text:
# The code in this file is a small port from scikit-learn's and numpy's
# library which is distributed under the 3-Clause BSD license.
# The rest of CausalForest.jl is released under the MIT license.

module treecausation
    include("../util.jl")

    import Random
    import PoissonRandom # to generate random mtry
    import StatsBase # to sample features
    export fit

    mutable struct NodeMeta{S}
        l           :: NodeMeta{S}  # right child
        r           :: NodeMeta{S}  # left child
        old_purity  :: Float64
        new_purity  :: Float64
        feature     :: Int          # feature used for splitting
        threshold   :: S            # threshold value
        is_leaf     :: Bool
        depth       :: Int
        region      :: UnitRange{Int} # a slice of the samples used to decide the split of the node
        features    :: Vector{Int}    # a list of features
        split_at    :: Int            # index of samples
        function NodeMeta{S}(
                features :: Vector{Int},
                region   :: UnitRange{Int},
                depth    :: Int,
                pur      :: Float64) where S
            node = new{S}()
            node.depth = depth
            node.old_purity = pur
            node.region = region
            node.features = features
            node.is_leaf = false
            node
        end
    end

    struct Tree{S}
        root   :: NodeMeta{S}
        inds   :: Vector{Int}
    end

    # find an optimal split that satisfy the given constraints
    # (max_depth, min_samples_split)
    function _split!(
            X                   :: AbstractMatrix{S}, # the feature array
            Y                   :: AbstractVector{Float64}, # the label array
            W                   :: AbstractVector{Int}, # the treatment array
            node                :: NodeMeta{S}, # the node to split
            const_mtry          :: Bool,
            m_pois              :: Int, # hyperparameter for poisson for random mtry
            max_depth           :: Int, # the maximum depth of the resultant tree
            min_samples_leaf    :: Int, # the minimum number of samples each leaf needs to have
            min_samples_split   :: Int, # the minimum number of samples in needed for a split
            indX                :: AbstractVector{Int}, # an array of sample indices,
                                                # we split using samples in indX[node.region]
            # the three arrays below are given for optimization purposes
            Xf                  :: AbstractVector{S},
            Yf                  :: AbstractVector{Float64},
            Wf                  :: AbstractVector{Int},
            rng                 :: Random.AbstractRNG) where {S}

        region = node.region
        n_samples = length(region)
        r_start = region.start - 1

        @inbounds @simd for i in 1:n_samples
            Yf[i] = Y[indX[i + r_start]]
            Wf[i] = W[indX[i + r_start]]
        end

        t_w_sum = zero(Int) # sum of W in node
        t_y_sum = zero(Float64) # sum of Y in node
        t_w_ssq = zero(Int) # sum of W squares in node
        t_wy_sum = zero(Float64) # sum of WY in node
        @inbounds @simd for i in 1:n_samples
            t_w_sum += Wf[i]
            t_y_sum += Yf[i]
            t_w_ssq += Wf[i]*Wf[i]
            t_wy_sum += Wf[i]*Yf[i]
        end

        if (min_samples_leaf * 2 >  n_samples
         || min_samples_split    >  n_samples
         || max_depth            <= node.depth
         )
            node.is_leaf = true
            return
        end

        features = node.features
        n_features = length(features)
        best_purity = typemin(Float64)
        best_feature = -1
        threshold_lo = X[1]
        threshold_hi = X[1]

        # true if every feature is constant
        unsplittable = true
        total_features = size(X, 2)
        if const_mtry
            mtry = m_pois
        else
            mtry = min(max(PoissonRandom.pois_rand(m_pois), 1), total_features)
        end
        random_features = StatsBase.sample(1:total_features, mtry, replace=false)
        for feature in random_features

            r_w_sum = t_w_sum
            l_w_sum = zero(Int)
            r_y_sum = t_y_sum
            l_y_sum = zero(Float64)
            r_w_ssq = t_w_ssq
            l_w_ssq = zero(Int)
            r_wy_sum = t_wy_sum
            l_wy_sum = zero(Float64)

            @simd for i in 1:n_samples
                Xf[i] = X[indX[i + r_start], feature]
            end

            # sort Yf and indX by Xf
            util.q_bi_sort!(Xf, indX, 1, n_samples, r_start)
            @simd for i in 1:n_samples
                Yf[i] = Y[indX[i + r_start]]
                Wf[i] = W[indX[i + r_start]]
            end
            nl, nr = 0, n_samples
            # lo and hi are the indices of
            # the least upper bound and the greatest lower bound
            # of the left and right nodes respectively
            hi = 0
            last_f = Xf[1]
            while hi < n_samples
                lo = hi + 1
                curr_f = Xf[lo]
                hi = (lo < n_samples && curr_f == Xf[lo+1]
                    ? searchsortedlast(Xf, curr_f, lo, n_samples, Base.Order.Forward)
                    : lo)

                if lo-1 >= min_samples_leaf && n_samples - (lo-1) >= min_samples_leaf
                    unsplittable = false
                    difference = (l_w_ssq/nl - (l_w_sum/nl)^2)/(l_wy_sum/nl - (l_w_sum/nl))-(r_w_ssq/nr - (r_w_sum/nr)^2)/(r_wy_sum/nr - (r_w_sum/nr))
                    purity = (nl*nr)/(n_samples*n_samples)*difference*difference
                    if purity > best_purity && !isapprox(purity, best_purity)
                        # will take average at the end, if possible
                        threshold_lo = last_f
                        threshold_hi = curr_f
                        best_purity  = purity
                        best_feature = feature
                    end
                end

                # update values
                # that would require the smaller number of iterations
                if (hi << 1) < n_samples + lo # i.e., hi - lo < n_samples - hi
                    @simd for i in lo:hi
                        nr   -= 1
                        r_w_sum -= Wf[i]
                        r_y_sum -= Yf[i]
                        r_w_ssq -= Wf[i]*Wf[i]
                        r_wy_sum -= Wf[i]*Yf[i]
                    end
                else
                    nr = r_w_sum = r_w_ssq = 0
                    r_y_sum = r_wy_sum = zero(Float64)
                    @simd for i in (hi+1):n_samples
                        nr   += 1
                        r_w_sum += Wf[i]
                        r_y_sum += Yf[i]
                        r_w_ssq += Wf[i]*Wf[i]
                        r_wy_sum += Wf[i]*Yf[i]
                    end
                end
                nl   = n_samples - nr
                l_w_sum = t_w_sum - r_w_sum
                l_y_sum = t_y_sum - r_y_sum
                l_w_ssq = t_w_ssq - r_w_ssq
                l_wy_sum = t_wy_sum - r_wy_sum

                last_f = curr_f
            end

        end

        # no splits honor min_samples_leaf
        node.new_purity = best_purity
        @inbounds if (unsplittable
                || best_feature == -1) # in case best pur = inf
            node.is_leaf = true
            return
        else

            @simd for i in 1:n_samples
                Xf[i] = X[indX[i + r_start], best_feature]
            end

            try
                node.threshold = (threshold_lo + threshold_hi) / 2.0
            catch
                node.threshold = threshold_hi
            end
            # split the samples into two parts: ones that are greater than
            # the threshold and ones that are less than or equal to the threshold
            #                                 ---------------------
            # (so we partition at threshold_lo instead of node.threshold)
            node.split_at = util.partition!(indX, Xf, threshold_lo, region)
            node.feature = best_feature
            node.features = features
        end

    end

    @inline function fork!(node :: NodeMeta{S}) where S
        ind = node.split_at
        region = node.region
        features = node.features
        # no need to copy because we will copy at the end
        node.l = NodeMeta{S}(features, region[    1:ind], node.depth + 1, node.new_purity)
        node.r = NodeMeta{S}(features, region[ind+1:end], node.depth + 1, node.new_purity)
    end

    function _fit(
            X                     :: AbstractMatrix{S},
            Y                     :: AbstractVector{Float64},
            W                     :: AbstractVector{Int},
            indX                  :: AbstractVector{Int},
            const_mtry            :: Bool,
            m_pois                :: Int,
            max_depth             :: Int,
            min_samples_leaf      :: Int,
            min_samples_split     :: Int,
            rng=Random.GLOBAL_RNG :: Random.AbstractRNG) where {S}

        n_samples, n_features = size(X)

        Yf  = Array{Float64}(undef, n_samples)
        Xf  = Array{S}(undef, n_samples)
        Wf  = Array{Int}(undef, n_samples)

        root = NodeMeta{S}(collect(1:n_features), 1:length(indX), 0, Inf)
        stack = NodeMeta{S}[root]

        @inbounds while length(stack) > 0
            node = pop!(stack)
            _split!(
                X, Y, W,
                node,
                const_mtry,
                m_pois,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                indX,
                Xf, Yf, Wf,
                rng)
            if !node.is_leaf
                fork!(node)
                push!(stack, node.r)
                push!(stack, node.l)
            end
        end
        return (root, indX)
    end

    function fit(;
            X                     :: AbstractMatrix{S},
            Y                     :: AbstractVector{Float64},
            W                     :: AbstractVector{Int},
            indX                  :: AbstractVector{Int},
            const_mtry            :: Bool,
            m_pois                :: Int,
            max_depth             :: Int,
            min_samples_leaf      :: Int,
            min_samples_split     :: Int,
            rng=Random.GLOBAL_RNG :: Random.AbstractRNG) where {S}

        n_samples, n_features = size(X)

        util.check_input(
            X,
            Y,
            W,
            indX,
            const_mtry,
            m_pois,
            max_depth,
            min_samples_leaf,
            min_samples_split)

        root, indX = _fit(
            X,
            Y,
            W,
            indX,
            const_mtry,
            m_pois,
            max_depth,
            min_samples_leaf,
            min_samples_split,
            rng)

        return Tree{S}(root, indX)

    end
end


text:
using GeometricIntegrators

set_config(:tab_compensated_summation, false)

const out_dir = "../../docs/src/integrators/runge_kutta/guiding_center_4d/"
const run_id  = "fast_deeply_passing_firk"

include("guiding_center_4d_fast_deeply_passing.jl")
include("../guiding_center_4d_firk.jl")

ext: jl


text:
module MagnitudeDistributions

using StatsBase, Statistics, Distributions, Random
using HypothesisTests # MannWhitneyUTest in mbass

export histmag # utils
export maxcurv, mbass, mbs # Magnitude of completeness
export GRclassic, fit_mle, fit_binned, fit_stderr, fit_stderrbinned # distributions
export bvalue, avalue

# utils
include("./utils/binning.jl")

# magnitude of completeness
include("./completeness/maxc.jl")
include("./completeness/mbass.jl")
include("./completeness/mbs.jl")

# frequency magnitude distributions
include("./distributions/GRclassic.jl")
include("./distributions/univariates.jl")

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_colorpicker

"""
    daq_colorpicker(;kwargs...)

A ColorPicker component.
A color picker.
Keyword arguments:
- `id` (String; optional): The ID used to identify the color picker in Dash callbacks
- `className` (String; optional): Class to apply to the root component element
- `disabled` (Bool; optional): If true, color cannot be picked.
- `label` (optional): Description to be displayed alongside the control. To control styling,
pass an object with label and style properties. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the indicator label is positioned
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
- `size` (Real; optional): Size (width) of the component in pixels
- `style` (Dict; optional): Style to apply to the root component element
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
- `value` (optional): Color value of the picker. value has the following type: lists containing elements 'hex', 'rbg'.
Those elements have the following types:
  - `hex` (String; optional): Hex string
  - `rbg` (optional): RGB/RGBA object. rbg has the following type: lists containing elements 'r', 'g', 'b', 'a'.
Those elements have the following types:
  - `r` (Real; optional)
  - `g` (Real; optional)
  - `b` (Real; optional)
  - `a` (Real; optional)
"""
function daq_colorpicker(; kwargs...)
        available_props = Symbol[:id, :className, :disabled, :label, :labelPosition, :persisted_props, :persistence, :persistence_type, :size, :style, :theme, :value]
        wild_props = Symbol[]
        return Component("daq_colorpicker", "ColorPicker", "dash_daq", available_props, wild_props; kwargs...)
end


text:
# Framework for
# Edward J. Xu, edxu96@outlook.com
# Aug 25, 2019

module MatrixOptimFramework

export

"MILP model for Benders decomposition in matrix form."
mutable struct Mod

    solution::Union{Solution, Missing}

    function Mod()
        new()
    end
end

mutable struct Mod()

    function Mod()
        model_ray = Model(with_optimizer(GLPK.Optimizer))
        @variable(model_ray, vec_u[1: n_cons] >= 0)
        @objective(model_ray, Max, 1)
        @constraint(model_ray, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1] == 1)
        @constraint(model_ray, vec_cons, transpose(mat_a) * vec_u .<= 0)
        new(expr=model)
    end
end

"Head of Linked List of Sub Model"
struct ModSubHead()
    expr
    ref

    function ModSubHead(n_cons, vec_yBar, vec_b, mat_b, mat_a, vec_c)
        expr = Model(with_optimizer(GLPK.Optimizer))
        @variable(expr, vec_u[1: n_cons] >= 0)
        @objective(expr, Max, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1])
        @constraint(expr, vec_cons, transpose(mat_a) * vec_u .<= vec_c)
        new(expr, ref=ModSub())
    end
end

function get_mod_sub(mod_sub_head::ModSubHead)
    return mod_sub
end

