text:
#========================================================================================#
"""
	Decoders

A Decoder of n-bit vectors to floating-point numbers.

Author: Niall Palfreyman, 8/3/2022.
"""
module Decoders

# Externally callable methods of Decoders
export Decoder, decode

#-----------------------------------------------------------------------------------------
# Module types:

"""
	Decoder

A Decoder of bit-vectors whose bit-entries encode ndims dimensions with precision nbits bits.
"""
struct Decoder
	lwb											# Lower bound of each domain dimension
	ndims										# Number of dimensions in domain
	nbits										# Number of bits encoding each dimension
	coeffs										# Matrix of coefficients for decoding bit-vectors

	"Construct a new Decoder"
	function Decoder( domain::Vector=[[0 1]], nbits::Int=10)
		# Set up decoding apparatus:
		span = ((x->x[2]-x[1]).(domain))'		# Row vector of domain dimension extents
		qlevels = 0.5.^(1:nbits)				# Column vector of successive halving fractions
		qlevels = qlevels / sum(qlevels)		# Normalise qlevels to cover entire unit interval

		# Initialise all fields of the decoding apparatus:
		new(
			(x->x[1]).(domain),					# lwb: Extract the lower bound of each dimension
			length(domain),						# ndims
			nbits,								# nbits
			qlevels*span						# coeffs
		)
	end
end

#-----------------------------------------------------------------------------------------
# Module methods:

"""
	decode( decoder, x)

Decode the bit-vector x into floating-point numbers. x is a data vector of ndims
floating-point values, each of which is encoded by nbits bits. Therefore, x contains
ndims*nbits bits that we must decode into a data vector of ndims floating-point values.
"""
function decode( decoder::Decoder, x::Vector)
	# Reshape this column of x into a data matrix of ndims columns, each of which contains
	# nbits rows. Then decode each column of data by multiplying it by components in the
	# corresponding column of coeffs, summing the column and adding the corresponding lwb:
	data = ones( decoder.nbits, decoder.ndims)		# reshape() x to the right size.
	sum(decoder.coeffs .* data,dims=1)'				# Add the sum of each column to lwb, then
end													# transpose it back into a column vector.

#-----------------------------------------------------------------------------------------
"""
	demo()

Provide a demonstration of using a Decoder to decode two bit-vectors into into the two
vectors [1.67742, 3.51613, -2.74194] and [1.32258, 4.41935,  1.51613].
"""
function demo()
	bitdata = [ [1,0,1,0,1,0,1,0,0,0,0,1,0,1,1], [0,1,0,1,0,1,0,1,1,0,1,0,1,1,0] ]
	decoder = Decoder( [[1 2],[3 5],[-7 5]], 5)		# Three dimensions, each of precision 5 bits

	println(
		"This bit-data contains two data vectors, each encoding 3 dimensions " *
		"with precision 5 bits:"
	)
	display( bitdata)
	println()
	println( "Here is the Decoder:")
	display( decoder)
	println()
	println( "Here are the two decoded vectors containing three floats:")
	display(
		map(bitdata) do x
			decode( decoder, x)
		end
	)
end

end		# of Decoders


text:
#=
old_code:
- Julia version: 1.0.3
- Author: rhdekker
- Date: 2019-01-11
=#

struct Container
    name::String
    id::Int64
end

function the_streaming_way_of_doing_things()
    tag_name = "test"
    L = []
    cb = XPCallbacks()
    cb.start_element = function (h, name, attrs)
        c = Container(name, 123)
        push!(L, c)
        return c
    end

    parse(xml, cb)
    print(L)

end

function bla(all_nodes)
    println("Called")
    vg = VariantGraphNode[]
    startNode = StartNode()
    push!(vg, startNode)

    vgedges = []
    last_node = startNode

    for node in all_nodes
        if typeof(node) == String
            println(node)
            # here we want to create a text node for in the variant hypergraph
            textNode = TextNode()
            push!(vg, startNode)
            # add an extra edge
            push!(vgedges, (last_node, textNode))
            last_node = textNode
        elseif typeof(node) == ETree
            println(node.name)
        end
    end

    # create end node
end

    # Ik wil over de nodes lopen en die op een andere manier bij elkaar groeperen zodat
    # het werken met mixed content xml niet meer zo onvoorspelbaar is.
    # Mixed content XML is in te delen in blokken met (tag, attr, textual content en tail textual content)
    # Deze constructie heb ik voor het eerst gezien in de XML parser van Python
    # Deze aanpak heeft voor en nadelen, maar in dit geval kan het zorgen voor een voorspelbare structuur
    # op basis waarvan we verdere bewerkingen kunnen doen.

    # OUDE comments
    # ik moet een methode hebben die gewoon over alle elements itereert.
    # maar het eerste etree element is al meteen de root..
    # Zijn er traversals ofzo?

    # wat we willen doen is alle text elementen uit de xml tree halen
    # en die dan in een variant graph achtige constructie plaatsen
    # Zoals elke variant graph is er een start en end vertex
    # Eigenlijk moeten we hier ook tokenizen
    # dus dan na de start node zijn er twee text nodes
    # dan een divergence node
    # met twee text nodes (well good) aan de ene kant
    # en twee texdt nodes (def bad) ana de andere kant.
    # had de java versie al divergence en convergence nodes?

    # moet ik zelf een graph struct aanmaken?
    # of een liberary daarvoor gebruiken
    # mijn eerste gevoel is om het gewoon even zelf te doen
    # even wat simpele structs te definieren

ext: jl


text:
#=
This script runs step-by-step through md-cRQA analysis in Julia. It uses DynamicalSystems.jl for delay embedding and analysis.
Please see 'example_embedding.jl' for more info
=#

# ----- this script relies on the following packages -----
using Pkg
Pkg.add("DataFrames");
Pkg.add("CSV");
Pkg.add("DelayEmbeddings");
Pkg.add("DynamicalSystems"); # RQA stuffs
Pkg.add("Statistics");
Pkg.add("Plots"); # plotting
Pkg.add("LinearAlgebra") # diag caculation
Pkg.add("Distances") # matrix pairwise distances
Pkg.add("ElasticArrays") # efficient matrices

using Distances, ElasticArrays, LinearAlgebra;
using DataFrames, CSV, DynamicalSystems, DelayEmbeddings, StatsBase, Statistics, Plots;

include("mdFnn.jl");
include("mdEmbed.jl");

# Step 1 ----- load in an example dataset
data = DataFrame(CSV.File("exampleData.csv"));

ts1 = data[1:end,2:3];
ts2 = data[1:end,5:6];

# get the dimensions (number of columns) of the data
ncol_ts1 = size(ts1,2);
ncol_ts2 = size(ts2,2);

if ncol_ts1 > 1
    ts1_names = names(ts1) |> permutedims; # use this for ami plotting
    ts1 = Matrix(ts1);
else
    ts1_names = "timeseries"
end

if ncol_ts2 > 1
    ts2_names = names(ts2) |> permutedims; # use this for ami plotting
    ts2 = Matrix(ts2);

else
    ts2_names = "timeseries"
end

# normalize the timeseries if from different signals or dimensions
# |> save as a Dataset (per DynamicalSystems.jl documentation): for now this does nothing,
# but futureproofing
ts1 = StatsBase.standardize(ZScoreTransform, ts1; dims=1, center=true, scale=true) |> Dataset
ts2 = StatsBase.standardize(ZScoreTransform, ts2; dims=1, center=true, scale=true) |> Dataset

# get the dimensions (number of columns) of the data
ncol_ts1 = size(ts1,2);
ncol_ts2 = size(ts2,2);

# Step 2 ----- find the delay (tau)

# parameters:
nbins = 10;
maxlag = (size(data,1) * .05) |> round |> Int;

# get delay for each dimension of ts1
lags_ts1 = zeros(ncol_ts1);

plotMI_ts1 = zeros(maxlag+1, ncol_ts1);
for dimension in 1:ncol_ts1
    plotMI_ts1[:,dimension] = selfmutualinfo(ts1[:,dimension],0:1:maxlag; nbins=nbins);
    lags_ts1[dimension] = estimate_delay(ts1[:,dimension], "mi_min", 0:1:maxlag; nbins=nbins);
end

# get delay for each dimension of ts2
lags_ts2 = zeros(ncol_ts2);

plotMI_ts2 = zeros(maxlag+1, ncol_ts2);
for dimension in 1:ncol_ts2
    plotMI_ts2[:,dimension] = selfmutualinfo(ts2[:,dimension],0:1:maxlag; nbins=nbins);
    lags_ts2[dimension] = estimate_delay(ts2[:,dimension], "mi_min", 0:1:maxlag; nbins=nbins);
end

τ = maximum([lags_ts1 lags_ts2]) |> Int;

# Step 3 ------ find the embedding dimension

# parameters
tau = τ;
maxEmb = 10;
numSamples = 500;
Rtol = 10;
Atol = 2;
diff_threshold = 3

fnnPerc_ts1 = mdFnn(ts1, tau, maxEmb, numSamples, Rtol, Atol);
fnnPerc_ts2 = mdFnn(ts2, tau, maxEmb, numSamples, Rtol, Atol);

# embedding dimension
D_ts1_below_threshold = findall(fnnPerc_ts1 .< 1)
D_ts2_below_threshold = findall(fnnPerc_ts2 .< 1)

if size(D_ts1_below_threshold)[1]>0
    D_ts1 = (findall(fnnPerc_ts1 .< 1) |> minimum);
else
    D_ts1 = (findall(abs.(diff(fnnPerc_ts1)) .< diff_threshold) |> minimum) + 1
end

if size(D_ts2_below_threshold)[1]>0
    D_ts2 = (findall(fnnPerc_ts2 .< 1) |> minimum);
else
    D_ts2 = (findall(abs.(diff(fnnPerc_ts2)) .< diff_threshold) |> minimum) + 1
end

D = maximum([D_ts1 D_ts2]);

# Step 4 ------ embed timeseries

embed_ts1 = mdEmbed(ts1, D, τ);
embed_ts2 = mdEmbed(ts2, D, τ);

# Step 5 ------ construct Recurrence Matrix

# cross recurrence matrix
CRP = CrossRecurrenceMatrix(embed_ts1, embed_ts2, .05; fixedrate = true);

# you can create a JRP using:
# JRP = JointRecurrenceMatrix(embed_ts1, embed_ts2, .05; fixedrate = true);

# Step 6 ------ Recurrence quantification

# note that NamedTuple will soon be depreciated. Need to adjust for Dict
crqaOUT = rqa(CRP; lmin = 20, thieller = 0);
crqaOUT = (; crqaOUT...) # convert Dict to NamedTuple

#jrqaOUT = rqa(JRP; lmin = 20, thieller = 0);

# diagonal crp profile
lags = 120;
diag_RR = [];

for i = (-1*lags):1:lags # caluculate diagonal recurrences at lag i
    push!(diag_RR, 100*sum(diag(CRP,i))/length(diag(CRP,i)));
end

# Step 7 ------- Construct output data frame

# get the obtained parameters
parameters = (delay = tau, embed = D);

# join parameters and rqaOUT
df_out = [merge(parameters, rqaOUT)] # be careful not to duplicate column names

#optionally save to csv
CSV.write("crqa_output.csv", df_out)

# Step 8 ------ Diagnostic Plots

## DO NOT ACCEPT ANY VALUES WITHOUT LOOKING AT THE PLOTS!!!!

# movement timeseries plots
ts1_plot = plot(ts1, labels = ts1_names);
ts2_plot = plot(ts2, labels = ts2_names);

# plot of the ami function
ami_plot_ts1 = plot(plotMI_ts1[:,:], label = ts1_names);
plot!(ami_plot_ts1, [τ], seriestype="vline", linestyle=:dash,label = "delay selected")

ami_plot_ts2 = plot(plotMI_ts2[:,:], label = ts2_names);
plot!(ami_plot_ts2, [τ], seriestype="vline", linestyle=:dash,label = "delay selected")

# false nearest neighbors plot
fnn_plot_ts1 = plot(fnnPerc_ts1, label = "%FNN");
plot!(fnn_plot_ts1, [D], seriestype="vline", linestyle=:dash,label = "embed selected")

fnn_plot_ts2 = plot(fnnPerc_ts2, label = "%FNN");
plot!(fnn_plot_ts2, [D], seriestype="vline", linestyle=:dash,label = "embed selected")

# recurrence plot
xs, ys = coordinates(CRP);
rec_plot = scatter(xs, ys, markersize = 0.1, markercolor = :black, legend = false);

# diagonal recurrence plot
dcrp = plot(-lags:1:lags,diag_RR, legend = false);

# combine plots
diagnostic_plots = plot(ts1_plot, ami_plot_ts1, fnn_plot_ts1, rec_plot, ts2_plot, ami_plot_ts2, fnn_plot_ts2, dcrp, layout = (2,4));
#diagnostic_plots = plot(rec_plot, ami_plot,fnn_plot, layout = l);
plot!(diagnostic_plots, size=(2000,1000));

diagnostic_plots_filename = string("cross_rec_diagnostic_plot.png")
png(diagnostic_plots, diagnostic_plots_filename)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Deprecated functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


