text:
module BouncingBeamsSimulation

using  ModiaLang
import Modia3D
using  Modia3D.ModiaInterface
using  Unitful

vmatRed   = Modia3D.VisualMaterial(color="Red")
vmatGreen = Modia3D.VisualMaterial(color="Green")
vmatBlue  = Modia3D.VisualMaterial(color="Blue")
vmatGrey  = Modia3D.VisualMaterial(color="Grey", transparency=0.5)

BouncingBeams = Model(
    boxHeigth = 0.1,
    world = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, 0, -1]),
                                   visualizeFrames=false,
                                   defaultFrameLength=0.2,
                                   visualizeBoundingBox = true,
                                   enableContactDetection=true,
                                   visualizeContactPoints=false)),
    worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
    ground = Object3D(parent=:world,
                      translation=:[0.0, 0.0, -boxHeigth/2],
                      feature=Solid(shape=Box(lengthX=5.0, lengthY=3.0, lengthZ=:boxHeigth),
                                    visualMaterial=vmatGrey,
                                    solidMaterial="DryWood",
                                    massProperties=MassPropertiesFromShape(),
                                    collision=true)),
    frameX = Object3D(parent=:world,
                      translation=:[-1.0, 0.0, 1.0],
                      rotation=:[-90*u"°", 0.0, -90*u"°"],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    beamX = Object3D(feature=Solid(shape=Beam(axis=1, length=1.0, width=0.4, thickness=0.2),
                                   visualMaterial=vmatRed,
                                   solidMaterial="DryWood",
                                   massProperties=MassPropertiesFromShape(),
                                   collision=true)),
    jointX = FreeMotion(obj1=:frameX, obj2=:beamX,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[0.0, -50*u"°", 0.0]),
                        v=Var(init=[0.0, 2.0, 0.0]),
                        w=Var(init=[2.0, 0.0, 0.0])),
    frameY = Object3D(parent=:world,
                      translation=:[0.0, 0.0, 1.0],
                      rotation=:[90*u"°", 90*u"°", 0.0],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    beamY = Object3D(feature=Solid(shape=Beam(axis=2, length=1.0, width=0.4, thickness=0.2),
                                   visualMaterial=vmatGreen,
                                   solidMaterial="DryWood",
                                   massProperties=MassPropertiesFromShape(),
                                   collision=true)),
    jointY = FreeMotion(obj1=:frameY, obj2=:beamY,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[0.0, 0.0, -50*u"°"]),
                        v=Var(init=[0.0, 0.0, 2.0]),
                        w=Var(init=[0.0, 2.0, 0.0])),
    frameZ = Object3D(parent=:world,
                      translation=:[1.0, 0.0, 1.0],
                      feature=Visual(shape=CoordinateSystem(length=0.5))),
    beamZ = Object3D(feature=Solid(shape=Beam(axis=3, length=1.0, width=0.4, thickness=0.2),
                                   visualMaterial=vmatBlue,
                                   solidMaterial="DryWood",
                                   massProperties=MassPropertiesFromShape(),
                                   collision=true)),
    jointZ = FreeMotion(obj1=:frameZ, obj2=:beamZ,
                        r=Var(init=[0.0, 0.0, 0.0]),
                        rot=Var(init=[-50*u"°", 0.0, 0.0]),
                        v=Var(init=[2.0, 0.0, 0.0]),
                        w=Var(init=[0.0, 0.0, 2.0]))
)

bouncingBeams = @instantiateModel(buildModia3D(BouncingBeams), unitless=true, log=false, logStateSelection=false, logCode=false)

#@show bouncingBeams.parameterExpressions
#@show bouncingBeams.parameters