mutable struct ModSub()
    bool
    obj
    vec_uBar
    vec_result_x

    ref

    function ModSub()
        new(missing, missing, missing, missing, ModSub())
    end
end

function solve_mod_sub!(mod::ModSub)
    mod.st_solution = optimize!(mod.expr)
    mod.vec_uBar = value_vec(vec_u)

    mod.bool = true
    mod.vec_result_x = zeros(length(vec_c))
    mod.obj = NaN
    if Int(primal_status(model)) == 1
        mod.vec_result_x = dual_vec(vec_cons)
        mod.obj = objective_value(model)
    elseif Int(primal_status(model)) == 4  # Unbounded
        println("    Not solved to optimality because feasible set is unbounded.")
        mod.bool = false
        mod.obj = objective_value(model)
        mod.vec_result_x = repeat([NaN], length(vec_c))
    else  # if Int(primal_status(model)) == 3  # Infeasible
        println(
            "    Not solved to optimality because infeasibility. Something is wrong. $(Int(primal_status(model)))"
        )
        mod.bool = false
        mod.vec_result_x = hcat(repeat([NaN], length(vec_c)))
    end
end

mutable struct ModRay()
    vec_uBar
    obj

    ref

    function ModRay()
        new(vec_uBar=value_vec(vec_u), obj_ray=objective_value(model_ray))
    end
end

function solve_mod_ray!(mod::ModRay)

    optimize!(model_ray)

    return (obj_ray, vec_uBar)
end

function check_whe_continue(boundUp, boundLow, epsilon, result_q, obj_sub,
        timesIteration, timesIterationMax)
    whe_continue = true
    if (boundUp - boundLow <= epsilon) & (boundUp - boundLow >= - epsilon)
        if (result_q - obj_sub <= epsilon) & (result_q - obj_sub >= - epsilon)
            whe_continue = false
        end
    end

    return (timesIteration <= timesIterationMax) && whe_continue
end

"""
Generic Benders Decomposition for Mixed Integer Linear Programming
"""
function solveModMixBD!(mod::ModBD, epsilon=1e-6, timesIterationMax=100)
    if mod.solution is not missing
        println("The model has been solved.")
    else
        println("The model has been solved.")
        mod_mas = ModMas(mod.n_y, mod.vec_min_y, mod.vec_max_y)
        mod_sub_head = ModSub()

        boundUp = Inf
        boundLow = - Inf
        # vec_uBar = zeros(n_cons, 1)  # initial value of master variables
        vec_yBar = zeros(mod.n_y, 1)
        vec_result_x = zeros(mod.n_x, 1)

        obj_sub = 0
        result_q = 0
        timesIteration = 1
        # Must make sure "result_q == obj_sub" in the final iteration
        # while ((boundUp - boundLow > epsilon) && (timesIteration <= timesIterationMax))  !!!
        while check_whe_continue(
                boundUp, boundLow, epsilon, result_q, obj_sub,
                timesIteration, timesIterationMax
            )
            mod_sub = get_mod_sub(mod_mod_sub_head)
            solve_mod_sub!(mod_sub)
            if mod_sub.bool
                println("obj_sub = $(obj_sub) ;")
                boundUp = min(boundUp, mod_sub.obj + (transpose(mod.vec_f) * vec_yBar)[1])
            else
                mod_ray = ModRay()
                solve_mod_ray!(mod_ray)
            end

            solve_mod_mas!(
                mod_mas, q, vec_y, mod.vec_b, mod.mat_b
                )
            vec_yBar = value_vec(mod_mas.vec_y)
            boundLow = max(boundLow, mod_mas.obj)

            timesIteration += 1
        end
    end
    println("End")
end

end


text:
using HDF5

struct ALPSComplex end

function read(g::HDF5Dataset, ::Type{ALPSComplex})
  data = read(g)

  # flip all dimensions since data is stored as row-major
  data = permutedims(data, reverse(1:ndims(data)))
  sz = size(data)

  # interpret last dimension as real/complex
  if exists(attrs(g), "__complex__")
    @assert sz[end] == 2
    R = CartesianIndices(sz[1:end-1])
    data = data[R,1] + 1.0im * data[R, 2]
  end

  return data
end

function write_alpscomplex(parent::Union{HDF5File, HDF5Group}, name::String, data::AbstractArray{Complex{T}}) where T <: HDF5.HDF5Scalar
  # flip all dimensions since data is stored as row-major
  data = permutedims(data, reverse(1:ndims(data)))

  # reinterpret real/complex as new axis
  data = reshape(reinterpret(T, data), 2, size(data)...)

  parent[name] = data
  attrs(parent[name])["__complex__"] = Int8(1)
end

function read(g::HDF5Group, ::Type{Contour})
  nb = read(g, "size")
  branches = Branch[]
  for i in 0:nb-1
    b = BranchEnum(read(g, "branch$i/type") + 1)
    l = read(g, "branch$i/len")
    push!(branches, Branch(b,l))
  end
  return Contour(branches)
end

function write(parent::Union{HDF5File, HDF5Group}, name::String, c::Contour)
  g = g_create(parent, name)
  g["size"] = UInt(nbranches(c))
  for (i, branch) in enumerate(c.branches)
    g["branch$(i-1)/type"] = Int32(Int(branch.domain) - 1)
    g["branch$(i-1)/len"] = length(branch)
  end
end

function read(g::HDF5Group, ::Type{TimeGrid})
  c = read(g["contour"], Contour)
  branch_enums = read(g["branch_enums"])
  values = read(g["values"], ALPSComplex)

  npts_fwd = count(branch_enums .== 0)
  npts_back = count(branch_enums .== 1)
  npts_imag = count(branch_enums .== 2)

  @assert npts_fwd == npts_back

  npts_real = npts_fwd
  grid = TimeGrid(c, npts_real=npts_real, npts_imag=npts_imag)

  @assert map(t -> t.val.val, grid) ≈ values

  return grid
end

function write(parent::Union{HDF5File, HDF5Group}, name::String, grid::TimeGrid)
  g = g_create(parent, name)

  write(g, "contour", grid.contour)

  values = map(t -> t.val.val, grid)
  refs = map(t -> t.val.ref, grid)
  branch_enums = map(t -> Int32(Int(t.val.domain) - 1), grid)

  write_alpscomplex(g, "values", values)
  write(g, "refs", refs)
  write(g, "branch_enums", branch_enums)
end

function read(g::HDF5Group, ::Type{TimeGF})
  data = read(g["data"], ALPSComplex)
  grid = read(g["mesh/1"], TimeGrid)
  return TimeGF(data, grid)
end

function write(parent::Union{HDF5File, HDF5Group}, name::String, gf::TimeGF)
  g = g_create(parent, name)
  write_alpscomplex(g, "data", gf.data)
  write(g, "mesh/1", gf.grid)
  write(g, "mesh/2", gf.grid)
  write(g, "mesh/N", 2)
end


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

using BAT
using Test

using Random, Statistics, LinearAlgebra
using Distributions, PDMats
using StableRNGs
using FillArrays
using ForwardDiff

@testset "standard_normal" begin
    stblrng() = StableRNG(789990641)

    @testset "StandardUvUniform" begin
        @test @inferred(BAT.StandardUvUniform()) isa BAT.StandardUvUniform{Float64}
        @test @inferred(BAT.StandardUvUniform{Float32}()) isa BAT.StandardUvUniform{Float32}
        @test @inferred(BAT.StandardUvUniform{Float64}()) isa BAT.StandardUvUniform{Float64}

        @test @inferred(Uniform(BAT.StandardUvUniform{Float32}())) isa Uniform{Float32}
        @test @inferred(Uniform(BAT.StandardUvUniform{Float64}())) isa Uniform{Float64}
        @test @inferred(Uniform(BAT.StandardUvUniform())) == Uniform()
        @test @inferred(convert(Uniform, BAT.StandardUvUniform())) == Uniform()

        d = BAT.StandardUvUniform()
        dref = Uniform()

        @test @inferred(minimum(d)) == minimum(dref)
        @test @inferred(maximum(d)) == maximum(dref)

        @test @inferred(params(d)) == params(dref)
        @test @inferred(partype(d)) == partype(dref)

        @test @inferred(location(d)) == location(dref)
        @test @inferred(scale(d)) == scale(dref)

        @test @inferred(eltype(typeof(d))) == eltype(typeof(dref))
        @test @inferred(eltype(d)) == eltype(dref)

        @test @inferred(length(d)) == length(dref)
        @test @inferred(size(d)) == size(dref)

        @test @inferred(mean(d)) == mean(dref)
        @test @inferred(median(d)) == median(dref)
        @test @inferred(mode(d)) == mode(dref)
        @test @inferred(modes(d)) ≈ modes(dref)

        @test @inferred(var(d)) == var(dref)
        @test @inferred(std(d)) == std(dref)
        @test @inferred(skewness(d)) == skewness(dref)
        @test @inferred(kurtosis(d)) == kurtosis(dref)

        @test @inferred(entropy(d)) == entropy(dref)

        for x in [-0.5, 0.0, 0.25, 0.75, 1.0, 1.5]
            @test @inferred(logpdf(d, x)) == logpdf(dref, x)
            @test @inferred(pdf(d, x)) == pdf(dref, x)
            @test @inferred(logcdf(d, x)) == logcdf(dref, x)
            @test @inferred(cdf(d, x)) == cdf(dref, x)
            @test @inferred(logccdf(d, x)) == logccdf(dref, x)
            @test @inferred(ccdf(d, x)) == ccdf(dref, x)
        end

        for p in [0.0, 0.25, 0.75, 1.0]
            @test @inferred(quantile(d, p)) == quantile(dref, p)
            @test @inferred(cquantile(d, p)) == cquantile(dref, p)
        end

        for t in [-3, 0, 3]
            @test @inferred(mgf(d, t)) == mgf(dref, t)
            @test @inferred(cf(d, t)) == cf(dref, t)
        end

        @test @inferred(rand(stblrng(), d)) == rand(stblrng(), d)
        @test @inferred(rand(stblrng(), d, 5)) == rand(stblrng(), d, 5)

        @test @inferred(truncated(BAT.StandardUvUniform{Float32}(), -0.5f0, 0.7f0)) isa Uniform{Float32}
        @test truncated(BAT.StandardUvUniform{Float32}(), -0.5f0, 0.7f0) == Uniform(0.0f0, 0.7f0)
        @test truncated(BAT.StandardUvUniform{Float32}(), 0.2f0, 0.7f0) == Uniform(0.2f0, 0.7f0)

        @test @inferred(product_distribution(fill(BAT.StandardUvUniform{Float32}(), 3))) isa BAT.StandardMvUniform{Float32}
        @test product_distribution(fill(BAT.StandardUvUniform{Float32}(), 3)) == BAT.StandardMvUniform{Float32}(3)
    end

    @testset "StandardMvUniform" begin
        @test @inferred(BAT.StandardMvUniform(3)) isa BAT.StandardMvUniform{Float64}
        @test @inferred(BAT.StandardMvUniform{Float32}(3)) isa BAT.StandardMvUniform{Float32}
        @test @inferred(BAT.StandardMvUniform{Float64}(3)) isa BAT.StandardMvUniform{Float64}

        @test @inferred(Distributions.Product(BAT.StandardMvUniform{Float32}(3))) isa Distributions.Product{Continuous,BAT.StandardUvUniform{Float32}}
        @test @inferred(Distributions.Product(BAT.StandardMvUniform{Float64}(3))) isa Distributions.Product{Continuous,BAT.StandardUvUniform{Float64}}
        @test @inferred(Distributions.Product(BAT.StandardMvUniform{Float64}(3))) == Distributions.Product(Fill(BAT.StandardUvUniform(), 3))
        @test @inferred(convert(Product, BAT.StandardMvUniform(3))) == Distributions.Product(Fill(BAT.StandardUvUniform(), 3))

        d = BAT.StandardMvUniform(3)
        dref = product_distribution(fill(Uniform(), 3))

        @test @inferred(eltype(typeof(d))) == eltype(typeof(dref))
        @test @inferred(eltype(d)) == eltype(dref)

        @test @inferred(length(d)) == length(dref)
        @test @inferred(size(d)) == size(dref)

        @test @inferred(view(BAT.StandardMvUniform{Float32}(7), 3)) isa BAT.StandardUvUniform{Float32}
        @test_throws BoundsError view(BAT.StandardMvUniform{Float32}(7), 9)
        @test @inferred(view(BAT.StandardMvUniform{Float32}(7), 2:4)) isa BAT.StandardMvUniform{Float32}
        @test view(BAT.StandardMvUniform{Float32}(7), 2:4) == BAT.StandardMvUniform{Float32}(3)
        @test_throws BoundsError view(BAT.StandardMvUniform{Float32}(7), 2:8)

        @test @inferred(mean(d)) == mean(dref)
        @test @inferred(var(d)) == var(dref)
        @test @inferred(cov(d)) == cov(dref)

        @test @inferred(mode(d)) == [0.5, 0.5, 0.5]
        @test @inferred(modes(d)) == Vector{Float64}[]

        @test @inferred(invcov(d)) == inv(cov(dref))
        @test @inferred(logdetcov(d)) == logdet(cov(dref))

        @test @inferred(entropy(d)) == entropy(dref)

        for x in fill.([-Inf, -1.3, 0.0, 1.3, +Inf], 3)
            @test @inferred(insupport(d, x)) == insupport(dref, x)
            @test @inferred(logpdf(d, x)) == logpdf(dref, x)
            @test @inferred(pdf(d, x)) == pdf(dref, x)
            @test @inferred(gradlogpdf(d, x)) == ForwardDiff.gradient(x -> logpdf(d, x), x)
        end

        @test @inferred(rand(stblrng(), d)) == rand(stblrng(), d)
        @test @inferred(rand!(stblrng(), d, zeros(3))) == rand!(stblrng(), d, zeros(3))
        @test @inferred(rand!(stblrng(), d, zeros(3, 10))) == rand!(stblrng(), d, zeros(3, 10))
    end