text:
"""
Create a 1Q GST sequence for a csv file of single qubit Clifford numbers

create_1Q_GST_seqs(file, q; Id_length=50e-9, num_cals=500)

file	 : path to csv file
q		 : name of qubit in scope
Id_length : length of the null {} pulse in seconds this is typically 0
num_cals	 : number of calibration sequences for 0 and 1

"""
function create_1Q_GST_seqs(file, q; Id_length=0, num_cals=500)

	# create pulses once in a library
	m = MEAS(q)
	pulse_lib = [DiAC(q,x) for x in 1:24]
	insert!(pulse_lib, 1, QGL.PulseBlock(Id(q,Id_length)))

	seqs = Vector{Vector{QGL.SequenceEntry}}()
	f = open(file, "r")
	for ln in eachline(f)
		# convert from string to Int and 0-indexed to 1-indexed
		pulse_nums = map(x -> parse(Int, x) + 1, split(ln, ','))
		seq = Vector{QGL.SequenceEntry}(length(pulse_nums) + 1)
		map!(x -> pulse_lib[x], seq, pulse_nums)
		seq[end] = m
		push!(seqs, seq)
	end
	cals = cal_seqs((q,),num_repeats=num_cals)
	seqs = vcat(seqs, cals)
	return seqs
end
"""
Create a 1Q GST sequence for a csv file of single qubit Clifford numbers
using diatomic pulses

create_1Q_GST_seqs(file::Array{String,1}, q; Id_length=50e-9, num_cals=500)

file	 : a string array with the paths to multiple csv files
q		 : name of qubit in scope
Id_length : length of the null {} pulse in seconds this is typically 0
num_cals	 : number of calibration sequences for 0 and 1
"""
function create_1Q_GST_seqs(file::Array{String, 1}, q; Id_length=0,
	num_cals=500)

	# create pulses once in a library
	m = MEAS(q)
	pulse_lib = [DiAC(q,x) for x in 1:24]
	insert!(pulse_lib, 1, QGL.PulseBlock(Id(q,Id_length)))

	seqs = Vector{Vector{QGL.SequenceEntry}}()
	cals = cal_seqs((q,),num_repeats=num_cals)
	for fname in file
		f = open(fname, "r")
		for ln in eachline(f)
			# convert from string to Int and 0-indexed to 1-indexed
			pulse_nums = map(x -> parse(Int, x) + 1, split(ln, ','))
			seq = Vector{QGL.SequenceEntry}(length(pulse_nums) + 1)
			map!(x -> pulse_lib[x], seq, pulse_nums)
			seq[end] = m
			push!(seqs, seq)
		end
		seqs = vcat(seqs, cals)
	end
	return seqs
end
"""
Create a 1Q GST sequence for a csv file of single qubit Clifford numbers using
diatomic pulses.  The experimental sequence specified by the csv file is
repeated 10 times including cals.

create_1Q_GST10_seqs(file, q; Id_length=50e-9, num_cals=500)

file	 : path to csv file
q		 : name of qubit in scope
Id_length : length of the null {} pulse in seconds this is typically 0
num_cals	 : number of calibration sequences for 0 and 1
"""
function create_1Q_GST10_seqs(file, q; Id_length=0, num_cals=500)

	# create pulses once in a library
	m = MEAS(q)
	pulse_lib = [DiAC(q,x) for x in 1:24]
	insert!(pulse_lib, 1, QGL.PulseBlock(Id(q,Id_length)))
	seqs = Vector{Vector{QGL.SequenceEntry}}()
	f = open(file, "r")
	for ln in eachline(f)
		# convert from string to Int and 0-indexed to 1-indexed
		pulse_nums = map(x -> parse(Int, x) + 1, split(ln, ','))
		seq = Vector{QGL.SequenceEntry}(length(pulse_nums) + 1)
		map!(x -> pulse_lib[x], seq, pulse_nums)
		seq[end] = m
		push!(seqs, seq)
	end
	cals = cal_seqs((q,),num_repeats=num_cals)
	seqs = vcat(seqs, cals)

	# x10
	seqs = repeat(seqs, outer=[10])

	return seqs
end


text:
#=

Chase Coleman (ccoleman@stern.nyu.edu)
Spencer Lyon (spencer.lyon@stern.nyu.edu)

References
----------

Judd, Kenneth L, Lilia Maliar, Serguei Maliar, and Rafael Valero. 2013.
    "Smolyak Method for Solving Dynamic Economic Models: Lagrange
    Interpolation, Anisotropic Grid and Adaptive Domain".

Krueger, Dirk, and Felix Kubler. 2004. "Computing Equilibrium in OLG
    Models with Stochastic Production." Journal of Economic Dynamics and
    Control 28 (7) (April): 1411-1436.

=#

# include helper file
include("smol_util.jl")

struct Smolyak <: BasisFamily end

struct SmolyakParams{T,Tmu<:IntSorV} <: BasisParams
    d::Int
    mu::Tmu
    a::Vector{T}
    b::Vector{T}

    # internal fields
    inds::Vector{Vector{Int}}  # Smolyak indices
    pinds::Matrix{Int64}  # Polynomial indices

    function SmolyakParams{T,Tmu}(
            d::Int, mu::Tmu, a::Vector{T}, b::Vector{T}
        ) where {T,Tmu}
        d < 2 && error("You passed d = $d. d must be greater than 1")
        if length(mu) > 1
            # working on building an anisotropic grid
            if length(mu) != d
                error("mu must have d elements. It has $(length(mu))")
            end

            if any(mu .< 1)
                error("All elements of mu must be greater than 1 (recieved $mu)")
            end
        else
            mu < 1 && error("You passed mu = $mu. mu must be greater than 1")
        end

        inds = smol_inds(d, mu)
        pinds = poly_inds(d, mu, inds)

        new{T,Tmu}(d, mu, a, b, inds, pinds)
    end
end

function SmolyakParams(d::Int, mu::Tmu, a::Vector{T}=fill(-1.0, d),
                       b::Vector{T}=fill(1.0, d)) where {T,Tmu}
    SmolyakParams{T,Tmu}(d, mu, a, b)
end

# add methods to helper routines from other files
smol_inds(sp::SmolyakParams) = sp.inds

function poly_inds(sp::SmolyakParams, inds::Vector{Vector{Int}}=sp.inds)
    poly_inds(sp.d, sp.mu, inds)
end

function build_grid(sp::SmolyakParams, inds::Vector{Vector{Int}}=sp.inds)
    build_grid(sp.d, sp.mu, inds)
end

function build_B!(out::AbstractMatrix{T}, sp::SmolyakParams,
                  pts::Matrix{Float64}, b_inds::Matrix{Int64}=sp.pinds) where T
    build_B!(out, sp.d, sp.mu, pts, b_inds)
end

function build_B(sp::SmolyakParams, pts::Matrix{Float64},
                 b_inds::Matrix{Int64}=sp.pinds)
    build_B!(Array{Float64}(size(pts, 1), size(b_inds, 1)), sp, pts, b_inds)
end

function dom2cube!(out::AbstractMatrix, pts::AbstractMatrix,
                   sp::SmolyakParams)
    dom2cube!(out, pts, sp.a, sp.b)
end
dom2cube(pts::AbstractMatrix, sp::SmolyakParams) = dom2cube(pts, sp.a, sp.b)

function cube2dom!(out::AbstractMatrix, pts::AbstractMatrix,
                   sp::SmolyakParams)
    cube2dom!(out, pts, sp.a, sp.b)
end
cube2dom(pts::AbstractMatrix, sp::SmolyakParams) = cube2dom(pts, sp.a, sp.b)

## BasisParams interface
# define these methods on the type, the instance version is defined over
# BasisParams
family(::Type{T}) where {T<:SmolyakParams} = Smolyak
family_name(::Type{T}) where {T<:SmolyakParams} = "Smolyak"
Base.eltype(::Type{SmolyakParams{T1,T2}}) where {T1,T2} = T1

# methods that only make sense for instances
Base.min(p::SmolyakParams) = p.a
Base.max(p::SmolyakParams) = p.b
Base.ndims(sp::SmolyakParams) = sp.d

function Base.show(io::IO, p::SmolyakParams)
    m = string("Smolyak interpolation parameters in $(p.d) dimensions",
               " from $(p.a) × $(p.b)")
    print(io, m)
end

# TODO: fix this
function Base.length{T,Ti<:Integer}(sp::SmolyakParams{T,Ti})::Int
    d, mu = sp.d, sp.mu
    mu == 1 ? 2d - 1 :
    mu == 2 ? Int(1 + 4d + 4d*(d-1)/2 ) :
    mu == 3 ? Int(1 + 8d + 12d*(d-1)/2 + 8d*(d-1)*(d-2)/6) :
    error("We only know the number of grid points for mu ∈ [1, 2, 3]")
end

function nodes(sp::SmolyakParams)
    dom_grid = build_grid(sp)
    cube2dom!(dom_grid, dom_grid, sp)
end

function evalbase(sp::SmolyakParams, x::AbstractVector, order=0)
    evalbase(sp, reshape(x, 1, :), order)
end

function evalbase(sp::SmolyakParams, x::AbstractArray, order=0)
    cube_pts = dom2cube(x, sp)
    build_B(sp, cube_pts, sp.pinds)
end


text:
using HeatEquation
using Test

@testset begin

    ncols, nrows = 1000, 1000
    nsteps = 500

    curr, prev = initialize(ncols, nrows)
    # initial average temperature
    @test average_temperature(curr) ≈ 59.76482

    # run simulation
    simulate!(curr, prev, nsteps)

    # final average temperature
    @test average_temperature(curr) ≈ 58.80755
end


text:
function mflows(cf, p, sval, rval, tracs, maxt; digits=ROUND_DIGITS)
    newtracs = [0., 0.]

    function update!(i, f)
        #newcf[i] = cf[i][tracs[i]] + f
        newtracs[i] = cvalarg(cf[i], cf[i][tracs[i]] + f, tracs[i], maxt[i])
    end

    if rval >= sum(sval)
        update!(1, sval[1]) #free
        update!(2, sval[2]) #free
    elseif sval[1] < 0.5*rval
        update!(1, sval[1]) #free
        update!(2, rval-sval[1]) #cong
    elseif sval[2] < 0.5*rval
        update!(1, rval-sval[2]) #free
        update!(2, sval[2]) #cong
    else
        update!(1, 0.5*rval) #cong
        update!(2, 0.5*rval) #cong
    end
    newtracs = round.(newtracs, digits=ROUND_DIGITS)

    newcf = round.([cf[i][newtracs[i]] for i in 1:2], digits=ROUND_DIGITS)

    newp = zeros(size(p[1])[2:end]...)
    for i in 1:2
        rng = tracs[i]:newtracs[i]
        if size(rng) > 0.
            newp .+= squeezesum(mucsum(cf[i], rng) .* p[i][indexify(rng),..], dims=1)
        end
    end
    if sum(newp) == 0.
        newp .= 1. /prod(size(newp))
    else
        newp ./= sum(newp)
    end

    return (newcf, newp, newtracs)
end

function dflows(cf, p, sr, sval, rval, trac, maxt; digits=ROUND_DIGITS)
    rng = trac:maxt
    fcut = [(mucsum(cf, rng) .* p[indexify(rng),..]) .* expand(sr[i], dims=1) for i in 1:2] # flow cut
    sva = valarg(sum(fcut), sval, 1)
    rva = min(valarg.(fcut, rval, 1)...)
    va = min(sva, rva)

    tmp = [squeezesum(fcut[i][0.:va,..], dims=1) for i in 1:2]
    newp = [(sum(tmptmp) == 0.) ? ones(size(tmptmp)...)*(1. /prod(size(tmptmp))) : (tmptmp ./ sum(tmptmp)) for tmptmp in tmp]

    newtrac = cvalarg(cf, cf[trac] + sum(sum(tmp)), trac, maxt)
    newtrac = round(newtrac, digits=ROUND_DIGITS)

    newcf = round(cf[newtrac], digits=ROUND_DIGITS)

    return (newcf, newp, newtrac, sum.(tmp))#[sum(tmptmp) for tmptmp in tmp])
end


text:
@inline function _countPotentialsSerial(potentials::Vector{Float64},
  ζcounts::Vector{Int64}, rpio::RPIO{Int64Particle}, p::Int64, G::F) where F<:Function
  v::Float64 = 0.0
  px::Int64Particle = rpio.scratchInt64Particle[1]
  for i in 1:length(ζcounts)
    @inbounds if ζcounts[i] > 0
      px.x = i-1
      @inbounds potentials[i] = ζcounts[i] * G(p, px)
      # @inbounds potentials[i] = ζcounts[i] * G(p, i-1)
      @inbounds v += potentials[i]
    else
      @inbounds potentials[i] = 0.0
    end
  end
  return v
end

@inline function _countPotentialsParallel(potentials::Vector{Float64},
  ζcounts::Vector{Int64}, rpio::RPIO{Int64Particle}, p::Int64, G::F) where F<:Function
  scratch::Vector{Float64} = rpio.scratch
  M::Int64 = length(ζcounts)
  nt::Int64 = rpio.threads
  b::Int64 = div(M, nt)
  Threads.@threads for i in 1:nt
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : M
    v::Float64 = 0.0
    px::Int64Particle = rpio.scratchInt64Particle[i]
    for j in start:finish
      @inbounds if ζcounts[j] > 0
        px.x = j-1
        @inbounds potentials[j] = ζcounts[j] * G(p, px)
        # @inbounds potentials[j] = ζcounts[j] * G(p, j-1)
        @inbounds v += potentials[j]
      else
        @inbounds potentials[j] = 0.0
      end
    end
    @inbounds scratch[i] = v
  end
  return sum(scratch)
