text:
function msfm2d(FSpeedImage, SourcePoints, usesecond, usecross)
# This function MSFM2D calculates the shortest distance from a list of
# points to all other pixels in an image, using the
# Multistencil Fast Marching Method (MSFM). This method gives more accurate
# distances by using second order derivatives and cross neighbours.
#
# T=msfm2d(F, SourcePoints, UseSecond, UseCross)
#
# inputs,
#   F: The speed image. The speed function must always be larger
#			than zero (min value 1e-8), otherwise some regions will
#			never be reached because the time will go to infinity.
#   SourcePoints : A list of starting points [2 x N] (distance zero)
#   UseSecond : Boolean Set to true if not only first but also second
#                order derivatives are used (default)
#   UseCross : Boolean Set to true if also cross neighbours
#                are used (default)
# outputs,
#   T : Image with distance from SourcePoints to all pixels
#
# Note:
#   Compile the c file "mex msfm2d.c"  for cpu-effective registration
#
# Literature : M. Sabry Hassouna et Al. Multistencils Fast Marching
#   Methods: A Highly Accurate Solution to the Eikonal Equation on
#   Cartesian Domains
#
# Example,
#   SourcePoint = [51; 51];
#   SpeedImage = ones([101 101]);
#   [X Y] = ndgrid(1:101, 1:101);
#   T1 = sqrt((X-SourcePoint(1)).^2 + (Y-SourcePoint(2)).^2);
#
#   # Run fast marching 1th order, 1th order multi stencil
#   # and 2th orde and 2th orde multi stencil
#
#   tic; T1_FMM1 = msfm2d(SpeedImage, SourcePoint, false, false); toc;
#   tic; T1_MSFM1 = msfm2d(SpeedImage, SourcePoint, false, true); toc;
#   tic; T1_FMM2 = msfm2d(SpeedImage, SourcePoint, true, false); toc;
#   tic; T1_MSFM2 = msfm2d(SpeedImage, SourcePoint, true, true); toc;
#
#   # Show results
#   fprintf('\nResults with T1 (Matlab)\n');
#   fprintf('Method   L1        L2        Linf\n');
#   Results = cellfun(@(x)([mean(abs(T1(:)-x(:))) mean((T1(:)-x(:)).^2) max(abs(T1(:)-x(:)))]), {T1_FMM1(:) T1_MSFM1(:) T1_FMM2(:) T1_MSFM2(:)}, 'UniformOutput',false);
#   fprintf('FMM1:   %9.5f %9.5f %9.5f\n', Results{1}(1), Results{1}(2), Results{1}(3));
#   fprintf('MSFM1:  %9.5f %9.5f %9.5f\n', Results{2}(1), Results{2}(2), Results{2}(3));
#   fprintf('FMM2:   %9.5f %9.5f %9.5f\n', Results{3}(1), Results{3}(2), Results{3}(3));
#   fprintf('MSFM2:  %9.5f %9.5f %9.5f\n', Results{4}(1), Results{4}(2), Results{4}(3));
#
# Example multiple starting points,
#   SourcePoint=rand(2,100)*255+1;
#   SpeedImage = ones([256 256]);
#   tic; T1_MSFM2 = msfm2d(SpeedImage, SourcePoint, true, true); toc;
#   figure, imshow(T1_MSFM2,[]); colormap(hot(256));
#
# Function is written by D.Kroon University of Twente (June 2009)

# Distance image, also used to store the index of narrowband pixels
# during marching process
T = zeros(size(F))-1;

# Augmented Fast Marching (For skeletonize)
Ed=nargout>1;

# Euclidian distance image
if(Ed)
	Y = zeros(size(F))
 end

# Pixels which are processed and have a final distance are frozen
Frozen   = zeros(size(F));

# Free memory to store neighbours of the (segmented) region
neg_free = 100000;
neg_pos=0;
if(Ed)
    neg_list = zeros(4,neg_free);
else
    neg_list = zeros(3,neg_free);
end

# (There are 3 pixel classes:
#   - frozen (processed)
#   - narrow band (boundary) (in list to check for the next pixel with smallest distance)
#   - far (not yet used)

# Neighbours
ne =[-1 0;
    1 0;
    0 -1;
    0 1];

SourcePoints=int32(floor(SourcePoints));

# set all starting points to distance zero and frozen
for z=1:size(SourcePoints,2)
    # starting point
    x= SourcePoints(1,z); y=SourcePoints(2,z);
    # Set starting point to frozen and distance to zero
    Frozen(x,y)=1; T(x,y)=0;
end

# Add all neighbours of the starting points to narrow list
for z=1:size(SourcePoints,2)
    # starting point
    x=SourcePoints(1,z);
	y=SourcePoints(2,z);
    for k=1:4,
        # Location of neighbour
        i=x+ne(k,1); j=y+ne(k,2);
        # Check if current neighbour is not yet frozen and inside the
        # picture
        if((i>0)&&(j>0)&&(i<=size(F,1))&&(j<=size(F,2))&&(Frozen(i,j)==0))
            Tt=1/max(F(i,j),eps);
	        Ty=1;
            # Update distance in neigbour list or add to neigbour list
            if(T(i,j)>0)
				if(neg_list(1,T(i,j))>Tt)
					neg_list(1,T(i,j))=Tt;
	            end
				if(Ed)
                    neg_list(4,T(i,j))=min(Ty,neg_list(4,T(i,j)));
                end
            else
	            neg_pos=neg_pos+1;
                # If running out of memory at a new block
                if(neg_pos>neg_free)
			 neg_free = neg_free +100000;
			 neg_list(1,neg_free)=0;
	         end
                if(Ed)
                    neg_list(:,neg_pos)=[Tt;i;j;Ty];
                else
                    neg_list(:,neg_pos)=[Tt;i;j];
                end
                T(i,j)=neg_pos;
            end
        end
    end
end

# Loop through all pixels of the image
for itt=1:numel(F)
    # Get the pixel from narrow list (boundary list) with smallest
    # distance value and set it to current pixel location
    [t,index]=min(neg_list(1,1:neg_pos));
    if(neg_pos==0)
	 break;
    end
    x=neg_list(2,index); y=neg_list(3,index);
    Frozen(x,y)=1;
    T(x,y)=neg_list(1,index);

    if(Ed)
	 Y(x,y)=neg_list(4,index);
    end

    # Remove min value by replacing it with the last value in the array
    if(index<neg_pos)
        neg_list(:,index)=neg_list(:,neg_pos);
        x2=neg_list(2,index); y2=neg_list(3,index);
        T(x2,y2)=index;
    end
    neg_pos =neg_pos-1;

    # Loop through all 4 neighbours of current pixel
    for k=1:4,
        # Location of neighbour
        i=x+ne(k,1); j=y+ne(k,2);

        # Check if current neighbour is not yet frozen and inside the
        # picture
        if((i>0)&&(j>0)&&(i<=size(F,1))&&(j<=size(F,2))&&(Frozen(i,j)==0))

            Tt=CalculateDistance(T,F(i,j),size(F),i,j,usesecond,usecross,Frozen);
            if(Ed)
                Ty=CalculateDistance(Y,1,size(F),i,j,usesecond,usecross,Frozen);
            end

            # Update distance in neigbour list or add to neigbour list
            if(T(i,j)>0)
                neg_list(1,T(i,j))=min(Tt,neg_list(1,T(i,j)));
                if(Ed)
                    neg_list(4,T(i,j))=min(Ty,neg_list(4,T(i,j)));
                end
            else
                neg_pos=neg_pos+1;
                # If running out of memory at a new block
                if(neg_pos>neg_free)
			 neg_free = neg_free +100000;
			 neg_list(1,neg_free)=0;
	        end
                if(Ed)
                    neg_list(:,neg_pos)=[Tt;i;j;Ty];
                else
                    neg_list(:,neg_pos)=[Tt;i;j];
                end
                T(i,j)=neg_pos;
            end
        end
    end
end

function CalculateDistance(T,Fij,sizeF,i,j,usesecond,usecross,Frozen)
# Boundary and frozen check -> current patch
Tpatch=inf(5,5);
for nx=-2:2
    for ny=-2:2
        in=i+nx; jn=j+ny;
        if((in>0)&&(jn>0)&&(in<=sizeF(1))&&(jn<=sizeF(2))&&(Frozen(in,jn)==1))
            Tpatch(nx+3,ny+3)=T(in,jn);
        end
    end
end
end

# The values in order is 0 if no neighbours in that direction
# 1 if 1e order derivatives is used and 2 if second order
# derivatives are used
Order=zeros(1,4);

# Make 1e order derivatives in x and y direction
Tm(1) = min( Tpatch(2,3) , Tpatch(4,3));
if(isfinite(Tm(1)))
 Order(1)=1;
 end
Tm(2) = min( Tpatch(3,2) , Tpatch(3,4));
if(isfinite(Tm(2)))
 Order(2)=1;
end
# Make 1e order derivatives in cross directions
if(usecross)
    Tm(3) = min( Tpatch(2,2) , Tpatch(4,4));
     if(isfinite(Tm(3)))
      Order(3)=1;
      end
    Tm(4) = min( Tpatch(2,4) , Tpatch(4,2));
    	  if(isfinite(Tm(4)))
		 Order(4)=1;
	   end
end

# Make 2e order derivatives
if(usesecond)
    Tm2=zeros(1,4);
    # pixels with a pixeldistance 2 from the center must be
    # lower in value otherwise use other side or first order
    ch1=(Tpatch(1,3)<Tpatch(2,3))&&isfinite(Tpatch(2,3)); ch2=(Tpatch(5,3)<Tpatch(4,3))&&isfinite(Tpatch(4,3));

    if(ch1&&ch2)
	Tm2(1) =min( (4*Tpatch(2,3)-Tpatch(1,3))/3 , (4*Tpatch(4,3)-Tpatch(5,3))/3);
	 Order(1)=2;
    elseif(ch1)
	 Tm2(1) =(4*Tpatch(2,3)-Tpatch(1,3))/3; Order(1)=2;
    elseif(ch2)
	 Tm2(1) =(4*Tpatch(4,3)-Tpatch(5,3))/3; Order(1)=2;
    end

    ch1=(Tpatch(3,1)<Tpatch(3,2))&&isfinite(Tpatch(3,2)); ch2=(Tpatch(3,5)<Tpatch(3,4))&&isfinite(Tpatch(3,4));

    if(ch1&&ch2)
	Tm2(2) =min( (4*Tpatch(3,2)-Tpatch(3,1))/3 , (4*Tpatch(3,4)-Tpatch(3,5))/3);
	Order(2)=2;
    elseif(ch1)
	 Tm2(2)=(4*Tpatch(3,2)-Tpatch(3,1))/3; Order(2)=2;
    elseif(ch2)
	 Tm2(2)=(4*Tpatch(3,4)-Tpatch(3,5))/3; Order(2)=2;
    end

    if(usecross)
        ch1=(Tpatch(1,1)<Tpatch(2,2))&&isfinite(Tpatch(2,2)); ch2=(Tpatch(5,5)<Tpatch(4,4))&&isfinite(Tpatch(4,4));
        if(ch1&&ch2)
		Tm2(3) =min( (4*Tpatch(2,2)-Tpatch(1,1))/3 , (4*Tpatch(4,4)-Tpatch(5,5))/3); Order(3)=2;
        elseif(ch1)
		 Tm2(3)=(4*Tpatch(2,2)-Tpatch(1,1))/3; Order(3)=2;
        elseif(ch2)
		 Tm2(3)=(4*Tpatch(4,4)-Tpatch(5,5))/3; Order(3)=2;
        end

        ch1=(Tpatch(1,5)<Tpatch(2,4))&&isfinite(Tpatch(2,4)); ch2=(Tpatch(5,1)<Tpatch(4,2))&&isfinite(Tpatch(4,2));
        if(ch1&&ch2)
		Tm2(4) =min( (4*Tpatch(2,4)-Tpatch(1,5))/3 , (4*Tpatch(4,2)-Tpatch(5,1))/3); Order(4)=2;
        elseif(ch1)
		 Tm2(4)=(4*Tpatch(2,4)-Tpatch(1,5))/3; Order(4)=2;
        elseif(ch2)
		 Tm2(4)=(4*Tpatch(4,2)-Tpatch(5,1))/3; Order(4)=2;
        end
    end