end


text:
mutable struct traceopt_type
    supp # support data
    coe # coefficient data
    constraint # "projection" or "unipotent"
    ptsupp # pure trace support
    wbasis # word basis
    tbasis # trace basis
    basis # non-trace basis
    ksupp # extending support at the k-th step
    sb # sizes of different blocks
    numb # numbers of different blocks
end

function get_sbasis(n, d, ptsupp)
    basis = Vector{UInt16}[[]]
    i = 0
    temp = UInt16[]
    while i < d+1
        if sum(temp) == n*i
            temp = ones(UInt16, i+1)
            if i < d && sum(length.(ptsupp[temp])) <= d
                push!(basis, temp)
            end
            i += 1
        else
            temp2 = copy(temp)
            j = temp[1]
            ind = findfirst(x->temp[x]!=j, 1:length(temp))
            if ind == nothing
                ind = length(temp)+1
            end
            if j != 1
                temp2[1:ind-2] = ones(UInt16, ind-2)
            end
            temp2[ind-1] = j+1
            temp = temp2
            if sum(length.(ptsupp[temp])) <= d
                push!(basis, temp)
            end
        end
    end
    return basis
end

function trace_basis(n, d, ptsupp, bsupp)
    ind = [length(bsupp[i]) <= d for i = 1:length(bsupp)]
    basis = bsupp[ind]
    tbasis = get_sbasis(length(ptsupp), d, ptsupp)
    wbasis = Vector{UInt16}[]
    for i = 1:length(tbasis), j = 1:length(basis)
        if sum(length.(ptsupp[tbasis[i]])) + length(basis[j]) <= d
            push!(wbasis, [i,j])
        end
    end
    return wbasis,tbasis,basis
end

function sym_cyclic(word)
    return min(_cyclic_canon(word), _cyclic_canon(reverse(word)))
end

function ptraceopt_first(tr_supp, coe, n, d; TS="block", monosquare=false, QUIET=false, constraint="unipotent", solve=true)
    println("***************************NCTSSOS***************************")
    println("NCTSSOS is launching...")
    bsupp = get_ncbasis(n, 2d)
    ind = [length(bsupp[i]) <= 1 || (bsupp[i][1] != bsupp[i][end] && findfirst(j -> bsupp[i][j]
    == bsupp[i][j+1], 1:length(bsupp[i])-1) == nothing) for i = 1:length(bsupp)]
    bsupp = bsupp[ind]
    ind = [sym_cyclic(bsupp[i])==bsupp[i] for i=1:length(bsupp)]
    ptsupp = bsupp[ind]
    ptsupp = ptsupp[2:end]
    sort!(ptsupp)
    supp = Vector{Vector{UInt16}}(undef, length(tr_supp))
    for i = 1:length(tr_supp)
        supp[i] = sort([ncbfind(ptsupp, length(ptsupp), tr_supp[i][j]) for j=1:length(tr_supp[i])])
    end
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    wbasis,tbasis,basis = trace_basis(n, d, ptsupp, bsupp)
    ksupp = copy(supp)
    if monosquare == true
        for i = 1:length(wbasis)
            bi1 = sort([tbasis[wbasis[i][1]]; tbasis[wbasis[i][1]]])
            bi2 = [reverse(basis[wbasis[i][2]]); basis[wbasis[i][2]]]
            constraint_reduce!(bi2, constraint=constraint)
            bi = trace_reduce(bi1, bi2, ptsupp)
            push!(ksupp, bi)
        end
    end
    sort!(ksupp)
    unique!(ksupp)
    blocks,cl,blocksize,sb,numb,_ = get_ncblocks(ksupp, ptsupp, wbasis, tbasis, basis, TS=TS, QUIET=QUIET, constraint=constraint)
    end
    if QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp = ptrace_SDP(supp, coe, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize, QUIET=QUIET, constraint=constraint, solve=solve)
    data = traceopt_type(supp, coe, constraint, ptsupp, wbasis, tbasis, basis, ksupp, sb, numb)
    return opt,data
end

function ptraceopt_higher!(data; TS="block", QUIET=false, solve=true)
    supp = data.supp
    coe = data.coe
    constraint = data.constraint
    ptsupp = data.ptsupp
    wbasis = data.wbasis
    tbasis = data.tbasis
    basis = data.basis
    ksupp = data.ksupp
    sb = data.sb
    numb = data.numb
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    blocks,cl,blocksize,sb,numb,status = get_ncblocks(ksupp, ptsupp, wbasis, tbasis, basis, sb=sb, numb=numb, TS=TS, QUIET=QUIET, constraint=constraint)
    opt = nothing
    if status == 1
        if QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure. The maximal size of blocks is $mb.")
        end
        opt,ksupp = ptrace_SDP(supp, coe, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize, QUIET=QUIET, constraint=constraint, solve=solve)
    end
    data.ksupp = ksupp
    data.sb = sb
    data.numb = numb
    return opt,data
end

function ptrace_SDP(supp, coe, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize; QUIET=false, constraint="unipotent", solve=true)
    ksupp = Vector{Vector{UInt16}}(undef, Int(sum(Int.(blocksize).^2+blocksize)/2))
    k = 1
    for i = 1:cl, j = 1:blocksize[i], r = j:blocksize[i]
        @inbounds bi1 = sort([tbasis[wbasis[blocks[i][j]][1]]; tbasis[wbasis[blocks[i][r]][1]]])
        @inbounds bi2 = [reverse(basis[wbasis[blocks[i][j]][2]]); basis[wbasis[blocks[i][r]][2]]]
        constraint_reduce!(bi2, constraint=constraint)
        @inbounds ksupp[k] = trace_reduce(bi1, bi2, ptsupp)
        k += 1
    end
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET == false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        cons = [AffExpr(0) for i=1:lksupp]
        for i = 1:cl
            bs = blocksize[i]
            if bs == 1
               @inbounds pos = @variable(model, lower_bound=0)
               @inbounds bi1 = sort([tbasis[wbasis[blocks[i][1]][1]]; tbasis[wbasis[blocks[i][1]][1]]])
               @inbounds bi2 = [reverse(basis[wbasis[blocks[i][1]][2]]); basis[wbasis[blocks[i][1]][2]]]
               constraint_reduce!(bi2, constraint=constraint)
               bi = trace_reduce(bi1, bi2, ptsupp)
               Locb = ncbfind(ksupp, lksupp, bi)
               @inbounds add_to_expression!(cons[Locb], pos)
            else
               @inbounds pos = @variable(model, [1:bs, 1:bs], PSD)
               for j = 1:blocksize[i], r = j:blocksize[i]
                   @inbounds bi1 = sort([tbasis[wbasis[blocks[i][j]][1]]; tbasis[wbasis[blocks[i][r]][1]]])
                   @inbounds bi2 = [reverse(basis[wbasis[blocks[i][j]][2]]); basis[wbasis[blocks[i][r]][2]]]
                   constraint_reduce!(bi2, constraint=constraint)
                   bi = trace_reduce(bi1, bi2, ptsupp)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   if Locb == 0
                       @error "The word does not exist!"
                       return nothing,nothing
                   end
                   if j == r
                       @inbounds add_to_expression!(cons[Locb], pos[j,r])
                   else
                       @inbounds add_to_expression!(cons[Locb], 2, pos[j,r])
                   end
               end
            end
        end
        bc = zeros(lksupp)
        for i = 1:length(supp)
            Locb = ncbfind(ksupp, lksupp, supp[i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] = coe[i]
            end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, cons.==bc)
        @objective(model, Max, lower)
        if QUIET == false
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end

function trace_reduce(word1, word2, ptsupp)
    if isempty(word2)
        ind = UInt16[]
    else
        ind = UInt16(ncbfind(ptsupp, length(ptsupp), sym_cyclic(word2)))
    end
    return sort([word1; ind])
end

function constraint_reduce!(word; constraint="unipotent")
    while length(word) > 1 && word[1] == word[end]
        deleteat!(word, 1)
        if constraint == "unipotent"
            deleteat!(word, length(word))
        end
    end
    i = 1
    while i < length(word)
        if word[i] == word[i+1]
            deleteat!(word, i)
            if constraint == "unipotent"
                deleteat!(word, i)
            end
        else
            i += 1
        end
    end
    return word
end

function get_ncgraph(ksupp, ptsupp, wbasis, tbasis, basis; constraint="unipotent")
    lb = length(wbasis)
    G = SimpleGraph(lb)
    lksupp = length(ksupp)
    for i = 1:lb, j = i+1:lb
        @inbounds bi1 = sort([tbasis[wbasis[i][1]]; tbasis[wbasis[j][1]]])
        @inbounds bi2 = [reverse(basis[wbasis[i][2]]); basis[wbasis[j][2]]]
        constraint_reduce!(bi2, constraint=constraint)
        bi = trace_reduce(bi1, bi2, ptsupp)
        if ncbfind(ksupp, lksupp, bi) != 0
            add_edge!(G, i, j)
        end
    end
    return G
end

function get_ncblocks(ksupp, ptsupp, wbasis, tbasis, basis; sb=[], numb=[], TS="block", QUIET=false, constraint="unipotent")
    if TS == false
        blocksize = [length(wbasis)]
        blocks = [[i for i=1:length(wbasis)]]
        cl = 1
    else
        G = get_ncgraph(ksupp, ptsupp, wbasis, tbasis, basis, constraint=constraint)
        if TS == "block"
            blocks = connected_components(G)
            blocksize = length.(blocks)
            cl = length(blocksize)
        else
            blocks,cl,blocksize = chordal_cliques!(G, method=TS)
        end
    end
    nsb = sort(unique(blocksize), rev=true)
    nnumb = [sum(blocksize.== i) for i in nsb]
    if isempty(sb) || nsb!=sb || nnumb!=numb
        status = 1
        if QUIET == false
            println("------------------------------------------------------")
            println("The sizes of PSD blocks:\n$nsb\n$nnumb")
            println("------------------------------------------------------")
        end
    else
        status = 0
        println("No higher NCTSSOS hierarchy!")
    end
    return blocks,cl,blocksize,nsb,nnumb,status
end

function Werner_witness_first(dY, sigma, n, d; TS="block", monosquare=false, QUIET=false, solve=true)
    println("***************************NCTSSOS***************************")
    println("NCTSSOS is launching...")
    bsupp = get_ncbasis(n, 2d)
    ind = [length(bsupp[i]) <= 1 || (bsupp[i][1] != bsupp[i][end] && findfirst(j -> bsupp[i][j]
    == bsupp[i][j+1], 1:length(bsupp[i])-1) == nothing) for i = 1:length(bsupp)]
    bsupp = bsupp[ind]
    ind = [sym_cyclic(bsupp[i])==bsupp[i] for i=1:length(bsupp)]
    ptsupp = bsupp[ind]
    ptsupp = ptsupp[2:end]
    sort!(ptsupp)
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    wbasis,tbasis,basis = trace_basis(n, d, ptsupp, bsupp)
    htrace = generate_htrace(n)
    supp = Vector{Vector{UInt16}}(undef, length(htrace))
    for i = 1:length(htrace)
        supp[i] = sort([ncbfind(ptsupp, length(ptsupp), sym_cyclic(htrace[i][j])) for j=1:length(htrace[i])])
    end
    if monosquare == true
        for i = 1:length(wbasis)
            bi1 = sort([tbasis[wbasis[i][1]]; tbasis[wbasis[i][1]]])
            bi2 = [reverse(basis[wbasis[i][2]]); basis[wbasis[i][2]]]
            constraint_reduce!(bi2, constraint="projection")
            bi = trace_reduce(bi1, bi2, ptsupp)
            push!(supp, bi)
        end
    end
    sort!(supp)
    unique!(supp)
    blocks,cl,blocksize,sb,numb,_ = get_ncblocks(supp, ptsupp, wbasis, tbasis, basis, TS=TS, QUIET=QUIET, constraint="projection")
    end
    if QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp = Werner_SDP(dY, sigma, htrace, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize, QUIET=QUIET, solve=solve)
    data = traceopt_type(htrace, dY, sigma, ptsupp, wbasis, tbasis, basis, ksupp, sb, numb)
    return opt,data
end

function Werner_witness_higher!(data; TS="block", QUIET=false, solve=true)
    htrace = data.supp
    dY = data.coe
    sigma = data.constraint
    ptsupp = data.ptsupp
    wbasis = data.wbasis
    tbasis = data.tbasis
    basis = data.basis
    ksupp = data.ksupp
    sb = data.sb
    numb = data.numb
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    blocks,cl,blocksize,sb,numb,status = get_ncblocks(ksupp, ptsupp, wbasis, tbasis, basis, sb=sb, numb=numb, TS=TS, QUIET=QUIET, constraint="projection")
    opt = nothing
    if status == 1
        if QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure. The maximal size of blocks is $mb.")
        end
        opt,ksupp = Werner_SDP(dY, sigma, htrace, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize, QUIET=QUIET, solve=solve)
    end
    data.ksupp = ksupp
    data.sb = sb
    data.numb = numb
    return opt,data