stopTime = 1.2
tolerance = 1e-8
requiredFinalStates = [-0.8823695490633129, 1.665391474522192, 0.04799938397293286, -1.0204351804781866, 0.43467213487238104, -0.05300623034434308, 4.594448165765738, 1.48458962376636, -1.5821822875254883, -8.616696059760667, -0.38077326260409977, 0.4606687533750722, 0.04448291563173141, -0.8635203337074655, 1.6627077970537474, -0.049200694422998424, -0.8870648971101465, 0.40291237892093623, 3.2801065817877686, 0.1345412028139165, -1.597965173809505, 0.541296269597139, -7.6554365382890746, -0.5338918826110308, 1.9235100250337558, 0.10167700504201196, -1.3991207821301785, 1.5242500783672102, 0.07788755980354922, -3.4667039460822022, -2.7882122968930343, 0.12494117610295302, 4.331872438738044, 1.9861239477070278, -2.898970671344765, 1.0516923739870323]
simulate!(bouncingBeams, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(bouncingBeams, [("jointX.r", "jointY.r", "jointZ.r") ("jointX.rot", "jointY.rot", "jointZ.rot")
                     ("jointX.v", "jointY.v", "jointZ.v") ("jointX.w", "jointY.w", "jointZ.w")], figure=1)

end


text:
export TrainingHistory
export @push!, @enumerate, @get, @length

using DataStructures

import Base: length, push!, get, enumerate

# ==========================================================================

if VERSION < v"0.4-"
  typealias TupleType{T,V} (T, V)
else
  typealias TupleType{T,V} Tuple{T, V}
end

type TrainingHistory{T<:Integer}
  _storage::Dict{Symbol, Queue}
end

function TrainingHistory{T<:Integer}(::Type{T} = Int64)
  TrainingHistory{T}(Dict{Symbol, Queue}())
end

# ==========================================================================
# Functions

function length{T<:Integer}(history::TrainingHistory{T}, key::Symbol)
  length(history._storage[key])
end

function push!{T<:Integer,V<:Any}(history::TrainingHistory{T},
                                  iteration::T,
                                  key::Symbol,
                                  value::V)
  lastiter = zero(T)
  if !haskey(history._storage, key)
    iteration >= lastiter || throw(ArgumentError("Iterations must be greater than or equal to 0"))
    history._storage[key] = Queue(TupleType{T,V})
  else
    lastiter, _ = back(history._storage[key])
    iteration > lastiter || throw(ArgumentError("Iterations must increase over time"))
  end
  enqueue!(history._storage[key], (iteration, value))
  value
end

function enumerate{T<:Integer}(history::TrainingHistory{T}, key::Symbol)
  history._storage[key].store
end

function get{T<:Integer}(history::TrainingHistory{T}, key::Symbol)
  l = length(history, key)
  k, v = front(history._storage[key])
  karray = zeros(T, l)
  varray = Array(typeof(v), l)
  i = 1
  for (k, v) in enumerate(history, key)
    karray[i] = k
    varray[i] = v
    i += 1
  end
  karray, varray
end

# ==========================================================================
# Convenience macros

macro push!(history, iteration, func)
  key = string(func)
  esc(:(push!($history, $iteration, symbol($key), $func)))
end

macro enumerate(history, func)
  key = string(func)
  esc(:(enumerate($history, symbol($key))))
end

macro get(history, func)
  key = string(func)
  esc(:(get($history, symbol($key))))
end

macro length(history, func)
  key = string(func)
  esc(:(length($history, symbol($key))))
end

ext: jl


text:
using ArviZ
using Test

include("helpers.jl")
include("test_dataset.jl")
include("test_data.jl")
include("test_stats.jl")
include("test_mcmcchains.jl")


text:
male_first_name()::String = executor(data["faker"]["name"]["male_first_name"])
female_first_name()::String = executor(data["faker"]["name"]["female_first_name"])
last_name()::String = executor(data["faker"]["name"]["last_name"])
prefixe()::String = executor(data["faker"]["name"]["prefix"])
job()::String = executor(data["faker"]["name"]["title"]["job"])
blood_type()::String= string(rand(["A", "B", "AB", "O"]),rand(["+", "-"]))

function first_name(genere::String="None")::String

 if(cmp(genere, "M") == 0)
 	return male_first_name()
 elseif(cmp(genere, "F") == 0)
 	return female_first_name()
 end

 return executor(data["faker"]["name"]["first_name"])
end

function seleccion_field(data::Dict, fields::Tuple)::Dict
	e = Dict()
	for (k, v) in data
		if (k in fields)
			push!(e, k => v)
		end
	end
	return e
end

function simple_profile(fields...="None")
	fields[1] == "None" && (fields = []);
	sex::String = rand(["M","F"])
	data = Dict(
		"username"=> user_name(),
		"name"=> first_name(sex),
		"sex"=> sex,
		"address"=> address(),
		"mail"=> free_email(),
		"birthdate"=> date() )

	if length(fields) > 0
		return  seleccion_field(data, fields)
	else
		return data
	end
end

function profile( fields...="None")
	fields[1] == "None" && (fields = []);
	sex::String = rand(["M","F"])
	d = Dict(
		"username"=> user_name(),
		"name"=> first_name(sex),
		"sex"=> sex,
		"address"=> address(),
		"mail"=> free_email(),
		"birthdate"=> date(),
		"job"=> job(),
		"company"=> company(),
		"residence"=> address(),
		"current_location"=> (latitude(), longitude()),
		"blood_group"=> blood_type(),
		"website"=> url()
		);
	if length(fields) > 0
		return  seleccion_field(d, fields)
	else
		return d
	end
end
#stringWithSpaces
#

user_profile(items::Tuple=(user_name(),first_name(),rand(["M","F"]),address(),free_email(),date(),job(),company(),address(),string(latitude(), longitude()),blood_type(),url()))::String = reduce(stringWithComa, items )


text:
cd("/home/kkaloudis/Documents/Mathematics/PhD/Julia/Prediction/fq pred")

using Distributions
using Polynomials
using DataFrames

include("toolBox.jl")
include("fqPredictionGSB.jl") # GSB
include("fqPredictionPar.jl") # Parametric

#########################################################
# Parameter Declaration
#########################################################

# map parameters
theta, x0, y0 = [1.38, 0., 0.27, 0., -1., 0.], -0.5, 1.2
# theta, x0, y0 = [1, 0., 0.3, 0., -1.4, 0.], 0.7, 0.7
# theta, x0, y0 = [1, 0., 0.23, 0., -1.31, 0.], 0.7, 0.5
# theta, x0, y0 = [1., 0., 0., 0., -1.8, 0.], 0.7, 0.5

# sample size
ss = 200
# prediction horizon
T = 25

# noise f_{2,l} precisions
w1, lam1 = 0.7, 1e07
lam2 = 1e02

# corresponding deterministic orbit
xdet = zeros(ss+T)
xdet[1] = henon(theta, x0, y0)
xdet[2] = henon(theta, xdet[1], x0)
for i in 3:(ss+T)
  xdet[i] = henon(theta, xdet[i-1], xdet[i-2])
end

# dataSeed = 122 #f_{1}
# dataSeed = 46512 #f_{2,1}
dataSeed = 202011 #f_{2,2}
# dataSeed = 4517 #f_{2,3}
# dataSeed = 20205 #f_{2,4}

# fulldata = copy(xdet) # no noise
# fulldata = genDataf1(ss + T, theta, x0, y0, lam2, dataSeed) #f_{1} noise
fulldata = genData(ss + T, theta, x0, y0, w1, lam1, lam2, dataSeed) #f_{2,l} noise

predValues = copy(fulldata[(ss+1):(ss+T)])
data = copy(fulldata[1:ss])

# GSBR reconstruction - Prediction
gibbsIter, burnIn, gEps1, gEps2, thLow, thUp, zLow, zUp, papr, pbpr, thin, samplerSeed, filename =
 500000, 100000, 1e-03, 1e-03, -10.0, 10.0, -2.0, 2.0, 0.5, 0.5, 5, 12345, "/Results";

savelocation = string(pwd(), filename, "/GSB/seed$samplerSeed")
mkpath(savelocation)
writedlm(string(savelocation,"/predValues.txt"), predValues)
writedlm(string(savelocation,"/data.txt"), data, '\n')
@time thetas = fqPredictionGSB(data, T, gibbsIter, burnIn, gEps1, gEps2, thLow, thUp, zLow, zUp, papr, pbpr, thin, samplerSeed);
thetas = 0
gc()

# Parametric reconstruction - Prediction
fulldata = genData(ss + T, theta, x0, y0, w1, lam1, lam2, dataSeed) #f_{2,l} noise
predValues = copy(fulldata[(ss+1):(ss+T)])
data = copy(fulldata[1:ss])
savelocation = string(pwd(), filename, "/Par/seed$samplerSeed")
mkpath(savelocation)
writedlm(string(savelocation,"/predValues.txt"), predValues)
writedlm(string(savelocation,"/data.txt"), data, '\n')
@time thetas = fqPredictionPar(data, T, gibbsIter, burnIn, gEps1, gEps2, thLow, thUp, zLow, zUp, thin, samplerSeed);
thetas = 0
gc()

max_stars_repo_path: Complete quadratic/runCode.jl

max_stars_repo_name: KosKal/GSBR-DPR

max_issues_repo_path: Complete quadratic/runCode.jl

max_issues_repo_name: KosKal/GSBR-DPR

max_forks_repo_path: Complete quadratic/runCode.jl

max_forks_repo_name: KosKal/GSBR-DPR


text:
include("rachford_rice.jl")
include("bubble_point.jl")
include("dew_point.jl")
include("azeotrope_point.jl")
include("LLE_point.jl")
include("VLLE.jl")
include("crit_mix.jl")
include("UCST_mix.jl")
include("UCEP.jl")
include("tp_flash.jl")
export bubble_pressure,    dew_pressure,    LLE_pressure,    azeotrope_pressure, VLLE_pressure
export bubble_temperature, dew_temperature, LLE_temperature, azeotrope_temperature, VLLE_temperature
export crit_mix, UCEP_mix, UCST_mix


text:
using CUJulia.CUDNN

@testset "functions" for i = 1:5

T = Float32

# activation
x = CuArray(rand(T,10,5,4,3))
for mode in (CUDNN_ACTIVATION_SIGMOID,CUDNN_ACTIVATION_RELU)
    y, df = CUDNN.activation(mode, x)
    @test true
end

# convolution
x = CuArray(rand(T,5,4,3,2))
w = CuArray(rand(T,2,2,3,4))
b = CuArray(rand(T,5,4))
CUDNN.convolution(x, w, b, (0,0), (1,1))
@test true

# dropout
x = CuArray(rand(T,10,5,4,3))
y, df = CUDNN.dropout(x, 0.5)
@test true

# pooling
x = CuArray(rand(T,5,4,3,10))
y, df = CUDNN.pooling(CUDNN_POOLING_MAX, (2,2), (1,1), (2,2), x)
@test true

# softmax
x = CuArray(rand(T,5,4,3,10))
y, df = CUDNN.softmax(x)
y, df = CUDNN.logsoftmax(x)

end


text:
module RangeIndexingArrays

using Statistics
using DataFrames
using CSV

export RangeIndexingArray

"""
RangeIndexingArray
"""
struct RangeIndexingArray{T, N}
    # the septa values defining the indexing ranges
    septaListTuple::NTuple{N}
    # the table to be indexed
    table::AbstractArray{T, N}
end

"""
    range_string_list2septa_list(rangeStringList::Vector)
transform list of string to a list of septas

Example:
    input:  3-element Array{Union{Missing, String},1}:
            [ "(0,0.75]",  "(0.75,1.5]",  "(1.5,2]" ]

    output: 4-element Array{Float32,1}:
            [  0.0 ,    0.75,    1.5,    2.0 ]
"""
function range_string_list2septa_list(rangeStringList::Vector)
    septaList = Vector{Float32}(undef, length(rangeStringList)+1)
    @inbounds for (i, distRangeString) in rangeStringList |> enumerate
        distRangeString = replace(distRangeString, "("=>"")
        distRangeString = replace(distRangeString, "]"=>"")
        start, stop = map(x->Meta.parse(x)|>Float32, split(distRangeString, ","))
        septaList[i] = start
        if i == length(rangeStringList)
            septaList[end] = stop
        end
    end
    # make sure that the septa were sorted assendingly
    # so we can use binary search algorithm directly
    @assert all(diff(septaList) .> zero(Float32))
    septaList
end

@inline function range_string_list2septa_list(rangeStringList::CSV.Column{String, String})
    range_string_list2septa_list(map(string, rangeStringList))
end

"""
    RangeIndexingArray{T}(df::DataFrame)

The dataframe is a matrix, so the result is 2D
"""
function RangeIndexingArray{T}(df::DataFrame) where {T}
    table = convert( Matrix{T}, df[1:end, 2:end] )

    septaList1 = range_string_list2septa_list(df[:Column1] |> Vector{String})
    rangeStringList = map(string, names(df)[2:end] )
    septaList2 = range_string_list2septa_list( rangeStringList )

    RangeIndexingArray{T, 2}(tuple(septaList1, septaList2), table)
end

"""
    RangeIndexingArray(; fileName::AbstractString=joinpath(@__DIR__, "../asset/zfish_score_table.csv"))

read default zfish score table
"""
function RangeIndexingArray{T}(; fileName::AbstractString=joinpath(@__DIR__, "../../asset/zfish_score_table.csv")) where {T}
    df = CSV.read(fileName)
    RangeIndexingArray{T}(df)
end

"""
    DataFrame(ria::RangeIndexingArray{T,N})
transform a RangeIndexingArray to DataFrame
"""
function DataFrame(ria::RangeIndexingArray{T,N}) where {T,N}
    distRangeList = Vector{Symbol}()
    @inbounds for i in 1:length(ria.septaListTuple[1])-1
        start = ria.septaListTuple[1][i]
        stop = ria.septaListTuple[1][i+1]
        range = Symbol("($(start),$(stop)]")
        push!(distRangeList, range)
    end

    df = DataFrame(Column1=distRangeList)

    @inbounds for i in 1:length(ria.septaListTuple[2])-1
        start = ria.septaListTuple[2][i]
        stop = ria.septaListTuple[2][i+1]
        range = Symbol("($(start),$(stop)]")
        df[range] = ria.table[:,i]
    end

    df
end

function Base.eltype(self::RangeIndexingArray{T,N}) where {T,N}
    return T
end

@inline function get_septa_list_tuple(self::RangeIndexingArray)
    self.septaListTuple
end

@inline function get_table_array(self::RangeIndexingArray)
    self.table
end

"""
    float_index2table_index(floatIndex::T, septaList::Vector{T})
transform index with type of floating point to the integer index in the table
To-Do: use the [Binary Search Algorithm](https://en.wikipedia.org/wiki/Binary_search_algorithm)
"""
function _float_index2table_index(floatIndex::T, septaList::Vector{T}) where T
    if floatIndex<septaList[1]
        return 1
    end
    @inbounds for i in 1:length(septaList)-1
        if floatIndex <= septaList[i+1]
            return i
        end
    end
    # this is even larger than the last septa value
    return length(septaList)-1
end

@inline function Base.fill!(self::RangeIndexingArray{T,N}, x::T) where {T,N}
    fill!(self.table, x)
end

@inline function Base.lastindex(self::RangeIndexingArray{T,N}, dim::Integer) where {T,N}
    table = get_table_array(self)
    lastindex(table, dim)
end

function Base.getindex(self::RangeIndexingArray{T,N}, idxes::Integer...) where {T,N}
    table = get_table_array(self)
    return table[idxes...]
end

"""
    Base.getindex(self::RangeIndexingArray{T,N}, idx::T...)

given a index which could be float type, return a value in the data array.
"""
function Base.getindex(self::RangeIndexingArray{T,N}, floatIndexList::AbstractFloat...) where {T,N}
    tableIndexList = map((x,y)->_float_index2table_index(x,y),
                         floatIndexList, get_septa_list_tuple(self))
    table = get_table_array(self)

#    for (i,index) in tableIndexList |> enumerate
#        if index<1
#            # if they are pretty similar, return score of 1
#            return one(T)
#        end
#    end
#    for (i,index) in tableIndexList |> enumerate
#        if index>size(table, i)
#            # if they are very different, return score of 0
#            return zero(T)
#        end
#    end

    return get_table_array(self)[tableIndexList...,]
end

function Base.setindex!(self::RangeIndexingArray{T,N}, value::T,
                       inds::AbstractFloat...) where {T,N}
    tableIndexList = map((x,y)->_float_index2table_index(x,y), inds, get_septa_list_tuple(self))
    table = get_table_array(self)
    table[tableIndexList...,] = value
end

"""
    to_position_only(self::RangeIndexingArray{T,2})

convert to position only range indexing array. This is used for dendrite query.
The dendrite seems only respect approximity rather than direction.
"""
function to_position_only(self::RangeIndexingArray{T,2}) where T
    ret = deepcopy(self)
    v = mean(ret.table; dims=2) |> vec
    for j in size(ret.table, 2)
        ret.table[:,j] = v
    end
    @assert size(ret.table, 2) == 10
    ret
end

end # end of module


text:
using Dates

struct TenorPeriod{P <: Dates.Period}
  period::P
  freq::Frequency
end

# Constructors
function TenorPeriod(f::Frequency)
  freq = value(f)
  if freq < 0
    period = Dates.Day(0)
  elseif freq <= 1
    period = Dates.Year(freq)
  elseif freq <= 12
    period = Dates.Month(12 / freq)
  elseif freq <= 52
    period = Dates.Week(52/freq)
  else
    period = Dates.Day(1)
  end

  return TenorPeriod(period, f)
end

TenorPeriod(p::Dates.Year) = TenorPeriod(p, Annual())

function TenorPeriod(p::Dates.Month)
  x = Dates.value(p)
  if x == 6
    return TenorPeriod(p, Semiannual())
  elseif x == 3
    return TenorPeriod(p, Quarterly())
  elseif x == 4
    return TenorPeriod(p, EveryFourthMonth())
  elseif x == 2
    return TenorPeriod(p, Bimonthly())
  else
    return TenorPeriod(p, Monthly())
  end
end

function TenorPeriod(p::Dates.Week)
  x = Dates.value(p)
  if x == 26
    return TenorPeriod(p, Biweekly())
  elseif x == 13
    return TenorPeriod(p, EveryFourthWeek())
  else
    return TenorPeriod(p, Weekly())
  end
end


text:
## coefficient tables with specialized show method

mutable struct CoefTable
    cols::Vector
    colnms::Vector
    rownms::Vector
    pvalcol::Int
    teststatcol::Int
    function CoefTable(cols::Vector,colnms::Vector,rownms::Vector,
                       pvalcol::Int=0,teststatcol::Int=0)
        nc = length(cols)
        nrs = map(length,cols)
        nr = nrs[1]
        length(colnms) in [0,nc] || throw(ArgumentError("colnms should have length 0 or $nc"))
        length(rownms) in [0,nr] || throw(ArgumentError("rownms should have length 0 or $nr"))
        all(nrs .== nr) || throw(ArgumentError("Elements of cols should have equal lengths, but got $nrs"))
        pvalcol in 0:nc || throw(ArgumentError("pvalcol should be between 0 and $nc"))
        teststatcol in 0:nc || throw(ArgumentError("teststatcol should be between 0 and $nc"))
        new(cols,colnms,rownms,pvalcol,teststatcol)
    end

    function CoefTable(mat::Matrix,colnms::Vector,rownms::Vector,
                       pvalcol::Int=0,teststatcol::Int=0)
        nc = size(mat,2)
        cols = Any[mat[:, i] for i in 1:nc]
        CoefTable(cols,colnms,rownms,pvalcol,teststatcol)
    end
end

Base.length(ct::CoefTable) = length(ct.cols[1])
function Base.eltype(ct::CoefTable)
    names = isempty(ct.rownms) ?
        tuple(Symbol.(ct.colnms)...) :
        tuple(Symbol("Name"), Symbol.(ct.colnms)...)
    types = isempty(ct.rownms) ?
        Tuple{eltype.(ct.cols)...} :
        Tuple{eltype(ct.rownms), eltype.(ct.cols)...}
    NamedTuple{names, types}
end

function Base.iterate(ct::CoefTable, i::Integer=1)
    if i in 1:length(ct)
        cols = getindex.(ct.cols, Ref(i))
        nt = isempty(ct.rownms) ?
            eltype(ct)(tuple(cols...)) :
            eltype(ct)(tuple(ct.rownms[i], cols...))
        (nt, i+1)
    else
        nothing
    end
end

"""
Show a p-value using 6 characters, either using the standard 0.XXXX
representation or as <Xe-YY.
"""
struct PValue <: Real
    v::Real
    function PValue(v::Real)
        0 <= v <= 1 || isnan(v) || error("p-values must be in [0; 1]")
        new(v)
    end
end
PValue(p::PValue) = p

function show(io::IO, pv::PValue)
    v = pv.v
    if isnan(v)
        @printf(io,"%d", v)
    elseif v >= 1e-4
        @printf(io,"%.4f", v)
    else
        @printf(io,"<1e%2.2d", ceil(Integer, max(nextfloat(log10(v)), -99)))
    end
end

"""Show a test statistic using 2 decimal digits"""
struct TestStat <: Real
    v::Real
end

show(io::IO, x::TestStat) = @printf(io, "%.2f", x.v)
TestStat(x::TestStat) = x

float(x::Union{TestStat, PValue}) = float(x.v)

for op in [:(==), :<, :≤, :>, :≥, :(isless), :(isequal)] # isless and < to place nice with NaN
    @eval begin
        Base.$op(x::Union{TestStat, PValue}, y::Real) = $op(x.v, y)
        Base.$op(y::Real, x::Union{TestStat, PValue}) = $op(y, x.v)
        Base.$op(x1::Union{TestStat, PValue}, x2::Union{TestStat, PValue}) = $op(x1.v, x2.v)
    end
end

Base.hash(x::Union{TestStat, PValue}, h::UInt) = hash(x.v, h)

# necessary to avoid a method ambiguity with isless(::TestStat, NaN)
Base.isless(x::Union{TestStat, PValue}, y::AbstractFloat) = isless(x.v, y)
Base.isless(y::AbstractFloat, x::Union{TestStat, PValue},) = isless(y, x.v)
Base.isequal(y::AbstractFloat, x::Union{TestStat, PValue}) = isequal(y, x.v)
Base.isequal(x::Union{TestStat, PValue}, y::AbstractFloat) = isequal(x.v, y)

Base.isapprox(x::Union{TestStat, PValue}, y::Real; kwargs...) = isapprox(x.v, y; kwargs...)
Base.isapprox(y::Real, x::Union{TestStat, PValue}; kwargs...) = isapprox(y, x.v; kwargs...)
Base.isapprox(x1::Union{TestStat, PValue}, x2::Union{TestStat, PValue}; kwargs...) = isapprox(x1.v, x2.v; kwargs...)

"""Wrap a string so that show omits quotes"""
struct NoQuote
    s::String
end

show(io::IO, n::NoQuote) = print(io, n.s)

function show(io::IO, ct::CoefTable)
    cols = ct.cols; rownms = ct.rownms; colnms = ct.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    mat = [j == 1 ? NoQuote(rownms[i]) :
           j-1 == ct.pvalcol ? NoQuote(sprint(show, PValue(cols[j-1][i]))) :
           j-1 in ct.teststatcol ? TestStat(cols[j-1][i]) :
           cols[j-1][i] isa AbstractString ? NoQuote(cols[j-1][i]) : cols[j-1][i]
           for i in 1:nr, j in 1:nc+1]
    # Code inspired by print_matrix in Base
    io = IOContext(io, :compact=>true, :limit=>false)
    A = Base.alignment(io, mat, 1:size(mat, 1), 1:size(mat, 2),
                       typemax(Int), typemax(Int), 3)
    nmswidths = pushfirst!(length.(colnms), 0)
    A = [nmswidths[i] > sum(A[i]) ? (A[i][1]+nmswidths[i]-sum(A[i]), A[i][2]) : A[i]
         for i in 1:length(A)]
    totwidth = sum(sum.(A)) + 2 * (length(A) - 1)
    println(io, repeat('─', totwidth))
    print(io, repeat(' ', sum(A[1])))
    for j in 1:length(colnms)
        print(io, "  ", lpad(colnms[j], sum(A[j+1])))
    end
    println(io, '\n', repeat('─', totwidth))
    for i in 1:size(mat, 1)
        Base.print_matrix_row(io, mat, A, i, 1:size(mat, 2), "  ")
        i != size(mat, 1) && println(io)
    end
    print(io, '\n', repeat('─', totwidth))
    nothing
end

function show(io::IO, ::MIME"text/markdown", ct::CoefTable)
    cols = ct.cols; rownms = ct.rownms; colnms = ct.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    mat = [j == 1 ? NoQuote(rownms[i]) :
           j-1 == ct.pvalcol ? NoQuote(sprint(show, PValue(cols[j-1][i]))) :
           j-1 in ct.teststatcol ? TestStat(cols[j-1][i]) :
           cols[j-1][i] isa AbstractString ? NoQuote(cols[j-1][i]) : cols[j-1][i]
           for i in 1:nr, j in 1:nc+1]
    # Code inspired by print_matrix in Base
    io = IOContext(io, :compact=>true, :limit=>false)
    A = Base.alignment(io, mat, 1:size(mat, 1), 1:size(mat, 2),
                       typemax(Int), typemax(Int), 3)
    nmswidths = pushfirst!(length.(colnms), 0)
    A = [nmswidths[i] > sum(A[i]) ? (A[i][1]+nmswidths[i]-sum(A[i]), A[i][2]) : A[i]
         for i in 1:length(A)]
    totwidth = sum(sum.(A)) + 2 * (length(A) - 1)

    # not using Markdown stdlib here because that won't give us nice decimal
    # alignment (even if that is lost when rendering to HTML, it's still nice
    # when looking at the markdown itself)

    print(io, '|', ' '^(sum(A[1])+1))
    for j in 1:length(colnms)
        print(io, " | ", lpad(colnms[j], sum(A[j+1])))
    end

    println(io, " |")
    print(io, '|', rpad(':', sum(A[1])+2, '-'))
    for j in 1:length(colnms)
        _pad = j-1 in [ct.teststatcol; ct.pvalcol] ? rpad : lpad
        print(io, '|', _pad(':', sum(A[j+1])+2, '-'))
    end
    println(io, '|')

    for i in 1:size(mat, 1)
        print(io, "| ")
        Base.print_matrix_row(io, mat, A, i, 1:size(mat, 2), " | ")
        print(io, " |")
        i != size(mat, 1) && println(io)
    end

    nothing
end

"""
    ConvergenceException(iters::Int, lastchange::Real=NaN, tol::Real=NaN)

The fitting procedure failed to converge in `iters` number of iterations,
i.e. the `lastchange` between the cost of the final and penultimate iteration was greater than
specified tolerance `tol`.
"""
struct ConvergenceException{T<:Real} <: Exception
    iters::Int
    lastchange::T
    tol::T
    msg::String
    function ConvergenceException{T}(iters, lastchange::T, tol::T, msg::String) where T<:Real
        if tol > lastchange
            throw(ArgumentError("Change must be greater than tol."))
        else
            new(iters, lastchange, tol, msg)
        end
    end
end

ConvergenceException(iters, lastchange::T=NaN, tol::T=NaN,
                     msg::AbstractString="") where {T<:Real} =
    ConvergenceException{T}(iters, lastchange, tol, String(msg))

function Base.showerror(io::IO, ce::ConvergenceException)
    print(io, "failure to converge after $(ce.iters) iterations.")
    if !isnan(ce.lastchange)
        print(io, " Last change ($(ce.lastchange)) was greater than tolerance ($(ce.tol)).")
    end
    if !isempty(ce.msg)
        print(io, ' ', ce.msg)
    end
end


text:
using Treesitter

mutable struct Cursor
    current_cursor::Union{TSTreeCursor,Nothing}
    current_node::TSNode
    next_sibling::Union{TSTreeCursor,Nothing}
    first_child::Union{TSTreeCursor,Nothing}
    text::Array{UInt8,1}
    current_text
end # struct

parser = Parser()
set_language(parser, tree_sitter_c())
tree = Treesitter.parse_filewith(
    parser,
    joinpath(@__DIR__, "corpus/easystruct.h"),
)

root_cursor = tree_getroot(tree)

root = cursor_new(root_cursor)

cursor = Cursor(
    root,
    cursor_current_node(root),
    cursor_next(root),
    cursor_firstchild(root),
    parser.buffer,
    nothing,
)
cursor.current_text = view(
    cursor.text,
    node_startbyte(cursor.current_node)+1:node_endbyte(cursor.current_node),
)

# TODO change to use nodes insea
function byterange(node::TSNode)
    start = node_startbyte(node) + 1
    stop = node_endbyte(node)
    return start:stop
end

# TODO change to use nodes instead of cursors.
function cursor_next!(cursor::Cursor)
    isnothing(cursor.next_sibling) && return

    cursor_next(cursor.current_cursor)

    cursor.current_cursor = cursor.next_sibling
    cursor.current_node = cursor_current_node(cursor.current_cursor)
    cursor.next_sibling = cursor_next(cursor.current_cursor)

    bytes = byterange(cursor.current_node)
    cursor.current_text = view(cursor.text, bytes)
end

# TODO change to use nodes instead of cursors.
function cursor_firstchild!(cursor::Cursor)
    isnothing(cursor.first_child) && return

    cursor.current_cursor = cursor.first_child

    cursor.current_node = cursor_current_node(cursor.current_cursor)

    cursor.next_sibling = cursor_next(cursor.current_cursor)
    cursor.first_child = cursor_firstchild(cursor.current_cursor)

    cursor.current_text = view
end


text:
module SpirvCross

include(joinpath(@__DIR__, "..", "gen", "LibSpirvCross.jl"))
using .LibSpirvCross
const api = LibSpirvCross

include("wrap_utils.jl")
include("resources.jl")

export
    Context,
    IR,
    Compiler,
    ShaderResources,
    get_resource_list

end


text:
#=
learn a bit more about MCMC. This explores using
a proposal density that takes into account the
correlations of the posterior, to get better
mixing

* play with scaling to see effects on acceptance rate
* examine chains with different acceptance rates
* look at scatter plots to see posterior correlations
=#
using Plots, LinearAlgebra, Statistics, Econometrics
include("QIVmodel.jl")
function main()
LNW, X, Z = getdata()
n = size(LNW,1)
mcmcreps = 100000
burnin = 100000
# use st. errs. from ordinary QR as quide to setting tuning.
basetuning = [0.1, 0.05, 0.03, 0.005, 0.3, 0.2, 0.2]
# adjust this up or down to achieve decent acceptance rate
scale = 0.1 # higher for higher taus
tuning = scale*basetuning
τ = 0.5
Σ = τ*(1.0-τ)*Z'Z/n
Σinv = inv(Σ)
θ = X\LNW  # OLS start values
# to do ordinary QR via MCMC, set Z=X (just to verify that MCMC works!)
lnL = θ -> likelihood(θ, LNW, X, Z, τ, Σinv)
Proposal = θ -> proposal(θ, tuning)
Prior = θ -> prior(θ)
chain = mcmc(θ, mcmcreps, burnin, Prior, lnL, Proposal)
d = dstats(chain)
# explore things:
p1 = plot(chain[:,1], legend=false, title="simple proposal: theta1") # chain is highly autocorrelated!
p2 = scatter(chain[:,1],chain[:,2], legend=false, title="simple proposal: theta1 vs. theta2") # there is significant posterior correlation
#= these results suggest that the proposal could be improved:
    * compute the covariance of the chain, and make draws from random walk MVN with that covariance
    * the following implements this idea
=#

# the following uses a different proposal, with correlations
# across the parameters, based on the previous chain.
V = cov(chain[:,1:end-1])
cholV = cholesky(V).U
tuning = 0.5
Proposal = θ -> proposal2(θ, tuning*cholV)
chain = mcmc(θ, mcmcreps, burnin, Prior, lnL, Proposal)
p3 = plot(chain[:,1], legend=false, title="correlated proposal: theta1") # chain is highly autocorrelated!
p4 = scatter(chain[:,1],chain[:,2], legend=false, title="correlated proposal: theta1 vs. theta2") # there is significant posterior correlation
plot(p1,p2,p3,p4,layout=(2,2))
gui()
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using CSV
using DataFrames
#include("clustering_ws.jl")
include("Network.jl")

#Number of Iterations Benders
K=1:20
#Sets
#Investemnt Options
set_opt_thermalgenerators = 1:5
set_opt_winds = 1:4
#Existing Units
set_thermalgenerators = 1:10
set_winds = 1:2
set_demands = 1:17
set_nodes = 1:24
set_nodes_ref=1
set_nodes_noref=2:24
#Operating conditions-scenarios-clsuters
set_scenarios= 1:100
#set_scenarios= 1:3
set_times= 1:48

#Penalty
P=10000
V=10000
#Assign the values of demands according to their locations in the network
#3 operating conditions below
#two time zones below#Maybe no key means and just montecarlo simulation? I do need to select, let's say 876 (10%)
max_demand=400
R=0.15
#max_demand=150
#max_wind=10

#Scenarios
a=1/length(set_scenarios)
γ=zeros(length(set_scenarios))
for s in set_scenarios
γ[s]=a
end

network=hcat(From,To,LineCapacity,Susceptance)

start_node=network[:,1]
end_node=network[:,2]
F_max=network[:,3]
B=network[:,4]

n_node=24
n_link=68

links = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links[i] = (start_node[i], end_node[i])
end

set_number_links=1:n_link
number_links=[i for i in set_number_links]

Z=[number_links,links]

links_rev = Array{Tuple{Int64, Int64},2}(undef, n_link, 1)

for i=1:n_link
  links_rev[i] = (end_node[i], start_node[i])
end

F_max_dict=Dict()
B_dict=Dict()

for i=1:n_link
  F_max_dict[(start_node[i], end_node[i])] = F_max[i]
  B_dict[(start_node[i], end_node[i])] = B[i]
end

MapG=[(1,2),(2,7),(3,13),(4,15),(5,15),(6,16),(7,21),(8,22),(9,23),(10,23)]
MapG_opt=[(1,16),(2,17),(3,19),(4,21),(5,22)]
MapD=[(1,1),(2,2),(3,3),(4,4),(5,5),(6,6),(7,7),(8,8),(9,9),(10,10),(11,13),(12,14),(13,15),(14,16),(15,18),(16,19),(17,20)]

MapW=[(1,1),(2,18)]
MapW_opt=[(1,19),(2,21),(3,23),(4,8)]

#Main difference with GAMS is that I calculate the variable cost beforehand
#NREL 2020 Data

#Nuclear and Coal in this case
tech_thermal=[        1500      9.20       7185655         10000            1              400              118988       0                 10.461        0.66        60
                      1500      9.20       7185655         10000            1              400              118988       0                 10.461        0.66        60
                      1500      9.20       7185655         10000            1              400              118988       0                 10.461        0.66        60
                      1000      21.76      4149988         10000            1              300              39695        0.0955             8.638        2.01        75
                      1000      21.76      4149988         10000            1              300              39695        0.0955             8.638        2.01        75
                      500       31.24      982962          10000            1              250              11395        0.0531            9.5145        2.81        55
                      500       31.24      982962          10000            1              250              11395        0.0531            9.5145        2.81        55
                      500       20.19      1088328         10000            1              250              12863        0.0531            6.4005        2.81        55
                      500       20.19      1088328         10000            1              250              12863        0.0531            6.4005        2.81        55
                      500       20.19      1088328         10000            1              250              12863        0.0531            6.4005        2.81        55]

tech_thermal_opt=[   1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
                     1000      35.84        6830748        10000            1              0               58242        0.0096            12.507        2.01        75
                     500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
                     500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
                     500       26.85        2778138        10000            1              0               26994        0.0053             7.525        2.81        55]

tech_wind=[   300         0.0               1604886        10000            1             200              43205        0                    0          0.00        30
              300         0.0               1604886        10000            1             200              43205        0                    0          0.00        30]

tech_wind_opt=[300         0.0              1604886        300               1             0               43205        0                    0          0.00        30
               300         0.0              1604886        300               1             0               43205        0                    0          0.00        30
               300         0.0              1604886        300               1             0               43205        0                    0          0.00        30
               300         0.0              1604886        300               1             0               43205        0                    0          0.00        30]

 #Indices for set_tech
capacity_per_unit=1
varcost=2
invcost=3
maxBuilds=4
ownership=5
capacity_existingunits=6
fixedcost=7
EmissionsRate=8
HeatRate=9
fuelprice=10
life=11

WACC=0.08

CRF_thermal_opt=zeros(5)

for i in 1:5
CRF_thermal_opt[i]=WACC/(1 - (1 / ( (1 + WACC)^tech_thermal_opt[i,life])))
end

CRF_wind_opt=zeros(4)

for i in 1:4
CRF_wind_opt[i]=WACC/(1 - (1 / ( (1 + WACC)^tech_wind_opt[i,life])))
end

α_down=-10000

Set_α=zeros(length(K))
Set_master_obj=zeros(length(K))
#Set_stage1_cost=zeros(length(K))
Set_invements_thermal=SharedArray(zeros(length(set_opt_thermalgenerators),length(K)))
Set_invements_wind=SharedArray(zeros(length(set_opt_winds),length(K)))

Set_investment_cost=zeros(length(K))
Set_lowerbound=zeros(length(K))

subproblem_obj_sce=SharedArray(zeros(length(set_scenarios),length(K)))
Set_stage2_cost=zeros(length(K))
Set_upperbound=zeros(length(K))
Set_Dual_constraint3_st2=SharedArray(zeros(length(set_thermalgenerators),length(set_times),length(set_scenarios),length(K)))
Set_Dual_constraint4_st2=SharedArray(zeros(length(set_opt_thermalgenerators),length(set_times),length(set_scenarios),length(K)))
Set_Dual_constraint5_st2=SharedArray(zeros(length(set_winds),length(set_times),length(set_scenarios),length(K)))
Set_Dual_constraint6_st2=SharedArray(zeros(length(set_opt_winds),length(set_times),length(set_scenarios),length(K)))
Set_Dual_constraint8_st2=SharedArray(zeros(n_link,length(set_times),length(set_scenarios),length(K)))
Set_Dual_constraint10_st2=SharedArray(zeros(length(set_nodes),length(set_times),length(set_scenarios),length(K)))

#Compute Var Cost
#2020 Data

#=

tech=[1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
      1000      35.84        6830748        10000            1              0               58242        0.0096            12.507        2.01        75
      500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
      500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
      500       26.85        2778138        10000            1              0               26994        0.0053             7.525        2.81        55
      300       0.0          1604886        10000            1              0               43205        0                    0          0.00        30
      300       0.0          1599902        10000            1              0               18760        0                    0          0.00        30
      1500      9.20         7185655        10000            1             10               118988       0                 10.461        0.66        60
      1000      21.76        4149988        10000            1              8               39695        0.0955             8.638        2.01        75]

TABLE Tech(q,GEN_PARAMS)
                    capacity   var_cost      invcost       maxBuilds      ownership   existingunits        fixedcost    EmissionsRate     HeatRate      FPPRICE     LIFE
*                    [MW]      [$/MWh]       [$/MW]   [Number of Units]                                   [$/MW-yr]    [ton/MMBTU]        [MMBTU/MWh] [$/MMBTU]   [years]
Coal_CCS_30_1        1000      26.50        5299096        10000            1              0               53114        0.0096            9.751         2.01        75
Coal_CCS_90_1        1000      35.84        6830748       10000            1              0               58242        0.0096            12.507        2.01        75
Gas_CT_1             500       31.24        982962         10000            1              0               11395        0.0531            9.5145        2.81        55
Gas_CC_1             500       20.19        1088328        10000            1              0               12863        0.0531            6.4005        2.81        55
Gas_CC_CCS_1         500       26.85       2778138        10000            1              0               26994        0.0053             7.525        2.81        55
Wind_onshore_1       300         0.0        1604886        10000            1              0               43205        0                    0          0.00        30
Solar_Utility_PV_1   300         0.0        1599902        10000            1              0               18760        0                    0          0.00        30
Nuclear_1            1500      9.20       7185655        10000            1             10               118988       0                 10.461        0.66        60
Coal_new_1           1000      21.76       4149988        10000            1              8                39695       0.0955             8.638        2.01        75
;

=#


text:
@recipe function plot(tracks::TracksAhead, camera::VideoCamera; sleeper_displace = 0.0, camera_grid = true)

    trackproperties = tracks.trackproperties

    # calculate tracks in the camera image
    uvLs = camera_image(camera, tracks.left_track);
    uvRs = camera_image(camera, tracks.right_track);

    uLs = [uv[1] for uv in uvLs];
    uRs = [uv[1] for uv in uvRs];
    vLs = [uv[2] for uv in uvLs];
    vRs = [uv[2] for uv in uvRs];

    uvL2s = camera_image(camera, [[t[1],t[2] - trackproperties.track_width, t[3]] for t in tracks.left_track]);
    uvR2s = camera_image(camera, [[t[1],t[2] + trackproperties.track_width, t[3]] for t in tracks.right_track]);

    uL2s = [uv[1] for uv in uvL2s];
    uR2s = [uv[1] for uv in uvR2s];
    vL2s = [uv[2] for uv in uvL2s];
    vR2s = [uv[2] for uv in uvR2s];

    Ymax = 0.8 * trackproperties.track_gauge;
    Xmax = maximum(t[1] for t in tracks.left_track) / 2 + maximum(t[1] for t in tracks.right_track) / 2
    Xmin = minimum(t[1] for t in tracks.left_track) / 2 + minimum(t[1] for t in tracks.right_track) / 2
    # Xmin = 0.0

    # calculate sleepers in the camera image
    sw = trackproperties.sleeper_width
    sd = trackproperties.sleeper_distance

    left_distances = cumsum(norm.(circshift(tracks.left_track,1)[2:end] - tracks.left_track[2:end]))
    left_distances = [0.0; left_distances]

    right_distances = cumsum(norm.(circshift(tracks.right_track,1)[2:end] - tracks.right_track[2:end]))
    right_distances = [0.0; right_distances]

    no = - floor(sleeper_displace/(sd + sw))
    sleeper_number = -1 + Int(min(
        - no + floor((left_distances[end] - sleeper_displace) / (sd + sw)),
        - no + floor((right_distances[end] - sleeper_displace) / (sd + sw))
    ))

    sleepers_uv = map(0:sleeper_number) do i
        t = findfirst((i+no) * (sd + sw) + sleeper_displace .<= left_distances)
        if t == length(left_distances) t = t - 1 end

        v = tracks.left_track[t+1] - tracks.left_track[t];
        v = v / norm(v);

        p1 =  tracks.left_track[t] + ((i+no) * (sd + sw) + sleeper_displace - left_distances[t]) .* v
        p2 =  tracks.left_track[t] + ((i+no) * (sd + sw) + sw + sleeper_displace - left_distances[t]) .* v

        t = findfirst((i+no) * (sd + sw) + sleeper_displace .< right_distances)
        if t == length(left_distances) t = t - 1 end

        v = tracks.right_track[t+1] - tracks.right_track[t];
        v = v / norm(v);

        p3 =  tracks.right_track[t] + ((i+no) * (sd + sw) + sw + sleeper_displace - right_distances[t]) .* v
        p4 =  tracks.right_track[t] + ((i+no) * (sd + sw) + sleeper_displace - right_distances[t]) .* v

        sleeper_XYZ = [p1,p2,p3,p4]

        camera_image(camera,sleeper_XYZ)
    end

    # Plot the camera grid
    if camera_grid
        @series begin
            Xmax --> Xmax
            Xmin --> Xmin
            Ymax --> Ymax
            camera
        end
    end

    # plot the sleepers
    for uvs in sleepers_uv
        @series begin
            linecolor --> :black
            fillcolor --> :brown
            fillalpha --> 0.6
            linewidth --> 0.1
            label --> ""
            Shape([uv[1] for uv in uvs], [uv[2] for uv in uvs])
        end
    end

    # plot the left tracks
    @series begin
        linecolor --> :black
        fillcolor --> :red
        fillalpha --> 0.8
        linewidth --> 0.1
        label --> ""
        Shape([uLs;reverse(uL2s)], [vLs;reverse(vL2s)])
    end

    # plot the right tracks
    @series begin
        yflip := true
        grid --> false
        xguide --> "u"
        yguide --> "v"
        aspect_ratio := 1.0
        label --> ""
        linecolor --> :black
        fillcolor --> :red
        fillalpha --> 0.8
        linewidth --> 0.1
        Shape([uRs;reverse(uR2s)], [vRs;reverse(vR2s)])
    end

end

@recipe function plot(camera::VideoCamera; Xdisplace = 0.0, Xmin = 1.0, Xmax = 20.0, Ymax = 1.5, Ymin = - Ymax, mesh_number = 6)

    Yspace = LinRange(Ymin,Ymax,mesh_number)
    dy = (Yspace[2] - Yspace[1])
    no = - floor(Xdisplace/dy)

    # Xspace = [reverse(min(Xmax,Xmin+Xdisplace):-(Yspace[2] - Yspace[1]):Xmin); (Xmin+Xdisplace):(Yspace[2] - Yspace[1]):Xmax ]
    Xspace = (Xmin + Xdisplace + no*dy):dy:Xmax

    hlines = [ [[x, y, 0.0] for y in Yspace] for x in Xspace];
    vlines = [ [[x, y, 0.0] for x in Xspace] for y in Yspace];

    cam_hlines = map(hlines) do h
        uvs = camera_image(camera, h)
        us = [uv[1] for uv in uvs]
        vs = [uv[2] for uv in uvs]
        [us, vs]
    end;
    cam_vlines = map(vlines) do v
        uvs = camera_image(camera, v)
        us = [uv[1] for uv in uvs]
        vs = [uv[2] for uv in uvs]
        [us, vs]
    end;

    for cam in [cam_hlines; cam_vlines]
        @series begin
            grid --> false
            xguide --> "u"
            yguide --> "v"
            aspect_ratio := 1.0
            label --> ""
            linecolor --> :gray
            linealpha --> 0.3
            linewidth --> 0.5
            yflip := true

            (cam[1],cam[2])
        end
    end
end

# pyplot(size = (400,300))
# plot()
# [plot!(cam[1],cam[2],grid=false, linecolor = :gray,label = "",linealpha = 0.3, linewidth = 0.5) for cam in [cam_hlines; cam_vlines]];
# gui()


text:
using Community
using PeaceVote
using Synchronizers: Ledger

setnamespace(@__MODULE__)

import PeaceVote.voters!

ledger = Ledger(tempdir() * "/ledger/")
sync!(ledger)

messages = Community.braidchain(ledger)

members = Community.members(messages)

voters = Set()
voters!(voters,messages)

for proposal in PeaceVote.proposals(messages)
    tally = count(proposal,messages)
    @show proposal,tally
end

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
# NFFT utilities

export NFFTLinOp, NFFTParametericLinOp, NFFTParametericDelayedEval, JacobianNFFTEvaluated
export nfft_linop, nfft, Jacobia, ∂, sparse_matrix_GaussNewton

## NFFT linear operator

struct NFFTLinOp{T}<:AbstractLinearOperator{Complex{T},3,2}
    X::RegularCartesianSpatialSampling{T}
    K::AbstractKSpaceSampling{T}
    tol::T
end

nfft_linop(X::RegularCartesianSpatialSampling{T}, K::AbstractKSpaceSampling{T}; tol::T=T(1e-6)) where {T<:Real} = NFFTLinOp{T}(X, K, tol)

nfft_linop(X::RegularCartesianSpatialSampling{T}; phase_encoding::NTuple{2,Integer}=(1,2), subsampling::Union{Nothing,AbstractVector{<:Integer}}=nothing, tol::T=T(1e-6)) where {T<:Real} = nfft_linop(X, kspace_Cartesian_sampling(X; phase_encoding=phase_encoding, subsampling=subsampling); tol=tol)

k_coord(F::NFFTLinOp) = coord(F.K)

AbstractLinearOperators.domain_size(F::NFFTLinOp) = F.X.n
AbstractLinearOperators.range_size(F::NFFTLinOp) = size(F.K)

function AbstractLinearOperators.matvecprod(F::NFFTLinOp{T}, u::AbstractArray{Complex{T},3}) where {T<:Real}
    Kh = reshape(k_coord(F).*reshape([F.X.h...], 1,1,3), :,3)
    phase_shift_origin = exp.(im*sum(Kh.*reshape([F.X.o...],1,3); dims=2)[:,1])
    return reshape(phase_shift_origin.*nufft3d2(Kh[:,1], Kh[:,2], Kh[:,3], -1, F.tol, u), range_size(F))/T(sqrt(prod(domain_size(F))))
end

function AbstractLinearOperators.matvecprod_adj(F::NFFTLinOp{T}, d::AbstractArray{Complex{T},2}) where {T<:Real}
    Kh = reshape(k_coord(F).*reshape([F.X.h...], 1,1,3), :,3)
    phase_shift_origin = exp.(-im*sum(Kh.*reshape([F.X.o...],1,3); dims=2)[:,1])
    return nufft3d1(Kh[:,1], Kh[:,2], Kh[:,3], phase_shift_origin.*vec(d), 1, F.tol, domain_size(F)...)[:,:,:,1]/T(sqrt(prod(domain_size(F))))
end

downscale(F::NFFTLinOp{T}; fact::Integer=1) where {T<:Real} = NFFTLinOp{T}(downscale(F.X; fact=fact), downscale(F.K; fact=fact), F.tol)

## Parameteric linear operator

struct NFFTParametericLinOp{T}
    X::RegularCartesianSpatialSampling{T}
    K::AbstractKSpaceSampling{T}
    tol::T
end

nfft(X::RegularCartesianSpatialSampling{T}, K::AbstractKSpaceSampling{T}; tol::T=T(1e-6)) where {T<:Real} = NFFTParametericLinOp{T}(X, K, tol)

nfft(X::RegularCartesianSpatialSampling{T}; phase_encoding::NTuple{2,Integer}=(1,2), subsampling::Union{Nothing,AbstractVector{<:Integer}}=nothing, tol::T=T(1e-6)) where {T<:Real} = nfft(X, kspace_Cartesian_sampling(X; phase_encoding=phase_encoding, subsampling=subsampling); tol=tol)

function (F::NFFTParametericLinOp{T})(θ::AbstractArray{T,2}) where {T<:Real}
    (size(θ,1) !== size(F.K)[1]) && error("Incompatible time dimension")
    Pτ = phase_shift_linop(F.K, θ[:,1:3])
    Rφ = rotation_linop(θ[:,4:end])
    Kφ = kspace_sampling(Rφ*coord(F.K))
    F = nfft_linop(F.X, Kφ; tol=F.tol)
    return Pτ*F
end

(F::NFFTParametericLinOp)() = F

downscale(F::NFFTParametericLinOp{T}; fact::Integer=1) where {T<:Real} = NFFTParametericLinOp{T}(downscale(F.X; fact=fact), downscale(F.K; fact=fact), F.tol)

## Parameteric delayed evaluation

struct NFFTParametericDelayedEval{T}
    F::NFFTParametericLinOp{T}
    u::AbstractArray{Complex{T},3}
end

Base.:*(F::NFFTParametericLinOp{T}, u::AbstractArray{Complex{T},3}) where {T<:Real} =  NFFTParametericDelayedEval{T}(F, u)

(Fu::NFFTParametericDelayedEval{T})(θ::AbstractArray{T,2}) where {T<:Real} = Fu.F(θ)*Fu.u

## Jacobian of nfft evaluated

struct JacobianNFFTEvaluated{T}<:AbstractLinearOperator{Complex{T},2,2}
    J::AbstractArray{Complex{T},3}
end

function Jacobian(Fu::NFFTParametericDelayedEval{T}, θ::AbstractArray{T,2}) where {T<:Real}

    # Simplifying notation
    u = Fu.u
    X = Fu.F.X
    K = Fu.F.K
    tol = Fu.F.tol
    τ = θ[:,1:3]
    φ = θ[:,4:end]
    P = phase_shift(K) # phase-shift
    R = rotation() # rotation

    # NFFT of u and derivatives thereof
    Kφ, ∂Kφ = ∂(R()*K, φ)
    Fφ = nfft_linop(X, kspace_sampling(Kφ); tol=tol)
    x, y, z = coord(X)
    Fu = Fφ*u
    ∇Fu = -im*cat(Fφ*(u.*x), Fφ*(u.*y), Fφ*(u.*z); dims=3)

    # Computing rigid-body motion Jacobian
    d, Pτ, ∂Pτu = ∂(P()*Fu, τ)
    J = cat(∂Pτu.J, Pτ*dot(∇Fu, ∂Kφ); dims=3)

    return d, Pτ*Fφ, JacobianNFFTEvaluated{T}(J)

end

∂(Fu::NFFTParametericDelayedEval{T}, θ::AbstractArray{T,2}) where {T<:Real} = Jacobian(Fu, θ)

AbstractLinearOperators.domain_size(∂Fu::JacobianNFFTEvaluated) = (size(∂Fu.J,1),6)
AbstractLinearOperators.range_size(∂Fu::JacobianNFFTEvaluated) = size(∂Fu.J)[1:2]
AbstractLinearOperators.matvecprod(∂Fu::JacobianNFFTEvaluated{T}, Δθ::AbstractArray{Complex{T},2}) where {T<:Real} = sum(∂Fu.J.*reshape(Δθ, :, 1, 6); dims=3)[:,:,1]
Base.:*(∂Fu::JacobianNFFTEvaluated{T}, Δθ::AbstractArray{T,2}) where {T<:Real} = ∂Fu*complex(Δθ)
AbstractLinearOperators.matvecprod_adj(∂Fu::JacobianNFFTEvaluated{T}, Δd::AbstractArray{Complex{T},2}) where {T<:Real} = sum(conj(∂Fu.J).*reshape(Δd, size(Δd)..., 1); dims=2)[:,1,:]

## Other utilities

function sparse_matrix_GaussNewton(J::JacobianNFFTEvaluated{T}; W::Union{Nothing,AbstractLinearOperator}=nothing, H::Union{Nothing,AbstractLinearOperator}=nothing) where {T<:Real}
    J = J.J
    GN = similar(J, T, size(J, 1), 6, 6)
    if ~isnothing(W)
        WJ = similar(J)
        @inbounds for i = 1:6
            WJ[:, :, i] = W*J[:,:,i]
        end
    else
        WJ = J
    end
    if ~isnothing(H)
        HWJ = similar(J)
        @inbounds for i = 1:6
            HWJ[:, :, i] = H*WJ[:,:,i]
        end
    else
        HWJ = J
    end
    @inbounds for i = 1:6, j = 1:6
        GN[:,i,j] = vec(real(sum(conj(WJ[:,:,i]).*HWJ[:,:,j]; dims=2)))
    end
    h(i,j) = spdiagm(0 => GN[:,i,j])
    return [h(1,1) h(1,2) h(1,3) h(1,4) h(1,5) h(1,6);
            h(2,1) h(2,2) h(2,3) h(2,4) h(2,5) h(2,6);
            h(3,1) h(3,2) h(3,3) h(3,4) h(3,5) h(3,6);
            h(4,1) h(4,2) h(4,3) h(4,4) h(4,5) h(4,6);
            h(5,1) h(5,2) h(5,3) h(5,4) h(5,5) h(5,6);
            h(6,1) h(6,2) h(6,3) h(6,4) h(6,5) h(6,6)]
end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
#=
    Network generation.
=#

"""
    Returns an Erdos-Renyi random graph according to the G(N, p) model.

    Args:
        N:
            Total number of nodes of the graph.
        p:
            Probability of existence between pair of nodes.
        is_directed:
            Boolean to signal whether the final graph is undirected or
            directed.
        self:
            {default=false} Boolean to signal whether self-loops are allowed.
    Return:
        graph:
            'Graph' structure representing the final graph according to G(N,p).
"""
function erdos_renyi_net(N::Int, p::Float64, is_directed::Bool, self=false)
    if p<=0 || p>=1
        return
    end

    graph = Graph(is_directed)
    for i in 1:N
        add_node(graph)
    end

    w = -1
    lp = log(1.0-p)

    if is_directed # ----- directed -----
        v = 0
        while v < N
            mu = rand(1)[1]
            lr =  log(1.0 - mu)
            w += (1 +  floor(Int, lr/lp))
            while v<N && N<=w
                w -= N
                v += 1
            end
            if v < N
                if v==w && self
                    add_edge(v+1, w+1, graph)
                elseif v!=w
                    add_edge(v+1, w+1, graph)
                end
            end
        end
    else # ----- undirected -----
        v = 1
        while v < N
            mu = rand(1)[1]
            lr = log(1.0 - mu)
            w += 1 + floor(Int, lr/lp)
            while w >= v && v < N
                w -= v
                v += 1
            end
            if v < N
                add_edge(v+1, w+1, graph)
            end
        end
    end
    return graph
end

"""
    Generate a random graph by parsing G(N, \$\\langle k \\rangle\$) where
    \$\\langle k \\rangle\$ is the desired mean degree of the final graph.

    Args:
        N:
            Total number of nodes of the graph.
        k_mean:
            Expected mean degree of the generated graph.
        is_directed:
            Boolean to signal whether the final graph is undirected or
            directed.
        self:
            {default=false} Boolean to signal whether self-loops are allowed.
    Return:
        graph:
            'Graph' structure representing the final ER graph with the expected
            mean degree.
"""
function ER_k(N::Int, k_mean::Float64, is_directed::Bool, self=false)
    p = k_mean/(N-1)
    graph = erdos_renyi_net(N, p, is_directed, self)
    return graph
end

"""
    Generate a random graph by parsing the total number of nodes, expected
    mean degree and number of edge types.

    Args:
        N:
            Total number of nodes of the graph.
        k_mean:
            Expected mean degree of the generated graph.
        n_types:
            Number of edge types for the final graph.
        is_directed:
            Boolean to signal whether the final graph is undirected or
            directed.
        self:
            {default=false} Boolean to signal whether self-loops are allowed.
    Return:
        graph:
            'Graph' structure representing the final ER graph with the expected
            mean degree and number of edge types. This final graph contains an
            edge property with name "edgetype".
"""
function ER_multi(N::Int, k_mean::Float64, n_types::Int, is_directed::Bool, self=false)
    graph = ER_k(N, k_mean, is_directed, self)
    possible_types = [j for j in 1:n_types]

    edgetypes_prop = [ rand(possible_types) for j in 1:length(graph.edges) ]
    set_edges_properties("edgetype", edgetypes_prop, graph)
    return graph
end

"""
    Specific for exporting the generated Erdos-Renyi random graphs.
    The only metadata exported together with the edgelist is the integer
    "edgetype" property

    Args:
        graph:
            'Graph' structure containing the "edgetype" integer property.
        fout:
            filename to where the edgelist should be exported.
"""
function export_edgefile_csv(graph::Graph, fout::String)
    edges = graph.edges
    edgetype = graph.int_eproperties["edgetype"]
    source_arr, target_arr = Int[], Int[]
    type_arr = Int[]
    for edge in edges
        push!(source_arr, edge.source)
        push!(target_arr, edge.target)
        push!(type_arr, edgetype[edge.index])
    end
    CSV.write(fout, DataFrame([source_arr, target_arr, type_arr], [:Source, :Target, :Type]))
end

# To delete
"""
    Function specific for exporting the generated random networks.
    The only metadata exported together is the "edgetype" property.
"""
function export_edgefile(graph::Graph, fout::String)
    edges = graph.edges
    edgetype = graph.int_eproperties["edgetype"]
    open("$fout.txt", "w") do f
        for edge in edges
            write(f, "$(edge.source)\t$(edge.target)\t$(edgetype[edge.index])\n")
        end
    end
end


text:
struct SimpleCurve
    start::Tuple{Number, Number}
    stop::Tuple{Number, Number}
    control::Tuple{Number, Number}

    SimpleCurve(start::Tuple{Number, Number}, stop::Tuple{Number, Number}, control::Tuple{Number, Number}=(start .+ stop) ./ 2) = new(start, stop, control)
end

function getPoint(curve::SimpleCurve, percent::Number)
    return (1 - percent)^2 .* curve.start .+ 2 .* (1 - percent) * percent .* curve.control .+ percent^2 .* curve.stop
end


text:
using GeometricIntegrators.Equations

ϑ₁(t, q, v) = g₁₁(t,q) * v[1] + A₁(t, q)
ϑ₂(t, q, v) = g₂₂(t,q) * v[2] + A₂(t, q)
ϑ₃(t, q, v) = g₃₃(t,q) * v[3] + A₃(t, q)

function ϑ(t, q, v, θ)
    θ[1] = ϑ₁(t,q,v)
    θ[2] = ϑ₂(t,q,v)
    θ[3] = ϑ₃(t,q,v)
    nothing
end

v¹(t, q, p) = g¹¹(t, q) * (p[1] - A₁(t, q))
v²(t, q, p) = g²²(t, q) * (p[2] - A₂(t, q))
v³(t, q, p) = g³³(t, q) * (p[3] - A₃(t, q))

hamiltonian(t,q,p) = (g₁₁(t,q) * v¹(t,q,p)^2 + g₂₂(t,q) * v²(t,q,p)^2 + g₃₃(t,q) * v³(t,q,p)^2) / 2 + φ(t,q)
toroidal_momentum(t,q,p) = p[3]

function charged_particle_3d_pᵢ(qᵢ, vᵢ, tᵢ=0)
    pᵢ = zero(qᵢ)

    if ndims(qᵢ) == 1
        ϑ(tᵢ, qᵢ, vᵢ, pᵢ)
    else
        for i in 1:size(qᵢ,2)
            @views ϑ(tᵢ, qᵢ[:,i], vᵢ[:,i], pᵢ[:,i])
        end
    end
    pᵢ
end

function charged_particle_3d_pode_v(t, q, p, v)
    v[1] = v¹(t,q,p)
    v[2] = v²(t,q,p)
    v[3] = v³(t,q,p)
    nothing
end

function charged_particle_3d_pode_f(t, q, p, f)
    f[1] = dA₁dx₁(t,q) * v¹(t,q,p) + dA₂dx₁(t,q) * v²(t,q,p) + dA₃dx₁(t,q) * v³(t,q,p) + E₁(t,q) +
           (dg₁₁dx₁(t,q) * v¹(t,q,p)^2 + dg₂₂dx₁(t,q) * v²(t,q,p)^2 + dg₃₃dx₁(t,q) * v³(t,q,p)^2) / 2
    f[2] = dA₁dx₂(t,q) * v¹(t,q,p) + dA₂dx₂(t,q) * v²(t,q,p) + dA₃dx₂(t,q) * v³(t,q,p) + E₂(t,q) +
           (dg₁₁dx₂(t,q) * v¹(t,q,p)^2 + dg₂₂dx₂(t,q) * v²(t,q,p)^2 + dg₃₃dx₂(t,q) * v³(t,q,p)^2) / 2
    f[3] = dA₁dx₃(t,q) * v¹(t,q,p) + dA₂dx₃(t,q) * v²(t,q,p) + dA₃dx₃(t,q) * v³(t,q,p) + E₃(t,q) +
           (dg₁₁dx₃(t,q) * v¹(t,q,p)^2 + dg₂₂dx₃(t,q) * v²(t,q,p)^2 + dg₃₃dx₃(t,q) * v³(t,q,p)^2) / 2
    nothing
end

charged_particle_3d_iode_ϑ(t, q, v, θ) = ϑ(t, q, v, θ)

function charged_particle_3d_iode_f(t, q, v, f)
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] + E₁(t,q) +
           (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] + E₂(t,q) +
           (dg₁₁dx₂(t,q) * v[1]^2 + dg₂₂dx₂(t,q) * v[2]^2 + dg₃₃dx₂(t,q) * v[3]^2) / 2
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] + E₃(t,q) +
           (dg₁₁dx₃(t,q) * v[1]^2 + dg₂₂dx₃(t,q) * v[2]^2 + dg₃₃dx₃(t,q) * v[3]^2) / 2
    nothing