else
    Tm2=zeros(1,4);
end

# Calculate the distance using x and y direction
Coeff = [0 0 -1/(max(Fij^2,eps))];
for t=1:2
    if Order(t) == 1
            Coeff=Coeff+[1 -2*Tm(t) Tm(t)^2];
    elseif Order(t) ==2
            Coeff=Coeff+[1 -2*Tm2(t) Tm2(t)^2]*(2.2500);
    end
end

Tt=roots(Coeff); Tt=max(Tt);
# Calculate the distance using the cross directions
if(usecross)
    Coeff = Coeff + [0 0 -1/(max(Fij^2,eps))];
    for t=3:4;
    	if Order(t) == 1
                Coeff=Coeff+0.5*[1 -2*Tm(t) Tm(t)^2];
	elseif Order(t) == 2
                Coeff=Coeff+0.5*[1 -2*Tm2(t) Tm2(t)^2]*(2.2500);
        end
    end
    Tt2=roots(Coeff); Tt2=max(Tt2);
    # Select minimum distance value of both stensils
    if(~isempty(Tt2))
	 Tt=min(Tt,Tt2);
    end
end

# Upwind condition check, current distance must be larger
# then direct neighbours used in solution
DirectNeigbInSol=Tm(isfinite(Tm));
if(nnz(DirectNeigbInSol>=Tt)>0) # Will this ever happen?
    Tt=min(DirectNeigbInSol)+(1/(max(Fij,eps)));
end

#function roots(Coeff)
#    a=Coeff(1);
#    b=Coeff(2);
#    c=Coeff(3);
#    d=max((b*b)-4.0*a*c,0);
#    if a != 0
#        z(1)= (-b - sqrt(d)) / (2.0*a);
#        z(2)= (-b + sqrt(d)) / (2.0*a);
#    else
#        z(1)= (2.0*c)/(-b - sqrt(d));
#        z(2)= (2.0*c)/(-b + sqrt(d));
#    end
#end

end
end


text:
using Test
using ABM4OSN
using Statistics

@test !(true in [((x > 1.0) | (x < -1.0)) for x in [ABM4OSN.generate_opinion() for i in 1:100]])
@test !(true in [(x < 0) for x in [ABM4OSN.generate_inclin_interact() for i in 1:100]])
@test !(true in [((x <= 0.9375) | (x > 1.0)) for x in [ABM4OSN.generate_check_regularity() for i in 1:100]])
@test mean([ABM4OSN.generate_desired_input_count(1) for i in 1:1000]) > 0

ext: jl


text:
function logit(theta, y, x)
    p = 1.0./(1.0 .+ exp.(-x*theta))
    obj = y.*log.(p) .+ (log.(1.0 .- p)).*(1.0 .- y)
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using Fermi.Options
using Formatting

export output

"""
    Fermi.output(str, x...; ending="\\n")

Writes a Python-style String `str` formatted with the values `x...`. The formatting is done through the package Formatting.jl.
The string is printed into the standard output or written into a file depending on the `printstyle` option. See below:

    printstyle:
    "repl"  Default - Returns the string into the Julia REPL (standard output)
    "file"    Write the string into a file. Path must be specified with the keyword `output`
    "both"    Returns the value into stdout and also writes to file
    "none"    Does not write or print anything

In order to write to a file, the keyword `output` must be the path to the destination file.
The default path is "output.jl". For help setting options see `@set`.

# Examples
```
julia> @set {
    printstyle file
    output hydrogen.jl
}
julia> output("The Hydrogen Atom")
julia> output("1s Energy", ending=":")
julia> output(" {:2.2f} hartrees",-0.5)
shell> cat hydrogen.jl
The Hydrogen Atom
1s Energy: -0.50 hartrees
```
"""
function output(str, x...; ending="\n")

    style = Options.get("printstyle")
    f = format(str*ending, x...)

    if style == "none"
        nothing
    elseif style == "repl"
        print(f)
    elseif style == "file"
        path = Options.get("output")
        open(path, "a") do io
            write(io, f)
            flush(io)
        end
    elseif style == "both"
        path = Options.get("output")
        open(path, "a") do io
            write(io, f)
            flush(io)
        end
        print(f)
    else
        throw(FermiException("printing style not recognized: $style. Accepted `printstyle` values: repl, file, both, and none"))
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
@enum msym_geometry begin
        MSYM_GEOMETRY_UNKNOWN = 0
        MSYM_GEOMETRY_SPHERICAL = 1
        MSYM_GEOMETRY_LINEAR = 2
        MSYM_GEOMETRY_PLANAR_REGULAR = 3
        MSYM_GEOMETRY_PLANAR_IRREGULAR = 4
        MSYM_GEOMETRY_POLYHEDRAL_PROLATE = 5
        MSYM_GEOMETRY_POLYHEDRAL_OBLATE = 6
        MSYM_GEOMETRY_ASSYMETRIC = 7
end

struct msym_symmetry_operation
    @enum msym_symmetry_operation_type begin
        MSYM_SYMMETRY_OPERATION_TYPE_IDENTITY = 0
        MSYM_SYMMETRY_OPERATION_TYPE_PROPER_ROTATION = 1
        MSYM_SYMMETRY_OPERATION_TYPE_IMPROPER_ROTATION = 2
        MSYM_SYMMETRY_OPERATION_TYPE_REFLECTION = 3
        MSYM_SYMMETRY_OPERATION_TYPE_INVERSION = 4
    end
    order::Cint                              # Order of proper/improper rotation
    power::Cint                              # Power (e.g. C2^2 = I)
    @enum msym_symmetry_operation_orientation begin
        MSYM_SYMMETRY_OPERATION_ORIENTATION_NONE = 0
        MSYM_SYMMETRY_OPERATION_ORIENTATION_HORIZONTAL = 1
        MSYM_SYMMETRY_OPERATION_ORIENTATION_VERTICAL = 2
        MSYM_SYMMETRY_OPERATION_ORIENTATION_DIHEDRAL = 3
    end
    v::NTuple{3,Cdouble}                            # Proper/improper rotation vector or reflection plane normal
    cla::Cint                                # Class of symmetry operation (point group dependant)
end

@enum msym_point_group_type begin
    MSYM_POINT_GROUP_TYPE_Kh  = 0
    MSYM_POINT_GROUP_TYPE_K   = 1
    MSYM_POINT_GROUP_TYPE_Ci  = 2
    MSYM_POINT_GROUP_TYPE_Cs  = 3
    MSYM_POINT_GROUP_TYPE_Cn  = 4
    MSYM_POINT_GROUP_TYPE_Cnh = 5
    MSYM_POINT_GROUP_TYPE_Cnv = 6
    MSYM_POINT_GROUP_TYPE_Dn  = 7
    MSYM_POINT_GROUP_TYPE_Dnh = 8
    MSYM_POINT_GROUP_TYPE_Dnd = 9
    MSYM_POINT_GROUP_TYPE_Sn  = 10
    MSYM_POINT_GROUP_TYPE_T   = 11
    MSYM_POINT_GROUP_TYPE_Td  = 12
    MSYM_POINT_GROUP_TYPE_Th  = 13
    MSYM_POINT_GROUP_TYPE_O   = 14
    MSYM_POINT_GROUP_TYPE_Oh  = 15
    MSYM_POINT_GROUP_TYPE_I   = 16
    MSYM_POINT_GROUP_TYPE_Ih  = 17
end

struct msym_subgroup
    type::msym_point_group_type
    n::Cint
    order::Cint
    primary::Ref{msym_symmetry_operation}
    #msym_symmetry_operation_t *primary;
    sops::Ref{Ref{msym_symmetry_operation}}
    #msym_symmetry_operation_t **sops;
    generators::NTuple{2,Ref{msym_subgroup}}
    name::NTuple{8,Cchar}
end

struct msym_thresholds
    zero::Cdouble                            # For determining if something is zero (e.g. vectors close to center of mass)
    geometry::Cdouble                        # For translating inertial tensor eigenvalues to geometric structures
    angle::Cdouble                           # For determining angles, (e.g. if vectors are parallel)
    equivalence::Cdouble                     # Equivalence test threshold
    eigfact::Cdouble                         # Jacobi eigenvalue algorithm threshold
    permutation::Cdouble                     # Equality test when determining permutation for symmetry operation
    orthogonalization::Cdouble               # For orthogonalizing orbital subspaces
end

struct msym_element
    id::Ref{Cvoid}                               # custom identifier
    m::Cdouble                              # Mass
    v::NTuple{3,Cdouble}                  # Position
    n::Cint                                 # Nuclear charge
    name::NTuple{4,Cchar}                 # Name
end

struct msym_equivalence_set
    elements::Ref{Ref{msym_element}}             # Pointers to elements
    err::Cdouble                            # Maximum error when detecting this equivalence set
    length::Cint                            # Number of elements
end

struct msym_real_spherical_harmonic
    n::Cint                                # Principal
    l::Cint                                # Azimuthal
    m::Cint                                # Liniear combination of magnetic quantum number (e.g. 2pz = 0, 2px = 1, 2py = -1)
end

struct msym_basis_function
    id::Ref{Cvoid}                               # custom identifier
    @enum msym_basis_type begin
        MSYM_BASIS_TYPE_REAL_SPHERICAL_HARMONIC = 0
        MSYM_BASIS_TYPE_CARTESIAN = 1
    end
    element::Ref{msym_element}
    rsh::msym_real_spherical_harmonic       # Atomic orbital basis
    name::NTuple{8,Cchar}
end

struct msym_partner_function
    i::Cint          # index of partner 0
    d::Cint          # component (dimension)
end

struct msym_salc
    d::Cint             # dimension of space (same as msym_character_table_t.s[msym_subrepresentation_space_t.s].d)
    fl::Cint            # number of basis functions
    pf::Ref{Cvoid}           # partner functions double[d][fl]
    f::Ref{Ref{msym_basis_function}}
end

struct msym_subrepresentation_space
    s::Cint      # symmetry species
    salcl::Cint  # nr of SALCs
    salc::Ref{msym_salc}
end

struct msym_symmetry_species
    d::Cint # dimensionality of a (ir)reducible representation of this species
    r::Cint # sum over all x ct->classc[x]*ct->table[i][x]^2/pg->order (can be decomposed into r irreducible representations)
    name::NTuple{8,Cchar}
end

struct msym_character_table
    d::Cint
    classc::Ref{Cint}
    sops::Ref{Ref{msym_symmetry_operation}}
    s::Ref{msym_symmetry_species}
    table::Ref{Cvoid}  #double[d][d]
end

struct msym_permutation_cycle_t
    l::Cint
    s::Cint
end

struct msym_permutation
    p::Ref{Cint}
    p_length::Cint
    c::Ref{msym_permutation_cycle_t}
    c_length::Cint
