text:
#!/bin/env julia

using CLI, JLD, Distributions, HypothesisTests

immutable Region
	start::Int32
	stop::Int32
end

# Print comment lines to stdout and return
# first non-comment line
function skip_to_header(vcf_file)
	while true
		line = readline(vcf_file);
		if !startswith( line, "#") return line; end
		println(line);
	end
end

drop_chr_prefix(chr) = startswith(chr, "chr") ? chr[4:end] : chr

is_protein_altering(effect) = ismatch(
	r"Missense|[Ff]rameshift|Stopgain|Stoploss|Splice", effect);

function mutation_class(ref, alt)
	if ref == "C"
		if alt == "A"; return 1; end
		if alt == "G"; return 2; end
		if alt == "T"; return 3; end
	elseif ref == "G"
		if alt == "T"; return 1; end
		if alt == "C"; return 2; end
		if alt == "A"; return 3; end
	elseif ref == "T"
		if alt == "A"; return 4; end
		if alt == "C"; return 5; end
		if alt == "G"; return 6; end
	elseif ref == "A"
		if alt == "T"; return 4; end
		if alt == "G"; return 5; end
		if alt == "C"; return 6; end
	end

	if length(ref) > 1 || ref == "-" || length(alt) > 1 || alt == "-"
		return 0
	else
		error("Mutation class could not be defined.")
	end
end

function somatic(vcf_file::IO, samples_file::IO; alt_reads=3, alt_frac=0.1,
	test_ref_ratio=5, test_bg_ratio=20, test_bg_alpha=1.0, ref_reads=0,
	keep_all=false, min_mapq=10, min_sidedness=6, keep_old_stars::Bool=false)

	header = skip_to_header(vcf_file); println(header)
	header_cols = split(header, '\t')
	notes_col = findfirst(header_cols .== "NOTES")
	alt_col = findfirst(header_cols .== "ALT")
	ref_col = findfirst(header_cols .== "REF")

	samples = header_cols[notes_col+1:end];
	n_samples = length(samples);

	# Read samples file and separate the row header
	d = readdlm(samples_file, '\t');
	sample_file_header = d[1, :][:];
	d = d[2:end, :];

	sample_test_col = findfirst(sample_file_header .== "TEST")
	sample_ref_col = findfirst(sample_file_header .== "REF")
	assert(sample_test_col != 0 && sample_ref_col != 0)

	# Check that all samples named in the design file are present in the VCF
	for s in unique(d[:, [sample_test_col, sample_ref_col]])
		if s != "" && !(s in samples)
			error("Sample $s is listed in the design file but missing from VCF.")
		end
	end

	# Background error is estimated based on paired *and* unpaired reference
	# samples.
	ref_samples = filter(s -> s != "", unique(d[:, sample_ref_col]))
	background = indexin(ref_samples, samples)

	# Get indices of test and ref sample rows in the samples file
	test_rows = d[:, sample_test_col] .!= ""
	test = indexin(d[test_rows, sample_test_col], samples)
	n_tests = length(test);
	ref = indexin(d[test_rows, sample_ref_col], samples)
	info("Analyzing $n_tests tumor-normal pairs...")

	if !isempty(findin(background, test))
		warn("Some background reference samples also used as test samples.")
	end
	info("Calculating background error rate based on $(length(background)) samples...")

	# Read alt_frac from design file or command line parameter.
	alt_frac_col = findfirst(sample_file_header .== "ALT_FRAC")
	if alt_frac_col == 0
		info("Using a mutant allele fraction threshold of $(alt_frac * 100)% for all samples...")
		alt_frac = ones(n_tests) * alt_frac
	else
		alt_frac = convert(Array{Float64}, d[test_rows, alt_frac_col])
	end

	# TODO: Allow sample-specific values to be specified in design file for
	# alt_reads, test_ref_ratio and test_bg_ratio as well.

	for line in eachline(vcf_file)

		cols = split(line, '\t')
		alt = zeros(Int32, n_samples);
		total = zeros(Int32, n_samples);
		# CALL2 support
		mapq = zeros(Int32, n_samples);
		baseq = zeros(Int32, n_samples);
		sidedness = zeros(Int32, n_samples);

		call2_format = false
		is_indel = length(cols[alt_col]) != length(cols[ref_col])

		for s in 1:n_samples
			col = cols[notes_col + s]
			gt = split(cols[notes_col + s], ':')
			alt[s] = parse(Int32, gt[1])
			total[s] = parse(Int32, gt[2])
			# CALL2 support
			if length(gt) >= 5
				mapq[s] = parse(Int32, gt[3])
				baseq[s] = parse(Int32, gt[4])
				sidedness[s] = parse(Int32, gt[5])
				call2_format = true
			end
		end

		fracs = alt ./ total
		bg_alt = sum(alt[background])
		bg_total = sum(total[background])
		avg_ref_frac = bg_alt / bg_total

		somatic = falses(n_samples)
		for k in 1:n_tests
			t = test[k]; # Index of test sample col
			r = ref[k];  # Index of ref sample col

			# Perform checks to determine somatic status of mutation
			if fracs[t] < alt_frac[k]; continue; end
			if r != 0 && fracs[t] < test_ref_ratio * fracs[r]; continue; end
			if fracs[t] < test_bg_ratio * avg_ref_frac; continue; end
			if alt[t] < alt_reads; continue; end
			if r != 0 && total[r] < ref_reads; continue; end
			if call2_format
				if !is_indel && mapq[t] < min_mapq; continue; end
				if sidedness[t] < min_sidedness; continue; end
			end

			#if test_bg_alpha < 1.0
			#	if pvalue(ChisqTest([alt[t] total[t] - alt[t]; bg_alt bg_total - bg_alt])) > test_ref_alpha; continue; end
			#end
			#chisq_p = pvalue(ChisqTest(
			#	[alt[t] total[t] - alt[t]; bg_alt bg_total - bg_alt]))
			#baseq[t] = round(Int32, min(-10 * log10(chisq_p), 100))

			somatic[t] = true
		end

		if keep_old_stars
			somatic |= map(c -> endswith(c, '*'), cols[notes_col+1:end])
		end

		if !any(somatic) && keep_all == false; continue; end

		print(join(cols[1:notes_col], '\t'))
		@printf("BG: %.3f%%. ", avg_ref_frac * 100)
		for s in 1:n_samples
			if call2_format
				print("\t$(alt[s]):$(total[s]):$(mapq[s]):$(baseq[s]):$(sidedness[s])")
			else
				print("\t$(alt[s]):$(total[s])")
			end

			if somatic[s]; print(":*"); end
		end
		println()
	end
end

function germline(vcf_file::IO, ref_regexps...; alt_reads=5, alt_frac=0.15, bg_ratio=20)
	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	samples = split(line, '\t')[sample_col:end]
	germline = falses(samples)
	for regex in ref_regexps; germline .|= ismatch.(Regex(regex), samples); end
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		gtypes = [split(c, ':') for c in cols[sample_col:end]]
		total = [parse(Int32, gt[2]) for gt in gtypes]
		alt = [parse(Int32, gt[1]) for gt in gtypes]
		fracs = alt ./ total
		alt_gt = (fracs .>= alt_frac) .& (alt .>= alt_reads) .& germline
		if !any(alt_gt); continue; end

		# There might be non-reference samples that are not in alt_gt because
		# of the alt_reads threshold. We don't want to include these
		# samples in the background, so we only use alt_frac when picking
		# background samples.
		avg_alt = sum(alt[alt_gt]) / sum(total[alt_gt])
		bg = (fracs .< alt_frac) .& germline
		avg_ref = sum(alt[bg]) / sum(total[bg])
		if avg_alt / avg_ref < bg_ratio; continue; end

		print(join(cols[1:sample_col-1], '\t'))
		for s in 1:length(total)
			print("\t$(alt[s]):$(total[s])")
			if alt_gt[s]; print(":*"); end
		end
		print('\n')
	end
end

function above_background(vcf_file::IO, ref_regex...; alt_reads=3, alt_frac=0.1, bg_ratio=20)
	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	samples = split(line, '\t')[sample_col:end]

	bg = falses(samples)
	for regex in ref_regex; bg .|= ismatch.(Regex(regex), samples); end
	info("Calculating background based on $(sum(bg)) samples...")

	for line in eachline(vcf_file)
		cols = split(line, '\t')
		gtypes = map(c -> split(c, ':'), cols[sample_col:end])
		total = map(gt -> parse(Int, gt[2]), gtypes)
		alt = map(gt -> parse(Int, gt[1]), gtypes)
		fracs = alt ./ total
		avg_bg_frac = sum(alt[bg]) / sum(total[bg])
		alt_gt = (alt .>= alt_reads) .&
			(fracs .>= max(alt_frac, avg_bg_frac * bg_ratio))
		if !any(alt_gt); continue; end

		print(join(cols[1:sample_col-1], '\t'))
		for s in 1:length(total)
			print("\t$(alt[s]):$(total[s])")
			if alt_gt[s]; print(":*"); end
		end
		println()
	end
end

function heterozygous_snps(vcf_file::IO, ref_regex...; min_depth=30, min_frac=0.3, bg_fraction=0.1)
	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	samples = split(line, '\t')[sample_col:end]; S = length(samples)

	ref = falses(samples)
	for regex in ref_regex; ref .|= ismatch.(Regex(regex), samples); end
	info("Searching for heterozygous SNPs in $(sum(ref)) germline samples...")

	for line in eachline(vcf_file)
		cols = split(line, '\t')
		alt = zeros(Int, S); total = zeros(Int, S); starred = falses(S);
		unclear = 0
		for (s, col) in enumerate(cols[sample_col:end])
			parts = split(col, ':')
			alt[s] = parse(Int, parts[1])
			total[s] = parse(Int, parts[2])
			if ref[s] == false; continue; end
			if total[s] < min_depth; continue; end

			if alt[s] / total[s] < min_frac
				error_lower = Binomial(total[s], 0.01)
				if 1 - cdf(error_lower, alt[s] - 1) <= 0.05; unclear += 1; end
			elseif alt[s] / total[s] > 1 - min_frac
				error_upper = Binomial(total[s], 0.99)
				if cdf(error_upper, alt[s]) <= 0.05; unclear += 1; end
			end

			if unclear / S >= bg_fraction
				starred[:] = false; break;
			end

			hetz = Binomial(total[s], 0.5)
			if 0.05 <= cdf(hetz, alt[s]) <= 0.95; starred[s] = true; end
		end

		if !any(starred); continue; end

		print(join(cols[1:sample_col-1], '\t'))
		for s in 1:S
			print("\t$(alt[s]):$(total[s])")
			if starred[s]; print(":*"); end
		end
		print('\n')
	end
end

annovar_valid_funcs = Set(["exonic", "splicing", "intergenic", "intronic",
	"upstream", "downstream", "UTR5", "UTR3", "ncRNA_exonic", "ncRNA_intronic",
	"ncRNA_splicing", ""])

function translate_annovar_effect(gene, func, exonic_func, aa_change)
	for f in split(func, ';')
		if !in(f, annovar_valid_funcs)
			warn("Unrecognized variant effect '$f'.")
		end
	end

	effects = []
	if contains(func, "splicing"); push!(effects, "Splice site"); end
	if contains(func, "exonic")
		details = unique(map(m -> m.captures[1],
			eachmatch(r":(p\..+?)(,|$)", aa_change)))
		if exonic_func == "synonymous SNV"
			push!(effects, "Synonymous $(join(details, ", "))")
		elseif exonic_func == "nonsynonymous SNV"
			push!(effects, "Missense $(join(details, ", "))")
		elseif exonic_func == "stopgain"
			push!(effects, "Stopgain $(join(details, ", "))")
		elseif exonic_func == "stoploss"
			push!(effects, "Stoploss $(join(details, ", "))")
		elseif contains(exonic_func, "nonframeshift")
			push!(effects, "Non-frameshift indel $(join(details, ", "))")
		elseif contains(exonic_func, "frameshift")
			push!(effects, "Frameshift $(join(details, ", "))")
		elseif contains(func, "ncRNA_exonic")
			push!(effects, "Exonic (ncRNA)")
		elseif exonic_func == "unknown"
			push!(effects, "Exonic (unknown)")
		else
			error("Unrecognized effect: $([func, exonic_func])")
		end
	end
	if contains(func, "UTR3"); push!(effects, "3'-UTR"); end
	if contains(func, "UTR5"); push!(effects, "5'-UTR"); end
	if contains(func, "upstream"); push!(effects, "Upstream"); end
	if contains(func, "downstream"); push!(effects, "Downstream"); end
	if contains(func, "intronic"); push!(effects, "Intronic"); end
	if func == "intergenic"; push!(effects, "Intergenic"); gene = ""; end

	return (gene, join(effects, ". "))
end