end

function Werner_SDP(dY, sigma, htrace, ptsupp, wbasis, tbasis, basis, blocks, cl, blocksize; QUIET=false, solve=true)
    ksupp = Vector{Vector{UInt16}}(undef, Int(sum(Int.(blocksize).^2+blocksize)/2))
    k = 1
    for i = 1:cl, j = 1:blocksize[i], r = j:blocksize[i]
        @inbounds bi1 = sort([tbasis[wbasis[blocks[i][j]][1]]; tbasis[wbasis[blocks[i][r]][1]]])
        @inbounds bi2 = [reverse(basis[wbasis[blocks[i][j]][2]]); basis[wbasis[blocks[i][r]][2]]]
        constraint_reduce!(bi2, constraint="projection")
        @inbounds ksupp[k] = trace_reduce(bi1, bi2, ptsupp)
        k += 1
    end
    sort!(ksupp)
    unique!(ksupp)
    lksupp = length(ksupp)
    if QUIET == false
        println("There are $lksupp affine constraints.")
    end
    objv = nothing
    if solve == true
        if QUIET == false
            println("Assembling the SDP...")
        end
        model = Model(optimizer_with_attributes(Mosek.Optimizer))
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        cons = [AffExpr(0) for i=1:lksupp]
        for i = 1:cl
            bs = blocksize[i]
            if bs == 1
               @inbounds pos = @variable(model, lower_bound=0)
               @inbounds bi1 = sort([tbasis[wbasis[blocks[i][1]][1]]; tbasis[wbasis[blocks[i][1]][1]]])
               @inbounds bi2 = [reverse(basis[wbasis[blocks[i][1]][2]]); basis[wbasis[blocks[i][1]][2]]]
               constraint_reduce!(bi2, constraint="projection")
               bi = trace_reduce(bi1, bi2, ptsupp)
               Locb = ncbfind(ksupp, lksupp, bi)
               @inbounds add_to_expression!(cons[Locb], pos)
            else
               @inbounds pos = @variable(model, [1:bs, 1:bs], PSD)
               for j = 1:blocksize[i], r = j:blocksize[i]
                   @inbounds bi1 = sort([tbasis[wbasis[blocks[i][j]][1]]; tbasis[wbasis[blocks[i][r]][1]]])
                   @inbounds bi2 = [reverse(basis[wbasis[blocks[i][j]][2]]); basis[wbasis[blocks[i][r]][2]]]
                   constraint_reduce!(bi2, constraint="projection")
                   bi = trace_reduce(bi1, bi2, ptsupp)
                   Locb = ncbfind(ksupp, lksupp, bi)
                   if Locb == 0
                       @error "The word does not exist!"
                       return nothing,nothing
                   end
                   if j == r
                       @inbounds add_to_expression!(cons[Locb], pos[j,r])
                   else
                       @inbounds add_to_expression!(cons[Locb], 2, pos[j,r])
                   end
               end
            end
        end
        bc = [AffExpr(0) for i=1:lksupp]
        coe = @variable(model, [1:length(htrace)])
        tcons = AffExpr(0)
        for i = 1:length(htrace)
            for k = 1:length(dY)
                tcons += sigma[k]*prod(tr(prod(dY[k][htrace[i][j]])) for j=1:length(htrace[i]))*coe[i]
            end
            temp = sort([ncbfind(ptsupp, length(ptsupp), sym_cyclic(htrace[i][j])) for j=1:length(htrace[i])])
            Locb = ncbfind(ksupp, lksupp, temp)
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing
            else
               bc[Locb] += coe[i]
            end
        end
        @constraint(model, tcons==-16)
        @variable(model, epsilon)
        cons[1] -= epsilon
        @constraint(model, cons.==bc)
        @objective(model, Min, epsilon)
        if QUIET == false
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
    end
    return objv,ksupp
end

function generate_htrace(n)
    return _generate_htrace(UInt16[i for i=1:n])
end

function _generate_htrace(var)
    if isempty(var)
         htrace = [Vector{UInt16}[]]
    else
        htrace = Vector{Vector{UInt16}}[]
        for i = 1:length(var)
            sset = _selete(var[2:end], i-1)
            pushfirst!.(sset, var[1])
            for mem in sset
                cbasis = _cyclic_basis(mem)
                for emem in cbasis
                    sub_htrace = _generate_htrace(UInt16.(setdiff(var, mem)))
                    for j = 1:length(sub_htrace)
                        pushfirst!(sub_htrace[j], emem)
                    end
                    append!(htrace, sub_htrace)
                end
            end
        end
    end
    return htrace
end

function _selete(var, d)
    if d > 0
        set = Vector{UInt16}[]
        for i = 1:length(var)-d+1
            iset = _selete(var[i+1:end], d-1)
            pushfirst!.(iset, var[i])
            append!(set, iset)
        end
    else
        set = [UInt16[]]
    end
    return set
end

function _cyclic_basis(var)
    basis = _permutation(var, ones(length(var)))
    return unique(_cyclic_canon.(basis))
end


text:
# This example seems to work

using Flux, Zygote
using Flux: @nograd, @epochs, onehotbatch, throttle, params
using Statistics: mean
using Plots

add_dim(x::Array) = reshape(x, (1,size(x)...))

TR = Float32
TC = ComplexF32

M = 64
constellation_dim = 2
N = 4 * M
SNR = 20
SNRlin = 10^(SNR / 10) |> TR;

encoder = Chain(Dense(M, 32, Flux.relu), Dense(32, 32, Flux.relu), Dense(32, constellation_dim))
decoder = Chain(Dense(constellation_dim, 32, Flux.relu), Dense(32, 32, Flux.relu), Dense(32, M))

function model(X)
    X_seed = Flux.onehotbatch(1:M,1:M)
    s_seed = encoder(X_seed)
    s_seed = add_dim(complex.(s_seed[1,:], s_seed[2,:]))
    norm_factor = sqrt(mean(abs.(s_seed).^2))
    s_seed = s_seed / norm_factor

    s = encoder(X)
    s = add_dim(complex.(s[1,:], s[2,:])) / norm_factor
    𝜎 = sqrt(1/SNRlin) |> TR

    r = s + 𝜎 * randn(TC, 1, N)
    r = [real(r); imag(r)]
    Y = decoder(r)
    return Y , s_seed
end

function loss(x)
    Y, s_seed = model(x)
    return Flux.logitcrossentropy(Y, x)
end

@nograd onehotbatch
X = onehotbatch(rand(1:M, N), 1:M)
@show loss(X)

opt = ADAM(0.001)
ps = params(encoder, decoder)
data = [[X]]

evalcb() = @show(loss(X));
@epochs 4000 Flux.train!(loss, ps, data, opt, cb = throttle(evalcb, 5));

Y, s_seed = model(X)

@show mean(abs.(s_seed).^2)

scatter(real(s_seed)[1,:], imag(s_seed)[1,:], aspect_ratio = :equal, markershape = :hexagon)
ylims!((-2,2))
xlims!((-2,2))

max_stars_repo_name: Rassibassi/jaude

max_issues_repo_name: Rassibassi/jaude

max_forks_repo_name: Rassibassi/jaude


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 8df38f30-b0e9-11eb-2f2a-1bb2f968d206
md"""
# [Day 20: Infinite Elves and Infinite Houses](https://adventofcode.com/2015/day/20)

To keep the Elves busy, Santa has them deliver some presents by hand, door-to-door. He sends them down a street with infinite houses numbered sequentially: `1`, `2`, `3`, `4`, `5`, and so on.

Each Elf is assigned a number, too, and delivers presents to houses based on that number:

-   The first Elf (number `1`) delivers presents to every house: `1`, `2`, `3`, `4`, `5`, ....
-   The second Elf (number `2`) delivers presents to every second house: `2`, `4`, `6`, `8`, `10`, ....
-   Elf number `3` delivers presents to every third house: `3`, `6`, `9`, `12`, `15`, ....

There are infinitely many Elves, numbered starting with `1`. Each Elf delivers presents equal to ten times his or her number at each house.

So, the first nine houses on the street end up like this:

```
House 1 got 10 presents.
House 2 got 30 presents.
House 3 got 40 presents.
House 4 got 70 presents.
House 5 got 60 presents.
House 6 got 120 presents.
House 7 got 80 presents.
House 8 got 150 presents.
House 9 got 130 presents.

```

The first house gets `10` presents: it is visited only by Elf `1`, which delivers `1 * 10 = 10` presents. The fourth house gets `70` presents, because it is visited by Elves `1`, `2`, and `4`, for a total of `10 + 20 + 40 = 70` presents.

**What is the lowest house number of the house to get at least as many presents as the number in your puzzle input?**
"""

# ╔═╡ 3508f641-2704-4bd5-9d44-987b29e87cf2
puzzle_input = "36000000"

# ╔═╡ 70139d68-0c2f-4d14-a753-fa6c5e573f35

# ╔═╡ Cell order:
# ╟─8df38f30-b0e9-11eb-2f2a-1bb2f968d206
# ╟─3508f641-2704-4bd5-9d44-987b29e87cf2
# ╠═70139d68-0c2f-4d14-a753-fa6c5e573f35


text:
export @threaded;

macro threaded(expr::Expr)
	local dict = Dict{Symbol, Symbol}(
		:foreach => :tforeach,
		:mapfoldl => :tmapfoldl,
		:mapfoldr => :tmapfoldr,
		:map => :tmap,
		:map! => :tmap!,
		:mapreduce => :tmaptreduce,
		:reduce => :treduce
	);

	local function replace_threaded(expr::Expr)
		if expr.head == :call
			try
				expr.args[1] = dict[expr.args[1]];
			catch
			end
		end
		for arg in expr.args
			if typeof(arg) == Expr
				replace_threaded(arg);
			end
		end
		return expr;
	end

	return quote
		$(esc(replace_threaded(expr)))
	end
end


text:
# Huber loss function

"""
  HuberLoss(rho::Real=1.0, mu::Real=1.0)

Returns the function `g(x) = (mu/2)||x||² if ||x|| ⩽ rho, and rho*mu*(||x||-rho/2) otherwise`.
"""

immutable HuberLoss{R <: Real} <: ProximableFunction
  rho::R
  mu::R
  function HuberLoss(rho::R, mu::R)
    if rho <= 0.0 || mu <= 0.0
      error("parameters rho and mu must be positive")
    else
      new(rho, mu)
    end
  end
end

HuberLoss{R <: Real}(rho::R=1.0, mu::R=1.0) = HuberLoss{R}(rho, mu)

function (f::HuberLoss){T <: Union{Real, Complex}}(x::AbstractArray{T})
  normx = vecnorm(x)
  if normx <= f.rho
    return (f.mu/2)*normx^2
  else
    return f.rho*f.mu*(normx-f.rho/2)
  end
end

function prox!{T <: Union{Real, Complex}}(f::HuberLoss, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  normx = vecnorm(x)
  mugam = f.mu*gamma
  scal = (1-min(mugam/(1+mugam), mugam*f.rho/(normx)))
  for k in eachindex(y)
    y[k] = scal*x[k]
  end
  normy = scal*normx
  if normy <= f.rho
    return (f.mu/2)*normy^2
  else
    return f.rho*f.mu*(normy-f.rho/2)
  end
end

fun_name(f::HuberLoss) = "Huber loss"
fun_dom(f::HuberLoss) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::HuberLoss) = "x ↦ (μ/2)||x||² if ||x||⩽ρ, μρ(||x||-ρ/2) otherwise"
fun_params(f::HuberLoss) = string("ρ = $(f.rho), μ = $(f.mu)")

function prox_naive{T <: Union{Real, Complex}}(f::HuberLoss, x::AbstractArray{T}, gamma::Real=1.0)
  y = (1-min(f.mu*gamma/(1+f.mu*gamma), f.mu*gamma*f.rho/(vecnorm(x))))*x
  if vecnorm(y) <= f.rho
    return y, (f.mu/2)*vecnorm(y)^2
  else
    return y, f.rho*f.mu*(vecnorm(y)-f.rho/2)
  end
end


text:
module MyFunRep

# Write your package code here.
"""
A simple add method for testing
    ```julia> my_fun1(2,3)
    ```
"""
my_fun1(x,y) = x+y

"""
A sample function as an experiment
    ```julia> my_fun2(1,2,3)
"""
function my_fun2(x,y,z)
    return x+y+z
end

"""
    func(x)
Returns double the number `x` plus `1`.
```julia-repl
julia > func(2)
```
"""
func(x) = 2x + 1

"""
    func2(x,y)
Returns something else double the number `x` plus `1`.
"""
func2(x,y) = 2x + 1 + 3y

export my_fun1, my_fun2, func,func2

end


text:
module TestTimescaleDB

using Test
using TSML

function test_timescaledb()
  tsdb = TimescaleDB()
  mpipeline = Pipeline(Dict(
         :transformers => [tsdb]
     )
  )
  fit!(mpipeline)
  respipe = transform!(mpipeline)
  @test sum(respipe.Value) ≈ 560997.0
  @test nrow(respipe) == 48349
  @test ncol(respipe) == 2