end

struct msym_context
    thresholds::Ref{msym_thresholds}
    elements::Ref{msym_element}
    pelements::Ref{Ref{msym_element}}
    basis::Ref{msym_basis_function}
    es::Ref{msym_equivalence_set}
    es_perm::Ref{Ref{msym_permutation}}
    srs::Ref{msym_subrepresentation_space}
    srsbf::Ref{Ref{msym_basis_function}}
    srs_span::Ref{Cint}
    flags::Culong
    elementsl::Cint
    basisl::Cint
    esl::Cint
    srsl::Cint
    es_perml::Cint
    sgl::Cint
    pg::Ref{msym_point_group_type}
    sg::Ref{msym_subgroup}
    cm::NTuple{3,Cdouble}
    geometry::msym_geometry
    eigval::NTuple{3,Cdouble}
    eigvec::NTuple{3,NTuple{Cdouble}}
    struct _external_data
        eesmap::Ref{Ref{msym_equivalence_set}}
        set_elements_ptr::Ref{msym_element}
        elements::Ref{msym_element}
        es::Ref{msym_equivalence_set}
    end
end

ctx = ccall((:msymCreateContext, "libmsym.so"), msym_context, ())


text:
using OodleAnalytics
using Documenter

makedocs(;
    modules=[OodleAnalytics],
    authors="tobi-lipede-oodle <tobi.lipede@oodlefinance.com> and contributors",
    repo="https://github.com/tobi-lipede-oodle/OodleAnalytics.jl/blob/{commit}{path}#L{line}",
    sitename="OodleAnalytics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tobi-lipede-oodle.github.io/OodleAnalytics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/tobi-lipede-oodle/OodleAnalytics.jl",
)


text:
# Rewriting the node graph to prepare it for translation.

# Auxiliary nodes.

# A SQL subquery with an undetermined SELECT args.
mutable struct BoxNode <: TabularNode
    over::Union{SQLNode, Nothing}
    type::BoxType
    handle::Int
    refs::Vector{SQLNode}

    BoxNode(; over = nothing, type = EMPTY_BOX, handle = 0, refs = SQLNode[]) =
        new(over, type, handle, refs)
end

Box(args...; kws...) =
    BoxNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Box), pats::Vector{Any}) =
    dissect(scr, BoxNode, pats)

function PrettyPrinting.quoteof(n::BoxNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Box))
    if !ctx.limit
        if n.type !== EMPTY_BOX
            push!(ex.args, Expr(:kw, :type, quoteof(n.type)))
        end
        if n.handle != 0
            push!(ex.args, Expr(:kw, :handle, n.handle))
        end
        if !isempty(n.refs)
            push!(ex.args, Expr(:kw, :refs, Expr(:vect, quoteof(n.refs, ctx)...)))
        end
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::BoxNode) =
    n.type.name

rebase(n::BoxNode, n′) =
    BoxNode(over = rebase(n.over, n′),
            type = n.type, handle = n.handle, refs = n.refs)

box_type(n::BoxNode) =
    n.type

box_type(n::SQLNode) =
    box_type(n[]::BoxNode)

# Get(over = Get(:a), name = :b) => NameBound(over = Get(:b), name = :a)
mutable struct NameBoundNode <: AbstractSQLNode
    over::SQLNode
    name::Symbol

    NameBoundNode(; over, name) =
        new(over, name)
end

NameBound(args...; kws...) =
    NameBoundNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(NameBound), pats::Vector{Any}) =
    dissect(scr, NameBoundNode, pats)

PrettyPrinting.quoteof(n::NameBoundNode, ctx::QuoteContext) =
    Expr(:call, nameof(NameBound), Expr(:kw, :over, quoteof(n.over, ctx)), Expr(:kw, :name, QuoteNode(n.name)))

# Get(over = q, name = :b) => HandleBound(over = Get(:b), handle = get_handle(q))
mutable struct HandleBoundNode <: AbstractSQLNode
    over::SQLNode
    handle::Int

    HandleBoundNode(; over, handle) =
        new(over, handle)
end

HandleBound(args...; kws...) =
    HandleBoundNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(HandleBound), pats::Vector{Any}) =
    dissect(scr, HandleBoundNode, pats)

PrettyPrinting.quoteof(n::HandleBoundNode, ctx::QuoteContext) =
    Expr(:call, nameof(NameBound), Expr(:kw, :over, quoteof(n.over, ctx)), Expr(:kw, :handle, n.handle))

# A generic From node is specialized to FromNothing, FromTable, or FromReference.
mutable struct FromNothingNode <: TabularNode
end

FromNothing(args...; kws...) =
    FromNothingNode(args...; kws...) |> SQLNode

PrettyPrinting.quoteof(::FromNothingNode, ::QuoteContext) =
    Expr(:call, nameof(FromNothing))

mutable struct FromTableNode <: TabularNode
    table::SQLTable

    FromTableNode(; table) =
        new(table)
end

FromTable(args...; kws...) =
    FromTableNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::FromTableNode, ctx::QuoteContext)
    tex = get(ctx.vars, n.table, nothing)
    if tex === nothing
        tex = quoteof(n.table, limit = true)
    end
    Expr(:call, nameof(FromTable), Expr(:kw, :table, tex))
end

mutable struct FromReferenceNode <: TabularNode
    over::SQLNode
    name::Symbol

    FromReferenceNode(; over, name) =
        new(over, name)
end

FromReference(args...; kws...) =
    FromReferenceNode(args...; kws...) |> SQLNode

PrettyPrinting.quoteof(n::FromReferenceNode, ctx::QuoteContext) =
    Expr(:call,
         nameof(FromReference),
         Expr(:kw, :over, quoteof(n.over, ctx)),
         Expr(:kw, :name, QuoteNode(n.name)))

# Annotated Bind node.
mutable struct IntBindNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    args::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}
    owned::Bool     # Did we find the outer query for this node?

    function IntBindNode(; over = nothing, args, label_map = nothing, owned = false)
        if label_map !== nothing
            return new(over, args, label_map, owned)
        end
        n = new(over, args, OrderedDict{Symbol, Int}(), owned)
        for (i, arg) in enumerate(n.args)
            n.label_map[label(arg)] = i
        end
        n
    end
end

IntBind(args...; kws...) =
    IntBindNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::IntBindNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(IntBind))
    push!(ex.args, Expr(:kw, :args, Expr(:vect, quoteof(n.args, ctx)...)))
    push!(ex.args, Expr(:kw, :owned, n.owned))
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

rebase(n::IntBindNode, n′) =
    IntBindNode(over = rebase(n.over, n′),
                     args = n.args, label_map = n.label_map, owned = n.owned)

# A recursive UNION ALL node.
mutable struct KnotNode <: TabularNode
    over::Union{SQLNode, Nothing}
    name::Symbol
    box::BoxNode
    iterator::SQLNode
    iterator_boxes::Vector{BoxNode}

    KnotNode(; over = nothing, iterator, name = label(iterator), iterator_boxes = SQLNode[], box) =
        new(over, name, box, iterator, iterator_boxes)
end

KnotNode(iterator; over = nothing, box) =
    KnotNode(over = over, iterator = iterator, box = box)

Knot(args...; kws...) =
    KnotNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::KnotNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Knot))
    if !ctx.limit
        push!(ex.args, quoteof(n.iterator, ctx))
    else
        push!(ex.args, :…)
    end
    push!(ex.args, Expr(:kw, :name, QuoteNode(n.name)))
    if !ctx.limit
        box_ex = Expr(:ref, quoteof(SQLNode(n.box), ctx))
        push!(ex.args, Expr(:kw, :box, box_ex))
        push!(ex.args, Expr(:kw, :iterator, quoteof(n.iterator, ctx)))
        iterator_boxes_ex =
            Expr(:vect, Any[Expr(:ref, quoteof(SQLNode(iterator_box), ctx))
                            for iterator_box in n.iterator_boxes]...)
        push!(ex.args, Expr(:kw, :iterator_boxes, iterator_boxes_ex))
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::KnotNode) =
    n.name

rebase(n::KnotNode, n′) =
    KnotNode(over = rebase(n.over, n′),
             name = n.name, box = n.box, iterator = n.iterator, iterator_boxes = n.iterator_boxes)

# Iterate node is split into Knot and IntIterate.
mutable struct IntIterateNode <: TabularNode
    over::Union{SQLNode, Nothing}
    name::Symbol            # Original label.
    iterator_name::Symbol   # Label of the iterator.

    IntIterateNode(; over = nothing, name, iterator_name) =
        new(over, name, iterator_name)
end

IntIterate(args...; kws...) =
    IntIterateNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::IntIterateNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(IntIterate))
    push!(ex.args, Expr(:kw, :name, QuoteNode(n.name)))
    push!(ex.args, Expr(:kw, :iterator_name, QuoteNode(n.iterator_name)))
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::IntIterateNode) =
    n.name

rebase(n::IntIterateNode, n′) =
    IntIterateNode(over = rebase(n.over, n′), name = n.name, iterator_name = n.iterator_name)

# Annotated Join node.
mutable struct IntJoinNode <: TabularNode
    over::Union{SQLNode, Nothing}
    joinee::SQLNode
    on::SQLNode
    left::Bool
    right::Bool
    type::BoxType               # Type of the product of `over` and `joinee`.
    lateral::Vector{SQLNode}    # References from `joinee` to `over` for JOIN LATERAL.

    IntJoinNode(; over, joinee, on, left, right, type = EMPTY_BOX, lateral = SQLNode[]) =
        new(over, joinee, on, left, right, type, lateral)
end

IntJoinNode(joinee, on; over = nothing, left = false, right = false, type = EMPTY_BOX, lateral = SQLNode[]) =
    IntJoinNode(over = over, joinee = joinee, on = on, left = left, right = right, type = type, lateral = lateral)

IntJoin(args...; kws...) =
    IntJoinNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::IntJoinNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(IntJoin))
    if !ctx.limit
        push!(ex.args, quoteof(n.joinee, ctx))
        push!(ex.args, quoteof(n.on, ctx))
        if n.left
            push!(ex.args, Expr(:kw, :left, n.left))
        end
        if n.right
            push!(ex.args, Expr(:kw, :right, n.right))
        end
        if n.type !== EMPTY_BOX
            push!(ex.args, Expr(:kw, :type, n.type))
        end
        if !isempty(n.lateral)
            push!(ex.args, Expr(:kw, :lateral, Expr(:vect, quoteof(n.lateral, ctx)...)))
        end
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

rebase(n::IntJoinNode, n′) =
    IntJoinNode(over = rebase(n.over, n′),
                     joinee = n.joinee, on = n.on, left = n.left, right = n.right, type = n.type, lateral = n.lateral)

label(n::Union{NameBoundNode, HandleBoundNode, IntBindNode, IntJoinNode}) =
    label(n.over)

# Annotation context.

# Maps a node in the annotated graph to a path in the original graph (for error reporting).
struct PathMap
    paths::Vector{Tuple{SQLNode, Int}}
    origins::IdDict{Any, Int}

    PathMap() =
        new(Tuple{SQLNode, Int}[], IdDict{Any, Int}())
end

function get_path(map::PathMap, idx::Int)
    path = SQLNode[]
    while idx != 0
        n, idx = map.paths[idx]
        push!(path, n)
    end
    path
end

get_path(map::PathMap, n) =
    get_path(map, get(map.origins, n, 0))