function predict_effect(vcf_file::IO; genome="~/tools/annovar-2016-02-01/humandb/hg38")
	genome_version = basename(genome)
	humandb_dir = expanduser(dirname(genome))
	tmp_prefix = "annovar-$(hex(rand(0:2^16-1), 4))"
	original_path = "$(tmp_prefix).original"
	reformatted_path = "$(tmp_prefix).avinput"

	# Reformat the input file to an ANNOVAR-compatible format. In particular,
	# indels must be changed so they use '-' as the ref or alt allele.
	# Since the input file could be a pipe, we must write the original VCF
	# into a separate file that will be used to build the final annotated VCF.
	original_file = open(original_path, "w")
	reformatted_file = open(reformatted_path, "w")
	headers = split(skip_to_header(vcf_file), '\t')

	for line in eachline(vcf_file)

		write(original_file, line, '\n')
		c = split(line, '\t')
		start = parse(Int, c[2]); stop = start; ref = c[3]; alt = c[4]

		if length(ref) == 1 && length(alt) > 1 && ref[1] == alt[1]
			# Simple insertion
			ref = "-"; alt = alt[2:end];
		elseif length(ref) > 1 && length(alt) == 1 && ref[1] == alt[1]
			# Simple deletion
			ref = ref[2:end]; alt = "-"; start += 1;
		elseif length(ref) > 1 && length(alt) > 1
			# Block substitution
		end
		stop = start + length(ref) - 1
		write(reformatted_file, "$(c[1])\t$(start)\t$(stop)\t$(ref)\t$(alt)\n")
	end
	close(original_file); close(reformatted_file); close(vcf_file)

	run(pipeline(`table_annovar.pl $(reformatted_path) $(humandb_dir) --buildver $(genome_version) --remove --outfile $(tmp_prefix) --operation g --protocol refGene`, stderr=DevNull))

	anno_path = "$(tmp_prefix).$(genome_version)_multianno.txt"
	original_file = open(original_path);
	anno_file = open(anno_path); readline(anno_file)
	print("CHROM\tPOSITION\tREF\tALT\tGENE\tEFFECT\tNOTES\t")
	println(join(headers[6:end], '\t'))
	for line in eachline(anno_file)
		c = split(line, '\t')
		if length(c) >= 7
			gene, effect = translate_annovar_effect(c[7], c[6],
				length(c) >= 9 ? c[9] : "", length(c) >= 10 ? c[10] : "")
		else
			gene = ""; effect = "";
		end

		orig_line = readline(original_file)
		tab = search(orig_line, '\t')
		for k in 1:3; tab = search(orig_line, '\t', tab + 1); end
		print(orig_line[1:tab])
		print("$(gene)\t$(effect)\t")
		println(orig_line[tab+1:end])
	end
	close(anno_file); rm(original_path); rm(reformatted_path); rm(anno_path);
end

type TextAnnotations
	position::Vector{Int32}
	change::Vector{Int16}
	annotation::Vector{String}
	TextAnnotations() = new([], [], [])
end

type NumericAnnotations
	position::Vector{Int32}
	change::Vector{Int16}
	annotation::Vector{Float32}
	NumericAnnotations() = new([], [], [])
end

base_num = Dict{String, Int16}("A" => 1, "C" => 2, "G" => 3, "T" => 4)
encode_ref_alt(ref, alt) = get(base_num, ref, 0) * 16 + get(base_num, alt, 0)

function build_annotation_database(anno_path, db_path; prefix="",
	format="text")
	db = jldopen(db_path, "w", compress=true)
	anno_file = zopen(anno_path)
	if format == "text"
		build_annotation_database_text(anno_file, db)
	elseif format == "numeric"
		build_annotation_database_numeric(anno_file, db)
	elseif format == "blank"
		build_annotation_database_blank(anno_file, db)
	else
		error("Unrecognized annotation database format '$format' requested.")
	end
	write(db, "prefix", prefix)
	write(db, "format", format)
	close(db)
end

function build_annotation_database_blank(anno_file, db)
	chr = ""; anno = TextAnnotations()
	for line in eachline(anno_file)
		c = split(rstrip(line, '\t'))
		if c[1] != chr
			if chr != ""
				assert(issorted(anno.position))
				write(db, chr, anno)
			end
			anno = TextAnnotations(); chr = String(c[1])
			println(STDERR, "Chromosome $(chr)...")
		end
		push!(anno.position, parse(Int32, c[2]))
		push!(anno.change, encode_ref_alt(c[3], c[4]))
	end
	write(db, chr, anno)
end

function build_annotation_database_text(anno_file, db)
	chr = ""; anno = TextAnnotations()
	for line in eachline(anno_file)
		c = split(rstrip(line, '\t'))
		if c[1] != chr
			if chr != ""
				assert(issorted(anno.position))
				write(db, chr, anno)
			end
			anno = TextAnnotations(); chr = String(c[1])
			println(STDERR, "Chromosome $(chr)...")
		end
		push!(anno.position, parse(Int32, c[2]))
		push!(anno.change, encode_ref_alt(c[3], c[4]))
		push!(anno.annotation, c[5])
	end
	write(db, chr, anno)
end

function build_annotation_database_numeric(anno_file, db)
	chr = ""; anno = NumericAnnotations()
	for line in eachline(anno_file)
		c = split(rstrip(line, '\t'))
		if c[1] != chr
			if chr != ""
				assert(issorted(anno.position))
				write(db, chr, anno)
			end
			anno = NumericAnnotations(); chr = String(c[1])
			println(STDERR, "Chromosome $(chr)...")
		end
		push!(anno.position, parse(Int32, c[2]))
		push!(anno.change, encode_ref_alt(c[3], c[4]))
		push!(anno.annotation, parse(Float32, c[5]))
	end
	write(db, chr, anno)
end

function row_to_key( row::String)
	# Firstly, compare non-numeric parts (case-insensitive)
	# Next, compare numeric parts as integer
	# Next, compare position as integer
	chrom, pos = split(row, '\t')[1:2]
	return (filter(x -> !isnumber(x), lowercase(chrom)),
		    parse(Int32, filter(x -> isnumber(x), '0'*chrom)),
		    parse(Int32, pos))
end

# Sorts and prints VCF file rows by their chromosomal position
# Does not allow comment lines anywhere else than the beginning
function sort_by_position(vcf_file::IO)

	println( skip_to_header(vcf_file));
	rows = readlines(vcf_file)

	sort!(rows, by=x -> row_to_key(x))
	for row in rows println(row); end
end

function annotate(vcf_file::IO, db_path)
	db = jldopen(db_path, "r")
	prefix = read(db, "prefix")
	format = read(db, "format")
	headers = split(skip_to_header(vcf_file), '\t')
	notes_col = findfirst(headers .== "NOTES")
	println(join(headers, '\t'))

	chr = ""; anno = nothing
	for line in eachline(vcf_file)

		c = split(line, '\t')
		if drop_chr_prefix(c[1]) != chr
			chr = String(drop_chr_prefix(c[1]))
			if any(x -> chr == x, names(db))   # in() does not work
				anno = read(db, chr)
			else
				warn("No annotations found for chromosome $(chr).")
				anno = NumericAnnotations()
			end
		end
		matches = searchsorted(anno.position, parse(Int32, c[2]))
		if !isempty(matches)
			change = encode_ref_alt(c[3], c[4])
			for m in matches
				if anno.change[m] == change
					if format == "numeric"
						# See https://github.com/JuliaLang/julia/issues/14331
						c[notes_col] *= rstrip(@sprintf("%s:%g", prefix,
							anno.annotation[m])) * ". "
					elseif format == "text"
						c[notes_col] *= "$(prefix):$(anno.annotation[m]). "
					else
						c[notes_col] *= "$(prefix). "
					end
				end
			end
		end
		println(join(c, '\t'))
	end
end

# Filter VCF file lines based on if they are protein altering or not
# ARG: invert: Set to false to output protein altering mutations"
#              Set to true to output non protein altering mutations"
function protein_altering(vcf_file::IO; invert=false)
	line = skip_to_header(vcf_file); println(line)
	func_col = findfirst(split(line, '\t') .== "EFFECT")
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		alters = ismatch(r"Missense|rameshift|Stopgain|Stoploss|Splice",
			cols[func_col])
		if !alters == invert; println(line); end
	end
end

function select_samples(invert::Bool, vcf_file::IO, regexps...)
	headers = split(skip_to_header(vcf_file), '\t')
	sample_col = findfirst(headers .== "NOTES") + 1
	keep = falses(length(headers))
	keep[1:sample_col-1] = true
	for col in sample_col:length(headers)
		keep[col] = any(rx -> ismatch(Regex(rx), headers[col]), regexps)
	end
	if invert; keep[sample_col:end] = .!keep[sample_col:end]; end
	println(join(headers[keep], '\t'))
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		println(join(cols[keep], '\t'))
	end
end
keep_samples(vcf_file::IO, regexps...) = select_samples(false, vcf_file, regexps...)
discard_samples(vcf_file::IO, regexps...) = select_samples(true, vcf_file, regexps...)

function discard_shallow(vcf_file::IO, min_median_depth::Float64)
	line = skip_to_header(vcf_file); println(line);
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	for line in eachline(vcf_file)
		total = map(c -> parse(Int, split(c, ':')[2]),
			split(line, '\t')[sample_col:end])
		if median(total) >= min_median_depth; println(line); end
	end
end

function discard_blacklisted(vcf_file::IO, blacklist::IO)
	d = readdlm(blacklist, '\t', String)
	blacklist = Set(map(r -> d[r, 1:4][:], 1:size(d, 1)))

	line = skip_to_header(vcf_file); println(line);
	for line in eachline(vcf_file)
		c = split(line, '\t')
		if in(c[1:4], blacklist); continue; end
		println(line)
	end
end

function discard_indels(vcf_file::IO)
	line = skip_to_header(vcf_file); println(line);
	for line in eachline(vcf_file)
		c = split(line, '\t')
		if length(c[3]) == 1 && length(c[4]) == 1 && c[3] != "-" && c[4] != "-"
			println(line)
		end
	end
end

function discard_if_frequency_above(vcf_file::IO, db_path, threshold::Float64)
	db = jldopen(db_path, "r")
	println(skip_to_header(vcf_file))
	chr = ""; anno = nothing
	for line in eachline(vcf_file)
		c = split(line, '\t')
		if drop_chr_prefix(c[1]) != chr
			chr = String(drop_chr_prefix(c[1])); anno = read(db, chr)
		end
		matches = searchsorted(anno.position, parse(Int32, c[2]))
		if !isempty(matches)
			change = encode_ref_alt(c[3], c[4])
			if any(m -> anno.change[m] == change &&
				anno.annotation[m] >= threshold, matches)
				continue;
			end
		end
		println(join(c, '\t'))
	end
end

function inside(vcf_file::IO, region)
	m = match(r"([a-zA-Z0-9]+):(\d+)-(\d+)", region)
	if m == nothing
		if isfile(region) || isfifo(region)
			inside_bed_regions(vcf_file, region)
			return
		else
			error("Invalid region specified.")
		end
	end

	chr = m.captures[1]
	start = parse(Int, m.captures[2])
	stop = parse(Int, m.captures[3])

	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	for line in eachline(vcf_file)
		c = split(line, '\t'); pos = parse(Int, c[2])
		if c[1] == chr && start <= pos <= stop; println(line); end
	end
end

function inside_bed_regions(vcf_file::IO, regions_path)
	regions = Dict{String, Array{Region}}()
	for line in eachline(open(regions_path))
		c = split(rstrip(line), '\t')
		push!(get!(() -> Array{Region}(0), regions, c[1]),
			Region(parse(Int, c[2])+1, parse(Int, c[3])))
	end

	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1 + 1
	for line in eachline(vcf_file)
		c = split(line, '\t'); pos = parse(Int, c[2])
		chr_regions = get(() -> Array{Region}(0), regions, c[1])
		if any(r -> r.start <= pos <= r.stop, chr_regions)
			println(line)
		end
	end
end