end

function charged_particle_3d_iode_g(t, q, v, f)
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] +
           (dg₁₁dx₁(t,q) * v[1]^2 + dg₂₂dx₁(t,q) * v[2]^2 + dg₃₃dx₁(t,q) * v[3]^2) / 2
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] +
           (dg₁₁dx₂(t,q) * v[1]^2 + dg₂₂dx₂(t,q) * v[2]^2 + dg₃₃dx₂(t,q) * v[3]^2) / 2
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] +
           (dg₁₁dx₃(t,q) * v[1]^2 + dg₂₂dx₃(t,q) * v[2]^2 + dg₃₃dx₃(t,q) * v[3]^2) / 2
    nothing
end

# function charged_particle_3d_pode(q₀=qᵢ, v₀=vᵢ)
#     PODE(charged_particle_3d_pode_v, charged_particle_3d_pode_f, q₀, charged_particle_3d_pᵢ(q₀, v₀))
# end
function charged_particle_3d_pode(q₀=qᵢ, p₀=pᵢ)
    PODE(charged_particle_3d_pode_v, charged_particle_3d_pode_f, q₀, p₀)
end

function charged_particle_3d_iode(q₀=qᵢ, p₀=pᵢ)
    IODE(charged_particle_3d_iode_ϑ, charged_particle_3d_iode_f, charged_particle_3d_iode_g,
            q₀, p₀;
            h=hamiltonian, v = (t, q, v) -> nothing)
