text:
# Matern 5/2 isotropic covariance function

@doc """
# Description
Constructor for the isotropic Matern 5/2 kernel (covariance)

k(x,x') = σ²(1+√5*d/ℓ + 5d²/3ℓ²)exp(-√5*d/ℓ), where d = |x-x'|
# Arguments:
* `ll::Float64`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Mat52Iso <: Isotropic
    ℓ::Float64      # Length scale
    σ2::Float64     # Signal std
    Mat52Iso(ll::Float64, lσ::Float64) = new(exp(ll), exp(2*lσ))
end

function set_params!(mat::Mat52Iso, hyp::Vector{Float64})
    length(hyp) == 2 || throw(ArgumentError("Matern 5/2 only has two parameters"))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2.0*hyp[2])
end
get_params(mat::Mat52Iso) = Float64[log(mat.ℓ), log(mat.σ2)/2.0]
get_param_names(mat::Mat52Iso) = [:ll, :lσ]
num_params(mat::Mat52Iso) = 2

metric(mat::Mat52Iso) = Euclidean()
cov(mat::Mat52Iso, r::Float64) = mat.σ2*(1+sqrt(5)*r/mat.ℓ+5*r^2/(3*mat.ℓ^2))*exp(-sqrt(5)*r/mat.ℓ)

function grad_kern(mat::Mat52Iso, x::Vector{Float64}, y::Vector{Float64})
    r = distance(mat,x,y)
    exp_r = exp(-sqrt(5)*r/mat.ℓ)

    g1 = mat.σ2*(5*r^2/mat.ℓ^2)*((1+sqrt(5)*r/mat.ℓ)/3)*exp_r   #dK_d(log ℓ)
    g2 = 2.0*mat.σ2*(1+sqrt(5)*r/mat.ℓ+5*r^2/(3*mat.ℓ^2))*exp_r #dK_d(log σ)
    return [g1,g2]
end

function grad_stack!(stack::AbstractArray, mat::Mat52Iso, X::Matrix{Float64}, data::IsotropicData)
    nobsv = size(X,2)
    R = distance(mat, X, data)
    exp_R = exp(-sqrt(5)*R/mat.ℓ)

    for i in 1:nobsv, j in 1:i
        @inbounds stack[i,j,1] = mat.σ2*(5*R[i,j]^2/mat.ℓ^2)*((1+sqrt(5)*R[i,j]/mat.ℓ)/3)*exp_R[i,j]      # dK_dℓ
        @inbounds stack[j,i,1] = stack[i,j,1]
        @inbounds stack[i,j,2] = 2.0*mat.σ2*(1+sqrt(5)*R[i,j]/mat.ℓ+5*R[i,j]^2/(3*mat.ℓ^2))*exp_R[i,j]   # dK_dσ
        @inbounds stack[j,i,2] = stack[i,j,2]
    end

    return stack
end


text:
using MyUtils
using Test

@testset "MyUtils.jl" begin
    # Write your tests here.
end

ext: jl


text:
# ---
# title: 982. Triples with Bitwise AND Equal To Zero
# id: problem982
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/triples-with-bitwise-and-equal-to-zero/description/>
# hidden: true
# ---
#
# Given an array of integers `A`, find the number of triples of indices (i, j,
# k) such that:
#
#   * `0 <= i < A.length`
#   * `0 <= j < A.length`
#   * `0 <= k < A.length`
#   * `A[i] & A[j] & A[k] == 0`, where `&` represents the bitwise-AND operator.
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,1,3]
#     Output: 12
#     Explanation: We could choose the following i, j, k triples:
#     (i=0, j=0, k=1) : 2 & 2 & 1
#     (i=0, j=1, k=0) : 2 & 1 & 2
#     (i=0, j=1, k=1) : 2 & 1 & 1
#     (i=0, j=1, k=2) : 2 & 1 & 3
#     (i=0, j=2, k=1) : 2 & 3 & 1
#     (i=1, j=0, k=0) : 1 & 2 & 2
#     (i=1, j=0, k=1) : 1 & 2 & 1
#     (i=1, j=0, k=2) : 1 & 2 & 3
#     (i=1, j=1, k=0) : 1 & 1 & 2
#     (i=1, j=2, k=0) : 1 & 3 & 2
#     (i=2, j=0, k=1) : 3 & 2 & 1
#     (i=2, j=1, k=0) : 3 & 1 & 2
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 1000`
#   2. `0 <= A[i] < 2^16`
#
#
## @lc code=start
using LeetCode

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


text:
using Parameters

function trankplot(model::SampleModel, param::Symbol; bins=40)
    nt = read_samples(model, :namedtuple)
    n_draws = model.method.num_samples
    n_chains = model.n_chains[1]

    if n_chains == 1
        @warn "trankplot() needs multiple chains."
        return
    end

    if !(param in keys(nt))
        @warn "Parameter $(param) not in model."
        return
    end

   if length(size(nt[param])) > 1
        params = [Symbol("$(param)[$i]") for i in 1:size(nt[param], 1)]
        n_pars = length(params)
    else
        params = [param]
        n_pars = 1
    end

    dfs = read_summary(model)
    n_eff = dfs[indexin(params, dfs[:, :parameters]), :ess]

    trankplot(nt, param; bins=bins, n_draws=n_draws, n_chains=n_chains, n_eff=n_eff)
end


text:
"""
$(TYPEDEF)

Binned point list structure allowing for fast
check for already existing points.

This provides better performance for indendifying
already inserted points than the naive linear search.

OTOH the implementation is still quite naive - it dynamically maintains
a cuboid bining region with a fixed number of bins.

Probably tree based adaptive methods (a la octree) will be more efficient,
however they will be harder to implement.

In an ideal world, we would maintain a dynamic Delaunay triangulation, which
at once could be the starting point of mesh generation which will follow here
anyway.
"""
mutable struct BinnedPointList{T}

    # Space dimension
    dim::Int32

    # Point distance tolerance. Points closer than tol
    # (in Euclidean distance) will be identified, i.e.
    # are collapsed to the first inserted.
    tol::T

    # The union of all bins is the binning region -
    # a cuboid given by two of its corners. It is calculated
    # dynamically depending on the inserted points.
    binning_region_min::Vector{T}
    binning_region_max::Vector{T}

    # Increase factor of binning region (with respect to
    # the cuboid defined by the coordinates of the binned points)
    binning_region_increase_factor::T

    # The actual point list
    points::ElasticArray{T,2}

    # The bins are vectors of indices of points in the point list
    # We store them in a dim-dimensional array  of length number_of_directional_bins^dim
    bins::Array{Vector{Int32}}

    # Number of bins in each space dimension
    number_of_directional_bins::Int32

    # Some points will fall outside of the binning region.
    # We collect them in vector of ubinned point indices
    unbinned::Vector{Int32}

    # Number of unbinned  points tolerated without rebinning
    num_allowed_unbinned_points::Int32

    # Maximum ratio of unbinned  points  in point list
    max_unbinned_ratio::T

    # Storage of current point bin
    current_bin::Vector{Int32}

    BinnedPointList{T}(::Nothing) where T=new()
end

"""
    $(SIGNATURES)

Create and initialize binned point list
"""
function BinnedPointList(::Type{T},dim;
                         tol=1.0e-12,
                         number_of_directional_bins=10,
                         binning_region_increase_factor=0.01,
                         num_allowed_unbinned_points=5,
                         max_unbinned_ratio=0.05
                         ) where T
    bpl=BinnedPointList{T}(nothing)

    bpl.dim=dim
    bpl.tol=tol
    bpl.number_of_directional_bins=number_of_directional_bins
    bpl.binning_region_increase_factor=binning_region_increase_factor
    bpl.num_allowed_unbinned_points=num_allowed_unbinned_points
    bpl.max_unbinned_ratio=max_unbinned_ratio

    bpl.points=ElasticArray{Cdouble}(undef,dim,0)

    bpl.binning_region_min=fill(floatmax(T),dim)
    bpl.binning_region_max=fill(-floatmax(T),dim)

    bpl.unbinned=Vector{Int32}(undef,0)
    bpl.bins=Array{Vector{Int32},dim}(undef,zeros(Int32,dim)...)

    bpl.current_bin=zeros(Int32,bpl.dim)
    bpl
end

"""
$(SIGNATURES)

Create and initialize binned point list
"""
BinnedPointList(dim; kwargs...)=BinnedPointList(Cdouble,dim; kwargs...)

#
# Find point in index list (by linear search)
# Return its index, or zero if not found
#
function _findpoint(bpl,index,p)
    for i=1:length(index)
        @views if norm(bpl.points[:,index[i]]-p)<bpl.tol
            return index[i]
        end
    end
    return 0
end

#
# Calculate the bin of the point. Result is stored
# in bpl.current_bin
#
function _bin_of_point!(bpl,p)
    for idim=1:bpl.dim
        # scaled value for particular dimension
        s=(p[idim]-bpl.binning_region_min[idim])/(bpl.binning_region_max[idim]-bpl.binning_region_min[idim])
        if s>1 || s<0
            # point lies outside of binning area
            bpl.current_bin[idim]=0
        else
            # calculate bin in particular dimension
            bpl.current_bin[idim]=ceil(Int32,bpl.number_of_directional_bins*s)
        end
    end
end

#
# Re-calculate binning if there are too many unbinned points
# This amounts to two steps:
# - Enlarge binning area in order to include all points
# - Re-calculate all point bins
function _rebin_all_points!(bpl)
    if length(bpl.unbinned)>max(bpl.num_allowed_unbinned_points,
                                bpl.max_unbinned_ratio*size(bpl.points,2))

        # Calculate extrema of unbinned points
        @views e=extrema(bpl.points[:,bpl.unbinned],dims=2)

        for i=1:bpl.dim
            # Increase binning region according to unbinned extrema
            bpl.binning_region_min[i]=min(bpl.binning_region_min[i],e[i][1])
            bpl.binning_region_max[i]=max(bpl.binning_region_max[i],e[i][2])

            # Slightly increase binning region further in order to
            # include all existing points with tolerance
            delta=max(bpl.binning_region_max[i]-bpl.binning_region_min[i],bpl.tol)
            bpl.binning_region_min[i]-=bpl.binning_region_increase_factor*delta
            bpl.binning_region_max[i]+=bpl.binning_region_increase_factor*delta
        end

        # Re-allocate all bins
        if bpl.dim==1
            bpl.bins=[zeros(Int32,0)  for i=1:bpl.number_of_directional_bins]
        elseif bpl.dim==2
            bpl.bins=[zeros(Int32,0)  for i=1:bpl.number_of_directional_bins,j=1:bpl.number_of_directional_bins]
        elseif bpl.dim==3
            bpl.bins=[zeros(Int32,0)  for i=1:bpl.number_of_directional_bins,j=1:bpl.number_of_directional_bins,k=1:bpl.number_of_directional_bins]
        end

        # Register all points in their respctive bins
        for i=1:size(bpl.points,2)
            @views  _bin_of_point!(bpl,bpl.points[:,i])
            push!(bpl.bins[bpl.current_bin...],i)
        end

        # Re-allocate unbinned index
        bpl.unbinned=zeros(Int32,0)
    end
end

"""
$(SIGNATURES)