function discard_contingent(vcf_path, alpha::Float64)
	vcf_file = zopen(vcf_path)
	line = skip_to_header(vcf_file); print(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	for line in eachline(vcf_file)
		gtypes = map(c -> split(c, ':'), split(line, '\t')[sample_col:end])
		total = map(gt -> parse(Int, gt[2]), gtypes)
		alt = map(gt -> parse(Int, gt[1]), gtypes)
		valid = total .>= 1
		p = any(valid) ? pvalue(ChisqTest([alt[valid] total[valid] - alt[valid]])) : 1
		if alpha == 1
			print(line[1:end-1]); print("\t$(p)\n")
		elseif p <= alpha
			print(line)
		end
	end
end

# Adds a warning to the Notes column if the mappability in the reference
# genome is below a specified treshold (-> low mappability)
function mappability(vcf_file::IO, track_path::String, threshold::Float64)
	tracks = jldopen(track_path, "r");
	headers = split(skip_to_header(vcf_file), '\t')
	notes_col = findfirst(headers .== "NOTES")
	println(join(headers, '\t'))

	# We assume that the VCF file is sorted by chromosome.
	curr_chr = ""; track = zeros(UInt8, 0);
	for line in eachline(vcf_file)
		c = split(line, '\t')
		if c[1] != curr_chr
			curr_chr = String(c[1])
			track = read(tracks, curr_chr)
		end
		mappability = track[parse(Int, c[2])]
		if mappability < threshold
			c[notes_col] = "$(c[notes_col])Mappability < $(threshold)%. "
		end
		println(join(c, '\t'))
	end
end

# Adds a warning to the NOTES column if there are indels within
# dist distance of the alt read (line in VCF file).
#
# ARG: dist : largest distance that is considered
#             being "nearby"
function nearby_indels(vcf_file::IO; dist=10)

	line = skip_to_header(vcf_file);
	println(line)
	headers = split(line, '\t')
	notes_col = findfirst(split(line, '\t') .== "NOTES")
	warn_msg = "INDEL within $dist bases. "

	#notes_col = findfirst(split(skip_to_header(vcf_file), '\t') .== "NOTES")

	prev_lines = Array{Array{String}}(0)

	for line in eachline(vcf_file)

		if startswith(line, "#") println(line); continue; end
		c = split(line, '\t')
		pos = parse(Int, c[2])
		starred = map(x -> endswith(x, '*'), c[notes_col+1:end])

		# Forget old lines that are so far behind they cannot be affected.
		keep = trues(length(prev_lines))
		for (k, prev_line) in enumerate(prev_lines)
			# On different chromosome or too distant
			if c[1] != prev_line[1] || parse(Int, prev_line[2]) < pos - dist
				keep[k] = false
				println(join(prev_line, '\t'))
			end
		end
		prev_lines = prev_lines[keep]

		# Add warning to old lines if current line is an indel.
		if length(c[3]) != length(c[4])
			for prev_line in prev_lines
				if any(s -> starred[s] && endswith(prev_line[notes_col+s], '*'), 1:length(starred))
					prev_line[notes_col] *= warn_msg
				end
			end
		end

		# Add warning to current line if any old lines were indels.
		for prev_line in prev_lines
			if length(prev_line[3]) != length(prev_line[4]) &&
				any(s -> starred[s] && endswith(prev_line[notes_col+s], '*'), 1:length(starred))
				c[notes_col] *= warn_msg
				break
			end
		end

		push!(prev_lines, c)
	end

	for prev_line in prev_lines; println(join(prev_line, '\t')); end
end

function statistics(vcf_path)
	vcf_file = zopen(vcf_path)
	headers = split(rstrip(skip_to_header(vcf_file)), '\t')
	sample_col = findfirst(headers .== "NOTES") + 1
	samples = headers[sample_col:end]; S = length(samples)
	total_subs = zeros(Int, S)
	possible_subs = zeros(Int, S, 6)
	total_indels = zeros(Int, S)

	for line in eachline(vcf_file)
		c = split(rstrip(line), '\t')
		alt = map(gt -> endswith(gt, '*'), c[sample_col:end])
		mut_class = mutation_class(c[3], c[4])
		if mut_class > 0
			total_subs += alt
			possible_subs[:, mut_class] += alt
		elseif mut_class == 0
			total_indels += alt
		end
	end

	println("SAMPLE\tTOTAL SUBSTITUTIONS\tTOTAL INDELS\tC>A\tC>G\tC>T\tT>A\tT>C\tT>G")
	for s in 1:S
		print("$(samples[s])\t$(total_subs[s])\t$(total_indels[s])")
		for k in 1:6; print("\t$(possible_subs[s, k])"); end
		println()
	end
end

# CALL2 output format: "alt_reads:total_reads:mapq:baseq:sidedness:significant"
function fractions(vcf_file::IO)

	line = skip_to_header(vcf_file); println(line)
	sample_col = findfirst(split(line, '\t') .== "NOTES") + 1
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		gtypes = [split(c, ':') for c in cols][sample_col:end]

		total = [parse(Int32, gt[2]) for gt in gtypes]
		alt = [parse(Int32, gt[1]) for gt in gtypes]

		print(join(cols[1:(sample_col-1)], '\t'))

		for s in 1:length(gtypes)
			frac = alt[s] / total[s]
			@printf("\t%.1f%% (%d)", isnan(frac) ? 0.0 : frac * 100, total[s])
			#CALL2 support
			if length(gtypes[s]) >= 5 && parse( Int32, gtypes[s][ 3]) > 0; @printf(" [mq:%s,bq:%s,sd:%s]", gtypes[s][3], gtypes[s][4], gtypes[s][5] ); end
			if length(gtypes[s]) >= 3 && gtypes[s][ length(gtypes[s])] == "*"; print(" *"); end
		end
		println()
	end
end

function error_rate(bam_path, genome_path; min_mapq=10)
	bam_path = expanduser(bam_path)
	genome_path = expanduser(genome_path)
	spileup = expanduser("~/tools/pypette/compiled/spileup")

	# 2D histogram, mutant read count on rows, mutant read fraction on columns
	# Histogram bins are left-exclusive, right-inclusive.
	read_thresholds = 0:100
	frac_thresholds = 0:0.001:0.10
	histogram = zeros(Int, length(read_thresholds), length(frac_thresholds))

	for line in eachline(pipeline(`samtools mpileup -d 100000 -A -x -R -sB -q$(min_mapq) -f $(genome_path) $(bam_path)`, `$(spileup) 0 $(min_mapq)`))
		cols = split(rstrip(line), '\t')
		if length(cols) < 4; continue; end
		#if cols[3] == "N"; print(line); continue; end
		assert(length(cols) == 4)
		assert(cols[3] != "N")

		pileup = cols[4]
		if pileup == ""; continue; end
		t = split(pileup, ' ')

		total_reads = 0
		alt_reads = 0

		for a in 1:3:length(t)
		 	count = parse(Int, t[a+1])
		 	total_reads += count
		 	if t[a] != "."; alt_reads += count; end
		end

		alt_frac = alt_reads / total_reads
		read_bin = Int(min(alt_reads + 1, 101))
		frac_bin = min(ceil(Int, alt_frac / 0.001) + 1, 101)
		histogram[read_bin, frac_bin] += 1

		# allele_reads = zeros(Int32, max_alleles)
		# alleles = Array{String}(0)

		# pileup = cols[4]
		# if pileup == ""; continue; end
		# t = split(pileup, ' ')
		# for a in 1:3:length(t)
		# 	count = int(t[a+1])
		# 	total_reads[s] += count
		# 	if t[a] == "."; continue; end
		# 	idx = findfirst(allele -> allele == t[a], alleles)
		# 	if idx == 0 && length(alleles) >= max_alleles; continue; end
		# 	if idx == 0; push!(alleles, t[a]); idx = length(alleles); end
		# 	allele_reads[s, idx] = count
		# end

		# for a in 1:length(alleles)
		# end
	end

	writedlm(STDOUT, histogram)
end

function group_by_snp_profile(hetz_snp_vcf::IO)
	headers = split(rstrip(skip_to_header(hetz_snp_vcf)), '\t')
	sample_col = findfirst(headers .== "NOTES") + 1
	samples = headers[sample_col:end]; S = length(samples);

	# The matrix stores genotypes for every sample:
	# 0 = unknown, 1 = homz ref, 2 = hetz, 3 = homz alt
	println("Constructing genotype matrix...")
	genotypes = zeros(Int8, 1_000_000, S); r = 0
	for line in eachline(hetz_snp_vcf)
		c = split(rstrip(line, '\n'), '\t')
		r += 1
		for (s, cell) in enumerate(c[sample_col:end])
			parts = split(cell, ':')
			alt = parse(Int, parts[1]); total = parse(Int, parts[2])
			if total < 30; genotypes[r, s] = 0; continue; end
			b = Binomial(total, 0.5)
			if 0.05 <= cdf(b, alt) <= 0.95
				genotypes[r, s] = 2
			elseif alt / total >= 0.8
				genotypes[r, s] = 3
			elseif alt / total <= 0.2
				genotypes[r, s] = 1
			end
		end
	end
	genotypes = genotypes[1:r, :]

	for i in 1:S
		any_matches = false
		for j in 1:S
			if j == i; continue; end
			matched = 0; total = 0;
			for r in 1:size(genotypes, 1)
				if genotypes[r, i] == 0 || genotypes[r, j] == 0; continue; end
				total += 1
				matched += (genotypes[r, i] == genotypes[r, j])
			end
			if matched / total >= 0.95 && total >= 50
				any_matches = true
				if j > i      # Don't print pairs twice
					println("$(samples[i]) <-> $(samples[j]): $matched / $total matched")
				end
			end
		end

		if any_matches == false
			println("$(samples[i]) does not match with any other samples.")
		end
	end
end

function mutation_rate(vcf_file::IO, coverage_histogram_dir; alt_reads=8)
	vcf = read_vcf(vcf_file); S = length(vcf.sample);
	human_chr = vcat(map(x -> "chr$(x)", 1:22), ["chrX", "chrY"]);
	println("SAMPLE\tPROTEIN ALTERING\tSILENT\tPROTEIN ALTERING (RATE)\tSILENT (RATE)")
	for s in 1:S
		info("Analyzing $(vcf.sample[s])...")
		protein_altering = 0; silent = 0;
		protein_altering_rate = 0; silent_rate = 0;
		histogram = zeros(5000)   # Bins are 0, 1, 2, ..., 4999
		genome_len = 0
		histogram_file = "$(coverage_histogram_dir)/$(vcf.sample[s]).tsv"
		if !isfile(histogram_file)
			warn("No coverage histogram found for $(vcf.sample[s])...")
			continue
		end
		d = readdlm(histogram_file, '\t')
		for chr in human_chr
			genome_len += d[findfirst(x -> x == chr, d[:, 1]), 4]
			for r in 1:size(d, 1)
				if d[r, 1] != chr || d[r, 2] >= 5000; continue; end
				histogram[d[r, 2] + 1] += d[r, 3]
			end
		end

		for r in 1:size(vcf, 1)
			if !vcf.star[r, s]; continue; end
			min_depth = ceil(Int, alt_reads / (vcf.alt[r, s] / vcf.total[r, s]))
			amenable_mbs = (genome_len - sum(histogram[1:min_depth-1])) / 1e6
			info("$(vcf.chromosome[r]):$(vcf.position[r]) - amenable $(amenable_mbs) Mb")
			if is_protein_altering(vcf.effect[r])
				protein_altering += 1
				protein_altering_rate += 1 / amenable_mbs
			else
				silent += 1
				silent_rate += 1 / amenable_mbs
			end
		end

		@printf("%s\t%d\t%d\t%.2f\t%.2f\n", vcf.sample[s], protein_altering, silent, protein_altering_rate, silent_rate)
	end
end

function discard_sketchy_silent(vcf_file::IO)
	line = skip_to_header(vcf_file); println(line)
	notes_col = findfirst(split(line, '\t') .== "NOTES")
	effect_col = findfirst(split(line, '\t') .== "EFFECT")
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		if !is_protein_altering(cols[effect_col])
			if contains(cols[notes_col], "INDEL within"); continue; end
			if contains(cols[notes_col], "Mappability < "); continue; end
		end
		println(line)
	end
end

function sort_samples_by_name(vcf_file::IO)
	line = skip_to_header(vcf_file)
	headers = split(line, '\t')
	notes_col = findfirst(headers .== "NOTES")
	samples = headers[notes_col+1:end]
	order = vcat(1:notes_col, notes_col .+ sortperm(samples))
	println(join(headers[order], '\t'))
	for line in eachline(vcf_file)
		cols = split(line, '\t')
		println(join(cols[order], '\t'))
	end
end

subcommands(somatic, germline, above_background, heterozygous_snps,
	keep_samples, discard_samples,
	discard_shallow, discard_if_frequency_above, discard_blacklisted,
	discard_indels, discard_contingent, discard_sketchy_silent,
	mappability, nearby_indels,
	inside, protein_altering, fractions, statistics,
	predict_effect, annotate, build_annotation_database, error_rate,
	group_by_snp_profile, mutation_rate, sort_by_position, sort_samples_by_name)


text:
# Single robot solver using Monte Carlo tree search

using LinearAlgebra

using MCTS
using POMDPs
using Random
import POMDPs.isterminal
import POMDPs.actions

export AbstractSingleRobotProblem, SingleRobotTargetTrackingProblem,
SingleRobotTargetCoverageProblem, MDPState, generate_solver, solve_single_robot,
isterminal

const default_solver_information_samples = 1

const default_horizon = 2

# This approximates the (half) the expected reward as prior tests indicated
# that MCTS performs best in this regime
#
# The exact values were tuned based on the default configuration (number of
# targets, sensor noise) for 4 and 8 robots and consideration of typical
# objective values per robot as obtained from the multi-robot example script
exploration_constant(horizon) = 0.18 * horizon^2 + 0.31

# We map the horizon length to the number of samples according to approximately
# where the average return crossed 97% in a prior test
#
# Note: the last test indicated slow convergence for a horizon length of 2.
const default_num_iterations = Dict(1=>512,
                                    2=>512,
                                    3=>1024,
                                    4=>1024,
                                    5=>2048,
                                    6=>4000)

#
# Define the MDP model for POMDPs.jl
#

# Simple struct for the state
# (defining the state trajectory in place probably ends up being faster than
# not)
#
# We will use the usual 0 indexing for depth (i.e. depth represents a number of
# actions taken)
struct MDPState
  state::State
  depth::Int64
  prev::Union{MDPState, Nothing}
end
# Constructor for initial states
MDPState(state) = MDPState(state, 0, nothing)
# Construct a new state in a trajectory
MDPState(m::MDPState, s::State) = MDPState(s, m.depth + 1, m)

# Abstract implementation of something like a single-robot tracking problem
# These single-robot problems are mdps with MDPStates for states
# and target tracking/coverage States for actions
abstract type AbstractSingleRobotProblem <: MDP{MDPState, State}
end

struct SingleRobotTargetTrackingProblem{F<:AbstractFilter} <: AbstractSingleRobotProblem
  grid::Grid
  sensor::RangingSensor
  horizon::Int64
  target_filters::Vector{F}
  prior_trajectories::Vector{Vector{State}}

  num_information_samples::Int64

  function SingleRobotTargetTrackingProblem(grid::Grid, sensor::RangingSensor,
                                            horizon::Integer,
                                            filters::Vector{T};
                                            prior_trajectories = Trajectory[],
                                            num_information_samples =
                                              default_solver_information_samples
                                           ) where T <: AbstractFilter
    new{T}(grid, sensor, horizon, filters, prior_trajectories,
        num_information_samples)
  end
end

struct SingleRobotTargetCoverageProblem <: AbstractSingleRobotProblem
  grid::Grid
  sensor::CoverageSensor
  horizon::Int64
  target_states::Vector{State}
  prior_trajectories::Vector{Vector{State}}

  function SingleRobotTargetCoverageProblem(grid, sensor, horizon,
                                            target_states;
                                            prior_trajectories = Trajectory[])
    new(grid, sensor, horizon, target_states, prior_trajectories)
  end
end

function generate_solver(depth;
                         n_iterations = default_num_iterations[depth],
                         exploration_constant = exploration_constant(depth),
                         enable_tree_vis = false)

  solver = MCTSSolver(n_iterations = n_iterations,
                      depth = depth,
                      exploration_constant = exploration_constant,
                      enable_tree_vis = enable_tree_vis
                     )
end

# Return the next state for the robot
function solve_single_robot(problem::AbstractSingleRobotProblem,
                            state::State;
                            n_iterations =
                              default_num_iterations[problem.horizon],
                            exploration_constant =
                              exploration_constant(problem.horizon))

  solver = generate_solver(horizon(problem), n_iterations = n_iterations,
                           exploration_constant = exploration_constant)
  policy = solve(solver, problem)

  action, info = action_info(policy, MDPState(state))
  tree = info[:tree]
  trajectory = extract_trajectory(problem, tree, state)

  (
   action=action,
   trajectory=trajectory,
   tree=tree
  )
end

# Pull the tree out of the info object
function extract_trajectory(problem::AbstractSingleRobotProblem,
                            tree::MCTS.MCTSTree,
                            initial_state::State)

  # Recall that our State objects are the actions in the MDP
  states = Vector{State}(undef, horizon(problem))

  mdp_state = MDPState(initial_state)

  # Allow for cases where the tree might not all be stored by sampling the rest
  # of the rollout
  for ii in 1:horizon(problem)
    state = try
      action(MCTS.best_sanode_Q(MCTS.StateNode(tree, mdp_state)))
    catch e
      target_dynamics(problem.grid, mdp_state.state)
    end

    states[ii] = state

    # Propagate the trajectory
    mdp_state = MDPState(mdp_state, state)
  end

  states
end

horizon(x::AbstractSingleRobotProblem) = x.horizon

# I probably dont need to define this, but it's good to have
function isterminal(model::AbstractSingleRobotProblem, state)
  state.depth == model.horizon
end

function actions(model::AbstractSingleRobotProblem, state::MDPState)
  neighbors(model.grid, state.state)
end
function actions(model::AbstractSingleRobotProblem)
  println("providing all actions")
  get_states(model.grid)
end
POMDPs.discount(model::AbstractSingleRobotProblem) = 1.0

# This tree implementation should be fine since everything should already be in
# cache once we reach the end.
# The alternative probably involves messy equality operators
#
# Note that the trajectory will *not* include the initial position
function state_trajectory(state::MDPState)
  ret = Array{State}(undef, state.depth)

  ret[1] = state.state

  curr_state = state
  for ii = 2:state.depth
    curr_state = curr_state.prev
    ret[ii] = curr_state.state
  end

  ret
end

# Information reward for target tracking problems
# * sample reward for all targets and sums
function sample_reward(model::SingleRobotTargetTrackingProblem,
                       trajectory::Trajectory; kwargs...)

  # Compute reward conditional on prior selections (trajectories)
  #
  # The more complex "mapreduce" method is necessary here instead of the simpler
  # "sum" since we may encounter empty arrays.
  mapreduce(+, model.target_filters, init=0.0) do filter
    finite_horizon_information(model.grid, filter, model.sensor,
                               trajectory;
                               prior_observations=model.prior_trajectories,
                               num_samples=model.num_information_samples,
                               kwargs...).reward
  end
end

# Coverage reward for target coverage problems
# * Computes *exact* reward for all targets and sums
function sample_reward(model::SingleRobotTargetCoverageProblem,
                       trajectory::Trajectory; rng = Nothing, kwargs...)

  # Compute reward conditional on prior selections (trajectories)
  #
  # The more complex "mapreduce" method is necessary here instead of the simpler
  # "sum" since we may encounter empty arrays.
  mapreduce(+, model.target_states, init=0.0) do target_state
    finite_horizon_coverage(model.grid, target_state, model.sensor, trajectory;
                            prior_trajectories=model.prior_trajectories,
                            kwargs...).reward
  end
end

# Generative model for the single-robot MDP
#
# We will sample the reward once we reach the end of the trajectory
function POMDPs.gen(model::AbstractSingleRobotProblem,
                    s::MDPState, a::State, rng)
  new_state = MDPState(s, a)
  reward = 0.0

  if isterminal(model, new_state)
    reward = sample_reward(model, state_trajectory(new_state); rng = rng)
  end

  if isnan(reward)
    error("Information reward (", reward, ") is NaN")
  end

  (sp=new_state , r=reward)
end

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
#using DataStructures
collections = pyimport("collections") # @pyimport collections

# import synapseclient.utils as utils
# from synapseclient.annotations import to_synapse_annotations, from_synapse_annotations, to_submission_status_annotations, from_submission_status_annotations, set_privacy
# from synapseclient.exceptions import *
using SynapseClient.annotations

@testset "annotations" begin
    # """Test string annotations"""
    a = Dict("foo"=>"bar", "zoo"=>["zing","zaboo"], "species"=>"Platypus")
    sa = to_synapse_annotations(a)
    # print(sa)
    @test sa["stringAnnotations"]["foo"] == ["bar"]
    @test sa["stringAnnotations"]["zoo"] == ["zing","zaboo"]
    @test sa["stringAnnotations"]["species"] == ["Platypus"]
end

@testset "annotation_name_collision" begin
    # """Test handling of a name collisions between typed user generated and untyped
       # system generated annotations, see SYNPY-203 and PLFM-3248"""

    ## order is important: to repro the erro, the key uri has to come before stringAnnotations
    sa = collections.OrderedDict()
    sa["uri"] = "/entity/syn47396/annotations"
    sa["doubleAnnotations"] = Dict()
    sa["longAnnotations"] = Dict()
    sa["stringAnnotations"] = Dict(
            "tissueType"=> ["Blood"],
            "uri"=> ["/repo/v1/dataset/47396"])
    sa["creationDate"] = "1321168909232"
    sa["id"] = "syn47396"

    a = from_synapse_annotations(sa)
    @test a["tissueType"] == ["Blood"]
##    assert a["uri"] == u"/entity/syn47396/annotations"
end

@testset "more_annotations" begin
    # """Test long, float and data annotations"""
    a = Dict("foo"=>1234,
             "zoo"=>[123.1, 456.2, 789.3],
             "species"=>"Platypus",
             # "birthdays"=>[Datetime(1969,4,28), Datetime(1973,12,8), Datetime(2008,1,3)],
             "birthdays"=>[Dates.DateTime(1969,4,28), Dates.DateTime(1973,12,8), Dates.DateTime(2008,1,3)],
             "test_boolean"=>true,
             "test_mo_booleans"=>[false, true, true, false])
    sa = to_synapse_annotations(a)
    println(sa)
    @test sa["longAnnotations"]["foo"] == [1234]
    @test sa["doubleAnnotations"]["zoo"] == [123.1, 456.2, 789.3]
    @test sa["stringAnnotations"]["species"] == ["Platypus"]
    @test sa["stringAnnotations"]["test_boolean"] == ["true"]
    @test lowercase.(sa["stringAnnotations"]["test_mo_booleans"]) == ["false", "true", "true", "false"]

    ## this part of the test is kinda fragile. It it breaks again, it should be removed
    bdays = [utils.from_unix_epoch_time(t) for t in sa["dateAnnotations"]["birthdays"]]
    @test all(Bool[t in bdays for t in [Dates.DateTime(1969,4,28), Dates.DateTime(1973,12,8), Dates.DateTime(2008,1,3)]]) == true
end
@testset "annotations_unicode" begin
    a = Dict("files"=>["tmp6y5tVr.txt"], "cacheDir"=>"/Users/chris/.synapseCache/python/syn1809087", "foo"=>1266)
    sa = to_synapse_annotations(a)
    @test sa["stringAnnotations"]["cacheDir"] == ["/Users/chris/.synapseCache/python/syn1809087"]
end
@testset "round_trip_annotations" begin
    # """Test that annotations can make the round trip from a simple dictionary
    # to the synapse format and back"""
    a = Dict("foo"=>1234, "zoo"=>[123.1, 456.2, 789.3], "species"=>"Moose", "birthdays"=>[Dates.DateTime(1969,4,28), Dates.DateTime(1973,12,8), Dates.DateTime(2008,1,3), Dates.DateTime(2013,3,15)])
    sa = to_synapse_annotations(a)
    # print(sa)
    a2 = from_synapse_annotations(sa)
    # print(a2)
    a = a2 # TODO: Shouldn't there be a test here?
end
@testset "mixed_annotations" begin
    # """test that to_synapse_annotations will coerce a list of mixed types to strings"""
    a = Dict("foo"=>[1, "a", Dates.DateTime(1969,4,28,11,47)])
    sa = to_synapse_annotations(a)
    # print(sa)
    a2 = from_synapse_annotations(sa)
    # print(a2)
    @test a2["foo"][1] == "1"
    @test a2["foo"][2] == "a"
    @test findfirst("1969",a2["foo"][3]) != nothing
end

@testset "idempotent_annotations" begin
    # """test that to_synapse_annotations won"t mess up a dictionary that"s already
    # in the synapse format"""
    a = Dict("species"=>"Moose", "n"=>42, "birthday"=>Dates.DateTime(1969,4,28))
    sa = to_synapse_annotations(a)
    a2 = Dict()
    merge!(a2,sa)#a2.update(sa)
    sa2 = to_synapse_annotations(a2)
    @test sa == sa2
end
@testset "submission_status_annotations_round_trip" begin
    april_28_1969 = Dates.DateTime(1969,4,28)
    a = Dict("screen_name"=>"Bullwinkle", "species"=>"Moose", "lucky"=>13, "pi"=>float(pi), "birthday"=>april_28_1969)
    sa = to_submission_status_annotations(a)
    println(sa)
    @test Set([kvp["key"]   for kvp in sa["stringAnnos"]]) == Set(["screen_name","species"])
    @test Set([kvp["value"] for kvp in sa["stringAnnos"]]) == Set(["Bullwinkle","Moose"])

    ## test idempotence
    @test to_submission_status_annotations(sa) == sa

    @test Set([kvp["key"] for kvp in sa["longAnnos"]]) == Set(["lucky", "birthday"])
    for kvp in sa["longAnnos"]
        key = kvp["key"]
        value = kvp["value"]
        key=="lucky" && @test value == 13

        key=="birthday" && @test utils.from_unix_epoch_time(value) == april_28_1969
    end

    @test Set([kvp["key"] for kvp in sa["doubleAnnos"]]) == Set(["pi"])
    @test Set([kvp["value"] for kvp in sa["doubleAnnos"]]) == Set([float(pi)])

    set_privacy(sa, key="screen_name", is_private=false)
    @test_pythrows PyKeyError set_privacy(sa, key="this_key_does_not_exist", is_private=false)

    for kvp in sa["stringAnnos"]
        kvp["key"] == "screen_name" && @test kvp["isPrivate"] == false

    end
    a2 = from_submission_status_annotations(sa)
    # TODO: is there a way to convert dates back from longs automatically?
    a2["birthday"] = utils.from_unix_epoch_time(a2["birthday"])
    @test a == a2

    ## test idempotence
    @test from_submission_status_annotations(a) == a
end
@testset "submission_status_double_annos" begin
    ssa = Dict("longAnnos"=>   [Dict("isPrivate"=>false, "value"=>13, "key"=>"lucky")],
           "doubleAnnos"=>[Dict("isPrivate"=>false, "value"=>3, "key"=>"three"), Dict("isPrivate"=>false, "value"=>float(pi), "key"=>"pi")])
    ## test that the double annotation "three":3 is interpretted as a floating
    ## point 3.0 rather than an integer 3
    annotations = from_submission_status_annotations(ssa)
    @test typeof(annotations["three"]) == Float64
    ssa2 = to_submission_status_annotations(annotations)
    @test Set([kvp["key"] for kvp in ssa2["doubleAnnos"]]) == Set(["three", "pi"])
    @test Set([kvp["key"] for kvp in ssa2["longAnnos"]]) == Set(["lucky"])
end


text:
import Base

const Iterable = Union{Set, AbstractSet, AbstractArray}

const NullableAbstractString = Union{AbstractString, Nothing}

Base.map(f, s::Set) = Set(Base.map(f, collect(s)))

Base.:-(a::AbstractSet, b::AbstractSet) = setdiff(a, b)

Base.:-(a::AbstractSet, b) = setdiff(a, Set([b]))


text:
function _try_extract_clickable_url_from_git_url(git_url::AbstractString)
    git_url_stripped = convert(String, strip(git_url))::String
    if startswith(lowercase(git_url_stripped), "https://")
        return git_url_stripped
    end
    return nothing
end

function _git_url_to_formatted_markdown(git_url::AbstractString)
    git_url_stripped = convert(String, strip(git_url))::String
    clickable_url_maybe = _try_extract_clickable_url_from_git_url(git_url_stripped)
    if clickable_url_maybe === nothing
        @warn "Could not convert Git URL into a clickable link" git_url
        return "`$(git_url_stripped)`"
    end
    return clickable_url_maybe
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D
using GLMakie

include("./models2D.jl")
include("./plot2D.jl")

#### calculation
elemsize = 0.01
n = 60

# m = model_circles_in_circle_rand_full(0.1, 1.8, elemsize)
# m = model_circles_in_circle_rand_half(0.1, 1.8, elemsize)
m = model_circles_in_circle_rand_quarter(0.1, 1.8, elemsize)
# m = model_circles_in_circle_cross(0.1, 1.8, 9, elemsize)

# tiles
dx, dy = create_tiles(m, n)
@time t_occ = check_tile_occupation(m, dx, dy, n, blockparts = 1:(m.npar-1))
tile_occ_analysis(t_occ, printit = true)

println("elem size / tile size: ", round(elemsize/((dx+dy)/2), digits = 2), " %")

# classic algo
vfmat = zeros(Float64, m.nelem, m.nelem)
existing_vf!(m, vfmat)
@time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
calculating_vf!(m, vfmat, normit = false)
vfmatp = vfmat_to_parts(m, vfmat, normit = false)

# simplified algo
vfmat_simp = zeros(Float64, m.nelem, m.nelem)
existing_vf!(m, vfmat_simp)
@time blocking_vf_with_tiles_simplified!(m, vfmat_simp, dx, dy, n, t_occ, elem_in_t = 3, skipped_t = 3)
calculating_vf!(m, vfmat_simp, normit = false)
vfmatp_simp = vfmat_to_parts(m, vfmat_simp, normit = false)

#### 2D plot
fig = Figure(resolution = (900, 900))
ax = fig[1, 1] = Axis(fig)
ax.xlabel = "X in m"
ax.ylabel = "Y in m"
ax.aspect = DataAspect()
plot_empty_tiles(fig, ax, dx, dy, n)
plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
plot_model(fig, ax, m, shownvec = true, shownodes = false, showcom = false, showleg = false)
display(fig)

#### compare results of classic algo with simplified algo

control_classic = sum(vfmatp, dims=2)
control = zeros(size(control_classic,1),3)
control[:,1] = control_classic

control_simp = sum(vfmatp_simp, dims=2)
control[:,2] = control_simp

diff = control_classic - control_simp
control[:,3] = diff

println("classic algo / simplified algo / diff")
for i = 1:size(control,1)
    println(control[i,:])
end

sum_classic = sum(control_classic)
sum_simp = sum(control_simp)
println("error rate: ", round((1-sum_simp/sum_classic)*100, digits = 2), " %")

# include("./test/run_algo_vs_simpalgo.jl")


text:
# ClusterManager for HTCondor

export LocalAffinityManager, AffinityMode, COMPACT, BALANCED

@enum AffinityMode COMPACT BALANCED

mutable struct LocalAffinityManager <: ClusterManager
    affinities::Array{Int}

    function LocalAffinityManager(;np=Sys.CPU_THREADS, mode::AffinityMode=BALANCED, affinities::Array{Int}=[])
        @assert(Sys.Sys.KERNEL == :Linux)

        if length(affinities) == 0
            if mode == COMPACT
                affinities = [i%Sys.CPU_THREADS for i in 1:np]
            else
                # mode == BALANCED
                if np > 1
                    affinities = [Int(floor(i)) for i in range(0, stop=Sys.CPU_THREADS - 1e-3, length=np)]
                else
                    affinities = [0]
                end
            end
        end

        return new(affinities)
    end
end

function launch(manager::LocalAffinityManager, params::Dict, launched::Array, c::Condition)
    dir = params[:dir]
    exename = params[:exename]
    exeflags = params[:exeflags]

    for core_id in manager.affinities
        io = open(detach(
            setenv(`taskset -c $core_id $(Base.julia_cmd(exename)) $exeflags $(worker_arg())`, dir=dir)), "r")
        wconfig = WorkerConfig()
        wconfig.process = io
        wconfig.io = io.out
        push!(launched, wconfig)
    end

    notify(c)
end

function manage(manager::LocalAffinityManager, id::Integer, config::WorkerConfig, op::Symbol)
    if op == :interrupt
        kill(get(config.process), 2)
    end
end


text:
export isAscii

function isascii(str::AbstractString)::Bool
    isAsciiReg = r"^[\x00-\x7F]*$"

    return contains(str, isAsciiReg)
end


text:
import MNIST

using Colors
using Images
using ImageView

"""
    mnist_img(row)

    Obtain the row of the MNIST data matrix as image.
"""
function mnist_img(row)
    img = Images.grayim(convert( Images.Image{Gray}, reshape(row, (28,28)) ))
    img["spatialorder"] = ["y", "x"]
    return img
end

"""
    mnist_view(img)

    View an image obtained by mnist_img.
"""
mnist_view(img) = ImageView.view(img) # xy=["y","x"])

"""
    mnist_saveview(filename, imgview)

    Save the image obtained by mnist_view(example)
"""
mnist_saveimg(filename, img) = Images.save(filename, img)

function drawmnist()
    X, y = MNIST.traindata()

    srand(1337)

    # select first three of each class, normalize and transpose
    X_0 = X[:,(y.==0.0)][:,1:3]' ./ 255.0
    X_1 = X[:,(y.==1.0)][:,1:3]' ./ 255.0
    X_2 = X[:,(y.==2.0)][:,1:3]' ./ 255.0
    X_3 = X[:,(y.==3.0)][:,1:3]' ./ 255.0
    X_4 = X[:,(y.==4.0)][:,1:3]' ./ 255.0
    X_5 = X[:,(y.==5.0)][:,1:3]' ./ 255.0
    X_6 = X[:,(y.==6.0)][:,1:3]' ./ 255.0
    X_7 = X[:,(y.==7.0)][:,1:3]' ./ 255.0
    X_8 = X[:,(y.==8.0)][:,1:3]' ./ 255.0
    X_9 = X[:,(y.==9.0)][:,1:3]' ./ 255.0

    # cleanup
    X = nothing
    y = nothing
    gc()

    mkdir("./img")
    for i in 1:3
        mnist_saveimg("./img/0_"*string(i)*".pdf", mnist_img(X_0[i,:]))
        mnist_saveimg("./img/1_"*string(i)*".pdf", mnist_img(X_1[i,:]))
        mnist_saveimg("./img/2_"*string(i)*".pdf", mnist_img(X_2[i,:]))
        mnist_saveimg("./img/3_"*string(i)*".pdf", mnist_img(X_3[i,:]))
        mnist_saveimg("./img/4_"*string(i)*".pdf", mnist_img(X_4[i,:]))
        mnist_saveimg("./img/5_"*string(i)*".pdf", mnist_img(X_5[i,:]))
        mnist_saveimg("./img/6_"*string(i)*".pdf", mnist_img(X_6[i,:]))
        mnist_saveimg("./img/7_"*string(i)*".pdf", mnist_img(X_7[i,:]))
        mnist_saveimg("./img/8_"*string(i)*".pdf", mnist_img(X_8[i,:]))
        mnist_saveimg("./img/9_"*string(i)*".pdf", mnist_img(X_9[i,:]))
    end

end

drawmnist()

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
function calc_dt_for_cleaning_speed(cfl::Real, mesh,
                                    equations::AbstractIdealGlmMhdEquations, dg::DG, cache)
# compute time step for GLM linear advection equation with c_h=1 for the DG discretization on
# Cartesian meshes
  max_scaled_speed_for_c_h = maximum(cache.elements.inverse_jacobian) * ndims(equations)
  # OBS! This depends on the implementation details of the StepsizeCallback and needs to be adapted
  #      as well if that callback changes.
  return cfl * 2 / (nnodes(dg) * max_scaled_speed_for_c_h)
end


text:
export nbc

import ..NBCModel

using mlpack._Internal.io

import mlpack_jll
const nbcLibrary = mlpack_jll.libmlpack_julia_nbc

# Call the C binding of the mlpack nbc binding.
function nbc_mlpackMain()
  success = ccall((:nbc, nbcLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module nbc_internal
  import ..nbcLibrary

import ...NBCModel

# Get the value of a model pointer parameter of type NBCModel.
function IOGetParamNBCModel(paramName::String)::NBCModel
  NBCModel(ccall((:IO_GetParamNBCModelPtr, nbcLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type NBCModel.
function IOSetParamNBCModel(paramName::String, model::NBCModel)
  ccall((:IO_SetParamNBCModelPtr, nbcLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeNBCModel(stream::IO, model::NBCModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeNBCModelPtr, nbcLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeNBCModel(stream::IO)::NBCModel
  buffer = read(stream)
  NBCModel(ccall((:DeserializeNBCModelPtr, nbcLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    nbc(; [incremental_variance, input_model, labels, test, training, verbose])

This program trains the Naive Bayes classifier on the given labeled training
set, or loads a model from the given model file, and then may use that trained
model to classify the points in a given test set.

The training set is specified with the `training` parameter.  Labels may be
either the last row of the training set, or alternately the `labels` parameter
may be specified to pass a separate matrix of labels.

If training is not desired, a pre-existing model may be loaded with the
`input_model` parameter.

The `incremental_variance` parameter can be used to force the training to use an
incremental algorithm for calculating variance.  This is slower, but can help
avoid loss of precision in some cases.

If classifying a test set is desired, the test set may be specified with the
`test` parameter, and the classifications may be saved with the
`predictions`predictions  parameter.  If saving the trained model is desired,
this may be done with the `output_model` output parameter.

Note: the `output` and `output_probs` parameters are deprecated and will be
removed in mlpack 4.0.0.  Use `predictions` and `probabilities` instead.

For example, to train a Naive Bayes classifier on the dataset `data` with labels
`labels` and save the model to `nbc_model`, the following command may be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> _, nbc_model, _, _, _ = nbc(labels=labels, training=data)
```

Then, to use `nbc_model` to predict the classes of the dataset `test_set` and
save the predicted classes to `predictions`, the following command may be used:

```julia
julia> using CSV
julia> test_set = CSV.read("test_set.csv")
julia> predictions, _, _, _, _ = nbc(input_model=nbc_model,
            test=test_set)
```

# Arguments

 - `incremental_variance::Bool`: The variance of each class will be
      calculated incrementally.  Default value `false`.

 - `input_model::NBCModel`: Input Naive Bayes model.
 - `labels::Array{Int, 1}`: A file containing labels for the training
      set.
 - `test::Array{Float64, 2}`: A matrix containing the test set.
 - `training::Array{Float64, 2}`: A matrix containing the training set.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Int, 1}`: The matrix in which the predicted labels for
      the test set will be written (deprecated).
 - `output_model::NBCModel`: File to save trained Naive Bayes model to.
 - `output_probs::Array{Float64, 2}`: The matrix in which the predicted
      probability of labels for the test set will be written (deprecated).
 - `predictions::Array{Int, 1}`: The matrix in which the predicted labels
      for the test set will be written.
 - `probabilities::Array{Float64, 2}`: The matrix in which the predicted
      probability of labels for the test set will be written.

"""
function nbc(;
             incremental_variance::Union{Bool, Missing} = missing,
             input_model::Union{NBCModel, Missing} = missing,
             labels = missing,
             test = missing,
             training = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, nbcLibrary), Nothing, ());

  IORestoreSettings("Parametric Naive Bayes Classifier")

  # Process each input argument before calling mlpackMain().
  if !ismissing(incremental_variance)
    IOSetParam("incremental_variance", convert(Bool, incremental_variance))
  end
  if !ismissing(input_model)
    nbc_internal.IOSetParamNBCModel("input_model", convert(NBCModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("output_model")
  IOSetPassed("output_probs")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  nbc_mlpackMain()

  return IOGetParamURow("output"),
         nbc_internal.IOGetParamNBCModel("output_model"),
         IOGetParamMat("output_probs", points_are_rows),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
export left_env

function LinearAlgebra.dot(ϕ::CuMPS, ψ::CuMPS)
    C = CUDA.ones(eltype(ψ), 1, 1)

    for i ∈ eachindex(ψ)
        M = ψ[i]
        M̃ = conj.(ϕ[i])
        @cutensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ)
    end
    return tr(C)
end

function left_env(ϕ::CuMPS, ψ::CuMPS)
    l = length(ψ)
    S = promote_type(eltype(ψ), eltype(ϕ))

    L = Vector{CuMatrix{S}}(undef, l+1)
    L[1] = CUDA.ones(eltype(ψ), 1, 1)

    for i ∈ 1:l
        M = ψ[i]
        M̃ = conj.(ϕ[i])

        C = L[i]
        @cutensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ)
        L[i+1] = C
    end
    return L
end

function right_env(ϕ::CuMPS, ψ::CuMPS)
    L = length(ψ)
    T = promote_type(eltype(ψ), eltype(ϕ))

    R = Vector{CuMatrix{T}}(undef, L+1)
    R[end] = ones(eltype(ψ), 1, 1)

    for i ∈ L:-1:1
        M = ψ[i]
        M̃ = conj.(ϕ[i])

        D = R[i+1]
        @cutensor D[x, y] := M[x, σ, α] * D[α, β] * M̃[y, σ, β] order = (β, α, σ)
        R[i] = D
    end
    return R
end

function LinearAlgebra.dot(ϕ::CuMPS, O::Vector{<:CuMatrix}, ψ::CuMPS)
    S = promote_type(eltype(ψ), eltype(ϕ), eltype(O[1]))
    C = CUDA.ones(S, 1, 1)

    for i ∈ eachindex(ψ)
        M = ψ[i]
        M̃ = conj.(ϕ[i])
        Mat = O[i]
        @cutensor C[x, y] := M̃[β, σ, x] * Mat[σ, η] * C[β, α] * M[α, η, y] order = (α, η, β, σ)
    end
    return tr(C)
end


text:
module SIWPD
export
    siwpd,
    makesiwpdtree

using
    Wavelets

using
    ..Utils,
    ..DWT

"""
    siwpd(x, wt[, L=maxtransformlevels(x), d=L])

Computes the Shift-Invariant Wavelet Packet Decomposition originally developed
by Cohen, Raz & Malah on the vector `x` using the discrete wavelet filter `wt`
for `L` levels with depth `d`.
"""
function siwpd(x::AbstractVector{T}, wt::OrthoFilter,
               L::S = maxtransformlevels(x), d::S = L) where {T<:AbstractFloat, S<:Integer}
    # Sanity check
    n = length(x)
    @assert 0 ≤ L ≤ maxtransformlevels(n)
    @assert 1 ≤ d ≤ L

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)
    # y = Matrix{T}(undef, (n, gettreelength(1<<(L+1))))
    y = zeros(T, (n, gettreelength(1<<(L+1))))
    @inbounds y[:,1] = x

    # Decomposition
    for i in axes(y,2)
        lvl = getdepth(i, :binary)
        len = nodelength(n, lvl)
        dₗ = (0 ≤ lvl ≤ L-d) ? d : L-lvl
        siwpd_subtree!(y, h, g, i, dₗ, 0, L=lvl, len=len)
    end
    return y
end

#=======================================================================================
Recursive calls for decomposition while d>0.
If shift is 0, ie. ss=0, the main root will not be further decomposed even if d>0.
=======================================================================================#
function siwpd_subtree!(y::AbstractMatrix{T₁}, h::Vector{T₂}, g::Vector{T₂},
                        i::S, d::S, ss::S;
                        L::S = getdepth(i, :binary),
                        len::S = nodelength(size(y,1), L)) where
                       {T₁<:AbstractFloat, T₂<:AbstractFloat, S<:Integer}
    # Sanity check
    n, k = size(y)
    @assert 0 ≤ L ≤ maxtransformlevels(n)
    @assert 0 ≤ ss < n
    @assert 0 ≤ d ≤ getdepth(k,:binary)-L
    @assert iseven(len) || len == 1

    # --- Base case ---
    d > 0 || return y

    # --- One level of decomposition ---
    v_start = ss*len+1                                      # Start index for parent node
    v_end = (ss+1)*len                                      # End index for parent node
    @inbounds v = @view y[v_start:v_end, i]                 # Parent node
    w_start = ss*(len÷2)+1                                  # Start index for child node
    w_end = w_start+(len÷2)-1                               # End index for child node
    @inbounds w₁ = @view y[w_start:w_end, getchildindex(i, :left)]    # Left child node
    @inbounds w₂ = @view y[w_start:w_end, getchildindex(i, :right)]   # Right child node
    sidwt_step!(w₁, w₂, v, h, g, false)                          # 1 step of non-shifted decomposition
    w_start = (ss+1<<L)*(len÷2)+1                           # Start index for shifted child node
    w_end = w_start+(len÷2)-1                               # End index for shifted child node
    @inbounds w₁ = @view y[w_start:w_end, getchildindex(i, :left)]    # Left child node
    @inbounds w₂ = @view y[w_start:w_end, getchildindex(i, :right)]   # Right child node
    sidwt_step!(w₁, w₂, v, h, g, true)                           # 1 step of shifted decomposition

    # --- Recursive decomposition while d>0 ---
    if ss > 0       # Decomposition of current shift
        siwpd_subtree!(y, h, g, getchildindex(i, :left), d-1, ss, L=L+1, len=len÷2)
        siwpd_subtree!(y, h, g, getchildindex(i, :right), d-1, ss, L=L+1, len=len÷2)
    end
    # Decomposition of shifted version
    siwpd_subtree!(y, h, g, getchildindex(i, :left), d-1, ss+1<<L, L=L+1, len=len÷2)
    siwpd_subtree!(y, h, g, getchildindex(i, :right), d-1, ss+1<<L, L=L+1, len=len÷2)

    return y
end

# Single step of decomposition
function sidwt_step!(w₁::AbstractVector{T}, w₂::AbstractVector{T},
                v::AbstractVector{T},
                h::Vector{S}, g::Vector{S},
                s::Bool) where {T<:AbstractFloat, S<:AbstractFloat}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)÷2
    @assert length(h) == length(g)

    # Setup
    n = length(v)           # Parent length
    n₁ = length(w₁)         # Child length
    filtlen = length(h)     # Filter length

    # One step of discrete transform
    for i in 1:n₁
        k₁ = mod1(2*i-1-s, n)   # Start index for low pass filtering
        k₂ = 2*i-s              # Start index for high pass filtering
        @inbounds w₁[i] = g[end] * v[k₁]
        @inbounds w₂[i] = h[1] * v[k₂]
        for j in 2:filtlen
            k₁ = k₁+1 |> k₁ -> k₁>n ? mod1(k₁,n) : k₁
            k₂ = k₂-1 |> k₂ -> k₂≤0 ? mod1(k₂,n) : k₂
            @inbounds w₁[i] += g[end-j+1] * v[k₁]
            @inbounds w₂[i] += h[j] * v[k₂]
        end
    end
    return w₁, w₂
end

"""
    makesiwpdtree(n, L, d)

Returns the multi-level, multi-depth binary tree corresponding to the Shift-
Invariant Wavelet Packet Decomposition.
"""
function makesiwpdtree(n::Integer, L::Integer, d::Integer)
    @assert 0 ≤ L ≤ maxtransformlevels(n)
    @assert 1 ≤ d ≤ L

    tree = Vector{BitVector}(undef, 2^(L+1)-1)
    for i in eachindex(tree)
        level = floor(Int, log2(i))
        len = nodelength(n, level)
        nshift = n ÷ len        # number of possible shifts for subspace Ω(i,j)
        node = falses(nshift)
        k = ceil(Int, nshift / 1<<d)  # "skips" between shifts
        node[1:k:end] .= true         # mark nodes present in tree
        tree[i] = node
    end
    return tree
end

end # end module


text:
module GameClient # LOLTools

using HTTP, JSON2

function http_action(verb, endpoint::HTTP.URI, path::String, headers::Vector{Pair{String,String}}=Vector{Pair{String,String}}(), query::Union{Nothing,Dict{String,String}}=nothing, body=nothing)
    url = string(merge(endpoint, path=path))
    sc = HTTP.Servers.SSLConfig(false)
    args = (query=query, sslconfig=sc)
    if body === nothing
        verb(url, headers; args...)
    else
        verb(url, headers; body=body, args...)
    end
end

function lol_replay_server(; host="127.0.0.1", port=2999)::HTTP.URI
    merge(HTTP.URI("https://"), host=host, port=port, path="/")
end

function call_api(endpoint::HTTP.URI, action::Function, path::String)
    resp = action(HTTP.get, endpoint, path)
    data = JSON2.read(IOBuffer(resp.body))
end

function call_api(endpoint::HTTP.URI, action::Function, path::String, data)
    headers = ["Content-Type" => "application/json"]
    body = JSON2.write(data)
    resp = action(HTTP.post, endpoint, path, headers, nothing, body)
    data = JSON2.read(IOBuffer(resp.body))
end

"""
         replay_game( ;
                     endpoint::HTTP.URI = lol_replay_server(),
                     action::Function = http_action)
"""
function replay_game( ;
                     endpoint::HTTP.URI = lol_replay_server(),
                     action::Function = http_action)
    data = call_api(endpoint, action, "/replay/game")
end

"""
         replay_particles( ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
"""
function replay_particles( ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
    data = call_api(endpoint, action, "/replay/particles")
end

"""
         replay_particles(particles ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
"""
function replay_particles(particles ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
    data = call_api(endpoint, action, "/replay/particles", particles)
end

"""
         replay_playback( ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
"""
function replay_playback( ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
    data = call_api(endpoint, action, "/replay/playback")
end

"""
         replay_playback(playback ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
"""
function replay_playback(playback ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
    data = call_api(endpoint, action, "/replay/playback", playback)
end

"""
         replay_render( ;
                       endpoint::HTTP.URI = lol_replay_server(),
                       action::Function = http_action)
"""
function replay_render( ;
                       endpoint::HTTP.URI = lol_replay_server(),
                       action::Function = http_action)
    data = call_api(endpoint, action, "/replay/render")
end

"""
         replay_render(render ;
                       endpoint::HTTP.URI = lol_replay_server(),
                       action::Function = http_action)
"""
function replay_render(render ;
                       endpoint::HTTP.URI = lol_replay_server(),
                       action::Function = http_action)
    data = call_api(endpoint, action, "/replay/render", render)
end

"""
         replay_recording( ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
"""
function replay_recording( ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
    data = call_api(endpoint, action, "/replay/recording")
end

"""
         replay_recording(recording ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
"""
function replay_recording(recording ;
                          endpoint::HTTP.URI = lol_replay_server(),
                          action::Function = http_action)
    data = call_api(endpoint, action, "/replay/recording", recording)
end

"""
         replay_sequence( ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
"""
function replay_sequence( ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
    data = call_api(endpoint, action, "/replay/sequence")
end

"""
         replay_sequence(sequence ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
"""
function replay_sequence(sequence ;
                         endpoint::HTTP.URI = lol_replay_server(),
                         action::Function = http_action)
    data = call_api(endpoint, action, "/replay/sequence", sequence)
end

"""
         swagger_v2_swagger( ;
                            endpoint::HTTP.URI = lol_replay_server(),
                            action::Function = http_action)
"""
function swagger_v2_swagger( ;
                            endpoint::HTTP.URI = lol_replay_server(),
                            action::Function = http_action)
    data = call_api(endpoint, action, "/swagger/v2/swagger.json")
end

"""
         swagger_v3_openapi( ;
                            endpoint::HTTP.URI = lol_replay_server(),
                            action::Function = http_action)
"""
function swagger_v3_openapi( ;
                            endpoint::HTTP.URI = lol_replay_server(),
                            action::Function = http_action)
    data = call_api(endpoint, action, "/swagger/v3/openapi.json")
end

end # module LOLTools.GameClient


text:
function evaluate(a::STaylor1{N,T}, dx::T) where {N, T<:Number}
    @inbounds suma = a[N-1]
    @inbounds for k in (N-1):-1:0
        suma = suma*dx + a[k]
    end
    suma
end
#=
function evaluate(a::STaylor1{N,T}, dx::S) where {N, T<:Number, S<:Number}
    suma = a[N-1]*one(dx)
    @inbounds for k in (N-1):-1:0
        suma = suma*dx + a[k]
    end
    suma
end
=#
evaluate(a::STaylor1{N,T}) where {N, T<:Number} = a[0]

evaluate(x::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}, δt::S) where
        {N, T<:Number, S<:Number} = evaluate.(x, δt)
evaluate(a::Union{Array{STaylor1{N,T}}, SubArray{STaylor1{N,T}}}) where
        {N, T<:Number} = evaluate.(a, zero(T))

(p::STaylor1)(x) = evaluate(p, x)
(p::STaylor1)() = evaluate(p)

(p::Array{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x)
(p::SubArray{STaylor1{N,T}})(x) where {N,T<:Number} = evaluate.(p, x)
(p::Array{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p)
(p::SubArray{STaylor1{N,T}})() where {N,T<:Number} = evaluate.(p)


text:
function title(s::AbstractString)
  join([ucfirst(w) for w in split(replace(s, "_", " "))], " ")::typeof(s)
end


text:
struct Interval{T}
    a::T
    b::T
end
Interval(a::T,b::U) where {T,U} = Interval(promote(a,b)...)

import Base: ∪, ∩, \, ⊆, ∈, isempty, show, isless

isless(i₁::Interval,i₂::Interval) = i₁.a < i₂.a
adjacent(i₁::Interval,i₂::Interval) = i₁.b == i₂.a
isempty(i::Interval) = i.a >= i.b

⊂(i₁::Interval,i₂::Interval) = i₂.a < i₁.a && i₁.b < i₂.b
⊆(i₁::Interval,i₂::Interval) = i₂.a <= i₁.a && i₁.b <= i₂.b

# This imposes right-continuity, as de Boor uses.
∈(e::Number, i::Interval) = i.a < e && e ≤ i.b

function ∩(i₁::Interval,i₂::Interval)
    a = max(i₁.a,i₂.a)
    b = min(i₁.b,i₂.b)
    Interval(a,b)
end

function ∪(i₁::Interval,i₂::Interval)
    a = min(i₁.a,i₂.a)
    b = max(i₁.b,i₂.b)
    if !isempty(i₁ ∩ i₂) || adjacent(i₁,i₂) || adjacent(i₂,i₁)
        Interval(a,b)
    else
        Interval(1,0)
    end
end

function ∪(i₁::Vector{Interval{T}}, i₂::Interval{T}) where T
    sort!(vcat([i ∪ i₂ for i ∈ i₁]...))
    # Should merge possibly adjacent intervals
end

function \(i₁::Interval,i₂::Interval)
    ii = i₁ ∩ i₂
    ia = Interval(i₁.a,i₂.a)
    ib = Interval(i₂.b,i₁.b)
    if isempty(ii)
        [i₁],[i₂]
    elseif i₂ ⊆ i₁
        if isempty(ia) && isempty(ib)
            Interval[],Interval[]
        elseif isempty(ia)
            [ib],Interval[]
        elseif isempty(ib)
            [ia],Interval[]
        else
            [ia,ib],Interval[]
        end
    elseif i₁ ⊂ i₂
        Interval[],[Interval(i₂.a,i₁.a),Interval(i₁.b,i₂.b)]
    else # i₁ and i₂ partially overlap
        if i₁.a <= i₂.a
            isempty(ia) ? Interval[] : [ia],[Interval(i₁.b,i₂.b)]
        else
            isempty(ib) ? Interval[] : [ib],[Interval(i₂.a,i₁.a)]
        end
    end
end

function \(i₁::Vector{Interval}, i₂::Interval)
    is = Interval[]
    for i ∈ i₁
        i′ = i \ i₂
        !isempty(i′) && append!(is, i′)
    end
    is
end

function collapse!(is::Vector{Interval{T}}) where T
    isempty(is) && return is
    sort!(is)
    nis = Interval{T}[is[1]]
    i = 1
    for j = 2:length(is)
        ij = nis[i] ∪ is[j]
        if !isempty(ij)
            nis[i] = ij
        else
            push!(nis, is[j])
            i += 1
        end
    end
    is[1:length(nis)] = nis
    resize!(is, length(nis))
    is
end

function show(io::IO, i::Interval)
    if isempty(i)
        write(io, "Empty Interval")
    else
        write(io, "Interval[$(i.a),$(i.b)]")
    end
end

export Interval, ⊂


text:
include("./1-indrnn.jl")
include("./2-indlstm.jl")
include("./3-rnn.jl")

export RNN
export RNNs
export IndRNN
export IndRNNs
export IndLSTM
export IndLSTMs

export resethidden
export PadSeqPackBatch
export PackedSeqPredict
export PackedSeqForward

global RNNLIST = [IndRNN,IndRNNs,
                  IndLSTM,IndLSTMs,
                  RNN,RNNs];
# global RNNLIST = [RNN, rin, lstm, IndRNN, IndLSTM, RNNs, RIN, LSTM, IndLSTMs];

"""
    PadSeqPackBatch(inputs::Vector; epsilon::Real=0.0) -> output
+ `inputs` <: AbstractArray{Real,2}
+ `output` <: AbstractArray{Real,3}
pad epsilon to align raw input features probably with different length
# Examples
    julia> PadSeqPackBatch([ones(2,1), 2ones(2,2), 3ones(2,3)])
    2×3×3 Array{Float64,3}:
    [:, :, 1] =
     1.0  0.0  0.0
     1.0  0.0  0.0

    [:, :, 2] =
     2.0  2.0  0.0
     2.0  2.0  0.0

    [:, :, 3] =
     3.0  3.0  3.0
     3.0  3.0  3.0
"""
function PadSeqPackBatch(inputs::Vector; epsilon::Real=0.0)
    # all Array of inputs shall have the same size in dim-1
    batchSize = length(inputs)
    lengths   = [size(inputs[i], 2) for i in 1:batchSize]
    featDims  = size(inputs[1], 1)
    maxSteps  = maximum(lengths)
    RNNBatch  = zeros(eltype(inputs[1]), featDims, maxSteps, batchSize)
    fill!(RNNBatch, epsilon)

    for i = 1:batchSize
        Tᵢ = lengths[i]
        RNNBatch[:,1:Tᵢ,i] .= inputs[i]
    end
    return RNNBatch
end

function PackedSeqForward(chain::Block, x::Variable{S}) where S
    T = size(x, 2)
    y = Vector{Variable{S}}(undef, T)
    resethidden(chain)
    for t = 1:T
        y[t] = forward(chain, x[:,t,:])
    end

    timeSteps = T
    featsDims = size(y[1], 1)
    batchSize = size(y[1], 2)
    RNNBatch  = Variable{S}(S(undef, featsDims, timeSteps, batchSize), x.backprop)

    for t = 1:timeSteps
        RNNBatch.value[:,t,:] .= y[t].value
    end

    if RNNBatch.backprop
        function PackSeqSlicesBackward()
            for t = 1:timeSteps
                if need2computeδ!(y[t])
                    y[t].delta .+= RNNBatch.delta[:,t,:]
                end
            end
            ifNotKeepδThenFreeδ!(RNNBatch)
        end
        push!(graph.backward, PackSeqSlicesBackward)
    end
    return RNNBatch
end

function PackedSeqPredict(chain::Block, x::AbstractArray{S}) where S
    T = size(x,2)
    y = Vector{AbstractArray{S}}(undef,T)
    resethidden(chain)
    for t = 1:T
        y[t] = predict(chain, x[:,t,:])
    end

    timeSteps = T
    featsDims = size(y[1], 1)
    batchSize = size(y[1], 2)
    RNNBatch  = typeof(x)(undef, featsDims, timeSteps, batchSize)
    for t = 1:timeSteps
        RNNBatch[:,t,:] .= y[t]
    end
    return RNNBatch
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
abstract type StandardBasisPolynomial{T,X} <: AbstractPolynomial{T,X} end

function showterm(io::IO, ::Type{<:StandardBasisPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T}

    if _iszero(pj) return false end

    pj = printsign(io, pj, first, mimetype)

    if !(_isone(pj) && !(showone(T) || j == 0))
        printcoefficient(io, pj, j, mimetype)
    end

    printproductsign(io, pj, j, mimetype)
    printexponent(io, var, j, mimetype)
    return true
end

"""
    evalpoly(x, p::StandardBasisPolynomial)
    p(x)

Evaluate the polynomial using [Horner's Method](https://en.wikipedia.org/wiki/Horner%27s_method), also known as synthetic division, as implemented in `evalpoly` of base `Julia`.

# Examples
```jldoctest
julia> using Polynomials

julia> p = Polynomial([1, 0, 3])
Polynomial(1 + 3*x^2)

julia> p(0)
1

julia> p.(0:3)
4-element Vector{Int64}:
  1
  4
 13
 28
```
"""
evalpoly(x, p::StandardBasisPolynomial) = EvalPoly.evalpoly(x, p.coeffs) # allows  broadcast  issue #209
constantterm(p::StandardBasisPolynomial) = p[0]

domain(::Type{<:StandardBasisPolynomial}) = Interval(-Inf, Inf)
mapdomain(::Type{<:StandardBasisPolynomial}, x::AbstractArray) = x

function Base.convert(P::Type{<:StandardBasisPolynomial}, q::StandardBasisPolynomial)
    if isa(q, P)
        return q
    else
        T = _eltype(P,q)
        X = indeterminate(P,q)
        return ⟒(P){T,X}([q[i] for i in 0:degree(q)])
    end
end

function Base.one(::Type{P}) where {P<:StandardBasisPolynomial}
    T,X = eltype(P), indeterminate(P)
    ⟒(P){T,X}(ones(T,1))
end
function variable(::Type{P}) where {P<:StandardBasisPolynomial}
    T,X = eltype(P), indeterminate(P)
    ⟒(P){T,X}([zero(T),one(T)])
end

## multiplication algorithms for computing p * q.
## default multiplication between same type.
## subtypes might relax to match T,S to avoid one conversion
function Base.:*(p::P, q::P) where {T,X, P<:StandardBasisPolynomial{T,X}}
    cs = ⊗(P, coeffs(p), coeffs(q))
    P(cs)
end

## put here, not with type defintion, in case reuse is possible
function ⊗(P::Type{<:StandardBasisPolynomial}, p::Vector{T}, q::Vector{S}) where {T,S}
    R = promote_type(T,S)
    fastconv(convert(Vector{R}, p), convert(Vector{R},q))
end

## Static size of product makes generated functions  a good choice
## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl
## convolution of two tuples
@generated function ⊗(::Type{<:StandardBasisPolynomial}, p1::NTuple{N,T}, p2::NTuple{M,S}) where {T,N,S,M}
    P = M + N - 1
    exprs = Any[nothing for i = 1 : P]
    for i in 1 : N
        for j in 1 : M
            k = i + j - 1
            if exprs[k] === nothing
                exprs[k] = :(p1[$i] * p2[$j])
            else
                exprs[k] = :(muladd(p1[$i], p2[$j], $(exprs[k])))
            end
        end
    end

    return quote
        Base.@_inline_meta
        tuple($(exprs...))
    end

end

function ⊗(P::Type{<:StandardBasisPolynomial}, p::Dict{Int,T}, q::Dict{Int,S}) where {T,S}
    R = promote_type(T,S)
    c = Dict{Int,R}()
    for (i,pᵢ) ∈ pairs(p)
        for (j,qⱼ) ∈ pairs(q)
            cᵢⱼ = get(c, i+j, zero(R))
            @inbounds c[i+j] = muladd(pᵢ, qⱼ, cᵢⱼ)
        end
    end
    c
end

## ---
function fromroots(P::Type{<:StandardBasisPolynomial}, r::AbstractVector{T}; var::SymbolLike = :x) where {T <: Number}
    n = length(r)
    c = zeros(T, n + 1)
    c[1] = one(T)
    for j in 1:n, i in j:-1:1
        c[(i + 1)] = c[(i + 1)] - r[j] * c[i]
    end
    #return P(c, var)
    if sort(r[imag(r).>0], lt = (x,y) -> real(x)==real(y) ? imag(x)<imag(y) : real(x)<real(y)) == sort(conj(r[imag(r).<0]), lt = (x,y) -> real(x)==real(y) ? imag(x)<imag(y) : real(x)<real(y))
        c = real(c)             # if complex poles come in conjugate pairs, the coeffs are real
    end
    return P(reverse(c), var)
end

function derivative(p::P, order::Integer = 1) where {T, X, P <: StandardBasisPolynomial{T, X}}
    order < 0 && throw(ArgumentError("Order of derivative must be non-negative"))

    # we avoid usage like Base.promote_op(*, T, Int) here, say, as
    # Base.promote_op(*, Rational, Int) is Any, not Rational in analogy to
    # Base.promote_op(*, Complex, Int)
    R = eltype(one(T)*1)
    Q = ⟒(P){R,X}

    order == 0 && return p
    hasnan(p) && return Q(R[NaN])
    order > length(p) && return zero(Q)
    d = degree(p)
    d <= 0 && return zero(Q)
    n = d + 1
    a2 = Vector{R}(undef, n - order)
    @inbounds for i in order:n - 1
        a2[i - order + 1] = reduce(*, (i - order + 1):i, init = p[i])
    end
    Q(a2)
end

function integrate(p::P) where {T, X, P <: StandardBasisPolynomial{T, X}}
    R = eltype(one(T)/1)
    Q = ⟒(P){R,X}

    hasnan(p) && return Q([NaN])
    iszero(p) && return zero(Q)

    n = length(p)
    as = Vector{R}(undef, n + 1)
    as[1] = zero(R)
    for (i, pᵢ) ∈ pairs(p)
        i′ = i + 1
        @inbounds as[i′+1] = pᵢ/i′
    end
    return Q(as)
end

function Base.divrem(num::P, den::Q) where {T, P <: StandardBasisPolynomial{T}, S, Q <: StandardBasisPolynomial{S}}

    assert_same_variable(num, den)
    X = indeterminate(num)

    n = degree(num)
    m = degree(den)

    m == -1 && throw(DivideError())
    if m == 0 && den[0] ≈ 0 throw(DivideError()) end

    R = eltype(one(T)/one(S))

    deg = n - m + 1

    if deg ≤ 0
        return zero(P), num
    end

    q_coeff = zeros(R, deg)
    r_coeff = R[ num[i-1] for i in 1:n+1 ]

    @inbounds for i in n:-1:m
        q = r_coeff[i + 1] / den[m]
        q_coeff[i - m + 1] = q
        @inbounds for j in 0:m
            elem = den[j] * q
            r_coeff[i - m + j + 1] -= elem
        end
    end
    resize!(r_coeff, min(length(r_coeff), m))

    return _convert(num, q_coeff), _convert(num, r_coeff)

end

"""
    gcd(p1::StandardBasisPolynomial, p2::StandardBasisPolynomial; method=:eculidean, kwargs...)

Find the greatest common divisor.

By default, uses the Euclidean division algorithm (`method=:euclidean`), which is susceptible to floating point issues.

Passing `method=:noda_sasaki` uses scaling to circumvent some of these.

Passing `method=:numerical` will call the internal method `NGCD.ngcd` for the numerical gcd. See the help page of [`Polynomials.NGCD.ngcd`](@ref) for details.
"""
function Base.gcd(p1::P, p2::Q, args...;
                  method=:euclidean,
                  kwargs...
                  ) where {T, P <: StandardBasisPolynomial{T}, Q <: StandardBasisPolynomial{T}}

    gcd(Val(method), p1, p2, args...; kwargs...)
end

function Base.gcd(::Val{:euclidean},
                  r₀::StandardBasisPolynomial{T}, r₁;
                  atol=zero(real(T)),
                  rtol=Base.rtoldefault(real(T)),
                  kwargs...) where {T}

    iter = 1
    itermax = length(r₁)

    while !iszero(r₁) && iter ≤ itermax
        _, rtemp = divrem(r₀, r₁)
        r₀ = r₁
        r₁ = truncate(rtemp; atol=atol, rtol=rtol)
        iter += 1
    end
    return r₀
end

Base.gcd(::Val{:noda_sasaki}, p, q; kwargs...) = gcd_noda_sasaki(p,q; kwargs...)

"""
     gcd_noda_sasaki(p,q; atol,  rtol)

Greatest common divisor of two polynomials.
Compute the greatest common divisor `d` of two polynomials `a` and `b` using
the Euclidian Algorithm with scaling as of [1].

References:

[1] M.-T. Noda and T. Sasaki. Approximate GCD and its application to ill-conditioned  algebraic equations. J. Comput. Appl. Math., 38:335–351, 1991.

Author: Andreas Varga

Note: requires Julia `v1.2` or greater.
"""
function  gcd_noda_sasaki(p::StandardBasisPolynomial{T}, q::StandardBasisPolynomial{S};
                          atol::Real=zero(real(promote_type(T,S))),
                          rtol::Real=Base.rtoldefault(real(promote_type(T,S)))
                          ) where {T,S}
    ⟒(typeof(p)) == ⟒(typeof(q)) ||  return gcd_noda_sasaki(promote(p,q);  atol=atol, rtol=rtol)
    ## check symbol
    a, b = coeffs(p), coeffs(q)
    as =  _gcd_noda_sasaki(a,b, atol=atol,  rtol=rtol)

    _convert(p, as)
end

function _gcd_noda_sasaki(a::Vector{T}, b::Vector{S};
              atol::Real=zero(real(promote_type(T,S))),
              rtol::Real=Base.rtoldefault(real(promote_type(T,S)))
              ) where {T,S}

    R = eltype(one(T)/one(S))

    na1 = findlast(!iszero,a) # degree(a) + 1
    na1 === nothing && return(ones(R, 1))

    nb1 = findlast(!iszero,b) # degree(b) + 1
    nb1 === nothing && return(ones(R, 1))

    a1 = R[a[i] for i in 1:na1]
    b1 = R[b[i] for i in 1:nb1]
    a1 ./= norm(a1)
    b1 ./= norm(b1)

    tol = atol + rtol

    # determine the degree of GCD as the nullity of the Sylvester matrix
    # this computation can be replaced by simply setting nd = 1, in which case the Sylvester matrix is not formed

    nd = na1 + nb1 - 2 - rank([NGCD.convmtx(a1,nb1-1) NGCD.convmtx(b1,na1-1)], atol = tol) # julia 1.1
    nd == 0 && (return [one(R)])

    sc = one(R)
    while na1 > nd
         na1 < nb1 && ((a1, b1, na1, nb1) = (b1, a1, nb1, na1))
         @inbounds for i in na1:-1:nb1
            s = -a1[i] / b1[nb1]
            sc = max(sc, abs(s))
            @inbounds for j in 1:nb1-1
                a1[i-nb1+j] += b1[j] * s
            end
        end
        a1 ./= sc
        na1 = findlast(t -> (abs(t) > tol),a1[1:nb1-1])
        na1 === nothing && (na1 = 0)
        resize!(a1, na1)
    end

    return nb1 == 1 ? [one(R)] : b1

end

Base.gcd(::Val{:numerical}, p, q, args...; kwargs...) = ngcd(p,q, args...; kwargs...).u

## --------------------------------------------------

function companion(p::P) where {T, P <: StandardBasisPolynomial{T}}
    d = length(p) - 1
    d < 1 && throw(ArgumentError("Series must have degree greater than 1"))
    d == 1 && return diagm(0 => [-p[0] / p[1]])

    R = eltype(one(T)/one(T))

    comp = diagm(-1 => ones(R, d - 1))
    ani = 1 / p[end]
    for j in  0:(degree(p)-1)
        comp[1,(d-j)] = -p[j] * ani # along top row has smaller residual than down column
    end
    return comp
end

function  roots(p::P; kwargs...)  where  {T, P <: StandardBasisPolynomial{T}}

    R = eltype(one(T)/one(T))
    d = degree(p)
    if d < 1
        return R[]
    end
    d == 1 && return R[-p[0] / p[1]]

    as = [p[i] for i in 0:d]
    K  = findlast(!iszero, as)
    if K == nothing
        return R[]
    end
    k =  findfirst(!iszero, as)

    k  == K && return zeros(R, k-1)

    # find eigenvalues of the  companion matrix
    comp  = companion(⟒(P)(as[k:K], indeterminate(p)))
    L = eigvals(comp; kwargs...)
    append!(L, zeros(eltype(L), k-1))

    L
end

function vander(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, n::Integer) where {T <: Number}
    A = Matrix{T}(undef, length(x), n + 1)
    A[:, 1] .= one(T)
    @inbounds for i in 1:n
        A[:, i + 1] = A[:, i] .* x
    end
    return A
end

## as noted at https://github.com/jishnub/PolyFit.jl, using method from SpecialMatrices is faster
## https://github.com/JuliaMatrices/SpecialMatrices.jl/blob/master/src/vandermonde.jl
## This is Algorithm 2 of https://www.maths.manchester.ac.uk/~higham/narep/narep108.pdf
## Solve V(αs)⋅x = y where V is (1+n) × (1+n) Vandermonde matrix (Vᵀ in the paper)
function solve_vander!(ys, αs)
    n = length(ys) - 1
    for k in 0:n-1
        for j in n:-1:k+1
            ys[1+j] = (ys[1+j] - ys[1+j-1])/(αs[1+j] - αs[1+j-k-1])
        end
    end
    for k in n-1:-1:0
        for j in k:n-1
            ys[1+j] = ys[1+j] - αs[1+k] * ys[1 + j + 1]
        end
    end

    nothing
end

# intercept one (typical) case for a faster variant
function fit(P::Type{<:StandardBasisPolynomial},
             x::AbstractVector{T},
             y::AbstractVector{T},
             deg::Integer = length(x) - 1;
             weights = nothing,
             var = :x,) where {T}

    if deg == length(x) -1 && weights == nothing
        _polynomial_fit(P, x, y; var=var)
    else
        _fit(P, x, y, deg; weights=weights, var=var)
    end
end

function _polynomial_fit(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, y; var=:x) where {T}
    R = float(T)
    coeffs = Vector{R}(undef, length(x))
    copyto!(coeffs, y)
    solve_vander!(coeffs, x)
    P(R.(coeffs), var)
end

## --------------------------------------------------
## More accurate fitting for higher degree polynomials

"""
    polyfitA(x, y, n=length(x)-1; var=:x)
    fit(ArnoldiFit, x, y, n=length(x)-1; var=:x)

Fit a degree ``n`` or less polynomial through the points ``(x_i, y_i)`` using Arnoldi orthogonalization of the Vandermonde matrix.

The use of a Vandermonde matrix to fit a polynomial to data is exponentially ill-conditioned for larger values of ``n``. The Arnoldi orthogonalization fixes this problem.

# Returns

Returns an instance of `ArnoldiFit`. This object can be used to evaluate the polynomial. To manipulate the polynomial, the object can be `convert`ed to other polynomial types, though there may be some loss in accuracy when doing polynomial evaluations afterwards for higher-degree polynomials.

# Citations:

The two main functions are translations from example code in:

*VANDERMONDE WITH ARNOLDI*;
PABLO D. BRUBECK, YUJI NAKATSUKASA, AND LLOYD N. TREFETHEN;
[arXiv:1911.09988](https://people.maths.ox.ac.uk/trefethen/vander_revised.pdf)

# Examples:

```
f(x) = 1/(1 + 25x^2)
N = 80; xs = [cos(j*pi/N) for j in N:-1:0];
p = fit(Polynomial, xs, f.(xs));
q = fit(ArnoldiFit, xs, f.(xs));
maximum(abs, p(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 3.304586010148457e16
maximum(abs, q(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 1.1939520722092922e-7

N = 250; xs = [cos(j*pi/N) for j in N:-1:0];
p = fit(Polynomial, xs, f.(xs));
q = fit(ArnoldiFit, xs, f.(xs));
maximum(abs, p(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 3.55318186254542e92
maximum(abs, q(x) - f(x) for x ∈ range(-1,stop=1,length=500)) # 8.881784197001252e-16

p = fit(Polynomial, xs, f.(xs), 10); # least-squares fit
q = fit(ArnoldiFit, xs, f.(xs), 10);
maximum(abs, q(x) - p(x) for x ∈ range(-1,stop=1,length=500)) # 4.6775083806238626e-14
Polynomials.norm(q-p, Inf) # 2.2168933355715126e-12 # promotes `q` to `Polynomial`
```

"""
function polyfitA(x, y, n=length(x)-1; var=:x)
    m = length(x)
    T = eltype(y)
    Q = ones(T, m, n+1)
    #H = UpperHessenberg(zeros(T, n+1, n))
    H = zeros(T, n+1, n)

    q = zeros(T, m)

    @inbounds for k = 1:n
        q .= x .* Q[:,k]
        for j in 1:k
            λ = dot(Q[:,j], q)/m
            H[j,k] = λ
            q .-= λ * Q[:,j]
        end
        H[k+1,k] = norm(q)/sqrt(m)
        Q[:,k+1] .= q/H[k+1,k]
    end
    d = Q \ y
    ArnoldiFit{eltype(d),typeof(H),Symbol(var)}(d, H)
end

function polyvalA(d, H::AbstractMatrix{S}, s::T) where {T, S}
    R = promote_type(T,S)
    n = length(d) - 1
    W = ones(R, n+1)
    @inbounds for k in 1:n
        w = s .* W[k]
        for j in 1:k
            w -= H[j,k] * W[j]
        end
        W[k+1] = w/H[k+1,k]
    end
    sum(W[i]*d[i] for i in eachindex(d))
end

# Polynomial Interface
"""
    ArnoldiFit

A polynomial type produced through fitting a degree ``n`` or less polynomial to data ``(x_1,y_1),…,(x_N, y_N), N ≥ n+1``, This uses Arnoldi orthogonalization to avoid the exponentially ill-conditioned Vandermonde polynomial. See [`Polynomials.polyfitA`](@ref) for details.
"""
struct ArnoldiFit{T, M<:AbstractArray{T,2}, X}  <: AbstractPolynomial{T,X}
    coeffs::Vector{T}
    H::M
end
export ArnoldiFit
@register ArnoldiFit
domain(::Type{<:ArnoldiFit}) = Interval(-Inf, Inf)

Base.show(io::IO, mimetype::MIME"text/plain", p::ArnoldiFit) = print(io, "ArnoldiFit of degree $(length(p.coeffs)-1)")

evalpoly(x, p::ArnoldiFit) = polyvalA(p.coeffs, p.H, x)

fit(::Type{ArnoldiFit}, x::AbstractVector{T}, y::AbstractVector{T}, deg::Int=length(x)-1;  var=:x, kwargs...) where{T} = polyfitA(x, y, deg; var=var)

Base.convert(::Type{P}, p::ArnoldiFit) where {P <: AbstractPolynomial} = p(variable(P,indeterminate(p)))

## --------------------------------------------------
"""
    compensated_horner(p::P, x)
    compensated_horner(ps, x)

Evaluate `p(x)` using a compensation scheme of S. Graillat, Ph. Langlois, N. Louve [Compensated Horner Scheme](https://cadxfem.org/cao/Compensation-horner.pdf). Either a `Polynomial` `p` or its coefficients may be passed in.

The Horner scheme has relative error given by

`|(p(x) - p̂(x))/p(x)| ≤ α(n) ⋅ u ⋅ cond(p, x)`, where `u` is the precision (`2⁻⁵³ = eps()/2`)).

The compensated Horner scheme has relative error bounded by

`|(p(x) - p̂(x))/p(x)| ≤  u +  β(n) ⋅ u² ⋅ cond(p, x)`.

As noted, this reflects the accuracy of a backward stable computation performed in doubled working precision `u²`. (E.g., polynomial evaluation of a `Polynomial{Float64}` object through `compensated_horner` is as accurate as evaluation of a `Polynomial{Double64}` object (using the `DoubleFloat` package), but significantly faster.

Pointed out in https://discourse.julialang.org/t/more-accurate-evalpoly/45932/5.
"""
function compensated_horner(p::P, x) where {T, P <: StandardBasisPolynomial{T}}
    compensated_horner(coeffs(p), x)
end

# rexpressed from paper to compute horner_sum in same pass
# sᵢ -- horner sum
# c -- compensating term
@inline function compensated_horner(ps, x)
    n, T = length(ps), eltype(ps)
    aᵢ = ps[end]
    sᵢ = aᵢ * EvalPoly._one(x)
    c = zero(T) * EvalPoly._one(x)
    for i in n-1:-1:1
	aᵢ = ps[i]
        pᵢ, πᵢ = two_product_fma(sᵢ, x)
	sᵢ, σᵢ = two_sum(pᵢ, aᵢ)
        c = fma(c, x, πᵢ + σᵢ)
    end
    sᵢ + c
end

function compensated_horner(ps::Tuple, x::S) where {S}
    ps == () && return zero(S)
    if @generated
        n = length(ps.parameters)
        sσᵢ =:(ps[end] * EvalPoly._one(x), zero(S))
        c = :(zero(S) * EvalPoly._one(x))
        for i in n-1:-1:1
            pπᵢ = :(two_product_fma($sσᵢ[1], x))
	    sσᵢ = :(two_sum($pπᵢ[1], ps[$i]))
            Σ = :($pπᵢ[2] + $sσᵢ[2])
            c = :(fma($c, x, $Σ))
        end
        s = :($sσᵢ[1] + $c)
        s
    else
        compensated_horner(ps, x)
    end
end

# error-free-transformations (EFT) of a∘b = x + y where x, y floating point, ∘ mathematical
# operator (not ⊕ or ⊗)
@inline function two_sum(a, b)
    x = a + b
    z = x - a
    y = (a - (x-z)) + (b-z)
    x, y
end

# return x, y, floats,  with x + y = a * b
@inline function two_product_fma(a, b)
    x = a * b
    y = fma(a, b, -x)
    x, y
end

# Condition number of a standard basis polynomial
# rule of thumb: p̂ a compute value
# |p(x) - p̃(x)|/|p(x)| ≤ α(n)⋅u ⋅ cond(p,x), where u = finite precision of compuation (2^-p)
function LinearAlgebra.cond(p::P, x) where {P <: StandardBasisPolynomial}
    p̃ = map(abs, p)
    p̃(abs(x))/ abs(p(x))
end