end

const PARALLEL_COUNTS_THRESHOLD = 256

function maximumI64(ps::Vector{Int64Particle})
  m::Int64 = 0
  for i in eachindex(ps)
    @inbounds m = max(m, ps[i].x)
  end
  return m
end

@inline function _countPotentials(potentials::Vector{Float64},
  ζcounts::Vector{Int64}, model::Model, rpio::RPIO{Int64Particle}, p::Int64)
  if rpio.threads > 1 && length(ζcounts) > PARALLEL_COUNTS_THRESHOLD
    return _countPotentialsParallel(potentials, ζcounts, rpio, p, model.G)
  else
    return _countPotentialsSerial(potentials, ζcounts, rpio, p, model.G)
  end
end

@inline function _getCountsSerial(rpio::RPIO{Int64Particle})
  particles::Vector{Int64Particle} = rpio.particles
  ζmax::Int64 = maximumI64(particles)
  ζcounts::Vector{Int64} = zeros(Int64, ζmax + 1)
  for i in 1:length(particles)
    @inbounds ζcounts[particles[i].x+1] += 1
  end
  return ζcounts
end

@inline function _threadCounts(threadCounts::Vector{Vector{Int64}},
  particles::Vector{Int64Particle}, ζmax::Int64)
  M::Int64 = length(particles)
  nt::Int64 = length(threadCounts)
  b::Int64 = div(M, nt)
  Threads.@threads for i in 1:nt
    for j in 1:ζmax + 1
      threadCounts[i][j] = 0
    end
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : M
    for j in start:finish
      @inbounds threadCounts[i][particles[j].x+1] += 1
    end
  end
end

@inline function _sumCounts(ζcounts::Vector{Int64},
  threadCounts::Vector{Vector{Int64}})
  M::Int64 = length(ζcounts)
  nt::Int64 = length(threadCounts)
  b::Int64 = div(M, nt)
  Threads.@threads for i in 1:nt
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : M
    for j in start:finish
      @inbounds ζcounts[j] = 0
      for k in 1:nt
        @inbounds ζcounts[j] += threadCounts[k][j]
      end
    end
  end
end

@inline function _maxParticlesParallel(rpio::RPIO{Int64Particle})
  nt::Int64 = rpio.threads
  particles::Vector{Int64Particle} = rpio.particles
  N::Int64 = rpio.N
  b::Int64 = div(N, nt)
  scratchInt::Vector{Int64} = rpio.scratchInt
  Threads.@threads for i in 1:nt
    start::Int64 = (i-1)*b + 1
    finish::Int64 = i != nt ? i*b : N
    v::Int64 = 0
    for j in start:finish
      @inbounds v = max(v, particles[j].x)
    end
    @inbounds scratchInt[i] = v
  end
  return maximum(scratchInt)
end

@inline function _getCountsParallel(rpio::RPIO{Int64Particle})
  ζmax::Int64 = _maxParticlesParallel(rpio)
  ζcounts::Vector{Int64} = Vector{Int64}(undef, ζmax + 1)
  threadCounts::Vector{Vector{Int64}} =
    Vector{Vector{Int64}}(undef, rpio.threads)
  for i in 1:rpio.threads
    threadCounts[i] = Vector{Int64}(undef, ζmax + 1)
  end
  _threadCounts(threadCounts, rpio.particles, ζmax)
  _sumCounts(ζcounts, threadCounts)
  return ζcounts
end

@inline function _getCounts(rpio::RPIO{Int64Particle})
  if rpio.threads > 1
    return _getCountsParallel(rpio)
  else
    return _getCountsSerial(rpio)
  end
end


text:
module API

import ..Accounts: Domain

using .Domain

export create

"""
    create(address_type::AddressType, address::String)::Address
    create(name::String, addresses::Array{Address,1})::Person
    create(name::String)::Person

Create an Address or a Person object.

# Example
```julia
julia> using Accounts

julia> address_email = create(EMAIL, "donald@duckcity.com")
Accounts.Domain.Address("7763679977726623090", Dates.DateTime("2020-09-28T10:56:29.461"), EMAIL, "donald@duckcity.com")

julia> donald = create("Donald Duck", [address_email])
Accounts.Domain.Person("14456927583164318539", Dates.DateTime("2020-09-28T10:56:29.997"), "Donald Duck", Accounts.Domain.Address[Accounts.Domain.Address("7763679977726623090", Dates.DateTime("2020-09-28T10:56:29.461"), EMAIL, "donald@duckcity.com")])

julia> donald.id
"14456927583164318539"

julia> donald.created
2020-09-28T10:56:29.997

julia> donald.name
"Donald Duck"

julia> donald.addresses
1-element Array{Accounts.Domain.Address,1}:
 Accounts.Domain.Address("7763679977726623090", Dates.DateTime("2020-09-28T10:56:29.461"), EMAIL, "donald@duckcity.com")
```
"""
function create end

create(address_type::AddressType, address::String)::Address = Address(address_type, address)

create(name::String, addresses::Array{Address,1})::Person = Person(name, addresses)

create(name::String)::Person = Person(name)

end


text:
@testset "Manopt.jl Error Measures" begin
    M = Sphere(2)
    N = PowerManifold(M, NestedPowerRepresentation(), 2)
    using Random: seed!
    seed!(42)
    d = Manifolds.uniform_distribution(M, [1.0, 0.0, 0.0])
    w = rand(d)
    x = rand(d)
    y = rand(d)
    z = rand(d)
    a = [w, x]
    b = [y, z]
    @test meanSquaredError(M, x, y) == distance(M, x, y)^2
    @test meanSquaredError(N, a, b) == 1 / 2 * (distance(M, w, y)^2 + distance(M, x, z)^2)
    @test meanAverageError(M, x, y) == distance(M, x, y)
    @test meanAverageError(N, a, b) == 1 / 2 * sum(distance.(Ref(M), a, b))
end


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter, CutPruners

makedocs()

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/JuliaPolyhedra/CutPruners.jl.git",
    julia  = "release"
)


text:
let
    using UnsteadyFlowSolvers

    # Set up a single vortex
    start_vc = 0.02 # Core radius
    start_vort = 1.
    re = 1e3
    nu = start_vort / re
    field = TwoDVFlowField(nu)
    vort = TwoDVVort(0., 0., start_vort, start_vc, 0., 0., 0.)
    push!(field.tev, vort)
    vort = TwoDVVort(0., 0.01, start_vort, start_vc, 0., 0., 0.)
    push!(field.tev, vort)

    # Solver settings
    dt = 0.015
    t0 = 0.
    tmax = 1.5
    t = t0
    nsteps = Int(ceil((tmax-t)/dt)) + 1

    for i = 1:nsteps
        t += dt
        UnsteadyFlowSolvers.wakeroll(field, dt)
        println("Step\t",i,":\tCore 1\t",field.tev[1].vc, "\tCore 2\t", field.tev[2].vc)
        println("\t\t", "VelX 1\t",field.tev[1].vx, "\tVelX 2\t", field.tev[2].vx)
        println("\t\t", "VelZ 1\t",field.tev[1].vz, "\tVelZ 2\t", field.tev[2].vz)
    end
end


text:
using juleva
using Test

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

max_stars_repo_name: denisbertini/juleva

max_issues_repo_name: denisbertini/juleva

max_forks_repo_name: denisbertini/juleva


text:
boundary_names(mesh::AbstractMesh) = boundary_names(domain(mesh))
coordinate_type(mesh::AbstractMesh) = coordinate_type(domain(mesh))

"""
    nelements(mesh::AbstractMesh)

The number of elements in the mesh.
"""
nelements(mesh::AbstractMesh) = length(elements(mesh))

"""
    i, ξs = split_refcoord(ξ, n)

Given a reference coordinate `ξ` in the interval `[-1, 1]``, divide the interval
into `n` evenly spaced subintervals, and return the index of the subinterval
(`i`), and the position in the subinterval (`ξs`), normalized to normalized `[-1, 1]`.
"""
function split_refcoord(ξ, n)
    ξn = ξ * n
    # we use Julia's range lookup here, which avoids intermediate rounding errors.
    i = min(searchsortedlast((-n):2:n, ξn), n)
    ξ = ξn - (2 * i - n - 1)
    return (i, ξ)
end

"""
    Meshes.SharedVertices(mesh, elem, vert)

An iterator over (element, vertex) pairs that are shared with `(elem,vert)`.
"""
struct SharedVertices{M <: AbstractMesh, E}
    mesh::M
    elem::E
    vert::Int
end
Base.IteratorSize(::Type{<:SharedVertices}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:SharedVertices}) = Base.HasEltype()
Base.eltype(::Type{SharedVertices{M, E}}) where {M, E} = Tuple{E, Int}

function Base.iterate(vertiter::SharedVertices)
    velem = vertiter.elem
    vvert = vertiter.vert
    ccw = false
    # return initial (element, vertex)
    return (velem, vvert), (velem, vvert, ccw)
end
function Base.iterate(vertiter::SharedVertices, (velem, vvert, ccw))
    # initially we go clockwise (ccw == false), we go to the face == vert
    # if ccw, then go to face = vert - 1
    vface = ccw ? mod1(vvert - 1, 4) : vvert
    if is_boundary_face(vertiter.mesh, velem, vface)
        if ccw
            # have already gone both directions: all done
            return nothing
        end
        # try counter-clockwise
        velem = vertiter.elem
        vvert = vertiter.vert
        ccw = true
        return Base.iterate(vertiter, (velem, vvert, ccw))
    end
    opelem, opface, reversed = opposing_face(vertiter.mesh, velem, vface)
    velem = opelem
    vvert = ccw ? opface : mod1(opface + 1, 4)
    if velem == vertiter.elem && vvert == vertiter.vert
        # we're back at where we started: all done
        return nothing
    end
    return (velem, vvert), (velem, vvert, ccw)
end

"""
    Meshes.linearindices(elemorder)

Given a data structure `elemorder[i] = elem` that orders elements, construct the
inverse map from `orderindex = linearindices(elemorder)` such that
`orderindex[elem] = i`.

This will try to use the most efficient structure available.
"""
linearindices(elemorder::CartesianIndices) = LinearIndices(elemorder)
function linearindices(elemorder::AbstractVector{CartesianIndex})
    cmax = maximum(elemorder)
    L = zeros(Int, cmax.I)
    for (i, c) in enumerate(elemorder)
        L[c] = i
    end
    return L
end
function linearindices(elemorder)
    orderindex = Dict{eltype(elemorder), Int}()
    for (i, elem) in elemorder
        orderindex[elem] = i
    end
    return orderindex
end

"""
    M = Meshes.face_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a `Bool`-valued `SparseCSCMatrix` containing the face connections of
`mesh`. Elements are indexed according to `elemorder`.

Note that `M[i,i] == true` only if two distinct faces of element `i` are connected.
"""
function face_connectivity_matrix(
    mesh::AbstractMesh,
    elemorder = elements(mesh),
    orderindex = linearindices(elemorder),
)
    m = n = length(elemorder)
    I = Int[]
    J = Int[]
    for (i, elem) in enumerate(elemorder)
        for face in 1:4
            if is_boundary_face(mesh, elem, face)
                continue
            end
            opelem, opface, reversed = Meshes.opposing_face(mesh, elem, face)
            j = orderindex[opelem]
            push!(I, i)
            push!(J, j)
        end
    end
    V = trues(length(I))
    return SparseArrays.sparse(I, J, V, m, n)
end

"""
    M = Meshes.vertex_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a `Bool`-valued `SparseCSCMatrix` containing the vertex connections of
`mesh`. Elements are indexed according to `elemorder`.

Note that `M[i,i] == true` only if two distinct vertices of element `i` are connected.
"""
function vertex_connectivity_matrix(
    mesh::AbstractMesh,
    elemorder = elements(mesh),
    orderindex = linearindices(elemorder),
)
    m = n = length(elemorder)
    I = Int[]
    J = Int[]
    for (i, elem) in enumerate(elemorder)
        for vert in 1:4
            for (velem, vvert) in SharedVertices(mesh, elem, vert)
                if velem == elem && vvert == vert
                    continue
                end
                j = orderindex[velem]
                push!(I, i)
                push!(J, j)
            end
        end
    end
    V = trues(length(I))
    return SparseArrays.sparse(I, J, V, m, n)
end

"""
    Meshes.coordinates(mesh, elem, vert::Int)
    Meshes.coordinates(mesh, elem, (ξ1, ξ2)::Tuple)

Return the physical coordinates of a point in an element `elem` of `mesh`. The
position of the point can either be a vertex number `vert` or a pair of coordinates
`(ξ1, ξ2)` in the reference element.
"""
function coordinates(mesh::AbstractMesh2D, elem, vert::Integer)
    FT = Domains.float_type(domain(mesh))
    ξ1 = (vert == 1 || vert == 4) ? FT(-1) : FT(1)
    ξ2 = (vert == 1 || vert == 2) ? FT(-1) : FT(1)
    coordinates(mesh, elem, (ξ1, ξ2))
end

"""
    elem, (ξ1, ξ2) = Meshes.containing_element(mesh::AbstractCubedSphere, coord)

Return the element (`elem`) in `mesh` containing the coordinate `coord`, and the position (`ξ1, ξ2`)
within the reference element.
"""
function containing_element end