If another point with distance less the tol from p is
in pointlist, return its index. Otherwise, insert point into pointlist.
"""
function Base.insert!(bpl::BinnedPointList{T},p) where T
    _rebin_all_points!(bpl)
    _bin_of_point!(bpl,p)
    if reduce(*,bpl.current_bin)>0
        i=_findpoint(bpl,bpl.bins[bpl.current_bin...],p)
        if i>0
            return i
        else
            append!(bpl.points,p)
            i=size(bpl.points,2)
            push!(bpl.bins[bpl.current_bin...],i)
            return i
        end
    else
        i=_findpoint(bpl,bpl.unbinned,p)
        if i>0
            return i
        else
            append!(bpl.points,p)
            i=size(bpl.points,2)
            push!(bpl.unbinned,i)
            return i
        end
    end
end

Base.insert!(bpl::BinnedPointList{T},x::Number) where T =insert!(bpl,(x))
Base.insert!(bpl::BinnedPointList{T},x::Number,y::Number) where T =insert!(bpl,(x,y))
Base.insert!(bpl::BinnedPointList{T},x::Number,y::Number,z::Number) where T =insert!(bpl,(x,y,z))

"""
$(SIGNATURES)

Return the array of points in the point list.
"""
points(bpl::BinnedPointList{T}) where T = bpl.points

# Just for being able to check of all of the above was worth the effort...
function naiveinsert!(bpl::BinnedPointList{T},p) where T
    for i=1:size(bpl.points,2)
        @views if norm(bpl.points[:,i]-p)<bpl.tol
            return i
        end
    end
    append!(bpl.points,p)
    size(bpl.points,2)
end


text:
"""
General Interface for Pricing

		Delta=delta(mcProcess,rfCurve,mcBaseData,payoff_)

Where:\n
		mcProcess          = Process to be simulated.
		rfCurve  = Zero Rate Data.
		mcBaseData = Basic properties of MonteCarlo simulation
		payoff_ = Payoff(s) to be priced

		Delta     = Delta of the derivative

"""
function delta(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoff::AbstractPayoff, dS0::Real = 1e-7)
    Price = pricer(mcProcess, rfCurve, mcConfig, abstractPayoff)
    mcProcess_up = deepcopy(mcProcess)
    mcProcess_up.underlying.S0 += dS0
    set_seed(mcConfig)
    PriceUp = pricer(mcProcess_up, rfCurve, mcConfig, abstractPayoff)
    Delta = (PriceUp - Price) / dS0

    return Delta
end

function delta(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoffs::Array{abstractPayoff}, dS0::Real = 1e-7) where {abstractPayoff <: AbstractPayoff}
    Prices = pricer(mcProcess, rfCurve, mcConfig, abstractPayoffs)
    mcProcess_up = deepcopy(mcProcess)
    mcProcess_up.underlying.S0 += dS0
    PricesUp = pricer(mcProcess_up, rfCurve, mcConfig, abstractPayoffs)
    Delta = (PricesUp .- Prices) ./ dS0

    return Delta
end

function delta(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoffs::Dict{AbstractPayoff, Number}, dS0::Real = 1e-7)
    Prices = pricer(mcProcess, rfCurve, mcConfig, abstractPayoffs)
    mcProcess_up = deepcopy(mcProcess)
    mcProcess_up.underlying.S0 += dS0
    set_seed(mcConfig)
    PricesUp = pricer(mcProcess_up, rfCurve, mcConfig, abstractPayoffs)
    Delta = (PricesUp .- Prices) ./ dS0

    return Delta
end

function delta(mcProcess::Dict{String, AbstractMonteCarloProcess}, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, dict_::Dict{String, Dict{AbstractPayoff, Number}}, underl_::String, dS::Real = 1e-7)
    if (!isnothing(findfirst("_", underl_)))
        error("deltas are defined on single name")
    end
    set_seed(mcConfig)
    underlyings_models = keys(mcProcess)
    underlyings_payoff = keys(dict_)
    price0 = pricer(mcProcess, rfCurve, mcConfig, dict_)
    price2 = 0.0
    set_seed(mcConfig)
    mcProcess_up = deepcopy(mcProcess)
    delta_ = 0.0
    keys_mkt = collect(keys(mcProcess_up))
    idx_supp = 0
    idx_1 = findfirst(y -> y == underl_, keys_mkt)
    if (isnothing(idx_1))
        idx_1 = findfirst(out_el -> any(out_el_sub -> out_el_sub == underl_, split(out_el, "_")), keys_mkt)
        if (isnothing(idx_1))
            return 0.0
        end
        multi_name = keys_mkt[idx_1]
        idx_supp = findfirst(x_ -> x_ == underl_, split(multi_name, "_"))
        model_ = mcProcess_up[multi_name]
        model_.models[idx_supp].underlying.S0 += dS
        delete!(mcProcess_up, multi_name)
        tmp_mkt = multi_name → model_
        mcProcess_up = mcProcess_up + tmp_mkt
        price2 = pricer(mcProcess_up, rfCurve, mcConfig, dict_)
    else
        model = mcProcess_up[keys_mkt[idx_1]]
        model.underlying.S0 += dS
        delete!(mcProcess_up, keys_mkt[idx_1])
        tmp_mkt = keys_mkt[idx_1] → model
        mcProcess_up = mcProcess_up + tmp_mkt
        price2 = pricer(mcProcess_up, rfCurve, mcConfig, dict_)
    end

    return (price2 - price0) / dS
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis

# Parameters
ex = [1.;0.;0.]

l1 = 1.0
l2 = sqrt(2) / 2
x, y = .1, .1

vert11 = [0.;0.;l1 / 2]
vert12 = -vert11

vert21 = [0.;0.;l2 / 2]

# Initial orientation
phi1 = pi / 4
q1 = UnitQuaternion(RotX(phi1))

# Links
origin = Origin{Float64}()
link1 = Box(x, y, l1, l1, color = RGBA(1., 1., 0.))
link2 = Box(x, y, l2, l2, color = RGBA(1., 1., 0.))

# Constraints
socket0to1 = EqualityConstraint(Spherical(origin, link1; p2=vert11))
socket1to2 = EqualityConstraint(Spherical(link1, link2; p1=vert12, p2=vert21))

links = [link1;link2]
constraints = [socket0to1;socket1to2]

mech = Mechanism(origin, links, constraints)
setPosition!(origin,link1,p2 = vert11,Δq = q1)
setPosition!(link1,link2,p1 = vert12,p2 = vert21,Δq = inv(q1)*UnitQuaternion(RotY(0.2)))

storage = simulate!(mech, 10., record = true)
visualize(mech, storage)


text:
# Plot and save the trained params
params = hcat(nodes[2:end], p[:,2:end])
CSV.write(string(ckpt_path,"/params.csv"), DataFrame(params), header=nodes)
if !parsed_args["disable-display"]
    if ispath(string(fig_path, "/p_inference_iter_tracking.png"))
        rm(string(fig_path, "/p_inference_iter_tracking.png"))
    end
    cbp(p, iter)
end

if "data" in keys(conf)
    cbrd(p)
end

# Plot for each conditions
if !parsed_args["disable-display"]
   for i in n_exp-2:n_exp
      plt = plot(size=(500, 300), legend=:right)
      cbi!(p, i, plt)
   end
end

# Plot comparison heatmaps
if !parsed_args["disable-display"]
    l = @layout [a{0.457w} b{0.457w}]
    pyplot_ticks = [string("#", g, "#") for g in nodes[2:end]]
    heatmap1 = heatmap(p_gold[:, 2:end], legend=false, size=(320, 320),
                       c=:ice, xtickfontrotation=90)
    title!("Ground truth")
    xticks!(1:length(nodes)-1, pyplot_ticks)
    yticks!(1:length(nodes)-1, pyplot_ticks)
    heatmap2 = heatmap(p[:, 2:end], size=(360, 320), c=:ice,
                       clim=(-1, 1), xtickfontrotation=90)
    title!("Training result")
    xticks!(1:length(nodes)-1, pyplot_ticks)
    yticks!(1:length(nodes)-1, pyplot_ticks)
    heatmaps = plot(heatmap1, heatmap2, layout = l, size=(700, 320))
    png(heatmaps, string(fig_path, "/heatmap.png"))
end

# Joint figure
if !parsed_args["disable-display"]
    pick_idx = n_exp
    plt = plot(size=(360, 320), legend=:right, leftmargin=40px, rightmargin=40px)
    cbi!(p, pick_idx, plt)
    title!("ODE trajectory")
    l = @layout [a{0.3w} b{0.3w} c{0.4w}]
    merged = plot(heatmap1, heatmap2, plt, layout = l, size=(1080, 320))
    png(merged, string(fig_path, "/merged.png"))
end


text:
using StatisticalRethinking, CmdStan
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df);

first(df2, 5)

struct Post
  mu::Float64
  sigma::Float64
  ll:: Float64
  prod::Float64
  prob::Float64
end

mu_list = repeat(range(140, 160, length=200), 200);
sigma_list = repeat(range(4, 9, length=200), inner=200);

ll = zeros(40000);
for i in 1:40000
    d1 = Normal(mu_list[i], sigma_list[i])
    ll[i] = sum(log.(pdf.(d1, df2[:height])))
end

d2 = Normal(178.0, 20.0)
d3 = Uniform(0, 50)
prod = ll + log.(pdf.(d2, mu_list)) + log.(pdf.(d3, sigma_list))
prob = exp.(prod .- maximum(prod))
post = DataFrame(mu=mu_list, sigma=sigma_list, ll=ll, prod=prod, prob=prob)
first(post, 10)

samples = post[sample(1:size(post, 1), Weights(post[:prob]), 10000, replace=true), :]

density(samples[:mu])

density(samples[:sigma])

#hpd(samples[:mu])

#hpd(samples[:sigma])

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# Types
include("QuadForm/Types.jl")
include("QuadForm/Quad/Types.jl")
include("QuadForm/Herm/Types.jl")

# Binary quadratic forms

include("QuadForm/QuadBin.jl")

# Basic
include("QuadForm/Basic.jl")
include("QuadForm/Lattices.jl")

include("QuadForm/Misc.jl")

# Quadratic
include("QuadForm/Quad/Spaces.jl")
include("QuadForm/Quad/Genus.jl")
include("QuadForm/Quad/GenusRep.jl")

# Hermitian
include("QuadForm/Herm/Spaces.jl")
include("QuadForm/Herm/Genus.jl")
include("QuadForm/Herm/GenusRep.jl")
include("QuadForm/Herm/Mass.jl")

include("QuadForm/Morphism.jl")
include("QuadForm/Database.jl")
include("QuadForm/Enumeration.jl")
include("QuadForm/LineOrbits.jl")
include("QuadForm/MassQuad.jl")


text:
const STATLIST = [:n, :posn, :mean, :var, :geom, :logvar, :sd, :se, :median, :min, :max, :range, :q1, :q3, :iqr, :kurt, :skew, :harmmean, :ses, :sek, :sum]

#=
function sortbyvec!(a, vec)
    sort!(a, by = x -> findfirst(y -> x == y, vec))
end
=#
#=
ispositive(::Missing) = false
ispositive(x::AbstractFloat) = isnan(x) ? false : x > zero(x)
ispositive(x) = x > zero(x)
################################################################################
struct SkipNonPositive{T}
    x::T
end
skipnonpositive(itr) = SkipNonPositive(itr)

Base.IteratorEltype(::Type{SkipNonPositive{T}}) where {T} = Base.IteratorEltype(T)
#Base.eltype(::Type{SkipNonPositive{T}}) where {T} = nonmissingtype(eltype(T))
function Base.iterate(itr::SkipNonPositive, state...)
    y = iterate(itr.x, state...)
    y === nothing && return nothing
    item, state = y
    while !ispositive(item)
        y = iterate(itr.x, state)
        y === nothing && return nothing
        item, state = y
    end
    item, state
end
Base.IndexStyle(::Type{<:SkipNonPositive{T}}) where {T} = IndexStyle(T)
Base.eachindex(itr::SkipNonPositive) =
    Iterators.filter(i -> ispositive(@inbounds(itr.x[i])), eachindex(itr.x))
Base.keys(itr::SkipNonPositive) =
    Iterators.filter(i -> ispositive(@inbounds(itr.x[i])), keys(itr.x))
Base.@propagate_inbounds function getindex(itr::SkipNonPositive, I...)
    v = itr.x[I...]
    !ispositive(v) && throw(ErrorException("the value at index $I is non positive"))
    v
end
function Base.length(itr::SkipNonPositive)
    n = 0
    for i in itr n+=1 end
    n
end
################################################################################
struct SkipNaNorMissing{T}
    x::T
end
skipnanormissing(itr) = SkipNaNorMissing(itr)

Base.IteratorEltype(::Type{SkipNaNorMissing{T}}) where {T} = Base.IteratorEltype(T)
#Base.eltype(::Type{SkipNaNorMissing{T}}) where {T} = nonmissingtype(eltype(T))
function Base.iterate(itr::SkipNaNorMissing, state...)
    y = iterate(itr.x, state...)
    y === nothing && return nothing
    item, state = y
    while isnanormissing(item)
        y = iterate(itr.x, state)
        y === nothing && return nothing
        item, state = y
    end
    item, state
end
Base.IndexStyle(::Type{<:SkipNaNorMissing{T}}) where {T} = IndexStyle(T)
Base.eachindex(itr::SkipNaNorMissing) =
    Iterators.filter(i -> isnanormissing(@inbounds(itr.x[i])), eachindex(itr.x))
Base.keys(itr::SkipNaNorMissing) =
    Iterators.filter(i -> isnanormissing(@inbounds(itr.x[i])), keys(itr.x))
Base.@propagate_inbounds function getindex(itr::SkipNaNorMissing, I...)
    v = itr.x[I...]
    !isnanormissing(v) && throw(ErrorException("The value at index $I is NaN or missing!"))
    v
end
function Base.length(itr::SkipNaNorMissing)
    n = 0
    for i in itr n+=1 end
    n
end
=#
################################################################################
length2(x) = length(x)
function length2(itr::Base.SkipMissing)
    n = 0
    for i in itr n+=1 end
    n
end
################################################################################
"""
    dataimport(data; vars, sort = nothing)

