text:
# Specifications for running multiple simulated environments

# For the Categorical task (Section 2.3.2)
mutable struct SpecsMultinomial
    stochasticities::Array{Float64, 1}
    changeprobabilities::Array{Float64, 1}
    nseeds::Int
    seedenvpolicy::Array{Any,1}
    seedlearner::Array{Any,1}
    ns::Int
    nsteps::Int
    envtype::Any
end
function SpecsMultinomial(; stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01],
                changeprobabilities = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0001],
                nseeds = 10,
                seedenvpolicy = collect(1:nseeds),
                # seedenvpolicy =  rand(1:typemax(UInt64)-1, nseeds),
                seedlearner =  collect(1:nseeds) * 99,
                # seedlearner =  rand(1:typemax(UInt64)-1, nseeds),
                ns = 5, nsteps = 10^5)
    envtype = Multi()
    SpecsMultinomial(stochasticities, changeprobabilities,
                    nseeds, seedenvpolicy, seedlearner,
                    ns, nsteps, envtype)
end

# For the Gaussian task (Section 2.3.1)
mutable struct SpecsGaussian
    mu_0::Float64
    sigma_0::Float64
    sigmas::Array{Float64, 1}
    changeprobabilities::Array{Float64, 1}
    nseeds::Int
    seedenvpolicy::Array{Any,1}
    seedlearner::Array{Any,1}
    nsteps::Int
    envtype::Any

end
function SpecsGaussian(; mu_0 = 0., sigma_0 = 1.,
                sigmas = [5., 2., 1., 0.5, 0.1, 0.01],
                changeprobabilities = [0.1, 0.05, 0.01, 0.005, 0.001, 0.0001],
                nseeds = 10,
                seedenvpolicy = collect(1:nseeds),
                #seedenvpolicy =  rand(1:typemax(UInt64)-1, nseeds),
                seedlearner =  collect(1:nseeds) * 99,
                # seedlearner =  rand(1:typemax(UInt64)-1, nseeds),
                nsteps = 10^5)
    envtype = Gaussian()
    SpecsGaussian(mu_0, sigma_0, sigmas, changeprobabilities,
                nseeds, seedenvpolicy, seedlearner,
                nsteps, envtype)
end

struct Multi end
Multi()
struct Gaussian end
Gaussian()


text:
module TestSpectralMixture
using CovarianceFunctions
using CovarianceFunctions: parameters, nparameters, isstationary
using Test
# IDEA: test on real data
# dataset:
# https://cdiac.ess-dive.lbl.gov/ftp/trends/co2/maunaloa.co2
@testset "spectral mixture kernel" begin
    n = 32
    x = sort!(randn(n))
    w = 1.
    μ = 0.
    l = 1.
    θ = [w, μ, l]
    k = CovarianceFunctions.Spectral(w, μ, l)
    @test k isa CovarianceFunctions.Product
    @test k(x[1], x[2]) isa Real
    @test isstationary(k)
    @test isstationary(k, x)

    ns = 3
    w = exp.(randn(ns))
    μ = exp.(randn(ns))
    l = exp.(randn(ns))
    θ = [w μ l]
    θ = θ'[:]
    sm = CovarianceFunctions.SpectralMixture(w, μ, l)
    @test isstationary(sm)
    @test isstationary(sm, x)

    # testing periodic
    # x = -1:.01:1
    # eq = CovarianceFunctions.EQ()
    # # eq = CovarianceFunctions.Lengthscale(eq, .5) # changes variation of periodic function without changing periodicity
    # per = CovarianceFunctions.Periodic(eq)
    # per = CovarianceFunctions.Lengthscale(per, .5) # changes periodicity
    # σ = .001
    # C = factorize(σ^2*I + CovarianceFunctions.gramian(per, x))
    # y = C.L * randn(length(x), 2)
    # using Plots
    # plot(x, y)
    # gui()
end

end # TestSpectralMixture

ext: jl


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    FunctionizeBridge{T}

`FunctionizeBridge` implements the following reformulations:

 * ``\\min \\{x\\}`` into ``\\min\\{1x + 0\\}``
 * ``\\max \\{x\\}`` into ``\\max\\{1x + 0\\}``

where `T` is the coefficient type of `1` and `0`.

## Source node

`FunctionizeBridge` supports:

 * [`MOI.ObjectiveFunction{MOI.VariableIndex}`](@ref)

## Target nodes

`FunctionizeBridge` creates:

 * One objective node: [`MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}`](@ref)
"""
struct FunctionizeBridge{T} <: AbstractBridge end

const Functionize{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{FunctionizeBridge{T},OT}

function bridge_objective(
    ::Type{FunctionizeBridge{T}},
    model::MOI.ModelLike,
    func::MOI.VariableIndex,
) where {T}
    F = MOI.ScalarAffineFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), convert(F, func))
    return FunctionizeBridge{T}()
end

function supports_objective_function(
    ::Type{<:FunctionizeBridge},
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:FunctionizeBridge},
)
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(::Type{<:FunctionizeBridge})
    return Tuple{Type,Type}[]
end

function MOI.Bridges.set_objective_function_type(
    ::Type{FunctionizeBridge{T}},
) where {T}
    return MOI.ScalarAffineFunction{T}
end

# Attributes, Bridge acting as a model
MOI.get(::FunctionizeBridge, ::MOI.NumberOfVariables)::Int64 = 0

function MOI.get(::FunctionizeBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# No variables or constraints are created in this bridge so there is nothing to
# delete.
MOI.delete(::MOI.ModelLike, ::FunctionizeBridge) = nothing

function MOI.set(
    ::MOI.ModelLike,
    ::MOI.ObjectiveSense,
    ::FunctionizeBridge,
    ::MOI.OptimizationSense,
)
    # `FunctionizeBridge` is sense agnostic, therefore, we don't need to change
    # anything.
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.Bridges.ObjectiveFunctionValue{MOI.VariableIndex},
    ::FunctionizeBridge{T},
) where {T}
    F = MOI.ScalarAffineFunction{T}
    attr_f = MOI.Bridges.ObjectiveFunctionValue{F}(attr.result_index)
    return MOI.get(model, attr_f)
end

function MOI.get(
    model::MOI.ModelLike,
    ::MOI.ObjectiveFunction{MOI.VariableIndex},
    ::FunctionizeBridge{T},
) where {T}
    F = MOI.ScalarAffineFunction{T}
    func = MOI.get(model, MOI.ObjectiveFunction{F}())
    return convert(MOI.VariableIndex, func)
end


text:
# some noise model (standard deviation)

"""
    LinearNoise

The noise increase linearly with the distance to the object

#Fields
- close_range: distance at which the error starts increasing
- min_noise: minimum error
- increase_rate: noise dimension/meter (e.g. position noise increase by 1m every meter)
"""
@with_kw struct LinearNoise <: AbstractNoiseModel
    close_range::Float64 = 2.0
    min_noise::Float64 = 0.5
    increase_rate::Float64 = 0.1
end

function noise(model::LinearNoise, dist::Float64)
    return max(model.min_noise, model.min_noise + dist*model.increase_rate)
end


text:
module Processing
import FFTW
using EllipsisNotation
import Glob: glob
import Luna: Maths, Fields, PhysData
import Luna.PhysData: wlfreq, c
import Luna.Grid: AbstractGrid, RealGrid, EnvGrid, from_dict
import Luna.Output: AbstractOutput, HDF5Output

"""
    Common(val)

Wrapper type to tell `scanproc` that `val` is the same for each simulation being processed,
and so only needs to be returned once rather than for each simulation in the scan.
"""
struct Common{dT}
    data::dT
end

"""
    VarLength(val)

Wrapper type to tell `scanproc` that the shape of `val` is different for each simulation being
processed. Return values wrapped in `VarLength` will be placed in an array of arrays.

!!! note
    While the **shape** of `val` can be different between simulations, the **type** must be
    the same, including the dimensionality and element type of arrays.
"""
struct VarLength{dT}
    data::dT
end

"""
    scanproc(f, scanfiles)
    scanproc(f, directory)
    scanproc(f, directory, pattern)
    scanproc(f)

Iterate over the scan output files, apply the processing function `f(o::AbstractOutput)`,
and collect the results in arrays.

The files can be given as:

- a `Vector` of `AbstractString`s containing file paths
- a directory to search for files according to the naming pattern of
    `Output.ScanHDF5Output`
- a directory and a `glob` pattern

If nothing is specified, `scanproc` uses the current working directory.

`f` can return a single value, an array, or a tuple/array of arrays/numbers. Arrays returned
by `f` must either be of the same size for each processed file, or wrapped in a `VarLength`.
Values returned by `f` which are guaranteed to be identical for each processed file can be
wrapped in a `Common`, and `scanproc` only returns these once.

# Example
```julia
Et, Eω = scanproc("path/to/scandir") do output
    t, Et = getEt(output)
    ω, Eω = getEω(output)
    energyout = energyout = Processing.VarLength(output["stats"]["energy"])
    Common(t), Et, Common(ω), Eω, energyout
end
```
"""
function scanproc(f, scanfiles::AbstractVector{<:AbstractString}; shape=nothing)
    local scanidcs, arrays
    scanfiles = sort(scanfiles)
    for (idx, fi) in enumerate(scanfiles)
        o = HDF5Output(fi)
        # wraptuple makes sure we definitely have a Tuple, even if f only returns one thing
        ret = wraptuple(f(o))
        if idx == 1 # initialise arrays
            isnothing(shape) && (shape = Tuple(o["meta"]["scanshape"]))
            scanidcs = CartesianIndices(shape)
            arrays = _arrays(ret, shape)
        end
        for (ridx, ri) in enumerate(ret)
            _addret!(arrays[ridx], scanidcs[idx], ri)
        end
    end
    unwraptuple(arrays) # if f only returns one thing, we also only return one array
end

wraptuple(x::Tuple) = x
wraptuple(x) = (x,)

unwraptuple(x::Tuple{<:Any}) = x[1] # single-element Tuple
unwraptuple(x) = x

function _addret!(array, aidcs, ri)
    array[aidcs] = ri
end

function _addret!(array, aidcs, ri::AbstractArray)
    idcs = CartesianIndices(ri)
    array[idcs, aidcs] .= ri
end

function _addret!(array, aidcs, ri::VarLength)
    array[aidcs] = ri.data
end

_addret!(array, aidcs, ri::Common) = nothing

# Default pattern for files named by ScanHDF5Output is [name]_[scanidx].h5 with 5 digits
defpattern = "*_[0-9][0-9][0-9][0-9][0-9].h5"

function scanproc(f, directory::AbstractString=pwd(), pattern::AbstractString=defpattern;
                  shape=nothing)
    scanfiles = glob(pattern, directory) # this returns absolute paths if directory given
    scanproc(f, scanfiles; shape=shape)
end

# Make array(s) with correct size to hold processing results
_arrays(ret, shape) = Array{typeof(ret)}(undef, shape)
_arrays(ret::AbstractArray, shape) = zeros(eltype(ret), (size(ret)..., shape...))
_arrays(ret::Tuple, shape) = Tuple([_arrays(ri, shape) for ri in ret])
_arrays(com::Common, shape) = com.data
_arrays(vl::VarLength, shape) = Array{typeof(vl.data), length(shape)}(undef, shape)

"""
    coherence(Eω; ndim=1)