struct AnnotateContext
    path_map::PathMap
    current_path::Vector{Int}
    handles::Dict{SQLNode, Int}
    boxes::Vector{BoxNode}
    with_nodes::Dict{Symbol, SQLNode}

    AnnotateContext() =
        new(PathMap(),
            Int[0],
            Dict{SQLNode, Int}(),
            BoxNode[],
            Dict{Symbol, SQLNode}())

    AnnotateContext(ctx::AnnotateContext; with_nodes = nothing) =
        new(ctx.path_map, ctx.current_path, ctx.handles, ctx.boxes, something(with_nodes, ctx.with_nodes))
end

function grow_path!(ctx::AnnotateContext, n::SQLNode)
    push!(ctx.path_map.paths, (n, ctx.current_path[end]))
    push!(ctx.current_path, length(ctx.path_map.paths))
end

function shrink_path!(ctx::AnnotateContext)
    pop!(ctx.current_path)
end

function mark_origin!(ctx::AnnotateContext, n::SQLNode)
    ctx.path_map.origins[n] = ctx.current_path[end]
end

mark_origin!(ctx::AnnotateContext, n::AbstractSQLNode) =
    mark_origin!(ctx, convert(SQLNode, n))

get_path(ctx::AnnotateContext) =
    get_path(ctx.path_map, ctx.current_path[end])

get_path(ctx::AnnotateContext, n::SQLNode) =
    get_path(ctx.path_map, n)

function make_handle!(ctx::AnnotateContext, n::SQLNode)
    get!(ctx.handles, n) do
        length(ctx.handles) + 1
    end
end

function get_handle(ctx::AnnotateContext, n::SQLNode)
    handle = 0
    idx = get(ctx.path_map.origins, n, 0)
    if idx > 0
        n = ctx.path_map.paths[idx][1]
        handle = get(ctx.handles, n, 0)
    end
    handle
end

get_handle(ctx::AnnotateContext, ::Nothing) =
    0

# Rewriting of the node graph.

function annotate(n::SQLNode, ctx)
    grow_path!(ctx, n)
    n′ = convert(SQLNode, annotate(n[], ctx))
    mark_origin!(ctx, n′)
    box = BoxNode(over = n′)
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    mark_origin!(ctx, n′)
    shrink_path!(ctx)
    n′
end

function annotate_scalar(n::SQLNode, ctx)
    grow_path!(ctx, n)
    n′ = convert(SQLNode, annotate_scalar(n[], ctx))
    mark_origin!(ctx, n′)
    shrink_path!(ctx)
    n′
end

annotate(ns::Vector{SQLNode}, ctx) =
    SQLNode[annotate(n, ctx) for n in ns]

annotate_scalar(ns::Vector{SQLNode}, ctx) =
    SQLNode[annotate_scalar(n, ctx) for n in ns]

function annotate(::Nothing, ctx)
    box = BoxNode()
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    mark_origin!(ctx, n′)
    n′
end

annotate_scalar(::Nothing, ctx) =
    nothing

annotate(n::AbstractSQLNode, ctx) =
    throw(IllFormedError(path = get_path(ctx)))

function annotate_scalar(n::TabularNode, ctx)
    n′ = convert(SQLNode, annotate(n, ctx))
    mark_origin!(ctx, n′)
    box = BoxNode(over = n′)
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    n′
end

function rebind(node, base, ctx)
    while @dissect(node, over |> Get(name = name))
        mark_origin!(ctx, base)
        base = NameBound(over = base, name = name)
        node = over
    end
    if node !== nothing
        handle = make_handle!(ctx, node)
        mark_origin!(ctx, base)
        base = HandleBound(over = base, handle = handle)
    end
    base
end

function annotate_scalar(n::AggregateNode, ctx)
    args′ = annotate_scalar(n.args, ctx)
    filter′ = annotate_scalar(n.filter, ctx)
    n′ = Agg(name = n.name, distinct = n.distinct, args = args′, filter = filter′)
    rebind(n.over, n′, ctx)
end

function annotate(n::AppendNode, ctx)
    over′ = annotate(n.over, ctx)
    args′ = annotate(n.args, ctx)
    Append(over = over′, args = args′)
end

function annotate(n::AsNode, ctx)
    over′ = annotate(n.over, ctx)
    As(over = over′, name = n.name)
end

