text:
#==============================================================================#
# preGS.jl
#==============================================================================#

# Author:    Austin Putz
# Created:   Sept 27, 2017
# Modified:  August 31, 2018
# License:   GPLv2

# set version number
preGS_VERSION = "0.0.3"

# developed in Julia version __
JULIA_VERSION = "0.7.0"

#==============================================================================#
# Usage
#==============================================================================#

# Usage:
#			$ julia preGS.jl input_file.jl

#==============================================================================#
# Will determine by itself
#==============================================================================#

# 1. Type of storage you have based on the first few lines it will read in

#==============================================================================#
# Laundry list of shit to do!
#==============================================================================#

# Stars link to importance, even if not in order

# *** Compile into a standalone program (people won't have to download Julia)

# ** Take missing values

#==============================================================================#
#
# PRINT HEADER
#
#==============================================================================#

# Images are https://www.text-image.com
# OR: http://patorjk.com/software/taag/#p=display&f=Graffiti&t=Type%20Something%20
# Simply select an image and upload it and it will convert it to ASCII Text

println("")
println("--------------------------------------------------------------------------------")
println("--                                                                            --")
println("--                                           `              -:.//-::.         --")
println("--                                   `.`:+/:o+//o+:.::.--oy+ss-o+/yss         --")
println("--       /:`                      .-:yss+so-:s:.oo:oysos+/:: / /oy//+         --")
println("--       sho+oo+/`             .-:yys-:/  +  +  + .oo:ossoso//oso+sss         --")
println("--       dmshysshy://:.      ..yys-`/` /` /` + ./oso/+/:o+ysss+/s/.           --")
println("--       -d-:oyyyoysoos/`..`-yso`.: `/  / -- + +sso+///::/+oyys+:             --")
println("--        +y   +/`+sssosso++/+.-- `:``: -- / +.ss+.         /ssyo-            --")
println("--          y/  `s`  -- /yyssso///:.`:.`: /`/`+ss-           .soyy:           --")
println("--            .h.  //  .: .:o+-+yso+++:./-./`/--os/             `::-          --")
println("--        `  /s   s`  :yhs/   :++/sho+//--/:-os+.                             --")
println("--        o   s:  -/`+oyo-        -osooy++/:+++:                              --")
println("--        y:  `y`  o/dyo.             /o++o+//:`                              --")
println("--        -s`  :o yhyh.                  .-.-`                                --")
println("--         /+  /hyhs/.                   _________     _ __                   --")
println("--           y+/ddyy.   ____   _______  / ___/ __/    ( ) /                   --")
println("--       `yhyyh-.      /  _ \\/  __/ -_) (_ / \\ \\_    / / /                    --")
println("--       hhdyo-       /  .__/ _/  \\__/\\___/____(_)_ / /_/                     --")
println("--       h/          /__/                        |___/                        --")
println("--                                                                            --")
println("--                                                                            --")
println("--                preGS.jl - A genomic pre-processor based on                 --")
println("--                        Ignacio Aguilar's 'preGSf90'                        --")
println("--                                                                            --")
println("--                               Austin Putz                                  --")
println("--                          Iowa State University                             --")
println("--                             License: GPLv2                                 --")
println("--                   Will be available on GitHub soon...                      --")
println("--                                                                            --")
println("--------------------------------------------------------------------------------")

#                      ________________ ______________
#  ______________________  ____/_  ___/ ______(_)__  /
#  ___  __ \_  ___/  _ \  / __ _____ \  _____  /__  /
#  __  /_/ /  /   /  __/ /_/ / ____/ /______  / _  /
#  _  .___//_/    \___/\____/  /____/_(_)__  /  /_/
#  /_/                                  /___/

#                 _________    _ __
#   ___  _______ / ___/ __/   (_) /
#  / _ \/ __/ -_) (_ /\ \_   / / /
# / .__/_/  \__/\___/___(_)_/ /_/
#/_/                     |___/

#println("||--------------------------------------------------------------------------------||")

# the following was taken from the website: www.network-science.de/ascii/
# using the 'speed' Font

# Print Header
#println("|                                       _________________                        |")
#println("|                 _____________________/   ____/    ____/                        |")
#println("|                 ___  __ \\_ / ___/  _ \\  /  ___ \\____ \\                         |")
#println("|                 __  /_/ / / /  /  __/  /_/ /  ____/  /                         |")
#println("|                 _  .___/ /_/   \\___/\\_____/  /______/                          |")
#println("|                 /_/                                                            |")

#==============================================================================#
# Packages and Functions
#==============================================================================#

# load Packages
using DataFrames
using UnicodePlots
using ProgressMeter
using CSV
using Printf
using Distributed
using Statistics
using DelimitedFiles
using LinearAlgebra
#using PyPlots

#==============================================================================#
# Read in input files to change defaults
#==============================================================================#

# set directory
current_directory = pwd()

# source input file
parameter_file = ARGS[1]

# read in file with inputs
include(string(current_directory, "/", parameter_file))

#==============================================================================#
#
# Begin Main Body
#
#==============================================================================#

#==============================================================================#
# Read in Map File
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Read in Map File ----------")
println("--------------------------------------------------------------------------------\n")

if map_file == ""

	println("--- No map file specified ----")

else

println("---- Read in Map File ----")

# read in Map File (SNP name, Chromosome, Position)
@time user_map = CSV.read(map_file,
	delim      = map_file_delimiter,
	header     = ["SNP", "Chr", "Pos"],
	datarow    = map_file_start_row,
	types      = [String, Int64, Int64])

println("\n---- Head of Map File ----")

# head map file
println(head(user_map))

# describe the Map File
println("\n---- Summary of the Map File ---- ")
@printf "-- %-30s: %d \n" "Chromosome Minimum" minimum(user_map[:, :Chr])
@printf "-- %-30s: %d \n" "Chromosome Maximum" maximum(user_map[:, :Chr])
@printf "-- %-30s: %d \n" "Number of SNPs"     length(user_map[:,  :Chr])

# describe the map file
#println(describe(user_map[:, :Chr]))

end

#==============================================================================#
# Read in genotypes
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Read in Genotypes ----------")
println("--------------------------------------------------------------------------------\n")

#------------------------------------------------------------------------------#
# Dense Storage
#------------------------------------------------------------------------------#

# Begin Dense Extraction
if SNP_format == "dense"

println("---- Reading in Genotypes From Compact Format ----")

# read in Map File (SNP name, Chromosome, Position)
#@time geno = CSV.read(genotype_file, delim = genotype_file_delimiter, header = ["ID", "SNPs"], datarow = genotype_file_start_row, types = [String, String])
@time geno = CSV.read(genotype_file, delim = ' ', header = ["ID", "SNPs"], datarow = 1, types = [String, String])

# set number of records
n_IDs_orig = size(geno, 1)

# set number of SNPs
n_SNPs_orig = length(geno[1, 2])

# print number of SNPs
@printf "-- %-50s: %d \n" "Number of SNPs in Raw Genotype File" n_SNPs_orig

# sample SNPs
if n_sample_M == 0
	# set the number to sample equal to all SNPs
	x = 1:n_SNPs_orig
	n_SNPs_sample = n_SNPs_orig
else
	# sample a number of SNPs from the total
	x = unique(rand(1:n_SNPs_orig, n_sample_M))
	n_SNPs_sample = length(x)
end

# print number of SNPs and IDs
@printf "-- %-50s: %d \n" "Number of SNPs (or sample) from Genotype File" n_SNPs_sample
@printf "-- %-50s: %d \n" "Number of IDs" n_IDs_orig

# function to fill the M matrix
function fillM(n_IDs::Int, n_SNPs::Int, geno)

	# initiate M with the size
	M = zeros(n_IDs, n_SNPs);

	# show progress meter (took out because it doesn't work well with @parallel)
	#n = n_IDs*n_SNPs
	#p = Progress(n, 1)

	# split string of genotypes into the M matrix
	@distributed for i = 1:n_IDs
		@distributed for j = 1:n_SNPs
			# test if you are sampling genotypes
			if n_sample_M == 0
				M[i, j] = parse(Int64, geno[i, 2][j])
			else
				M[i, j] = parse(Int64, geno[i, 2][x[j]])
			end
			# use progress bar (doesn't work well with parallel computing so I removed
			#next!(p)
		end
	end

	# convert M to an Array
	M = Array(M)

	# return M
	return M

end

println("-- Loop to fill M Matrix --")
@time M = fillM(n_IDs_orig, n_SNPs_sample, geno)

# print messages of the number of rows and columns
@printf "-- %-50s: %d \n" "Number of rows in M"    size(M, 1)
@printf "-- %-50s: %d \n" "Number of columns in M" size(M, 2)

# calculate the number of SNPs
n_SNPs_current = size(M, 2);

end  # end of "if SNP_format == "dense"

#------------------------------------------------------------------------------#
# Plot M
#------------------------------------------------------------------------------#

println("-- Reshape M --")

# reshape M into a 1-D Vector to Plot!
@time Mvec = vec(reshape(M, (length(M), 1)))

# print histogram of elements of M
println("-- Plot of elements in M --")
println(histogram(Mvec, bins = 5, title = "Plot of the Elements in M (raw)"))
println("-- Done Getting M --")

#==============================================================================#
# Allele Frequencies
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Allele Frequencies (prior to cleaning) ----------")
println("--------------------------------------------------------------------------------\n")

# Function to calculate allele frequencies with missing values
function calcAlleleFreqs(M::Array)

	# initiate vector for allele frequencies
	allele_freqs = zeros(size(M, 2))

	# calculate allele frequencies
	@distributed for i = 1:size(M,2)

		# subset each vector
		current_SNP = vec(M[:, i])

		# subset to 0, 1, 2
		current_SNP = current_SNP[current_SNP .< 2.0001]
		current_SNP = current_SNP[current_SNP .> -0.0001]

		# take mean and divide by 2
		allele_freqs[i] = Statistics.mean(current_SNP) / 2

	end

	# return allele freqs
	return allele_freqs

end

# Calculate allele freqs
allele_freqs_before_clean = calcAlleleFreqs(M)

# histogram of allele frequencies
println("\n-- Plot of Allele Frequenceis --\n")
println(histogram(vec(allele_freqs_before_clean), bins=20, title = "Allele Frequencies (prior to cleaning)"))

#==============================================================================#
# Call Rates
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Call Rates ----------")
println("--------------------------------------------------------------------------------\n")

# Initalize Data Frame for IDs
data_IDs  = DataFrame(ID_Number = 1:n_IDs_orig, ID_Orig = geno[:, 1]);

# Initalize Data Frame for SNPs
data_SNPs = DataFrame(SNP_Number = 1:n_SNPs_current);

# add allele freqs
data_SNPs[:AlleleFreqsBeforeClean] = allele_freqs_before_clean

# set new column if map file exits
#(TO-DO)

# set up vectors to fill
ID_call_rate_orig  = zeros(n_IDs_orig);
SNP_call_rate_orig = zeros(n_SNPs_current);

# loop to check ID call rate
MissingMat = M .== genotype_missing_value

# sum to across row to get animal call rate
data_IDs[:CallRate]  = 1 .- (vec(sum(MissingMat, dims=2)) / size(MissingMat, 2))

# sum down row and divide to get SNP call rate
data_SNPs[:CallRate] = 1 .- (vec(sum(MissingMat, dims=1)) / size(MissingMat, 1))

# set true/false call rates
data_IDs[:CallRateTF]  = data_IDs[:CallRate]  .> call_rate_animals
data_SNPs[:CallRateTF] = data_SNPs[:CallRate] .> call_rate_SNPs

# Print
println(head(data_IDs))
println(head(data_SNPs))

# plot Call Rate IDs
println("\n---- Call Rate IDs ----")

# convert call rate to a vector to plot
IDs_call_rate_vec = vec(data_IDs[:, :CallRate])

# plot histogram of the Animal Call Rates
#println(histogram(IDs_call_rate_vec, bins=5, title="Call Rate - Animals"))
println(histogram(IDs_call_rate_vec, bins=10, title="Call Rate - Animals"))