text:
"""
    show_collection(io, collection, ::Type{H<:ObjectHandle})

Given a collection-like object, (`Symbols`, `DynamicLinks`, ``)
"""
function show_collection(io::IO, stuff::ST, ::Type{H}) where {ST} where {H <: ObjectHandle}
    h_str = format_string(H)

    # Get the type name, removing the `ELF`/`MachO`/`COFF` at the beginning
    type_name = if isa(ST, UnionAll)
        string(ST.body.name.name)
    else
        string(ST.name.name)
    end

    if startswith(type_name, h_str)
        type_name = type_name[length(h_str)+1:end]
    end

    # Start by outputting the container type
    print(io, h_str," ", type_name)

    # Next, print either the entire thing, or a limited subsection of it
    limited = get(io, :limit, false)
    if limited && length(stuff) > 20
        for idx in 1:10
            print(io, "\n  [$(idx)] ")
            show(IOContext(io, :compact => true), stuff[idx])
        end
        print(io, "\n   \u2026")
        for idx in length(stuff)-10:length(stuff)
            print(io, "\n  [$(idx)] ")
            show(IOContext(io, :compact => true), stuff[idx])
        end
    else
        for idx in 1:length(stuff)
            print(io, "\n  [$(idx)] ")
            show(IOContext(io, :compact => true), stuff[idx])
        end
    end
end


text:
function rk4(f, t, y, dt, P, k)
    k1 = f(t, y, P)
    k2 = f(t + 1/2*dt, y + 1/2*dt*k1, P)
    k3 = f(t + 1/2*dt, y + 1/2*dt*k2, P)
    k4 = f(t + dt, y + dt*k3, P)
    y + dt*(k1+2*k2+2*k3+k4)/6
end

struct Rk4 <: ODESolverType end

#=
        FOUR components
=#
function update(::Rk4, fs, t, A, B, C, D, dt, P, tableau=NaN)
    kA1 = update(fs[1], t, A, B, C, D, P)
    kB1 = update(fs[2], t, A, B, C, D, P)
    kC1 = update(fs[3], t, A, B, C, D, P)
    kD1 = update(fs[4], t, A, B, C, D, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1
    Ci = C + 1/2*dt*kC1
    Di = D + 1/2*dt*kD1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kC2 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kD2 = update(fs[4], t + 1/2*dt, Ai, Bi, Ci, Di, P)

    Ai = A + 1/2*dt*kA2
    Bi = B + 1/2*dt*kB2
    Ci = C + 1/2*dt*kC2
    Di = D + 1/2*dt*kD2

    kA3 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kB3 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kC3 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, Di, P)
    kD3 = update(fs[4], t + 1/2*dt, Ai, Bi, Ci, Di, P)

    Ai = A + dt*kA3
    Bi = B + dt*kB3
    Ci = C + dt*kC3
    Di = D + dt*kD3

    kA4 = update(fs[1], t + dt, Ai, Bi, Ci, Di, P)
    kB4 = update(fs[2], t + dt, Ai, Bi, Ci, Di, P)
    kC4 = update(fs[3], t + dt, Ai, Bi, Ci, Di, P)
    kD4 = update(fs[4], t + dt, Ai, Bi, Ci, Di, P)

    (A + dt*(kA1 + 2*kA2 + 2*kA3 + kA4)/6, B + dt*(kB1 + 2*kB2 + 2*kB3 + kB4)/6,
     C + dt*(kC1 + 2*kC2 + 2*kC3 + kC4)/6, D + dt*(kD1 + 2*kD2 + 2*kD3 + kD4)/6)
end

#=
        THREE components
=#
function update(::Rk4, fs, t, A, B, C, dt, P, tableau=NaN)
    kA1 = update(fs[1], t, A, B, C, P)
    kB1 = update(fs[2], t, A, B, C, P)
    kC1 = update(fs[3], t, A, B, C, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1
    Ci = C + 1/2*dt*kC1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, P)
    kC2 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, P)

    Ai = A + 1/2*dt*kA2
    Bi = B + 1/2*dt*kB2
    Ci = C + 1/2*dt*kC2

    kA3 = update(fs[1], t + 1/2*dt, Ai, Bi, Ci, P)
    kB3 = update(fs[2], t + 1/2*dt, Ai, Bi, Ci, P)
    kC3 = update(fs[3], t + 1/2*dt, Ai, Bi, Ci, P)

    Ai = A + dt*kA3
    Bi = B + dt*kB3
    Ci = C + dt*kC3

    kA4 = update(fs[1], t + dt, Ai, Bi, Ci, P)
    kB4 = update(fs[2], t + dt, Ai, Bi, Ci, P)
    kC4 = update(fs[3], t + dt, Ai, Bi, Ci, P)

    (A + dt*(kA1 + 2*kA2 + 2*kA3 + kA4)/6, B + dt*(kB1 + 2*kB2 + 2*kB3 + kB4)/6,
     C + dt*(kC1 + 2*kC2 + 2*kC3 + kC4)/6)
end

#=
        TWO components
=#
function update(::Rk4, fs, t, A, B, dt, P, tableau=NaN)
    kA1 = update(fs[1], t, A, B, P)
    kB1 = update(fs[2], t, A, B, P)

    Ai = A + 1/2*dt*kA1
    Bi = B + 1/2*dt*kB1

    kA2 = update(fs[1], t + 1/2*dt, Ai, Bi, P)
    kB2 = update(fs[2], t + 1/2*dt, Ai, Bi, P)

    Ai = A + 1/2*dt*kA2
    Bi = B + 1/2*dt*kB2

    kA3 = update(fs[1], t + 1/2*dt, Ai, Bi, P)
    kB3 = update(fs[2], t + 1/2*dt, Ai, Bi, P)

    Ai = A + dt*kA3
    Bi = B + dt*kB3

    kA4 = update(fs[1], t + dt, Ai, Bi, P)
    kB4 = update(fs[2], t + dt, Ai, Bi, P)

    (A + dt*(kA1 + 2*kA2 + 2*kA3 + kA4)/6, B + dt*(kB1 + 2*kB2 + 2*kB3 + kB4)/6)
end


text:
export iauEe06a
"""
Equation of the equinoxes, compatible with IAU 2000 resolutions and
IAU 2006/2000A precession-nutation.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double    TT as a 2-part Julian Date (Note 1)

Returned (function value):
               double    equation of the equinoxes (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The result, which is in radians, operates in the following sense:

      Greenwich apparent ST = GMST + equation of the equinoxes

Called:
   iauAnpm      normalize angle into range +/- pi
   iauGst06a    Greenwich apparent sidereal time, IAU 2006/2000A
   iauGmst06    Greenwich mean sidereal time, IAU 2006

Reference:

   McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
   IERS Technical Note No. 32, BKG

This revision:  2008 May 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEe06a(double date1, double date2)

function iauEe06a(date1::Real, date2::Real)
   return ccall((:iauEe06a, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, date1), convert(Float64, date2))
end


text:
using BenchmarkTools
using CSV
using Distributions
using LACosmic
using ProgressLogging
using PSFModels: Gaussian
using PyCall
using Random
using Statistics

ENV["OMP_NUM_THREADS"] = 1
@info "OMP NUM THREADS=$(ENV["OMP_NUM_THREADS"])"

scrap = pyimport("astroscrappy")

# generate data
rng = Random.MersenneTwister(444219)

function make_data(rng, T, N; N_sources=N ÷ 10, N_cosmics=N ÷ 10)
    imdata = fill(T(200), (N, N))

    # Add some fake sources
    for _ in 1:N_sources
        x = rand(rng, Uniform(1, N + 1))
        y = rand(rng, Uniform(1, N + 1))
        brightness = rand(rng, Uniform(1000, 30000))
		model = Gaussian(T; x=x, y=y, fwhm=3.5, amp=brightness)
        imdata .+= model[axes(imdata)...]
	end

    # Add the poisson noise
    imdata .= rand.(rng, Poisson.(imdata))

    # Add readnoise
    imdata .+= rand(rng, Normal(0, 10), (N, N))

	clean_image = T.(imdata)

    # Add Nc fake cosmic rays
    crmask = falses((N, N))
	for _ in 1:N_cosmics
    	cr_x = round(Int, rand(rng, Uniform(6, N - 5)))
    	cr_y = round(Int, rand(rng, Uniform(6, N - 5)))
    	cr_brightnesses = rand(rng, Uniform(1000, 30000))
    	imdata[cr_y, cr_x] += T(cr_brightnesses)
    	crmask[cr_y, cr_x] = true
	end

    # Make a mask where the detected cosmic rays should be
    return (image=imdata, clean_image, mask=crmask)
end

data_sizes = 101:100:1001

results = []
@progress for N_data in data_sizes
    data = make_data(rng, Float64, N_data)

    bench_jl = @benchmark lacosmic($(data.image);
        sigma_clip=6,
        contrast=5,
        neighbor_thresh=1)
    bench_py = @benchmark scrap.detect_cosmics($(data.image);
        sigclip=6,
        objlim=5,
        sigfrac=1,
        readnoise=0,
        sepmed=false)

    res = (N_data, t_python=median(bench_py).time/1e9, t_julia=median(bench_jl).time/1e9)
    @info "" res...
    push!(results, res)
end

CSV.write(joinpath(@__DIR__, "benchmark_results.csv"), results)


text:
#
# Test for ArgMin.jl
#
# author: Atsushi Sakai
#

using ArgMin
using Test

@testset "ArgMin.jl" begin
    include("test_least_square_solver.jl")
    include("test_convex_optimization_solver.jl")
    include("test_projected_newton_qp_solver.jl")
end


text:
import Oceananigans: tupleit

#####
##### Some utilities for tupling
#####

tupleit(::Nothing) = ()
tupleit(t::NamedTuple) = t
tupleit(t::Tuple) = t
tupleit(nt) = tuple(nt)
tupleit(nt::Vector) = tuple(nt...)

parenttuple(obj) = Tuple(f.data.parent for f in obj)

@inline datatuple(obj::Nothing) = nothing
@inline datatuple(obj::AbstractArray) = obj
@inline datatuple(obj::Tuple) = Tuple(datatuple(o) for o in obj)
@inline datatuple(obj::NamedTuple) = NamedTuple{propertynames(obj)}(datatuple(o) for o in obj)
@inline datatuples(objs...) = (datatuple(obj) for obj in objs)


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 5/3
equations = CompressibleEulerEquations2D(gamma)

"""
    initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)