function annotate_scalar(n::AsNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    As(over = over′, name = n.name)
end

function annotate(n::BindNode, ctx)
    over′ = annotate(n.over, ctx)
    args′ = annotate_scalar(n.args, ctx)
    IntBind(over = over′, args = args′, label_map = n.label_map)
end

annotate_scalar(n::BindNode, ctx) =
    annotate(n, ctx)

function annotate(n::DefineNode, ctx)
    over′ = annotate(n.over, ctx)
    args′ = annotate_scalar(n.args, ctx)
    Define(over = over′, args = args′, label_map = n.label_map)
end

function annotate(n::FromNode, ctx)
    source = n.source
    if source isa SQLTable
        FromTable(table = source)
    elseif source isa Symbol
        over = get(ctx.with_nodes, source, nothing)
        if over === nothing
            throw(ReferenceError(REFERENCE_ERROR_TYPE.UNDEFINED_TABLE_REFERENCE,
                                 name = source,
                                 path = get_path(ctx)))
        end
        FromReference(over = over, name = source)
    else
        FromNothing()
    end
end

function annotate_scalar(n::FunctionNode, ctx)
    args′ = annotate_scalar(n.args, ctx)
    Fun(name = n.name, args = args′)
end

function annotate_scalar(n::GetNode, ctx)
    rebind(n.over, Get(name = n.name), ctx)
end

function annotate(n::GroupNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    Group(over = over′, by = by′, label_map = n.label_map)
end

function annotate(n::HighlightNode, ctx)
    over′ = annotate(n.over, ctx)
    Highlight(over = over′, color = n.color)
end

function annotate_scalar(n::HighlightNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    Highlight(over = over′, color = n.color)
end

function annotate(n::IterateNode, ctx)
    over′ = annotate(n.over, ctx)
    knot_box = BoxNode()
    knot = KnotNode(over = over′, iterator = n.iterator, box = knot_box)
    mark_origin!(ctx, knot)
    knot_box.over = knot
    push!(ctx.boxes, knot_box)
    over′ = convert(SQLNode, knot_box)
    mark_origin!(ctx, over′)
    with_nodes′ = copy(ctx.with_nodes)
    with_nodes′[knot.name] = over′
    ctx′ = AnnotateContext(ctx, with_nodes = with_nodes′)
    range_start = length(ctx.boxes) + 1
    iterator′ = annotate(n.iterator, ctx′)
    range_stop = length(ctx.boxes)
    knot.iterator = iterator′
    knot.iterator_boxes = ctx.boxes[range_start:range_stop]
    IntIterateNode(over = over′, name = label(n.over), iterator_name = knot.name)
end

function annotate(n::JoinNode, ctx)
    over′ = annotate(n.over, ctx)
    joinee′ = annotate(n.joinee, ctx)
    on′ = annotate_scalar(n.on, ctx)
    IntJoin(over = over′, joinee = joinee′, on = on′, left = n.left, right = n.right)
end

function annotate(n::LimitNode, ctx)
    over′ = annotate(n.over, ctx)
    Limit(over = over′, offset = n.offset, limit = n.limit)
end

annotate_scalar(n::LiteralNode, ctx) =
    n

function annotate(n::OrderNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    Order(over = over′, by = by′)
end

function annotate(n::PartitionNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    order_by′ = annotate_scalar(n.order_by, ctx)
    Partition(over = over′, by = by′, order_by = order_by′, frame = n.frame)
end

function annotate(n::SelectNode, ctx)
    over′ = annotate(n.over, ctx)
    args′ = annotate_scalar(n.args, ctx)
    Select(over = over′, args = args′, label_map = n.label_map)
end

function annotate_scalar(n::SortNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    Sort(over = over′, value = n.value, nulls = n.nulls)
end

annotate_scalar(n::VariableNode, ctx) =
    n

function annotate(n::WhereNode, ctx)
    over′ = annotate(n.over, ctx)
    condition′ = annotate_scalar(n.condition, ctx)
    Where(over = over′, condition = condition′)
end

function annotate(n::WithNode, ctx)
    args′ = annotate(n.args, ctx)
    with_nodes′ = copy(ctx.with_nodes)
    for (name, i) in n.label_map
        with_nodes′[name] = args′[i]
    end
    ctx′ = AnnotateContext(ctx, with_nodes = with_nodes′)
    over′ = annotate(n.over, ctx′)
    With(over = over′, args = args′, label_map = n.label_map)
end

# Type resolution.

resolve!(ctx::AnnotateContext) =
    resolve!(ctx.boxes, ctx)

function resolve!(boxes::AbstractVector{BoxNode}, ctx)
    for box in boxes
        over = box.over
        if over !== nothing
            h = get_handle(ctx, over)
            t = resolve(over[], ctx)
            t = add_handle(t, h)
            box.handle = h
            box.type = t
        end
    end
end

function resolve(n::AppendNode, ctx)
    t = box_type(n.over)
    for arg in n.args
        t = intersect(t, box_type(arg))
    end
    t
end

function resolve(n::Union{AsNode, KnotNode}, ctx)
    t = box_type(n.over)
    fields = FieldTypeMap(n.name => t.row)
    row = RowType(fields)
    BoxType(n.name, row, t.handle_map)
end

function resolve(n::DefineNode, ctx)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for (f, ft) in t.row.fields
        if f in keys(n.label_map)
            ft = ScalarType()
        end
        fields[f] = ft
    end
    for f in keys(n.label_map)
        if !haskey(fields, f)
            fields[f] = ScalarType()
        end
    end
    row = RowType(fields, t.row.group)
    BoxType(t.name, row, t.handle_map)
end

resolve(n::FromNothingNode, ctx) =
    EMPTY_BOX

function resolve(n::FromReferenceNode, ctx)
    t = box_type(n.over)
    ft = get(t.row.fields, n.name, nothing)
    if !(ft isa RowType)
        throw(ReferenceError(REFERENCE_ERROR_TYPE.INVALID_TABLE_REFERENCE,
                             name = n.name,
                             path = get_path(ctx, n.over)))
    end
    BoxType(n.name, ft)
end

function resolve(n::FromTableNode, ctx)
    fields = FieldTypeMap()
    for f in n.table.columns
        fields[f] = ScalarType()
    end
    row = RowType(fields)
    BoxType(n.table.name, row)
end

function resolve(n::GroupNode, ctx)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for name in keys(n.label_map)
        fields[name] = ScalarType()
    end
    row = RowType(fields, t.row)
    BoxType(t.name, row)
end

resolve(n::Union{HighlightNode, IntBindNode, LimitNode, OrderNode, WhereNode, WithNode}, ctx) =
    box_type(n.over)

resolve_knot!(n::SQLNode, ctx) =
    resolve_knot!(n[], ctx)

function resolve_knot!(n::BoxNode, ctx)
    knot = n.over[]::KnotNode
    iterator_t = box_type(knot.iterator)
    while !issubset(n.type.row, iterator_t.row)
        n.type = intersect(n.type, iterator_t)
        resolve!(knot.iterator_boxes, ctx)
        iterator_t = box_type(knot.iterator)
    end
    over = n.over
    n.type = add_handle(n.type, n.handle)
end

function resolve(n::IntIterateNode, ctx)
    resolve_knot!(n.over, ctx)
    t = box_type(n.over)
    row = t.row.fields[n.iterator_name]::RowType
    BoxType(n.name, row)
end

function resolve(n::IntJoinNode, ctx)
    lt = box_type(n.over)
    rt = box_type(n.joinee)
    t = union(lt, rt)
    n.type = t
    t
end

function resolve(n::PartitionNode, ctx)
    t = box_type(n.over)
    row = RowType(t.row.fields, t.row)
    BoxType(t.name, row, t.handle_map)
end

function resolve(n::SelectNode, ctx)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for name in keys(n.label_map)
        fields[name] = ScalarType()
    end
    row = RowType(fields)
    BoxType(t.name, row)
end

# Collecting references.

gather!(refs::Vector{SQLNode}, n::SQLNode) =
    gather!(refs, n[])

function gather!(refs::Vector{SQLNode}, ns::Vector{SQLNode})
    for n in ns
        gather!(refs, n)
    end
end

gather!(refs::Vector{SQLNode}, ::Union{AbstractSQLNode, Nothing}) =
    nothing

gather!(refs::Vector{SQLNode}, n::Union{AsNode, BoxNode, HighlightNode, SortNode}) =
    gather!(refs, n.over)

function gather!(refs::Vector{SQLNode}, n::IntBindNode)
    gather!(refs, n.over)
    gather!(refs, n.args)
    n.owned = true
end

gather!(refs::Vector{SQLNode}, n::FunctionNode) =
    gather!(refs, n.args)

function gather!(refs::Vector{SQLNode}, n::Union{AggregateNode, GetNode, HandleBoundNode, NameBoundNode})
    push!(refs, n)
end

# Validating references.

function validate(t::BoxType, ref::SQLNode, ctx)
    if @dissect(ref, over |> HandleBound(handle = handle))
        if handle in keys(t.handle_map)
            ht = t.handle_map[handle]
            if ht isa AmbiguousType
                throw(ReferenceError(REFERENCE_ERROR_TYPE.AMBIGUOUS_HANDLE,
                                     path = get_path(ctx, ref)))
            end
            validate(ht, over, ctx)
        else
            throw(ReferenceError(REFERENCE_ERROR_TYPE.UNDEFINED_HANDLE,
                                 path = get_path(ctx, ref)))
        end
    else
        validate(t.row, ref, ctx)
    end
end

function validate(t::RowType, ref::SQLNode, ctx)
    while @dissect(ref, over |> NameBound(name = name))
        ft = get(t.fields, name, EmptyType())
        if !(ft isa RowType)
            type =
                ft isa EmptyType ? REFERENCE_ERROR_TYPE.UNDEFINED_NAME :
                ft isa ScalarType ? REFERENCE_ERROR_TYPE.UNEXPECTED_SCALAR_TYPE :
                ft isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_NAME : error()
            throw(ReferenceError(type, name = name, path = get_path(ctx, ref)))
        end
        t = ft
        ref = over
    end
    if @dissect(ref, nothing |> Get(name = name))
        ft = get(t.fields, name, EmptyType())
        if !(ft isa ScalarType)
            type =
                ft isa EmptyType ? REFERENCE_ERROR_TYPE.UNDEFINED_NAME :
                ft isa RowType ? REFERENCE_ERROR_TYPE.UNEXPECTED_ROW_TYPE :
                ft isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_NAME : error()
            throw(ReferenceError(type, name = name, path = get_path(ctx, ref)))
        end
    elseif @dissect(ref, nothing |> Agg(name = name))
        if !(t.group isa RowType)
            type =
                t.group isa EmptyType ? REFERENCE_ERROR_TYPE.UNEXPECTED_AGGREGATE :
                t.group isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_AGGREGATE : error()
            throw(ReferenceError(type, path = get_path(ctx, ref)))
        end
    else
        error()
    end
end

function gather_and_validate!(refs::Vector{SQLNode}, n, t::BoxType, ctx)
    start = length(refs) + 1
    gather!(refs, n)
    for k in start:length(refs)
        validate(t, refs[k], ctx)
    end
end

function route(lt::BoxType, rt::BoxType, ref::SQLNode)
    if @dissect(ref, over |> HandleBound(handle = handle))
        if get(lt.handle_map, handle, EmptyType()) isa EmptyType
            return 1
        else
            return -1
        end
    end
    return route(lt.row, rt.row, ref)
end

function route(lt::RowType, rt::RowType, ref::SQLNode)
    while @dissect(ref, over |> NameBound(name = name))
        lt′ = get(lt.fields, name, EmptyType())
        if lt′ isa EmptyType
            return 1
        end
        rt′ = get(rt.fields, name, EmptyType())
        if rt′ isa EmptyType
            return -1
        end
        @assert lt′ isa RowType && rt′ isa RowType
        lt = lt′
        rt = rt′
        ref = over
    end
    if @dissect(ref, Get(name = name))
        if name in keys(lt.fields)
            return -1
        else
            return 1
        end
    elseif @dissect(ref, over |> Agg(name = name))
        if lt.group isa RowType
            return -1
        else
            return 1
        end
    else
        error()
    end
end

# Linking references through box nodes.

function link!(ctx::AnnotateContext)
    root_box = ctx.boxes[end]
    for (f, ft) in root_box.type.row.fields
        if ft isa ScalarType
            push!(root_box.refs, Get(f))
        end
    end
    link!(reverse(ctx.boxes), ctx)
end

function link!(boxes::AbstractVector{BoxNode}, ctx)
    for box in boxes
        box.over !== nothing || continue
        refs′ = SQLNode[]
        for ref in box.refs
            if @dissect(ref, over |> HandleBound(handle = handle)) && handle == box.handle
                push!(refs′, over)
            else
                push!(refs′, ref)
            end
        end
        link!(box.over[], refs′, ctx)
    end
end

function link!(n::AppendNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    for arg in n.args
        box = arg[]::BoxNode
        append!(box.refs, refs)
    end
end

function link!(n::AsNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        if @dissect(ref, over |> NameBound(name = name))
            @assert name == n.name
            push!(box.refs, over)
        elseif @dissect(ref, HandleBound())
            push!(box.refs, ref)
        else
            error()
        end
    end
end

function link!(n::DefineNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    seen = Set{Symbol}()
    for ref in refs
        if @dissect(ref, nothing |> Get(name = name)) && name in keys(n.label_map)
            !(name in seen) || continue
            push!(seen, name)
            col = n.args[n.label_map[name]]
            gather_and_validate!(box.refs, col, box.type, ctx)
        else
            push!(box.refs, ref)
        end
    end
end

link!(::Union{FromNothingNode, FromTableNode}, ::Vector{SQLNode}, ctx) =
    nothing

function link!(n::FromReferenceNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        push!(box.refs, NameBound(over = ref, name = n.name))
    end
end

function link!(n::GroupNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    gather_and_validate!(box.refs, n.by, box.type, ctx)
    for ref in refs
        if @dissect(ref, nothing |> Agg(args = args, filter = filter))
            gather_and_validate!(box.refs, args, box.type, ctx)
            if filter !== nothing
                gather_and_validate!(box.refs, filter, box.type, ctx)
            end
        end
    end
end

function link!(n::Union{HighlightNode, LimitNode, WithNode}, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
end

function link!(n::IntBindNode, refs::Vector{SQLNode}, ctx)
    if !n.owned
        gather_and_validate!(SQLNode[], n.args, EMPTY_BOX, ctx)
    end
    box = n.over[]::BoxNode
    append!(box.refs, refs)
end

function link!(n::IntIterateNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        push!(box.refs, NameBound(over = ref, name = n.iterator_name))
    end
end

function link!(n::IntJoinNode, refs::Vector{SQLNode}, ctx)
    lbox = n.over[]::BoxNode
    rbox = n.joinee[]::BoxNode
    gather_and_validate!(n.lateral, n.joinee, lbox.type, ctx)
    append!(lbox.refs, n.lateral)
    refs′ = SQLNode[]
    gather_and_validate!(refs′, n.on, n.type, ctx)
    append!(refs′, refs)
    for ref in refs′
        turn = route(lbox.type, rbox.type, ref)
        if turn < 0
            push!(lbox.refs, ref)
        else
            push!(rbox.refs, ref)
        end
    end
end

function link!(n::KnotNode, ::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    iterator_box = n.iterator[]::BoxNode
    refs = SQLNode[]
    seen = Set{SQLNode}()
    while true
        repeat = false
        for ref in n.box.refs
            @dissect(ref, over |> NameBound(name = name)) || error()
            @assert name == n.name ref
            if !(over in seen)
                push!(refs, ref)
                push!(seen, over)
                repeat = true
            end
        end
        empty!(n.box.refs)
        append!(n.box.refs, refs)
        repeat || break
        for ibox in n.iterator_boxes
            empty!(ibox.refs)
        end
        append!(iterator_box.refs, refs)
        link!(reverse(n.iterator_boxes), ctx)
    end
    for ref in refs
        @dissect(ref, over |> NameBound(name = name)) || error()
        push!(box.refs, over)
    end
end

function link!(n::OrderNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    gather_and_validate!(box.refs, n.by, box.type, ctx)
end

function link!(n::PartitionNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        if @dissect(ref, nothing |> Agg(args = args, filter = filter))
            gather_and_validate!(box.refs, args, box.type, ctx)
            if filter !== nothing
                gather_and_validate!(box.refs, filter, box.type, ctx)
            end
        else
            push!(box.refs, ref)
        end
    end
    gather_and_validate!(box.refs, n.by, box.type, ctx)
    gather_and_validate!(box.refs, n.order_by, box.type, ctx)
end

function link!(n::SelectNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    gather_and_validate!(box.refs, n.args, box.type, ctx)
end

function link!(n::WhereNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    gather_and_validate!(box.refs, n.condition, box.type, ctx)
end


text:
## --- Bootstrap resampling

    # Kernel functions for bsr!
    @inline function uniform(rng::AbstractRNG, mu, halfwidth)
        mu + halfwidth * (2*rand(rng)-1)
    end
    export uniform
    @inline function triangular(rng::AbstractRNG, mu, halfwidth)
        mu + halfwidth * (rand(rng) - rand(rng))
    end
    export triangular
    @inline function gaussian(rng::AbstractRNG, mu, sigma)
        mu + sigma * randn(rng)
    end
    export gaussian

    """
    ```julia
    bsr!([f::Function], resampled::Array, index::Vector{Int}, data, sigma, p;
        \trng::AbstractRNG=MersenneTwister()
    )
    ```

    Fill `resampled` with data boostrap resampled from a (sample-per-row / element-per-column)
    dataset `data` with uncertainties `sigma` and resampling probabilities `p`, optionally using
    random numbers generated by `f` where `f` is a function of the form `f(rng, data[i], sigma[i])`
    """
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::Number, p::Number; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma)
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::Number, p::AbstractVector; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma)
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::AbstractArray, p::Number; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[row,j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function bsr!(f::Function, resampled::Array, index::Vector{Int}, data::AbstractArray, sigma::AbstractArray, p::AbstractVector; rng::AbstractRNG=MersenneTwister())
        # Prepare
        ndata = size(data,1)
        nrows = size(resampled,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        nrows_accepted = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    index[nrows_accepted] = i
                    if nrows_accepted == nrows
                        break
                    end
                end
            end
            nrows_new = min(nrows_accepted - n, nrows - n)

            # Columns go in outer loop because of column major indexing
            for j=1:ncolumns
                # Optimized inner loop
                @inbounds for i = 1:nrows_new
                    row = index[n+i]
                    resampled[n+i,j] = f(rng, data[row,j], sigma[row,j])
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    # default method if f not specified
    bsr!(resampled::Array, index::Vector{Int}, data::AbstractArray, sigma, p; rng=MersenneTwister()) = bsr!(gaussian, resampled, index, data, sigma, p, rng=rng)
    export bsr!

    """
    ```julia
    resampled = bsresample(data::AbstractArray, sigma, nrows, [p];
        \t kernel = gaussian,
        \t rng = MersenneTwister(),
        \t return_index = false
    )
    ```
    Bootstrap resample a (sample-per-row / element-per-column) array of `data`
    with uncertainties `sigma` and resampling probabilities `p`
    """
    function bsresample(data::AbstractArray, sigma, nrows, p=min(0.2,nrows/size(data,1));
            kernel = gaussian,
            rng = MersenneTwister(),
            return_index = false
        )
        index = Array{Int}(undef, nrows)
        resampled = Array{float(eltype(data))}(undef, nrows, size(data,2))
        bsr!(kernel, resampled, index, data, sigma, p, rng=rng)
        if return_index
            return resampled, index
        else
            return resampled
        end
    end
    """
    ```julia
    resampled = bsresample(dataset::Dict, nrows, [elements], [p];
        \t kernel = gaussian,
        \t rng = MersenneTwister()
    )
    ```
    Bootstrap resample a dictionary-based `dataset` with uncertainties stored either
    in `dataset["err"]` or `dataset["[variable]_sigma"]`
    """
    function bsresample(dataset::Dict, nrows, elements=dataset["elements"], p=min(0.2,nrows/length(elements[1]));
            kernel = gaussian,
            rng = MersenneTwister()
        )
        # 2d array of nominal values
        data = unelementify(dataset, elements, floatout=true)

        # 2d array of absolute 1-sigma uncertainties
        if haskey(dataset, "err") && isa(dataset["err"], Dict)
            sigma = unelementify(dataset["err"], elements, floatout=true)
        else
            sigma = unelementify(dataset, elements.*"_sigma", floatout=true)
        end

        # Resample
        sdata = bsresample(data, sigma, nrows, p, kernel=kernel, rng=rng)
        return elementify(sdata, elements, skipstart=0)
    end
    export bsresample

    function randsample!(resampled::AbstractArray, data::AbstractArray, nrows::Integer, p::Number, rng::AbstractRNG=MersenneTwister(), buffer::Vector{Int}=Array{Int}(undef,size(data,1)))
        # Prepare
        ndata = size(data,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            nrows_accepted = 0
            @inbounds for i=1:ndata
                if rand(rng) < p
                    nrows_accepted += 1
                    buffer[nrows_accepted] = i
                end
            end
            nrows_new = min(nrows_accepted, nrows - n)

            # Columns go in outer loop because of column major indexing
            @inbounds @simd for j=1:ncolumns
                # Optimized inner loop
                for i = 1:nrows_new
                    resampled[n+i,j] = data[buffer[i],j]
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end
    function randsample!(resampled::AbstractArray, data::AbstractArray, nrows::Integer, p::AbstractVector, rng::AbstractRNG=MersenneTwister(), buffer::Vector{Int}=Array{Int}(undef,size(data,1)))
        # Prepare
        ndata = size(data,1)
        ncolumns = size(resampled,2)

        # Resample
        n = 0
        while n < nrows

            # Compare acceptance probability p against Unif(0,1)
            nrows_accepted = 0
            @inbounds for i=1:ndata
                if rand(rng) < p[i]
                    nrows_accepted += 1
                    buffer[nrows_accepted] = i
                end
            end
            nrows_new = min(nrows_accepted, nrows - n)

            # Columns go in outer loop because of column major indexing
            @inbounds @simd for j=1:ncolumns
                # Optimized inner loop
                for i = 1:nrows_new
                    resampled[n+i,j] = data[buffer[i],j]
                end
            end

            # Keep track of current filled rows
            n += nrows_new
        end

        return resampled
    end

    """
    ```julia
    randsample(data::Array, nrows, [p])
    ```
    Bootstrap resample (without uncertainty) a `data` array to length `nrows`.
    Optionally provide weights `p` either as a vector (one-weight-per-sample) or scalar.
    """
    function randsample(data::AbstractArray, nrows::Integer, p=min(0.2,nrows/size(data,1));
            rng::AbstractRNG=MersenneTwister(),
            buffer::Vector{Int}=Array{Int}(undef,size(data,1))
        )
        resampled = Array{eltype(data)}(undef,nrows,size(data,2))
        return randsample!(resampled, data, nrows, p, rng, buffer)
    end
    # Second method for randsample that takes a dictionary as input
    """
    ```julia
    randsample(dataset::Dict, nrows, [elements], [p])
    ```
    Bootstrap resample (without uncertainty) a `dataset` dict to length `nrows`.
    Optionally provide weights `p` either as a vector (one-weight-per-sample) or scalar.
    """
    function randsample(dataset::Dict, nrows::Integer, elements=in["elements"], p=min(0.2,nrows/length(in[elements[1]])))
        data = unelementify(dataset, elements, floatout=true)
        sdata = randsample(data, nrows, p)
        return elementify(sdata, elements, skipstart=0)
    end
    export randsample

## --- Bin a dataset by a given independent variable

    """
    ```julia
    (bincenters, N) = bincounts(x::AbstractArray, xmin::Number, xmax::Number, nbins::Integer)
    ```
    Tally the number of samples that fall into each of `nbins` equally spaced `x`
    bins between `xmin` and `xmax`, aligned with bin edges as
    `xmin:(xmax-xmin)/nbins:xmax`
    """
    function bincounts(x::AbstractArray, xmin::Number, xmax::Number, nbins::Integer)
        # Tally the number of samples (either resampled or corrected/original) that fall into each bin
        binwidth = (xmax-xmin)/nbins
        bincenters = (xmin+binwidth/2):binwidth:(xmax-binwidth/2)

        # Calculate index from x value
        scalefactor = nbins / (xmax - xmin)
        index = (x .- xmin) .* scalefactor

        # Add up the results
        N = fill(0,nbins)
        @inbounds for i in index
            if 0 < i <= nbins
                N[ceil(Int, i)] += 1
            end
        end
        return bincenters, N
    end
    export bincounts

    """
    ```julia
    (c,m,e) = binmeans(x, y, xmin, xmax, nbins, [weight]; resamplingratio::Number=1)
    ```
    The means (ignoring NaNs) of `y` values binned by `x`, into each of `nbins`
    equally spaced `x` bins between `xmin` and `xmax`, returning bincenters,
    means, and standard errors of the mean.

    To calculate binned medians only (without uncertainties), see `nanmean`
    """
    function binmeans(x::AbstractArray, y::AbstractArray, xmin::Number, xmax::Number, nbins::Integer; resamplingratio::Number=1)
        binwidth = (xmax-xmin)/nbins
        bincenters = (xmin+binwidth/2):binwidth:(xmax-binwidth/2)

        # Calculate index from x value
        scalefactor = nbins / (xmax - xmin)
        index_float = (x .- xmin) .* scalefactor

        # Calculate the nanmeans and nansems
        N = fill(0,nbins)
        mu = fill(0.0,nbins)
        s2 = fill(0.0,nbins)
        for i = 1:length(x)
            if (0 < index_float[i] <= nbins) && !isnan(y[i])
                index = ceil(Int, index_float[i])
                N[index] += 1
                delta = y[i] - mu[index]
                mu[index] += delta / N[index]
                s2[index] += delta * (y[i] - mu[index])
            end
        end

        # Calculate standard errors
        se = Array{Float64}(undef, nbins)
        for i=1:nbins
            if N[i] > 0
                s2[i] /= N[i] - 1 # Variance
                se[i] = sqrt(s2[i] * resamplingratio / N[i])
            else
                mu[i] = NaN
                se[i] = NaN
            end
        end
        return bincenters, mu, se
    end
    function binmeans(x::AbstractArray, y::AbstractArray, min::Number, max::Number, nbins::Integer, weight::AbstractArray; resamplingratio::Number=1)
        binwidth = (max-min)/nbins
        binedges = range(min,max,length=nbins+1)
        bincenters = (min+binwidth/2):binwidth:(max-binwidth/2)

        means = Array{Float64}(undef,nbins)
        errors = Array{Float64}(undef,nbins)
        for i = 1:nbins
            t = (binedges[i] .< x .<= binedges[i+1]) .& (y.==y)
            w = ProbabilityWeights(weight[t])
            means[i] = mean(y[t], w)
            errors[i] = std(y[t], w, corrected=true) * sqrt(resamplingratio) / sqrt(count(t))
        end

        return bincenters, means, errors
    end
    export binmeans

    """
    ```julia
    (c,m,e) = binmedians(x, y, xmin, xmax, nbins; resamplingratio::Number=1)
    ```

    The medians (ignoring NaNs) of `y` values binned by `x`, into each of `nbins`
    equally spaced `x` bins between `xmin` and `xmax`, returning bincenters, medians,
    and equivalent standard errors of the mean (1.4828 * median abolute deviation)

    To calculate binned medians only (without uncertainties), see `nanmedian`
    """
    function binmedians(x::AbstractArray, y::AbstractArray, min::Number, max::Number, nbins::Integer; resamplingratio::Number=1)
        binwidth = (max-min)/nbins
        binedges = range(min,max,length=nbins+1)
        bincenters = (min+binwidth/2):binwidth:(max-binwidth/2)

        medians = Array{Float64}(undef,nbins)
        errors = Array{Float64}(undef,nbins)
        for i = 1:nbins
            t = (binedges[i] .< x .<= binedges[i+1]) .& (y.==y)
            medians[i] = median(y[t])
            errors[i] = 1.4826 * nanmad(y[t]) * sqrt(resamplingratio / count(t))
        end

        return bincenters, medians, errors
    end
    export binmedians

## --- Bin bootstrap resampled data

    """
    ```julia
    bin_bsr([f!::Function], x::Vector, y::VecOrMat, xmin, xmax, nbins, [w];
        \tx_sigma = zeros(size(x)),
        \ty_sigma = zeros(size(y)),
        \tnresamplings = 1000,
        \tsem = :sigma,
        \tp = 0.2
    )
    ```
    Returns the bincenters `c`, means or medians `m`, and uncertainties of the
    mean or median for a variable `y` binned by independent variable `x` into
    `nbins` equal bins between `xmin` and `xmax`, after `nresamplings` boostrap
    resamplings with acceptance probability `p`.

    If a 2-d array (matrix) of `y` values is provided, each column will be treated
    as a separate variable, means and uncertainties will be returned column-wise.

    Optional keyword arguments and defaults:

        x_sigma = zeros(size(x))

    A vector representing the uncertainty (standard deviation) of each x value

        y_sigma = zeros(size(y))

    A vector representing the uncertainty (standard deviation) of each y value

        nresamplings = 1000

    The number of resamplings to conduct

        sem = :sigma

    Format of the uncertainty estimate of the distribution of the mean.
    If `:sigma` is chosen, a tuple of three vectors `(c, m, e)` will be returned,
    where `e` is the standard error of the mean.
    If `:CI` or `:pctile` is chosen, a tuple of four vectors `(c, m, el, eu)`
    will be returned, where `el` and `eu` are the lower and upper bounds of the 95%
    confidence interval.

        p = 0.2

    Resampling probabilities, either as a scalar or a vector of the same length as `x`

    # Examples:
    ```julia
    (c,m,e) = bin_bsr(nanbinmedian!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:sigma)
    ```
    ```julia
    (c,m,el,eu) = bin_bsr(nanbinmean!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:pctile)
    ```
    """
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractVector, xmin, xmax, nbins::Integer;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(y)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y)
        sigma = hcat(x_sigma, y_sigma)
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins)
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i), N, view(dbs,:,1), view(dbs,:,2), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        if sem == :sigma
            # Standard deviation of means (sem)
            e = nanstd(means,dim=2)
            return c, m, e
        elseif sem == :credibleinterval || sem == :CI || sem == :pctile
            # Lower bound of central 95% CI of means
            el = m .- nanpctile(means,2.5,dim=2)
            # Upper bound of central 95% CI of means
            eu = nanpctile(means,97.5,dim=2) .- m
            return c, m, el, eu
        else
            return c, means
        end
    end
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractMatrix, xmin, xmax, nbins::Integer;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(y)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y)
        sigma = hcat(x_sigma, y_sigma)
        dtype = float(eltype(data))
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{dtype}(undef, nrows, ncols)
        means = Array{dtype}(undef, nbins, nresamplings, size(y,2))
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins, size(y,2))
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i,:), N, view(dbs,:,1), view(dbs,:,2:1+size(y,2)), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        if sem == :sigma
            m = Array{dtype}(undef, nbins, size(y,2))
            e = Array{dtype}(undef, nbins, size(y,2))
            for j = 1:size(y,2)
                m[:,j] .= nanmean(view(means,:,:,j),dim=2) # Mean-of-means
                e[:,j] .= nanstd(view(means,:,:,j),dim=2) # Standard deviation of means (sem)
            end
            return c, m, e
        elseif sem == :credibleinterval || sem == :CI || sem == :pctile
            m = Array{dtype}(undef, nbins, size(y,2))
            el = Array{dtype}(undef, nbins, size(y,2))
            eu = Array{dtype}(undef, nbins, size(y,2))
            for j = 1:size(y,2)
                m[:,j] .= nanmean(view(means,:,:,j),dim=2) # Mean-of-means
                el[:,j] .= m[:,j] .- nanpctile(view(means,:,:,j),dim=2)
                eu[:,j] .= nanpctile(view(means,:,:,j),dim=2) .- m[:,j]
            end
            return c, m, el, eu
        else
            return c, means
        end
    end
    function bin_bsr(f!::Function, x::AbstractVector, y::AbstractVector, xmin, xmax, nbins::Integer, w::AbstractVector;
            x_sigma = zeros(size(x)),
            y_sigma = zeros(size(x)),
            nresamplings = 1000,
            sem = :credibleinterval,
            p = 0.2
        )

        data = hcat(x, y, w)
        sigma = hcat(x_sigma, y_sigma, zeros(size(w)));

        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        means = Array{Float64}(undef, nbins, nresamplings)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins)
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            f!(view(means,:,i), N, view(dbs,:,1), view(dbs,:,2), view(dbs,:,3), xmin, xmax, nbins)
        end

        # Return summary of results
        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        if sem == :sigma
            # Standard deviation of means (sem)
            e = nanstd(means,dim=2)
            return c, m, e
        elseif sem == :credibleinterval || sem == :CI || sem == :pctile
            # Lower bound of central 95% CI of means
            el = m .- nanpctile(means,2.5,dim=2)
            # Upper bound of central 95% CI of means
            eu = nanpctile(means,97.5,dim=2) .- m
            return c, m, el, eu
        else
            return c, means
        end
    end
    bin_bsr(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:sigma, p=0.2) =
        bin_bsr(nanbinmean!,x,y,xmin,xmax,nbins,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    bin_bsr(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer, w::AbstractVector; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:sigma, p=0.2) =
        bin_bsr(nanbinmean!,x,y,xmin,xmax,nbins,w,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    export bin_bsr

    bin_bsr_means(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr(nanbinmean!,x,y,xmin,xmax,nbins,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    bin_bsr_means(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer, w::AbstractVector; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr(nanbinmean!,x,y,xmin,xmax,nbins,w,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    export bin_bsr_means

    bin_bsr_medians(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr(nanbinmedian!,x,y,xmin,xmax,nbins,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    bin_bsr_medians(x::AbstractVector, y::AbstractVecOrMat, xmin, xmax, nbins::Integer, w::AbstractVector; x_sigma=zeros(size(x)), y_sigma=zeros(size(y)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr(nanbinmedian!,x,y,xmin,xmax,nbins,w,x_sigma=x_sigma,y_sigma=y_sigma,nresamplings=nresamplings,sem=sem,p=p)
    export bin_bsr_medians

    """
    ```julia
    (c, m, el, eu) = bin_bsr_ratios([f!::Function], x::Vector, num::Vector, denom::Vector, xmin, xmax, nbins, [w];
        \tx_sigma = zeros(size(x)),
        \tnum_sigma = zeros(size(num)),
        \tdenom_sigma = zeros(size(denom)),
        \tnresamplings = 1000,
        \tp::Union{Number,Vector} = 0.2
    )
    ```

    Returns the bincenters `c`, means `m`, as well as upper (`el`) and lower (`eu`) 95% CIs of the mean
    for a ratio `num`/`den` binned by `x` into `nbins` equal bins between `xmin` and `xmax`,
    after `nresamplings` boostrap resamplings with acceptance probability `p`.
    """
    function bin_bsr_ratios(f!::Function, x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer;
            x_sigma::AbstractVector=zeros(size(x)),
            num_sigma::AbstractVector=zeros(size(num)),
            denom_sigma::AbstractVector=zeros(size(denom)),
            nresamplings=1000,
            p::Union{Number,AbstractVector}=0.2
        )

        data = hcat(x, num, denom)
        sigma = hcat(x_sigma, num_sigma, denom_sigma)
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        fractions = Array{Float64}(undef, nrows)
        fraction_means = Array{Float64}(undef, nbins)
        rng = MersenneTwister()
        N = Array{Int}(undef, nbins) # Array of bin counts -- Not used but preallocated for speed
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            @views @avx @. fractions = dbs[:,2] / (dbs[:,2] + dbs[:,3])
            f!(fraction_means, N, view(dbs,:,1), fractions, xmin, xmax, nbins)
            @. means[:,i] = fraction_means / (1 - fraction_means)
        end

        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        el = m .- nanpctile(means,2.5,dim=2) # Lower bound of central 95% CI
        eu = nanpctile(means,97.5,dim=2) .- m # Upper bound of central 95% CI

        return c, m, el, eu
    end
    function bin_bsr_ratios(f!::Function, x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer, w::AbstractVector;
            x_sigma::AbstractVector=zeros(size(x)),
            num_sigma::AbstractVector=zeros(size(num)),
            denom_sigma::AbstractVector=zeros(size(denom)),
            nresamplings=1000,
            p::Union{Number,AbstractVector}=0.2
        )

        data = hcat(x, num, denom, w)
        sigma = hcat(x_sigma, num_sigma, denom_sigma, zeros(size(w)))
        binwidth = (xmax-xmin)/nbins
        nrows = size(data,1)
        ncols = size(data,2)

        # Preallocate
        dbs = Array{Float64}(undef, nrows, ncols)
        index = Array{Int}(undef, nrows) # Must be preallocated even if we don't want it later
        means = Array{Float64}(undef, nbins, nresamplings)
        fractions = Array{Float64}(undef, nrows)
        fraction_means = Array{Float64}(undef, nbins)
        rng = MersenneTwister()
        W = Array{Float64}(undef, nbins) # Array of bin weights -- Not used but preallocated for speed
        # Resample
        for i=1:nresamplings
            bsr!(dbs, index, data, sigma, p, rng=rng) # Boostrap Resampling
            @views @avx @. fractions = dbs[:,2] / (dbs[:,2] + dbs[:,3])
            f!(fraction_means, W, view(dbs,:,1), fractions, view(dbs,:,4), xmin, xmax, nbins)
            @. means[:,i] = fraction_means / (1 - fraction_means)
        end

        c = (xmin+binwidth/2):binwidth:(xmax-binwidth/2) # Bin centers
        m = nanmean(means,dim=2) # Mean-of-means
        el = m .- nanpctile(means,2.5,dim=2) # Lower bound of central 95% CI
        eu = nanpctile(means,97.5,dim=2) .- m # Upper bound of central 95% CI

        return c, m, el, eu
    end
    bin_bsr_ratios(x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer; x_sigma=zeros(size(x)), num_sigma=zeros(size(num)), denom_sigma=zeros(size(denom)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr_ratios(nanbinmean!,x,num,denom,xmin,xmax,nbins,x_sigma=x_sigma,num_sigma=num_sigma,denom_sigma=denom_sigma,nresamplings=nresamplings,p=p)
    bin_bsr_ratios(x::AbstractVector, num::AbstractVector, denom::AbstractVector, xmin, xmax, nbins::Integer, w::AbstractVector; x_sigma=zeros(size(x)), num_sigma=zeros(size(num)), denom_sigma=zeros(size(denom)), nresamplings=1000, sem=:pctile, p=0.2) =
        bin_bsr_ratios(nanbinmean!,x,num,denom,xmin,xmax,nbins,w,x_sigma=x_sigma,num_sigma=num_sigma,denom_sigma=denom_sigma,nresamplings=nresamplings,p=p)
    export bin_bsr_ratios

## --- Quick Monte Carlo binning/interpolation functions

    function mcfit(x::AbstractVector, σx::AbstractVector, y::AbstractVector, σy::AbstractVector,
            xmin::Number, xmax::Number, nbins::Integer=10;
            binwidth::Number=(xmax-xmin)/(nbins-1),
            minrows::Number=100000
        )
        # Run a simplified Monte Carlo fit with nbins of witdth binwidth between xmin and xmax

        # Remove missing data
        hasdata = .!(isnan.(x) .| isnan.(y))
        x′ = x[hasdata]
        y′ = y[hasdata]
        σx′ = σx[hasdata]
        σy′ = σy[hasdata]

        # Fill in variances where not provided explicitly
        σx′[isnan.(σx′)] .= nanstd(x′)
        σy′[isnan.(σy′)] .= nanstd(y′)

        # Increase x uncertainty if x sampling is sparse
        xsorted = [xmin; sort(x′); xmax]
        minerr = maximum(xsorted[2:end] - xsorted[1:end-1]) / 2
        σx′[σx′ .< minerr] .= minerr

        # Bin centers
        c = xmin:(xmax-xmin)/(nbins-1):xmax
        halfwidth = binwidth / 2

        # Run the Monte Carlo
        N = fill(0, nbins)
        m = fill(zero(float(eltype(y′))), nbins)
        xresampled = similar(x′, float(eltype(x′)))
        yresampled = similar(y′, float(eltype(y′)))
        @inbounds for n = 1:ceil(Int, minrows/length(x′))
            randn!(xresampled)
            xresampled .= x′ .+ σx′ .* xresampled
            randn!(yresampled)
            yresampled .= y′ .+ σy′ .* yresampled
            for j = 1:nbins
                l = (c[j] - halfwidth)
                u = (c[j] + halfwidth)
                for i ∈ eachindex(xresampled)
                    if l < xresampled[i] <= u
                        m[j] += yresampled[i]
                        N[j] += 1
                    end
                end
            end
        end
        m ./= N

        return c, m
    end
    export mcfit

## --- Downsample an image / array

    function downsample(A::AbstractArray, factor::Integer, jfactor=factor::Integer)
        if ndims(A)==2
            rows = size(A,1) ÷ factor
            cols = size(A,2) ÷ jfactor

            result = Array{eltype(A)}(undef,rows,cols)
            @avx for i=1:rows
                for j=1:cols
                    result[i,j]=A[i*factor,j*jfactor]
                end
            end
        else
            result = A[factor:factor:end]
        end
        return result
    end
    export downsample

## --- Spatiotemporal sample weighting

    """
    ```julia
    k = invweight(lat::AbstractArray, lon::AbstractArray, age::AbstractArray;
        \tlp::Number=2,
        \tspatialscale=1.8,
        \tagescale=38.0
    )
    ```

    Find the inverse weights `k` (proportional to spatiotemporal sample density) for
    a set of geological samples with specified latitude (`lat`), logitude (`lon`),
    and `age` (of crystallization, deposition, etc.).

    The default `spatialscale` and `agescale` are taken from Keller and Schoene 2012.
    However, alternative scalings can be supplied. If an array is supplied for either
    `spatialscale`, `agescale`, or both, a 3-d matrix of `k` values will be returned,
    with dimensions length(`spatialscale`)*length(`agescale`)*nrows.
    """
    function invweight(lat::AbstractArray, lon::AbstractArray, age::AbstractArray;
        lp::Number=2, spatialscale=1.8, agescale=38.0)

        # Check if there is lat, lon, and age data
        nodata = vec(isnan.(lat) .| isnan.(lon) .| isnan.(age))

        # Convert lat and lon to radians
        latr = vec(lat/180*pi)
        lonr = vec(lon/180*pi)
        spatialscalr = spatialscale/180*pi

        # Precalculate some sines and cosines
        latsin = sin.(latr)
        latcos = cos.(latr)

        # Allocate and fill ks
        if isa(spatialscale, Number) && isa(agescale, Number)
            k = Array{Float64}(undef,length(lat))
            @showprogress 1 "Calculating weights: " for i=1:length(lat)
                if nodata[i] # If there is missing data, set k=inf for weight=0
                    k[i] = Inf
                else # Otherwise, calculate weight
                    k[i] = nansum( @avx @. 1.0 / ( (acos(min( latsin[i] * latsin + latcos[i] * latcos * cos(lonr[i] - lonr), 1.0 ))/spatialscalr)^lp + 1.0) + 1.0 / ((abs(age[i] - age)/agescale)^lp + 1.0) )
                end
            end
        else
            k = Array{Float64}(undef,length(spatialscale),length(agescale),length(lat))
            spatialdistr = similar(lat)
            @showprogress 1 "Calculating weights: " for i=1:length(lat)
                if nodata[i] # If there is missing data, set k=inf for weight=0
                    k[:,:,i] .= Inf
                else # Otherwise, calculate weight
                    @avx @. spatialdistr = acos(min( latsin[i] * latsin + latcos[i] * latcos * cos(lonr[i] - lonr), 1.0 ))
                    Threads.@threads for g = 1:length(spatialscale)
                        for h = 1:length(agescale)
                            k[g,h,i] = nansum( @avx @. 1.0 / ( (spatialdistr/spatialscalr[g])^lp + 1.0) + 1.0 / ((abs(age[i] - age)/agescale[h])^lp + 1.0) )
                        end
                    end
                end
            end
        end
        return k
    end
    export invweight

    """
    ```julia
    k = invweight_location(lat::AbstractArray, lon::AbstractArray;
        \tlp::Number=2,
        \tspatialscale::Number=1.8
    )
    ```

    Find the inverse weights `k` (proportional to spatial sample density) for
    a set of geological samples with specified latitude (`lat`), and logitude (`lon`).
    """
    function invweight_location(lat::AbstractArray, lon::AbstractArray;
        lp::Number=2, spatialscale::Number=1.8)

        # Check if there is lat, lon data
        nodata = vec(isnan.(lat) .| isnan.(lon))

        # Convert lat and lon to radians
        latr = vec(lat/180*pi)
        lonr = vec(lon/180*pi)
        spatialscalr = spatialscale/180*pi

        # Precalculate some sines and cosines
        latsin = sin.(latr)
        latcos = cos.(latr)

        # Allocate and fill ks
        k = Array{Float64}(undef,length(lat))
        @showprogress 1 "Calculating weights: " for i=1:length(lat)
            if nodata[i] # If there is missing data, set k=inf for weight=0
                k[i] = Inf
            else # Otherwise, calculate weight
                k[i] = nansum( @avx @. 1.0 / ( (acos(min( latsin[i] * latsin + latcos[i] * latcos * cos(lonr[i] - lonr), 1.0 ))/spatialscalr)^lp + 1.0) )
            end
        end
        return k
    end
    export invweight_location

    """
    ```julia
    k = invweight(nums::AbstractArray, scale::Number; lp=2)
    ```

    Find the inverse weights for a single array `nums` for a given `scale`, and
    exponent `lp` (default lp = 2).

    Returns an array k where k[i] is the "inverse weight" for element i of the
    input array.
    """
    function invweight(nums::AbstractArray, scale::Number; lp=2)
        # Check if there is data
        nodata = isnan.(nums)

        k = Array{Float64}(undef,length(nums))
        @showprogress 1 "Calculating weights: " for i=1:length(nums)
            if nodata[i] # If there is missing data, set k=inf for weight=0
                k[i] = Inf
            else # Otherwise, calculate weight
                k[i] = nansum( @avx @. 1.0 / ( (abs(nums[i] - nums)/scale)^lp + 1.0) )
            end
        end
        return k
    end
    export invweight

    """
    ```julia
    k = invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
    ```

    Find the inverse weights `k` (proportional to temporal sample density) for
    a set of geological samples with specified `age` (of crystallization, deposition, etc.).
    """
    function invweight_age(age::AbstractArray; lp::Number=2, agescale::Number=38.0)
        return invweight(age, agescale, lp=lp)
    end
    export invweight_age

## --- End of file


text:
null = FixedTimeZone("", 10800)
fixed = FixedTimeZone("UTC+01:00")
est = FixedTimeZone("EST", -18000)
warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)
apia = resolve("Pacific/Apia", tzdata["australasia"]...)
honolulu = resolve("Pacific/Honolulu", tzdata["northamerica"]...)  # Uses cutoff
ulyanovsk = resolve("Europe/Ulyanovsk", tzdata["europe"]...)  # No named abbreviations
dt = DateTime(1942,12,25,1,23,45)

buffer = IOBuffer()

# TimeZones as a string
@test string(null) == "UTC+03:00"
@test string(fixed) == "UTC+01:00"
@test string(est) == "EST"
@test string(warsaw) == "Europe/Warsaw"
@test string(apia) == "Pacific/Apia"
@test string(honolulu) == "Pacific/Honolulu"
@test string(ulyanovsk) == "Europe/Ulyanovsk"

showcompact(buffer, null)
@test takebuf_string(buffer) == "UTC+03:00"
showcompact(buffer, fixed)
@test takebuf_string(buffer) == "UTC+01:00"
showcompact(buffer, est)
@test takebuf_string(buffer) == "EST"
showcompact(buffer, warsaw)
@test takebuf_string(buffer) == "Europe/Warsaw"
showcompact(buffer, apia)
@test takebuf_string(buffer) == "Pacific/Apia"
showcompact(buffer, honolulu)
@test takebuf_string(buffer) == "Pacific/Honolulu"
showcompact(buffer, ulyanovsk)
@test takebuf_string(buffer) == "Europe/Ulyanovsk"

show(buffer, null)
@test takebuf_string(buffer) == "UTC+03:00"
show(buffer, fixed)
@test takebuf_string(buffer) == "UTC+01:00"
show(buffer, est)
@test takebuf_string(buffer) == "EST (UTC-5)"
show(buffer, warsaw)
@test takebuf_string(buffer) == "Europe/Warsaw (UTC+1/UTC+2)"
show(buffer, apia)
@test takebuf_string(buffer) == "Pacific/Apia (UTC+13/UTC+14)"
show(buffer, honolulu)
@test takebuf_string(buffer) == "Pacific/Honolulu (UTC-10)"
show(buffer, ulyanovsk)
@test takebuf_string(buffer) == "Europe/Ulyanovsk (UTC+4)"

# UTC and GMT are special cases
show(buffer, FixedTimeZone("UTC"))
@test takebuf_string(buffer) == "UTC"
show(buffer, FixedTimeZone("GMT", 0))
@test takebuf_string(buffer) == "GMT"
show(buffer, FixedTimeZone("FOO", 0))
@test takebuf_string(buffer) == "FOO (UTC+0)"

# ZonedDateTime as a string
zdt = ZonedDateTime(dt, warsaw)
@test string(zdt) == "1942-12-25T01:23:45+01:00"

show(buffer, zdt)
@test takebuf_string(buffer) == "1942-12-25T01:23:45+01:00"

# TimeZone parsing
df = Dates.DateFormat("z")
@test !isempty(df.slots)  # Ensure that 'z' slot character is recognized (issue #24)
@test Dates.parse("+0100", df) == Any[FixedTimeZone("+01:00")]
@test_throws ArgumentError Dates.parse("EST", df)
@test_throws ArgumentError Dates.parse("UTC", df)
@test_throws ArgumentError Dates.parse("Europe/Warsaw", df)

df = Dates.DateFormat("Z")
@test !isempty(df.slots)  # Ensure that 'Z' slot character is recognized (issue #24)
@test_throws ArgumentError Dates.parse("+0100", df)
@test_throws ArgumentError Dates.parse("EST", df)
@test Dates.parse("UTC", df) == Any[FixedTimeZone("UTC")]
@test Dates.parse("Europe/Warsaw", df) == Any[warsaw]

# ZonedDateTime parsing.
# Note: uses compiled time zone information. If these tests are failing try to rebuild
# the TimeZones package.
@test ZonedDateTime("1942-12-25T01:23:45.0+01:00") == ZonedDateTime(dt, fixed)
@test ZonedDateTime("1942-12-25T01:23:45+0100", "yyyy-mm-ddTHH:MM:SSzzz") == ZonedDateTime(dt, fixed)
@test ZonedDateTime("1942-12-25T01:23:45 Europe/Warsaw", "yyyy-mm-ddTHH:MM:SS ZZZ") == ZonedDateTime(dt, warsaw)

# Note: CET here represents the FixedTimeZone used in Europe/Warsaw and not the
# VariableTimeZone CET.
@test_throws ArgumentError ZonedDateTime("1942-12-25T01:23:45 CET", "yyyy-mm-ddTHH:MM:SS ZZZ")

# Creating a ZonedDateTime requires a TimeZone to be present.
@test_throws ArgumentError ZonedDateTime("1942-12-25T01:23:45", "yyyy-mm-ddTHH:MM:SSzzz")

# ZonedDateTime formatting
f = "yyyy/m/d H:M:S ZZZ"
@test Dates.format(ZonedDateTime(dt, fixed), f) == "1942/12/25 1:23:45 UTC+01:00"
@test Dates.format(ZonedDateTime(dt, warsaw), f) == "1942/12/25 1:23:45 CET"
@test Dates.format(ZonedDateTime(dt, ulyanovsk), f) == "1942/12/25 1:23:45 UTC+04:00"

f = "yyyy/m/d H:M:S zzz"
@test Dates.format(ZonedDateTime(dt, fixed), f) == "1942/12/25 1:23:45 +01:00"
@test Dates.format(ZonedDateTime(dt, warsaw), f) == "1942/12/25 1:23:45 +01:00"
@test Dates.format(ZonedDateTime(dt, ulyanovsk), f) == "1942/12/25 1:23:45 +04:00"

# The "Z" slot displays the time zone abbreviation for VariableTimeZones. It is fine to use
# the abbreviation for display purposes but not fine for parsing. This means that we
# currently cannot parse all strings produced by format.
f = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS ZZZ")
@test_throws ArgumentError Dates.parse(Dates.format(ZonedDateTime(dt, warsaw), f), f)