Import data.

* data - tabular data;
* vars - variables;
* sort - sort by categories.
"""
function dataimport(data; vars, sort = nothing)
    if isa(vars, Symbol) vars = [vars] end

    if !isnothing(sort) && isa(sort, Symbol) sort = [sort] end

    dataimport_(data, vars, sort)
end

function dataimport_(data, vars::AbstractVector, sort::AbstractVector)
    cols   = Tables.columns(data)
    cdata  = Tuple(Tables.getcolumn(cols, y) for y in sort)
    d      = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
    indsdict!(d, cdata)
    varlabel = :Variable
    while varlabel in sort
        varlabel = Symbol(varlabel, '_')
    end

    sdata = Vector{ObsData}(undef, length(d)*length(vars))

    i = one(Int)
    @inbounds for (k, v) in d
        for c in vars
            dict = Dict{Symbol, Any}(sort .=> k)
            dict[varlabel] = c
            sdata[i] = ObsData(view(Tables.getcolumn(data, c), v), varlabel,  dict)
            i += one(Int)
        end
    end
    return DataSet(identity.(sdata))
end
function dataimport_(data, vars, sort::Nothing)
    sdata = Vector{ObsData}(undef, length(vars))
    for i in 1:length(vars)
        sdata[i] = ObsData(Tables.getcolumn(data, vars[i]), vars[i],  Dict(:Variable=>vars[i]))
    end
    DataSet(identity.(sdata))
end
"""
    descriptives(data, vars, sort = nothing; kwargs...)
"""
function descriptives(data, vars = nothing, sort = nothing; kwargs...)
    if isa(vars, String) vars = [Symbol(vars)] end
    if isa(vars, Symbol) vars = [vars] end
    if isa(sort, String) sort = [Symbol(sort)] end
    if isa(sort, Symbol) sort = [sort] end

    descriptives(dataimport_(data, vars, sort); kwargs...)
end

function descriptives(data; vars = nothing, sort = nothing, kwargs...)
    if isnothing(vars)
        vars = Vector{Symbol}(undef, 0)
        coln   = Tables.columnnames(data)
        for i in coln
            if eltype(Tables.getcolumn(data, i)) <: Number push!(vars, i) end
        end
        if length(vars) == 0 error("No column found for descriptive statistics!") end
    end
    descriptives(data, vars, sort; kwargs...)
end
"""
    descriptives(data::DataSet{T}; kwargs...) where T <: ObsData

* kwargs:
* `skipmissing`
* `skipnonpositive`
"""
function descriptives(data::DataSet{T}; kwargs...) where T <: ObsData
    kwargs = Dict{Symbol, Any}(kwargs)
    k = keys(kwargs)

    if !(:corrected in k)
        kwargs[:corrected] = true
    end
    if !(:skipmissing in k)
        kwargs[:skipmissing] = false
    end
    if !(:skipnonpositive in k)
        kwargs[:skipnonpositive] = false
    end

    if !(:stats in k)
        kwargs[:stats] = [:n, :mean, :sd, :se, :median, :min, :max]
    else
        if isa(kwargs[:stats], Symbol)  kwargs[:stats] = [kwargs[:stats]] end
        if isa(kwargs[:stats], String)  kwargs[:stats] = [Symbol(kwargs[:stats])] end
    end

     kwargs[:stats] ⊆ STATLIST || error("Some statistics not known!")

    ds = Vector{Descriptives}(undef, length(data))
    i  = 1
    for d in data
        # skipmissing
        if kwargs[:skipmissing]
            logvec = vec = skipnanormissing(d.obs)
        else
            vec = d.obs
        end
        n_  = logn_ = length2(vec)
        # skipnonpositive
        if kwargs[:skipnonpositive]
            logvec = skipnonpositive(d.obs)
            logn_  = length2(logvec)
        elseif !kwargs[:skipmissing]
            logvec = d.obs
        end

        result = Dict{Symbol, Float64}()

        for s in kwargs[:stats]
            if !(n_ > 0)
                result[s] = NaN
                continue
            end
            if s == :n
                result[s] = n_
            elseif s == :posn
                result[s] = logn_
            elseif s == :mean
                result[s] = sum(vec) / n_
            elseif s == :sd
                Base.ht_keyindex(result, :mean) > 0 || begin result[:mean] =  sum(vec) / n_ end
                result[s] = std(vec; corrected = kwargs[:corrected], mean = result[:mean])
            elseif s == :var
                Base.ht_keyindex(result, :mean) > 0 || begin result[:mean] = sum(vec) / n_ end
                result[s] = var(vec; corrected = kwargs[:corrected], mean = result[:mean])
            elseif s == :se
                Base.ht_keyindex(result, :mean) > 0 || begin result[:mean] = sum(vec) / n_ end
                Base.ht_keyindex(result, :sd) > 0 || begin result[:sd] = std(vec; corrected = kwargs[:corrected], mean = result[:mean]) end
                result[s] = result[:sd] / sqrt(n_)
            elseif s == :median
                result[:median] = median(vec)
            elseif s == :min
                result[s] = minimum(vec)
            elseif s == :max
                result[s] = maximum(vec)
            elseif s == :q1
                result[:q1] = quantile(vec, 0.25)
            elseif s == :q3
                result[:q3] = quantile(vec, 0.75)
            elseif s == :iqr
                result[s] = abs(quantile(vec, 0.75) - quantile(vec, 0.25))
            elseif s == :range
                result[s] = abs(maximum(vec) - minimum(vec))
            elseif s == :kurt
                result[s] = kurtosis_(vec, sum(vec) / n_)
            elseif s == :skew
                result[s] = skewness_(vec, sum(vec) / n_)
            elseif s == :harmmean
                result[s] = harmmean(vec)
            elseif s == :ses
                result[s] = sesvec(vec)
            elseif s == :sek
                result[s] = sekvec(vec)
            elseif s == :sum
                result[s] = sum(vec)
            end
            if !(logn_ > 0)
                result[s] = NaN
                continue
            end
            if s == :geom
                result[:logmean] = sum(log, logvec) / logn_
                result[:geom] = exp(result[:logmean])
            elseif s == :logvar
                Base.ht_keyindex(result, :logmean) > 0 || begin result[:logmean] = sum(log, logvec) / logn_ end
                result[:logvar] = var(logvec; corrected = kwargs[:corrected], mean = result[:logmean])
            end
        end
        filter!(x -> x.first in kwargs[:stats], result)
        ds[i] =  Descriptives(d, result)
        i += 1
    end
    DataSet(identity.(ds))
end

function sesvec(vec)
    ses(length(vec))
end
function ses(n)
    return sqrt(6 * n *(n - 1) / ((n - 2) * (n + 1) * (n + 3)))
end

function sekvec(vec)
    n   = length(vec)
    sek(n; ses = ses(n))
end
function sek(n; ses = ses(n))
    return 2 * ses * sqrt((n * n - 1)/((n - 3) * (n + 5)))
end

function kurtosis_(v, m)
    n = length(v)
    cm2 = 0.0  # empirical 2nd centered moment (variance)
    cm4 = 0.0  # empirical 4th centered moment
    for i in v
        z = i - m
        z2 = z * z
        cm2 += z2
        cm4 += z2 * z2
    end
    cm4 /= n
    cm2 /= n
    return (cm4 / (cm2 * cm2)) - 3.0
end

function skewness_(v, m)
    n = length(v)
    cm2 = 0.0   # empirical 2nd centered moment (variance)
    cm3 = 0.0   # empirical 3rd centered moment
    for i in v
        z = i - m
        z2 = z * z
        cm2 += z2
        cm3 += z2 * z
    end
    cm3 /= n
    cm2 /= n
    return cm3 / sqrt(cm2 * cm2 * cm2)  # this is much faster than cm2^1.5
end

################################################################################
#
################################################################################

function MetidaBase.metida_table(obj::DataSet{DS}; sort = STATLIST, stats = nothing, id = nothing) where DS <: Descriptives
    idset  = Set(keys(first(obj).data.id))
    resset = Set(keys(first(obj).result))
    if length(obj) > 1
        for i = 2:length(obj)
            union!(idset,  Set(keys(obj[i].data.id)))
            union!(resset, Set(keys(obj[i].result)))
        end
    end
    if !isnothing(stats)
        stats ⊆ STATLIST || error("Some statistics not known!")
        if isa(stats, Symbol) stats = [stats] end
        ressetl = sortbyvec!(collect(intersect(resset, stats)), sort)
    else
        ressetl = sortbyvec!(collect(resset), sort)
    end
    if !isnothing(id)
        if isa(id, Symbol) id = [id] end
        id ⊆ idset || error("Some id not in dataset!")
        idset = intersect(idset, id)
    end
    mt1 = metida_table((getid(obj, :, c) for c in idset)...; names = idset)
    mt2 = metida_table((obj[:, c] for c in ressetl)...; names = ressetl)
    MetidaTable(merge(mt1.table, mt2.table))
end

################################################################################
# SHOW
################################################################################

function Base.show(io::IO, obj::DataSet{DS}) where DS <: Descriptives
    show(io, metida_table(obj))
end

function Base.show(io::IO, obj::DataSet{OD}) where OD <: ObsData
    println(io, "DataSet: observations")
    for i in obj
        println(io, "  Var: $(i.var); ID: $(i.id); Length: $(length(i.obs))")
    end
end

################################################################################
# To MetidaBase

function cvfromsd(σ)
    return sqrt(exp(σ^2)-1)
end


text:
module AttractorNetworksBase
using LinearAlgebra
using OrdinaryDiffEq
using Distributions, Random
# numerical defaults for the constructors
const default_values = Dict(
    :muw=>3.0,
    :Wmeanf => [ 2.5  -1.3
                 2.4   -1.  ],
    :Wmeanf_alt => [ 1.08  -0.8
                       0.8   -2.25 ],
    :external_current => 7.0,
    :tau_ei => (0.02,0.01),
    :gainalpha => 0.02 )

abstract type IOFunction{R} end
Base.Broadcast.broadcastable(g::IOFunction)=Ref(g)
Base.copy(g::IOFunction) = g

struct IOQuad{T} <: IOFunction{T}
    α::T
end
@inline function (gf::IOQuad{T})(x::T) where T
     return gf.α*(log(1.0+exp(x)))^2
end
@inline function ioprime(x::T,gf::IOQuad{T}) where T
    ex= exp(x)
    return gf.α * 2.0ex*log(1. + ex)/(1. + ex)
end
@inline function ioprimeprime(x::T,gf::IOQuad{T}) where T
    ex= exp(x)
    oex= 1. + ex
    return 2.0*gf.α*ex*(log(oex)+ex) / (oex*oex)
end
@inline function ioinv(y::T,gf::IOQuad{T}) where T
    return log(exp( sqrt(y/gf.α) )-1.0)
end

struct IOId{T} <:IOFunction{T}
end
@inline function (gf::IOId{T})(x::T) where T
     return x
end
@inline function ioprime(x::T,gf::IOId{T}) where T
    return one(T)
end
@inline function ioprimeprime(x::T,gf::IOId{T}) where T
    return zero(T)
end
@inline function ioinv(x::T,gf::IOId{T}) where T
    return x
end

# struct IORelu{T} <:IOFunction{T}
# end
#

"""
        populations(ne::I,ni::I) where I<:Integer -> (pe::BitArray,pi::BitArray)