end


text:
# Uses the Newton method to find the analytic center of a polytope.

function polytope_analcenter(A::AbstractMatrix,
                             b::AbstractVector,
                             x0::AbstractVector,
                             ϵ::Real;
                             α::Real=0.01, β::Real=0.05)
  all(A*x0.<b) || throw(
    DomainError("x0 must lie strictly inside the polytope"))
  x = copy(x0)
  H = zeros(size(x0, 1), size(x0, 1))
  polytope_log(x) = -sum(log.(b-A*x))
  while true
    d = 1./(b-A*x)
    ∇ = A'*d
    H = A'*Diagonal(d)^2*A
    Δxnt = -H\∇
    λ² = -∇⋅Δxnt
    λ²/2 < ϵ && break
    t = 1
    xtest = x + t*Δxnt
    while(any(A*xtest .≥ b)
          || polytope_log(xtest) > polytope_log(x)-t*α*λ²)
      t *= β
      xtest = x + t*Δxnt
    end
    # numeric limit
    polytope_log(xtest) == polytope_log(x) && break
    x = xtest
  end
  x, H
end


text:
using Pkg
Pkg.activate(".")
using DataFrames,Distances,  BenchmarkTools
x = rand(1:4,100)
y = rand(100,5) .+ (x .* (x .==3))