Calculate the first-order coherence function g₁₂ of the set of fields `Eω`. The ensemble
average is taken over the last `ndim` dimensions of `Eω`, other dimensions are preserved.

See J. M. Dudley and S. Coen, Optics Letters 27, 1180 (2002).
"""
function coherence(Eω; ndim=1)
    dimsize = size(Eω)[end-ndim+1:end]
    outsize = size(Eω)[1:end-ndim]
    prodidcs = CartesianIndices(dimsize)
    restidcs = CartesianIndices(outsize)
    coherence(Eω, prodidcs, restidcs)
end

# function barrier for speedup
function coherence(Eω, prodidcs, restidcs)
    num = zeros(ComplexF64, size(restidcs))
    den1 = zeros(ComplexF64, size(restidcs))
    den2 = zeros(ComplexF64, size(restidcs))
    it = Iterators.product(prodidcs, prodidcs)
    for (idx1, idx2) in it
        Eω1 = Eω[restidcs, idx1]
        Eω2 = Eω[restidcs, idx2]
        @. num += conj(Eω1)*Eω2
        @. den1 += abs2(Eω1)
        @. den2 += abs2(Eω2)
    end
    @. abs(num/sqrt(den1*den2))
end

"""
    arrivaltime(grid, Eω; bandpass=nothing, method=:moment, oversampling=1)

Extract the arrival time of the pulse in the wavelength limits `λlims`.

# Arguments
- `bandpass` : method to bandpass the field if required. See [`window_maybe`](@ref)
- `method::Symbol` : `:moment` to use 1st moment to extract arrival time, `:peak` to use
                    the time of peak power
- `oversampling::Int` : If >1, oversample the time-domain field before extracting delay
- `sumdims` : Single `Int` or `Tuple` of `Int`s. The time-domain power will be summed over
            these dimensions (e.g. modes) before extracting the arrival time.
"""
function arrivaltime(grid::AbstractGrid, Eω;
                     bandpass=nothing, method=:moment, oversampling=1, sumdims=nothing)
    to, Eto = getEt(grid, Eω; oversampling=oversampling, bandpass=bandpass)
    Pt = abs2.(Eto)
    if !isnothing(sumdims)
        Pt = dropdims(sum(Pt; dims=sumdims); dims=sumdims)
    end
    arrivaltime(to, Pt; method=method)
end

function arrivaltime(t::AbstractVector, It::AbstractVector; method)
    if method == :moment
        Maths.moment(t, It)
    elseif method == :peak
        t[argmax(It)]
    else
        error("Unknown arrival time method $method")
    end
end

function arrivaltime(t::AbstractVector, It::AbstractArray; method)
    out = Array{Float64, ndims(It)-1}(undef, size(It)[2:end])
    cidcs = CartesianIndices(size(It)[2:end])
    for ii in cidcs
        out[ii] = arrivaltime(t, It[:, ii]; method=method)
    end
    out
end

"""
    time_bandwidth(grid, Eω; bandpass=nothing, oversampling=1)

Extract the time-bandwidth product, after bandpassing if required. The TBP
is defined here as ΔfΔt where Δx is the FWHM of x. (In this definition, the TBP of
a perfect Gaussian pulse is ≈0.44). If `oversampling` > 1, the time-domain field is
oversampled before extracting the FWHM.
"""
function time_bandwidth(grid, Eω; bandpass=nothing, oversampling=1, sumdims=nothing)
    fwt = fwhm_t(grid, Eω; bandpass=bandpass, oversampling=oversampling, sumdims=nothing)
    fwf = fwhm_f(grid, Eω; bandpass=bandpass)
    fwt.*fwf
end

"""
    fwhm_t(grid::AbstractGrid, Eω; bandpass=nothing, oversampling=1)

Extract the temporal FWHM. If `bandpass` is given, bandpass the fieldaccording to
[`window_maybe`](@ref). If `oversampling` > 1, the  time-domain field is oversampled before
extracting the FWHM. If `sumdims` is given, the time-domain power is summed over these
dimensions (e.g. modes) before extracting the FWHM.
"""
function fwhm_t(grid::AbstractGrid, Eω; bandpass=nothing, oversampling=1, sumdims=nothing)
    to, Eto = getEt(grid, Eω; oversampling=oversampling, bandpass=bandpass)
    Pt = abs2.(Eto)
    if !isnothing(sumdims)
        Pt = dropdims(sum(Pt; dims=sumdims); dims=sumdims)
    end
    fwhm(to, Pt)
end

function fwhm_t(output::AbstractOutput; kwargs...)
    grid = makegrid(output)
    fwhm_t(grid, output["Eω"]; kwargs...)
end

"""
    fwhm_f(grid, Eω::Vector; bandpass=nothing, oversampling=1)

Extract the frequency FWHM. If `bandpass` is given, bandpass the field according to
[`window_maybe`](@ref). If `sumdims` is given, the energy density is summed over these
dimensions (e.g. modes) before extracting the FWHM.
"""
function fwhm_f(grid::AbstractGrid, Eω; bandpass=nothing, oversampling=1, sumdims=nothing)
    Eω = window_maybe(grid.ω, Eω, bandpass)
    f, If = getIω(getEω(grid, Eω)..., :f)
    if !isnothing(sumdims)
        If = dropdims(sum(If; dims=sumdims); dims=sumdims)
    end
    fwhm(f, If)
end

function fwhm(x, I)
    out = Array{Float64, ndims(I)-1}(undef, size(I)[2:end])
    cidcs = CartesianIndices(size(I)[2:end])
    for ii in cidcs
        out[ii] = fwhm(x, I[:, ii])
    end
    out
end

fwhm(x::Vector, I::Vector) = Maths.fwhm(x, I)

"""
    peakpower(grid, Eω; bandpass=nothing, oversampling=1)

Extract the peak power. If `bandpass` is given, bandpass the field according to
[`window_maybe`](@ref).
"""
function peakpower(grid, Eω; bandpass=nothing, oversampling=1)
    to, Eto = getEt(grid, Eω; oversampling=oversampling, bandpass=bandpass)
    dropdims(maximum(abs2.(Eto); dims=1); dims=1)
end

"""
    energy(grid, Eω; bandpass=nothing)

Extract energy. If `bandpass` is given, bandpass the field according to
[`window_maybe`](@ref).
"""
function energy(grid, Eω; bandpass=nothing)
    Eω = window_maybe(grid.ω, Eω, bandpass)
    _, energyω = Fields.energyfuncs(grid)
    _energy(Eω, energyω)
end

function energy(output::AbstractOutput; bandpass=nothing)
    grid = makegrid(output)
    energy(grid, output["Eω"]; bandpass=bandpass)
end

_energy(Eω::Vector, energyω) = energyω(Eω)

function _energy(Eω, energyω)
    out = Array{Float64, ndims(Eω)-1}(undef, size(Eω)[2:end])
    cidcs = CartesianIndices(size(Eω)[2:end])
    for ii in cidcs
        out[ii] = _energy(Eω[:, ii], energyω)
    end
    out
end

"""
    field_autocorrelation(Et; dims=1)

Calculate the field autocorrelation of `Et`.
"""
function field_autocorrelation(Et, grid::EnvGrid; dims=1)
    FFTW.fftshift(FFTW.ifft(abs2.(FFTW.fft(Et, dims)), dims), dims)
end

function field_autocorrelation(Et, grid::RealGrid; dims=1)
    fac = FFTW.fftshift(FFTW.irfft(abs2.(FFTW.rfft(Et, dims)), length(grid.t), dims), dims)
    Maths.hilbert(fac, dim=dims)
end

"""
    intensity_autocorrelation(Et, grid)

Calculate the intensity autocorrelation of `Et` over `grid`.
"""
function intensity_autocorrelation(Et, grid; dims=1)
    real.(FFTW.fftshift(FFTW.irfft(abs2.(FFTW.rfft(Fields.It(Et, grid), dims)), length(grid.t), dims), dims))
end

"""
    coherence_time(grid, Et; dims=1)

Get the coherence time of a field `Et` over `grid`.
"""
function coherence_time(grid, Et; dims=1)
    Maths.fwhm(grid.t, abs2.(field_autocorrelation(Et, grid, dims=dims)))
end

"""
    specres(ω, Iω, specaxis, resolution, specrange; window=nothing, nsamples=10)

Smooth the spectral energy density `Iω(ω)` to account for the given `resolution`
on the defined `specaxis` and `specrange`. The `window` function to use defaults
to a Gaussian function with FWHM of `resolution`, and by default we sample `nsamples=10`
times within each `resolution`.

Note that you should prefer the `resolution` keyword of [`getIω`](@ref) instead of calling
this function directly.

The input `ω` and `Iω` should be as returned by [`getIω`](@ref) with `specaxis = :ω`.

Returns the new specaxis grid and smoothed spectrum.
"""
function specres(ω, Iω, specaxis, resolution, specrange; window=nothing, nsamples=10)
    if isnothing(window)
        window = let ng=Maths.gaussnorm(fwhm=resolution), σ=resolution/(2*(2*log(2))^(1/2))
            (x,x0) -> exp(-0.5*((x - x0)/σ)^2)/ng
        end
    end
    if specaxis == :λ
        xg, Ix = _specres(ω, Iω, resolution, specrange, window, nsamples, wlfreq, wlfreq)
    elseif specaxis == :f
        xg, Ix = _specres(ω, Iω, resolution, specrange, window, nsamples, x -> x/(2π), x -> x*(2π))
    else
        error("`specaxis` must be one of `:λ` or `:f`")
    end
    xg, Ix
end

function _specres(ω, Iω, resolution, xrange, window, nsamples, ωtox, xtoω)
    # build output grid and array
    x = ωtox.(ω)
    fxrange = extrema(x[(x .> 0) .& isfinite.(x)])
    if isnothing(xrange)
        xrange = fxrange
    else
        xrange = extrema(xrange)
        xrange = (max(xrange[1], fxrange[1]), min(xrange[2], fxrange[2]))
    end
    nxg = ceil(Int, (xrange[2] - xrange[1])/resolution*nsamples)
    xg = collect(range(xrange[1], xrange[2], length=nxg))
    rdims = size(Iω)[2:end]
    Ix = Array{Float64, ndims(Iω)}(undef, ((nxg,)..., rdims...))
    fill!(Ix, 0.0)
    cidcs = CartesianIndices(rdims)
    # we find a suitable nspan
    nspan = 1
    while window(nspan*resolution, 0.0)/window(0.0, 0.0) > 1e-8
        nspan += 1
    end
    # now we build arrays of start and end indices for the relevant frequency
    # band for each output. For a frequency grid this is a little inefficient
    # but for a wavelength grid, which has varying index ranges, this is essential
    # and I think having a common code is simpler/cleaner.
    istart = Array{Int,1}(undef,nxg)
    iend = Array{Int,1}(undef,nxg)
    δω = ω[2] - ω[1]
    i0 = argmin(abs.(ω))
    ωs = ω[i0]
    for i in 1:nxg
        i1 = i0 + round(Int, (xtoω(xg[i] + resolution*nspan) - ωs)/δω)
        i2 = i0 + round(Int, (xtoω(xg[i] - resolution*nspan) - ωs)/δω)
        # we want increasing indices
        if i1 > i2
            i1,i2 = i2,i1
        end
        # handle boundaries
        if i2 > length(ω)
            i2 = length(ω)
        end
        if i1 < i0
            i1 = i0
        end
        istart[i] = i1
        iend[i] = i2
    end
    # run the convolution kernel - the function barrier massively improves performance
    _specres_kernel!(Ix, cidcs, istart, iend, Iω, window, x, xg, δω)
    xg, Ix
end

"""
Convolution kernel for each output point. We simply loop over all outer indices
and output points. The inner loop adds up the contributions from the specified window
around the target point. Note that this works without scaling also for wavelength ranges
because the integral is still over a frequency grid (with appropriate frequency dependent
integration bounds).
"""
function _specres_kernel!(Ix, cidcs, istart, iend, Iω, window, x, xg, δω)
    @inbounds @fastmath for ii in cidcs
        for j in 1:size(Ix, 1)
            for k in istart[j]:iend[j]
                Ix[j,ii] += Iω[k,ii] * window(x[k], xg[j]) * δω
            end
        end
    end
    Ix[Ix .<= 0.0] .= minimum(Ix[Ix .> 0.0])
end

function _specrangeselect(x, Ix; specrange=nothing, sortx=false)
    cidcs = CartesianIndices(size(Ix)[2:end])
    if !isnothing(specrange)
        specrange = extrema(specrange)
        idcs = ((x .>= specrange[1]) .& (x .<= specrange[2]))
        x = x[idcs]
        Ix = Ix[idcs, cidcs]
    end
    if sortx
        idcs = sortperm(x)
        x = x[idcs]
        Ix = Ix[idcs, cidcs]
    end
    x, Ix
end

"""
    ωwindow_λ(ω, λlims; winwidth=:auto)