The blob test case taken from
- Agertz et al. (2006)
  Fundamental differences between SPH and grid methods
  [arXiv: astro-ph/0610051](https://arxiv.org/abs/astro-ph/0610051)
"""
function initial_condition_blob(x, t, equations::CompressibleEulerEquations2D)
  # blob test case, see Agertz et al. https://arxiv.org/pdf/astro-ph/0610051.pdf
  # other reference: https://arxiv.org/pdf/astro-ph/0610051.pdf
  # change discontinuity to tanh
  # typical domain is rectangular, we change it to a square
  # resolution 128^2, 256^2
  # domain size is [-20.0,20.0]^2
  # gamma = 5/3 for this test case
  R = 1.0 # radius of the blob
  # background density
  dens0 = 1.0
  Chi = 10.0 # density contrast
  # reference time of characteristic growth of KH instability equal to 1.0
  tau_kh = 1.0
  tau_cr = tau_kh/1.6 # crushing time
  # determine background velocity
  velx0 = 2*R*sqrt(Chi)/tau_cr
  vely0 = 0.0
  Ma0 = 2.7 # background flow Mach number Ma=v/c
  c = velx0/Ma0 # sound speed
  # use perfect gas assumption to compute background pressure via the sound speed c^2 = gamma * pressure/density
  p0 = c*c*dens0/equations.gamma
  # initial center of the blob
  inicenter = [-15,0]
  x_rel = x-inicenter
  r = sqrt(x_rel[1]^2 + x_rel[2]^2)
  # steepness of the tanh transition zone
  slope = 2
  # density blob
  dens = dens0 + (Chi-1) * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  # velocity blob is zero
  velx = velx0 - velx0 * 0.5*(1+(tanh(slope*(r+R)) - (tanh(slope*(r-R)) + 1)))
  return prim2cons(SVector(dens, velx, vely0, p0), equations)
end
initial_condition = initial_condition_blob

surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha
basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.05,
                                         alpha_min=0.0001,
                                         alpha_smooth=true,
                                         variable=pressure)

volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-32.0, -32.0)
coordinates_max = ( 32.0,  32.0)
refinement_patches = (
  (type="box", coordinates_min=(-40.0, -5.0), coordinates_max=(40.0, 5.0)),
  (type="box", coordinates_min=(-40.0, -5.0), coordinates_max=(40.0, 5.0)),
)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                refinement_patches=refinement_patches,
                n_cells_max=100_000,)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 16.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.7)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# @ToDo:
# implement delete operation

# abstract type AbstractBinaryTree end # to be removed
# abstract type AbstractBinaryTreeNode end # to be removed

mutable struct RedBlackTreeNode{T} <: AbstractBinaryTreeNode
    color::Bool # 0 black, 1 red
    parent::RedBlackTreeNode
    left::RedBlackTreeNode
    right::RedBlackTreeNode
    data::T

    function RedBlackTreeNode{T}(content::T) where {T}
        obj = new()
        obj.data = content
        obj.color = true
        return obj
    end
    function RedBlackTreeNode{T}() where {T}
        obj = new()
        obj.color = true
        return obj
    end
end

function RedBlackTreeNode(content::T) where {T}
    return RedBlackTreeNode{T}(content::T)
end

function RedBlackTreeNode()
    return RedBlackTreeNode(undef)
end

function Base.eltype( ::Type{ <: RedBlackTreeNode{T}} ) where {T}
    return T
end

function Base.isempty(node::RedBlackTreeNode)
    return node.data == undef
end

function is_red(node::RedBlackTreeNode)
    return node.color # == true
end
function is_black(node::RedBlackTreeNode)
    return !node.color # == false
end

function set_red!(node::RedBlackTreeNode)
    node.color = true
    return node
end
function set_black!(node::RedBlackTreeNode)
    node.color = false
    return node
end

function is_nil(tree::T, node::N) where {T <: AbstractBinaryTree, N <: AbstractBinaryTreeNode}
    return tree.nil == node
end

mutable struct RedBlackTree{T} <: AbstractBinaryTree
    count::Int64
    node::RedBlackTreeNode{T}
    nil::RedBlackTreeNode{Int}

    function RedBlackTree{T}(content::T) where {T}
        count_       = 1
        node_        = RedBlackTreeNode{T}(content)
        nil_         = RedBlackTreeNode{Int}(-1)
        node_.parent = nil_
        node_.left   = nil_
        node_.right  = nil_
        set_black!(nil_)
        set_black!(node_)
        return new(count_, node_, nil_)
    end
    function RedBlackTree{T}() where {T}
        count_       = 1
        node_        = RedBlackTreeNode{T}()
        nil_         = RedBlackTreeNode{Int}(-1)
        node_.parent = nil_
        node_.left   = nil_
        node_.right  = nil_
        set_black!(nil_)
        set_black!(node_)
        return new(count_, node_, nil_)
    end
end

function RedBlackTree(content::T) where {T}
    return RedBlackTree{T}(content::T)
end

function right_rotate!(tree::T, x::N) where {T <: AbstractBinaryTree, N <: AbstractBinaryTreeNode}

    if isdefined(x, :left)
        alpha = x.left

        if isdefined(x, :parent)
            p = x.parent
            alpha.parent = p
            if isdefined(p, :left) && x == p.left
                p.left = alpha
            elseif isdefined(p, :right)
                p.right = alpha
            end
        end
        if isdefined(alpha, :right)
            x.left             = alpha.right
            alpha.right.parent = x
        end

        alpha.right = x
        x.parent    = alpha

        if is_root(tree, x)
            tree.node = alpha
        end
    end

    return tree
end

function left_rotate!(tree::T, x::N) where {T <: AbstractBinaryTree, N <: AbstractBinaryTreeNode}

    if isdefined(x, :right)
        beta = x.right

        if isdefined(x, :parent)
            p = x.parent
            beta.parent = p
            if isdefined(p, :left) && x == p.left
                p.left = beta
            elseif isdefined(p, :right)
                p.right = beta
            end
        end
        if isdefined(beta, :left)
            x.right          = beta.left
            beta.left.parent = x
        end

        beta.left = x
        x.parent  = beta

        if is_root(tree, x)
            tree.node = beta
        end
    end

    return tree
end

function Base.insert!(tree::T, new_node::N) where {T <: RedBlackTree, N <: RedBlackTreeNode}
    node     = tree.node
    parent   = node
    to_left  = false
    to_right = false

    while true
        if new_node.data < node.data
            if isdefined(node, :left) && !is_nil(tree, node.left)
                child = node.left
            else
                parent  = node
                to_left = true
                break
            end
        else
            if isdefined(node, :right) && !is_nil(tree, node.right)
                child = node.right
            else
                parent   = node
                to_right = true
                break
            end
        end
        parent = node
        node = child
    end

    new_node.parent = parent
    if to_left
        parent.left = new_node
    end
    if to_right
        parent.right = new_node
    end

    new_node.left  = tree.nil
    new_node.right = tree.nil

    set_red!(new_node)
    red_black_insert_fixup!(tree, new_node)

    tree.count += 1
    return tree
end

function red_black_insert_fixup!(tree::T, node::N) where {T <: RedBlackTree, N <: RedBlackTreeNode}
    while is_red(node.parent)

        # parent  = node.parent
        # grandpa = node.parent.parent
        # isdefined(grandpa, :left) &&
        if node.parent == node.parent.parent.left

            y = node.parent.parent.right     # uncle
            if is_red(y)
                set_black!(node.parent)      # parent is black
                set_black!(y)                # uncle is black
                set_red!(node.parent.parent) # grandparent is red
                node = node.parent.parent    # move 2 levels up the hierarchy
            else
                if node == node.parent.right
                    node = node.parent
                    left_rotate!(tree, node)
                end
                set_black!(node.parent)
                set_red!(node.parent.parent)
                right_rotate!(tree, node.parent.parent)
            end

        else
            # same with left and right exchanged

            y = node.parent.parent.left
            if is_red(y)
                set_black!(node.parent)
                set_black!(y)
                set_red!(node.parent.parent)
                node = node.parent.parent
            else
                if node == node.parent.left
                    node = node.parent
                    right_rotate!(tree, node)
                end
                set_black!(node.parent)
                set_red!(node.parent.parent)
                left_rotate!(tree, node.parent.parent)
            end

        end
    end

    set_black!(tree.node)
    set_black!(tree.nil)

    return tree
end

function transplant!(tree::T, old_node::N, new_node::N) where {T <: RedBlackTree, N <: RedBlackTreeNode}

    if isdefined(old_node, :parent)
        parent = old_node.parent
        if old_node == parent.left
            parent.left = new_node
        else
            parent.right = new_node
        end
        new_node.parent = parent
    else
        tree.node = new_node
    end

    return tree
end

function remove!(tree::T, node::N) where {T <: RedBlackTree, N <: RedBlackTreeNode}

    return tree
end

function Base.show(io::IO, tree::RedBlackTree{T}) where {T}
    println("RedBlackTree{$(T)} with properties")
    println(" \t count: $(tree.count)")
    print(" \t node: {")
    show_node(tree.node)
    print("}\n\n")
end

function show_node(node::RedBlackTreeNode)
    print(" color: ")
    if is_red(node)
        print("red, ")
    else
        print("black, ")
    end
    show(node.data)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export get_Jacobian

"""
$(SIGNATURES)

Obtain the symbolic Jacobian matrix of `eom` (either a `HarmonicEquation` or a `DifferentialEquation`).
This is the linearised left-hand side of F(u) = du/dT.

"""
function get_Jacobian(eom::HarmonicEquation)

    rearr = !HarmonicBalance.is_rearranged(eom) ? HarmonicBalance.rearrange_standard(eom) : eom
    lhs = _remove_brackets(rearr)
    vars = _remove_brackets.(eom.variables)

    get_Jacobian(lhs, vars)

end

" Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`. "
function get_Jacobian(diff_eom::DifferentialEquation)
    @variables T
    harmonic_eq = get_harmonic_equations(diff_eom, slow_time=T, fast_time=first(get_independent_variables(diff_eom)))
    get_Jacobian(harmonic_eq)
end

" Get the Jacobian of a set of equations `eqs` with respect to the variables `vars`. "
function get_Jacobian(eqs::Vector{Num}, vars::Vector{Num})
    length(eqs) == length(vars) || error("Jacobians are only defined for square systems!")
    M = Matrix{Num}(undef, length(vars), length(vars))

    for idx in CartesianIndices(M)
        M[idx] = expand_derivatives(d(eqs[idx[1]], vars[idx[2]]))
    end
    M
end

get_Jacobian(eqs::Vector{Equation}, vars::Vector{Num}) = get_Jacobian(Num.(getfield.(eqs, :lhs) .- getfield.(eqs, :rhs)), vars)

# the Jacobian of some equations again
function get_Jacobian_steady(eom::HarmonicEquation; differential_order=0)
    Hopf_vars = first.(getfield.(filter(x -> x.type == "Hopf", eom.variables), :symbol))
    Hopf_idx = findall(x -> any(isequal.(x, Hopf_vars)) , get_variables(eom))
    nonsingular = filter( x -> x ∉ Hopf_idx, 1:length(get_variables(eom)))

    vars_simp = Dict([var => HarmonicBalance.declare_variable(var_name(var)) for var in get_variables(eom)])
    T = get_independent_variables(eom)[1]
    J = get_Jacobian(eom.equations, d(get_variables(eom), T, differential_order))[nonsingular, nonsingular]

    expand_derivatives.(HarmonicBalance.substitute_all(J, vars_simp))
end

# COMPILE THIS!
function get_implicit_Jacobian(eom::HarmonicEquation)::Function
    M = get_Jacobian_steady(eom, differential_order=0)
    Mp = get_Jacobian_steady(eom, differential_order=1)

    function J(soln::OrderedDict)
        -inv(ComplexF64.(substitute_all(Mp, soln))) * ComplexF64.(substitute_all(M, soln))
    end
    J
end

###
# STUFF BELOW IS FOR THE CORRECTED JACOBIAN METHOD
###

"""
    get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)

Obtain the symbolic linear response matrix of a `diff_eq` corresponding to a perturbation frequency `freq`.
This routine cannot accept a `HarmonicEquation` since there, some time-derivatives are already dropped.
`order` denotes the highest differential order to be considered.