preds = [["a","b","c","d"][i] for i in x]
preds2 = [rand(["a","b","c","d"]) for i in x]
preds3 = [rand(["a","b","c","d"]) for i in x]
df = DataFrame([preds,preds2,preds3],[:X,:Z,:W])

using SnoopCompile,StatsModels,ProfileView
tinf1 =@snoopi_deep begin
    using PERMANOVA
end
ProfileView.view(flamegraph(tinf1))
tinf2 =@snoopi_deep begin
    Hyde = hydra(df,y,BrayCurtis,@formula(1~X+Z+W),10000)
end
ProfileView.view(flamegraph(tinf2))
ttot, pcs = SnoopCompile.parcel(tinf2)
SnoopCompile.write("./precomp/precompiles", pcs)


text:
# an implementation of Uniform Manifold Approximation and Projection
# for Dimension Reduction, L. McInnes, J. Healy, J. Melville, 2018.

# NOTE: unused for now
struct UMAP_{S}
    graph::AbstractMatrix{S}
    embedding::AbstractMatrix{S}

    function UMAP_(graph::AbstractMatrix{S}, embedding::AbstractMatrix{S}) where {S<:AbstractFloat}
        issymmetric(graph) || throw(MethodError("UMAP_ constructor expected graph to be a symmetric matrix"))
        new{S}(graph, embedding)
    end