Create a ω-axis filtering window to filter in `λlims`. `winwidth`, if a `Number`, sets
the smoothing width of the window in rad/s.
"""
function ωwindow_λ(ω, λlims; winwidth=:auto)
    ωmin, ωmax = extrema(wlfreq.(λlims))
    winwidth == :auto && (winwidth = 64*abs(ω[2] - ω[1]))
    window = Maths.planck_taper(ω, ωmin-winwidth, ωmin, ωmax, ωmax+winwidth)
end

"""
    getIω(ω, Eω, specaxis; specrange=nothing, resolution=nothing)

Get spectral energy density and x-axis given a frequency array `ω` and frequency-domain field
`Eω`, assumed to be correctly normalised (see [`getEω`](@ref)). `specaxis` determines the
x-axis:

- :f -> x-axis is frequency in Hz and Iω is in J/Hz
- :ω -> x-axis is angular frequency in rad/s and Iω is in J/(rad/s)
- :λ -> x-axis is wavelength in m and Iω is in J/m

# Keyword arguments
- `specrange::Tuple` can be set to a pair of limits on the spectral range (in `specaxis` units).
- `resolution::Real` is set, smooth the spectral energy density as defined by [`specres`](@ref).

Note that if `resolution` and `specaxis=:λ` is set it is highly recommended to also set `specrange`.
"""
function getIω(ω, Eω, specaxis; specrange=nothing, resolution=nothing)
    sortx = false
    if specaxis == :ω || !isnothing(resolution)
        specx = ω
        Ix = abs2.(Eω)
        if !isnothing(resolution)
            return specres(ω, Ix, specaxis, resolution, specrange)
        end
    elseif specaxis == :f
        specx = ω./2π
        Ix = abs2.(Eω)*2π
    elseif specaxis == :λ
        specx = wlfreq.(ω)
        Ix = @. ω^2/(2π*c) * abs2.(Eω)
        sortx = true
    else
        error("Unknown specaxis $specaxis")
    end
    if !isnothing(specrange) || sortx
        specx, Ix = _specrangeselect(specx, Ix, specrange=specrange, sortx=sortx)
    end
    return specx, Ix
end

"""
    getIω(output, specaxis[, zslice]; kwargs...)

Calculate the correctly normalised frequency-domain field and convert it to spectral
energy density on x-axis `specaxis` (`:f`, `:ω`, or `:λ`). If `zslice` is given,
returs only the slices of `Eω` closest to the given distances. `zslice` can be a single
number or an array. `specaxis` determines the
x-axis:

- :f -> x-axis is frequency in Hz and Iω is in J/Hz
- :ω -> x-axis is angular frequency in rad/s and Iω is in J/(rad/s)
- :λ -> x-axis is wavelength in m and Iω is in J/m

# Keyword arguments
- `specrange::Tuple` can be set to a pair of limits on the spectral range (in `specaxis` units).
- `resolution::Real` is set, smooth the spectral energy density as defined by [`specres`](@ref).

Note that `resolution` is set and `specaxis=:λ` it is highly recommended to also set `specrange`.
"""
getIω(output::AbstractOutput, specaxis; kwargs...) = getIω(getEω(output)..., specaxis; kwargs...)

function getIω(output::AbstractOutput, specaxis, zslice; kwargs...)
    ω, Eω, zactual = getEω(output, zslice)
    specx, Iω = getIω(ω, Eω, specaxis; kwargs...)
    return specx, Iω, zactual
end

"""
    getEω(output[, zslice])