# print stats
Printf.@printf "-- %-40s: %0.3f, %0.3f, %0.3f, %0.3f \n" "Minimum, Maximum, Average, Median" minimum(IDs_call_rate_vec) maximum(IDs_call_rate_vec) mean(IDs_call_rate_vec) median(IDs_call_rate_vec)

# print summary statistics for Animal Call Rate
describe(data_IDs[:, :CallRate])

# plot Call Rate IDs
println("\n---- Call Rate SNPs ----")
SNPs_call_rate_vec = vec(data_SNPs[:, :CallRate])

# plot histogram of the Animal Call Rates
println(histogram(SNPs_call_rate_vec, bins=10, title="Call Rate - SNPs"))

# print stats
Printf.@printf "-- %-40s: %0.3f, %0.3f, %0.3f, %0.3f \n" "Minimum, Maximum, Average, Median" minimum(SNPs_call_rate_vec) maximum(SNPs_call_rate_vec) mean(SNPs_call_rate_vec) median(SNPs_call_rate_vec)

# print summary statistics for Animal Call Rate
describe(data_SNPs[:, :CallRate])

#==============================================================================#
# Subset M matrix
#==============================================================================#

# Remove Rows for Call Rate on IDs
M = M[vec(data_IDs[:CallRateTF]), :];

# Remove Columns for Call Rate on SNPs
M = M[:, vec(data_SNPs[:CallRateTF])];

# Print size of M after cleaning for Call Rate
println("-- Size of M after Call Rates Removed --")
Printf.@printf "-- %-40s: %d, %d \n" "Size of new M matrix after subsetting for outliers" size(M,1) size(M,2)

#==============================================================================#
# Create Main Matrices
#==============================================================================#

println("\n\n\n--------------------------------------------------------------------------------")
println("---------- Create Matrices for G ----------")
println("--------------------------------------------------------------------------------\n")

# check for anything put 0/1/2 genotypes

# ToDo: need to write check for 0/1/2 genotypes

#------------------------------------------------------------------------------#
# Calculate column means of M
#------------------------------------------------------------------------------#

println("---- Column Means of M ----")

# get column means of M matrix
@time col_means = vec(Statistics.mean(M, dims=1))

#------------------------------------------------------------------------------#
# Calculate P
#------------------------------------------------------------------------------#

println("---- Calculate P ----")

# repeat column means (2*p) for P matrix
@time P = reshape(repeat(col_means, inner=size(M, 1)), (size(M,1), size(M,2)));

#------------------------------------------------------------------------------#
# Calculate Z
#------------------------------------------------------------------------------#

println("---- Calculate Z ----")

# subtract off column means
@time Z = M - P;

# calculate the size of G
Gdim = size(Z, 1);

#------------------------------------------------------------------------------#
# Calculate G
#------------------------------------------------------------------------------#

#----------------------------------------#
# Calculate G
#----------------------------------------#

println("---- Calculate G (unscaled) ----")

# calculate G
@time G = Z*Z'

println("---- Size of G ----")

# get size of G
@printf "-- %-40s: %d \n" "Number of rows/columns in G" size(G, 2)

#------------------------------------------------------------------------------#
# Allele Freqs & Scale Factor
#------------------------------------------------------------------------------#

println("---- Calculate Allele Frequencies (p & q) ----")

# get allele frequencies
@time p = vec(col_means ./ 2)  # column means already calculated
@time q = vec(1 .- p)

println("---- Calculate Scale Factor ----")

# set sum 2pq
@time sum2pq = 2*sum(p.*q)

println("Sum2pq is: $sum2pq")

#------------------------------------------------------------------------------#
# Scale G
#------------------------------------------------------------------------------#

println("---- Calculate G (scaled) ----")

# calculate G
@time G = Array(G / sum2pq)

# write G matrix
open("G.txt", "w") do io
	writedlm(io, G)
end

println("")

#==============================================================================#
# Plot components
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Plot Elements of G ----------")
println("--------------------------------------------------------------------------------\n")

#----------------------------------------#
# Diagonals of G
#----------------------------------------#

println("---- Print Diagonal Elements of G ----\n")
# get diags
@time Gdiags = vec(LinearAlgebra.diag(G))

# histogram of my diagonals
println(histogram(Gdiags, bins = 20, title = "G Matrix Diagonals"))

#----------------------------------------#
# Off-Diagonals of G
#----------------------------------------#

println("---- Print Off-Diagonal Elements of G ----\n")

function extractOffdiags(G::Array)

	# calculate the size of G
	sizeG = size(G, 1)

	# initiate off-diagonal Vector
	# offdiagsG = [];
	#offdiagsG = zeros((((1+sizeG)*sizeG)/2), 1);  # equation for upper triangular part
	#offdiagsG = zeros(Float64, convert(Float64, (((1+sizeG)*sizeG)/2), 1));
	offdiag_length = convert(BigInt, ((1+sizeG)*sizeG) / 2)
	println("\nOffdiagonal length is : $offdiag_length\n")
	#offdiagsG = ones(offdiag_length, 1)
	offdiagsG = Array{Float64,1}(undef, offdiag_length)

	# set count
	count = 1;

	# start progress bar
	n = length(offdiagsG)
	p = Progress(n, 1)

	# loop to fill off-diag vector
	for i = 1:sizeG
		for j = i+1:sizeG

			# add off-diag to list
			if i < j
				offdiagsG[count, 1] = G[i, j]
			end

			# advance count
			count = count + 1

			# advance progress bar
			next!(p)

		end
	end

	# convert to vector
	offdiagsG = vec(offdiagsG)

	# return off-diagonals
	return offdiagsG

end # END function

# extract the off-diagonals of G
@time offdiagsG = extractOffdiags(G)

# print offdiags of G
println(histogram(offdiagsG, bins = 20, title = "G Matrix Off-Diagonals"))

#----------------------------------------#
# Relationships of G
#----------------------------------------#

println("---- Print Relationships of G ----\n")

function extractRelationships(G::Array)

	# size of G
	sizeG = size(G, 1)

	# initalize the vector
	Grels = Float64[];

	# loop through to standardize the G matrix to relationships
	for i = 1:sizeG
		for j = i+1:sizeG

			# Skip if it's a diagonal or lower triangular
			if i == j || i > j
				continue
			end

			# fill G relationship matrix
			push!(Grels, (G[i, j] / sqrt( G[i,i] * G[j,j] )))

		end
	end

	# convert to vector
	Grels = vec(Grels)

	# return relationships
	return Grels

end # END Function

@time Grels = extractRelationships(G)

# print histogram of relationships
println(histogram(Grels, bins = 20, title = "G relationships"))

#==============================================================================#
# PCA Plot
#==============================================================================#

println("\n--------------------------------------------------------------------------------")
println("---------- Plot PCA's ----------")
println("--------------------------------------------------------------------------------\n")

println("---- Singular Value Decomposition for PCA's ----")
# Singular Value Decomposition
@time U, S, V = LinearAlgebra.svd(G)

# sum of s^2 values is the total variance
S2 = S.^2

# get each percent explained
S_pct_exp = vec((S2 ./ sum(S2)).*100)

# pull out the first few PCAs
S_1 = S_pct_exp[1]
S_2 = S_pct_exp[2]
S_3 = S_pct_exp[3]
S_4 = S_pct_exp[4]
S_5 = S_pct_exp[5]

# print how much they explain
@printf "\n-- %-40s: %0.2f, %0.2f, %0.2f, %0.2f, %0.2f \n " "Variances Explained of 1st 5 PCAs" S_1 S_2 S_3 S_4 S_5

# extract first 10 variance explained
S_pct_exp_first_10 = S_pct_exp[1:10]

# get vector list 1 to 10
first_10 = collect(1:10)

println("\n")

# plot
println(lineplot(first_10, S_pct_exp_first_10, title = "Scree Plot of First 10 PCAs (var exp)"))

println("\n\n")

# if you set only one PC, set to 2 to plot between 1 and 2
if n_PCs == 1
	n_PCs = 2
end

	# Loop through to make Plots
	for i = 1:n_PCs
		for j = 2:n_PCs

		if i == j || j < i || i == n_PCs
			continue
		end

		# plot the 1,2 principal components
		println(scatterplot(vec(U[:, i]), vec(U[:, j]),
			color = :red,
			title = "PC $i vs PC $j"))

		end
	end # end loop for PCA plots

# print space
println("\n\n")