end

"""
    umap(X::AbstractMatrix[, n_components=2]; <kwargs>) -> embedding

Embed the data `X` into a `n_components`-dimensional space. `n_neighbors` controls
how many neighbors to consider as locally connected.

# Keyword Arguments
- `n_neighbors::Integer = 15`: the number of neighbors to consider as locally connected. Larger values capture more global structure in the data, while small values capture more local structure.
- `metric::SemiMetric = Euclidean()`: the metric to calculate distance in the input space
- `n_epochs::Integer = 300`: the number of training epochs for embedding optimization
- `learning_rate::AbstractFloat = 1.`: the initial learning rate during optimization
- `init::Symbol = :spectral`: how to initialize the output embedding; valid options are `:spectral` and `:random`
- `min_dist::AbstractFloat = 0.1`: the minimum spacing of points in the output embedding
- `spread::AbstractFloat = 1.0`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::AbstractFloat = 1.0`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::AbstractFloat = 1.0`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a::AbstractFloat = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b::AbstractFloat = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
"""
function umap(args...; kwargs...)
    # this is just a convenience function for now
    return UMAP_(args...; kwargs...).embedding
end

function UMAP_(X::AbstractMatrix{S},
               n_components::Integer = 2;
               n_neighbors::Integer = 15,
               metric::SemiMetric = Euclidean(),
               n_epochs::Integer = 300,
               learning_rate::AbstractFloat = 1.,
               init::Symbol = :spectral,
               min_dist::AbstractFloat = 0.1,
               spread::AbstractFloat = 1.0,
               set_operation_ratio::AbstractFloat = 1.0,
               local_connectivity::Integer = 1,
               repulsion_strength::AbstractFloat = 1.0,
               neg_sample_rate::Integer = 5,
               a::Union{AbstractFloat, Nothing} = nothing,
               b::Union{AbstractFloat, Nothing} = nothing
               ) where {S <: AbstractFloat, V <: AbstractVector}
    # argument checking
    size(X, 2) > n_neighbors > 0|| throw(ArgumentError("size(X, 2) must be greater than n_neighbors and n_neighbors must be greater than 0"))
    size(X, 1) > n_components > 1 || throw(ArgumentError("size(X, 1) must be greater than n_components and n_components must be greater than 1"))
    min_dist > 0. || throw(ArgumentError("min_dist must be greater than 0"))
    #n_epochs > 0 || throw(ArgumentError("n_epochs must be greater than 1"))

    # main algorithm
    graph = fuzzy_simplicial_set(X, n_neighbors, metric, local_connectivity, set_operation_ratio)

    embedding = initialize_embedding(graph, n_components, Val(init))

    embedding = optimize_embedding(graph, embedding, n_epochs, learning_rate, min_dist, spread, repulsion_strength, neg_sample_rate)
    # TODO: if target variable y is passed, then construct target graph
    #       in the same manner and do a fuzzy simpl set intersection

    return UMAP_(graph, embedding)
end