end
@testset "TimescaleDB" begin
  test_timescaledb()
end

end


text:
##############################################################################
#
# REFERENCES: "Statistical Distributions"
#
##############################################################################

immutable EmpiricalUnivariateDistribution <: ContinuousUnivariateDistribution
	values::Vector{Float64}
    support::Vector{Float64}
	cdf::Function
	entropy::Float64
	kurtosis::Float64
	mean::Float64
	median::Float64
	modes::Vector{Float64}
	skewness::Float64
	var::Float64
end

function EmpiricalUnivariateDistribution(x::Vector)
    sx = sort(x)
	EmpiricalUnivariateDistribution(sx,
                                    unique(sx),
	                                ecdf(x),
	                                NaN,
	                                NaN,
	                                mean(x),
	                                median(x),
	                                Float64[],
	                                NaN,
	                                var(x))
end

cdf(d::EmpiricalUnivariateDistribution, q::Real) = d.cdf(q)

entropy(d::EmpiricalUnivariateDistribution) = d.entropy

kurtosis(d::EmpiricalUnivariateDistribution) = d.kurtosis

mean(d::EmpiricalUnivariateDistribution) = d.mean

median(d::EmpiricalUnivariateDistribution) = d.median

modes(d::EmpiricalUnivariateDistribution) = Float64[]

function pdf(d::EmpiricalUnivariateDistribution, x::Real)
    ## TODO: Create lookup table for discrete case
    1.0 / length(d.values)
end

function quantile(d::EmpiricalUnivariateDistribution, p::Real)
    n = length(d.values)
    index = ifloor(p * n) + 1
    index > n ? d.values[n] : d.values[index]
end

function rand(d::EmpiricalUnivariateDistribution)
    d.values[rand(1:length(d.values))]
end

skewness(d::EmpiricalUnivariateDistribution) = NaN

var(d::EmpiricalUnivariateDistribution) = d.var

### handling support

insupport(d::EmpiricalUnivariateDistribution, x::Real) = contains(d.support, x)

isupperbounded(::Union(EmpiricalUnivariateDistribution, Type{EmpiricalUnivariateDistribution})) = true
islowerbounded(::Union(EmpiricalUnivariateDistribution, Type{EmpiricalUnivariateDistribution})) = true
isbounded(::Union(EmpiricalUnivariateDistribution, Type{EmpiricalUnivariateDistribution})) = true

hasfinitesupport(d::Union(EmpiricalUnivariateDistribution, Type{EmpiricalUnivariateDistribution})) = true
minimum(d::EmpiricalUnivariateDistribution) = d.values[1]
maximum(d::EmpiricalUnivariateDistribution) = d.values[end]
support(d::EmpiricalUnivariateDistribution) = d.support

### fit model

function fit_mle{T <: Real}(::Type{EmpiricalUnivariateDistribution},
	                    x::Vector{T})
	EmpiricalUnivariateDistribution(x)
end


text:
module ArgoTools

using Dates, YAML, NCDatasets, CSV, DataFrames, Interpolations, Printf
import ArgoData.ProfileNative

"""
    meta_init(fil::String)

Get parameters to call `MITprof.format` from yaml file (`fil`, e.g. "../examples/ArgoToMITprof.yml").
"""
function meta_init(fil="")

    if !isempty(fil)
        meta=YAML.load(open(fil))
    else
        z_std=[[5:10:185]... ; [200:20:500]... ; [550:50:1000]...; [1100:100:6000]...]
        meta=Dict("name" => "argo",
        "variables"   => ["depth", "T", "S"],
        "subset"      => Dict{Any, Any}("year"=>2016, "basin"=>"indian"),
        "depthLevels" => z_std,
        "dirOut"      => joinpath(tempdir(),"argo_example"),
        "depthRange"  => [0, 2000],
        "dirIn"       =>  joinpath(tempdir(),"argo_example"),
        )
    end

    d=meta["dirIn"]
    b=meta["subset"]["basin"]
    y=meta["subset"]["year"]

    list0=Array{Array,1}(undef,12)
    for m=1:12
        sd="$b"*"_ocean/$y/"*Printf.@sprintf("%02d/",m)
        if isdir(d*sd)
            tmp=readdir(d*sd)
            list0[m]=[sd*tmp[i] for i=1:length(tmp)]
        else
            list0[m]=[]
        end
    end

    nf=sum(length.(list0))
    list1=Array{String,1}(undef,nf)
    f=0
    for m=1:12
        for ff=1:length(list0[m])
            f+=1
            list1[f]=list0[m][ff]
        end
    end

    meta["fileInList"]=list1;

    #2. coordinate

    z_std=meta["depthLevels"]
    if length(z_std)>1
        tmp1=(z_std[2:end]+z_std[1:end-1])/2
        z_top=[z_std[1]-(z_std[2]-z_std[1])/2;tmp1]
        z_bot=[tmp1;z_std[end]+(z_std[end]-z_std[end-1])/2]
    else
        z_top=0.9*z_std
        z_bot=1.1*z_std
    end

    meta["z_std"]=z_std
    meta["z_top"]=z_top
    meta["z_bot"]=z_bot

    #3. various other specs

    meta["inclZ"] = false
    meta["inclT"] = true
    meta["inclS"] = true
    meta["inclU"] = false
    meta["inclV"] = false
    meta["inclPTR"] = false
    meta["inclSSH"] = false
    meta["TPOTfromTINSITU"] = true

    meta["doInterp"] = true
    meta["addGrid"] = true
    meta["outputMore"] = false
    meta["method"] = "interp"
    meta["fillval"] = -9999.0
    meta["buffer_size"] = 10000

    meta["fileOut"]=meta["name"]*"_"*meta["subset"]["basin"]*"_"*string(meta["subset"]["year"])*".nc"
    meta["var_out"]=meta["variables"]

    return meta
end

"""
    meta(input_file,output_file)

Get parameters to call `MITprof.format` which will read from `input_file` to create `output_file`.
"""
function meta(input_file,output_file)
    meta=ArgoTools.meta_init()
    #f=1
    #input_file=meta["dirIn"]*meta["fileInList"][f]
    meta["fileOut"]=output_file

    meta
end

"""
    GetOneProfile(m)

Get one profile from netcdf file.
"""
function GetOneProfile(ds,m)

    #
    t=ds["JULD"][m]
    if !ismissing(t)
        ymd=Dates.year(t)*1e4+Dates.month(t)*1e2+Dates.day(t)
        hms=Dates.hour(t)*1e4+Dates.minute(t)*1e2+Dates.second(t)

        prof_date=t-DateTime(0)
        prof_date=prof_date.value/86400/1000+1 #+1 is to match Matlab's datenum
    else
        ymd=0
        hms=0
        prof_date=0.0
    end

    lat=ds["LATITUDE"][m]
    lon=ds["LONGITUDE"][m]
    (!ismissing(lon))&&(lon < 0.0) ? lon=lon+360.0 : nothing

    direction=ds["DIRECTION"][m]
    direc=0
    direction=='A' ? direc=1 : nothing
    direction=='D' ? direc=2 : nothing

    #
    pnum_txt=ds["PLATFORM_NUMBER"][:,m]
    #pnum_txt=pnum_txt[findall((!ismissing).(pnum_txt))]
    ii=findall(skipmissing(in.(pnum_txt,"0123456789")))
    ~isempty(ii) ? pnum_txt=String(vec(Char.(pnum_txt[ii]))) : pnum_txt="9999"
    pnum=parse(Int,pnum_txt)

    #
    p=ds["PRES_ADJUSTED"][:,m]
    p_QC=ds["PRES_ADJUSTED_QC"][:,m]
    if sum((!ismissing).(p))==0
        p=ds["PRES"][:,m]
        p_QC=ds["PRES_QC"][:,m]
    end

    #set qc to 5 if missing
    p_QC[findall(ismissing.(p))].='5'
    #avoid potential duplicates
    for n=1:length(p)-1
        if ~ismissing(p[n])
            tmp1=( (!ismissing).(p[n+1:end]) ).&( p[n+1:end].==p[n] )
            tmp1=findall(tmp1)
            p[n.+tmp1].=missing
            p_QC[n.+tmp1].='5'
        end
    end

    #position and date
    isBAD=0
    ismissing(ds["POSITION_QC"][m])||~in(ds["POSITION_QC"][m],"1258") ? isBAD=1 : nothing
    ismissing(ds["JULD_QC"][m])||~in(ds["JULD_QC"][m],"1258") ? isBAD=1 : nothing

    #pressure
    tmp1=findall( (!in).(p_QC,"1258") )
    if (length(tmp1)<=5)&&(length(tmp1)>0)
        #omit these few bad points but keep the profile
        p[tmp1].=missing
    elseif length(tmp1)>5
        #flag the profile (will be omitted later)
        #but keep the bad points (for potential inspection)
        isBAD=1
    end

    #temperature
    t=ds["TEMP_ADJUSTED"][:,m]
    t_QC=ds["TEMP_ADJUSTED_QC"][:,m]
    t_ERR=ds["TEMP_ADJUSTED_ERROR"][:,m]
    t_ERR[findall( (ismissing).(t_ERR) )].=0.0

    if sum((!ismissing).(t))==0
        t=ds["TEMP"][:,m]
        t_QC=ds["TEMP_QC"][:,m]
    end

    #salinity
    if haskey(ds,"PSAL")
        s=ds["PSAL_ADJUSTED"][:,m]
        s_QC=ds["PSAL_ADJUSTED_QC"][:,m]
        s_ERR=ds["PSAL_ADJUSTED_ERROR"][:,m]
        s_ERR[findall( (ismissing).(t_ERR) )].=0.0
        if sum((!ismissing).(t))==0
            s=ds["PSAL"][:,m]
            s_QC=ds["PSAL_QC"][:,m]
        end
    else
        s=fill(missing,size(t))
        s_QC=Char.(32*ones(size(t_QC)))
    end

    if ismissing(t[1]) #this file does not contain temperature data...
        t=fill(missing,size(t))
        t_ERR=fill(0.0,size(t))
    else #apply QC
        tmp1=findall(skipmissing((!in).(t_QC,"1258")))
        t[tmp1].=missing
    end

    if ismissing(s[1]) #this file does not contain salinity data...
        s=fill(missing,size(s))
        s_ERR=fill(0.0,size(s))
    else #apply QC
        tmp1=findall(skipmissing((!in).(s_QC,"1258")))
        s[tmp1].=missing
    end

    prof=ProfileNative(
        convert(Array{Union{Float64,Missing}},[lon]),
        convert(Array{Union{Float64,Missing}},[lat]),
        convert(Union{Float64,Missing},prof_date),
        convert(Array{Union{Int,Missing}},[ymd]),
        convert(Array{Union{Int,Missing}},[hms]),
        convert(Array{Union{Float64,Missing}},t),
        convert(Array{Union{Float64,Missing}},s),
        convert(Array{Union{Float64,Missing}},p),
        convert(Array{Union{Float64,Missing}},p), #place holder for depth
        convert(Array{Union{Float64,Missing}},t_ERR),
        convert(Array{Union{Float64,Missing}},s_ERR),
        convert(Array{Union{String,Missing}},[pnum_txt]),
        convert(Array{Union{Int,Missing}},[direc]),
        convert(Array{Union{Char,Missing}},[ds["DATA_MODE"][m]]),
        convert(Array{Union{Int,Missing}},[isBAD])
        )

    return prof
end

## ArgoProfile data structure

"""
    prof_PtoZ!(prof)

Convert prof["p"] to prof["depth"]
"""
function prof_PtoZ!(prof)
    l=prof.lat[1]
    k=findall((!ismissing).(prof.pressure))
    prof.depth[k].=[sw_dpth(Float64(prof.pressure[kk]),Float64(l)) for kk in k]
end

"""
    prof_TtoΘ!(prof)

Convert prof["T"] to potential temperature
"""
function prof_TtoΘ!(prof)
    T=prof.T
    P=0.981*1.027*prof.depth
    S=35.0*ones(size(T))
    k=findall( ((!ismissing).(T)).*((!ismissing).(P)) )
    T[k]=[sw_ptmp(Float64(S[kk]),Float64(T[kk]),Float64(P[kk])) for kk in k]
end

"""
    prof_interp!(prof,meta)

Interpolate from prof["depth"] to meta["z_std"]
"""
function prof_interp!(prof,prof_std,meta)
    z_std=meta["z_std"]
    for ii=2:length(meta["var_out"])
        v=meta["var_out"][ii]
        v_e=v*"_ERR"

        t_std=similar(z_std,Union{Missing,Float64})
        e_std=similar(z_std,Union{Missing,Float64})

        z=prof.depth
        v=="T" ? t=prof.T : t=prof.S
        do_e=true #haskey(prof,v_e)
        v=="T" ? e=prof.T_ERR : e=prof.S_ERR

        kk=findall((!ismissing).(z.*t))
        if (meta["doInterp"])&&(length(kk)>1)
            z_in=z[kk]; t_in=t[kk]
            do_e ? e_in=e[kk] : nothing
            k=sort(1:length(kk),by= i -> z_in[i])
            z_in=z_in[k]; t_in=t_in[k]
            do_e ? e_in=e_in[k] : nothing
            #omit values outside observed range:
            D0=minimum(skipmissing(z_in))
            D1=maximum(skipmissing(z_in))
            msk1=findall( (z_std.<D0).|(z_std.>D1) )
            #avoid duplicates:
            msk2=findall( ([false;(z_in[1:end-1]-z_in[2:end]).==0.0]).==true )
            if length(kk)>5
                interp_linear_extrap = LinearInterpolation(Float64.(z_in), Float64.(t_in), extrapolation_bc=Line())
                t_std[:] = interp_linear_extrap(z_std)
                t_std[msk1].=missing
                t_std[msk2].=missing
                if do_e
                    e_in[findall(ismissing.(e_in))].=0.0
                    interp_linear_extrap = LinearInterpolation(Float64.(z_in), Float64.(e_in), extrapolation_bc=Line())
                    e_std[:] = interp_linear_extrap(z_std)
                    e_std[msk1].=missing
                    e_std[msk2].=missing
                end
            else
                t_std = []
                e_std = []
            end
        end
        if v=="T"&&!isempty(t_std)
            prof_std.T .=t_std
            prof_std.T_ERR .=e_std
        elseif v=="S"&&!isempty(t_std)
            prof_std.S .=t_std
            prof_std.S_ERR .=e_std
        end
    end