"""
function get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)::Matrix
    @variables T, i
    time = get_independent_variables(diff_eq)[1]

    eom = HarmonicBalance.harmonic_ansatz(diff_eq, time)

    # replace the time-dependence of harmonic variables by slow time BUT do not drop any derivatives
    eom = HarmonicBalance.slow_flow(eom, fast_time=time, slow_time=T, degree=order+1)

    eom = HarmonicBalance.fourier_transform(eom, time)

    M = get_Jacobian(eom.equations, get_variables(eom))
    for n in 1:order
        M += (i * freq)^n * get_Jacobian(eom.equations, d(get_variables(eom), T, n))
    end
    M = substitute_all(M, [var => HarmonicBalance.declare_variable(var_name(var)) for var in get_variables(eom)])
    substitute_all(expand_derivatives.(M), i => im)
end


text:
# FWI on the 2D Overthrust model using spectral projected gradient descent
# Author: mlouboutin3@gatech.edu
# Date: February 2021
#

using DrWatson
@quickactivate :TimeProbeSeismic
import TimeProbeSeismic: qr_data

# Load starting model
# Load starting model
~isfile(datadir("models", "overthrust_model.h5")) && run(`curl -L ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_model_2D.h5 --create-dirs -o $(datadir("models", "overthrust_model.h5"))`)
n, d, o, m0, m = h5read(datadir("models", "overthrust_model.h5"), "n", "d", "o", "m0", "m")

n = Tuple(n)
o = Tuple(o)
d = Tuple(d)

# Setup info and model structure
nsrc = 1	# number of sources
model = Model(n, d, o, m)

# Set up receiver geometry
nxrec = n[1]
xrec = range(0f0, stop=(n[1] - 1)*d[1], length=nxrec)
yrec = 0f0
zrec = range(12.5f0, stop=12.5f0, length=nxrec)

# receiver sampling and recording time
timeD = 3000f0   # receiver recording time [ms]
dtD = 4f0  # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtD, t=timeD, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
xsrc = (n[1] - 1)*d[1]/2
ysrc = 0f0
zsrc = 12.5f0

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtD, t=timeD)

# setup wavelet
f0 = 0.010f0     # kHz
wavelet = ricker_wavelet(timeD, dtD, f0)
q = judiVector(srcGeometry, wavelet)
q_dist = generate_distribution(q)

###################################################################################################
# Write shots as segy files to disk
opt = Options(space_order=16)

# Setup operators
F = judiModeling(model, srcGeometry, recGeometry; options=opt)

# data
# Nonlinear modeling
dobs = F*q

# Probing vectors
figure();imshow(dobs.data[1]*dobs.data[1]', cmap="seismic", vmin=-10, vmax=10)
title(L"$d_{obs} d_{obs}^\top$")

timea = 0:4f0:3000f0

fig, axs = subplots(5, 3, sharex=:col, sharey=:row, gridspec_kw=Dict(:hspace=> 0, :wspace=> 0))
for (i, ps)=enumerate([4, 16, 32, 64, 256])
    @show i
    # QR
    Q = qr_data(dobs.data[1], ps)
    axs[i, 1].imshow(Q*Q',vmin=-.05, vmax=.05, cmap="seismic", aspect="auto")
    i == 1 && axs[i, 1].set_title(L"$\mathbf{Q}\mathbf{Q}^\top$")
    axs[i, 1].set_ylabel("r=$(ps)", rotation=0, labelpad=20)
    # gaussian
    E = rand([-1f0, 1f0], length(timea), ps)/sqrt(ps)
    axs[i, 2].imshow(E*E',vmin=-.25, vmax=.25, cmap="seismic", aspect="auto")
    i == 1 && axs[i, 2].set_title(L"$\mathbf{Z}\mathbf{Z}^\top$")
    # DFT
    freq = select_frequencies(q_dist; fmin=0.002, fmax=0.04, nf=ps)
    F = exp.(-2*im*pi*timea*freq')/sqrt(ps)
    axs[i, 3].imshow(real.(F*F'),vmin=-.5, vmax=.5, cmap="seismic", aspect="auto")
    i == 1 && axs[i, 3].set_title(L"$\mathbf{F}\mathbf{F}^\top$")
end
plt.setp(plt.gcf().get_axes(), xticks=[], yticks=[])

wsave(plotsdir("overthrust_single", "Zi.png"), gcf())

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
if abspath(PROGRAM_FILE) == @__FILE__
    @assert(b"foo" != b"bar")
    @assert(b"\"" == b"\"")
    @assert(b"'" == b"'")
    @assert(b"\xbbfoo" == b"\xbbfoo")
    println("OK")
end


text:
α = 0.025

# Hepers

@testset "correlations" begin
  @testset "copula functions" begin
    @test Ccl([1., 0.5], 3.) == 0.5
    @test Cg([1., 0.5], 3.) == 0.5
  end
  @testset "Kendall cor" begin
    @test Debye(0.5, 1) ≈ 0.8819271567906056
    @test τ2θ(0.6, "frank") ≈ 7.929642284264058
    @test τ2θ(0.5, "gumbel") ≈ 2.
    @test τ2θ(1/3, "clayton") ≈ 1.
    @test τ2θ(1/4, "amh") ≈ 0.8384520912688538
    @test_throws AssertionError τ2θ(1/3, "frak")
    @test_throws DomainError τ2θ(0., "frank")
    @test frankτ2θ(0.6) ≈ 7.929642284264058
    @test AMHτ2θ(1/4) ≈ 0.8384520912688538
    @test AMHτ2θ(0.) ≈ 0. atol=1.0e-6
    @test AMHτ2θ(0.28) ≈ 0.9999
    @test AMHτ2θ(-2/11) ≈ -0.9999
  end
  @testset "Spearman cor" begin
    @test dilog(0.5) ≈ 0.5822405264650125
    @test claytonρ2θ(1/3) ≈ 0.58754 atol=1.0e-5
    @test claytonρ2θ(0.01) ≈ 0. atol=1.0e-1
    @test gumbelρ2θ(0.5) ≈ 1.5410704204332681
    @test gumbelρ2θ(0.0001) == 1.
    @test ρ2θ(1/3, "frank") ≈ 2.116497 atol=1.0e-5
    @test ρ2θ(0.2, "amh") ≈ 0.5168580913147318
    @test_throws AssertionError ρ2θ(1/3, "frak")
    @test_throws DomainError ρ2θ(0., "frank")
    @test frankρ2θ(1/3) ≈ 2.116497 atol=1.0e-5
    @test AMHρ2θ(0.2) ≈ 0.5168580913147318
    @test claytonθ2ρ(3.) ≈ 0.78645 atol=1.0e-4
    @test claytonθ2ρ(0.001) ≈ 0. atol=1.0e-2
    @test gumbelθ2ρ(3.) ≈ 0.8489 atol=1.0e-4
    @test AMHρ2θ(0.) ≈ 0. atol=1.0e-4
    @test AMHρ2θ(0.49) ≈ 1 atol=1.0e-4
    @test AMHρ2θ(-0.273) ≈ -1 atol=1.0e-4
  end
  @testset "negative cor" begin
    @test ρ2θ(-0.2246, "amh") ≈ -0.8 atol=1.0e-3
    @test τ2θ(-0.1505, "amh") ≈ -0.8 atol=1.0e-3
    @test ρ2θ(-0.5572, "frank") ≈ -4. atol=1.0e-3
    @test τ2θ(-0.3881, "frank") ≈ -4. atol=1.0e-3
    @test τ2θ(-1/3, "clayton") ≈ -.5 atol=1.0e-5
    @test ρ2θ(-0.4668, "clayton") ≈ -.5 atol=1.0e-3
  end
end
@testset "Archimedean copulas axiliary functions" begin
  c = arch_gen("clayton", [0.2 0.4 0.8; 0.2 0.8 0.6; 0.3 0.9 0.6], 1.; rng = Random.GLOBAL_RNG)
  @test c ≈ [0.5 0.637217; 0.362783 0.804163; 0.432159 0.896872] atol=1.0e-5
  @test useτ(0.5, "clayton") == 2.
  @test useρ(0.75, "gumbel") ≈ 2.285220798876495
  @test getθ4arch(0.5, "gumbel", SpearmanCorrelation) ≈ 1.541070420842913
  @test getθ4arch(0.5, "gumbel", KendallCorrelation) ≈ 2.0
end

@testset "Axiliary function exceptions" begin
  @test_throws DomainError testθ(0.5, "gumbel")
  @test_throws DomainError useρ(0.6, "amh")
  @test_throws DomainError useτ(0.45, "amh")
end

@testset "Gumbel copula" begin
  @testset "exceptions" begin
    @test_throws DomainError Gumbel_cop(3, 0.3)
    @test_throws DomainError Gumbel_cop(3, 1.1, KendallCorrelation)
    @test_throws DomainError Gumbel_cop(3, 1.1, SpearmanCorrelation)
    @test_throws DomainError Gumbel_cop(1, 1.1)
    @test_throws DomainError Gumbel_cop(1, 1.1, KendallCorrelation)

    @test_throws DomainError Gumbel_cop_rev(3, 0.3)
    @test_throws DomainError Gumbel_cop_rev(3, 1.1, KendallCorrelation)
    @test_throws DomainError Gumbel_cop_rev(3, 1.1, SpearmanCorrelation)
    @test_throws DomainError Gumbel_cop_rev(1, 1.1)
    @test_throws DomainError Gumbel_cop_rev(1, 1.1, KendallCorrelation)

    u = zeros(3,3)
    @test_throws AssertionError simulate_copula!(u, Gumbel_cop(2, 2.))
    @test_throws AssertionError simulate_copula!(u, Gumbel_cop_rev(2, 2.))
  end

  @testset "small example" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Gumbel_cop(2, 2.); rng = rng) ≈ [0.76229 0.812431] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Gumbel_cop_rev(2, 2.); rng = rng) ≈ [0.237703 0.1875681] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    simulate_copula!(u, Gumbel_cop(2, 2.); rng = rng)
    @test u ≈ [0.76229 0.812431] atol=1.0e-5

    Random.seed!(43)
    u = zeros(1,2)
    rng = StableRNG(123)
    simulate_copula!(u, Gumbel_cop_rev(2, 2.); rng = rng)
    @test u ≈ [0.237703 0.1875681] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Gumbel_cop(2, 2.); rng = rng) ≈ [0.76229 0.81243] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    simulate_copula!(u, Gumbel_cop(2, 2.); rng = rng)
    @test u ≈ [0.7622966 0.8124318] atol=1.0e-5
  end

  @testset "tests on larger data" begin
    Random.seed!(1234)
    x = simulate_copula(100_000, Gumbel_cop(3, 2.))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "r") ≈ 2-2^(1/2) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r") ≈ 2-2^(1/2) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0.
    @test tail(x[:,1], x[:,3], "l", 0.00001) ≈ 0.
    @test corkendall(x) ≈ [1. 1/2 1/2; 1/2 1. 1/2; 1/2 1/2 1.] atol=1.0e-2

    Random.seed!(43)
    x = simulate_copula(50000, Gumbel_cop_rev(2, 1.5))
    @test tail(x[:,1], x[:,2], "l") ≈ 2-2^(1/1.5) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r", 0.00001) ≈ 0.
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    Random.seed!(43)
    x = simulate_copula(50000, Gumbel_cop(2, 0.5, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.5; 0.5 1.] atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(50000, Gumbel_cop(2, 0.5, SpearmanCorrelation))
    @test corspearman(x) ≈ [1. 0.5; 0.5 1.] atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(50000, Gumbel_cop_rev(2, 0.5, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.5; 0.5 1.] atol=1.0e-2
  end
end
@testset "Clayton copula" begin
  @testset "exceptions" begin
    @test_throws DomainError Clayton_cop(3, -0.3)
    @test_throws DomainError Clayton_cop(2, -2.3)
    @test_throws DomainError Clayton_cop(3, 1.1, KendallCorrelation)
    @test_throws DomainError Clayton_cop(3, 1.1, SpearmanCorrelation)
    @test_throws DomainError Clayton_cop(1, 1.1)
    @test_throws DomainError Clayton_cop(1, 1.1, KendallCorrelation)

    @test_throws DomainError Clayton_cop_rev(3, -0.3)
    @test_throws DomainError Clayton_cop_rev(2, -2.3)
    @test_throws DomainError Clayton_cop_rev(3, 1.1, KendallCorrelation)
    @test_throws DomainError Clayton_cop_rev(3, 1.1, SpearmanCorrelation)
    @test_throws DomainError Clayton_cop_rev(1, 1.1)
    @test_throws DomainError Clayton_cop_rev(1, 1.1, KendallCorrelation)

    u = zeros(3,3)
    @test_throws AssertionError simulate_copula!(u, Clayton_cop(2, 2.))
    @test_throws AssertionError simulate_copula!(u, Clayton_cop_rev(2, 2.))
  end
  @testset "small example" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Clayton_cop(2, 2.); rng = rng) ≈ [0.46542 0.56630] atol=1.0e-5
    @test simulate_copula(1, Clayton_cop(2, -0.5); rng = rng) ≈ [0.04273 0.78114] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Clayton_cop_rev(2, 2.); rng = rng) ≈ 1. .- [0.46542 0.56630]  atol=1.0e-5
    @test simulate_copula(1, Clayton_cop_rev(2, -.5); rng = rng) ≈ 1. .- [0.04273 0.78114] atol=1.0e-5

    Random.seed!(43)
    u = zeros(1,2)
    u1 = zeros(1,2)
    rng = StableRNG(123)
    simulate_copula!(u, Clayton_cop(2, 2.); rng = rng)
    simulate_copula!(u1, Clayton_cop(2, -0.5); rng = rng)
    @test u ≈ [0.46542 0.56630] atol=1.0e-5
    @test u1 ≈ [0.04273 0.78114]  atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    u1 = zeros(1,2)
    simulate_copula!(u, Clayton_cop_rev(2, 2.); rng = rng)
    simulate_copula!(u1, Clayton_cop_rev(2, -0.5); rng = rng)
    @test u ≈ 1 .- [0.46542 0.56630] atol=1.0e-5
    @test u1 ≈ 1 .- [0.04273 0.78114] atol=1.0e-5
  end
  @testset "test on larger data" begin
    Random.seed!(43)
    x = simulate_copula(50000, Clayton_cop(3, 1.))

    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l") ≈ 2.0^(-1) atol=1.0e-1
    @test tail(x[:,1], x[:,3], "l") ≈ 2.0^(-1) atol=1.0e-1
    @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
    @test corkendall(x) ≈ [1. 1/3 1/3; 1/3 1. 1/3; 1/3 1/3 1.] atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(50000, Clayton_cop(2, 0.5, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.5; 0.5 1.] atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(50000, Clayton_cop(2, -0.9))
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test corkendall(x)[1,2] ≈ -0.9/(2-0.9) atol=1.0e-2

    Random.seed!(43)
    x = simulate_copula(50000, Clayton_cop_rev(2, 0.5, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.5; 0.5 1.] atol=1.0e-2

    Random.seed!(43)
    x = simulate_copula(50000, Clayton_cop_rev(2, -0.9))
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test corkendall(x)[1,2] ≈ -0.9/(2-0.9) atol=1.0e-2

  end
end

@testset "Ali-Mikhail-Haq copula" begin
  @testset "exceptions" begin
    @test_throws DomainError AMH_cop(3, -0.3)
    @test_throws DomainError AMH_cop(2, -1.2)
    @test_throws DomainError AMH_cop(3, 0.5, KendallCorrelation)
    @test_throws DomainError AMH_cop(3, 0.6, SpearmanCorrelation)
    @test_throws DomainError AMH_cop(1, 1.1)
    @test_throws DomainError AMH_cop(1, 1.1, KendallCorrelation)

    @test_throws DomainError AMH_cop_rev(3, -0.3)
    @test_throws DomainError AMH_cop_rev(2, -1.2)
    @test_throws DomainError AMH_cop_rev(3, 0.5, KendallCorrelation)
    @test_throws DomainError AMH_cop_rev(3, 0.55, SpearmanCorrelation)
    @test_throws DomainError AMH_cop_rev(1, 1.1)
    @test_throws DomainError AMH_cop_rev(1, 1.1, KendallCorrelation)

    u = zeros(3,3)
    @test_throws AssertionError simulate_copula!(u, AMH_cop(2, .5))
    @test_throws AssertionError simulate_copula!(u, AMH_cop_rev(2, .5))
  end
  @testset "small example" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, AMH_cop(2, 0.5); rng = rng) ≈ [0.270221 0.434935] atol=1.0e-5
    @test simulate_copula(1, AMH_cop(2, -0.5); rng = rng) ≈ [0.042730 0.5317979] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, AMH_cop_rev(2, 0.5); rng = rng) ≈ 1 .- [0.270221 0.434935] atol=1.0e-5
    @test simulate_copula(1, AMH_cop_rev(2, -0.5); rng = rng) ≈ 1 .- [0.042730 0.5317979] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    u1 = zeros(1,2)
    simulate_copula!(u, AMH_cop(2, 0.5); rng = rng)
    simulate_copula!(u1, AMH_cop(2, -0.5); rng = rng)
    @test u ≈ [0.270221 0.434935] atol=1.0e-5
    @test u1 ≈ [0.042730 0.5317979] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    u1 = zeros(1,2)
    simulate_copula!(u, AMH_cop_rev(2, 0.5); rng = rng)
    simulate_copula!(u1, AMH_cop_rev(2, -0.5); rng = rng)
    @test u ≈ 1 .- [0.270221 0.434935] atol=1.0e-5
    @test u1 ≈ 1 .- [0.042730 0.5317979] atol=1.0e-5
  end
  @testset "test on larger data" begin

    Random.seed!(43)
    x = simulate_copula(100_000, AMH_cop(3, 0.8))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
    @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0
    @test corkendall(x)[1:2, 1:2] ≈ [1. 0.233; 0.233 1.] atol=1.0e-2
    x = simulate_copula(100_000, AMH_cop(2, 0.25, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.25; 0.25 1.] atol=1.0e-3

    Random.seed!(43)
    x = simulate_copula(100000, AMH_cop_rev(3, 0.8))

    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
    @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0

    Random.seed!(43)
    x = simulate_copula(100000, AMH_cop_rev(2, -0.4))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
    @test tail(x[:,1], x[:,2], "r", 0.0001) ≈ 0

    Random.seed!(43)
    x = simulate_copula(100000, AMH_cop_rev(2, 0.2, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.2; 0.2 1.] atol=1.0e-2
  end
end

@testset "Frank copula" begin
  @testset "exceptions" begin
    @test_throws DomainError Frank_cop(3, -0.3)
    @test_throws DomainError Frank_cop(2, 0.)
    @test_throws DomainError Frank_cop(3, 1.1, KendallCorrelation)
    @test_throws DomainError Frank_cop(3, 1.1, SpearmanCorrelation)
    @test_throws DomainError Frank_cop(1, 1.1)
    @test_throws DomainError Frank_cop(1, 1.1, KendallCorrelation)

    u = zeros(3,3)
    @test_throws AssertionError simulate_copula!(u, Frank_cop(2, .5))
  end
  @testset "small example" begin
    Random.seed!(43)
    rng = StableRNG(123)
    @test simulate_copula(1, Frank_cop(2, 2.); rng = rng) ≈ [0.33582 0.48292] atol=1.0e-5
    @test simulate_copula(1, Frank_cop(2, -2.); rng = rng) ≈ [0.04273 0.64937] atol=1.0e-5

    Random.seed!(43)
    rng = StableRNG(123)
    u = zeros(1,2)
    u1 = zeros(1,2)
    simulate_copula!(u, Frank_cop(2, 2.); rng= rng)
    simulate_copula!(u1, Frank_cop(2, -2.); rng = rng)

    @test u ≈ [0.33582 0.48292] atol=1.0e-5
    @test u1 ≈ [0.04273 0.64937] atol=1.0e-5

  end
  @testset "test on larger data" begin
    Random.seed!(43)
    x = simulate_copula(100000, Frank_cop(3, 0.8))

    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
    @test tail(x[:,2], x[:,3], "r", 0.0001) ≈ 0
    Random.seed!(43)
    x = simulate_copula(100000, Frank_cop(2, 0.2, KendallCorrelation))
    @test corkendall(x) ≈ [1. 0.2; 0.2 1.] atol=1.0e-2
    Random.seed!(43)
    x = simulate_copula(100000, Frank_cop(2, -2.))
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
    @test tail(x[:,1], x[:,2], "l", 0.0001) ≈ 0
  end
end

@testset "tests on Big Float" begin
  Random.seed!(1234)
  θ = BigFloat(2.)
  x = simulate_copula(1000, Gumbel_cop(3, θ))

  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α

  Random.seed!(1234)
  θ = BigFloat(-.5)
  x = simulate_copula(1000, Clayton_cop(2, θ))

  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

  Random.seed!(1234)
  θ = BigFloat(2.5)
  if false
    x = simulate_copula(1000, Clayton_cop(3, θ))

    @test typeof(x) == Array{BigFloat,2}
    @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
    @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  end

  Random.seed!(1234)
  θ = BigFloat(2.)
  x = simulate_copula(1000, Frank_cop(2, θ))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

  Random.seed!(1234)
  θ = BigFloat(.5)
  x = simulate_copula(25, AMH_cop(3, θ))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

  Random.seed!(1234)
  θ = BigFloat(-.3)
  x = simulate_copula(10, AMH_cop(2, θ))
  @test typeof(x) == Array{BigFloat,2}
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α

end


text:
"""
    open_xlsx_template(filepath::AbstractString) :: XLSXFile