Get frequency-domain modal field from `output` with correct normalisation (i.e.
`abs2.(Eω)`` gives angular-frequency spectral energy density in J/(rad/s)).
"""
getEω(output::AbstractOutput, args...) = getEω(makegrid(output), output, args...)
getEω(grid, output) = getEω(grid, output["Eω"])

function getEω(grid::RealGrid, Eω::AbstractArray)
    ω = grid.ω[grid.sidx]
    Eω = Eω[grid.sidx, CartesianIndices(size(Eω)[2:end])]
    return ω, Eω*fftnorm(grid)
end

function getEω(grid::EnvGrid, Eω::AbstractArray)
    idcs = FFTW.fftshift(grid.sidx)
    Eωs = FFTW.fftshift(Eω, 1)
    ω = FFTW.fftshift(grid.ω)[idcs]
    Eω = Eωs[idcs, CartesianIndices(size(Eω)[2:end])]
    return ω, Eω*fftnorm(grid)
end

function getEω(grid, output, zslice)
    zidx = nearest_z(output, zslice)
    ω, Eω = getEω(grid, output["Eω", .., zidx])
    return ω, Eω, output["z"][zidx]
end

fftnorm(grid::RealGrid) = Maths.rfftnorm(grid.t[2] - grid.t[1])
fftnorm(grid::EnvGrid) = Maths.fftnorm(grid.t[2] - grid.t[1])

"""
    getEt(output[, zslice]; kwargs...)

Get the envelope time-domain electric field (including the carrier wave) from the `output`.
If `zslice` is given, returs only the slices of `Eω` closest to the given distances. `zslice`
can be a single number or an array.
"""
getEt(output::AbstractOutput, args...; kwargs...) = getEt(
    makegrid(output), output, args...; kwargs...)

"""
    getEt(grid, Eω; trange=nothing, oversampling=4, bandpass=nothing, FTL=false)

Get the envelope time-domain electric field (including the carrier wave) from the frequency-
domain field `Eω`. The field can be cropped in time using `trange`, it is oversampled by
a factor of `oversampling` (default 4) and can be bandpassed with `bandpass`
(see [`window_maybe`](@ref)). If `FTL` is `true`, return the Fourier-transform limited pulse,
i.e. remove any spectral phase.

If `zslice` is given, returs only the slices of `Eω` closest to the given distances. `zslice`
can be a single number or an array.
"""
function getEt(grid::AbstractGrid, Eω::AbstractArray;
               trange=nothing, oversampling=4, bandpass=nothing,
               FTL=false, propagate=nothing)
    t = grid.t
    Eω = window_maybe(grid.ω, Eω, bandpass)
    if FTL
        τ = length(grid.t) * (grid.t[2] - grid.t[1])/2
        Eω .= abs.(Eω) .* exp.(-1im .* grid.ω .* τ)
    end
    Eω = prop_maybe(grid, Eω, propagate)
    Etout = envelope(grid, Eω)
    if isnothing(trange)
        idcs = 1:length(t)
    else
        idcs = @. (t < max(trange...)) & (t > min(trange...))
    end
    to, Eto = Maths.oversample(t[idcs], Etout[idcs, ..], factor=oversampling)
    return to, Eto
end

getEt(grid::AbstractGrid, output::AbstractOutput; kwargs...) = getEt(grid, output["Eω"]; kwargs...)

function getEt(grid::AbstractGrid, output::AbstractOutput, zslice;
               kwargs...)
    t = grid.t
    zidx = nearest_z(output, zslice)
    to, Eto = getEt(grid, output["Eω", .., zidx]; kwargs...)
    return to, Eto, output["z"][zidx]
end

"""
    AutoWindow(width, λmin, λmax, ω0fun; relative=false, ndims=1)

Window function generator which automatically tracks the central frequency in the spectral
region given by `λmin` and `λmax` and applies a window of a specific `width` around the peak.
The central frequency is found using the function `ω0fun(ω, Iω::AbstractVector)`, where
`ω` and `Iω` are already cropped to within the wavelength limits given.
If `relative` is `true`, `width` is relative bandwidth instead of the wavelength width.
`ndims` determines how many dimensions of the array to sum over. For a field array with size
`(Nω, N1, N2, ...)`, the first dimension is always assumed to be frequency. `ndim=1` means
each field to be analysed is 1-dimensional, so the window iterates over all of `(N1, N2, ...)`.
`ndim=2` means each field to be analysed is 2-dimensional, `(Nω, N1)` in size, and will be
summed over its second dimension before finding the central frequency. The window iterates
over all other dimensions, `(N2, ...)`.

A `AutoWindow` automatically stores the limits of the windows it applies in the field `lims`.
"""
mutable struct AutoWindow
    width::Float64
    λmin::Float64
    λmax::Float64
    ω0fun
    relative::Bool
    ndims
    lims
end

function AutoWindow(width, λmin, λmax, ω0fun; relative=false, ndims=1)
    AutoWindow(width, λmin, λmax, ω0fun, relative, ndims, nothing)
end

function (pw::AutoWindow)(ω, Eω)
    cidcs = CartesianIndices(size(Eω)[(pw.ndims+1):end])
    out = similar(Eω)
    cropidcs = (ω .> wlfreq(pw.λmax)) .& (ω .< wlfreq(pw.λmin))
    cropω = ω[cropidcs]
    Iω = abs2.(Eω)
    limsA = zeros((2, size(Eω)[(pw.ndims+1):end]...))
    for cidx in cidcs
        Iω_this = Iω[.., Tuple(cidx)...]
        Iωsum = sum(Iω_this; dims=2:ndims(Iω_this))
        λ0 = wlfreq(pw.ω0fun(cropω, Iωsum[cropidcs]))
        lims = pw.relative ? λ0.*(1 .+ (-0.5, 0.5).*pw.width) : λ0 .+ (-0.5, 0.5).*pw.width
        window = ωwindow_λ(ω, lims)
        limsA[:, Tuple(cidx)...] .= lims
        out[.., Tuple(cidx)...] .= Eω[.., Tuple(cidx)...] .* window
    end
    pw.lims = limsA
    out
end

"""
    PeakWindow(width, λmin, λmax; relative=false, ndims=1)

An [`AutoWindow`](@ref) which uses the peak of the spectral energy density as the central
frequency.
"""
function PeakWindow(width, λmin, λmax; relative=false, ndims=1)
    ω0fun = (ω, Iω) ->  ω[argmax(Iω)]
    AutoWindow(width, λmin, λmax, ω0fun; relative=relative, ndims=ndims)
end

"""
    CentroidWindow(width, λmin, λmax; relative=false, ndims=1, power=1)

An [`AutoWindow`](@ref) which uses the centroid (centre of mass or first moment) of the
spectral energy density as the central frequency. Before calculating the centroid, the
SED is raised to the `power` given.
"""
function CentroidWindow(width, λmin, λmax; relative=false, ndims=1, power=1)
    ω0fun = (ω, Iω) -> Maths.moment(ω, Iω.^power)
    AutoWindow(width, λmin, λmax, ω0fun; relative=relative, ndims=ndims)
end

"""
    window_maybe(ω, Eω, win)

Apply a frequency window to the field `Eω` if required. Possible values for `win`:

- `nothing` : no window is applied
- 4-`Tuple` of `Number`s : the 4 parameters for a `Maths.planck_taper` in **wavelength**
- 3-`Tuple` of `Number`s : minimum, maximum **wavelength**, and smoothing in **radial frequency**
- 2-`Tuple` of `Number`s : minimum and maximum **wavelength** with automatically chosen smoothing
- `Vector{<:Real}` : a pre-defined window function (shape must match `ω`)
- `PeakWindow` : automatically track the peak in a given range and apply the window around it
- `window(ω, Eω)` : an arbitrary user-supplied window function
"""
window_maybe(ω, Eω, ::Nothing) = Eω
window_maybe(ω, Eω, win::NTuple{4, Number}) = Eω.*Maths.planck_taper(
    ω, sort(wlfreq.(collect(win)))...)
window_maybe(ω, Eω, win::NTuple{2, Number}) = Eω .* ωwindow_λ(ω, win)
window_maybe(ω, Eω, win::NTuple{3, Number}) = Eω .* ωwindow_λ(ω, win[1:2]; winwidth=win[3])
window_maybe(ω, Eω, window) = window(ω, Eω)
window_maybe(ω, Eω, window::AbstractVector) = Eω.*window

prop_maybe(grid, Eω, ::Nothing) = Eω
prop_maybe(grid, Eω, propagator) = propagator(grid, Eω)

"""
    envelope(grid, Eω)

Get the envelope electric field including the carrier wave from the frequency-domain field
`Eω` sampled on `grid`.
"""
envelope(grid::RealGrid, Eω) = Maths.hilbert(FFTW.irfft(Eω, length(grid.t), 1))
envelope(grid::EnvGrid, Eω) = FFTW.ifft(Eω, 1) .* exp.(im.*grid.ω0.*grid.t)

"""
    makegrid(output)

Create an `AbstractGrid` from the `"grid"` dictionary saved in `output`.
"""
function makegrid(output)
    if output["simulation_type"]["field"] == "field-resolved"
        from_dict(RealGrid, output["grid"])
    else
        from_dict(EnvGrid, output["grid"])
    end
end

"""
    nearest_z(output, z)

Return the index of saved z-position(s) closest to the position(s) `z`. Output is always
an array, even if `z` is a number.
"""
nearest_z(output, z::Number) = [argmin(abs.(output["z"] .- z))]
nearest_z(output, z) = [argmin(abs.(output["z"] .- zi)) for zi in z]

end


text:
module SpecialUpdates

using Random
using Statistics
using Parameters
using Printf
using LinearAlgebra
using SparseArrays
using StatsBase

using ..Utilities: get_index, reshaped, swap!
using ..Models: AbstractModel, HolsteinModel, SSHModel, update_model!, mulM!, muldMdx!, mulMᵀ!, construct_M
using ..PhononAction: calc_Sb
using ..HMC: HybridMonteCarlo, refresh_ϕ!, calc_S, calc_Sf, calc_O⁻¹Λϕ!, update_Λ!

export SpecialUpdate, NullUpdate, ReflectionUpdate, SwapUpdate, special_update!

"""
Abstract type to represent special updates.
"""
abstract type SpecialUpdate end

"""
Represent Null Special Update, meant to act as a placeholder in the code
when no special update is being used.
"""
struct NullUpdate <: SpecialUpdate

    """
    Whether the updater is active (always false).
    """
    active::Bool

    """
    Frequency of updates (set to 1 because never active).
    """
    freq::Int

    function NullUpdate()

        return new(false,1)
    end
end

"""
NullUpdate behavior.
"""
function special_update!(model::AbstractModel{T},hmc::HybridMonteCarlo{T},nu::NullUpdate,preconditioner)::T where {T<:AbstractFloat}

    return 0.0
end

"""
Performs Reflection update on Holstein model where (xᵢ)⟶(-xᵢ) on some set of sites i.
On all other models it is a null operation.
"""
mutable struct ReflectionUpdate <: SpecialUpdate

    """
    Whether the updater is turned on.
    """
    active::Bool

    """
    Frequency of special updates.
    """
    freq::Int

    """
    Number of sites the reflection will be applied to.
    """
    nsites::Int

    """
    Which sites the reflection will be applied to.
    """
    sites::Vector{Int}
end

function ReflectionUpdate(model::HolsteinModel,freq::Int,nsites::Int)

    nsites = min(model.Nph,nsites)
    sites = zeros(Int,nsites)
    return ReflectionUpdate(true,freq,nsites,sites)
end

function ReflectionUpdate(model::AbstractModel,freq::Int,nsites::Int)

    sites = Vector{Int}(undef,0)
    return ReflectionUpdate(false,freq,0,sites)
end

"""
Apply reflection updates to Holstein model.
"""
function special_update!(model::HolsteinModel{T},hmc::HybridMonteCarlo{T},ru::ReflectionUpdate,preconditioner)::T where {T<:AbstractFloat}

    @unpack Nph, Lτ = model
    @unpack nsites, sites = ru
    @unpack ϕ₊, ϕ₋ = hmc

    # counts number of accepted reflections
    accepted = 0.0

    # if updater is active
    if ru.active

        # get all phonon fields
        x = reshaped(model.x,Lτ,Nph)

        # randomly sample sites
        sample!(model.rng,1:Nph,sites,replace=true)

        # update exp{-Δτ⋅V[x]}
        update_model!(model)

        # iterate over sites to apply reflection operation to
        for i in sites

            # get phonon fields associated with site
            xᵢ = @view x[:,i]

            # CALCULATE FORWARD TRANSITION PROBABILITY

            # resample ϕ and calculate initial action
            S₀ = refresh_ϕ!(hmc,model,sample_R=true)

            # reflect phonon fields
            @. xᵢ = -xᵢ

            # update exp{-Δτ⋅V[x]}
            update_model!(model)

            # calculate O⁻¹⋅Λ⋅ϕ₊ and O⁻¹⋅Λ⋅ϕ₋
            iters, flag = calc_O⁻¹Λϕ!(hmc,model,preconditioner,2.0)

            # get final action
            S₁ = calc_S(hmc,model)

            # forward transition probability
            Pf = min( 1.0 , exp(-(S₁-S₀)) )

            # accept/reject decision
            if rand(model.rng) < Pf && iszero(flag)

                accepted += 1.0
            else

                # reflect phonon fields
                @. xᵢ = -xᵢ

                # update exp{-Δτ⋅V[x]}
                update_model!(model)
            end
        end
    end

    return accepted/nsites
end

function special_update!(model::AbstractModel{T},hmc::HybridMonteCarlo{T},ru::ReflectionUpdate,preconditioner)::T where {T<:AbstractFloat}

    return 0.0
end

"""
Swap update in Holstein model where adjacent sites swap phonon worldlines.
"""
mutable struct SwapUpdate <: SpecialUpdate

    """
    Whether is turned on.
    """
    active::Bool

    """
    Frequency of special update.
    """
    freq::Int

    """
    Number bonds to swap phonon position across.
    """
    nbonds::Int

    """
    The sites to apply the swap updates to.
    """
    bonds::Vector{Int}
end

function SwapUpdate(model::HolsteinModel,freq::Int,nbonds::Int)

    @unpack Nsites, Nbonds, neighbor_table = model

    if Nbonds==0
        active = false
    else
        active = true
    end
    nbonds = min(model.Nbonds,nbonds)
    bonds  = zeros(Int,nbonds)
    return SwapUpdate(active,freq,nbonds,bonds)
end

function SwapUpdate(model::AbstractModel,freq::Int,nbonds::Int)

    nbonds    = 0
    bonds     = Vector{Int}(undef,nbonds)
    return SwapUpdate(active,freq,nbonds,bonds)
end

"""
Apply swap updates to Holstein model.
"""
function special_update!(model::HolsteinModel{T},hmc::HybridMonteCarlo{T},su::SwapUpdate,preconditioner)::T where {T<:AbstractFloat}

    @unpack Nbonds, Nsites, Lτ, neighbor_table = model
    @unpack nbonds, bonds = su

    # counts number of accepted reflections
    accepted = 0.0

    # if updater is active
    if su.active

        # get all phonon fields
        x = reshaped(model.x,Lτ,Nsites)

        # randomly sample sites
        sample!(model.rng,1:Nbonds,bonds,replace=true)

        # update exp{-Δτ⋅V[x]}
        update_model!(model)

        # iterate over sites to apply reflection operation to
        for b in bonds

            # get a randomly selected neighbor table
            i = neighbor_table[1,b]
            j = neighbor_table[2,b]

            # get phonon fields associated with each sites
            xᵢ = @view x[:,i]
            xⱼ = @view x[:,j]

            # CALCULATE FORWARD TRANSITION PROBABILITY

            # resample ϕ and calculate initial action
            S₀ = refresh_ϕ!(hmc,model,sample_R=true)

            # swap mean phonon positions
            swap!(xᵢ,xⱼ)

            # update exp{-Δτ⋅V[x]}
            update_model!(model)

            # calculate O⁻¹⋅Λ⋅ϕ₊ and O⁻¹⋅Λ⋅ϕ₋
            iters, flag = calc_O⁻¹Λϕ!(hmc,model,preconditioner,2.0)

            # get final action
            S₁ = calc_S(hmc,model)

            # forward transition proability
            Pf = min( 1.0 , exp(-(S₁-S₀)) )

            # accept/reject decision
            if rand(model.rng) < Pf && iszero(flag)

                accepted += 1.0
            else

                # swap phonon positions
                swap!(xᵢ,xⱼ)

                # update exp{-Δτ⋅V[x]}
                update_model!(model)
            end
        end
    end

    return accepted/nbonds
end

function special_update!(model::AbstractModel{T},hmc::HybridMonteCarlo{T},ru::SwapUpdate,preconditioner)::T where {T<:AbstractFloat}

    return 0.0
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
# test memory allocation

using MicrochipRTDM
using FTD2XX

projectdirectory = "c:/Users/Chris/MPLABXProjects/lbcmcp/PWM_RTDM.x"
mapdict = microchip_parsemap(projectdirectory,"production")
uartconfig1 = UARTConfiguration(921600,8,1,"n",readtimeout = 1000,writetimeout = 1000)
io = open(FT_DeviceIndex(0),uartconfig1)
setlatencytimer(io,2); #Should speed things up a little
rtdminterface = RTDMInterface(mapdict,io)

const address = mapdict["offsetP3"]
const buffer = Array(UInt16,1)
function f()
  for i in 1:1000
    rtdm_read!(rtdminterface, buffer, address)
    # x = rtdm_read(rtdminterface, UInt16, address)
  end
  return nothing
end

f()

Profile.clear_malloc_data()
f()
#rtdm_read!(io, buffer, address)

# close(rtdminterface)
quit()


text:
"""
Covalent radius in Å of each element from the Additional file 1 of PICCOLO [1].
Hydrogen was updated using the value on Table 2 from Cordero et. al. [2].

1. Bickerton, G. R., Higueruelo, A. P., & Blundell, T. L. (2011).
Comprehensive, atomic-level characterization of structurally
characterized protein-protein interactions: the PICCOLO database.
BMC bioinformatics, 12(1), 313.
2. Cordero, B., Gómez, V., Platero-Prats, A. E., Revés, M.,
Echeverría, J., Cremades, E., ... & Alvarez, S. (2008).
Covalent radii revisited. Dalton Transactions, (21), 2832-2838.
"""
const covalentradius = Dict{ASCIIString,Float64}("C" => 0.77,
                                                 "N" => 0.70,
                                                 "O" => 0.66,
                                                 "S" => 1.04,
                                                 "H" => 0.31 )

const _3_letter_aa = ASCIIString[ "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL" ]

"""
van der Waals radius in Å from the Additional file 1 of Bickerton et. al. 2011

- Bickerton, G. R., Higueruelo, A. P., & Blundell, T. L. (2011).
Comprehensive, atomic-level characterization of structurally characterized protein-protein interactions: the PICCOLO database.
BMC bioinformatics, 12(1), 313.
"""
const vanderwaalsradius = Dict{Tuple{ASCIIString,ASCIIString},Float64}(("ALA","C") => 1.61,
                                                                       ("ALA","CA") => 1.88,
                                                                       ("ALA","CB") => 1.88,
                                                                       ("ALA","N") => 1.64,
                                                                       ("ALA","O") => 1.42,
                                                                       ("ARG","C") => 1.61,
                                                                       ("ARG","CA") => 1.88,
                                                                       ("ARG","CB") => 1.88,
                                                                       ("ARG","CD") => 1.88,
                                                                       ("ARG","CG") => 1.88,
                                                                       ("ARG","CZ") => 1.61,
                                                                       ("ARG","N") => 1.64,
                                                                       ("ARG","NE") => 1.64,
                                                                       ("ARG","NH1") => 1.64,
                                                                       ("ARG","NH2") => 1.64,
                                                                       ("ARG","O") => 1.42,
                                                                       ("ASN","C") => 1.61,
                                                                       ("ASN","CA") => 1.88,
                                                                       ("ASN","CB") => 1.88,
                                                                       ("ASN","CG") => 1.61,
                                                                       ("ASN","N") => 1.64,
                                                                       ("ASN","ND2") => 1.64,
                                                                       ("ASN","O") => 1.42,
                                                                       ("ASN","OD1") => 1.42,
                                                                       ("ASP","C") => 1.61,
                                                                       ("ASP","CA") => 1.88,
                                                                       ("ASP","CB") => 1.88,
                                                                       ("ASP","CG") => 1.61,
                                                                       ("ASP","N") => 1.64,
                                                                       ("ASP","O") => 1.42,
                                                                       ("ASP","OD1") => 1.42,
                                                                       ("ASP","OD2") => 1.42,
                                                                       ("CYS","C") => 1.61,
                                                                       ("CYS","CA") => 1.88,
                                                                       ("CYS","CB") => 1.88,
                                                                       ("CYS","N") => 1.64,
                                                                       ("CYS","O") => 1.42,
                                                                       ("CYS","SG") => 1.77,
                                                                       ("GLN","C") => 1.61,
                                                                       ("GLN","CA") => 1.88,
                                                                       ("GLN","CB") => 1.88,
                                                                       ("GLN","CD") => 1.61,
                                                                       ("GLN","CG") => 1.88,
                                                                       ("GLN","N") => 1.64,
                                                                       ("GLN","NE2") => 1.64,
                                                                       ("GLN","O") => 1.42,
                                                                       ("GLN","OE1") => 1.42,
                                                                       ("GLU","C") => 1.61,
                                                                       ("GLU","CA") => 1.88,
                                                                       ("GLU","CB") => 1.88,
                                                                       ("GLU","CD") => 1.61,
                                                                       ("GLU","CG") => 1.88,
                                                                       ("GLU","N") => 1.64,
                                                                       ("GLU","O") => 1.42,
                                                                       ("GLU","OE1") => 1.42,
                                                                       ("GLU","OE2") => 1.42,
                                                                       ("GLY","C") => 1.61,
                                                                       ("GLY","CA") => 1.88,
                                                                       ("GLY","N") => 1.64,
                                                                       ("GLY","O") => 1.42,
                                                                       ("HIS","C") => 1.61,
                                                                       ("HIS","CA") => 1.88,
                                                                       ("HIS","CB") => 1.88,
                                                                       ("HIS","CD2") => 1.76,
                                                                       ("HIS","CE1") => 1.76,
                                                                       ("HIS","CG") => 1.61,
                                                                       ("HIS","N") => 1.64,
                                                                       ("HIS","ND1") => 1.64,
                                                                       ("HIS","NE2") => 1.64,
                                                                       ("HIS","O") => 1.42,
                                                                       ("ILE","C") => 1.61,
                                                                       ("ILE","CA") => 1.88,
                                                                       ("ILE","CB") => 1.88,
                                                                       ("ILE","CD1") => 1.88,
                                                                       ("ILE","CG1") => 1.88,
                                                                       ("ILE","CG2") => 1.88,
                                                                       ("ILE","N") => 1.64,
                                                                       ("ILE","O") => 1.42,
                                                                       ("LEU","C") => 1.61,
                                                                       ("LEU","CA") => 1.88,
                                                                       ("LEU","CB") => 1.88,
                                                                       ("LEU","CD1") => 1.88,
                                                                       ("LEU","CD2") => 1.88,
                                                                       ("LEU","CG") => 1.88,
                                                                       ("LEU","N") => 1.64,
                                                                       ("LEU","O") => 1.42,
                                                                       ("LYS","C") => 1.61,
                                                                       ("LYS","CA") => 1.88,
                                                                       ("LYS","CB") => 1.88,
                                                                       ("LYS","CD") => 1.88,
                                                                       ("LYS","CE") => 1.88,
                                                                       ("LYS","CG") => 1.88,
                                                                       ("LYS","N") => 1.64,
                                                                       ("LYS","NZ") => 1.64,
                                                                       ("LYS","O") => 1.42,
                                                                       ("MET","C") => 1.61,
                                                                       ("MET","CA") => 1.88,
                                                                       ("MET","CB") => 1.88,
                                                                       ("MET","CE") => 1.88,
                                                                       ("MET","CG") => 1.88,
                                                                       ("MET","N") => 1.64,
                                                                       ("MET","O") => 1.42,
                                                                       ("MET","SD") => 1.77,
                                                                       ("PHE","C") => 1.61,
                                                                       ("PHE","CA") => 1.88,
                                                                       ("PHE","CB") => 1.88,
                                                                       ("PHE","CD1") => 1.76,
                                                                       ("PHE","CD2") => 1.76,
                                                                       ("PHE","CE1") => 1.76,
                                                                       ("PHE","CE2") => 1.76,
                                                                       ("PHE","CG") => 1.61,
                                                                       ("PHE","CZ") => 1.76,
                                                                       ("PHE","N") => 1.64,
                                                                       ("PHE","O") => 1.42,
                                                                       ("PRO","C") => 1.61,
                                                                       ("PRO","CA") => 1.88,
                                                                       ("PRO","CB") => 1.88,
                                                                       ("PRO","CD") => 1.88,
                                                                       ("PRO","CG") => 1.88,
                                                                       ("PRO","N") => 1.64,
                                                                       ("PRO","O") => 1.42,
                                                                       ("SER","C") => 1.61,
                                                                       ("SER","CA") => 1.88,
                                                                       ("SER","CB") => 1.88,
                                                                       ("SER","N") => 1.64,
                                                                       ("SER","O") => 1.42,
                                                                       ("SER","OG") => 1.46,
                                                                       ("THR","C") => 1.61,
                                                                       ("THR","CA") => 1.88,
                                                                       ("THR","CB") => 1.88,
                                                                       ("THR","CG2") => 1.88,
                                                                       ("THR","N") => 1.64,
                                                                       ("THR","O") => 1.42,
                                                                       ("THR","OG1") => 1.46,
                                                                       ("TRP","C") => 1.61,
                                                                       ("TRP","CA") => 1.88,
                                                                       ("TRP","CB") => 1.88,
                                                                       ("TRP","CD1") => 1.76,
                                                                       ("TRP","CD2") => 1.61,
                                                                       ("TRP","CE2") => 1.61,
                                                                       ("TRP","CE3") => 1.76,
                                                                       ("TRP","CG") => 1.61,
                                                                       ("TRP","CH2") => 1.76,
                                                                       ("TRP","CZ2") => 1.76,
                                                                       ("TRP","CZ3") => 1.76,
                                                                       ("TRP","N") => 1.64,
                                                                       ("TRP","NE1") => 1.64,
                                                                       ("TRP","O") => 1.42,
                                                                       ("TYR","C") => 1.61,
                                                                       ("TYR","CA") => 1.88,
                                                                       ("TYR","CB") => 1.88,
                                                                       ("TYR","CD1") => 1.76,
                                                                       ("TYR","CD2") => 1.76,
                                                                       ("TYR","CE1") => 1.76,
                                                                       ("TYR","CE2") => 1.76,
                                                                       ("TYR","CG") => 1.61,
                                                                       ("TYR","CZ") => 1.61,
                                                                       ("TYR","N") => 1.64,
                                                                       ("TYR","O") => 1.42,
                                                                       ("TYR","OH") => 1.46,
                                                                       ("VAL","C") => 1.61,
                                                                       ("VAL","CA") => 1.88,
                                                                       ("VAL","CB") => 1.88,
                                                                       ("VAL","CG1") => 1.88,
                                                                       ("VAL","CG2") => 1.88,
                                                                       ("VAL","N") => 1.64,
                                                                       ("VAL","O") => 1.42 )

function _add_CTER_O!(dict)
    for aa in _3_letter_aa
        push!(dict, (aa, "OXT"))
        push!(dict, (aa, "OT2"))
        push!(dict, (aa, "OT1"))
    end
    dict
end

function _add_CTER_O!(dict, value)
    for aa in _3_letter_aa
        push!(dict, (aa, "OXT") => value)
        push!(dict, (aa, "OT2") => value)
        push!(dict, (aa, "OT1") => value)
    end
    dict
end

_add_CTER_O!(vanderwaalsradius, 1.42) # Using 1.42 because OXT is the terminal oxigen I assume same vdw radii than O

const _hydrophobic = Set{Tuple{ASCIIString, ASCIIString}}( [ ("ALA","CB"),
                                                            ("ARG","CB"),
                                                            ("ARG","CG"),
                                                            ("ASN","CB"),
                                                            ("ASP","CB"),
                                                            ("CYS","CB"),
                                                            ("GLN","CB"),
                                                            ("GLN","CG"),
                                                            ("GLU","CB"),
                                                            ("GLU","CG"),
                                                            ("HIS","CB"),
                                                            ("ILE","CB"),
                                                            ("ILE","CD1"),
                                                            ("ILE","CG1"),
                                                            ("ILE","CG2"),
                                                            ("LEU","CB"),
                                                            ("LEU","CD1"),
                                                            ("LEU","CD2"),
                                                            ("LEU","CG"),
                                                            ("LYS","CB"),
                                                            ("LYS","CD"),
                                                            ("LYS","CG"),
                                                            ("MET","CB"),
                                                            ("MET","CE"),
                                                            ("MET","CG"),
                                                            ("MET","SD"),
                                                            ("PHE","CB"),
                                                            ("PHE","CD1"),
                                                            ("PHE","CD2"),
                                                            ("PHE","CE1"),
                                                            ("PHE","CE2"),
                                                            ("PHE","CG"),
                                                            ("PHE","CZ"),
                                                            ("PRO","CB"),
                                                            ("PRO","CG"),
                                                            ("THR","CG2"),
                                                            ("TRP","CB"),
                                                            ("TRP","CD2"),
                                                            ("TRP","CE3"),
                                                            ("TRP","CG"),
                                                            ("TRP","CH2"),
                                                            ("TRP","CZ2"),
                                                            ("TRP","CZ3"),
                                                            ("TYR","CB"),
                                                            ("TYR","CD1"),
                                                            ("TYR","CD2"),
                                                            ("TYR","CE1"),
                                                            ("TYR","CE2"),
                                                            ("TYR","CG"),
                                                            ("VAL","CB"),
                                                            ("VAL","CG1"),
                                                            ("VAL","CG2") ] )

const _aromatic_res = Set{ASCIIString}( [ "HIS", "PHE", "TRP", "TYR" ] )

const _aromatic = Set{Tuple{ASCIIString, ASCIIString}}( [ ("HIS","CD2"),
                                                         ("HIS","CE1"),
                                                         ("HIS","CG"),
                                                         ("HIS","ND1"),
                                                         ("HIS","NE2"),
                                                         ("PHE","CD1"),
                                                         ("PHE","CD2"),
                                                         ("PHE","CE1"),
                                                         ("PHE","CE2"),
                                                         ("PHE","CG"),
                                                         ("PHE","CZ"),
                                                         ("TRP","CD1"),
                                                         ("TRP","CD2"),
                                                         ("TRP","CE2"),
                                                         ("TRP","CE3"),
                                                         ("TRP","CG"),
                                                         ("TRP","CH2"),
                                                         ("TRP","CZ2"),
                                                         ("TRP","CZ3"),
                                                         ("TRP","NE1"),
                                                         ("TYR","CD1"),
                                                         ("TYR","CD2"),
                                                         ("TYR","CE1"),
                                                         ("TYR","CE2"),
                                                         ("TYR","CG"),
                                                         ("TYR","CZ") ] )

const _cationic = Set{Tuple{ASCIIString, ASCIIString}}( [ ("ARG","CZ"),
                                                         ("ARG","NE"),
                                                         ("ARG","NH1"),
                                                         ("ARG","NH2"),
                                                         ("HIS","CD2"),
                                                         ("HIS","CE1"),
                                                         ("HIS","CG"),
                                                         ("HIS","ND1"),
                                                         ("HIS","NE2"),
                                                         ("LYS","NZ") ] )

const _anionic = Set{Tuple{ASCIIString, ASCIIString}}( [ ("ASP","CG"),
                                                        ("ASP","OD1"),
                                                        ("ASP","OD2"),
                                                        ("GLU","CD"),
                                                        ("GLU","OE1"),
                                                        ("GLU","OE2") ] )

_add_CTER_O!(_anionic)

"""Keys come from Table 1 of Bickerton et. al. 2011,
The hydrogen names of the donor comes from: http://biomachina.org/courses/modeling/download/topallh22x.pro
Synonyms come from: http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl"""
const _hbond_donor = Dict{Tuple{ASCIIString, ASCIIString}, Vector{ASCIIString}}( ("ALA","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("ARG","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("ARG","NE") => ["HE", "HNE"],
                                                                                ("ARG","NH1") => ["HH11","HH12", "1HH1","2HH1"],
                                                                                ("ARG","NH2") => ["HH22","HH21", "2HH1", "1HH2"],
                                                                                ("ASN","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("ASN","ND2") => ["HD21","HD22", "HN21", "HN22", "1HD2","2HD2"],
                                                                                ("ASP","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("CYS","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("CYS","SG") => ["HG1", "HG", "HSG"],
                                                                                ("GLN","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("GLN","NE2") => ["HE21","HE22", "1HE2","2HE2", "HN21","HN22"],
                                                                                ("GLU","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("GLY","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("HIS","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("HIS","ND1") => ["HD1"],
                                                                                ("HIS","NE2") => [ "HE2","HNE2"],
                                                                                ("ILE","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("LEU","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("LYS","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("LYS","NZ") => ["HZ1","HZ2","HZ3", "1HZ","2HZ","3HZ", "HNZ1","HNZ2","HNZ3"],
                                                                                ("MET","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("PHE","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("SER","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("SER","OG") => ["HG1", "HG", "HOG"],
                                                                                ("THR","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("THR","OG1") => ["HG1", "HOG1"],
                                                                                ("TRP","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("TRP","NE1") => ["HE1", "HNE1"],
                                                                                ("TYR","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("TYR","OH") => ["HH", "HOH"],
                                                                                ("VAL","N") => ["HN", "H", "HN1","H1","1H", "HN2","H2","2H", "HN3","H3","3H", "HT1","HT2","HT3"],
                                                                                ("PRO","N") => ["HN1","H1","1H", "HN2","H2","2H", "HT1","HT2"] )
# Proline N-Terminal (RESIDUE PROP)

"""
Keys come from Table 1 of Bickerton et. al. 2011,
Antecedents come from come from: http://biomachina.org/courses/modeling/download/topallh22x.pro
Synonyms come from: http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl
"""
const _hbond_acceptor = Dict{Tuple{ASCIIString, ASCIIString}, Vector{ASCIIString}}( ("ALA","O") => ["C"], ("ALA","OT1") => ["C"], ("ALA","OXT") => ["C"], ("ALA","OT2") => ["C"],
                                                                                ("ARG","O") => ["C"], ("ARG","OT1") => ["C"], ("ARG","OXT") => ["C"], ("ARG","OT2") => ["C"],
                                                                                ("ASN","O") => ["C"], ("ASN","OT1") => ["C"], ("ASN","OXT") => ["C"], ("ASN","OT2") => ["C"],
                                                                                ("ASN","OD1") => ["CG"],
                                                                                ("ASP","O") => ["C"], ("ASP","OT1") => ["C"], ("ASP","OXT") => ["C"], ("ASP","OT2") => ["C"],
                                                                                ("ASP","OD1") => ["CG"],
                                                                                ("ASP","OD2") => ["CG"],
                                                                                ("CYS","O") => ["C"], ("CYS","OT1") => ["C"], ("CYS","OXT") => ["C"], ("CYS","OT2") => ["C"],
                                                                                ("CYS","SG") => ["CB"],
                                                                                ("GLN","O") => ["C"], ("GLN","OT1") => ["C"], ("GLN","OXT") => ["C"], ("GLN","OT2") => ["C"],
                                                                                ("GLN","OE1") => ["CD"],
                                                                                ("GLU","O") => ["C"], ("GLU","OT1") => ["C"], ("GLU","OXT") => ["C"], ("GLU","OT2") => ["C"],
                                                                                ("GLU","OE1") => ["CD"],
                                                                                ("GLU","OE2") => ["CD"],
                                                                                ("GLY","O") => ["C"], ("GLY","OT1") => ["C"], ("GLY","OXT") => ["C"], ("GLY","OT2") => ["C"],
                                                                                ("HIS","ND1") => ["CG", "CE1"],
                                                                                ("HIS","NE2") => ["CD2", "CE1"],
                                                                                ("HIS","O") => ["C"], ("HIS","OT1") => ["C"], ("HIS","OXT") => ["C"], ("HIS","OT2") => ["C"],
                                                                                ("ILE","O") => ["C"], ("ILE","OT1") => ["C"], ("ILE","OXT") => ["C"], ("ILE","OT2") => ["C"],
                                                                                ("LEU","O") => ["C"], ("LEU","OT1") => ["C"], ("LEU","OXT") => ["C"], ("LEU","OT2") => ["C"],
                                                                                ("LYS","O") => ["C"], ("LYS","OT1") => ["C"], ("LYS","OXT") => ["C"], ("LYS","OT2") => ["C"],
                                                                                ("MET","O") => ["C"], ("MET","OT1") => ["C"], ("MET","OXT") => ["C"], ("MET","OT2") => ["C"],
                                                                                ("MET","SD") => ["CG", "CE"],
                                                                                ("PHE","O") => ["C"], ("PHE","OT1") => ["C"], ("PHE","OXT") => ["C"], ("PHE","OT2") => ["C"],
                                                                                ("PRO","O") => ["C"], ("PRO","OT1") => ["C"], ("PRO","OXT") => ["C"], ("PRO","OT2") => ["C"],
                                                                                ("SER","O") => ["C"], ("SER","OT1") => ["C"], ("SER","OXT") => ["C"], ("SER","OT2") => ["C"],
                                                                                ("SER","OG") => ["CB"],
                                                                                ("THR","O") => ["C"], ("THR","OT1") => ["C"], ("THR","OXT") => ["C"], ("THR","OT2") => ["C"],
                                                                                ("THR","OG1") => ["CB"],
                                                                                ("TRP","O") => ["C"], ("TRP","OT1") => ["C"], ("TRP","OXT") => ["C"], ("TRP","OT2") => ["C"],
                                                                                ("TYR","O") => ["C"], ("TYR","OT1") => ["C"], ("TYR","OXT") => ["C"], ("TYR","OT2") => ["C"],
                                                                                ("VAL","O") => ["C"], ("VAL","OT1") => ["C"], ("VAL","OXT") => ["C"], ("VAL","OT2") => ["C"] )

function _generate_dict!(dict, input_dict)
    for (res, atom) in keys(input_dict)
        if haskey(dict, res)
            push!(dict[res], atom)
        else
            dict[res] = Set{ASCIIString}(ASCIIString[ atom ])
        end
    end
    dict
end

function _generate_dict!(dict, input_set::Set{Tuple{ASCIIString,ASCIIString}})
    for (res, atom) in input_set
        if haskey(dict, res)
            push!(dict[res], atom)
        else
            dict[res] = Set{ASCIIString}(ASCIIString[ atom ])
        end
    end
    dict
end

function _generate_interaction_keys(vdw, hyd, aro, cat, ani)
    dict = Dict{ASCIIString, Set{ASCIIString}}()
    _generate_dict!(dict, vdw)
    _generate_dict!(dict, hyd)
    _generate_dict!(dict, aro)
    _generate_dict!(dict, cat)
    _generate_dict!(dict, ani)
    dict
end

const _interaction_keys = _generate_interaction_keys(vanderwaalsradius, _hydrophobic, _aromatic, _cationic, _anionic)

_generate_atoms_set(res::PDBResidue) = ASCIIString[ atom.atom for atom in res.atoms[findheavy(res)] ]

"""
This function takes a `PDBResidue` and returns `true` only if all the atoms can be used for checking interactions.
"""
function check_atoms_for_interactions(res::PDBResidue)
    atoms = _generate_atoms_set(res)
    if haskey(_interaction_keys, res.id.name)
        used = _interaction_keys[res.id.name]
    else
        warn(string("RESIDUE: ", res.id.name, " is unknown for MIToS.PDB (AtomsData.jl)"))
        return(false)
    end
    for atom in atoms
        if !( atom in used )
            warn(string("RESIDUE ", res.id.name," ATOM ", atom, " is unknown for MIToS.PDB (AtomsData.jl)"))
            return( false )
        end
    end
    true
end


text:
module ReadSensors

include("../Types.jl")

using Dates

export readsensors

nowepoch = () -> floor(Int, datetime2unix(now())*1000)

active = true

start = () -> active = true
stop = () -> active = false

function init(output:: Channel)
    #while active
    for i in range(1,600; step=1)
        x = parse(Float64, readlines("/sys/bus/iio/devices/iio:device2/in_accel_x_raw")[1])
        y = parse(Float64, readlines("/sys/bus/iio/devices/iio:device2/in_accel_y_raw")[1])
        z = parse(Float64, readlines("/sys/bus/iio/devices/iio:device2/in_accel_z_raw")[1])
        magnitude = x^2 + y^2 + z^2
        dp = Types.Datapoint(nowepoch(), magnitude)
        put!(output, dp)
        sleep(0.1)
    end
end

end

max_stars_repo_name: eliabiraschi/pedometer

max_issues_repo_name: eliabiraschi/pedometer

max_forks_repo_name: eliabiraschi/pedometer


text:
struct DistanceMatrix{T<:Real} <: AbstractArray{T,2}
    dm::AbstractArray{T,2}
    samples::AbstractVector{S} where S
    distance::PreMetric
end

struct PCoA{T<:AbstractFloat} <: AbstractArray{T,2}
    eigenvectors::Array{T,2}
    eigenvalues::Array{T,1}
    variance_explained::Array{T,1}
end

DistanceMatrix(dm::AbstractArray{T,2}, distance) where {T<:Real} = DistanceMatrix(dm, Vector(1:size(dm,1)), distance)

@forward_func DistanceMatrix.dm Base.getindex, Base.setindex, Base.length, Base.size
@forward_func PCoA.eigenvectors Base.getindex, Base.setindex, Base.length, Base.size

function getdm(t::AbstractComMatrix, distance::PreMetric)
    dm = pairwise(distance, t.occurrences)
    for i in eachindex(dm); if isnan(dm[i]); dm[i] = 1; end; end
    return DistanceMatrix(
            dm,
            samplenames(t),
            distance)
end

function getdm(t::AbstractArray, distance::PreMetric)
    dm = pairwise(distance, t)
    for i in eachindex(dm); if isnan(dm[i]); dm[i] = 1; end; end
    return DistanceMatrix(
            dm,
            Vector(1:size(t,2)),
            distance)
end

function getdm(df::DataFrame, distance::PreMetric)
    dm = pairwise(distance, convert(Matrix, df[2:end]))
    for i in eachindex(dm); if isnan(dm[i]); dm[i] = 1; end; end
    return DistanceMatrix(
            dm,
            Vector(names(df[2:end])),
            distance)
end

function getrowdm(abt::AbstractComMatrix, distance::PreMetric)
    m = abt.occurrences'
    return DistanceMatrix(
            pairwise(distance, m),
            samplenames(abt),
            distance)
end

function getrowdm(arr::AbstractArray, distance::PreMetric)
    m = arr'
    return DistanceMatrix(
            pairwise(distance, m),
            Vector(1:size(arr,1)),
            distance)
end

function getrowdm(df::DataFrame, distance::PreMetric)
    m = convert(Matrix, df[2:end])'
    return DistanceMatrix(
            pairwise(distance, m),
            Vector(df[:,1]),
            distance)
end

function pcoa(D::DistanceMatrix; correct_neg::Bool=false)
    n = size(D,1)
    A = -1/2 * D.^2
    Δ1 = getdelta(A)

    f = sortedeig(Δ1)

    if correct_neg && f.values[end] < 0
        c = abs(f.values[end])
        for (i,h) in [(i, h) for i in 1:n for h in 1:n if i != h]
            A[h,i] = -1/2 * D[h,i]^2 - c
        end
        Δ1 = getdelta(A)
        f = sortedeig(Δ1)
    end

    vals = f.values[1:n-1]
    return PCoA(f.vectors[:,1:n-1], vals, [v/sum(vals) for v in vals])
end

function sortedeig(M::Array{Float64,2})
    f = eigen(M, scale=true, permute=true)
    v = real.(f.values)
    p = sortperm(v, rev = true)
    return LinearAlgebra.Eigen(v[p], real.(f.vectors[:,p]))
end

function getdelta(A::AbstractArray{T,2}) where T <: AbstractFloat
    n = size(A,1)
    return reshape(
        [A[h,i] - mean(A[h,:]) - mean(A[:,i]) + mean(A) for i in 1:n for h in 1:n],n,n)
end

@inline eigenvalue(p::PCoA, inds...) = p.eigenvalues[inds...]
@inline variance(p::PCoA, inds...) = p.variance_explained[inds...]
@inline principalcoord(p::PCoA, inds...) = p[:,inds...]

## Diversity Measures

function shannon(v::AbstractVector{T}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return -sum([log(x^x) for x in relab])
end

function ginisimpson(v::AbstractVector{T}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return 1 - sum([x^2 for x in relab])
end


text:
@testset "confident tests" begin
    a = 1.0
    b = 3
    ci1 = ConfidenceInterval(a, b)
    @test ci1.E == 2.0
    @test ci1.min === 1.0
    @test ci1.max === 3.0

    @test_throws ErrorException ConfidenceInterval(b, a)
end

@testset "not so confident tests" begin
    a = 9.7
    b = 153.9

    nc1 = notsoconfident(b, a)
    nc2 = notsoconfident(a, b)

    @test nc1.min == 9.7
    @test nc1.max == 153.9
    @test nc1.E == 81.8

    @test nc1.min == nc2.min
    @test nc1.max == nc2.max
    @test nc1.E == nc2.E
end


text:
using PyCall
using DataFrames
using Pandas: query, read_csv, read_excel, read_pickle
using RCall
using Pipe
#using Winston
using MLDataUtils
using Gadfly
using TSne
using Clustering
using Plots
using StatsBase

plotlyjs()

@pyimport re
@pyimport pandas as pd
@pyimport numpy as np

@pyimport scipy.stats as sps

unshift!(PyVector(pyimport("sys")["path"]),
            "/home/benmo/OneDrive/GitHub/PyPackages/Algos")
@pyimport ModelsML

include("./juliaAlgo.jl")

stocks_df = pd.read_pickle("/home/benmo/Data/PyObjects/stocks_df.gz",
                            compression="gzip")

stocks_df = pd.pivot_table(stocks_df, values=collect(stocks_df[:columns])[1:end-1],
                       index=stocks_df[:index], columns="Ticker")

stocksPrice = stocks_df[:Adj_Close]
stocksVol = stocks_df[:Volume]

sp500 = pd.read_pickle("/home/benmo/Data/PyObjects/SP500.pkl")
sp500long = pd.read_pickle("/home/benmo/Data/PyObjects/SP500long.pkl")

sp500[:DATE] = pd.to_datetime(sp500[:DATE])
sp500 = sp500[:set_index]("DATE")
sp500[:SP500] = sp500[:SP500][:apply](pd.to_numeric, errors="coerce")

stocksRtn = stocksPrice[:apply](np.log)[:diff]()
sp500 = sp500[:apply](np.log)[:diff]()

sp500μ = sp500[:mean]()
sp500σ = sp500[:std](ddof=0)

stocksRtnσ = stocksRtn[:std](ddof=0)

βⱼ = pd.read_pickle("/home/benmo/Data/PyObjects/stocks_betaj1.pkl")
βⱼ = py"$βⱼ.drop(columns='T').drop('T')"
#βⱼ = βⱼ[:dropna](axis=1,how="all")[:dropna](how="all")

stocksRtn, stocksVol, sVol_pd, mVol_pd, lVol_pd = clean_Vol(
                                                    βⱼ, stocksRtn, stocksVol;
                                                    startY=2013
                                                    )

stockcols = py"list(filter(lambda x: x in $βⱼ.columns.tolist(), $stocksRtn.columns.tolist()))"

βⱼ = py"$βⱼ[$stockcols].loc[$stockcols]" |> pandas2julia

#stocksRtn = stocksRtn[:join](sp500[:apply](
        #np.log)[:diff](), how="inner")
#stocksRtn[:SP500] = stocksRtn[:SP500][:fillna](method="ffill")
#stocksRtn[:SP500] = stocksRtn[:SP500][:fillna](method="bfill")

Rσ = stocksRtn[:std](ddof=0)

for stock in stocksRtn[:columns] stocksRtn[stock] =
    stocksRtn[stock] ./ Rσ[stock] end

stockCats, dims, km1 = classifyRⱼ(βⱼ; perplx=50, nclust=30)
Gadfly.plot(x=dims[:,1], y=dims[:,2], color = km1.assignments, Geom.point)

tech_idx = stockCats[stockCats[:,2] .== "AAPL",1][1]
namesTech = stockCats[stockCats[:,1] .== tech_idx,2]
techVol = py"$stocksVol[$namesTech]['2013':]"
techPrice = py"$stocksPrice[$namesTech]['2013':]"

techRtn = techPrice[:apply](np.log)[:diff]()[:fillna](0)

techRtnσ, techRtnμ = techRtn[:std](ddof=0), techRtn[:mean]()
techRtn_max, techRtn_min = techRtn[:max](), techRtn[:min]()

scaler = ModelsML.scale(py"$techRtn")

techRtn = scaler[:tfm](py"$techRtn")

techRtn_bxcx = techPrice[:apply](np.log)[:diff]()[:fillna](0)
scaler_bxcx = ModelsML.scale(py"$techRtn_bxcx", how="boxcox")
techRtn_bxcx = scaler_bxcx[:tfm](py"$techRtn_bxcx")

techRtnJ = techRtn |> pandas2julia
techRtnJ_bxcx = techRtn_bxcx |> pandas2julia

dates = techRtn[:index][:tolist]()

#histogram(techRtnJ[3], opacity=0.45)
#histogram!(techRtnJ_bxcx[3], opacity=0.45)

stockᵢ = try sample(filter((x) -> x ∈ sVol_pd[:columns][:tolist](), namesTech))
            catch sample(filter((x) -> x ∈ mVol_pd[:columns][:tolist](), namesTech)) end

stocksⱼ = sampleStocks(namesTech[namesTech .!= stockᵢ], sVol_pd, mVol_pd, lVol_pd, 0.1, 0.4, 0.5; N=10)

#mdl = ModelsML.ts_LSTM(techRtn, stockᵢ)
#ŷ, y_test = py"$mdl.yhat.astype('float64'), $mdl.y_test.astype('float64')"

#ŷ , y_test = (ŷ, y_test) .|> ((x) -> x * techRtnσ[stockᵢ][:tolist]() +
                                                        #techRtnμ[stockᵢ][:tolist]()) .|>
                                                        #transpose .|> DataFrame

yᵣ = DataFrame([])
yᵣ[:time] = 1:4

yᵣ = hcat(yᵣ, ŷ[1], y_test[1]; makeunique=true)

Gadfly.plot( ŷ, x=:time, y=Col.value(:x1,:x1_1),
            color=Col.index(:x1,:x1_1), Geom.line)

y₁ = DataFrame([])
y₁[:time] = 1:251
y₁ = hcat(y₁, reduce(vcat, DataFrames.columns( ŷ[1,:])),
                    reduce(vcat, DataFrames.columns( y_test[1,:])))

Gadfly.plot( y₁, x=:time, y=Col.value(:x1,:x1_1),
            color=Col.index(:x1,:x1_1), Geom.line)

max_stars_repo_name: benzomo/PyPackages

max_issues_repo_name: benzomo/PyPackages

max_forks_repo_name: benzomo/PyPackages


text:
using Serialization
using Test

using ComparativeAutograder
using ComparativeAutograder: RunnerInput, runner_main

@testset "Test runner_main" begin
    func_body = """
        function add(a, b)
            a + b
        end
        """
    suite = TestSuite([
        TestCase(1, 2),
        TestCase("foo", "bar"),
    ])
    subfile = tempname()
    outfile = tempname()
    runner_input = RunnerInput(
        suite,
        subfile,
        :add,
        outfile,
    )

    routput, foutput = nothing, nothing
    try
        open(io -> write(io, func_body), subfile, "w")
        runner_output = runner_main(runner_input; exit=false, write=true)

        # The result is returned but it is also written to file.
        # We will check both.
        routput = runner_output
        foutput = open(io -> deserialize(io), outfile)
    finally
        # pass
        # rm.((subfile, outfile); force=true)
    end

    @test routput !== nothing
    @test foutput !== nothing

    for output in (routput, foutput)
        cases = output.result.results
        @test cases[1].result == 3
        @test occursin("MethodError", cases[2].exception)
    end
end

@testset "Test subprocess runner invalid invokation" begin
    @testset "Test subprocess runner with invalid usage" begin
        @test_throws Exception run(`$(Base.julia_cmd()) -e $(ComparativeAutograder.RUNNER_CMD)`)
    end

    @testset "Test subprocess runner with nonexistant input file" begin
        @test_throws Exception run(`$(Base.julia_cmd()) -e $(ComparativeAutograder.RUNNER_CMD) $(tempname())`)
    end
end


text:
# # 2D Brusselator
#
#
#
# This example solves the 2D Brusselator.
#
# ```math
# \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + E - (B+1)u + v u^2, \\
# \frac{\partial v}{\partial t} = \frac{\partial^2 v}{\partial x^2} + B u - v u^2.
# ```
#

module Brusselator

using Printf,
      Random,
      Reexport

@reexport using FourierFlows

using LinearAlgebra: mul!, ldiv!

export
  set_uv!,
  updatevars!,
  get_righthandside

# ## Coding up the equations
# ### A demonstration of FourierFlows.jl framework
#
# What follows is a step-by-step tutorial showing how you can create your own
# solver for an equation of your liking.

# The basic building blocks for a `FourierFlows.Problem()` are:
# - `Grid` struct containining the physical and wavenumber grid for the problem,
# - `Params` struct containining all the parameters of the problem,
# - `Vars` struct containining arrays with the variables used in the problem,
# - `Equation` struct containining the coefficients of the linear operator $\mathcal{L}$ and the function that computes the nonlinear terms, usually named `calcN!()`.
#
# The `Grid` structure is provided by FourierFlows.jl. One simply has to call one of
# the `OneDGrid()`,  `TwoDGrid()`, or `ThreeDGrid()` grid constructors, depending
# on the dimensionality of the problem. All other structs mentioned above are problem-specific
# and need to be constructed for every set of equations we want to solve.

# First lets construct the `Params` struct that contains all parameters of the problem

mutable struct Params{T} <: AbstractParams
   B :: T         # B parameter
   D :: T         # D parameter
   E :: T         # E parameter
end

# Now the `Vars` struct that contains all variables used in this problem. For this
# problem `Vars` includes the represenations of the flow fields in physical space
# `u` and `v` and their Fourier transforms `uh` and `vh`.

struct Vars{Aphys, Atrans} <: AbstractVars
     u :: Aphys
     v :: Aphys
   vu² :: Aphys
    uh :: Atrans
    vh :: Atrans
  vu²h :: Atrans
end

# A constructor populates empty arrays based on the dimension of the `grid`
# and then creates `Vars` struct.
"""
    Vars!(dev, grid)
Constructs Vars based on the dimensions of arrays of the `grid`.
"""
function Vars(::Dev, grid::AbstractGrid) where Dev
  T = eltype(grid)
  @devzeros Dev T grid.nx u v vu²
  @devzeros Dev Complex{T} grid.nkr uh vh vu²h
  return Vars(u, v, vu², uh, vh, vu²h)
end

# In Fourier space, the 1D linear shallow water dynamics are:
#
# ```math
# \frac{\partial \hat{u}}{\partial t} = \underbrace{- D |\boldsymbol{k}|^2 -(B+1) }_{\mathcal{L}_u} \hat{u} + \underbrace{ E + \widehat{v u^2} }_{\mathcal{N}_u} \; , \\
# \frac{\partial \hat{v}}{\partial t} = \underbrace{- D |\boldsymbol{k}|^2 }_{\mathcal{L}_v} \hat{v} } + \underbrace{ B \hat{u} - \widehat{v u^2} }_{\mathcal{N}_v} \; . \\
# ```
#
# With these in mind, we construct function `calcN!` that computes the nonlinear terms.
#
"""
    calcN!(N, sol, t, clock, vars, params, grid)
The function that computes the nonlinear terms for our problem.
"""
function calcN!(N, sol, t, clock, vars, params, grid)
  @. vars.uh = sol[:, 1]
  @. vars.vh = sol[:, 2]

  ldiv!(vars.u, grid.rfftplan, vars.uh)
  ldiv!(vars.v, grid.rfftplan, vars.vh)

  @. vars.vu² = vars.v * vars.u^2

  mul!(vars.vu²h, grid.rfftplan, vars.vu²)

  @views @. N[:, 1] = + vars.vu²h - (params.B + 1) * sol[:, 1]    # + v*u² - (B+1) u
  N[1, 1] += params.E * grid.nx # note that fft(constant) = constant * nx

  @views @. N[:, 2] = - vars.vu²h + params.B * sol[:, 1]          # - v*u² + B*u

  return nothing
end

# Next we construct the `Equation` struct:

"""
    Equation!(prob)
Construct the equation: the linear part and the nonlinear part, which is
computed by `caclN!` function.
"""
function Equation(dev, params, grid::AbstractGrid)
  T = eltype(grid)
  L = zeros(dev, T, (grid.nkr, 2))

  diffusion = @. - grid.kr^2

  @. L[:, 1] = params.D * diffusion  # D*∂²u/∂x²
  @. L[:, 2] = diffusion             #   ∂²v/∂x²

  return FourierFlows.Equation(L, calcN!, grid)
end

# We now have all necessary building blocks to construct a `FourierFlows.Problem`.
# It would be useful, however, to define some more "helper functions". For example,
# a function that updates all variables given the solution `sol` which comprises $\hat{u}$
# and $\hat{v}$:

"""
    updatevars!(prob)
Update the variables in `prob.vars` using the solution in `prob.sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  @. vars.uh = sol[:, 1]
  @. vars.vh = sol[:, 2]

  ldiv!(vars.u, grid.rfftplan, deepcopy(sol[:, 1])) # use deepcopy() because irfft destroys its input
  ldiv!(vars.v, grid.rfftplan, deepcopy(sol[:, 2])) # use deepcopy() because irfft destroys its input
  return nothing
end

# Another useful function is one that prescribes an initial condition to the state variable `sol`.

"""
    set_uv!(prob, u0, v0)
Sets the state variable `prob.sol` as the Fourier transforms of `u0` and `v0`
and update all variables in `prob.vars`.
"""
function set_uv!(prob, u0, v0)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  @. vars.u = u0
  @. vars.v = v0

  mul!(vars.uh, grid.rfftplan, vars.u)
  mul!(vars.vh, grid.rfftplan, vars.v)

  @. sol[:, 1] = vars.uh
  @. sol[:, 2] = vars.vh

  updatevars!(prob)
  return nothing
end

"""
    get_righthandside!(prob)
Returns the right-hand-side of the equation. In this case, the function returns
a 2-column array with the time-tendencies for uh and vh.
"""
function get_righthandside(prob)

  vars, params, grid, sol, clock = prob.vars, prob.params, prob.grid, prob.sol, prob.clock

  L = prob.eqn.L

  N = zeros(eltype(sol), size(sol))

  prob.eqn.calcN!(N, sol, clock.t, clock, vars, params, grid)

  return @. L*sol + N
end

"""
    timestep_up_to_t!(prob, t)
Timesteps the problem u the right-hand-side of the equation. In this case, the function returns
a 2-column array with the time-tendencies for uh and vh.
"""
function timestep_up_to_t(prob, t_final)
  dt = prob.clock.dt

  t_initial = prob.clock.t

  t_final = t

  nsteps = Int(floor((t_final - t_initial) / dt))

  if nsteps < 0; error("t_final < prob.t"); end

  stepforward!(prob, nsteps)

  t_remaining = t_final - t_initial - nsteps*dt

  prob.clock.dt = t_remaining

  stepforward!(prob)

  prob.clock.dt = dt
end

end # module