end

function prof_test_set1!(prof,prof_std,meta)

    #test for 'not enough data near standard level'
    z=prof.depth
    nz1=length(z)
    nz2=length(meta["z_top"])
    tmp1=[meta["z_top"][j]-z[i] for i in 1:nz1, j in 1:nz2]
    tmp2=[meta["z_bot"][j]-z[i] for i in 1:nz1, j in 1:nz2]
    tmp1=1.0*(tmp1.<0.0)
    tmp2=1.0*(tmp2.>=0.0)
    tmp3=tmp1.*tmp2
    if isa(tmp3,Matrix{Missing})
        tmp3=convert(Array{Union{Float64,Missing}},tmp3)
    end
    tmp3[findall(ismissing.(tmp3))].=0.0
    tmp3=sum(tmp3,dims=1)

    prof_std.Ttest.=0.0
    prof_std.Ttest[findall( (tmp3[:].<=0.0).&((!ismissing).(prof_std.T[:])) )].=1.0
    prof_std.Stest.=0.0
    prof_std.Stest[findall( (tmp3[:].<=0.0).&((!ismissing).(prof_std.S[:])) )].=1.0

    #test for "absurd" salinity values :
    prof_std.Stest[findall( (prof_std.S[:].>42.0).&((!ismissing).(prof_std.S[:])) )].=
    10*prof_std.Stest[findall( (prof_std.S[:].>42.0).&((!ismissing).(prof_std.S[:])) )] .+ 2
    prof_std.Stest[findall( (prof_std.S[:].<15.0).&((!ismissing).(prof_std.S[:])) )].=
    10*prof_std.Stest[findall( (prof_std.S[:].<15.0).&((!ismissing).(prof_std.S[:])) )] .+ 2

    #bad pressure flag:
    if prof.isBAD[1]>0
        prof_std.Ttest.=10*prof_std.Ttest .+ 6
        prof_std.Stest.=10*prof_std.Stest .+ 6
    end;

    #Argo grey list:
    if haskey(meta,"greylist") && !isempty(meta["greylist"])
        test1=!(prof.DATA_MODE[1].=='D') #true = real time profile ('R' or 'A')
        test2=sum(parse(Int,prof.pnum_txt[1]).==meta["greylist"][:,"PLATFORM_CODE"]) #is in grey list
        if test1&(test2>0)
            II=findall(parse(Int,prof.pnum_txt[1]).==meta["greylist"][:,"PLATFORM_CODE"])
            timeP=prof.ymd[1]
            for ii in II
                time0=meta["greylist"][ii,"START_DATE"]
                time1=meta["greylist"][ii,"END_DATE"]
                if (time0<timeP)&&(ismissing(time1)||(time1>timeP))
                    prof_std.Ttest.=10*prof_std.Ttest .+ 4
                    prof_std.Stest.=10*prof_std.Stest .+ 4
                end
            end
        end
    end

end

function prof_test_set2!(prof_std,meta)
    max_cost=50

    tmp_cost=prof_std.Tweight.*((prof_std.T-prof_std.Testim).^2)
    tmp_cost=convert(Array{Union{Float64,Missing}},tmp_cost)
    tmp_cost[findall(ismissing.(tmp_cost))].=0.0
    ii=findall( (tmp_cost.>max_cost) )
    prof_std.Ttest[ii].=10*prof_std.Ttest[ii] .+5

    tmp_cost=prof_std.Sweight.*((prof_std.S-prof_std.Sestim).^2)
    tmp_cost=convert(Array{Union{Float64,Missing}},tmp_cost)
    tmp_cost[findall(ismissing.(tmp_cost))].=0.0
    ii=findall( (tmp_cost.>max_cost) )
    prof_std.Stest[ii].=10*prof_std.Stest[ii] .+5

    ii=findall(
        ((ismissing).(prof_std.T)).+((ismissing).(prof_std.Testim)).+(prof_std.Ttest.>0).>0
        )
    prof_std.Tweight[ii].=0

    ii=findall(
        ((ismissing).(prof_std.S)).+((ismissing).(prof_std.Sestim)).+(prof_std.Stest.>0).>0
        )
    prof_std.Sweight[ii].=0

end

"""
    prof_convert!(prof,meta)

Appply conversions to variables (lon,lat,depth,temperature) in `prof` if specified in `meta`.
"""
function prof_convert!(prof,meta)
    lonlatISbad=false
    ismissing(prof.lat[1])||(prof.lat[1]<-90.0)||(prof.lat[1]>90.0) ? lonlatISbad=true : nothing
    ismissing(prof.lon[1])||(prof.lon[1]<-180.0)||(prof.lon[1]>360.0) ? lonlatISbad=true : nothing

    #if needed then reset lon,lat after issuing a warning
    lonlatISbad==true ? println("warning: out of range lon/lat was reset to 0.0,-89.99") : nothing
    lonlatISbad ? (prof.lon[1],prof.lat[1])=(0.0,-89.99) : nothing

    #if needed then fix longitude range to 0-360
    (~lonlatISbad)&(prof.lon[1]>180.0) ? prof.lon[1]-=360.0 : nothing

    #if needed then convert pressure to depth
    (~meta["inclZ"])&(~lonlatISbad) ? ArgoTools.prof_PtoZ!(prof) : nothing

    #if needed then convert T to potential temperature θ
    meta["TPOTfromTINSITU"] ? ArgoTools.prof_TtoΘ!(prof) : nothing
end

function interp_z(x,y,xi)
    jj=findall(isfinite.(y))
    interp_linear_extrap = LinearInterpolation(Float64.(x[jj]), Float64.(y[jj]), extrapolation_bc=Flat())
    return interp_linear_extrap(xi)
end

"""
    sw_dpth(P,LAT)

Calculate depth in meters from pressure (P; in decibars) and
latitude (LAT; in °N)

```
d = sw_dpth(100.0,20.0)
```
"""
function sw_dpth(P,LAT)
    # Original author:  Phil Morgan 92-04-06  (morgan@ml.csiro.au)
    # Reference: Unesco 1983. Algorithms for computation of fundamental properties of
    # seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp. Eqn 25, p26

    c1 = 9.72659
    c2 = -2.2512E-5
    c3 = 2.279E-10
    c4 = -1.82E-15
    gam_dash = 2.184e-6

    LAT = abs.(LAT)
    X   = sin.(deg2rad.(LAT))
    X   = X.*X
    bot_line = 9.780318*(1.0+(5.2788E-3+2.36E-5*X).*X) + gam_dash*0.5*P
    top_line = (((c4*P+c3).*P+c2).*P+c1).*P
    DEPTHM   = top_line./bot_line

    return DEPTHM
end

"""
    sw_ptmp(S,T,P,PR)

Calculate potential temperature as per UNESCO 1983 report from salinity (S;
in psu), in situ temperature (T; in °C), and pressure (P; in decibar)
relative to PR (in decibar; 0 by default).

```
ptmp = sw_ptmp(S,T,P,PR=missing)
```
"""
function sw_ptmp(S,T,P,PR=0.0)
# Original author:  Phil Morgan
# Reference: Fofonoff, P. and Millard, R.C. Jr
#    Unesco 1983. Algorithms for computation of fundamental properties of
#    seawater. _Unesco Tech. Pap. in Mar. Sci._, No. 44. Eqn.(31) p.39
#    Bryden, H. 1973.
#    "New Polynomials for thermal expansion, adiabatic temperature gradient
#    and potential temperature of sea water." DEEP-SEA RES., 1973, Vol 20

# theta1
del_P  = PR - P
del_th = del_P.*sw_adtg(S,T,P);
th     = T + 0.5*del_th;
q      = del_th;

# theta2
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
th     = th + (1 - 1/sqrt(2))*(del_th - q);
q      = (2-sqrt(2))*del_th + (-2+3/sqrt(2))*q;

# theta3
del_th = del_P.*sw_adtg(S,th,P+0.5*del_P);
th     = th + (1 + 1/sqrt(2))*(del_th - q);
q      = (2 + sqrt(2))*del_th + (-2-3/sqrt(2))*q;

# theta4
del_th = del_P.*sw_adtg(S,th,P+del_P);
PT     = th + (del_th - 2*q)/6;

return PT
end

"""
    sw_adtg(S,T,P)

Calculate adiabatic temperature gradient as per UNESCO 1983 routines from salinity
(S; in psu), in situ temperature (T; in °C), and pressure (P; in decibar)
```
adtg = sw_adtg(S,T,P)
```
"""
function sw_adtg(S,T,P)
    # Original author:  Phil Morgan
    # Reference: Fofonoff, P. and Millard, R.C. Jr
    #    Unesco 1983. Algorithms for computation of fundamental properties of
    #    seawater. _Unesco Tech. Pap. in Mar. Sci._, No. 44. Eqn.(31) p.39
    #    Bryden, H. 1973.
    #    "New Polynomials for thermal expansion, adiabatic temperature gradient
    #    and potential temperature of sea water." DEEP-SEA RES., 1973, Vol 20

    a0 =  3.5803E-5;
    a1 =  8.5258E-6;
    a2 = -6.836E-8;
    a3 =  6.6228E-10;

    b0 =  1.8932E-6;
    b1 = -4.2393E-8;

    c0 =  1.8741E-8;
    c1 = -6.7795E-10;
    c2 =  8.733E-12;
    c3 = -5.4481E-14;

    d0 = -1.1351E-10;
    d1 =  2.7759E-12;

    e0 = -4.6206E-13;
    e1 =  1.8676E-14;
    e2 = -2.1687E-16;

    ADTG =      a0 + (a1 + (a2 + a3.*T).*T).*T + (b0 + b1.*T).*(S-35.0) +
    ( (c0 + (c1 + (c2 + c3.*T).*T).*T) + (d0 + d1.*T).*(S-35.0) ).*P +
    (  e0 + (e1 + e2.*T).*T ).*P.*P

    return ADTG
end

"""
    monthly_climatology_factors(date)

if `date` is a date in `days since DateTime(0)`

and `ff(rec)` returns a value for month `rec` in `1:12`

then compute the factors to interpolate from `rec[1],rec[2]` to `date`.

For example :

```
ff(x)=sin((x-0.5)/12*2pi)
fac,rec=monthly_climatology_factors(prof["date"])

gg=fac[1]*ff(rec[1])+fac[2]*ff(rec[2])
(ff(rec[1]),gg,ff(rec[2]))
```
"""
function monthly_climatology_factors(date)

    tmp2=ones(13,1)*[1991 1 1 0 0 0]; tmp2[1:12,2].=(1:12); tmp2[13,1]=1992.0;
    tmp2=[DateTime(tmp2[i,:]...) for i in 1:13]

    tim_fld=tmp2 .-DateTime(1991,1,1);
    tim_fld=1/2*(tim_fld[1:12]+tim_fld[2:13])
    tim_fld=[tim_fld[i].value for i in 1:12]/86400/1000

    tim_fld=[tim_fld[12]-365.0;tim_fld...;tim_fld[1]+365.0]
    rec_fld=[12;1:12;1]

    year0=year(DateTime(0,1,1)+Day(Int(floor(date))))
    date0=DateTime(year0,1,1)-DateTime(0)

    date0=date0.value/86400/1000+1 #+1 is to match Matlab's datenum
    tim_prof=date-date0
    tim_prof>365.0 ? tim_prof=365.0 : nothing

    tt=maximum(findall(tim_fld.<=tim_prof))
    a0=(tim_prof-tim_fld[tt])/(tim_fld[tt+1]-tim_fld[tt])

    return (1-a0,a0),(rec_fld[tt],rec_fld[tt+1])
end

end #module ArgoTools

module GriddedFields

using MeshArrays, OceanStateEstimation, Statistics

function NaN_mask(Γ)
    msk=write(Γ.hFacC)
    msk[findall(msk.>0.0)].=1.0
    msk[findall(msk.==0.0)].=NaN
    msk=MeshArrays.read(msk,Γ.hFacC)
end

function MonthlyClimatology(fil,msk)
    fid = open(fil)
    tmp = Array{Float32,4}(undef,(90,1170,50,12))
    read!(fid,tmp)
    tmp = hton.(tmp)
    close(fid)

    T=Array{MeshArray,1}(undef,12)
    for tt=1:12
        T[tt]=msk*MeshArrays.read(tmp[:,:,:,tt],msk.grid)
    end

    return T