Open an Excel file as template for editing and saving to another file with `XLSX.writexlsx`.

The returned `XLSXFile` instance is in closed state.
"""
@inline open_xlsx_template(filepath::AbstractString) :: XLSXFile = open_or_read_xlsx(filepath, true, true, true)

const EMPTY_EXCEL_TEMPLATE = joinpath(@__DIR__, "..", "data", "blank.xlsx")

"""
    open_empty_template(sheetname::AbstractString="") :: XLSXFile

Returns an empty, writable `XLSXFile` with 1 worksheet.

`sheetname` is the name of the worksheet, defaults to `Sheet1`.
"""
function open_empty_template(sheetname::AbstractString="") :: XLSXFile
    @assert isfile(EMPTY_EXCEL_TEMPLATE) "Couldn't find template file $EMPTY_EXCEL_TEMPLATE."
    xf = open_xlsx_template(EMPTY_EXCEL_TEMPLATE)

    if sheetname != ""
        rename!(xf[1], sheetname)
    end

    return xf
end

"""
    writexlsx(output_filepath, xlsx_file; [overwrite=false])

Writes an Excel file given by `xlsx_file::XLSXFile` to file at path `output_filepath`.

If `overwrite=true`, `output_filepath` will be overwritten if it exists.
"""
function writexlsx(output_filepath::AbstractString, xf::XLSXFile; overwrite::Bool=false)

    @assert is_writable(xf) "XLSXFile instance is not writable."
    @assert !isopen(xf) "Can't save an open XLSXFile."
    @assert all(values(xf.files)) "Some internal files were not loaded into memory. Did you use `XLSX.open_xlsx_template` to open this file?"
    if !overwrite
        @assert !isfile(output_filepath) "Output file $output_filepath already exists."
    end

    update_worksheets_xml!(xf)

    xlsx = ZipFile.Writer(output_filepath)

    # write XML files
    for f in keys(xf.files)
        if f == "xl/sharedStrings.xml"
            # sst will be generated below
            continue
        end

        io = ZipFile.addfile(xlsx, f)
        EzXML.print(io, xf.data[f])
    end

    # write binary files
    for f in keys(xf.binary_data)
        io = ZipFile.addfile(xlsx, f)
        ZipFile.write(io, xf.binary_data[f])
    end

    if !isempty(get_sst(xf))
        io = ZipFile.addfile(xlsx, "xl/sharedStrings.xml")
        print(io, generate_sst_xml_string(get_sst(xf)))
    end

    close(xlsx)

    # fix libuv issue on windows (#42)
<<<<<<< HEAD
    Sys.iswindows ? gc() : nothing
=======
    @static Sys.iswindows() ? GC.gc() : nothing
>>>>>>> upstream/master
end

get_worksheet_internal_file(ws::Worksheet) = "xl/" * get_relationship_target_by_id(get_workbook(ws), ws.relationship_id)
get_worksheet_xml_document(ws::Worksheet) = get_xlsxfile(ws).data[ get_worksheet_internal_file(ws) ]

function set_worksheet_xml_document!(ws::Worksheet, xdoc::EzXML.Document)
    xf = get_xlsxfile(ws)
    filename = get_worksheet_internal_file(ws)
    @assert haskey(xf.data, filename) "Internal file not found for $(ws.name)."
    xf.data[filename] = xdoc
end

function generate_sst_xml_string(sst::SharedStrings) :: String
    @assert sst.is_loaded "Can't generate XML string from a Shared String Table that is not loaded."
    buff = IOBuffer()

    # TODO: <sst count="89"
    print(buff, """
<?xml version="1.0" encoding="UTF-8" standalone="yes"?><sst uniqueCount="$(length(sst))" xmlns="http://schemas.openxmlformats.org/spreadsheetml/2006/main">
""")

    for s in sst.formatted_strings
        print(buff, s)
    end

    print(buff, "</sst>")
    return String(take!(buff))
end

function update_worksheets_xml!(xl::XLSXFile)
    buff = IOBuffer()

    wb = get_workbook(xl)
    for i in 1:sheetcount(wb)
        sheet = getsheet(wb, i)
        doc = get_worksheet_xml_document(sheet)
        xroot = EzXML.root(doc)

        # check namespace and root node name
        @assert get_default_namespace(xroot) == SPREADSHEET_NAMESPACE_XPATH_ARG[1][2] "Unsupported Spreadsheet XML namespace $(get_default_namespace(xroot))."
        @assert EzXML.nodename(xroot) == "worksheet" "Malformed Excel file. Expected root node named `worksheet` in worksheet XML file."

        # forces a document copy to avoid crash: munmap_chunk(): invalid pointer
        EzXML.print(buff, doc)
        doc_copy = EzXML.parsexml(String(take!(buff)))

        # deletes all elements under sheetData
        child_nodes = EzXML.findall("/xpath:worksheet/xpath:sheetData/xpath:row", EzXML.root(doc_copy), SPREADSHEET_NAMESPACE_XPATH_ARG)
        for c in child_nodes
            EzXML.unlink!(c)
        end
        c = nothing
        child_nodes = nothing

        # updates sheetData
        sheetData_node = EzXML.findfirst("/xpath:worksheet/xpath:sheetData", EzXML.root(doc_copy), SPREADSHEET_NAMESPACE_XPATH_ARG)
        spans_str = string(column_number(get_dimension(sheet).start), ":", column_number(get_dimension(sheet).stop))

        # iterates over WorksheetCache cells and write the XML
        for r in eachrow(sheet)
            ordered_column_indexes = sort(collect(keys(r.rowcells)))

            row_node = EzXML.addelement!(sheetData_node, "row")
            row_node["r"] = string(row_number(r))
            row_node["spans"] = spans_str

            # add cells to row
            for c in ordered_column_indexes
                cell = getcell(r, c)
                c_element = EzXML.addelement!(row_node, "c")

                c_element["r"] = cell.ref.name

                if cell.datatype != ""
                    c_element["t"] = cell.datatype
                end

                if cell.style != ""
                    c_element["s"] = cell.style
                end

                if cell.formula != ""
                    f_node = EzXML.addelement!(c_element, "f")
                    EzXML.setnodecontent!(f_node, cell.formula)
                end

                if cell.value != ""
                    v_node = EzXML.addelement!(c_element, "v")
                    EzXML.setnodecontent!(v_node, cell.value)
                end
            end
        end

        # updates worksheet dimension
        dimension_node = EzXML.findfirst("/xpath:worksheet/xpath:dimension", EzXML.root(doc_copy), SPREADSHEET_NAMESPACE_XPATH_ARG)
        dimension_node["ref"] = string(get_dimension(sheet))

        set_worksheet_xml_document!(sheet, doc_copy)
    end

    nothing
end

function setdata!(ws::Worksheet, cell::Cell)
    @assert is_writable(get_xlsxfile(ws)) "XLSXFile instance is not writable."
    @assert ws.cache != nothing "Can't write data to a Worksheet with empty cache."
    cache = ws.cache

    r = row_number(cell)
    c = column_number(cell)

    if !haskey(cache.cells, r)
        push!(cache.rows_in_cache, r)
        sort!(cache.rows_in_cache)
        cache.cells[r] = Dict{Int, Cell}()

        for i in 1:length(cache.rows_in_cache)
            if cache.rows_in_cache[i] == r
                cache.row_index[r] = i
            end
        end
    end
    cache.cells[r][c] = cell

    # update worksheet dimension
    ws_dimension = get_dimension(ws)

    top = row_number(ws_dimension.start)
    left = column_number(ws_dimension.start)

    bottom = row_number(ws_dimension.stop)
    right = column_number(ws_dimension.stop)

    if r < top || c < left
        top = min(r, top)
        left = min(c, left)
        set_dimension!(ws, CellRange(top, left, bottom, right))
    elseif r > bottom || c > right
        bottom = max(r, bottom)
        right = max(c, right)
        set_dimension!(ws, CellRange(top, left, bottom, right))
    end

    nothing
end

function xlsx_escape(str::AbstractString)
    if isempty(str)
        return str
    end

    buffer = IOBuffer()

    for c in str
        if c == '&'
            write(buffer, "&")
        elseif c == '"'
            write(buffer, """)
        elseif c == '<'
            write(buffer, "<")
        elseif c == '>'
            write(buffer, ">")
        elseif c == '\''
            write(buffer, "'")
        else
            write(buffer, c)
        end
    end

    return String(take!(buffer))
end

"""
Returns the datatype and value for `val` to be inserted into `ws`.
"""
function xlsx_encode(ws::Worksheet, val::AbstractString)
    if isempty(val)
        return ("", "")
    end
    sst_ind = add_shared_string!(get_workbook(ws), xlsx_escape(val))
    return ("s", string(sst_ind))
end