For `ne` excitatory an `ni` inhibitory neurons, returns the positions of each in the form of BitArrays. The convention is that E neurons come first.  The indexes of the connections in `W` *from* poulation `pa` *to* population `pb` can then be selected as: `view(W,pb,pa)`
"""
function populations(ne::I,ni::I) where I<:Integer
    ntot = ne+ni
    pe = BitArray(falses(ntot))
    pi = copy(pe)
    pe[1:ne] .= true
    pi[ne+1:end] .= true
    return (pe,pi)
end
"""
        populations(ne::I,ni::I, ne_aux::I,ni_aux::I) where I<:Integer
            -> ((pe,pi,pe_pri,pe_aux,pi_pri,pi_aux)

Like `populations(ne,ni)`, but now E and I populations can have a subset of auxiliary neurons, with `ne = ne_pri+ne_aux`. Returns the indexing for each population in the form of BitArrays. The convetion is that excitatory neurons go first, and principal neurons go first.
"""
function populations(ne::I,ni::I, ne_aux::I,ni_aux::I) where I<:Integer
    ne_pri = ne - ne_aux
    ni_pri = ni - ni_aux
    @assert (ne_pri >= 0 ) && (ni_pri >= 0)
    ntot = ne+ni
    pe,pi = populations(ne,ni)
    pe_pri,pe_aux,pi_pri,pi_aux = (BitArray(falses(ntot)) for _ in 1:4)
    pe_pri[1:ne_pri] .= true
    pe_aux[ne_pri+1:ne_pri+ne_aux] .= true
    pi_pri[ne+1:ne+ni_pri] .= true
    pi_aux[ne+ni_pri+1:end] .= true
    return (pe,pi,pe_pri,pe_aux,pi_pri,pi_aux)
end

"""
        diagtozero!(M::AbstractMatrix{T}) where T
Replaces the matrix diagonal with zeros
"""
function diagtozero!(M::AbstractMatrix{T}) where T<:Real
    ms = minimum(size(M))
    for i in 1:ms
        @inbounds M[i,i] = zero(T)
    end
    return M
end

"""
        norm_sum_rows!(mat)
Rescales the matrix by row so that the sum of each row is 1.0
"""
function norm_sum_rows!(mat::AbstractMatrix{T}) where T<:Real
    normf=abs.(inv.(sum(mat,dims=2)))
    broadcast!(*,mat,mat,normf)
    return mat
end

 # two populations means Dale!
function make_wmat(ne::I,ni::I,wmeanf::M ; noautapses=true) where
            {I<: Integer,M <: AbstractMatrix{<:Real}}
    nw = size(wmeanf,1)
    @assert (nw == 2) || (nw==ne+ni)
    if nw == ne+ni #nothing to do
        return wmeanf
    end
    # else, expand the 2D
    Wμ = wmeanf
    d = Exponential(default_values[:muw])
    W=rand(d,ne+ni,ne+ni)
    noautapses && diagtozero!(W)
    # now normalize by block
    pe,pi = populations(ne,ni)
    let wee =  view(W,pe,pe)
        norm_sum_rows!(wee)
        wee .*= Wμ[1,1]
    end
    let wii =  view(W,pi,pi)
        norm_sum_rows!(wii)
        wii .*= Wμ[2,2]
    end
    let wei =  view(W,pe,pi) # inhibitory to excitatory
        norm_sum_rows!(wei)
        wei .*= Wμ[1,2]
    end
    let wie =  view(W,pi,pe) # excitatory to inhibitory
        norm_sum_rows!(wie)
        wie .*= Wμ[2,1]
    end
    return W
end

abstract type RecurrentNetwork{R} end

struct BaseNetwork{R} <: RecurrentNetwork{R}
    weights::Matrix{R}
    iofunction::IOFunction{R}
    membrane_taus::Vector{R}
    external_input::Vector{R}
end
struct AttractorNetwork{R} <: RecurrentNetwork{R}
    weights::Matrix{R}
    iofunction::IOFunction{R}
    attractors_u::Matrix{R}
    membrane_taus::Vector{R}
    external_input::Vector{R}
end

n_neurons(rn::RecurrentNetwork) = size(rn.weights,1)
n_attractors(rn::AttractorNetwork) = size(rn.attractors_u,2)
Base.ndims(rn::RecurrentNetwork)=n_neurons(rn)

#Base.Broadcast.broadcastable(g::RecurrentNetwork)=Ref(g)

Base.copy(ntw::BaseNetwork{R}) where R =
    BaseNetwork{R}( (copy(getfield(ntw,n))
        for n in fieldnames(BaseNetwork) )...)
Base.copy(ntw::AttractorNetwork{R}) where R =
    AttractorNetwork{R}( (copy(getfield(ntw,n))
                for n in fieldnames(AttractorNetwork) )...)

@inline function iofunction!(dest::AbstractVector{R},u::AbstractVector{R},
    nt::RecurrentNetwork{R}) where R<:Real
  for i in eachindex(dest)
      dest[i]=nt.iofunction(u[i])
  end
  return dest
end
@inline function ioprime!(dest::Vector{R},u::AbstractVector{R},
    nt::RecurrentNetwork{R}) where R<:Real
  for i in eachindex(dest)
      dest[i]=ioprime(u[i],nt.iofunction)
  end
  return dest
end
@inline function ioprime(u::AbstractVector{R},nt::RecurrentNetwork{R}) where R<:Real
    return ioprime!(similar(u),u,nt)
end
@inline function ioprimeprime!(dest::AbstractVector{R},
    u::AbstractVector{R},nt::RecurrentNetwork{R}) where R<:Real
  for i in eachindex(dest)
      dest[i]=ioprimeprime(u[i],nt.iofunction)
  end
  return dest
end
@inline function ioprimeprime(u::AbstractVector{R},nt::RecurrentNetwork{R}) where R<:Real
    return ioprimeprime!(similar(u),u,nt)
end

@inline function ioinv!(dest::Vector{R},u::AbstractVector{R},
    nt::RecurrentNetwork{R}) where R<:Real
  for i in eachindex(dest)
      dest[i]=ioinv(u[i],nt.iofunction)
  end
  return dest
end
@inline function ioinv(u::AbstractVector{R},nt::RecurrentNetwork{R}) where R<:Real
    return ioinv!(similar(u),u,nt)
end

function normalize_membrane_taus!(v::Vector{R},ntw::RecurrentNetwork{R}) where R
    for (i,τ) in enumerate(ntw.membrane_taus)
        v[i] /= τ
    end
    return v
end

function BaseNetwork(ne::I,ni::I;
        gfun::Union{Nothing,G}=nothing,
        taus::Union{Nothing,V}=nothing,
        external_current::Union{Nothing,V}=nothing,
        W::Union{Nothing,M}=nothing,
        noautapses = true) where
                {I<:Integer, G<:IOFunction,M<:AbstractMatrix,V<:AbstractVector}

   W = if isnothing(W)
       make_wmat(ne,ni,default_values[:Wmeanf];noautapses=noautapses)
   # if full D , just use it
   elseif size(W,1)==ne+ni
       W
   # if 2D use as mean field
   elseif (size(W,1) == 2 ) && (size(W,2) == 2 )
       make_wmat(ne,ni,W;noautapses=noautapses)
   else
       @warn "something might be wrong with the weight matrix"
       W
   end
   # now the taus
   taus =
       if isnothing(taus)
           τe,τi =default_values[:tau_ei]
           vcat(fill(τe,ne),fill(τi,ni))
       elseif length(taus)==2
           τe,τi =taus
           vcat(fill(τe,ne),fill(τi,ni))
       else
           taus
    end
    # currents
    h  = something(external_current,
            fill(default_values[:external_current],ne+ni) )
    # gain
    gfun = something(gfun, IOQuad(default_values[:gainalpha]))
   return BaseNetwork(W,gfun,taus,h)
end
function AttractorNetwork(bn::BaseNetwork{R},attr::Matrix{R}) where R
    return AttractorNetwork{R}(bn.weights,bn.iofunction,attr,
        bn.membrane_taus,bn.external_input)
end

function velocity!(dest::AbstractVector{R},u::AbstractVector{R},
        gu::AbstractVector{R},rn::RecurrentNetwork{R}) where R<:Real
    copy!(dest, rn.external_input)
    dest .-= u  #  v_out <-  - u +  h
    LinearAlgebra.BLAS.gemv!('N',1.0,rn.weights,gu,1.0,dest) # W*g(u) - u + h
    normalize_membrane_taus!(dest,rn) # ( W*g(u) - u + h) / taus
    return dest
end
velocity(u,rn) = velocity!(similar(u),u,rn.iofunction.(u), rn)

# Jacobian stuff ! Compute and derivatives!
# let's define (and test) the derivatives of the jacobian here!

struct JAlloc{R}
    dgu::Vector{R}
    ddgu::Vector{R}
    inv_taus::Vector{R}
end
function JAlloc(ntw::RN) where {R<:Real,RN<:RecurrentNetwork{R}}
    dgu = similar(ntw.membrane_taus)
    ddgu = similar(ntw.membrane_taus)
    JAlloc{R}(dgu,ddgu,inv.(ntw.membrane_taus))
end

function jacobian!(Jdest::AbstractMatrix{R},
        Wgrad::Union{Nothing,Matrix{R}},ugrad::Union{Nothing,Matrix{R}},
        u::AbstractVector{R},ntw::RecurrentNetwork{R},
        Jallocs::Union{Nothing,JAlloc{R}}=nothing) where R
  Jallocs=something(Jallocs,JAlloc(ntw))
  ioprime!(Jallocs.dgu,u,ntw)
  broadcast!(*,Jdest, ntw.weights, transpose(Jallocs.dgu)) # multiply columnwise
  for i in 1:n_neurons(ntw)
    Jdest[i,i] -= 1.0 #subtract diagonal
  end
  #normalize by taus, rowwise
  broadcast!(/,Jdest,Jdest,ntw.membrane_taus)
  # GRADIENTS ! with respect to  W first
  if !isnothing(Wgrad)
    broadcast!(*,Wgrad,Jallocs.inv_taus,transpose(Jallocs.dgu))
  end
  # now with respect to u
  if !isnothing(ugrad)
    ioprimeprime!(Jallocs.ddgu,u,ntw)
    #broadcast!(/,Jallocs.weights,transpose(dgu),ntw.membrane_taus)
    broadcast!(*,ugrad,Jallocs.inv_taus,
      ntw.weights, transpose(Jallocs.ddgu))
  end
  return Jdest
end

"""
        jacobian(u,rn::RecurrentNetwork)
Jacobian matrix of the
"""
@inline function jacobian(u::AbstractVector{R},rn::RecurrentNetwork{R}) where R
    return jacobian!(similar(rn.weights),nothing,nothing,u,rn)
end

function spectral_abscissa(u,rn::RecurrentNetwork)
    J=jacobian(u,rn)
    return maximum(real.(eigvals(J)))
end

# Dynamics
function run_network(ntw::RecurrentNetwork{R},r_start::Vector{R},
    t_end::Real; verbose::Bool=false,stepsize=0.05) where R<:Real
  u0=ioinv(r_start,ntw)
  iofun_alloc=similar(r_start)
  f(du,u,p,t) = velocity!(du,u,iofunction!(iofun_alloc,u,ntw),ntw)
  prob = ODEProblem(f,u0,(0.,t_end))
  solv = solve(prob,Tsit5();verbose=verbose,saveat=stepsize)
  ret_u = hcat(solv.u...)
  ret_r = ntw.iofunction.(ret_u)
  return solv.t,ret_u,ret_r
end

function run_network_noise(ntw::RecurrentNetwork{R},
    r_start::Vector{R},noiselevel::Real,t_end::Real;
    verbose::Bool=false,stepsize=0.05) where R<:Real
  u0=ioinv(r_start,ntw)
  iofun_alloc=similar(r_start)
  f(du,u,p,t) = velocity!(du,u,iofunction!(iofun_alloc,u,ntw),ntw)
  σ_f(du,u,p,t) = fill!(du,noiselevel)
  prob = SDEProblem(f,σ_f,u0,(0.,t_end))
  solv =  solve(prob,EM();verbose=verbose,dt=stepsize)
  ret_u = hcat(solv.u...)
  ret_r = ntw.iofunction.(ret_u)
  return solv.t,ret_u,ret_r
end

"""
        run_network_to_convergence(rn::RecurrentNetwork,r_start::Vector{<:Real} ;
                t_end=80. , veltol=1E-1)
Runs the network as described in [`run_network_nonoise`](@ref), but stops as soon as
`norm(v) / n < veltol` where `v` is the velocity at time `t`.
If this condition is not satisfied (no convergence to attractor), it runs until `t_end` and prints a warning.
# Arguments
- `rn::RecurrentNetwork`
- `r_start::Vector` : initial conditions
- `t_end::Real`: the maximum time considered
- `veltol::Real` : the norm (divided by num. dimensions) for velocity at convergence
# Outputs
- `u_end`::Vector : the final state at convergence
- `r_end`::Vector : the final state at convergence as rate
"""
function run_network_to_convergence(ntw::RecurrentNetwork{R},
     r_start::Vector{R};t_end::Real=50. , veltol::Real=1E-1) where R
    function  condition(u,t,integrator)
        v = get_du(integrator)
        return norm(v) / length(v) < veltol
    end
    function affect!(integrator)
        savevalues!(integrator)
        return terminate!(integrator)
    end
    u0=ioinv(r_start,ntw)
    cb=DiscreteCallback(condition,affect!)
    iofun_alloc=similar(r_start)
    f(du,u,p,t) = velocity!(du,u,iofunction!(iofun_alloc,u,ntw),ntw)
    prob = ODEProblem(f,u0,(0.,t_end))
    out = solve(prob,Tsit5();verbose=false,callback=cb)
    u_out = out.u[end]
    t_out = out.t[end]
    if isapprox(t_out,t_end; atol=0.05)
        vel = velocity(u_out,ntw)
        @warn "no convergence after max time $t_end"
        @info "the norm (divided by n) of the velocity is $(norm(vel)/length(vel)) "
    end
    ret_r = ntw.iofunction.(u_out)
    return u_out,ret_r
end

# here I make the AttractorNetwork, the one with attractors

"""
        lognorm_reparametrize(m,std) -> distr::LogNormal
# parameters
  + `m`   sample mean
  + `std` sample std
"""
function lognorm_reparametrize(m,std)
    vm2= (std/m)^2
    μ = log(m / sqrt(1. + vm2))
    σ = sqrt(log( 1. + vm2))
    return LogNormal(μ,σ)
end

function make_attractors(ntot,natt,gainf::IOFunction{R};
        mu_r=5.0,std_r=2.0) where R
    distr = lognorm_reparametrize(mu_r,std_r)
    attr_r = rand(distr,(ntot,natt))
    attr_u = ioinv.(attr_r,gainf)
    return (attr_r,attr_u)
end
make_attractors(ntot,natt,ntw::RecurrentNetwork;mu_r=5.0,std_r=2.0) =
    make_attractors(ntot,natt,ntw.iofunction;mu_r=mu_r,std_r=std_r)

function AttractorNetwork(natt::I,ne::I,ni::I;
        gfun::Union{Nothing,G}=nothing,
        taus::Union{Nothing,V}=nothing,
        external_current::Union{Nothing,V}=nothing,
        W::Union{Nothing,M}=nothing,
        noautapses::Bool = true ,
        mu_attr::Real=5.0,
        std_attr::Real=2.0) where
                {I<:Integer, G<:IOFunction,M<:AbstractMatrix,V<:AbstractVector}
  base=BaseNetwork(ne,ni;gfun=gfun,taus=taus,external_current=external_current,
    W=W,noautapses=noautapses)
  _,attr_u = make_attractors(ne+ni,natt,base;mu_r=mu_attr,std_r=std_attr)
  return AttractorNetwork(base,attr_u)
end

end # module


text:
#=
	Author @maurxeugenio
	Algoritmo de ordenação buble sort, feito em Julia [juja]
	[os comentários em muitas linhas é bem mais lindo que em python.]
=#

lista_desordenada = [10, 3, 4, 5, 6, 1, 4, 5,6, 30, 1, 500, 200]

function odernar(lista)
    for i in 1:length(lista)
		for j in 1:(length(lista) - 1)
			if lista[j] > lista[j+1]
				aux = lista[j]
				lista[j] = lista[j+1]
				lista[j+1] = aux
			end
		end
	end
	return lista
end

print(odernar(lista_desordenada))

max_stars_repo_name: maurxeugenio/juja

max_issues_repo_name: maurxeugenio/juja

max_forks_repo_name: maurxeugenio/juja


text:
lane() = 1

@spmd lane() = vect(ntuple(identity, length(__mask__))...)

lanewidth() = 1

@spmd lanewidth() = length(__mask__)

lanesum(x) = lanewidth() * x

@spmd lanesum(x::AbstractVec) = sum(x)
@spmd lanesum(x) = length(__mask__) * x # TODO remove this

function shuffle(x, i)
  i == 1 || error("Tried to shuffle from lane $i of 1")
  return x
end

@spmd shuffle(x, i) = x
@spmd shuffle(x::AbstractVec, i) = x[i]


text:
using CuArrays, CUDAnative, CUDAdrv
using Test

# Kernel
function add!(a,b,c, index)
    c[index] = a[index] + b[index]
    return nothing
end

function main()

    # Initialise a, b and c
    a = CuArrays.CuArray(fill(0, 10))
    b = CuArrays.CuArray(fill(0, 10))
    c = CuArrays.CuArray(fill(0, 10))

    # Put values in a and b
    for i in 1:10
        a[i] = -i
        b[i] = i * i
    end

    # Create two streams
    s1 = CuStream()
    s2 = CuStream()

    # Call add! asynchronously in two streams
    for i in 1:2:min(length(a), length(b), length(c))
        @cuda threads = 1 stream = s1 add!(a,b,c, i)
        @cuda threads = 1 stream = s2 add!(a,b,c, i+1)
    end

    # Copy arrays back to host (CPU)
    a=Array(a)
    b=Array(b)
    c=Array(c)

    # Check the addition worked
    for i in 1:length(a)
        @test a[i] + b[i] ≈ c[i]
    end

end

main()


text:
######################
### Run simulation ###
######################

"""
Run one simulation iteration

**Arguments**
* `sim_data` : `SimData` object
* `buffer` : list of `OpenStreetMapXSim.Road` objects containing informations about routes selected during the simulation run.
* `nodes_stats` : dictionary of `OpenStreetMapXSim.NodeStat` objects containing informations about each intersection in simulation.
* `destination_selector` : string determining a way how the destination (workplace) will be selected (based on journey matrix, business data or on both)
* `weight_var` : weighting variable name (or nothing)
* `google` : boolean variable; if true simulation will generates routes based on Google Distances API
"""
function run_once!(sim_data::OpenStreetMapXSim.SimData,
                            buffer::Array{OpenStreetMapXSim.Road,1},
                            nodes_stats::Dict{Int,OpenStreetMapXSim.NodeStat},
                            destination_selector::String,
							agentid::Int,demographic_profile::Function,
							additional_activity::Function;
							weight_var:: Union{Symbol,Nothing} = nothing,
							google::Bool = false
							)
    loc = OpenStreetMapXSim.start_location(sim_data.demographic_data, weight_var = weight_var)
    agent = demographic_profile(loc, sim_data.demographic_data[loc])
    agent[:id]=agentid
    if destination_selector == "flows"
        OpenStreetMapXSim.destination_location!(agent,sim_data.DAs_flow_dictionary,sim_data.DAs_flow_matrix)
    elseif destination_selector == "business"
        OpenStreetMapXSim.destination_location!(agent,sim_data.business_data)
	else
		OpenStreetMapXSim.destination_location!(agent,sim_data)
    end
    #before work
    activity = additional_activity(agent,true,sim_data)
    if isa(activity,Nothing)
        routebefore = OpenStreetMapXSim.select_route(agent.DA_home[1], agent.DA_work[1],sim_data, buffer, google = google)
    else
        routebefore = OpenStreetMapXSim.select_route(agent.DA_home[1], agent.DA_work[1], activity,sim_data, buffer,google = google)
    end
    OpenStreetMapXSim.stats_aggregator!(nodes_stats, agent, routebefore)
    local routeafter
    #after work
    activity = additional_activity(agent,false,sim_data)
    if isa(activity,Nothing)
        routeafter = OpenStreetMapXSim.select_route(agent.DA_work[1], agent.DA_home[1], sim_data, buffer,google = google)
    else
        routeafter = OpenStreetMapXSim.select_route(agent.DA_work[1], agent.DA_home[1], activity, sim_data, buffer,google = google)
    end
    OpenStreetMapXSim.stats_aggregator!(nodes_stats, agent, routeafter)
    return (routebefore,routeafter)
end

"""
Run simulation

Run simulation for data stored in `SimData` object

**Arguments**
* `sim_data` : `SimData` object
* `destination_selector` : string determining a way how the destination (workplace) will be selected (based on journey matrix, business data or on both)
* `N` : number of iterations
* `demographic_profile` : function generating a demographic profile on the base of DA_home and DA_demostat
* `additional_activity` : function generating an additional_activity on the base agents profile and direction
* `weight_var` : weighting variable name (or nothing)
* `google` : boolean variable; if true simulation will generates routes based on Google Distances API
"""
function run_simulation(sim_data::OpenStreetMapXSim.SimData,
            destination_selector::String,
			job::Int,
            N::Int,
			demographic_profile::Function,
			additional_activity::Function,
			weight_var::Union{Symbol,Nothing};
			google::Bool = false)
	if !in(destination_selector,["flows","business","both"])
		error("destination_selector not declared properly! It can only takes flows, business or both values!")
	end
    nodes_stats = OpenStreetMapXSim.node_statistics(sim_data)
    buffer = Array{OpenStreetMapXSim.Road,1}()
    routes = Dict()
	startt = Dates.now()
	@info "Worker $(myid()) Starting simulation for seed $job at $startt"
    Random.seed!(job);
    for i = 1:N
        agentid = (job*1000000) + i
        routes[agentid] = OpenStreetMapXSim.run_once!(sim_data,buffer,nodes_stats,destination_selector,agentid,demographic_profile,additional_activity, weight_var = weight_var, google = google)
		i == 1 && @info "Worker: $(Distributed.myid()) First out of $N simulation completed"
    end
    #$(Distributed.myid())
	@info "Worker: $(Distributed.myid()) All $N sims completed with time per sim = $((Dates.now()-startt).value/N)ms"
    return nodes_stats,buffer,routes
end

function run_dist_sim(resultspath,version::String,master_ID::Int,N::Int,max_jobs_worker::Int,
					  sim_data::OpenStreetMapXSim.SimData, mode::String,
					  demographic_profile,additional_activity,weight_var;
					  s3action::Union{Function,Nothing}=nothing)
	for ii in 1:max_jobs_worker
	    d = master_ID*max_jobs_worker+ii;
	    nodes, buffer,routes = run_simulation(sim_data, mode, d, N,
		        demographic_profile,additional_activity,weight_var);

		# add nodeids and distributed.myid to nodes statistics
		nodeids = collect(keys(nodes));
		for i in nodeids
			if nodes[i].agents_data != nothing
				insert!(nodes[i].agents_data, 1, nodes[i].latitude, :latitude)
				insert!(nodes[i].agents_data, 1, nodes[i].longitude, :longitude)
				insert!(nodes[i].agents_data, 1, i, :NODE_ID)
				insert!(nodes[i].agents_data, 1, Distributed.myid(), :DISTRIBUTED_ID)
			else
				delete!(nodes, i)
			end
		end

	    # merge results
		results = collect(values(nodes))[1].agents_data
		for df in collect(values(nodes))[2:end]
	       append!(results, df.agents_data)
	    end
	    filenamebase="res_V$(version)_M$(mode)_W$(@sprintf("%04d", Distributed.myid()))_S$(@sprintf("%05d",d))"
	    filenodes  = "$(filenamebase)_nodes.csv"
	    targetfile=joinpath(resultspath,filenodes)
		Nanocsv.write_csv(targetfile,results, delim = ';')
		s3action != nothing && s3action(resultspath, filenodes)
	    fileroutes = "$(filenamebase)_routes.csv"
	    f = open(joinpath(resultspath, fileroutes),"w")
	    for agentid in keys(routes)
	        mode = "towork"
	        for rset in routes[agentid]
	            print(f,"$(agentid);$(mode)")
	            for nn in rset
	               print(f,";$nn")
	            end
	            println(f,"")
	            mode = "tohome"
	        end
	    end
	    close(f)
	    s3action != nothing && s3action(resultspath, fileroutes)
		@info "Results exported for worker $(myid()) seed $d to: $targetfile"
	    #1
	end
end


text:
@testset "rotations" begin
    qq = Quat(1.0, 0.01, 0.02, 0.03)
    rv = RodriguesVec(qq)
    rv_sv = SVector(rv.sx, rv.sy, rv.sz)
    rv_cheap = cheapRV(qq)
    rel_err = (rv_sv - rv_cheap) ./ rv_sv
    @test all(rel_err .< 1.0e-3)

    t = rand(SVector{4,Float64})
    q = Quat(t..., false)
    @test components(q) == t

    t = rand(SVector{3,Float64})
    s = SPQuat(t...)
    @test components(s) == t

    q = rand(Quat{Float64})
    s = SPQuat(q)
    @test cheapRV(q) ≈ cheapRV(s)

    q_0 = Quat(0.1, 0.2, 0.3, 0.4)
    for k = 1:100
        q_1 = rand(Quat{Float64})
        q_2 = Quat(-components(q_1)..., false)
        q_err_1 = quatErr(q_1, q_0)
        q_err_2 = quatErr(q_2, q_0)
        @test q_err_1 ≈ q_err_1
        q_err_check = components(q_1 * inv(q_0))
        if q_err_check[1] <= 0.0
            @test q_err_check[2:4] ≈ -q_err_1
        else
            @test q_err_check[2:4] ≈ q_err_1
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "raptor"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/tbronzwaer/raptor.git", "6cf2a5903af340c3cbf07fc9e13d43364c3187b4")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/raptor/
make CPU="1" CFLAGS="-fopenmp -std=c99 -I${includedir} -Wno-unused-result" LDFLAGS="-lm -lgsl -lgslcblas"
cp RAPTOR ${bindir}/.
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("RAPTOR", :RAPTOR),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@info "Loading packages"
import FoldsChainRules
import FoldsChainRulesBenchmarks

@info "Start warmup"
run(FoldsChainRulesBenchmarks.setup_smoke(); verbose = true)

@info "Start benchmarks"
SUITE = FoldsChainRulesBenchmarks.setup()
RESULTS = run(SUITE; verbose = true)

if abspath(PROGRAM_FILE) == @__FILE__
    display(RESULTS)
    println()
end

RESULTS


text:
using LinearAlgebra: diag, eigen, eigvals, eigvecs, Symmetric, Diagonal, tr, norm, cholesky, LowerTriangular, UpperTriangular

@doc doc"""
    CholeskySpace{N} <: Manifold

the manifold of lower triangular matrices with positive diagonal and
a metric based on the cholesky decomposition. The formulae for this manifold
are for example summarized in Table 1 of
> Lin, Zenhua: "Riemannian Geometry of Symmetric Positive Definite Matrices via
> Cholesky Decomposition", arXiv: [1908.09326](https://arxiv.org/abs/1908.09326).
"""
struct CholeskySpace{N} <: Manifold end
CholeskySpace(n::Int) = CholeskySpace{n}()

# two small helper for strictly lower and upper triangulars
strictlyLowerTriangular(x) = LowerTriangular(x) - Diagonal(diag(x))
strictlyUpperTriangular(x) = UpperTriangular(x) - Diagonal(diag(x))

@generated representation_size(::CholeskySpace{N}) where N = (N,N)

@generated manifold_dimension(::CholeskySpace{N}) where N = div(N*(N+1), 2)

@doc doc"""
    inner(M,x,v,w)

computes the inner product on the [`CholeskySpace`](@ref) `M` at the
lower triangular matric with positive diagonal `x` and the two tangent vectors
`v`,`w`, i.e they are both lower triangular matrices with arbitrary diagonal.
The formula reads

````math
    g_{x}(v,w) = \sum_{i>j} v_{ij}w_{ij} + \sum_{j=1}^m v_{ii}w_{ii}x_{ii}^{-2}
````
"""
inner(::CholeskySpace{N},x,v,w) where N = sum(strictlyLowerTriangular(v).*strictlyLowerTriangular(w)) + sum(diag(v).*diag(w)./( diag(x).^2 ))

@doc doc"""
    distance(M,x,y)

computes the Riemannian distance on the [`CholeskySpace`](@ref) `M` between two
matrices `x`, `y` that are lower triangular with positive diagonal. The formula
reads

````math
d_{\mathcal M}(x,y) = \sqrt{
\sum_{i>j} (x_{ij}-y_{ij})^2 +
\sum_{j=1}^m (\log x_{jj} - \log y_{jj})^2
}
````
"""
distance(::CholeskySpace{N},x,y) where N = sqrt(
  sum( (strictlyLowerTriangular(x) - strictlyLowerTriangular(y)).^2 ) + sum( (log.(diag(x)) - log.(diag(y))).^2 )
)

@doc doc"""
    exp!(M,y,x,v)

compute the exponential map on the [`CholeskySpace`](@ref) `M` eminating from
the lower triangular matrix with positive diagonal `x` towards the lower triangular
matrix `v` and return the result in `y`. The formula reads

````math
\exp_x v = \lfloor x \rfloor + \lfloor v \rfloor
+\operatorname{diag}(x)\operatorname{diag}(x)\exp\bigl( \operatorname{diag}(v)\operatorname{diag}(x)^{-1}\bigr)
````
where $\lfloor x\rfloor$ denotes the strictly lower triangular matrix of $x$ and
$\operatorname{diag}(x)$ the diagonal matrix of $x$
"""
function exp!(::CholeskySpace{N},y,x,v) where N
    y .= strictlyLowerTriangular(x) + strictlyLowerTriangular(v) + Diagonal(x)*Diagonal(exp.(diag(v)./diag(x)))
    return y
end
@doc doc"""
    log!(M,v,x,y)

compute the exponential map on the [`CholeskySpace`](@ref) `M` eminating from
the lower triangular matrix with positive diagonal `x` towards the lower triangular
matrx `v` and return the result in `y`. The formula reads

````math
\log_x v = \lfloor x \rfloor - \lfloor y \rfloor
+\operatorname{diag}(x)\log\bigl(\operatorname{diag}(y)\operatorname{diag}(x)^{-1}\bigr)
````
where $\lfloor x\rfloor$ denotes the strictly lower triangular matrix of $x$ and
$\operatorname{diag}(x)$ the diagonal matrix of $x$
"""
function log!(::CholeskySpace{N},v,x,y) where N
    v .= strictlyLowerTriangular(y) - strictlyLowerTriangular(x) + Diagonal(diag(x))*Diagonal(log.(diag(y)./diag(x)))
    return v
end

@doc doc"""
    zero_tangent_vector!(M,v,x)

returns the zero tangent vector on the [`CholeskySpace`](@ref) `M` at `x` in
the variable `v`.
"""
function zero_tangent_vector!(M::CholeskySpace{N},v,x) where N
    fill!(v,0)
    return v
end

@doc doc"""
    vector_transport!(M,vto,x,v,y)

parallely transport the tangent vector `v` at `x` along the geodesic to `y`
on respect to the [`CholeskySpace`](@ref) manifold `M`. The formula reads

````math
    \mathcal P_{x\to y}(v) = \lfloor v \rfloor  + \operatorname{diag}(y)\operatorname{diag}(x)^{-1}\operatorname{diag}(v),
````
where $\lfloor\cdot\rfloor$ denotes the strictly lower triangular matrix,
and $\operatorname{diag}$ extracts the diagonal matrix.
"""
function vector_transport_to!(::CholeskySpace{N}, vto, x, v, y, ::ParallelTransport) where N
    vto .= strictlyLowerTriangular(x) + Diagonal(diag(y))*Diagonal(1 ./ diag(x))*Diagonal(v)
    return vto
end

@doc doc"""
    check_manifold_point(M,x;kwargs...)

check whether the matrix `x` lies on the [`CholeskySpace`](@ref) `M`, i.e.
it's size fits the manifold, it is a lower triangular matrix and has positive
entries on the diagonal.
The tolerance for the tests can be set using the `kwargs...`.
"""

function check_manifold_point(M::CholeskySpace{N}, x; kwargs...) where N
    if size(x) != representation_size(M)
        return DomainError(size(x),"The point $(x) does not lie on $(M), since its size is not $(representation_size(M)).")
    end
    if !isapprox( norm(strictlyUpperTriangular(x)), 0.; kwargs...)
        return DomainError(norm(UpperTriangular(x) - Diagonal(x)), "The point $(x) does not lie on $(M), since it strictly upper triangular nonzero entries")
    end
    if any( diag(x) .<= 0)
        return DomainError(min(diag(x)...), "The point $(x) does not lie on $(M), since it hast nonpositive entries on the diagonal")
    end
    return nothing
end
"""
    check_tangent_vector(M,x,v; kwargs... )

checks whether `v` is a tangent vector to `x` on the [`CholeskySpace`](@ref) `M`, i.e.
atfer [`check_manifold_point`](@ref)`(M,x)`, `v` has to be of same dimension as `x`
and a symmetric matrix.
The tolerance for the tests can be set using the `kwargs...`.
"""
function check_tangent_vector(M::CholeskySpace{N}, x,v; kwargs...) where N
    mpe = check_manifold_point(M, x)
    if mpe !== nothing
        return mpe
    end
    if size(v) != representation_size(M)
        return DomainError(size(v),
            "The vector $(v) is not a tangent to a point on $(M) since its size does not match $(representation_size(M)).")
    end
    if !isapprox(norm(v-transpose(v)), 0.; kwargs...)
        return DomainError(size(v),
            "The vector $(v) is not a tangent to a point on $(M) (represented as an element of the Lie algebra) since its not symmetric.")
    end
    return nothing
end


text:
# primitives, strings and arrays

writeobj(io::IO, n::Integer) = write(io, hton(n))
readobj(io::IO, ::Type{T}) where T<:Integer = ntoh(read(io, T))

function writeobj(io::IO, s::String)
    len = Int16(length(s))
    writeobj(io, len > 0 ? len : -1)
    write(io, s)
end
function readobj(io::IO, ::Type{String})
    len = readobj(io, Int16)
    return len > 0 ? String(read(io, len)) : ""
end

function writeobj(io::IO, arr::Vector{T}) where T
    len = Int32(length(arr))
    writeobj(io, len > 0 ? len : -1)
    for x in arr
        writeobj(io, x)
    end
end
function readobj(io::IO, ::Type{Vector{T}}) where T
    len = readobj(io, Int32)
    if len <= 0
        return T[]
    end
    arr = Array{T}(undef, len)
    for i=1:len
        arr[i] = readobj(io, T)
    end
    return arr
end

# tuples
function writeobj(io::IO, tp::Tuple)
    for x in tp
        writeobj(io, x)
    end
end
# there should be a better way to do it, but overloading for
# 1-4 element tuples works too
function readobj(io::IO, ::Type{Tuple{T}}) where T
    return readobj(io, T)
end
function readobj(io::IO, ::Type{Tuple{T1,T2}}) where {T1,T2}
    v1 = readobj(io, T1)
    v2 = readobj(io, T2)
    return (v1, v2)
end
function readobj(io::IO, ::Type{Tuple{T1,T2,T3}}) where {T1,T2,T3}
    v1 = readobj(io, T1)
    v2 = readobj(io, T2)
    v3 = readobj(io, T3)
    return (v1, v2, v3)
end
function readobj(io::IO, ::Type{Tuple{T1,T2,T3,T4}}) where {T1,T2,T3,T4}
    v1 = readobj(io, T1)
    v2 = readobj(io, T2)
    v3 = readobj(io, T3)
    v4 = readobj(io, T4)
    return (v1, v2, v3, v4)
end

# composite types
function writeobj(io::IO, o)
    for f in fieldnames(typeof(o))
        writeobj(io, getfield(o, f))
    end
end
function readobj(io::IO, ::Type{T}) where T
    vals = Array{Any}(undef, length(T.types))
    for (i, t) in enumerate(T.types)
        vals[i] = readobj(io, t)
    end
    return T(vals...)
end

# Weird part: OffsetMessage, MessageSet and FetchResponsePartitionData
# Unlike other data types, OffsetMessage and MessageSet can't
# be encoded/decoded as a sum of its parts, but require special rules:
#
# 1. OffsetMessage contains field message_size with actual size (in bytes)
#    of encoded message. In most cases it is the same as if message were decoded
#    normally. But in some specific cases (e.g. when value is null), message
#    may contain additional bytes (for which I couldn't understand the meaning).
#    So we need to read message_size bytes from a stream to a buffer and
#    then parse actual message from that buffer.

function readobj(io::IO, ::Type{OffsetMessage})
    offset = readobj(io, Int64)
    message_size = readobj(io, Int32)
    data = read(io, message_size)
    buf = IOBuffer()
    write(buf, data)
    seek(buf, 0)
    msg = readobj(buf, Message)
    return OffsetMessage(offset, message_size, msg)
end

# 2. MessageSet is essentially an array of OffsetMessages, but unlike other
#    arrays it doesn't contain Int32 prefix of length. Instead, its
#    parent struct (during reading) - FetchResponsePartitionData - contains a field
#    message_set_size. Just as with OffsetMessage, we need first to read
#    that much bytes into a buffer and then deserialize all full instances of
#    OffsetMessage. Note that data is copied from brokers as byte arrays
#    so the resulting buffer may contain some partial messages which we dismiss

function writeobj(io::IO, message_set::MessageSet)
    # skip array header
    for msg_elem in message_set.elements
        writeobj(io, msg_elem)
    end
end
function readobj(io::IO, ::Type{MessageSet}, message_set_size::Int32)
    message_set_bytes = read(io, message_set_size)
    buf = IOBuffer()
    write(buf, message_set_bytes)
    seek(buf, 0)
    elements = OffsetMessage[]
    while !eof(buf)
        try
            elem = readobj(buf, OffsetMessage)
            push!(elements, elem)
        catch EOFError
            # incomplete message
        end
    end
    return MessageSet(elements)
end
function readobj(io::IO, ::Type{FetchResponsePartitionData})
    partition = readobj(io, Int32)
    error_code = readobj(io, Int16)
    highwater_mark_offset = readobj(io, Int64)
    message_set_size = readobj(io, Int32)
    # we don't know how many messeges there are, so reading that much bytes
    # and trying to parse as many messages as we can
    message_set = readobj(io, MessageSet, message_set_size)
    return FetchResponsePartitionData(partition, error_code, highwater_mark_offset,
                                message_set_size, message_set)
end

# requests and responses

# TODO: server waits for some more bytes
#       one possible reason is that the length of objects is calculated incorrectly
function obj_size(objs...)
    buf = IOBuffer()
    for obj in objs
        writeobj(buf, obj)
    end
    return Int32(buf.size)
end

function send_request(io::IO, header::RequestHeader, req)
    len = obj_size(header, req)
    writeobj(io, len)
    writeobj(io, header)
    writeobj(io, req)
end
function send_request(io::IO, header::RequestHeader)
    len = obj_size(header)
    writeobj(io, len)
    writeobj(io, header)
end
function recv_header(io::IO)
    readobj(io, Int32) # length, do we really need it?
    return readobj(io, ResponseHeader)
end
function recv_response_with_header(io::IO, ::Type{T}) where T
    header = recv_header(io)
    resp = readobj(io, T)
    return header, resp
end
recv_response(io::IO, ::Type{T}) where T = recv_response_with_header(io, T)[2]


text:
using MatchCore
using Julog
#
# # example assuming * operation is always binary

struct JulogAnalysis <: AbstractAnalysis
    egraph::EGraph
    facts::Vector{Clause}
    data::Dict{Int64, Vector{Clause}}
end

# Mandatory for AbstractAnalysis
JulogAnalysis(g::EGraph, facts::Vector{Clause}) =
    JulogAnalysis(g, facts, Dict{Int64, Vector{Clause}}())

# This should be auto-generated by a macro
function EGraphs.make(an::JulogAnalysis, n)
    !(n isa Expr) && return []

    m = @matcher begin
        l * r |> begin
            get(an, l, [])
        end

        _ => nothing
    end

    return m(n)
end

EGraphs.join(analysis::NumberFold, from, to) = from ∪ to

function EGraphs.modify!(an::NumberFold, id::Int64)
    g = an.egraph
    # !haskey(an, id) && return nothing
    if an[id] isa Number
        newclass = EGraphs.add!(g, an[id])
        merge!(g, newclass.id, id)
    end
end

Base.setindex!(an::JulogAnalysis, value, id::Int64) = setindex!(an.data, value, id)
Base.getindex(an::JulogAnalysis, id::Int64) = an.data[id]
Base.haskey(an::JulogAnalysis, id::Int64) = haskey(an.data, id)
Base.delete!(an::JulogAnalysis, id::Int64) = delete!(an.data, id)

g = EGraph(:(x * (cos(z)/cos(z))))
facts = @julog [
    iszero(z) <<= true,
    iszero(X) <<= fail
]
extractor = addanalysis!(g, ExtractionAnalysis, astsize)
# ja = addanalysis!(g, JulogAnalysis, facts)

macro res(goals, facts)
    quote first(resolve(@julog($goals), $facts)) end
end

example = @theory begin
    x/x |> (@res(not(iszero(x)), facts) ? :($x) : :($x/$x))
end

display(g.M); println()
saturate!(g, example; mod=@__MODULE__)
println(extract!(g, extractor))


text:
# Non-transducer reducing function combinators

# Some common definitions for TeeRF and ProductRF
abstract type AbstractMultiCastingRF <: _Function end
# TOOD: Use `AbstractMultiCastingRF <: AbstractReduction` so that
# extra `BottomRF` wrapper can be removed.

start(rf::AbstractMultiCastingRF, init) = map(f -> start(f, init), rf.fs)

# Trying to detect bugs by narrowing method signatures when defining
# `next` for concrete reducing functions. Otherwise, it falls back to
# this method; e.g., foldl(TeeRF(min, max), [5, 2, 6, 8, 3])
next(rf::AbstractMultiCastingRF, accs, x) = error(
    "$rf received an incompatible accumulator of type $(typeof(accs))\n",
    "It requires a $(length(rf.fs))-tuple.\n",
    "Accumulator (first argument): $(accs)\n",
    "Input (second argument): $x",
)

@inline function complete(rf::AbstractMultiCastingRF, accs)
    results = map((f, a) -> complete(f, a), rf.fs, accs)
    if any(r -> r isa DefaultInitOf, results)  # is `any` too strict?
        return DefaultInitOf{typeof(rf)}()
    end
    return results
end

@inline combine(rf::AbstractMultiCastingRF, lefts, rights) =
    map((f, l, r) -> combine(f, l, r), rf.fs, lefts, rights)

# # Something like these are needed for `AbstractMultiCastingRF <: AbstractReduction`:
# Transducer(::AbstractMultiCastingRF) = IdentityTransducer()
# has(rf::AbstractMultiCastingRF, ::Type{<:Transducer}) = false  # should be multicast?
# _realbottomrf(rf::AbstractMultiCastingRF) = rf
# _show_impl(io, m, rf::AbstractMultiCastingRF) = invoke(Base.show, Tuple{IO,Any}, io, rf)

# Adjoining init has to happen at "leafs":
_asmonoid(rf::AbstractMultiCastingRF) =
    ConstructionBase.constructorof(typeof(rf))(map(_asmonoid, rf.fs)::Tuple)

# See: [[./core.jl::TODOs for `Completing` and `skipcomplete`]]
Completing(rf::AbstractMultiCastingRF) =
    ConstructionBase.constructorof(typeof(rf))(map(Completing, rf.fs)::Tuple)

function Base.show(io::IO, rf::AbstractMultiCastingRF)
    @nospecialize
    T = ConstructionBase.constructorof(typeof(rf))
    if rf === T(rf.fs)
        print(io, T, rf.fs)
    else
        invoke(show, Tuple{IO,Any}, io, rf)
    end
end

"""
    TeeRF(reducing_functions::Tuple)
    TeeRF(reducing_functions...)

Combine multiple reducing functions into a new reducing function that
"multicast" the input to multiple reducing functions.

Roughly speaking, `TeeRF(op₁, op₂, ..., opₙ)` is equivalent to

```julia
((a₁, a₂, ..., aₙ), x) -> (op₁(a₁, x), op₂(a₂, x), ..., opₙ(aₙ, x))
```

For [`combine`](@ref), it behaves like [`ProductRF`](@ref).

!!! compat "Transducers.jl 0.4.32"

    New in version 0.4.32.

# Examples
```jldoctest
julia> using Transducers

julia> extrema′(xs, xf = Map(identity)) = foldl(TeeRF(min, max), xf, xs);

julia> extrema′([5, 2, 6, 8, 3])
(2, 8)
```

Note that the input is considered empty unless _all_ reducing
functions call their bottom reducing functions.  Specify `init` to
obain results even when the input collection is empty or all filtered
out.

```jldoctest; setup = :(using Transducers), filter = r"EmptyResultError: .*"
julia> filtering_max = reducingfunction(Filter(isodd), max);

julia> foldl(TeeRF(min, filtering_max), Map(identity), [5, 2, 6, 8, 3])
(2, 5)

julia> foldl(TeeRF(min, filtering_max), Map(identity), 2:2:8)
ERROR: EmptyResultError: ...

julia> foldl(TeeRF(min, filtering_max), Map(identity), 2:2:8; init = Init)
(2, Init(max))
```
"""
struct TeeRF{T<:Tuple} <: AbstractMultiCastingRF
    fs::T
    TeeRF{T}(fs) where {T} = new{T}(fs)
end

TeeRF(fs::Tuple) = TeeRF{typeof(fs)}(fs)
TeeRF(f, fs...) = TeeRF((f, fs...))

@inline next(rf::TeeRF{<:NTuple{N,Any}}, accs::NTuple{N,Any}, x) where {N} =
    map((f, a) -> next(f, a, x), rf.fs, accs)
@inline (rf::TeeRF)(accs, x) = next(rf, accs, x)

"""
    ProductRF(reducing_functions::Tuple)
    ProductRF(reducing_functions...)

Combine `N` reducing functions into a new reducing function that work
on `N`-tuple.  The `i`-th reducing function recieves the `i`-th
element of the input tuple.

Roughly speaking, `ProductRF(op₁, op₂, ..., opₙ)` is equivalent to

```julia
((a₁, a₂, ..., aₙ), (b₁, b₂, ..., bₙ)) -> (op₁(a₁, b₁), op₂(a₂, b₂), ..., opₙ(aₙ, bₙ))
```

!!! compat "Transducers.jl 0.4.32"

    New in version 0.4.32.

# Examples

Like [`TeeRF`](@ref), `ProductRF` can be used to drive multiple reducing
functions.  `ProductRF` is more "low-level" in the sense that `TeeRF`
can be defined in terms of `ProductRF` (other direction is much harder):

```jldoctest
julia> using Transducers

julia> TeeRF′(fs...) = reducingfunction(
           Map(x -> ntuple(_ -> x, length(fs))),
           ProductRF(fs...),
       );

julia> foldl(TeeRF′(min, max), Map(identity), [5, 2, 6, 8, 3])
(2, 8)
```

`ProductRF` may be useful for handling pre-existing stream whose item
type is already a tuple:

```julia; setup = (using Transducers)
julia> foldl(ProductRF(&, +), Map(x -> (isodd(x), x)), [5, 2, 6, 8, 3])
(false, 24)

julia> foldl(TeeRF(reducingfunction(Map(isodd), &), +), Map(identity), [5, 2, 6, 8, 3])
(false, 24)
```
"""
struct ProductRF{T<:Tuple} <: AbstractMultiCastingRF
    fs::T
    ProductRF{T}(fs) where {T} = new{T}(fs)
end

ProductRF(fs::Tuple) = ProductRF{typeof(fs)}(fs)
ProductRF(f, fs...) = ProductRF((f, fs...))

# Default `start` is defined by AbstractMultiCastingRF
function start(rf::ProductRF, inits::Tuple)
    _check_multirf_inits(rf.fs, inits)
    return map((f, i) -> start(f, i), rf.fs, inits)
end

_check_multirf_inits(::NTuple{N,Any}, ::NTuple{N,Any}) where {N} = nothing
function _check_multirf_inits(fs, inits)
    @nospecialize
    throw(ArgumentError("incompatible length of tuple is provided"))
end

@inline next(rf::ProductRF{<:NTuple{N,Any}}, accs::NTuple{N,Any}, xs) where {N} =
    map(next, rf.fs, accs, xs)
@inline (rf::ProductRF)(accs, x) = next(rf, accs, x)


text:
# ---
# title: 971. Flip Binary Tree To Match Preorder Traversal
# id: problem971
# author: Indigo
# date: 2021-02-20
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/flip-binary-tree-to-match-preorder-traversal/description/>
# hidden: true
# ---
#
# Given a binary tree with `N` nodes, each node has a different value from `{1,
# ..., N}`.
#
# A node in this binary tree can be _flipped_  by swapping the left child and
# the right child of that node.
#
# Consider the sequence of `N` values reported by a preorder traversal starting
# from the root.  Call such a sequence of `N` values the  _voyage_  of the tree.
#
# (Recall that a _preorder traversal_  of a node means we report the current
# node's value, then preorder-traverse the left child, then preorder-traverse
# the right child.)
#
# Our goal is to flip the **least number** of nodes in the tree so that the
# voyage of the tree matches the `voyage` we are given.
#
# If we can do so, then return a list of the values of all nodes flipped.  You
# may return the answer in any order.
#
# If we cannot do so, then return the list `[-1]`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-01.png)**
#
#
#
#     Input: root = [1,2], voyage = [2,1]
#     Output: [-1]
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-02.png)**
#
#
#
#     Input: root = [1,2,3], voyage = [1,3,2]
#     Output: [1]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/02/1219-02.png)**
#
#
#
#     Input: root = [1,2,3], voyage = [1,2,3]
#     Output: []
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 100`
#
#
## @lc code=start
using LeetCode

function flip_match_voyage(root::TreeNode, voyage::Vector{Int})
    res = Int[]
    i = 1
    dfs(::Nothing) = nothing
    function dfs(node::TreeNode)
        if node.val != voyage[i]
            res = [-1]
            return
        end
        i += 1
        if i ≤ length(voyage) && node.left !== nothing && node.left.val != voyage[i]
            push!(res, node.val)
            dfs(node.right)
            dfs(node.left)
        else
            dfs(node.left)
            dfs(node.right)
        end
    end
    dfs(root)
    return (!isempty(res) && res[1] == -1) ? [-1] : res
end
## @lc code=end