end

function AnnualClimatology(fil,msk)
    fid = open(fil)
    tmp=Array{Float32,3}(undef,(90,1170,50))
    read!(fid,tmp)
    tmp = hton.(tmp)
    close(fid)

    T=MeshArray(msk.grid,Float64,50)
    T=msk*MeshArrays.read(convert(Array{Float64},tmp),T)
    return T
end

"""
    load()

Load gridded fields from files (and download the files if needed).
Originally this function returned `(Γ=Γ,msk,T=T,S=S,σT=σT,σS=σS)`.
"""
function load()
    γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
    Γ=GridLoad(γ,option="full")

    msk=NaN_mask(Γ)

    pth=MITPROFclim_path
    OceanStateEstimation.MITPROFclim_download()

    T=MonthlyClimatology(pth*"T_OWPv1_M_eccollc_90x50.bin",msk)
    S=MonthlyClimatology(pth*"S_OWPv1_M_eccollc_90x50.bin",msk)
    σT=AnnualClimatology(pth*"sigma_T_nov2015.bin",msk)
    σS=AnnualClimatology(pth*"sigma_S_nov2015.bin",msk)

    tmp=σT.grid.write(σT)
    for kk in 1:size(tmp,3)
        tmp1=tmp[:,:,kk]; tmp2=tmp1[findall((!isnan).(tmp1))];
        tmp3=max(quantile(tmp2,0.05),1e-3)
        tmp1[findall(tmp1.<tmp3)].=tmp3
        tmp[:,:,kk].=tmp1
    end
    σT=σT.grid.read(tmp,σT.grid)

    tmp=σS.grid.write(σS)
    for kk in 1:size(tmp,3)
        tmp1=tmp[:,:,kk]; tmp2=tmp1[findall((!isnan).(tmp1))];
        tmp3=max(quantile(tmp2,0.05),1e-3)
        tmp1[findall(tmp1.<tmp3)].=tmp3
        tmp[:,:,kk].=tmp1
    end
    σS=σS.grid.read(tmp,σS.grid)

    (Γ=Γ,msk,T=T,S=S,σT=σT,σS=σS)
end

function interp_h(z_in::MeshArray,f,i,j,w,z_out)
    for k in 1:50
        z_out[k]=NaN
        if !isnan(sum(w[1,:]))
            x=[z_in[f[1,ii],k][i[1,ii],j[1,ii]] for ii=1:4]
            kk=findall(isfinite.(x))
            ~isempty(kk) ? z_out[k]=sum(w[1,kk].*x[kk])/sum(w[1,kk]) : nothing
        end
    end
end

end #module GriddedFields


text:
import Base: isnull, get

export BestEntity, isnull, get, reset!, set!, update!, release!

type BestEntity{T<:AbstractComparer}
    comparer::T
    value::Nullable{Entity}
end

BestEntity(comparer) = BestEntity(comparer, Nullable{Entity}())

release!(best::BestEntity) = reset!(best)

isnull(best::BestEntity) = isnull(best.value)

get(best::BestEntity) = get(best.value)

reset!(best::BestEntity) = set!(best, Nullable{Entity}())

set!(best::BestEntity, value::Entity) = set!(best, Nullable{Entity}(value))

function set!(best::BestEntity, value::Nullable{Entity})
    if !isnull(best)
        release!(get(best))
    end
    if isnull(value)
        best.value = Nullable{Entity}()
    else
        best.value = copy(get(value))
    end
end

update!(best::BestEntity, value::Entity) = update!(best, Nullable{Entity}(value))

function update!(best::BestEntity, other::Entity)
    if isnull(best) || fn(best.comparer)(other.body, get(best).body)
        set!(best, other)
        return true
    end
    false
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
const δceil_const = true
const δcolon_const = true
const δfloor_const = true
const δlength_const = true
const δones_const = true
const δrand_const = true
const δrandn_const = true
const δref_const2 = true
const δround_const = true
const δsign_const = true
const δsize_const = true
const δsrand_const = true
const δtrunc_const = true
const δtuple = true
const δzeros_const = true

function fit{T}(::Type{T}, x, sz::Tuple{Int, Int})::T
    if size(x) == sz
        x
    elseif sz[1] == 1 && sz[2] != 1
        sum(x, 1)
    elseif sz[2] == 1 && sz[1] != 1
        sum(x, 2)
    else
        fill(sum(x), sz)
    end
end
function fit{T}(::Type{T}, x, sz::Tuple{Int})::T
    if size(x) == sz
        vec(x)
    elseif sz[1] == 1
        fill(sum(x), sz)
    else
        vec(sum(x, 2))
    end
end
fit{T}(::Type{T}, x, sz::Tuple{})::T = sum(x)

safediv{T}(x::T, y) = y == 0 ? T(0) : x / y
δabs(x) = (abs(x), z->z*sign(x))
δacos(x) = (acos(x), z->-z/sqrt(1-x*x))
δasin(x) = (asin(x), z->z/sqrt(1-x*x))
δatan(x) = (atan(x), z->z/(1+x*x))
δcos(x) = (cos(x), z->-z*sin(x))
δcosh(x) = (cosh(x), z->z*sinh(x))
δdivide(x::AbstractArray, y::Real) = (t = x/y; (t, z->(z/y, -sum(z.*t)/y)))
δdivide(x::Real, y::Real) = (t = x/y; (t, z->(z/y, -z*t/y)))
δdivide_const1(x, y) = (t = x/y; (t, z->(-z*t/y)))
δdivide_const2(x, y) = (x/y, z->(z/y))
δdot(x, y) = (dot(x, y), z->(z.*y, z.*x))
δdot_divide_const1{T}(x, y::T) = (sy = size(y); t = x./y; (t, z->-fit(T, z.*t./y, sy)))
δdot_divide_const2{T}(x::T, y) = (sx = size(x); (x./y, z->fit(T, z./y, sx)))
δdot_divide{TX,TY}(x::TX, y::TY) = (sx = size(x); sy = size(y); t = x./y; (t, z->(fit(TX, z./y, sx), -fit(TY, z.*t./y, sy))))
δdot_minus_const1{T}(x, y::T) = (sy = size(y); (x.-y, z->-fit(T, z, sy)))
δdot_minus_const2{T}(x::T, y) = (sx = size(x); (x.-y, z->fit(T, z, sx)))
δdot_minus{TX,TY}(x::TX, y::TY) = (sx = size(x); sy = size(y); (x.-y, z->(fit(TX, z, sx), -fit(TY, z, sy))))
δdot_plus_const1{T}(x, y::T) = (sy = size(y); (x.+y, z->fit(T, z, sy)))
δdot_plus_const2{T}(x::T, y) = (sx = size(x); (x.+y, z->fit(T, z, sx)))
δdot_plus{TX,TY}(x::TX, y::TY) = (sx = size(x); sy = size(y); (x.+y, z->(fit(TX, z, sx), fit(TY, z, sy))))
δdot_power_const1{T}(x, y::T) = (sy = size(y); t = x.^y; (t, z->fit(T, z.*t.*log.(x), sy)))
δdot_power_const2{T}(x::T, y) = (sx = size(x); t = x.^y; (t, z->fit(T, y == 2 ? z.*2x : safediv.(z.*y.*t, x), sx)))
function δdot_power{TX,TY}(x::TX, y::TY)
    sx = size(x); sy = size(y)
    t = x.^y
    t, z->(fit(TX, safediv.(z.*y.*t, x), sx), fit(TY, z.*t.*log.(x), sy))