text:
# Preprocessor
function EnterMainSourceFile(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_EnterMainSourceFile(x)
end

function EnterSourceFile(x::Preprocessor, id::FileID, loc::SourceLocation=SourceLocation())
    @check_ptrs x id
    return clang_Preprocessor_EnterSourceFile(x, id, loc)
end

function EndSourceFile(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_EndSourceFile(x)
end

function getHeaderSearchInfo(x::Preprocessor)
    @check_ptrs x
    return HeaderSearch(clang_Preprocessor_getHeaderSearchInfo(x))
end

function PrintStats(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_PrintStats(x)
end

function InitializeBuiltins(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_InitializeBuiltins(x)
end

function isIncrementalProcessingEnabled(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_isIncrementalProcessingEnabled(x)
end

function enableIncrementalProcessing(x::Preprocessor)
    @check_ptrs x
    return clang_Preprocessor_enableIncrementalProcessing(x)
end

function DumpToken(x::Preprocessor, tok::Token, flag=false)
    @check_ptrs x
    return clang_Preprocessor_DumpToken(x, tok, flag)
end

function DumpLocation(x::Preprocessor, loc::SourceLocation)
    @check_ptrs x
    return clang_Preprocessor_DumpLocation(x, loc)
end

ext: jl


text:
function initialrelease!(node::Int64, params::Parameters)
  for i in 1:6
    params.rdel[i, 1] = copy(params.del[i,1, node])
  end

  if 1 < node < params.noofnodes
    params.fvec[1,1,node] = params.f[1,1,node] - (
      params.sa[1,1,node]*params.del[1,1,node-1] + params.sc[1,1,node]*params.del[1,1,node+1] +
      params.sa[1,5,node]*params.del[5,1,node-1] + params.sc[1,5,node]*params.del[5,1,node+1])

    params.fvec[2,1,node] = params.f[2,1,node] - (
      params.sa[2,2,node]*params.del[2,1,node-1] + params.sc[2,2,node]*params.del[2,1,node+1] +
      params.sa[2,4,node]*params.del[4,1,node-1] + params.sc[2,4,node]*params.del[4,1,node+1])

    params.fvec[3,1,node] = params.f[3,1,node] - (
      params.sa[3,3,node]*params.del[3,1,node-1] + params.sc[3,3,node]*params.del[3,1,node+1])

    params.fvec[4,1,node] = params.f[4,1,node] - (
      params.sa[4,2,node]*params.del[2,1,node-1] + params.sc[4,2,node]*params.del[2,1,node+1] +
      params.sa[4,4,node]*params.del[4,1,node-1] + params.sc[4,4,node]*params.del[4,1,node+1])

    params.fvec[5,1,node] = params.f[5,1,node] - (
      params.sa[5,1,node]*params.del[1,1,node-1] + params.sc[5,1,node]*params.del[1,1,node+1] +
      params.sa[5,5,node]*params.del[5,1,node-1] + params.sc[5,5,node]*params.del[5,1,node+1])

    params.fvec[6,1,node] = params.f[6,1,node] - (
      params.sa[6,6,node]*params.del[6,1,node-1] + params.sc[6,6,node]*params.del[6,1,node+1])
  else
    if node == 1
      params.fvec[1,1,node] = params.f[1,1,node] - (
        params.sc[1,1,node]*params.del[1,1,node+1] + params.sc[1,5,node]*params.del[5,1,node+1])

      params.fvec[2,1,node] = params.f[2,1,node] - (
        params.sc[2,2,node]*params.del[2,1,node+1] + params.sc[2,4,node]*params.del[4,1,node+1])

      params.fvec[3,1,node] = params.f[3,1,node] - params.sc[3,3,node]*params.del[3,1,node+1]

      params.fvec[4,1,node] = params.f[4,1,node] - (
        params.sc[4,2,node]*params.del[2,1,node+1] + params.sc[4,4,node]*params.del[4,1,node+1])

      params.fvec[5,1,node] = params.f[5,1,node] - (
        params.sc[5,1,node]*params.del[1,1,node+1] + params.sc[5,5,node]*params.del[5,1,node+1])

      params.fvec[6,1,node] = params.f[6,1,node] - params.sc[6,6,node]*params.del[6,1,node+1]
    else
      params.fvec[1,1,node] = params.f[1,1,node] - (
        params.sa[1,1,node]*params.del[1,1,node-1] + params.sa[1,5,node]*params.del[5,1,node-1])

      params.fvec[2,1,node] = params.f[2,1,node] - (
        params.sa[2,2,node]*params.del[2,1,node-1] + params.sa[2,4,node]*params.del[4,1,node-1])

      params.fvec[3,1,node] = params.f[3,1,node] - params.sa[3,3,node]*params.del[3,1,node-1]

      params.fvec[4,1,node] = params.f[4,1,node] - (
        params.sa[4,2,node]*params.del[2,1,node-1] + params.sa[4,4,node]*params.del[4,1,node-1])

      params.fvec[5,1,node] = params.f[5,1,node] - (
        params.sa[5,1,node]*params.del[1,1,node-1] + params.sa[5,5,node]*params.del[5,1,node-1])

      params.fvec[6,1,node] = params.f[6,1,node] - params.sa[6,6,node]*params.del[6,1,node-1]
    end
  end

  params.del[1, 1, node] = params.si[1, 1, node] * params.fvec[1, 1, node] +
    params.si[1, 5, node] * params.fvec[5, 1, node]

  params.del[2, 1, node] = params.si[2, 2, node] * params.fvec[2, 1, node] +
    params.si[2, 4, node] * params.fvec[4, 1, node]

  params.del[3, 1, node] = params.si[3, 3, node] * params.fvec[3, 1, node]

  params.del[4, 1, node] = params.si[4, 2, node] * params.fvec[2, 1, node] +
    params.si[4, 4, node] * params.fvec[4, 1, node]

  params.del[5, 1, node] = params.si[5, 1, node] * params.fvec[1, 1, node] +
    params.si[5, 5, node] * params.fvec[5, 1, node]

  params.del[6, 1, node] = params.si[6, 6, node] * params.fvec[6, 1, node]
  params
end


text:
# Autogenerated wrapper script for LicenseSpring_jll for i686-w64-mingw32
export libLicenseSpring

using LibCURL_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("LicenseSpring")
JLLWrappers.@declare_library_product(libLicenseSpring, "LicenseSpring.dll")
function __init__()
    JLLWrappers.@generate_init_header(LibCURL_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libLicenseSpring,
        "bin\\LicenseSpring.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module AtomicData

export AtomicStage
export get_atomic_stage, read_NIST
export partition_function, partition_function_interpolator

using DelimitedFiles
using Interpolations
using PeriodicTable
using Unitful
import PhysicalConstants.CODATA2018: h, k_B, c_0

const element_symbols = [el.symbol for el in elements]

include("types.jl")
include("read_utils.jl")
include("partition_function.jl")

end


text:
export approxConvCircular

"""
    $SIGNATURES

Construct (new) subgraph and draw the subgraph associated with clique `frontalSym::Symbol`.

Notes
- See `drawGraphCliq`/`writeGraphPdf` for details on keyword options.

Related

drawGraphCliq, spyCliqMat, drawTree, buildCliqSubgraphUp, buildSubgraphFromLabels!
"""
function drawCliqSubgraphUpMocking(fgl::G,
                                   treel::AbstractBayesTree,
                                   frontalSym::Symbol;
                                   show::Bool=true,
                                   filepath::String="/tmp/caesar/random/cliq_sfg.dot",
                                   engine::AS1="sfdp",
                                   viewerapp::AS2="xdot"  ) where {G <: AbstractDFG, AS1 <: AbstractString, AS2 <: AbstractString}
  #
  sfg = buildCliqSubgraphUp(fgl, treel, frontalSym)
  drawGraph(sfg, show=show, viewerapp=viewerapp, engine=engine, filepath=filepath)
  nothing
end

"""
    $SIGNATURES

Draw and show the factor graph `<:AbstractDFG` via system graphviz and xdot app.

Notes
- Requires system install on Linux of `sudo apt-get install xdot`
- Should not be calling outside programs.
- Need long term solution
- DFG's `toDotFile` a better solution -- view with `xdot` application.
- also try `engine={"sfdp","fdp","dot","twopi","circo","neato"}`

Notes:
- Calls external system application `xdot` to read the `.dot` file format
  - ```toDot(fg,file=...); @async run(`xdot file.dot`)```

Related

drawGraphCliq, [`drawTree`](@ref), printCliqSummary, spyCliqMat
"""
function drawGraph( fgl::AbstractDFG;
                    viewerapp::AbstractString="xdot",
                    filepath::AbstractString="/tmp/caesar/random/fg.dot",
                    engine::AbstractString="neato", #sfdp
                    show::Bool=true )
  #
  mkpath(joinpath( "/", (split(filepath, '/')[1:(end-1)])...) )

  @debug "Writing factor graph file"
  fext = split(filepath, '.')[end]
  fpwoext = filepath[1:(end-length(fext)-1)] # split(filepath, '.')[end-1]
  dotfile = fpwoext*".dot"

  # create the dot file
  DFG.toDotFile(fgl, dotfile)

  try
    # run(`$(engine) $(dotfile) -T$(fext) -o $(filepath)`)
    show ? (@async run(`$(viewerapp) $(dotfile)`)) : nothing
  catch e
    @warn "not able to show $(filepath) with viewerapp=$(viewerapp). Exception e=$(e)"
  end
  nothing
end

"""
    $SIGNATURES

Draw the factor graph from a clique state machine history at a particular step as pdf and show.

Related

drawCliqSubgraphUpMocking, drawGraph, drawTree
"""
function drawGraphCliq( hists::Dict{Int, <: Tuple},
                        step::Int,
                        tree::AbstractBayesTree,
                        frontal::Symbol;
                        show::Bool=true  )
  #
  cid = getId(getClique(tree, frontal))
  cfg = hists[cid][step][4].cliqSubFg
  drawGraph(cfg, show=show)
end

"""
    $SIGNATURES

Print basic statistics about a clique variables and factors.

Related

printCliqHistorySummary
"""
function printCliqSummary(dfg::G,
                          cliq::TreeClique,
                          logger=ConsoleLogger() ) where G <: AbstractDFG
  #
  frtl = getCliqFrontalVarIds(cliq)
  seps = getCliqSeparatorVarIds(cliq)
  fcts = getCliqFactorIdsAll(cliq)

  isinit = map(x->isInitialized(dfg,x), [frtl;seps])
  infdim = map(x->getVariableInferredDim(dfg, x), [frtl;seps])

  with_logger(logger) do
    @info "Clique $(getId(cliq)) summary:"
    @info "  num frontals:    $(length(frtl))"
    @info "  num separators:  $(length(seps))"
    @info "  num factors:     $(length(fcts))"
    @info "  num initialized: $(sum(isinit)) of $(length(isinit))"
    @info ""
    @info "  frontals:  $(frtl)"
    @info "  separator: $(seps)"
    @info "  factors:   $(fcts)"
    @info "  init'ed:   $(Int.(isinit))"
    @info "  infr'dims: $(infdim)"
  end
  nothing
end

"""
    $SIGNATURES

Print basic summary of graph to `logger=ConsoleLogger()`.
"""
function printGraphSummary(dfg::G, logger=ConsoleLogger())::Nothing where {G <: AbstractDFG}
    vars = ls(dfg)
    fcts = lsf(dfg)

    prio = lsfPriors(dfg)

    isinit = map(x->isInitialized(dfg,x), vars)
    infdim = map(x->getVariableInferredDim(dfg, x), vars)
    numedges = map(v->length(ls(dfg, v)), vars)
    numfed = map(fc->length(ls(dfg, fc)), fcts)
    vardims = map(v->getDimension(getVariable(dfg, v)), vars)
    fctdims = map(v->getDimension(getFactor(dfg, v)), fcts)
    priodims = map(v->getDimension(getFactor(dfg, v)), prio)

    with_logger(logger) do
      @info "Distributed Factor Graph summary:"
      @info "  num variables:    $(length(vars))"
      @info "  num factors:      $(length(fcts)), w/ $(length(prio)) priors"
      @info "  var initialized:  $(sum(isinit))"
      @info ""
      @info "  var num edges: min. $(minimum(numedges)) | mean $(round(Statistics.mean(numedges),digits=2)) | 90% $(round(quantile(numedges,0.9),digits=2)) | max. $(maximum(numedges))"
      @info "  fct num edges: min. $(minimum(numfed)) | mean $(round(Statistics.mean(numfed),digits=2)) | 90% $(round(quantile(numfed,0.9),digits=2)) | max. $(maximum(numfed))"
      @info "  Variable dims: min. $(minimum(vardims)) | mean $(round(Statistics.mean(vardims),digits=2)) | 90% $(round(quantile(vardims,0.9),digits=2)) | max. $(maximum(vardims))"
      @info "  Factor dims:   min. $(minimum(fctdims)) | mean $(round(Statistics.mean(fctdims),digits=2)) | 90% $(round(quantile(fctdims,0.9),digits=2)) | max. $(maximum(fctdims))"
      @info "  Prior dimens:  min. $(minimum(priodims)) | mean $(round(Statistics.mean(priodims),digits=2)) | 90% $(round(quantile(priodims,0.9),digits=2)) | max. $(maximum(priodims))"
      @info "  var infr'dims: min. $(minimum(infdim)) | mean $(round(Statistics.mean(infdim),digits=2)) | 90% $(round(quantile(infdim,0.9),digits=2)) | max. $(maximum(infdim))"
    end
    nothing
end

"""
    $SIGNATURES

Print basic summary of graph to `logger=ConsoleLogger()`.
"""
function printSummary(dfg::G, logger=ConsoleLogger()) where G <: AbstractDFG
    printGraphSummary(dfg, logger)
end

"""
    $SIGNATURES

Build an approximate density `[Y|X,DX,.]=[X|Y,DX][DX|.]` as proposed by the conditional convolution.

Notes
- Assume both are on circular manifold, `manikde!(pts, (:Circular,))`
"""
function approxConvCircular(pX::BallTreeDensity, pDX::BallTreeDensity; N::Int=100)
  #

  # building basic factor graph
  tfg = initfg()
  addVariable!(tfg, :s1, Sphere1)
  addVariable!(tfg, :s2, Sphere1)
  addFactor!(tfg, [:s1;:s2], Sphere1Sphere1(pDX), graphinit=false)
  initManual!(tfg,:s1, pX)

  # solve for outgoing proposal value
  approxConv(tfg,:s1s2f1,:s2)
end

function approxConvCircular(pX::BallTreeDensity, pDX::SamplableBelief; N::Int=100)
  pts = reshape(rand(pDX, N), 1, :)
  pC = manikde!(pts, Sphere1)
  approxConvCircular(pX, pC)
end

function approxConvCircular(pX::SamplableBelief, pDX::BallTreeDensity; N::Int=100)
  pts = reshape(rand(pX, N), 1, :)
  pC = manikde!(pts, Sphere1)
  approxConvCircular(pC, pDX)
end

#


text:
# This file was generated, do not modify it. # hide
y2_oos = vcat(y21, y22, y23);
fit!(y2_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y2_oos(), ls="none", marker="o", label="knn oos")
legend()

savefig(joinpath(@OUTPUT, "s3.svg")) # hide


text:
x = range(0, stop=0.2, length=10)
y = zeros(size(x))
z = zeros(size(x))
r = [x'; y'; z']
a = randn(10).^2  # squared to make sure the're all positive

# Different types to check promotion
h = 1.02 * ones(Float32, length(a))
g = 1.04 * ones(Float16, length(a))

weird_zoop = Scatterer(r, a, h, g)
werid_zoop = Scatterer(r, a, 1.02, 1.04)
werid_zoop = Scatterer(r, a, 1.02, 1)

# Make sure all models load
krill1 = Models.krill_mcgeehee
krill2 = Models.krill_conti
cope = Models.calanoid_copepod
sandeel = Models.sandeel
daphnia = Models.daphnia

@assert length(rotate(krill1, roll=90, tilt=30, yaw=12)) ≈ length(krill1)
@assert length(resize(krill1, 0.03)) == 0.03


text:
include("benchtools.jl")
module TensorKitTimers
    using TensorKit
    import ..Timers

    function mpo_timer(f = randn, T = Float64; Vmpo, Vmps, Vphys)
        A = Tensor(f, T, Vmps ⊗ Vphys ⊗ Vmps')
        M = Tensor(f, T, Vmpo ⊗ Vphys ⊗ Vphys' ⊗ Vmpo')
        FL = Tensor(f, T, Vmps ⊗ Vmpo' ⊗ Vmps')
        FR = Tensor(f, T, Vmps ⊗ Vmpo ⊗ Vmps')

        return Timers.Timer(A, M, FL, FR) do A, M, FL, FR
            @tensor C = FL[4,2,1]*A[1,3,6]*M[2,5,3,7]*conj(A[4,5,8])*FR[6,7,8]
            return C
        end
    end

    function pepo_timer(f = randn, T = Float64; Vpepo, Vpeps, Venv, Vphys)
        A = Tensor(f, T, Vpeps ⊗ Vpeps ⊗ Vphys ⊗ Vpeps' ⊗ Vpeps')
        P = Tensor(f, T, Vpepo ⊗ Vpepo ⊗ Vphys ⊗ Vphys' ⊗ Vpepo' ⊗ Vpepo')
        FL = Tensor(f, T, Venv ⊗ Vpeps ⊗ Vpepo' ⊗ Vpeps' ⊗ Venv')
        FD = Tensor(f, T, Venv ⊗ Vpeps ⊗ Vpepo' ⊗ Vpeps' ⊗ Venv')
        FR = Tensor(f, T, Venv ⊗ Vpeps ⊗ Vpepo ⊗ Vpeps' ⊗ Venv')
        FU = Tensor(f, T, Venv ⊗ Vpeps ⊗ Vpepo ⊗ Vpeps' ⊗ Venv')
        return Timers.Timer(A, P, FL, FD, FR, FU) do A, P, FL, FD, FR, FU
            @tensor C = FL[18,7,4,2,1]*FU[1,3,6,9,10]*
                        A[2,17,5,3,11]*P[4,16,8,5,6,12]*conj(A[7,15,8,9,13])*
                        FR[10,11,12,13,14]*FD[14,15,16,17,18]
            return C
        end
    end

    function mera_timer(f = randn, T = Float64; Vmera)
        u = Tensor(f, T, Vmera ⊗ Vmera ⊗ Vmera' ⊗ Vmera')
        w = Tensor(f, T, Vmera ⊗ Vmera ⊗ Vmera')
        ρ = Tensor(f, T, Vmera ⊗ Vmera ⊗ Vmera ⊗ Vmera' ⊗ Vmera' ⊗ Vmera')
        h = Tensor(f, T, Vmera ⊗ Vmera ⊗ Vmera ⊗ Vmera' ⊗ Vmera' ⊗ Vmera')
        return Timers.Timer(u, w, ρ, h) do u, w, ρ, h
            @tensor C = (((((( (h[9,3,4,5,1,2]*u[1,2,7,12]) * conj(u[3,4,11,13]) ) *
                                (u[8,5,15,6]*w[6,7,19]) ) *
                                (conj(u[8,9,17,10])*conj(w[10,11,22])) ) *
                               ( (w[12,14,20] * conj(w[13,14,23])) * ρ[18,19,20,21,22,23])) * w[16,15,18]) * conj(w[16,17,21]))
            return C
        end
    end
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test, Arrow, Tables, Dates, PooledArrays, TimeZones, UUIDs, CategoricalArrays, DataAPI

include(joinpath(dirname(pathof(Arrow)), "ArrowTypes/test/tests.jl"))
include(joinpath(dirname(pathof(Arrow)), "../test/testtables.jl"))
include(joinpath(dirname(pathof(Arrow)), "../test/testappend.jl"))
include(joinpath(dirname(pathof(Arrow)), "../test/integrationtest.jl"))
include(joinpath(dirname(pathof(Arrow)), "../test/dates.jl"))

struct CustomStruct
    x::Int
    y::Float64
    z::String
end

struct CustomStruct2{sym}
    x::Int
end

@testset "Arrow" begin

@testset "table roundtrips" begin

for case in testtables
    testtable(case...)
end

end # @testset "table roundtrips"

@testset "table append" begin

    for case in testtables
        testappend(case...)
    end

    testappend_partitions()

    for compression_option in (:lz4, :zstd)
        testappend_compression(compression_option)
    end

end # @testset "table append"

@testset "arrow json integration tests" begin

for file in readdir(joinpath(dirname(pathof(Arrow)), "../test/arrowjson"))
    jsonfile = joinpath(joinpath(dirname(pathof(Arrow)), "../test/arrowjson"), file)
    println("integration test for $jsonfile")
    df = ArrowJSON.parsefile(jsonfile);
    io = Arrow.tobuffer(df)
    tbl = Arrow.Table(io; convert=false);
    @test isequal(df, tbl)
end

end # @testset "arrow json integration tests"

@testset "misc" begin

# multiple record batches
t = Tables.partitioner(((col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],), (col1=Union{Int64, Missing}[missing,11],)))
io = Arrow.tobuffer(t)
tt = Arrow.Table(io)
@test length(tt) == 1
@test isequal(tt.col1, vcat([1,2,3,4,5,6,7,8,9,missing], [missing,11]))
@test eltype(tt.col1) === Union{Int64, Missing}

# Arrow.Stream
seekstart(io)
str = Arrow.Stream(io)
state = iterate(str)
@test state !== nothing
tt, st = state
@test length(tt) == 1
@test isequal(tt.col1, [1,2,3,4,5,6,7,8,9,missing])

state = iterate(str, st)
@test state !== nothing
tt, st = state
@test length(tt) == 1
@test isequal(tt.col1, [missing,11])

@test iterate(str, st) === nothing

@test isequal(collect(str)[1].col1, [1,2,3,4,5,6,7,8,9,missing])
@test isequal(collect(str)[2].col1, [missing,11])

# dictionary batch isDelta
t = (
    col1=Int64[1,2,3,4],
    col2=Union{String, Missing}["hey", "there", "sailor", missing],
    col3=NamedTuple{(:a, :b), Tuple{Int64, Union{Missing, NamedTuple{(:c,), Tuple{String}}}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b=(c="sailor",)), (a=Int64(4), b=(c="jo-bob",))]
)
t2 = (
    col1=Int64[1,2,5,6],
    col2=Union{String, Missing}["hey", "there", "sailor2", missing],
    col3=NamedTuple{(:a, :b), Tuple{Int64, Union{Missing, NamedTuple{(:c,), Tuple{String}}}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(5), b=(c="sailor2",)), (a=Int64(4), b=(c="jo-bob",))]
)
tt = Tables.partitioner((t, t2))
tt = Arrow.Table(Arrow.tobuffer(tt; dictencode=true, dictencodenested=true))
@test tt.col1 == [1,2,3,4,1,2,5,6]
@test isequal(tt.col2, ["hey", "there", "sailor", missing, "hey", "there", "sailor2", missing])
@test isequal(tt.col3, vcat(NamedTuple{(:a, :b), Tuple{Int64, Union{Missing, NamedTuple{(:c,), Tuple{String}}}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b=(c="sailor",)), (a=Int64(4), b=(c="jo-bob",))], NamedTuple{(:a, :b), Tuple{Int64, Union{Missing, NamedTuple{(:c,), Tuple{String}}}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(5), b=(c="sailor2",)), (a=Int64(4), b=(c="jo-bob",))]))

t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
meta = Dict("key1" => "value1", "key2" => "value2")
Arrow.setmetadata!(t, meta)
meta2 = Dict("colkey1" => "colvalue1", "colkey2" => "colvalue2")
Arrow.setmetadata!(t.col1, meta2)
tt = Arrow.Table(Arrow.tobuffer(t))
@test length(tt) == length(t)
@test tt.col1 == t.col1
@test eltype(tt.col1) === Int64
@test Arrow.getmetadata(tt) == meta
@test Arrow.getmetadata(tt.col1) == meta2

# custom compressors
lz4 = Arrow.CodecLz4.LZ4FrameCompressor(; compressionlevel=8)
Arrow.CodecLz4.TranscodingStreams.initialize(lz4)
t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
tt = Arrow.Table(Arrow.tobuffer(t; compress=lz4))
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

zstd = Arrow.CodecZstd.ZstdCompressor(; level=8)
Arrow.CodecZstd.TranscodingStreams.initialize(zstd)
t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
tt = Arrow.Table(Arrow.tobuffer(t; compress=zstd))
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# custom alignment
t = (col1=Int64[1,2,3,4,5,6,7,8,9,10],)
tt = Arrow.Table(Arrow.tobuffer(t; alignment=64))
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# 53
s = "a" ^ 100
t = (a=[SubString(s, 1:10), SubString(s, 11:20)],)
tt = Arrow.Table(Arrow.tobuffer(t))
@test tt.a == ["aaaaaaaaaa", "aaaaaaaaaa"]

# 49
@test_throws ArgumentError Arrow.Table("file_that_doesnt_exist")

# 52
t = (a=Arrow.DictEncode(string.(1:129)),)
tt = Arrow.Table(Arrow.tobuffer(t))

# 60: unequal column lengths
io = IOBuffer()
@test_throws ArgumentError Arrow.write(io, (a = Int[], b = ["asd"], c=collect(1:100)))

# nullability of custom extension types
t = (a=['a', missing],)
tt = Arrow.Table(Arrow.tobuffer(t))
@test isequal(tt.a, ['a', missing])

# automatic custom struct serialization/deserialization
t = (col1=[CustomStruct(1, 2.3, "hey"), CustomStruct(4, 5.6, "there")],)

Arrow.ArrowTypes.arrowname(::Type{CustomStruct}) = Symbol("JuliaLang.CustomStruct")
Arrow.ArrowTypes.JuliaType(::Val{Symbol("JuliaLang.CustomStruct")}, S) = CustomStruct
tt = Arrow.Table(Arrow.tobuffer(t))
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# 76
t = (col1=NamedTuple{(:a,),Tuple{Union{Int,String}}}[(a=1,), (a="x",)],)
tt = Arrow.Table(Arrow.tobuffer(t))
@test length(tt) == length(t)
@test all(isequal.(values(t), values(tt)))

# 89 etc. - test deprecation paths for old UUID autoconversion + UUID FixedSizeListKind overloads
u = 0x6036fcbd20664bd8a65cdfa25434513f
@test Arrow.ArrowTypes.arrowconvert(UUID, (value=u,)) === UUID(u)
@test Arrow.ArrowTypes.arrowconvert(UUID, u) === UUID(u)
@test Arrow.ArrowTypes.gettype(Arrow.ArrowTypes.ArrowKind(UUID)) == UInt8
@test Arrow.ArrowTypes.getsize(Arrow.ArrowTypes.ArrowKind(UUID)) == 16

# 98
t = (a = [Nanosecond(0), Nanosecond(1)], b = [uuid4(), uuid4()], c = [missing, Nanosecond(1)])
tt = Arrow.Table(Arrow.tobuffer(t))
@test copy(tt.a) isa Vector{Nanosecond}
@test copy(tt.b) isa Vector{UUID}
@test copy(tt.c) isa Vector{Union{Missing,Nanosecond}}

# copy on DictEncoding w/ missing values
x = PooledArray(["hey", missing])
x2 = Arrow.toarrowvector(x)
@test isequal(copy(x2), x)

# some dict encoding coverage

# signed indices for DictEncodedKind #112 #113 #114
av = Arrow.toarrowvector(PooledArray(repeat(["a", "b"], inner = 5)))
@test isa(first(av.indices), Signed)

av = Arrow.toarrowvector(CategoricalArray(repeat(["a", "b"], inner = 5)))
@test isa(first(av.indices), Signed)

av = Arrow.toarrowvector(CategoricalArray(["a", "bb", missing]))
@test isa(first(av.indices), Signed)
@test length(av) == 3
@test eltype(av) == Union{String, Missing}

av = Arrow.toarrowvector(CategoricalArray(["a", "bb", "ccc"]))
@test isa(first(av.indices), Signed)
@test length(av) == 3
@test eltype(av) == String

# 120
x = PooledArray(["hey", missing])
x2 = Arrow.toarrowvector(x)
@test eltype(DataAPI.refpool(x2)) == Union{Missing, String}
@test eltype(DataAPI.levels(x2)) == String
@test DataAPI.refarray(x2) == [1, 2]

# 121
a = PooledArray(repeat(string.('S', 1:130), inner=5), compress=true)
@test eltype(a.refs) == UInt8
av = Arrow.toarrowvector(a)
@test eltype(av.indices) == Int16

# 123
t = (x = collect(zip(rand(10), rand(10))),)
tt = Arrow.Table(Arrow.tobuffer(t))
@test tt.x == t.x

# 144
t = Tables.partitioner(((a=Arrow.DictEncode([1,2,3]),), (a=Arrow.DictEncode(fill(1, 129)),)))
tt = Arrow.Table(Arrow.tobuffer(t))
@test length(tt.a) == 132

# 126
t = Tables.partitioner(
    (
        (a=Arrow.toarrowvector(PooledArray([1,2,3  ])),),
        (a=Arrow.toarrowvector(PooledArray([1,2,3,4])),),
        (a=Arrow.toarrowvector(PooledArray([1,2,3,4,5])),),
    )
)
tt = Arrow.Table(Arrow.tobuffer(t))
@test length(tt.a) == 12
@test tt.a == [1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 4, 5]

t = Tables.partitioner(
    (
        (a=Arrow.toarrowvector(PooledArray([1,2,3  ], signed=true, compress=true)),),
        (a=Arrow.toarrowvector(PooledArray(collect(1:129))),),
    )
)
io = IOBuffer()
@test_throws ErrorException Arrow.write(io, t)

# 75
tbl = Arrow.Table(Arrow.tobuffer((sets = [Set([1,2,3]), Set([1,2,3])],)))
@test eltype(tbl.sets) <: Set

# 85
tbl = Arrow.Table(Arrow.tobuffer((tups = [(1, 3.14, "hey"), (1, 3.14, "hey")],)))
@test eltype(tbl.tups) <: Tuple

# Nothing
tbl = Arrow.Table(Arrow.tobuffer((nothings=[nothing, nothing, nothing],)))
@test tbl.nothings == [nothing, nothing, nothing]

# arrowmetadata
t = (col1=[CustomStruct2{:hey}(1), CustomStruct2{:hey}(2)],)
ArrowTypes.arrowname(::Type{<:CustomStruct2}) = Symbol("CustomStruct2")
tbl = Arrow.Table(Arrow.tobuffer(t))
# test we get the warning about deserializing
@test eltype(tbl.col1) <: NamedTuple
ArrowTypes.arrowmetadata(::Type{CustomStruct2{sym}}) where {sym} = sym
ArrowTypes.JuliaType(::Val{:CustomStruct2}, S, meta) = CustomStruct2{Symbol(meta)}
tbl = Arrow.Table(Arrow.tobuffer(t))
@test eltype(tbl.col1) == CustomStruct2{:hey}

# 170
tbl = Arrow.Table(Arrow.tobuffer((x = [1,2,3],)))
m = Dict("a" => "b")
Arrow.setmetadata!(tbl, m)
@test Arrow.getmetadata(tbl) === m

# 166
t = (
    col1=[zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.NANOSECOND, nothing})],
)
tbl = Arrow.Table(Arrow.tobuffer(t))
@test tbl.col1[1] == Dates.DateTime(1970)

# 95; Arrow.ToTimestamp
x = [ZonedDateTime(Dates.DateTime(2020), tz"Europe/Paris")]
c = Arrow.ToTimestamp(x)
@test eltype(c) == Arrow.Timestamp{Arrow.Flatbuf.TimeUnitModule.MILLISECOND, Symbol("Europe/Paris")}
@test c[1] == Arrow.Timestamp{Arrow.Flatbuf.TimeUnitModule.MILLISECOND, Symbol("Europe/Paris")}(1577836800000)

# 158
# arrow ipc stream generated from pyarrow with no record batches
bytes = UInt8[0xff, 0xff, 0xff, 0xff, 0x78, 0x00, 0x00, 0x00, 0x10, 0x00, 0x00, 0x00, 0x00, 0x00, 0x0a, 0x00, 0x0c, 0x00,
    0x06, 0x00, 0x05, 0x00, 0x08, 0x00, 0x0a, 0x00, 0x00, 0x00, 0x00, 0x01, 0x04, 0x00, 0x0c, 0x00, 0x00, 0x00, 0x08, 0x00,
    0x08, 0x00, 0x00, 0x00, 0x04, 0x00, 0x08, 0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x14, 0x00,
    0x00, 0x00, 0x10, 0x00, 0x14, 0x00, 0x08, 0x00, 0x06, 0x00, 0x07, 0x00, 0x0c, 0x00, 0x00, 0x00, 0x10, 0x00, 0x10, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x01, 0x02, 0x10, 0x00, 0x00, 0x00, 0x1c, 0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x01, 0x00, 0x00, 0x00, 0x61, 0x00, 0x00, 0x00, 0x08, 0x00, 0x0c, 0x00, 0x08, 0x00, 0x07, 0x00, 0x08, 0x00,
    0x00, 0x00, 0x00, 0x00, 0x00, 0x01, 0x40, 0x00, 0x00, 0x00, 0xff, 0xff, 0xff, 0xff, 0x00, 0x00, 0x00, 0x00]
tbl = Arrow.Table(bytes)
@test length(tbl.a) == 0
@test eltype(tbl.a) == Union{Int64, Missing}

# 181
tbl = (x = [Dict()],)
d = tbl.x[];
for i in 1:20
    d[i] = Dict()
    d = d[i]
end
msg = "reached nested serialization level (42) deeper than provided max depth argument (41); to increase allowed nesting level, pass `maxdepth=X`"
@test_throws ErrorException(msg) Arrow.tobuffer(tbl; maxdepth=41).x
@test Arrow.Table(Arrow.tobuffer(tbl; maxdepth=42)).x == tbl.x

# 167
t = (
    col1=[["boop", "she"], ["boop", "she"], ["boo"]],
)
tbl = Arrow.Table(Arrow.tobuffer(t))
@test eltype(tbl.col1) == Vector{String}

end # @testset "misc"

end


text: base_global = "hello world"


text:
using LinearAlgebraStudy
using Test

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

max_stars_repo_name: gjunqueira-sys/LinearAlgebraStudy


text:
# we don't use Some because it does not behave like a container, but more like `Ref`
# for details/update see https://github.com/JuliaLang/julia/issues/35911

"""
    Option{T} = Union{Const{Nothing}, Identity{T}}
    Option(nothing) == Const(nothing)
    Option("anything but nothing") == Identity("anything but nothing")
    Option() == Const(nothing)

Like `Union{T, Nothing}`, however with container semantics. While `Union{T, Nothing}` can be thought of like
a value which either exists or not, `Option{T} = Union{Identity{T}, Const{Nothing}}` is a container which
is either empty or has exactly one element.

We reuse [`Identity`](@ref) as representing the single-element-container and `Const(nothing)` as the empty container.
"""
const Option{T} = Union{Const{Nothing}, Identity{T}}

Option(::Nothing) = Const(nothing)
Option(a) = Identity(a)
Option() = Const(nothing)

"""
    iftrue(bool_condition, value)
    iftrue(func, bool_condition)
    iftrue(bool_condition) do
      # only executed if bool_condition is true
    end

Helper to create an Option based on some condition. If `bool_condition` is
true, the function `func` is called with no arguments, and its result wrapped into `Identity`.
If `bool_condition` is false, `Option()` is returned.
"""
function iftrue(func::Function, b::Bool)
  b ? Identity(func()) : Const(nothing)
end
function iftrue(b::Bool, t)
  b ? Identity(t) : Const(nothing)
end

"""
    iffalse(bool_condition, value)
    iffalse(func, bool_condition)
    iffalse(bool_condition) do
      # only executed if bool_condition is true
    end

Helper to create an Option based on some condition. If `bool_condition` is
false, the function `func` is called with no arguments, and its result wrapped into `Identity`.
If `bool_condition` is true, `Option()` is returned.

Precisely the opposite of [`iftrue`](@ref).
"""
function iffalse(func::Function, b::Bool)
  !b ? Identity(func()) : Const(nothing)
end
function iffalse(b::Bool, t)
  !b ? Identity(t) : Const(nothing)
end

"""
    isoption(::Const{Nothing}) = true
    isoption(::Identity) = true
    isoption(other) = false

check whether something is an [`Option`](@ref)
"""
isoption(::Const{Nothing}) = true
isoption(::Identity) = true
isoption(other) = false

"""
    issome(::Identity) = true
    issome(::Const{Nothing}) = false

Similar to [`isright`](@ref), but only defined for `Const{Nothing}`. Will
throw MethodError when applied on other `Const`.
"""
issome(::Identity) = true
issome(::Const{Nothing}) = false

"""
    isnone(::Identity) = false
    isnone(::Const{Nothing}) = true

Similar to [`isleft`](@ref), but only defined for `Const{Nothing}`. Will
throw MethodError when applied on other `Const`.
"""
isnone(::Identity) = false
isnone(::Const{Nothing}) = true

Base.eltype(::Type{Option{T}}) where {T} = T
Base.eltype(::Type{Option}) = Any


text:
function luhn(s::AbstractString)
    length(s) < 2 && return false
    tonumber(c::Char) = Int(c - codepoint('0'))
    s = replace(s, ' ' => "")
    !all(isdigit,s) && return false
    checksum = sum(
        tuple -> begin
            (i,c) = tuple
            rem(i,2) == 0 ? rem(tonumber(c)*2,9) : tonumber(c)
        end,
        s|>reverse|>enumerate
    )
    rem(checksum,10) == 0
end

# Slower but cooler
function luhn2(s::AbstractString)
    length(s) < 2 && return false
    s = replace(s, ' ' => "")
    !all(isdigit,s) && return false
    digits = parse.(Int,s|>collect)
    digits[end-1:-2:1] = rem.(digits[end-1:-2:1].*2,9)
    rem(sum(digits),10) == 0
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ClientIPConfig <: SwaggerModel
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiCoreV1ClientIPConfig(;timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ClientIPConfig, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiCoreV1ClientIPConfig

const _property_map_IoK8sApiCoreV1ClientIPConfig = Dict{Symbol,Symbol}(Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiCoreV1ClientIPConfig = Dict{Symbol,String}(Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiCoreV1ClientIPConfig }) = collect(keys(_property_map_IoK8sApiCoreV1ClientIPConfig))
Swagger.property_type(::Type{ IoK8sApiCoreV1ClientIPConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ClientIPConfig[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ClientIPConfig }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ClientIPConfig[property_name]

function check_required(o::IoK8sApiCoreV1ClientIPConfig)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ClientIPConfig }, name::Symbol, val)
end


text:
module BudanTables

using DocStringExtensions:SIGNATURES
using Requires

using RecipesBase
using Polynomials # Polynomial type, roots function
# using Plots

export BudanTable, budantable, set_getting_roots

const _ROOTS = Ref{Function}(Polynomials.roots)

"""
$(SIGNATURES)

Sets the function used to compute the roots of polynomials building a Budan's table.
The input is a function `f` computing the roots of a polynomial represented as an object of type `T`, for which `degree(p::T)` and `derivative(p::T)` must be defined.

See also [`budantable`](@ref)
"""
function set_getting_roots(f)
    _ROOTS[] = f
end

function __init__()
        @require PolynomialRoots="3a141323-8675-5d76-9d11-e1df1406c778" begin
            # Use Polynomials api for poly arith. and PR's roots for roots
            _BigFloat(A::AbstractArray) = BigFloat.(A)
            set_getting_roots(PolynomialRoots.roots∘_BigFloat∘Polynomials.coeffs)
        end
end

include("BudanTablesBase.jl")
include("Recipes.jl")

end


text:
using AbstractDomains

a = HyperBox([0.0 0.0 0.0
              1.0 1.0 1.0])
b = HyperBox([0.0 0.0 0.0 0.0
              1.0 1.0 1.0 1.0])
@test mid(a) == [0.5,0.5,0.5]
@test ndims(a) == 3
@test a != b
@test isrelational(HyperBox) == false

# Mid
c = HyperBox([0.0 0.0
              1.0 1.0])
prev = prevfloat(0.5)
next = nextfloat(0.5)
c1 = HyperBox([0.0 0.0
              0.5 0.5])
c2 = HyperBox([next 0.0
              1.0 0.5])
c3 = HyperBox([0.0 next
               0.5 1.0])
c4 = HyperBox([next next
               1.0  1.0])
q = mid_split(c)

# Partial Split
d = HyperBox([0.0 0.0 0.0
              10.0 20.0 1.0])
split_points = Dict(1 => 5.0, 3 => 0.5)
@test length(partial_split_box(d,split_points)) == 4
@test length(mid_partial_split(d,[1,2,3])) == 8
split_points2 = Dict(1 => 5.0, 2=> 0.1, 3 => 0.5)
@test length(partial_split_box(d,split_points2)) == 8

@test isequal(q[1],c1)
@test isequal(q[2],c2)
@test isequal(q[3],c3)
@test isequal(q[4],c4)

@test all(Bool[0 <= rb <= 1 for rb in rand(b)])


text:
# ---
# title: 1290. Convert Binary Number in a Linked List to Integer
# id: problem1290
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List, Bit Manipulation
# link: <https://leetcode.com/problems/convert-binary-number-in-a-linked-list-to-integer/description/>
# hidden: true
# ---
#
# Given `head` which is a reference node to a singly-linked list. The value of
# each node in the linked list is either 0 or 1. The linked list holds the
# binary representation of a number.
#
# Return the _decimal value_ of the number in the linked list.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/05/graph-1.png)
#
#
#
#     Input: head = [1,0,1]
#     Output: 5
#     Explanation: (101) in base 2 = (5) in base 10
#
#
# **Example 2:**
#
#
#
#     Input: head = [0]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: head = [1]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: head = [1,0,0,1,0,0,1,1,1,0,0,0,0,0,0]
#     Output: 18880
#
#
# **Example 5:**
#
#
#
#     Input: head = [0,0]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * The Linked List is not empty.
#   * Number of nodes will not exceed `30`.
#   * Each node's value is either `0` or `1`.
#
#
## @lc code=start
using LeetCode

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


text:
module DSPHelp
export ev, od

function ev(x::Array{Float64,2})
  return (x + fliplr(x))/2
end

function od(x::Array{Float64,2})
  return (x - fliplr(x))/2
end

end # module


text:
function invariants!(W1::Vector{Float64},W2::Vector{Float64},lf::Vector{Float64},
    lb::Vector{Float64},A::Vector{Float64},Q::Vector{Float64},
    c0::Float64,beta::Float64,rho::Float64,μ::Float64,diff::Float64,
    Qs::Float64,As::Float64)

    # characteristic speeds
    lb[1] = (Q[1]./A[1].*Qs./As .- sqrt.(0.5.*beta./rho).*(A[1].*As).^0.25).*As./Qs;
    lf[1] = (Q[end]./A[end].*Qs./As .+ sqrt.(0.5.*beta./rho).*(A[end].*As).^0.25).*As./Qs;

    # invariants
    W1[:] .= (Q[end-1:end]./A[end-1:end].*Qs./As .+
        4.*(sqrt.(0.5.*beta./rho).*(A[end-1:end].*As).^0.25 .- c0)).*As./Qs;
    W2[:] .= (Q[1:2]./A[1:2].*Qs./As .-
        4.*(sqrt.(0.5.*beta./rho).*(A[1:2].*As).^0.25 .- c0)).*As./Qs;

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
#=
epigraph of Euclidean norm (AKA second-order cone)
(u, w) : u ≥ ‖w‖
self-dual

extended formulation
∃ λ, u ≥ 2 Σᵢ λᵢ, (u, λᵢ, wᵢ) ∈ EpiPerSquare
i.e. λᵢ ≥ 0, 2 u λᵢ ≥ wᵢ²
=#

mutable struct EpiNormEucl{E <: NatExt} <: Cache
    oa_s::Vector{AE}
    d::Int
    λ::Vector{VR}
    EpiNormEucl{E}() where {E <: NatExt} = new{E}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.EpiNormEuclCone{RealF},
    opt::Optimizer,
)
    dim = MOI.dimension(cone)
    @assert dim == length(oa_s)
    d = dim - 1
    E = nat_or_ext(opt, d)
    cache = EpiNormEucl{E}()
    cache.oa_s = oa_s
    cache.d = d
    return cache
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::EpiNormEucl,
    opt::Optimizer,
)
    return _get_cuts(z[2:end], cache, opt)
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::EpiNormEucl,
    opt::Optimizer,
)
    us = s[1]
    @views ws = s[2:end]
    ws_norm = LinearAlgebra.norm(ws)
    if us - ws_norm > -opt.tol_feas
        return AE[]
    end

    # gradient cut is (1, -ws / ‖ws‖)
    r = ws / -ws_norm
    return _get_cuts(r, cache, opt)
end

# unextended formulation

function MOIPajarito.Cones.add_init_cuts(cache::EpiNormEucl{Nat}, opt::Optimizer)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add cuts u ≥ |wᵢ|
    @views w = cache.oa_s[2:end]
    JuMP.@constraints(opt.oa_model, begin
        u .>= w
        u .>= -w
    end)
    return
end

function _get_cuts(r::Vector{RealF}, cache::EpiNormEucl{Nat}, opt::Optimizer)
    # strengthened cut is (‖r‖, r)
    clean_array!(r) && return AE[]
    p = LinearAlgebra.norm(r)
    u = cache.oa_s[1]
    @views w = cache.oa_s[2:end]
    cut = JuMP.@expression(opt.oa_model, p * u + JuMP.dot(r, w))
    return [cut]
end

# extended formulation

MOIPajarito.Cones.num_ext_variables(cache::EpiNormEucl{Ext}) = cache.d

function MOIPajarito.Cones.extend_start(
    cache::EpiNormEucl{Ext},
    s_start::Vector{RealF},
    opt::Optimizer,
)
    u_start = s_start[1]
    w_start = s_start[2:end]
    if u_start < 1e-9
        return zeros(cache.d)
    end
    return [w_i / 2u_start * w_i for w_i in w_start]
end

function MOIPajarito.Cones.setup_auxiliary(cache::EpiNormEucl{Ext}, opt::Optimizer)
    @assert cache.d >= 2
    λ = cache.λ = JuMP.@variable(opt.oa_model, [1:(cache.d)], lower_bound = 0)
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 2 * sum(λ))
    return λ
end

function MOIPajarito.Cones.add_init_cuts(cache::EpiNormEucl{Ext}, opt::Optimizer)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add disaggregated cuts (1, 2, ±2) on (u, λᵢ, wᵢ), implying u ≥ |wᵢ|
    @views w = cache.oa_s[2:end]
    d = cache.d
    λ = cache.λ
    JuMP.@constraints(opt.oa_model, begin
        [i in 1:d], u + 2λ[i] + 2w[i] >= 0
        [i in 1:d], u + 2λ[i] - 2w[i] >= 0
    end)
    return
end

function _get_cuts(r::Vector{RealF}, cache::EpiNormEucl{Ext}, opt::Optimizer)
    clean_array!(r) && return AE[]
    p = LinearAlgebra.norm(r)
    u = cache.oa_s[1]
    @views w = cache.oa_s[2:end]
    λ = cache.λ
    cuts = AE[]
    for i in 1:(cache.d)
        r_i = r[i]
        iszero(r_i) && continue
        # strengthened disaggregated cut on (u, λᵢ, wᵢ) is (rᵢ² / 2‖r‖, ‖r‖, rᵢ)
        cut = JuMP.@expression(opt.oa_model, r_i^2 / 2p * u + p * λ[i] + r_i * w[i])
        push!(cuts, cut)
    end
    return cuts
end


text:
"""
    setExperimentSummary(genj::GeneticJuliaStructure = GenJ; batchSize::UInt16 = -1)

documentation
"""
function setExperimentSummary(genj::GeneticJuliaStructure = GenJ; batchSize::Integer = 0,
                              displayFitness::Bool = true, displayBestFitness::Bool = true,
                              displayFitnessMean::Bool = true, displayFitnessSTD::Bool = true,
                              printDuringExperiment::Bool = false, outputFile::String = "")

    if batchSize < 0
        GenJ._experimentInfo._experimentSummary = ExperimentSummary(outputFile,
                                                                    -1, false,
                                                                    Array{Float64}(undef, 0),
                                                                    Array{Float64}(undef, 0),
                                                                    Array{Float64}(undef, 0),
                                                                    Array{Float64}(undef, 0),
                                                                    Array{Any}(undef, 0))
    else
        if !(isdefined(genj, :_generator))
            error("Generator must be defined first")
        elseif !(isdefined(genj, :_evaluator))
            error("Evaluator must be defined first")
        elseif !(isdefined(genj, :_stopCondition))
            error("StopCondition must be defined first")
        end

        popSize = getPopSize(genj._generator)
        nFitness = getNumberFitness(genj._evaluator)
        maxIter = getMaxIterations(genj._stopCondition)

        if batchSize == 0
            batchSize = maxIter
        end

        intDiv = div(maxIter, batchSize)

        # Using column-first order for higher performance
        if displayFitness
            if hasGlobalFitnessFunction(genj._evaluator)
                fitnessValues = Array{Float64}(undef, (nFitness+1, popSize, intDiv))
            else
                fitnessValues = Array{Float64}(undef, (nFitness, popSize, intDiv))
            end
        else
            fitnessValues = Array{Float64}(undef, 0)
        end

        if displayBestFitness
            if hasGlobalFitnessFunction(genj._evaluator)
                bestFitnessValues = Array{Float64}(undef, (nFitness+1, intDiv))
            else
                bestFitnessValues = Array{Float64}(undef, (nFitness, intDiv))
            end
            bestIndividuals = Array{Any}(undef, intDiv)
        else
            bestFitnessValues = Array{Float64}(undef, 0)
            bestIndividuals = Array{Any}(undef, 0)
        end

        if displayFitnessMean
            if hasGlobalFitnessFunction(genj._evaluator)
                meanFitness = Array{Float64}(undef, (nFitness+1, intDiv))
            else
                meanFitness = Array{Float64}(undef, (nFitness, intDiv))
            end
        else
            meanFitness = Array{Float64}(undef, 0)
        end

        if displayFitnessSTD
            if hasGlobalFitnessFunction(genj._evaluator)
                stdFitness = Array{Float64}(undef, (nFitness+1, intDiv))
            else
                stdFitness = Array{Float64}(undef, (nFitness, intDiv))
            end
        else
            stdFitness = Array{Float64}(undef, 0)
        end

        GenJ._experimentInfo._experimentSummary = ExperimentSummary(outputFile,
                                                                    batchSize,
                                                                    printDuringExperiment,
                                                                    fitnessValues,
                                                                    bestFitnessValues,
                                                                    meanFitness,
                                                                    stdFitness,
                                                                    bestIndividuals)
    end

    return nothing
end # function

"""
    printInformation_(summary::ExperimentSummary)

Prints all the information in the summary.
"""
function printInformation_(summary::ExperimentSummary; printFitness::Bool = true,
                           printBestFitness::Bool = true, printMeanFitness::Bool = true,
                           printSTDFitness::Bool = true, outputFile::String = "")
    if somethingToDisplay(summary) && (printFitness || printBestFitness ||
        printMeanFitness || printSTDFitness)

        if outputFile == ""
            outputFile = getOutputFile(summary)
        elseif outputFile == "Base.stdout"
            outputFile = ""
        end
        io = outputFile != "" ? open(outputFile, "a") : Base.stdout
        flush(io)

        batchSize = getBatchSize(summary)
        anyDisplayable = displayAnything(summary)
        nGens = axes(anyDisplayable)[end]

        for gen = nGens
            println(io, "GENERATION ", batchSize * gen, ":")
            println(io)

            if displayFitness(summary) && printFitness
                println(io, "  RESULTS OF ALL THE INDIVIDUALS")
                println(io, "  ------------------------------")
                for ind = eachindex(view(summary._fitnessValues, 1, :, 1))
                    println(io, "\tIndividual ", ind, ":")
                    for fit = eachindex(view(summary._fitnessValues, :, 1, 1))
                        println(io, "\t  Fitness Value ", fit, ": ", summary._fitnessValues[fit, ind, gen])
                    end
                    println(io)
                end
                println(io)
            end

            if displayBestFitness(summary) && printBestFitness
                println(io, "  RESULTS OF THE BEST INDIVIDUAL")
                println(io, "  ------------------------------")
                for fit = eachindex(view(summary._bestFitnessValues, :, 1))
                    println(io, "\tFitness ", fit, " of the Best Individual: ", summary._bestFitnessValues[fit, gen])
                end
                println(io, "\tBest Individual: ", summary._bestIndividuals[gen])
                println(io)
            end

            if displayMeanFitness(summary) && printMeanFitness
                println(io, "  MEAN OF THE RESULTS")
                println(io, "  -------------------")
                for fit = eachindex(view(summary._meanFitness, :, 1))
                    println(io, "\tMean of Fitness ", fit, ": ", summary._meanFitness[fit, gen])
                end
                println(io)
            end

            if displayBestFitness(summary) && printSTDFitness
                println(io, "  STANDARD DEVIATION OF THE RESULTS")
                println(io, "  ---------------------------------")
                for fit = eachindex(view(summary._stdFitness, :, 1))
                    println(io, "\tStandard Deviation of Fitness ", fit, ": ", summary._stdFitness[fit, gen])
                end
                println(io)
            end

            println(io, "============================")
            println(io)
            println(io)

            io == Base.stdout || close(io)
        end
    elseif !(somethingToDisplay(summary))
        println("Nothing to display! Set the information that you want to display by calling setExperimentSummary function")
    else
        println("Nothing to display!")
    end

    return nothing
end # function

"""
    printInformation_(summary::ExperimentSummary)

Prints all the information in the summary.
"""
function printInformationWithGlobal_(summary::ExperimentSummary; printFitness::Bool = true,
                                     printBestFitness::Bool = true, printMeanFitness::Bool = true,
                                     printSTDFitness::Bool = true, outputFile::String = "")
    if somethingToDisplay(summary) && (printFitness || printBestFitness ||
        printMeanFitness || printSTDFitness)

        if outputFile == ""
            outputFile = getOutputFile(summary)
        elseif outputFile == "Base.stdout"
            outputFile = ""
        end
        io = outputFile != "" ? open(outputFile, "a") : Base.stdout
        flush(io)

        batchSize = getBatchSize(summary)
        anyDisplayable = displayAnything(summary)
        nGens = axes(anyDisplayable)[end]

        for gen = nGens
            println(io, "GENERATION ", batchSize * gen, ":")
            println(io)

            if displayFitness(summary) && printFitness
                println(io, "  RESULTS OF ALL THE INDIVIDUALS")
                println(io, "  ------------------------------")
                for ind = eachindex(view(summary._fitnessValues, 1, :, 1))
                    println(io, "\tIndividual ", ind, ":")
                    for fit = eachindex(view(summary._fitnessValues, :, 1, 1))
                        if fit != lastindex(view(summary._fitnessValues, :, 1, 1))
                            println(io, "\t  Fitness Value ", fit, ": ", summary._fitnessValues[fit, ind, gen])
                        else
                            println(io, "\t  Global Fitness Value: ", summary._fitnessValues[fit, ind, gen])
                        end
                    end
                    println(io)
                end
                println(io)
            end

            if displayBestFitness(summary) && printBestFitness
                println(io, "  RESULTS OF THE BEST INDIVIDUAL")
                println(io, "  ------------------------------")
                for fit = eachindex(view(summary._bestFitnessValues, :, 1))
                    if fit != lastindex(view(summary._bestFitnessValues, :, 1))
                        println(io, "\tFitness Value ", fit, " of the Best Individual: ", summary._bestFitnessValues[fit, gen])
                    else
                        println(io, "\tGlobal Fitness Value of the Best Individual: ", summary._bestFitnessValues[fit, gen])
                    end
                end
                println(io, "\tBest Individual: ", summary._bestIndividuals[gen])
                println(io)
            end

            if displayMeanFitness(summary) && printMeanFitness
                println(io, "  MEAN OF THE RESULTS")
                println(io, "  -------------------")
                for fit = eachindex(view(summary._meanFitness, :, 1))
                    if fit != lastindex(view(summary._meanFitness, :, 1))
                        println(io, "\tMean of Fitness ", fit, ": ", summary._meanFitness[fit, gen])
                    else
                        println(io, "\tMean of Global Fitness: ", summary._meanFitness[fit, gen])
                    end
                end
                println(io)
            end

            if displayBestFitness(summary) && printSTDFitness
                println(io, "  STANDARD DEVIATION OF THE RESULTS")
                println(io, "  ---------------------------------")
                for fit = eachindex(view(summary._stdFitness, :, 1))
                    if fit != lastindex(view(summary._stdFitness, :, 1))
                        println(io, "\tStandard Deviation of Fitness ", fit, ": ", summary._stdFitness[fit, gen])
                    else
                        println(io, "\tStandard Deviation of Global Fitness: ", summary._stdFitness[fit, gen])
                    end
                end
                println(io)
            end

            println(io, "============================")
            println(io)
            println(io)

            io == Base.stdout || close(io)
        end
    elseif !(somethingToDisplay(summary))
        println("Nothing to display! Set the information that you want to display by calling setExperimentSummary function")
    else
        println("Nothing to display!")
    end

    return nothing
end # function

"""
    printLastInformation_(summary::ExperimentSummary)

Prints the last bit of information collected of the summary.
"""
function printLastInformation_(summary::ExperimentSummary, currGeneration::Integer)
    if somethingToDisplay(summary)

        outputFile = getOutputFile(summary)
        io = outputFile != "" ? open(outputFile, "a") : Base.stdout
        flush(io)
        currGeneration = div(currGeneration, getBatchSize(summary))

        if displayFitness(summary)
            println(io, "  RESULTS OF ALL THE INDIVIDUALS")
            println(io, "  ------------------------------")
            for ind = eachindex(view(summary._fitnessValues, 1, :, 1))
                println(io, "\tIndividual ", ind, ":")
                for fit = eachindex(view(summary._fitnessValues, :, 1, 1))
                    println(io, "\t  Fitness Value ", fit, ": ", summary._fitnessValues[fit, ind, currGeneration])
                end
                println(io)
            end
            println(io)
        end

        if displayBestFitness(summary)
            println(io, "  RESULTS OF THE BEST INDIVIDUAL")
            println(io, "  ------------------------------")
            for fit = eachindex(view(summary._bestFitnessValues, :, 1))
                println(io, "\tFitness Value ", fit, " of the Best Individual: ", summary._bestFitnessValues[fit, currGeneration])
            end
            println(io, "\tBest Individual: ", summary._bestIndividuals[currGeneration])
            println(io)
        end

        if displayMeanFitness(summary)
            println(io, "  MEAN OF THE RESULTS")
            println(io, "  -------------------")
            for fit = eachindex(view(summary._meanFitness, :, 1))
                println(io, "\tMean of Fitness ", fit, ": ", summary._meanFitness[fit, currGeneration])
            end
            println(io)
        end

        if displayBestFitness(summary)
            println(io, "  STANDARD DEVIATION OF THE RESULTS")
            println(io, "  ---------------------------------")
            for fit = eachindex(view(summary._stdFitness, :, 1))
                println(io, "\tStandard Deviation of Fitness ", fit, ": ", summary._stdFitness[fit, currGeneration])
            end
            println(io)
        end

        io == Base.stdout || close(io)
    end

    return nothing
end # function

"""
    printLastInformation_(summary::ExperimentSummary)

Prints the last bit of information collected of the summary.
"""
function printLastInformationWithGlobal_(summary::ExperimentSummary, currGeneration::Integer)
    if somethingToDisplay(summary)

        outputFile = getOutputFile(summary)
        io = outputFile != "" ? open(outputFile, "a") : Base.stdout
        flush(io)
        currGeneration = div(currGeneration, getBatchSize(summary))

        if displayFitness(summary)
            println(io, "  RESULTS OF ALL THE INDIVIDUALS")
            println(io, "  ------------------------------")
            for ind = eachindex(view(summary._fitnessValues, 1, :, 1))
                println(io, "\tIndividual ", ind, ":")
                for fit = eachindex(view(summary._fitnessValues, :, 1, 1))
                    if fit != lastindex(view(summary._fitnessValues, :, 1, 1))
                        println(io, "\t  Fitness Value ", fit, ": ", summary._fitnessValues[fit, ind, currGeneration])
                    else
                        println(io, "\t  Global Fitness Value: ", summary._fitnessValues[fit, ind, currGeneration])
                    end
                end
                println(io)
            end
            println(io)
        end

        if displayBestFitness(summary)
            println(io, "  RESULTS OF THE BEST INDIVIDUAL")
            println(io, "  ------------------------------")
            for fit = eachindex(view(summary._bestFitnessValues, :, 1))
                if fit != lastindex(view(summary._bestFitnessValues, :, 1))
                    println(io, "\tFitness Value", fit, " of the Best Individual: ", summary._bestFitnessValues[fit, currGeneration])
                else
                    println(io, "\tGlobal Fitness Value of the Best Individual: ", summary._bestFitnessValues[fit, currGeneration])
                end
            end
            println(io, "\tBest Individual: ", summary._bestIndividuals[currGeneration])
            println(io)
        end

        if displayMeanFitness(summary)
            println(io, "  MEAN OF THE RESULTS")
            println(io, "  -------------------")
            for fit = eachindex(view(summary._meanFitness, :, 1))
                if fit != lastindex(view(summary._meanFitness, :, 1))
                    println(io, "\tMean of Fitness ", fit, ": ", summary._meanFitness[fit, currGeneration])
                else
                    println(io, "\tMean of Global Fitness: ", summary._meanFitness[fit, currGeneration])
                end
            end
            println(io)
        end

        if displayBestFitness(summary)
            println(io, "  STANDARD DEVIATION OF THE RESULTS")
            println(io, "  ---------------------------------")
            for fit = eachindex(view(summary._stdFitness, :, 1))
                if fit != lastindex(view(summary._stdFitness, :, 1))
                    println(io, "\tStandard Deviation of Fitness ", fit, ": ", summary._stdFitness[fit, currGeneration])
                else
                    println(io, "\tStandard Deviation of Global Fitness: ", summary._stdFitness[fit, currGeneration])
                end
            end
            println(io)
        end

        io == Base.stdout || close(io)
    end

    return nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
@doc (@doc NodalData)
mutable struct NodalChannel <: GphysChannel
  ox::Float64                         # origin x
  oy::Float64                         # origin y
  oz::Float64                         # origin z
  id::String                          # id
  name::String                        # name
  loc::InstrumentPosition             # loc
  fs::Float64                         # fs
  gain::Float64                       # gain
  resp::InstrumentResponse            # resp
  units::String                       # units
  src::String                         # src
  misc::Dict{String,Any}              # misc
  notes::Array{String,1}              # notes
  t::Array{Int64,2}                   # time
  x::FloatArray                       # data

  function NodalChannel()
    return new(
                default_fs, default_fs, default_fs,
                "",
                "",
                NodalLoc(),
                default_fs,
                default_gain,
                PZResp(),
                "",
                "",
                Dict{String,Any}(),
                Array{String,1}(undef,0),
                Array{Int64,2}(undef,0,2),
                Array{Float32,1}(undef,0)
              )
  end

  function NodalChannel(
            ox::Float64,
            oy::Float64,
            oz::Float64,
            id::String                          , # id
            name::String                        , # name
            loc::InstrumentPosition             , # loc
            fs::Float64                         , # fs
            gain::Float64                       , # gain
            resp::InstrumentResponse            , # resp
            units::String                       , # units
            src::String                         , # src
            misc::Dict{String,Any}              , # misc
            notes::Array{String,1}              , # notes
            t::Array{Int64,2}                   , # time
            x::FloatArray
            )

    return new(ox, oy, oz, id, name, loc, fs, gain, resp, units, src, misc, notes, t, x)
  end
end

function sizeof(S::NodalChannel)
  s = 120
  for f in nodalfields
    v = getfield(S, f)
    s += sizeof(v)
    if f == :notes
      if !isempty(v)
        s += sum([sizeof(j) for j in v])
      end
    elseif f == :misc
      k = collect(keys(v))
      s += sizeof(k) + 64 + sum([sizeof(j) for j in k])
      for p in values(v)
        s += sizeof(p)
        if typeof(p) == Array{String,1}
          s += sum([sizeof(j) for j in p])
        end
      end
    end
  end
  return s
end

function getindex(S::NodalData, j::Int)
  C = NodalChannel()
  [setfield!(C, f, getfield(S,f)[j]) for f in nodalfields]
  C.x = copy(S.data[:,j])
  return C
end

function setindex!(S::NodalData, C::NodalChannel, j::Int)
  [(getfield(S, f))[j] = getfield(C, f) for f in nodalfields]
  S.data[:,j] .= C.x
  S.x[j] = view(S.data, :, j)
  return nothing
end

function isempty(C::NodalChannel)
  q::Bool = C.gain == default_gain
  for f in (:ox, :oy, :oz, :fs)
    q = min(q, getfield(C, f) == default_fs)
    (q == false) && return q
  end
  for f in (:id, :loc, :misc, :name, :notes, :resp, :src, :units, :t, :x)
    q = min(q, isempty(getfield(C, f)))
    (q == false) && return q
  end
  return q
end

function push!(S::NodalData, C::NodalChannel)
 for f in nodalfields
   push!(getfield(S, f), getfield(C, f))
 end
 S.data = hcat(S.data, C.x)
 S.n += 1
 resize!(S.x, S.n+1)
 S.x[S.n] = view(S.data, :, S.n)
 return nothing
end

function write(io::IO, S::NodalChannel)
  write(io, S.ox)                                                     # ox
  write(io, S.oy)                                                     # oy
  write(io, S.oz)                                                     # oz
  write(io, Int64(sizeof(S.id)))
  write(io, S.id)                                                     # id
  write(io, Int64(sizeof(S.name)))
  write(io, S.name)                                                   # name
  write(io, loctyp2code(S.loc))
  write(io, S.loc)                                                    # loc
  write(io, S.fs)                                                     # fs
  write(io, S.gain)                                                   # gain
  write(io, resptyp2code(S.resp))
  write(io, S.resp)                                                   # resp
  write(io, Int64(sizeof(S.units)))
  write(io, S.units)                                                  # units
  write(io, Int64(sizeof(S.src)))
  write(io, S.src)                                                    # src
  write_misc(io, S.misc)                                              # misc
  write_string_vec(io, S.notes)                                       # notes
  write(io, Int64(size(S.t,1)))
  write(io, S.t)                                                      # t
  write(io, typ2code(eltype(S.x)))
  write(io, Int64(length(S.x)))
  write(io, S.x)                                                      # x
  return nothing
end

read(io::IO, ::Type{NodalChannel}) = NodalChannel(
  fastread(io, Float64),                                                # ox
  fastread(io, Float64),                                                # oy
  fastread(io, Float64),                                                # oz
  String(fastread(io, fastread(io, Int64))),                            # id
  String(fastread(io, fastread(io, Int64))),                            # name
  read(io, code2loctyp(fastread(io))),                                  # loc
  fastread(io, Float64),                                                # fs
  fastread(io, Float64),                                                # gain
  read(io, code2resptyp(fastread(io))),                                 # resp
  String(fastread(io, fastread(io, Int64))),                            # units
  String(fastread(io, fastread(io, Int64))),                            # src
  read_misc(io, getfield(BUF, :buf)),                                   # misc
  read_string_vec(io, getfield(BUF, :buf)),                             # notes
  read!(io, Array{Int64, 2}(undef, fastread(io, Int64), 2)),              # t
  read!(io, Array{code2typ(read(io, UInt8)), 1}(undef, read(io, Int64))), # x
  )


text:
using Visualize, GeometryTypes, ModernGL
using Visualize: GLRasterization, Uniforms, VertexCS, get_texture_atlas, Sprite2
using Visualize: orthographicprojection, perspectiveprojection, lookat, normalmesh, ColorBuffer, Area
using Visualize.JLRasterization: JLCanvas, JLRasterizer
resolution = (1024, 1024)
canvas = JLCanvas((Area => SimpleRectangle(0, 0, resolution...),))

proj = orthographicprojection(SimpleRectangle(0, 0, resolution...), -10_000f0, 10_000f0)

uniforms = Uniforms(
    proj,
    Vec4f0(1, 0, 0, 1),
    Vec4f0(1, 0, 1, 1),
    Visualize.circle
)

N = 20
middle = Vec2f0(resolution) / 2f0
radius = min(resolution...) / 2f0
vertices = [(VertexCS(
    Vec2f0((sin(2pi * (i / N)) , cos(2pi * (i / N))) .* radius) .+ middle,
    Vec4f0(1, i/N, 0, 1),
    Vec2f0(40, 40)
),) for i = 1:N]
Visualize.visualize(canvas, x)

visualize(canvas, )
using FileIO, Images
save(joinpath(homedir(), "Desktop/test.png"), clamp01.(canvas[Visualize.ColorBuffer][1]))


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Polytope{K,Dim,T}

We say that a geometry is a K-polytope when it is a collection of "flat" sides
that constitue a `K`-dimensional subspace. They are called polygon and polyhedron
respectively for 2D (`K=2`) and 3D (`K=3`) subspaces, embedded in a `Dim`-dimensional
space. The parameter `K` is also known as the rank or parametric dimension of the
polytope: https://en.wikipedia.org/wiki/Abstract_polytope.

The term polytope expresses a particular combinatorial structure. A polyhedron,
for example, can be decomposed into faces. Each face can then be decomposed into
edges, and edges into vertices. Some conventions act as a mapping between vertices
and higher dimensional features (edges, faces, cells...), removing the need to
store all features.

Additionally, the following property must hold in order for a geometry to be considered
a polytope: the boundary of a (K+1)-polytope is a collection of K-polytopes, which may
have (K-1)-polytopes in common. See https://en.wikipedia.org/wiki/Polytope.

### Notes

- Type aliases are `Polygon`, `Polyhedron`.
"""
abstract type Polytope{K,Dim,T} <: Geometry{Dim,T} end

(::Type{PL})(vertices::Vararg{P}) where {PL<:Polytope,P<:Point} = PL(SVector(vertices))
(::Type{PL})(vertices::AbstractVector{TP}) where {PL<:Polytope,TP<:Tuple} = PL(Point.(vertices))
(::Type{PL})(vertices::Vararg{TP}) where {PL<:Polytope,TP<:Tuple} = PL(collect(vertices))

"""
    paramdim(polytope)

Return the parametric dimension or rank of the polytope.
"""
paramdim(::Type{<:Polytope{K}}) where {K} = K

"""
    vertices(polytope)

Return the vertices of a `polytope`.
"""
vertices(p::Polytope) = p.vertices

"""
    nvertices(polytope)

Return the number of vertices in the `polytope`.
"""
nvertices(p::Polytope) = length(vertices(p))

"""
    p1 == p2

Tells whether or not polytopes `p1` and `p2` are equal.
"""
==(p1::Polytope, p2::Polytope) = vertices(p1) == vertices(p2)

"""
    centroid(polytope)

Return the centroid of the `polytope`.
"""
centroid(p::Polytope) = Point(sum(coordinates, vertices(p)) / length(vertices(p)))

"""
    unique(polytope)

Return a new `polytope` without duplicate vertices.
"""
Base.unique(p::Polytope) = unique!(deepcopy(p))

"""
    unique!(polytope)

Remove duplicate vertices in `polytope`.
"""
function Base.unique!(::Polytope) end

function Base.show(io::IO, p::Polytope)
  kind = prettyname(p)
  vert = join(vertices(p), ", ")
  print(io, "$kind($vert)")
end

function Base.show(io::IO, ::MIME"text/plain", p::Polytope{K,Dim,T}) where {K,Dim,T}
  kind = prettyname(p)
  lines = ["  └─$v" for v in vertices(p)]
  println(io, "$kind{$Dim,$T}")
  print(io, join(lines, "\n"))
end

prettyname(p::Polytope) = prettyname(typeof(p))
function prettyname(PL::Type{<:Polytope})
  n = string(PL)
  i = findfirst('{', n)
  isnothing(i) ? n : n[1:i-1]
end

# --------
# POLYGON
# --------

"""
    Polygon{Dim,T}

A polygon is a 2-polytope, i.e. a polytope with parametric dimension 2.
"""
const Polygon = Polytope{2}

"""
    area(polygon)

Return the area of the `polygon`.
"""
area(p::Polygon) = measure(p)

"""
    chains(polygon)

Return the outer and inner chains of the polygon.
"""
function chains(::Polygon) end

"""
    hasholes(polygon)

Tells whether or not the `polygon` contains holes.
"""
function hasholes(::Polygon) end

"""
    issimple(polygon)

Tells whether or not the `polygon` is simple.
See [https://en.wikipedia.org/wiki/Simple_polygon]
(https://en.wikipedia.org/wiki/Simple_polygon).
"""
issimple(p::Polygon) = issimple(typeof(p))

"""
    windingnumber(point, polygon)

Winding number of `point` with respect to the `polygon`.
"""
function windingnumber(::Point, ::Polygon) end

"""
    orientation(polygon)

Returns the orientation of the rings of the `polygon`
as either counter-clockwise (CCW) or clockwise (CW).
"""
orientation(p::Polygon) = orientation(p, WindingOrientation())

function orientation(p::Polygon, algo)
  o = [orientation(c, algo) for c in chains(p)]
  hasholes(p) ? o : first(o)
end

"""
    boundary(polygon)

Returns the boundary of the `polygon`.
"""
boundary(p::Polygon) = hasholes(p) ? Multi(chains(p)) : first(chains(p))

"""
    isconvex(polygon)

Tells whether or not the `polygon` is convex.
"""
isconvex(p::Polygon{Dim,T}) where {Dim,T} =
  issimple(p) && all(≤(T(π)), innerangles(boundary(p)))

"""
    bridge(polygon; width=0)

Transform `polygon` with holes into a single outer chain
via bridges of given `width` as described in Held 1998.
Return the outer chain and a vector with pairs of indices
for duplicate vertices. These indices can be used to undo
the bridges.

## References

* Held. 1998. [FIST: Fast Industrial-Strength Triangulation of Polygons]
  (https://link.springer.com/article/10.1007/s00453-001-0028-4)
"""
function bridge(p::Polygon{Dim,T}; width=zero(T)) where {Dim,T}
  if hasholes(p)
    bridge(chains(p), width=width)
  else
    first(chains(p)), []
  end
end

# -----------
# POLYHEDRON
# -----------

"""
    Polyhedron{Dim,T}

A polyhedron is a 3-polytope, i.e. a polytope with parametric dimension 3.
"""
const Polyhedron = Polytope{3}

"""
   volume(polyhedron)

Return the volume of the `polyhedron`.
"""
volume(p::Polyhedron) = measure(p)

# ----------------
# IMPLEMENTATIONS
# ----------------

include("polytopes/segment.jl")
include("polytopes/ngon.jl")
include("polytopes/chain.jl")
include("polytopes/polyarea.jl")
include("polytopes/pyramid.jl")
include("polytopes/tetrahedron.jl")
include("polytopes/hexahedron.jl")