"""
    fuzzy_simplicial_set(X, n_neighbors) -> graph::SparseMatrixCSC

Construct the local fuzzy simplicial sets of each point in `X` by
finding the approximate nearest `n_neighbors`, normalizing the distances
on the manifolds, and converting the metric space to a simplicial set.
"""
function fuzzy_simplicial_set(X::AbstractMatrix,
                              n_neighbors,
                              metric,
                              local_connectivity,
                              set_operation_ratio) where {V <: AbstractVector}
    if size(X, 2) < 4096
        # compute all pairwise distances
        knns, dists = pairwise_knn(X, n_neighbors, metric)
    else
        knngraph = DescentGraph(X, n_neighbors, metric)
        knngraph.graph::Matrix{Tuple{S,T}} where {S<:Integer,T<:AbstractFloat}
        knns = Array{typeof(knngraph.graph[1][1])}(undef, size(knngraph.graph))
        dists = Array{typeof(knngraph.graph[1][2])}(undef, size(knngraph.graph))
        for index in eachindex(knngraph.graph)
            @inbounds knns[index] = knngraph.graph[index][1]
            @inbounds dists[index] = knngraph.graph[index][2]
        end
    end

    σs, ρs = smooth_knn_dists(dists, n_neighbors, local_connectivity)

    rows, cols, vals = compute_membership_strengths(knns, dists, σs, ρs)
    fs_set = sparse(rows, cols, vals, size(knns, 2), size(knns, 2))
                                      # sparse matrix M[i, j] = vᵢⱼ where
                                      # vᵢⱼ is the probability that j is in the
                                      # simplicial set of i
    return dropzeros(fs_set .+ fs_set' .- fs_set .* fs_set')
end

"""
    smooth_knn_dists(dists, k; <kwargs>) -> knn_dists, nn_dists

Compute the distances to the nearest neighbors for a continuous value `k`. Returns
the approximated distances to the kth nearest neighbor (`knn_dists`)
and the nearest neighbor (nn_dists) from each point.

# Keyword Arguments
...
"""
function smooth_knn_dists(knn_dists::AbstractMatrix{S}, k::Integer, local_connectivity::Integer;
                          niter::Integer=64,
                          bandwidth::AbstractFloat=1.,
                          ktol = 1e-5) where {S <: Real}
    @inline minimum_nonzero(dists) = minimum(dists[dists .> 0.])
    ρs = S[minimum_nonzero(knn_dists[:, i]) for i in 1:size(knn_dists, 2)]
    σs = Array{S}(undef, size(knn_dists, 2))

    for i in 1:size(knn_dists, 2)
        @inbounds σs[i] = smooth_knn_dist(knn_dists[:, i], ρs[i], k, local_connectivity, bandwidth, niter, ktol)
    end
    return ρs, σs
end

@fastmath function smooth_knn_dist(dists::AbstractVector, ρ, k, local_connectivity, bandwidth, niter, ktol)
    target = log2(k)*bandwidth
    lo, mid, hi = 0., 1., Inf
    #psum(dists, ρ) = sum(exp.(-max.(dists .- ρ, 0.)/mid))
    for n in 1:niter
        psum = sum(exp.(-max.(dists .- ρ, 0.)./mid))
        if abs(psum - target) < ktol
            break
        end
        if psum > target
            hi = mid
            mid = (lo + hi)/2.
        else
            lo = mid
            if hi == Inf
                mid *= 2.
            else
                mid = (lo + hi) / 2.
            end
        end
    end
    # TODO: set according to min k dist scale
    return mid
end

"""
    compute_membership_strengths(knns, dists, σ, ρ) -> rows, cols, vals

Compute the membership strengths for the 1-skeleton of each fuzzy simplicial set.
"""
function compute_membership_strengths(knns::AbstractMatrix{S},
                                      dists::AbstractMatrix{T},
                                      ρs::Vector{T},
                                      σs::Vector{T}) where {S <: Integer, T}
    # set dists[i, j]
    rows = sizehint!(S[], length(knns))
    cols = sizehint!(S[], length(knns))
    vals = sizehint!(T[], length(knns))
    for i in 1:size(knns, 2), j in 1:size(knns, 1)
        @inbounds if i == knns[j, i] # dist to self
            d = 0.
        else
            @inbounds d = exp(-max(dists[j, i] - ρs[i], 0.)/σs[i])
        end
        append!(cols, i)
        append!(rows, knns[j, i])
        append!(vals, d)
    end
    return rows, cols, vals
end

function initialize_embedding(graph, n_components, ::Val{:spectral})
    embed = spectral_layout(graph, n_components)
    # expand
    expansion = 10. / maximum(embed)
    @. embed = (embed*expansion) + randn(size(embed))*0.0001
    return embed
end

function initialize_embedding(graph, n_components, ::Val{:random})
    return 20. .* rand(n_components, size(graph, 1)) .- 10.
end

"""
    optimize_embedding(graph, embedding, n_epochs, initial_alpha, min_dist, spread) -> embedding

Optimize an embedding by minimizing the fuzzy set cross entropy between the high and low
dimensional simplicial sets using stochastic gradient descent.

# Arguments
- `graph`: a sparse matrix of shape (n_samples, n_samples)
- `embedding`: a dense matrix of shape (n_components, n_samples)
- `neg_sample_rate::Integer`: the number of negative samples per positive sample
"""
function optimize_embedding(graph,
                            embedding,
                            n_epochs,
                            initial_alpha,
                            min_dist,
                            spread,
                            gamma,
                            neg_sample_rate)
    a, b = fit_ϕ(min_dist, spread)

    alpha = initial_alpha
    for e in 1:n_epochs

        @fastmath @inbounds for i in 1:size(graph, 2)
            for ind in nzrange(graph, i)
                j = rowvals(graph)[ind]
                p = nonzeros(graph)[ind]
                if rand() <= p
                    @views sdist = evaluate(SqEuclidean(), embedding[:, i], embedding[:, j])
                    if sdist > 0.
                        delta = (-2. * a * b * sdist^(b-1))/(1. + a*sdist^b)
                    else
                        delta = 0.
                    end
                    @simd for d in 1:size(embedding, 1)
                        grad = clamp(delta * (embedding[d,i] - embedding[d,j]), -4., 4.)
                        embedding[d,i] += alpha * grad
                        embedding[d,j] -= alpha * grad
                    end

                    for _ in 1:neg_sample_rate
                        k = rand(1:size(graph, 2))
                        @views sdist = evaluate(SqEuclidean(),
                                                embedding[:, i], embedding[:, k])
                        if sdist > 0
                            delta = (2. * gamma * b) / ((0.001 + sdist)*(1. + a*sdist^b))
                        elseif i == k
                            continue
                        else
                            delta = 0.
                        end
                        @simd for d in 1:size(embedding, 1)
                            if delta > 0.
                                grad = clamp(delta * (embedding[d, i] - embedding[d, k]), -4., 4.)
                            else
                                grad = 4.
                            end
                            embedding[d, i] += alpha * grad
                        end
                    end

                end
            end
        end
        alpha = initial_alpha*(1. - e/n_epochs)
    end

    return embedding
end

"""
    fit_ϕ(min_dist, spread) -> a, b

Find a smooth approximation to the membership function of points embedded in ℜᵈ.
This fits a smooth curve that approximates an exponential decay offset by `min_dist`.
"""
function fit_ϕ(min_dist, spread)
    ψ(d) = d >= min_dist ? exp(-(d - min_dist)/spread) : 1.
    xs = LinRange(0., spread*3, 300)
    ys = map(ψ, xs)
    @. curve(x, p) = (1. + p[1]*x^(2*p[2]))^(-1)
    result = curve_fit(curve, xs, ys, [1., 1.], lower=[0., -Inf])
    a, b = result.param
    return a, b
end