end
δdot_times_const1{T}(x, y::T) = (sy = size(y); (x.*y, z->fit(T, z.*x, sy)))
δdot_times_const2{T}(x::T, y) = (sx = size(x); (x.*y, z->fit(T, z.*y, sx)))
δdot_times{TX,TY}(x::TX, y::TY) = (sx = size(x); sy = size(y); (x.*y, z->(fit(TX, z.*y, sx), fit(TY, z.*x, sy))))
δerf(x) = (erf(x), y->y*2/sqrt(π)*exp(-x*x))
δerfc(x) = (erfc(x), y->-y*2/sqrt(π)*exp(-x*x))
δexp(x) = (t = exp(x); (t, z->z*t))
δexpm1(x) = (t = expm1(x); (t, z->z*(1 + t)))
δfanout(x) = (x..., (z...) -> z)
δgamma(x) = (t=gamma(x); (t, y->y*polygamma(0,x)*t))
δlgamma(x) = (lgamma(x), y->y*polygamma(0,x))
δlog(x) = (log(x), z->z/x)
δlog1p(x) = (log1p(x), z->z/(1 + x))
δmax(x, y) = (max(x, y), z->(z*(x>y),z*(x<y)))
δmax_const1(x, y) = (max(x, y), z->z*(x<y))
δmax_const2(x, y) = (max(x, y), z->z*(x>y))
δmaximum(x) = (t=maximum(x); (t, y->(t.==x).*y))
δmin(x, y) = (min(x, y), z->(z*(x<y),z*(x>y)))
δmin_const1(x, y) = (min(x, y), z->z*(x>y))
δmin_const2(x, y) = (min(x, y), z->z*(x<y))
δminimum(x) = (t=minimum(x); (t, y->(t.==x).*y))
δminus(x) = (-x, z->-z)
δminus(x::AbstractArray, y::Real) = (x-y, z->(z, -sum(z)))
δminus(x::Real, y::AbstractArray) = (x-y, z->(sum(z), -z))
δminus_const1(x, y) = (x-y, z->-z)
δminus_const2(x, y) = (x-y, z->z)
δminus_const1(x, y::Real) = (x-y, z->-sum(z))
δminus_const2(x::Real, y) = (x-y, z->sum(z))
δminus{T}(x::T, y::T) = (x-y, z->(z, -z))
δmod2pi(x::Real) = (mod2pi(x), y->y)
δplus(x::AbstractArray, y::Real) = (x+y, z->(z, sum(z)))
δplus(x::Real, y::AbstractArray) = (x+y, z->(sum(z), z))
δplus_const1(x, y) = (x+y, z->z)
δplus_const2(x, y) = (x+y, z->z)
δplus_const1(x, y::Real) = (x+y, z->sum(z))
δplus_const2(x::Real, y) = (x+y, z->sum(z))
δplus{T}(x::T, y::T) = (x+y, z->(z, z))
δpower(x::Real, y::Real) = (t = x^y; (t, z->(safediv(z*y*t, x), z*t*log(x))))
δpower_const1(x, y) = (t = x^y; (t, z->z*t*log(x)))
δpower_const2(x, y) = (t = x^y; (t, z->safediv(z*y*t, x)))
δsin(x) = (sin(x), z->z*cos(x))
δsinh(x) = (sinh(x), z->z*cosh(x))
δsqrt(x) = (t = sqrt(x); (t, z->0.5*z/t))
δsum(x::AbstractArray) = (t = size(x); (sum(x), z->fill(z, t)))
δsum(x::Real) = (x, z->z)
δtan(x) = (t = tan(x); (t, z->z*(1+t*t)))
δtanh(x) = (t = tanh(x); (t, z->z*(1-t*t)))
δtimes(x::AbstractVector, y::AbstractMatrix) = (x*y, z->(vec(z*y'), x'*z))
δtimes(x::AbstractMatrix, y::AbstractVector) = (x*y, z->(z*y', vec(x'*z)))
δtimes(x::AbstractMatrix, y::AbstractMatrix) = (x*y, z->(z*y', x'*z))
δtimes(x::AbstractArray, y::Real) = (x*y, z->(z.*y, sum(z.*x)))
δtimes(x::Real, y::AbstractArray) = (x*y, z->(sum(z.*y), z.*x))
δtimes(x::Real, y::Real) = (x*y, z->(z*y, z*x))
δtimes_const1(x, y) = (x*y, z->(x'*z))
δtimes_const1(x, y::Real) = (x*y, z->sum(x.*z))
δtimes_const2(x, y) = (x*y, z->(z*y'))
δtimes_const2(x::Real, y) = (x*y, z->sum(z.*y))
δtranspose(x::AbstractVector) = (x', y->vec(y'))
δtranspose(x) = (x', y->y')
δzeros(x::AbstractArray) = zeros(x)
δzeros{T}(x::T)::T = 0.
δmean(x::Real) = (x, z->z)
δmean(x::AbstractArray) = (t = size(x); (mean(x), z->fill(z/prod(t), t)))


text:
sp = SplineParams(range(0, stop = 5, length = 10), 0, 3)
cp = ChebParams(6, 2.0, 5.0)
lp = LinParams(range(0, stop = 5, length = 10), 0)

@testset "Test original API compat" begin

    s = Spline()
    c = Cheb()
    l = Lin()

    @testset "old_name for basis family" begin
        @test CompEcon.old_name(s) == :spli
        @test CompEcon.old_name(c) == :cheb
        @test CompEcon.old_name(l) == :lin
    end

    @testset "old_name for params" begin
        @test CompEcon.old_name(sp) == :spli
        @test CompEcon.old_name(cp) == :cheb
        @test CompEcon.old_name(lp) == :lin
    end

    @testset "old_params" begin
        @test CompEcon.old_params(sp) == Any[sp.breaks, sp.evennum, sp.k]
        @test CompEcon.old_params(cp) == Any[cp.n, cp.a, cp.b]
        @test CompEcon.old_params(lp) == Any[lp.breaks, lp.evennum]
    end

    @testset "convert(Basis, Dict) and revert(Basis) --> Dict" begin
        b = Basis(sp, cp)
        d = CompEcon.revert(b)
        @test convert(Basis, d) == b
    end

end


text:
"""
# SAC.Stack

SAC.Stack is a submodule of SAC which implements simple stacking routines
to create stacks of SAC traces.

### Using

Pass an array of `SACtr` objects to the `stack` routine, with a range of times
over which to sum, and optionally times on which to align the traces, weights
to apply and the method to use.

### Available stacking methods

- `:linear` [default].  A linear (optionally weighted) stack.
"""
module Stack

import SAC

export
    stack

"List of available stacking routines"
const available_stack_methods = [:linear]

"""
    stack(S::AbstractArray{SACtr}, time_range, align=zeros(length(S)), weight; method=:linear) -> s::SACtr

Return the linear stack of all traces in `S`, aligned in time on the
value in `align`, between `times[1]` and `times[end]`.

`times` can be a Tuple, Range or Array; in all cases, only the first and last values
are used.  **N.B. If using a non-integer Range, this may not include the end value
you expect.  E.g., `0.1:0.9` by default contains only the value `0.1` because the
default step is `1`.**

`align` can be a header symbol, or an array of equal length to `S`.  If it is a Symbol,
then the values in the headers are used.

All traces must have the same sampling frequency, and an error is thrown
if not all traces can be included in the stack.

If provided, traces are weighted by the values in `weight`, which are normalised
to sum to unity.
"""
function stack(S::AbstractArray{SAC.SACtr}, time_range,
               align::Union{Symbol,AbstractArray}=zeros(SAC.SACFloat, length(S)),
               weight=ones(SAC.SACFloat, length(S));
               method=:linear)
    method in available_stack_methods ||
        throw(ArgumentError("Stacking method '$(method)' is not valid.  " *
                            "Choose from: $(available_stack_methods)"))
    # Convert to a vector of times if given as a symbol
    if typeof(align) == Symbol
        align in SAC.sac_float_hdr || error("Header '$(align)' is not a float header")
        any(S[align] .== SAC.sac_rnull) && error("Not all traces have header $align")
        align = S[align]
    end
    all(S[1].delta .== S[:delta]) || error("All traces must have same delta")
    delay = -align
    t1, t2 = time_range[1], time_range[end]
    weight /= sum(weight) # Normalise weights
    if lowercase(string(method)) == "linear"
        stack_trace = _stack_linear([w.*@view(s.t[:]) for (s, w) in zip(S, weight)],
            S[1].delta, S[:b], delay, t1, t2)
    end
    out = SAC.SACtr(S[1].delta, length(stack_trace), t1)
    out.t[:] .= stack_trace[:]
    SAC.update_headers!(out)
    out
end

"""
    _stack_linear(a, delta, b, delay, t1, t2) -> stacked_trace::Vector{eltype(a)}

Return an array containing the linear stack of the vector of vectors `a`, whose
start times are given in `b`, and where the stack is returned between times `t1` and `t2`
(which can be negative if delay is also negative).  `delta` is the sampling interval
for all the traces.

The trace starts at time `t1` and has sample spacing `delta`.
"""
function _stack_linear(a::Vector{A}, delta::Real, b::Vector,
        delay::Vector, t1::Real, t2::Real) where {A<:AbstractArray}
    N, npts, b_shift, stack_trace = _stack_array(a, delta, b, delay, t1, t2)
    # stack_trace checks that we are always in array bounds
    @inbounds for j in 1:N
        ip1 = round(Int, (t1 - b_shift[j])/delta) + 1
        @simd for i in ip1:(ip1+npts-1)
            stack_trace[i-ip1+1] += a[j][i]
        end
    end
    stack_trace # Do not normalise as weightings are applied in `stack`
end

"""
    _stack_array(a::Vector{Vector{<:Real}}, delta, b::Vector, delay::Vector, t1, t2) -> N, npts, b_shifted, stack_trace

Return the number of traces `N`, the number of points in the output stack trace,
`npts`, an array of shifted trace start times `b_shifted` and a zeroed array
for the stack itself, `stack_trace`.  This routine also checks that the input
arrays are the required length and no points are missing from the desired stack.
"""
function _stack_array(a::Vector{A}, delta, b, delay, t1, t2) where {A<:AbstractArray}
    N = length(a)
    length(b) == length(delay) == N ||
        throw(ArgumentError("Lengths of `a`, `b` and `delay` must all be the same"))
    t1 < t2 || throw(ArgumentError("`t2` is before or the same as `t1`"))
    b_shift = b .+ delay
    e_shift = b_shift .+ (length.(a) .- 1).*delta
    all(b_shift .<= t1) && all(t2 .<= e_shift) ||
        error("Not all traces have data within time range $(t1):$(t2) s
        range in b: $(extrema(b_shift))
        range in e: $(extrema(e_shift))")
    npts = round(Int, (t2 - t1)/delta) + 1
    stack_trace = zeros(eltype(eltype(a)), npts)
    N, npts, b_shift, stack_trace
end

end # module


text:
function step!(x,bkhs)
       nε = size(bkhs.D)[2]
       ε  = rand(Normal(),nε)
       x .= bkhs.C*x + bkhs.D*ε
       return x
end

function timeseries(bkhs,n)
       m  = size(bkhs.C)[1] # number of equations
       x  = zeros(m)
       ts = zeros(n,m)
       for τ in 1:n
           ts[τ,:] = step!(x,bkhs)
       end
       return ts
end

function impulse_response(bkhs,n,iε)
       m  = size(bkhs.C)[1] # number of equations
       x  = zeros(m)
       ts = zeros(n,m)
       nε = size(bkhs.D)[2]
       ε  = I(nε)[:,iε]
       x .= bkhs.C*x + bkhs.D*ε
       for τ in 1:n
           ts[τ,:] = x
           x      .= bkhs.C*x
       end
       return ts
end


text:
@testset "Logitnormal numerical moments" begin
    DN = LogitNormal(1.3)
    x = rand(DN, 1_000_000);
    m = mean(DN)
    @test abs(m - mean(x))/mean(x) <= 1e-3
    s2 = var(DN)
    @test abs(s2 - var(x))/var(x) <= 1e-2 #1e-4 too strong for random numbers
    s2m = var(DN; mean=m) # specify so that not need to recompute
    @test abs(s2m - var(x))/var(x) <= 1e-2 #1e-4 too strong for random numbers
    sd = std(DN) # specify so that not need to recompute
    @test abs(sd - std(x))/std(x) <= 1e-2 #1e-4 too strong for random numbers
    # sd2 = std(DN, mean=m) # specify so that not need to recompute
    # @test abs(sd2 - std(x))/std(x) <= 1e-2 #1e-4 too strong for random numbers
end

function test_lognormal_mode(g)
    mo = mode(g)
    @test pdf(g, mo-1e-6) < pdf(g,mo) && pdf(g, mo+1e-6) < pdf(g,mo)
end

###### numerical estimation of moments
@testset "Logitnormal seek mode" begin
    #plot(g); vline!([mode(g)])
    # median
    g = LogitNormal(); test_lognormal_mode(g)
    g = LogitNormal(0,1.4); test_lognormal_mode(g)
    # larger one of two modes, infer the lower one
    g = LogitNormal(0,1.6); test_lognormal_mode(g)
    # two modes, larger is on the right
    g = LogitNormal(0.1,1.6); test_lognormal_mode(g)
    # two modes, larger is on the left
    g = LogitNormal(-0.1,1.6); test_lognormal_mode(g)
end

@testset "fit by mode and upper quantile" begin
    #plot(g); vline!([mode(g)])
    # median
    #g = fit_mode_quantile(LogitNormal, 0.68, @qp(0.71,0.99))
    g = fit(LogitNormal, 0.68, @qp_u(0.9), Val(:mode))
    @test mode(g) ≈ 0.68
    @test quantile(g, 0.95) ≈ 0.9
    #plot(g)
end

function is_logit_slope_monotone(d, upper=0.5, lower=0.0, decreasing = false)
  x = range(lower,upper,41)[2:40]	# plotting grid
  dx = pdf.(d, x)	#density function
  if decreasing
    all(diff(dx) .<= 0)
  else
    all(diff(dx) .>= 0)
  end
end

@testset "fit by single mode and flat" begin
    d9 = d = fit_mode_flat(LogitNormal, 0.9)
    @test isapprox( mode(d), 0.9, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.9)
    @test is_logit_slope_monotone(d, 0.0, 0.9, true)
    d1 = d = fit_mode_flat(LogitNormal, 0.1)
    @test isapprox( mode(d), 0.1, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.1)
    @test is_logit_slope_monotone(d, 0.0, 0.1, true)
    @test d1.σ == d9.σ
    @test d1.μ == -d9.μ
    d5 = d = fit_mode_flat(LogitNormal, 0.5)
    @test isapprox( mode(d), 0.5, atol=1e-4)
    @test is_logit_slope_monotone(d, 0.5)
    @test is_logit_slope_monotone(d, 0.0, 0.5, true)
    #plot(d1)
end


text:
using HybridArrays
using StaticArrays
using Test

@testset "Linear algebra" begin

    @testset "HybridMatrix as a (mathematical) vector space" begin
        c = 2
        a1 = [2 4; 6 8]
        s1 = SMatrix{2,2}(a1)
        m1 = HybridMatrix{2,StaticArrays.Dynamic()}(a1)
        a2 = [4 3; 2 1]
        s2 = SMatrix{2,2}(a2)
        m2 = HybridMatrix{2,StaticArrays.Dynamic()}(a2)

        @test @inferred(c * m1)::HybridMatrix == @SMatrix [4 8; 12 16]
        @test @inferred(m1 * c)::HybridMatrix == @SMatrix [4 8; 12 16]
        @test @inferred(m1 / c)::HybridMatrix == @SMatrix [1.0 2.0; 3.0 4.0]
        @test @inferred(c \ m1)::HybridMatrix ≈ @SMatrix [1.0 2.0; 3.0 4.0]

        @test @inferred(m1 + m2) == @SMatrix [6  7;  8 9]
        @test @inferred(m1 - m2) == @SMatrix [-2  1; 4 7]

        @test @inferred(a1 + m2) == @SMatrix [6  7;  8 9]
        @test @inferred(a1 - m2) == @SMatrix [-2  1; 4 7]
        @test @inferred(m1 + a2) == @SMatrix [6  7;  8 9]
        @test @inferred(m1 - a2) == @SMatrix [-2  1; 4 7]

        @test @inferred(s1 + m2) == @SMatrix [6  7;  8 9]
        @test @inferred(s1 - m2) == @SMatrix [-2  1; 4 7]
        @test @inferred(m1 + s2) == @SMatrix [6  7;  8 9]
        @test @inferred(m1 - s2) == @SMatrix [-2  1; 4 7]

        @test vcat(m1) == m1
        @test hcat(m1) == m1
        @test @inferred(vcat(m1, m2)) == vcat(a1, a2)
        @test @inferred(hcat(m1, m2)) == hcat(a1, a2)
        @test isa(vcat(m1, m2), HybridMatrix{4,StaticArrays.Dynamic()})
        @test isa(vcat(m1, m2, m2), HybridMatrix{6,StaticArrays.Dynamic()})
        @test isa(hcat(m1, m2), HybridMatrix{2,StaticArrays.Dynamic()})
        @test isa(hcat(m1, m2, m2), HybridMatrix{2,StaticArrays.Dynamic()})
        @test isa(vcat(HybridVector{2}([1, 2]), HybridVector{3}([1, 2, 3])), HybridVector{5, Int})

        @test_throws DimensionMismatch hcat(HybridVector{2}([1, 2]), HybridVector{3}([1, 2, 3]))
        @test_throws ArgumentError vcat(m1, hcat(m1, m2))
    end
end


text:
#=
timer.jl
Utilities used in the Literate example
=#

using Statistics: median
using Base: gc_num, gc_alloc_count, GC_Diff, timev_print
using Base: cumulative_compile_time_ns_before, cumulative_compile_time_ns_after

# cf https://github.com/JuliaLang/julia/blob/6aaedecc447e3d8226d5027fb13d0c3cbfbfea2a/base/timing.jl#L327-L362

"""
    @timeo
A version of `@time` that shows the timing only, not the computed values.
Returns a string so that Literate will capture the output.
"""
macro timeo(ex) # @timeo macro to show only the timing, not results
    quote
        while false; end # compiler heuristic:
        local stats = gc_num()
#       local compile_time = cumulative_compile_time_ns_before()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
#       compile_time = cumulative_compile_time_ns_after() - compile_time
        local diff = GC_Diff(gc_num(), stats)
        local allocs = gc_alloc_count(diff)
#       timev_print(elapsedtime, diff, compile_elapsedtime)
#       time_print(elapsedtime, diff.allocd, diff.total_time, allocs, compile_time, true)
        string(
            "time=", round(elapsedtime/10^6, digits=1), # time in ms
            "ms mem=", diff.allocd,
            " alloc=", allocs,
        )
    end
end

"""
    btime()
Pretty-print the @benchmark output for non-interactive use with Literate.
Returns a string so that Literate will capture the output.
"""
btime(t) = string(
    "time=", round(median(t.times)/10^6, digits=1), # median time in ms
    "ms mem=", t.memory,
    " alloc=", t.allocs,
)