xlsx_encode(::Worksheet, val::Missing) = ("", "")
xlsx_encode(::Worksheet, val::Bool) = ("b", val ? "1" : "0")
xlsx_encode(::Worksheet, val::Union{Int, Float64}) = ("", string(val))
xlsx_encode(ws::Worksheet, val::Dates.Date) = ("", string(date_to_excel_value(val, isdate1904(get_xlsxfile(ws)))))
xlsx_encode(ws::Worksheet, val::Dates.DateTime) = ("", string(datetime_to_excel_value(val, isdate1904(get_xlsxfile(ws)))))
xlsx_encode(::Worksheet, val::Dates.Time) = ("", string(time_to_excel_value(val)))

function setdata!(ws::Worksheet, ref::CellRef, val::CellValue)
    t, v = xlsx_encode(ws, val.value)
    cell = Cell(ref, t, id(val.styleid), v, "")

    setdata!(ws, cell)
end

setdata!(ws::Worksheet, ref::CellRef, val::CellValueType) = setdata!(ws, ref, CellValue(ws, val))
setdata!(ws::Worksheet, ref_str::AbstractString, value) = setdata!(ws, CellRef(ref_str), value)

setdata!(ws::Worksheet, ref::CellRef, value) = error("Unsupported datatype $(typeof(value)) for writing data to Excel file. Supported data types are $(CellValueType) or $(CellValue).")

Base.setindex!(ws::Worksheet, v, ref) = setdata!(ws, ref, v)

function writetable!(sheet::Worksheet, data, columnnames; anchor_cell::CellRef=CellRef("A1"))

    # read dimensions
    col_count = length(data)
    @assert col_count == length(columnnames) "Column count moccursin between `data` ($col_count columns) and `columnnames` ($(length(columnnames)) columns)."
    @assert col_count > 0 "Can't write table with no columns."
    row_count = length(data[1])
    if col_count > 1
        for c in 2:col_count
            @assert length(data[c]) == row_count "Row count moccursin between column 1 ($row_count rows) and column $c ($(length(data[c])) rows)."
        end
    end

    anchor_row = row_number(anchor_cell)
    anchor_col = column_number(anchor_cell)

    # write table header
    for c in 1:col_count
        target_cell_ref = CellRef(anchor_row, c + anchor_col - 1)
        sheet[target_cell_ref] = string(columnnames[c])
    end

    # write table data
    for r in 1:row_count, c in 1:col_count
        target_cell_ref = CellRef(r + anchor_row, c + anchor_col - 1)
        sheet[target_cell_ref] = data[c][r]
    end
end

function Base.setindex!(sheet::Worksheet, value, row::Integer, col::Integer)
    target_cell_ref = CellRef(row, col)
    sheet[target_cell_ref] = value
end

function Base.setindex!(sheet::Worksheet, data::AbstractVector, row::Integer, cols::UnitRange{<:Integer})
    col_count = length(data)

    @assert col_count == length(cols) "Column count moccursin between `data` ($col_count columns) and column range $cols ($(length(cols)) columns)."

    for c in 1:col_count
        target_cell_ref = CellRef(row, c + first(cols) - 1)
        sheet[target_cell_ref] = data[c]
    end
end

function Base.setindex!(sheet::Worksheet, data::AbstractVector, row::Integer, c::Colon)
    col_count = length(data)

    for c in 1:col_count
        target_cell_ref = CellRef(row, c)
        sheet[target_cell_ref] = data[c]
    end
end

Base.setindex!(sheet::Worksheet, data::AbstractVector, ref_str::AbstractString) = setindex!(sheet, data, CellRef(ref_str))

function Base.setindex!(sheet::Worksheet, data::AbstractVector, index::CellRef)
    col_count = length(data)
    anchor_row = row_number(index)
    anchor_col = column_number(index)

    for c in 1:col_count
        target_cell_ref = CellRef(anchor_row, c + anchor_col - 1)
        sheet[target_cell_ref] = data[c]
    end
end

function rename!(ws::Worksheet, name::AbstractString)
    xf = get_xlsxfile(ws)
    @assert is_writable(xf) "XLSXFile instance is not writable."

    # updates XML
    xroot = xmlroot(xf, "xl/workbook.xml")
    for node in EzXML.eachelement(xroot)
        if EzXML.nodename(node) == "sheets"

            for sheet_node in EzXML.eachelement(node)
                if sheet_node["name"] == ws.name
                    # assign new name
                    sheet_node["name"] = name
                    break
                end
            end

            break
        end
    end

    # updates the new name in the worksheet instance
    ws.name = name
    nothing
end

const FILEPATH_SHEET_TEMPLATE = joinpath(@__DIR__, "..", "data", "sheet_template.xml")

addsheet!(xl::XLSXFile, name::AbstractString="") :: Worksheet = addsheet!(get_workbook(xl), name)

"""
    addsheet!(workbook, [name]) :: Worksheet

Create a new worksheet with named `name`.
If `name` is not provided, a unique name is created.

"""
function addsheet!(wb::Workbook, name::AbstractString="") :: Worksheet

    xf = get_xlsxfile(wb)
    @assert is_writable(xf) "XLSXFile instance is not writable."

    @assert isfile(FILEPATH_SHEET_TEMPLATE) "Couldn't find template file $FILEPATH_SHEET_TEMPLATE."

    if name == ""
        # name was not provided. Will find a unique name.
        i = 1
        current_sheet_names = sheetnames(wb)
        while true
            name = "Sheet$i"
            if !in(name, current_sheet_names)
                # found a unique name
                break
            end
            i += 1
        end
    end

    @assert name != ""

    # generate sheetId
    current_sheet_ids = [ ws.sheetId for ws in wb.sheets ]
    sheetId = max(current_sheet_ids...) + 1

    xdoc = EzXML.readxml(FILEPATH_SHEET_TEMPLATE)

    # generate a unique name for the XML
    local xml_filename::String
    i = 1
    while true
        xml_filename = "xl/worksheets/sheet$i.xml"
        if !in(xml_filename, keys(xf.files))
            break
        end
        i += 1
    end

    # adds doc do XLSXFile
    xf.files[xml_filename] = true # is read
    xf.data[xml_filename] = xdoc

    # adds workbook-level relationship
    # <Relationship Id="rId1" Target="worksheets/sheet1.xml" Type="http://schemas.openxmlformats.org/officeDocument/2006/relationships/worksheet"/>
    rId = add_relationship!(wb, xml_filename[4:end], "http://schemas.openxmlformats.org/officeDocument/2006/relationships/worksheet")

    # creates Worksheet instance
    ws = Worksheet(xf, sheetId, rId, name, CellRange("A1:A1"))

    # creates a mock WorksheetCache
    # because we can't write to sheet with empty cache (see setdata!(ws::Worksheet, cell::Cell))
    # and the stream should be closed
    # to indicate that no more rows will be fetched from SheetRowStreamIterator in Base.iterate(ws_cache::WorksheetCache, row_from_last_iteration::Int)
    itr = SheetRowStreamIterator(ws)
    zip_io, reader = open_internal_file_stream(xf, "[Content_Types].xml") # could be any file
    state = SheetRowStreamIteratorState(zip_io, reader, true, 0)
    close(state)
    ws.cache = WorksheetCache(CellCache(), Vector{Int}(), Dict{Int, Int}(), itr, state)

    # adds the new sheet to the list of sheets in the workbook
    push!(wb.sheets, ws)

    # updates workbook xml
    xroot = xmlroot(xf, "xl/workbook.xml")
    for node in EzXML.eachelement(xroot)
        if EzXML.nodename(node) == "sheets"

            #<sheet name="Sheet1" r:id="rId1" sheetId="1"/>
            sheet_element = EzXML.addelement!(node, "sheet")
            sheet_element["name"] = name
            sheet_element["r:id"] = rId
            sheet_element["sheetId"] = string(sheetId)

            break
        end
    end

    return ws
end

#
# Helper Functions
#

"""
    writetable(filename, data, columnnames; [overwrite], [sheetname])

`data` is a vector of columns.
`columnames` is a vector of column labels.
`overwrite` is a `Bool` to control if `filename` should be overwritten if already exists.
`sheetname` is the name for the worksheet.

Example using `DataFrames.jl`:

```julia
import DataFrames, XLSX
df = DataFrames.DataFrame(integers=[1, 2, 3, 4], strings=["Hey", "You", "Out", "There"], floats=[10.2, 20.3, 30.4, 40.5])
XLSX.writetable("df.xlsx", DataFrames.columns(df), DataFrames.names(df))
```
"""
function writetable(filename::AbstractString, data, columnnames; overwrite::Bool=false, sheetname::AbstractString="", anchor_cell::Union{String, CellRef}=CellRef("A1"))

    if !overwrite
        @assert !isfile(filename) "$filename already exists."
    end

    xf = open_empty_template(sheetname)
    sheet = xf[1]

    if isa(anchor_cell, String)
        anchor_cell = CellRef(anchor_cell)
    end

    writetable!(sheet, data, columnnames; anchor_cell=anchor_cell)

    # write output file
    writexlsx(filename, xf, overwrite=overwrite)
    nothing
end

"""
    writetable(filename::AbstractString; overwrite::Bool=false, kw...)
    writetable(filename::AbstractString, tables::Vector{Tuple{String, Vector{Any}, Vector{String}}}; overwrite::Bool=false)

Write multiple tables.

`kw` is a variable keyword argument list. Each element should be in this format: `sheetname=( data, column_names )`,
where `data` is a vector of columns and `column_names` is a vector of column labels.

Example:

```julia
import DataFrames, XLSX

df1 = DataFrames.DataFrame(COL1=[10,20,30], COL2=["Fist", "Sec", "Third"])
df2 = DataFrames.DataFrame(AA=["aa", "bb"], AB=[10.1, 10.2])

XLSX.writetable("report.xlsx", REPORT_A=( DataFrames.columns(df1), DataFrames.names(df1) ), REPORT_B=( DataFrames.columns(df2), DataFrames.names(df2) ))
```
"""
function writetable(filename::AbstractString; overwrite::Bool=false, kw...)

    if !overwrite
        @assert !isfile(filename) "$filename already exists."
    end

    xf = open_empty_template()
    is_first = true

    for (sheetname, (data, column_names)) in kw
        if is_first
            # first sheet already exists in template file
            sheet = xf[1]
            rename!(sheet, string(sheetname))
            writetable!(sheet, data, column_names)

            is_first = false
        else
            sheet = addsheet!(xf, string(sheetname))
            writetable!(sheet, data, column_names)
        end
    end

    # write output file
    writexlsx(filename, xf, overwrite=overwrite)
    nothing
end

function writetable(filename::AbstractString, tables::Vector{Tuple{String, Vector{Any}, Vector{T}}}; overwrite::Bool=false) where {T<:Union{String, Symbol}}

    if !overwrite
        @assert !isfile(filename) "$filename already exists."
    end

    xf = open_empty_template()

    is_first = true

    for (sheetname, data, column_names) in tables
        if is_first
            # first sheet already exists in template file
            sheet = xf[1]
            rename!(sheet, string(sheetname))
            writetable!(sheet, data, column_names)

            is_first = false
        else
            sheet = addsheet!(xf, string(sheetname))
            writetable!(sheet, data, column_names)
        end
    end

    # write output file
    writexlsx(filename, xf, overwrite=overwrite)
end


text:
"""
Query a function pointer for an API function.

One of `get_instance_proc_addr` or `get_device_proc_addr` will be used depending on the argument types.
"""
function function_pointer end

function_pointer(name::AbstractString) = get_instance_proc_addr(name)
function_pointer(instance::Instance, name::AbstractString) = get_instance_proc_addr(name; instance)
function_pointer(device::Device, name::AbstractString) = get_device_proc_addr(device, name)
function_pointer(x, name::AbstractString) = function_pointer(handle(x), name)


text:
using LightQuery: @>, @_, By, Group, over, make_columns, @name_str, Rows, value

using CSV: File
using DataFrames: DataFrame, groupby, combine
using Tables: columntable

cd("/home/brandon/benchmark")
download("http://rapidsai-data.s3-website.us-east-2.amazonaws.com/notebook-mortgage-data/mortgage_2000.tgz", "mortgage_2000.tgz")
run(`tar --gzip --extract --file=mortgage_2000.tgz`)
cd("perf")

file = File(
    "Performance_2000Q1.txt",
    delim = '|',
    header = Symbol.(string.("Column", 1:31)),
    missingstrings = ["NULL", ""],
    dateformat = "mm/dd/yyyy",
    truestrings = ["Y"],
    falsestrings = ["N"],
)

# just the type stable ones
columns = (name"Column1", name"Column2", name"Column4", name"Column6", name"Column10", name"Column11", name"Column12")(columntable(file))

function process_with_lightquery(columns)
    @> Rows(; columns...) |>
    Group(By(_, name"Column1")) |>
    Iterators.map((@_ (Count = length(value(_)),)), _) |>
    make_columns
end

@time process_with_lightquery(columns)


text:
const Maybe{T} = Union{T,Nothing}
const TupleN{T,N} = NTuple{N,T}

const AbsArr{T,N} = AbstractArray{T,N}
const AbsMat{T} = AbstractMatrix{T}
const AbsVec{T} = AbstractVector{T}

macro mustimplement(sig)
    :($(esc(sig)) = error("must implement ", $(string(sig))))
end

@noinline argerror(msg::AbstractString) = throw(ArgumentError(msg))
@noinline dimerror(msg::AbstractString) = throw(DimensionMismatch(msg))
@noinline internalerror(msg::AbstractString) = error("Internal error: $msg\nPlease file a bug.")
@noinline internalerror() = error("Internal error. Please file a bug.")