"""
    spectral_layout(graph, embed_dim) -> embedding

Initialize the graph layout with spectral embedding.
"""
function spectral_layout(graph::SparseMatrixCSC{T},
                         embed_dim::Integer) where {T<:AbstractFloat}
    D_ = Diagonal(dropdims(sum(graph; dims=2); dims=2))
    D = inv(sqrt(D_))
    # normalized laplacian
    # TODO: remove sparse() when PR #30018 is merged
    L = sparse(Symmetric(I - D*graph*D))

    k = embed_dim+1
    num_lanczos_vectors = max(2k+1, round(Int, sqrt(size(L, 1))))
    local layout
    try
        # get the 2nd - embed_dim+1th smallest eigenvectors
        eigenvals, eigenvecs = eigs(L; nev=k,
                                       ncv=num_lanczos_vectors,
                                       which=:SM,
                                       tol=1e-4,
                                       v0=ones(T, size(L, 1)),
                                       maxiter=size(L, 1)*5)
        layout = permutedims(eigenvecs[:, 2:k])::Array{T, 2}
    catch e
        print("\n", e, "\n")
        print("Error occured in spectral_layout;
               falling back to random layout.\n")
        layout = 20 .* rand(T, embed_dim, size(L, 1)) .- 10
    end
    return layout
end


text:
module UPPE

export CPU, GPU, adapt_array,
       Grid, GridT, GridRT,
       Field, field2intensity, intensity2field,
       refractive_index, k_func, chi3_func, group_velocity, diffraction_length,
       dispersion_length,
       Response, Model,
       Simulation, run!,
       RK2, RK3, RK4, Tsit5, ATsit5   # ODEIntegrators algorithms

using AnalyticSignals: rsig2asig!, rsig2aspec!
using CUDA: CuVector, CuMatrix, CuArray, @cuda, launch_configuration,
            threadIdx, blockIdx, blockDim, gridDim
using FFTW: fftfreq, ifftshift, plan_fft!
using HankelTransforms: dhtcoord, dhtfreq, plan_dht
using HDF5
using ODEIntegrators
using Printf: @printf, @sprintf
using TimerOutputs: @timeit, reset_timer!, get_defaulttimer

using PhysicalConstants.CODATA2018
const C0 = SpeedOfLightInVacuum.val
const EPS0 = VacuumElectricPermittivity.val
const MU0 = VacuumMagneticPermeability.val
const QE = ElementaryCharge.val
const ME = ElectronMass.val
const HBAR = ReducedPlanckConstant.val

abstract type ARCH{T} end
struct CPU{T} <: ARCH{T} end
struct GPU{T} <: ARCH{T} end
CPU() = CPU{Float64}()
GPU() = GPU{Float32}()

float_type(::ARCH{T}) where T = T

adapt_array(::CPU, a::AbstractArray) = a
adapt_array(::GPU, a::AbstractArray) = CuArray(a)

include("grids.jl")
include("fields.jl")
include("medium.jl")
include("guards.jl")
include("models.jl")
include("field_analyzers.jl")
include("outputs.jl")
include("simulations.jl")

end


text:
const PRIMARY_COLOR = :blue
const SECONDARY_COLOR = :red

# this function is required to plot integration areas
function get_left_right_points(
    xs::AbstractArray{T},
    ys::AbstractArray{T},
    xₗ::T,
    xᵣ::T,
    b::Union{Nothing, UncertainBound}
    ;
    baseline_handling=nothing
) where {T<:AbstractFloat}

    xₗ, xᵣ = xₗ < xᵣ ? (xₗ, xᵣ) : (xᵣ, xₗ)

    i = searchsortedfirst(xs, xₗ)
    j = searchsortedfirst(xs, xᵣ)
    any(
        [
            i < 2,
            i > length(xs) - 1,
            j < 3,
            j > length(xs)
        ]
    ) && throw(error("At least one integration bound is outside the support region ($(minimum(x)), $(maximum(x)))."))

    if baseline_handling == "local"
        _, _, yₗ, yᵣ = _local_baseline(xs, ys, xₗ, xᵣ, b)
    else
        # bound `b` is unused in this case
        _, _, yₗ, yᵣ = _endpoint_to_endpoint_baseline(xs, ys, xₗ, xᵣ)
    end

    return i, j, xₗ, xᵣ, yₗ, yᵣ
end

@recipe function plot_recipe(crv::AbstractCurve)
    return crv.x, crv.y
end

# --------------------------------------
# enable plotting of sample draws
# --------------------------------------

@recipe function plot_recipe(c::AbstractCurve, uc::UncertainCurve, draws=3)
    draws < 0 && throw(ArgumentError("Number of samples must be > 0."))

    legend := :none
    layout := (draws + 1, 1)
    link := :both
    size --> (500, 600)

    delete!(plotattributes, :draws)

    for i ∈ 0:draws
        @series begin
            subplot := i + 1
            if i == 0
                seriescolor := SECONDARY_COLOR
                yguide := "input"
                c
            else
                yguide := "sample $(i)"
                get_draw(i, uc)
            end
        end
    end
end

@recipe function plot_recipe(crv::Curve{T}, left::T, right::T; subtract_baseline=false, local_baseline=false, bound=nothing) where T
    fillrange := 0
    fillalpha --> 0.5
    fillcolor --> :orange
    linewidth --> 0
    label     --> nothing

    (local_baseline && bound == nothing) && error("You have to provide a bound if local_baseline == true.") |> throw
    (subtract_baseline && local_baseline) && error("local_baseline and subtract_baseline cannot both be true.") |> throw

    left = T(left)
    right = T(right)

    if subtract_baseline
        bh = "end-to-end"
    elseif local_baseline
        bh = "local"
    else
        bh = nothing
    end

    l, r, xl, xr, yl, yr = get_left_right_points(crv.x, crv.y, left, right, bound; baseline_handling=bh)
    if !(local_baseline || subtract_baseline)
        x = [xl; crv.x[l+1:r-1]; xr; xr; xl]
        y = [yl; crv.y[l+1:r-1]; yr; zero(T); zero(T)]
    elseif subtract_baseline
        x = [xl; crv.x[l+1:r-1]; xr; xl]
        y = [yl; crv.y[l+1:r-1]; yr; yl]
    else
        _, _, ycl, ycr = _endpoint_to_endpoint_baseline(crv.x, crv.y, left, right) # to get extra points located on the curve
        x = [xl; crv.x[l+1:r-1]; xr; xr; xl]
        y = [ycl; crv.y[l+1:r-1]; ycr; yr; yl]

    end

    return x, y
end

@recipe function plot_recipe(crv::Curve{T}, bnds::Vector{UncertainBound{T, N}}, draw::Int) where {T, N}
    @series begin
        crv
    end
    for (j, b) ∈ enumerate(bnds)
        @series begin
            fillcolor := j % 2 == 1 ? :red : :orange
            left, right = get_draw(draw, b)
            bound := b
            crv, left, right
        end
    end
end

@recipe function plot_recipe(crv::Curve{T}, bnd::UncertainBound{T, N}, draw::Int) where {T, N}
    @series begin
        crv, [bnd], draw
    end
end

# plot draws of curves alongside with draws of bounds
@recipe function plot_recipe(
    uc::UncertainCurve{T, N},
    bnds::Vector{UncertainBound{T, N}}
    ;
    draws=3,
    subtract_baseline=true
) where {T, N}

    legend := :none
    layout := (draws + 1, 1)
    link := :both
    size --> (500, 600)

    mean_uc = mean(uc)

    for i ∈ 0:draws
        for (j, b) in enumerate(bnds)
            @series begin
                fillcolor := j % 2 == 1 ? :red : :orange
                subplot := i + 1
                bound := b
                if i == 0
                    mean_uc, mean(b)...
                else
                    get_draw(i, uc), get_draw(i, b)...
                end
            end
        end
        @series begin
            subplot := i + 1
            if i == 0
                # mean spectrum
                seriescolor := SECONDARY_COLOR
                yguide := "mean"
                mean_uc
            else
                seriescolor := PRIMARY_COLOR
                yguide := "sample $(i)"
                get_draw(i, uc)
            end
        end
    end
end

@recipe function plot_recipe(
    uc::UncertainCurve{T, N},
    bnd::UncertainBound{T, N}
    ;
    draws=3
) where {T, N}
    draws := draws
    uc, [bnd]
end

# --------------------------------------
# enable plotting of noise sample draws
# --------------------------------------

@recipe function plot_recipe(x::Vector{T}, nm::AbstractNoiseModel; draws=3, subplot_offset=0) where {T <: Real}
    draws < 0 && throw(ArgumentError("Number of samples must be > 0."))

    layout --> (draws, 1)
    legend --> :none
    link --> :both
    size --> (500, 600)

    delete!(plotattributes, :grid_points)
    delete!(plotattributes, :draws)

    S = generate_noise(nm, length(x), draws)
    for i ∈ 1:draws
        @series begin
            subplot := i + subplot_offset # the offset does only apply if plotting together with a noise sample
            yguide := "sample $(i)"
            x, get_draw.(i, S)
        end
    end
end
@recipe function plot_recipe(nm::AbstractNoiseModel; gridpoints=1000, draws=3)
    @series begin
        x = collect(eltype(nm), 1:gridpoints)
        x, nm
    end
end

@recipe function plot_recipe(ns::NoiseSample, nm::AbstractNoiseModel; draws=3)

    layout --> (draws+1, 1)
    legend --> :none
    link --> :both
    size --> (500, 600)

    @series begin
            subplot := 1
            yguide := "input"
            seriescolor --> SECONDARY_COLOR
            ns
    end
    @series begin
            subplot_offset := 1
            draws := draws
            ns.x, nm
    end
end

# -----------------------------------
# enable plotting of auto covariance
# -----------------------------------

"""
    plotautocovfit(ns::NoiseSample, nm::MvGaussianNoiseModel; kw...)

Plot results of autocovariance fit.
"""
@userplot PlotAutoCovFit
@recipe function plot_repice(pac::PlotAutoCovFit)

    ns, nm = pac.args

    xguide := "lag"
    yguide := "auto-covariance"

    lags, acov = estimate_autocov(ns)

    @series begin
        label := "estimate"
        seriescolor --> PRIMARY_COLOR
        lags, acov
    end

    @series begin
        label := @sprintf "fit (α = %.3e, λ = %.3e)" nm.α nm.λ
        seriescolor --> SECONDARY_COLOR
        lags, gauss_kernel(lags, [nm.α, nm.λ])
    end
end

MonteCarloMeasurements.mcplot(uc::UncertainCurve; draws=10, alpha=0.5, kw...) = MonteCarloMeasurements.mcplot(uc.x, uc.y, draws; alpha=0.5, kw...)

# --------------------------------------------
# enable plotting of UncertainBound histograms
# --------------------------------------------

@recipe plot_repice(::Type{T}, ub::T) where {T <: UncertainBound} = [ub.left.particles, ub.right.particles]


text:
using Transformers
using Transformers.Pretrain
using Transformers.BidirectionalEncoder: WordPiece, bert_cased_tokenizer

using TextEncodeBase
using TextEncodeBase: NestedTokenizer, BaseTokenization, Sentence, Word, SubWord, getvalue, Splittable

struct BertCasedTokenization <: BaseTokenization
    wordpiece::WordPiece
end

# split sentence with `bert_cased_tokenizer` (define with WordTokenizers.jl's `TokenBuffer`)
TextEncodeBase.splitting(::BertCasedTokenization, s::Sentence) = bert_cased_tokenizer(getvalue(s))

# word is splittable with WordPiece
TextEncodeBase.splittability(::BertCasedTokenization, w::Word) = Splittable()

# split word with `WordPiece`
TextEncodeBase.splitting(t::BertCasedTokenization, w::Word) = t.wordpiece(getvalue(w))

tokenizer = pretrain"bert-cased_L-12_H-768_A-12:tokenizer" # this is just `bert_cased_tokenizer`
wordpiece = pretrain"bert-cased_L-12_H-768_A-12:wordpiece"

tkr = NestedTokenizer(BertCasedTokenization(wordpiece))

text1 = "Peter Piper picked a peck of pickled peppers"
single_without_TEB = text1 |> tokenizer |> wordpiece
single_with_TEB = tkr(Sentence(text1))

# `NestedTokenizer` return vector of vector
@assert single_without_TEB == map(getvalue, single_with_TEB[])

#=
julia> single_without_TEB
11-element Vector{String}:
 "Peter"
 "Piper"
 "picked"
 "a"
 "p"
 "##eck"
 "of"
 "pick"
 "##led"
 "pepper"
 "##s"

julia> single_with_TEB
1-element Vector{Vector{TextEncodeBase.TokenStage}}:
 [Token("Peter"), Token("Piper"), Token("picked"), Token("a"), Token("p"), Token("##eck"), Token("of"), Token("pick"), Token("##led"), Token("pepper"), Token("##s")]

julia> single_without_TEB == map(getvalue, single_with_TEB[])
true

=#

# define stage for batch of data
TextEncodeBase.@stage BatchSentence{A<:AbstractVector, M} TextEncodeBase.DocumentStage

# struct BatchSentence{A<:AbstractVector, M} <: TextEncodeBase.DocumentStage
#     x::A
#     meta::M
# end
# BatchSentence(x) = BatchSentence(x, nothing)
# TextEncodeBase.setmeta(x::BatchSentence, meta) = BatchSentence(x.x, meta)
# TextEncodeBase.setvalue(x::BatchSentence, y) = BatchSentence(y, x.meta)

# splittability and split behavior for `BatchSentence`
TextEncodeBase.splittability(::BertCasedTokenization, ::BatchSentence) = Splittable()
TextEncodeBase.splitting(::BertCasedTokenization, s::BatchSentence) = s.x

text2 = "Fuzzy Wuzzy was a bear"
texts = [text1, text2]

batch_without_TEB = map(wordpiece∘tokenizer, texts)
batch_with_TEB = tkr(BatchSentence(texts))

@assert batch_without_TEB == TextEncodeBase.nestedcall(getvalue, batch_with_TEB)

#=
julia> batch_without_TEB
2-element Vector{Vector{String}}:
 ["Peter", "Piper", "picked", "a", "p", "##eck", "of", "pick", "##led", "pepper", "##s"]
 ["Fu", "##zzy", "Wu", "##zzy", "was", "a", "bear"]

julia> batch_with_TEB
2-element Vector{Vector{TextEncodeBase.TokenStage}}:
 [Token("Peter"), Token("Piper"), Token("picked"), Token("a"), Token("p"), Token("##eck"), Token("of"), Token("pick"), Token("##led"), Token("pepper"), Token("##s")]
 [Token("Fu"), Token("##zzy"), Token("Wu"), Token("##zzy"), Token("was"), Token("a"), Token("bear")]

julia> batch_without_TEB == TextEncodeBase.nestedcall(getvalue, batch_with_TEB)
true

=#

using TextEncodeBase: IndexedTokenization

itkr = NestedTokenizer(IndexedTokenization(BertCasedTokenization(wordpiece)))

ibatch_with_TEB = itkr(BatchSentence(texts))

#=
# subword from same word having the same `word_id`
julia> ibatch_with_TEB[1]
11-element Vector{TextEncodeBase.TokenStage}:
 Token("Peter", (sentence_id = 1, word_id = 1, token_id = 1))
 Token("Piper", (sentence_id = 1, word_id = 2, token_id = 2))
 Token("picked", (sentence_id = 1, word_id = 3, token_id = 3))
 Token("a", (sentence_id = 1, word_id = 4, token_id = 4))
 Token("p", (sentence_id = 1, word_id = 5, token_id = 5))
 Token("##eck", (sentence_id = 1, word_id = 5, token_id = 6))
 Token("of", (sentence_id = 1, word_id = 6, token_id = 7))
 Token("pick", (sentence_id = 1, word_id = 7, token_id = 8))
 Token("##led", (sentence_id = 1, word_id = 7, token_id = 9))
 Token("pepper", (sentence_id = 1, word_id = 8, token_id = 10))
 Token("##s", (sentence_id = 1, word_id = 8, token_id = 11))

julia> ibatch_with_TEB[2]
7-element Vector{TextEncodeBase.TokenStage}:
 Token("Fu", (sentence_id = 2, word_id = 1, token_id = 1))
 Token("##zzy", (sentence_id = 2, word_id = 1, token_id = 2))
 Token("Wu", (sentence_id = 2, word_id = 2, token_id = 3))
 Token("##zzy", (sentence_id = 2, word_id = 2, token_id = 4))
 Token("was", (sentence_id = 2, word_id = 3, token_id = 5))
 Token("a", (sentence_id = 2, word_id = 4, token_id = 6))
 Token("bear", (sentence_id = 2, word_id = 5, token_id = 7))

=#

using TextEncodeBase: nestedcall, with_head_tail, trunc_and_pad, nested2batch

# construct `Vocab` with `WordPiece`
vocab = Vocab(wordpiece.vocab, wordpiece.vocab[wordpiece.unk_idx])

# define encoder with `TextEncoder`
enc = TextEncoder(
    itkr, vocab,
    nested2batch ∘ trunc_and_pad(nothing, vocab.unk) ∘ with_head_tail("[CLS]", "[SEP]") ∘ nestedcall(getvalue)
)

#=
julia> encode(enc, BatchSentence(texts))
28996x13x2 OneHotArray{28996, 3, Matrix{OneHot{0x00007144}}}:
[...]

julia> decode(enc, ans)
13×2 Matrix{String}:
 "[CLS]"   "[CLS]"
 "Peter"   "Fu"
 "Piper"   "##zzy"
 "picked"  "Wu"
 "a"       "##zzy"
 "p"       "was"
 "##eck"   "a"
 "of"      "bear"
 "pick"    "[SEP]"
 "##led"   "[UNK]"
 "pepper"  "[UNK]"
 "##s"     "[UNK]"
 "[SEP]"   "[UNK]"

=#