text:
# ADAPTER TRIMMING AND ALIGNMENT
echo *_1.fq.gz | parallel -n1 '[ ! -e ../alignments/${x/_1.fq.gz/.bam} ] &&
mkfifo ${x/.fq.gz/.reads.pipe} &&
mkfifo ${x/_1.fq.gz/_2.reads.pipe} &&
(cutadapt -f fastq -m 20 -a AGATCGGAAGAGC -A AGATCGGAAGAGC $x ${x/_1.fq/_2.fq} -o >(fasta trim by quality - 30 | fasta mask by quality - 20 | fasta to raw - > ${x/.fq.gz/.reads.pipe}) -p >(fasta trim by quality - 30 | fasta mask by quality - 20 | fasta to raw - > ${x/_1.fq.gz/_2.reads.pipe}) &) &&
bowtie2 -r -p8 -X 1000 --score-min L,0,-0.6 -x ~/tools/bowtie2-indexes/homo_sapiens/hg38 -1 ${x/.fq.gz/.reads.pipe} -2 ${x/_1.fq.gz/_2.reads.pipe} | samblaster | samtools view -u - | samtools sort -@8 -o ../alignments/${x/_1.fq.gz/.bam} -;
rm ${x/.fq.gz/.reads.pipe} ${x/_1.fq.gz/_2.reads.pipe}'
echo *.bam | parallel -n8 '[ ! -e ${x}.bai ] && samtools index $x'

# ESTIMATE PANEL SIZE BY CALCULATING COVERAGE HISTOGRAM FOR ENTIRE GENOME
echo *-UMI.bam | parallel -n8 'samtools view -b -q10 -F 0x400 $x | bedtools genomecov -ibam stdin -g ~/homo_sapiens/hg38.chrom.sizes > ../coverage_histograms/${x/.bam/.tsv}'

tsv_files = readdir();
histogram = zeros(2000, length(tsv_files));
total = zeros(1, length(tsv_files));
human_chr = vcat(map(x -> "chr$(x)", 1:22), ["chrX", "chrY"]);
for (s, tsv_file) in enumerate(tsv_files)
	d = readtsv(tsv_file)
	for chr in human_chr
		total[s] += d[findfirst(x -> x == chr, d[:, 1]), 4]
		for r in find((d[:, 1] .== chr) & (d[:, 2] .< 2000))
			histogram[d[r, 2] + 1, s] += d[r, 3]
		end
	end
end

threshold = 100;
for s in 1:length(tsv_files)
	captured = round(Int, total[s] - sum(histogram[1:threshold, s]))
	println("$(tsv_files[s]): $(captured)")
end
median(map(s -> total[s] - sum(histogram[1:threshold, s]),
	1:length(tsv_files)))
total[1] - sum(median(histogram, 2)[1:threshold])

# IDENTIFY SOMATIC MUTATIONS AND GERMLINE VARIANTS IN TARGETED PANEL
echo X `seq 22` Y | parallel -n5 'mutato call2 --region=chr${x} --alt-reads=8 --alt-frac=0.005 --min-mapq=0 ~/homo_sapiens/hg38.fa *.bam > ../mutations/chr${x}.vcf'
cat chr1.vcf <(cat chr{2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y}.vcf | grep -v 'CHROM\|##') > variants.vcf
variant nearby indels variants.vcf | variant somatic --alt-reads=10 --alt-frac=0.01 --test-ref-ratio=3 --test-bg-ratio=25 --ref-reads=20 --min-sidedness=15 --min-mapq=10 - ../tumor_normal_pairs.txt | variant predict effect - | variant mappability - ~/homo_sapiens/mappability_170bp_hg38.jld 90 | variant discard sketchy silent - | variant annotate - ~/homo_sapiens/cosmic_77_hg38.jld | variant annotate - ~/homo_sapiens/exac_0.3_hg38.jld | variant discard if frequency above - ~/homo_sapiens/exac_0.3_hg38.jld 0.005 > somatic.vcf

# SUPPLEMENTARY TABLE 1: LIQUID BIOPSIES
using Distributions
d = readtsv("table_s4_mutations.tsv");
headers = d[1, :][:]; d = d[2:end, :];
samples = unique(d[:, 1]);
chromosome = d[:, findone(headers .== "Chrom")];
total_reads = d[:, findone(headers .== "Total reads")];
maf = [isempty(x) ? NaN : parse(Float64, x[1:end-1]) / 100
	for x in d[:, findone(headers .== "Allele-%")]];
effect = d[:, findone(headers .== "Effect")];

# Estimate ctDNA fractions using code borrowed from "clonality" tool
cancer_frac = zeros(length(samples));
for s in 1:length(samples)
	for r in find(d[:, 1] .== samples[s])
		if isnan(maf[r]); continue; end

		ploidy = 2
		if endswith(chromosome[r], 'X') || endswith(chromosome[r], 'Y')
			ploidy = 1
		end

		# The observed mutant allele fraction is drawn from a binomial
		# distribution Bin(n, p) where p = true mutant allele fraction.
		# Here we see what the true MAF could be, assuming that the
		# observed MAF is a 95% quantile anomaly.
		true_maf = NaN;
		for p in 0:0.005:1
			dist = Binomial(total_reads[r], p)
			if cdf(dist, round(Int, maf[r] * total_reads[r])) < 0.95
				true_maf = p; break;
			end
		end

		if ploidy == 1
			this_cancer_frac = true_maf
		elseif ploidy == 2
			# Since we cannot detect copy number changes when cancer-%
			# is low, we conservatively always assume LOH for autosomes.
			this_cancer_frac = 2 / (1 / true_maf + 1)
		end

		cancer_frac[s] = max(cancer_frac[s], this_cancer_frac)
	end
end

# Output as file "ctdna_fractions.tsv"
out = open("ctdna_fractions.tsv", "w");
for s in 1:length(samples)
	@printf(out, "%s\t%.3f\n", samples[s], cancer_frac[s])
end
close(out);

# Count somatic mutations of different types
total_somatic = zeros(length(samples));
total_somatic_cds = zeros(length(samples));
total_somatic_5utr = zeros(length(samples));
total_somatic_3utr = zeros(length(samples));
for s in 1:length(samples)
	for r in find(d[:, 1] .== samples[s])
		total_somatic[s] += 1
		if is_coding_region(effect[r])
			total_somatic_cds[s] += 1
		elseif rx(effect[r], "5'-UTR")
			total_somatic_5utr[s] += 1
		elseif rx(effect[r], "3'-UTR")
			total_somatic_3utr[s] += 1
		end
	end
end

# Perform some sanity checks on the table
for s in 1:length(samples)
	if cancer_frac[s] == 0 && total_somatic[s] > 0
		info("$(samples[s]) has rearrangements but no mutations.")
	end
end

# How many total in coding regions, UTRs, and intronic/intergenic?
sum(total_somatic_cds)
sum(total_somatic_5utr) + sum(total_somatic_3utr)
2000 - sum(total_somatic_cds) - sum(total_somatic_5utr) - sum(total_somatic_3utr)

# SHARED CODE FOR COUNTING MUTATION STATISTICS
struct Region
	chromosome::String
	start::Int
	stop::Int
	len::Int
	covered_len::Int
	gene::String
	class::String
end
d = readtsv("table_s5_gene_regions.tsv"); d = d[2:end, :];
d = d[d[:, 5] .!= "N/A", :];
regions = [Region(d[r, 4], d[r, 5], d[r, 6], d[r, 6] - d[r, 5] + 1,
	d[r, 8], d[r, 1], d[r, 3]) for r in 1:size(d, 1)];

d = readtsv("sample_name_to_patient.tsv");
sample_to_patient = Dict(d[r, 1] => d[r, 2] for r in 1:size(d, 1));

struct Mutation
	chromosome::String
	position::Int
	ref_allele::String
	alt_allele::String
	alt_frac::Float64
	mut_len::Int
	in_cds::Bool
	in_utr3::Bool
	in_utr5::Bool
	gene::String
	sample::String
	patient::String
	indel::Bool
	primary::Bool   # If mutation found at multiple TP, only one is primary
end

using Variant
d = readtsv("table_s4_mutations.tsv"); d = d[2:end, :];
mutations = Vector{Mutation}(0);
for m in 1:size(d, 1)
	ref_allele = d[m, 4]; alt_allele = d[m, 5]
	position = Int(d[m, 3])
	mut_len = length(ref_allele) - length(alt_allele) + 1
	if mut_len < 1; mut_len = 1; end

	in_cds = false; in_utr3 = false; in_utr5 = false; gene = ""
	for r in regions
		if r.chromosome != d[m, 2]; continue; end
		if r.start <= position <= r.stop || r.start <= position + mut_len - 1 <= r.stop
			if r.class == "CDS"; in_cds = true; end
			if r.class == "5'-UTR"; in_utr5 = true; end
			if r.class == "3'-UTR"; in_utr3 = true; end
			gene = r.gene
		end
	end

	patient = sample_to_patient[d[m, 1]]
	alt_frac = d[m, 8] == "" ? NaN : parse(Float64, rstrip(d[m, 8], '%')) / 100
	primary = !any(m -> m.position == position && m.ref_allele == ref_allele && m.alt_allele == alt_allele && m.patient == patient, mutations)

	push!(mutations, Mutation(d[m, 2], position, ref_allele, alt_allele,
		alt_frac, mut_len, in_cds, in_utr3, in_utr5, gene,
		d[m, 1], patient, is_indel(ref_allele, alt_allele), primary))
end

# RESULTS: CALCULATE SOMATIC MUTATION RATES

# Figure 1B barplot
with_ctdna = length(unique([m.sample for m in mutations]));
total_cds_len = sum([r.len for r in regions if r.class == "CDS"]) / 1e6;
total_utr5_len = sum([r.len for r in regions if r.class == "5'-UTR"]) / 1e6;
total_utr3_len = sum([r.len for r in regions if r.class == "3'-UTR"]) / 1e6;

count(m.in_cds for m in mutations) / total_cds_len / with_ctdna
count(m.in_utr5 for m in mutations) / total_utr5_len / with_ctdna
count(m.in_utr3 for m in mutations) / total_utr3_len / with_ctdna
(utr5_mutations + utr3_mutations) / (total_utr5_len + total_utr3_len) / with_ctdna

using Distributions
1 - cdf(Binomial(cds_mutations + utr5_mutations + utr3_mutations, total_cds_len / (total_cds_len + total_utr5_len + total_utr3_len)), cds_mutations)

# FIGURE 1C AND 1D: BARPLOTS OF MUTATION COUNT
cd ~/datasets/foxa1_utr
using Plot2;
with_ctdna = 439;

# Figure 1C
d = readtsv("table_s4_mutations.tsv"); d = d[2:end, :];
cds_truncating = [fill("", 0) for g in 1:G];
cds_missense = [fill("", 0) for g in 1:G];
utr5_mutations = [fill("", 0) for g in 1:G];
utr3_mutations = [fill("", 0) for g in 1:G];
for r in 1:size(d, 1)
	sample = d[r, 1]
	g = in_region(d[r, :], "5'-UTR")
	if g > 0; push!(utr5_mutations[g], sample); end

	g = in_region(d[r, :], "3'-UTR")
	if g > 0; push!(utr3_mutations[g], sample); end

	g = findone(genes .== d[r, 6])
	if g > 0
		if is_truncating(d[r, 7])
			push!(cds_truncating[g], sample)
		elseif is_missense(d[r, 7])
			push!(cds_missense[g], sample)
		end
	end
end

order = sortperm(length.(cds_truncating) + length.(cds_missense) + length.(utr5_mutations) + length.(utr3_mutations), rev=true)[1:15];
stacked_bar_plot(hcat(length.(cds_truncating), length.(cds_missense), length.(utr5_mutations), length.(utr3_mutations))[order, :])

for g in order
	@printf("%s\t%.1f%%\n", genes[g], length(unique(vcat(cds_truncating[g], cds_missense[g], utr5_mutations[g], utr3_mutations[g]))) / with_ctdna * 100)
end

# Figure 1D
with_ctdna = 439;
genes = unique([r.gene for r in regions]);
utr3_len = [sum(r.covered_len for r in regions if r.gene == g && r.class == "3'-UTR") for g in genes];
utr3_subs = [count(m -> m.gene == g && m.in_utr3 && !m.indel, mutations)
	for g in genes];
utr3_indels = [count(m -> m.gene == g && m.in_utr3 && m.indel, mutations)
	for g in genes];

utr3_subs_rate = utr3_subs ./ (utr3_len ./ 1e6) ./ with_ctdna;
utr3_indels_rate = utr3_indels ./ (utr3_len ./ 1e6) ./ with_ctdna;

order = sortperm(utr3_subs_rate + utr3_indels_rate, rev=true)[1:15];
using Plot2; stacked_bar_plot(hcat(utr3_subs_rate, utr3_indels_rate)[order, :])

for g in order
	@printf("%s\t%.1f%%\n", genes[g], length(unique(vcat(utr3_subs[g], utr3_indels[g]))) / with_ctdna * 100)
end

# FIGURE 1E: MUTATION LOLLIPOP PLOT
cd ~/datasets/foxa1_utr
# Mutations from our own cohort
foxa1_muts = [m for m in mutations if m.gene == "FOXA1" && m.primary];
position = [round(Int, m.position + m.mut_len / 2) for m in foxa1_muts];
ref_allele = [m.ref_allele for m in foxa1_muts];
alt_allele = [m.alt_allele for m in foxa1_muts];

# Mutations from cBioPortal
d = readtsv("cbioportal_mutations.tsv"); d = d[2:end, :];
d = d[d[:, 2] .== "chr14", :];
position = Float64.(d[:, 3]);
ref_allele = d[:, 4]; alt_allele = d[:, 5];
for k in 1:length(ref_allele)
	if ref_allele[k] == "-"
		ref_allele[k] = "N"; alt_allele[k] = "N$(alt_allele[k])";
	elseif alt_allele[k] == "-"
		ref_allele[k] = "N$(ref_allele[k])"; alt_allele[k] = "N";
	end
end

# Calculate the dot positions in the scatter plot
using Plot2
mut_mrna_pos = zeros(Int, 0); mut_color = Vector{RGB}(0);
for r in 1:length(position)
	# Calculate mutation's position inside the mRNA transcript
	mrna_pos = 0
	if 37_594_901 <= position[r] <= 37_595_120         # First exon (220 bp)
		mrna_pos = 37_595_120 - position[r] + 1
	elseif 37_589_552 <= position[r] <= 37_592_711     # Second exon (3160 bp)
		mrna_pos = 37_592_711 - position[r] + 220 + 1
	else
		continue
	end

	push!(mut_mrna_pos, mrna_pos)
	if is_indel(ref_allele[r], alt_allele[r])
		if mod(length(alt_allele[r]) - length(ref_allele[r]), 3) == 0
			if length(alt_allele[r]) > length(ref_allele[r])
				push!(mut_color, RGB(255, 128, 128))
			else
				push!(mut_color, RGB(128, 128, 255))
			end
		else
			if length(alt_allele[r]) > length(ref_allele[r])
				push!(mut_color, RGB(255, 0, 0))
			else
				push!(mut_color, RGB(0, 0, 255))
			end
		end
	else
		push!(mut_color, RGB(0, 255, 0))
	end
end

order = sortperm(mut_mrna_pos);
mut_mrna_pos = mut_mrna_pos[order];
mut_color = mut_color[order];

# Stagger the mutations on the y-axis to make them not overlap
mrna_len = 3160 + 220;
y_level = fill(NaN, length(mut_mrna_pos));
y_space = falses(50, mrna_len);
for k in 1:length(mut_mrna_pos)
	level = findfirst(y_space[:, mut_mrna_pos[k]] .== false)
	window = max(mut_mrna_pos[k] - 20, 1):min(mut_mrna_pos[k] + 20, mrna_len)
	y_space[level, window] = true
	y_level[k] = level
end

figure("~/plot.pdf", size=(10, 3)) do
	scatter_plot(mut_mrna_pos, y_level, color=mut_color, marker="o", size=11)
	xlim(0, 3380); ylim(0, 40);
end

# FIGURE 2B: BASE COMPOSITION AROUND FOXA1 3'-UTR INDELS
using BioSequences, Variant
genome_path = expanduser("~/homo_sapiens/hg38.fa");
genome = FASTA.Reader(open(genome_path), index="$(genome_path).fai");
chr14 = sequence(genome["chr14"]);
foxa1_utr3_start = 37589551; foxa1_utr3_end = 37591364;
d = readtsv("table_s4_mutations.tsv"); d = d[2:end, :];
junctions = fill("", 0);
for r in 1:size(d, 1)
	if d[r, 6] != "FOXA1"; continue; end
	if !is_indel(d[r, 4], d[r, 5]); continue; end
	if any(k -> d[k, 6] == "FOXA1" && d[k, 3] == d[r, 3], 1:r-1); continue; end
	pos = Int(d[r, 3])
	if foxa1_utr3_start <= pos <= foxa1_utr3_end
		junction = "$(chr14[pos-19:pos])|$(chr14[pos+1:pos+20])"
		push!(junctions, junction)
		#println("$(junction)\t$(d[r, 4])\t$(d[r, 5])")
	end
end

gc_frac = fill(NaN, 20);
for k in 1:20
	gc_frac[k] = count(j -> j[21 - k] in "GC", junctions) + count(j -> j[21 + k] in "GC", junctions)
end
gc_frac /= 2 * length(junctions)
using Plot2; line_plot(1:length(gc_frac), gc_frac)

# Generate list of random 10 bp sequences from within the 3'-UTR
junctions = fill("", 0);
for k in 1:1000
	pos = rand(foxa1_utr3_start+50:foxa1_utr3_end-50)
	push!(junctions, "$(chr14[pos-19:pos])|$(chr14[pos+1:pos+20])")
end

# FIGURE 2C: SEQUENCE MOTIFS AROUND FOXA1 3'-UTR INDELS
using BioSequences, Variant, Plot2
genome_path = expanduser("~/homo_sapiens/hg38.fa");
genome = FASTA.Reader(open(genome_path), index="$(genome_path).fai");
chr14 = ReferenceSequence(sequence(genome["chr14"]));
foxa1_utr3_start = 37589551; foxa1_utr3_end = 37591364;

K = 4;
kmers = [DNAKmer{K}(UInt64(k)) for k in 0:(4^K-1)];

# Regions near somatic indels
kmer_count_near_indels = zeros(length(kmers)); J = 0;
d = readtsv("table_s4_mutations.tsv"); d = d[2:end, :];
for r in 1:size(d, 1)
	if d[r, 6] != "FOXA1"; continue; end
	if !is_indel(d[r, 4], d[r, 5]); continue; end
	if any(k -> d[k, 6] == "FOXA1" && d[k, 3] == d[r, 3], 1:r-1); continue; end
	pos = Int(d[r, 3])
	if !(foxa1_utr3_start <= pos <= foxa1_utr3_end); continue; end
	junction = chr14[pos-19:pos+20]
	for (_, kmer) in each(DNAKmer{K}, junction)
		kmer_count_near_indels[convert(UInt64, canonical(kmer)) + 1] += 1
		J += 1
	end
end
kmer_count_near_indels /= J;

# All regions inside FOXA1 3'-UTR
kmer_count_everywhere = zeros(length(kmers)); J = 0;
for pos in foxa1_utr3_start:(foxa1_utr3_end - 39)
	junction = chr14[pos:pos+39]
	for (_, kmer) in each(DNAKmer{K}, junction)
		kmer_count_everywhere[convert(UInt64, canonical(kmer)) + 1] += 1
		J += 1
	end
end
kmer_count_everywhere /= J;

figure("~/plot.pdf", size=(3.5,3)) do
	valid = (kmer_count_everywhere .> 0) | (kmer_count_near_indels .> 0)
	scatter_plot(kmer_count_everywhere[valid], kmer_count_near_indels[valid], marker="o", point_size=20)
	if K == 3; scale = 0.12; xlim(0, scale); ylim(0, scale); end
	if K == 4; scale = 0.035; xlim(0, scale); ylim(0, scale); end
end

order = sortperm(kmer_count_near_indels, rev=true);
hcat(kmers, kmer_count_near_indels, kmer_count_everywhere)[order[1:15], :]

# ARE FOXA1 MUTATIONS SUBCLONAL OR CLONAL?
d = readtsv("ctdna_fractions.tsv");
ctdna_frac = Dict{String, Float64}(d[r, 1] => d[r, 2] for r in 1:size(d, 1));
# Run first lines of code from Figure 1C code here...

foxa1_cds_mut = zeros(0); foxa1_utr3_mut = zeros(0);
other_cds_mut = zeros(0); other_utr3_mut = zeros(0);
for m in [m for m in mutations if m.chromosome != "chrX"]
	ctdna = ctdna_frac[m.sample]
	caf = m.alt_frac / ctdna
	if m.in_cds
		push!(m.gene == "FOXA1" ? foxa1_cds_mut : other_cds_mut, caf)
	end
	if m.in_utr3
		push!(m.gene == "FOXA1" ? foxa1_utr3_mut : other_utr3_mut, caf)
	end
end

using Plot2, Statistics;
figure("~/plot.pdf", size=(3,3)) do
	beeswarm_plot(foxa1_cds_mut, other_cds_mut, foxa1_utr3_mut, other_utr3_mut, point_size=8)
end
ranksum_test(foxa1_cds_mut, other_cds_mut)
ranksum_test(foxa1_utr3_mut, other_utr3_mut)

ranksum_test(vcat(foxa1_cds_mut, other_cds_mut), vcat(foxa1_utr3_mut, other_utr3_mut))

# SUPPLEMENTARY FIGURE: CO-OCCURRENCE OF FOXA1 AND OTHER RECURRENT MUTATIONS
using Variant
d = readtsv("table_s4_mutations.tsv"); d = d[2:end, :];
samples = unique(d[:, 1]); S = length(samples);

foxa1_cds_mutation = falses(S);
foxa1_utr3_mutation = falses(S);
tp53_mutation = falses(S);
ctnnb1_mutation = falses(S);
apc_mutation = falses(S);
pik3ca_mutation = falses(S);
pten_mutation = falses(S);
akt1_mutation = falses(S);
spop_mutation = falses(S);
ar_mutation = falses(S);

for r in 1:size(d, 1)
	gene = d[r, 6]; effect = d[r, 7];
	s = findone(samples .== d[r, 1])
	if gene == "FOXA1" && is_coding_region(effect)
		foxa1_cds_mutation[s] = true
	elseif gene == "FOXA1" && contains(effect, "3'-UTR")
		foxa1_utr3_mutation[s] = true
	elseif gene == "TP53" && is_protein_altering(effect)
		tp53_mutation[s] = true
	elseif gene == "CTNNB1" && is_protein_altering(effect)
		ctnnb1_mutation[s] = true
	elseif gene == "APC" && is_truncating(effect)
		apc_mutation[s] = true
	elseif gene == "PIK3CA" && is_protein_altering(effect)
		pik3ca_mutation[s] = true
	elseif gene == "PTEN" && is_protein_altering(effect)
		pten_mutation[s] = true
	elseif gene == "AKT1" && is_protein_altering(effect)
		akt1_mutation[s] = true
	elseif gene == "SPOP" && is_protein_altering(effect)
		spop_mutation[s] = true
	elseif gene == "AR" && is_protein_altering(effect) && d[r, 3] >= 67711620
		ar_mutation[s] = true
	end
end

using Plot2
order = hierarchical_order(foxa1_utr3_mutation, foxa1_cds_mutation, tp53_mutation, ctnnb1_mutation, apc_mutation, pik3ca_mutation, pten_mutation, akt1_mutation, spop_mutation, ar_mutation);
matrix = MatrixLayer(10, S, glyph=GLYPH_TILE);
shade(x) = x ? RGB(0, 0, 0) : RGB(230, 230, 230);
matrix.color[1, :] = map(shade, foxa1_utr3_mutation);
matrix.color[2, :] = map(shade, foxa1_cds_mutation);
matrix.color[3, :] = map(shade, tp53_mutation);
matrix.color[4, :] = map(shade, ctnnb1_mutation);
matrix.color[5, :] = map(shade, apc_mutation);
matrix.color[6, :] = map(shade, pik3ca_mutation);
matrix.color[7, :] = map(shade, pten_mutation);
matrix.color[8, :] = map(shade, akt1_mutation);
matrix.color[9, :] = map(shade, spop_mutation);
matrix.color[10, :] = map(shade, ar_mutation);
matrix_plot([matrix[:, order]], path="~/plot.svg", cell_width=26)

using Statistics
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, foxa1_cds_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, tp53_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, ctnnb1_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, apc_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, pik3ca_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, pten_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, akt1_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, spop_mutation));
@printf("%.3f\n", fisher_exact(foxa1_utr3_mutation, ar_mutation));

ref = foxa1_utr3_mutation | foxa1_cds_mutation;
@printf("%.3f\n", fisher_exact(ref, tp53_mutation));
@printf("%.3f\n", fisher_exact(ref, ctnnb1_mutation));
@printf("%.3f\n", fisher_exact(ref, apc_mutation));
@printf("%.3f\n", fisher_exact(ref, pik3ca_mutation));
@printf("%.3f\n", fisher_exact(ref, pten_mutation));
@printf("%.3f\n", fisher_exact(ref, akt1_mutation));
@printf("%.3f\n", fisher_exact(ref, spop_mutation));
@printf("%.3f\n", fisher_exact(ref, ar_mutation));

sum(foxa1_utr3_mutation | foxa1_cds_mutation)
sum(foxa1_cds_mutation)


text:
using QuantumOptics
using PyPlot

b = SpinBasis(1//2)
sx = 1/2*sigmax(b)
sy = 1/2*sigmay(b)
sz = 1/2*sigmaz(b)
XX = sx ⊗ sx
YY = sy ⊗ sy
ZZ = sz ⊗ sz

Jxy=1
Jz=1

function basis(N)
    spinchain=b
    for ispin in 2:N
        spinchain = spinchain ⊗ b

    end
    return spinchain # same local
end

function CrtModelH(N)
    Htwoint=Jxy*(XX+YY)+Jz*ZZ
    HXXZ = embed(basis(N), [1,2] , Htwoint)
    for s in 2:N-1
        HXXZ += embed(basis(N), [s,s+1] , Htwoint)
    end
    return HXXZ
end

s5z = embed(basis(10) , [5] , sz)
s6z = embed(basis(10) , [6] , sz)
sup = spinup(b)
sdo = spindown(b)
psi0 = sup
for ispin in 2:10
    global psi0
    if ispin == 5
       psi0 = psi0 ⊗ sdo
    else
       psi0 = psi0 ⊗ sup
    end
end

times = 0:.1:12
tout, psi_t = timeevolution.schroedinger(times, psi0, CrtModelH(10))

m5z = expect(s5z, psi_t)
m6z = expect(s6z, psi_t)

ext: jl

max_stars_repo_name: BlackBodyRadiation/JupyterTNN

max_issues_repo_name: BlackBodyRadiation/JupyterTNN

max_forks_repo_name: BlackBodyRadiation/JupyterTNN


text:
function raindrops(number)

end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    return readdlm(fname)
end

function Part1(instructions::AbstractArray)
    x, y = 0, 0
    for i = 1:length(instructions[:,1])
        if instructions[i,1] == "forward"
            x += instructions[i,2]
        elseif instructions[i,1] == "up"
            y -= instructions[i,2]
        elseif instructions[i,1] == "down"
            y += instructions[i,2]
        end
    end

    return x*y
end

function Part2(instructions::AbstractArray)
    x, y, aim = 0, 0, 0
    for i = 1:length(instructions[:,1])
        if instructions[i,1] == "forward"
            x += instructions[i,2]
            y += instructions[i,2]*aim
        elseif instructions[i,1] == "up"
            aim -= instructions[i,2]
        elseif instructions[i,1] == "down"
            aim += instructions[i,2]
        end
    end

    return x*y
end

instructions = ReadData("../data/day2.txt")
println(Part1(instructions))
println(Part2(instructions))


text:
@testset "BiMADS" begin
    function test(x)
        f1(x) = (x[1] + 1) .^ 2 + (x[2] - 1) .^ 2 - 10
        f2(x) = (x[1] + 12) .^ 2 + (x[2] - 3) .^ 2 + 20.0
        return f1, f2
    end

    @testset "B_points" begin
        cost = [1.0, 2.0]
        coordinates = [3.0, 4.0]
        points = DS.B_points(cost, coordinates)
        @test points.cost == cost
        @test points.x_now == coordinates
        @test points.weight == 0
    end

    @testset "BiMADS_status" begin
        s = DS.BiMADS_status()
        @test s.iteration == 0
        @test s.func_evaluation == 0
        @test s.total_time == 0.0
        @test s.hypervolume == 0.0
        @test s.opt_status == DS.Unoptimized
        @test s.opt_string == "Unoptimized"
        @test typeof(s.start_time) == Float64
    end

    @testset "Pareto Evaluation" begin
        p = DSProblem(2; objective = test)
        result = Optimize!(p)
        hv1 = hvIndicator(result)
        p = DSProblem(2; objective = test)
        SetIterationLimit(p, 100)
        result = Optimize!(p)
        hv2 = hvIndicator(result)
        @test typeof(DS.pareto_front(result)) == Vector{DS.B_points}
        @test hv1 >= hv2
    end

    @testset "Initialization" begin
        @testset "p_split" begin
            p = DSProblem(2; objective = test)
            @test p_dim(p) == 2
            # if successfully splitting BOP into two SOPs
            p1, p2 = DS.p_split(p)
            x = zeros(Float64, p.N)
            p.objective(x)
            f1, f2 = p.objective(x)
            @test p1.objective == f1
            @test p2.objective == f2
        end
        @testset "Optimization for p1 and P2" begin
            SetIterationLimit(p1, 1)
            @test_throws ErrorException DS.get_Bpoints(p1, p2, 1)
            SetFunctionEvaluationLimit(p2, 1)
            @test_throws ErrorException DS.get_Bpoints(p2, p1, 2)
            #reset p1 p2 for other tests
            p1, p2 = DS.p_split(DSProblem(2; objective = test))
            @test length(DS.get_Bpoints(p1, p2, 1)) > 0
            @test length(DS.get_Bpoints(p2, p1, 2)) > 0
        end
        p1, p2 = DS.p_split(DSProblem(2; objective = test))
        @test typeof(DS.initial_X_L(p1, p2)) == Vector{DS.B_points}
    end
    @testset "Main Iteration" begin
        @testset "Reference Point Determination" begin
            undominated_points = Vector{DS.B_points}()
            # one undominated point
            push!(undominated_points, DS.B_points([0.1, 0.2], [0.3, 0.4]))
            j, δ, ref_point = DS.ReferencePointDetermination(undominated_points)
            @test j == 1.0 && δ == 1.0 && ref_point == [0.1, 0.2]
            # two undominated point
            push!(undominated_points, DS.B_points([0.5, 0.6], [0.7, 0.8]))
            j, δ, ref_point = DS.ReferencePointDetermination(undominated_points)
            @test j == 1.0 &&
                  isapprox(0.32, δ, atol = 0.1) &&
                  ref_point == [0.5, 0.2]
        end

        @testset "Auxiliary function phi" begin
            p1, p2 = DS.p_split(DSProblem(2; objective = test))
            f1 = p1.objective
            f2 = p2.objective
            @test DS.phi(f1, f2, [16.0, 61.0], [-6.0, 1.0]) == -1.0
            @test DS.phi(f1, f2, [15.0, 61.0], [-6.0, 1.0]) ==
                  DS.phi(f1, f2, [16.0, 60.0], [-6.0, 1.0]) ==
                  0.0
            @test DS.phi(f1, f2, [14.0, 61.0], [-6.0, 1.0]) == 1.0
        end

        @testset "Stopping Conditions" begin
            @testset "Iteration Limit" begin
                p = DSProblem(2; objective = test)
                SetIterationLimit(p, 62)
                @test_throws ErrorException Optimize!(p)
                status = DS.BiMADS_status()
                p.status.optimization_status = DS.IterationLimit
                undominated_points = Vector{DS.B_points}()
                @test DS.checkBiMADSStopping(p, status, undominated_points) ==
                      false
            end
            @testset "Function Evalution Limit" begin
                p = DSProblem(2; objective = test)
                SetFunctionEvaluationLimit(p, 234)
                @test_throws ErrorException Optimize!(p)
                status = DS.BiMADS_status()
                p.status.optimization_status = DS.FunctionEvaluationLimit
                undominated_points = Vector{DS.B_points}()
                @test DS.checkBiMADSStopping(p, status, undominated_points) ==
                      false
            end
            @testset "Run time Limit" begin
                p = DSProblem(2; objective = test)
                AddStoppingCondition(p, RuntimeStoppingCondition(0.001))
                @test_throws ErrorException Optimize!(p)
                status = DS.BiMADS_status()
                p.status.optimization_status = DS.RuntimeLimit
                undominated_points = Vector{DS.B_points}()
                @test DS.checkBiMADSStopping(p, status, undominated_points) ==
                      false
            end

            @testset "Hyper-Volume Limit" begin
                p = DSProblem(2; objective = test)
                AddStoppingCondition(p, HypervolumeStoppingCondition(1.0))
                @test typeof(Optimize!(p)) == Vector{DS.B_points}
                status = DS.BiMADS_status()
                p.status.optimization_status = DS.HypervolumeLimit
                undominated_points = Vector{DS.B_points}()
                @test DS.checkBiMADSStopping(p, status, undominated_points) ==
                      false
            end

        end
    end

end


text:
using NextGP
using Test

@testset "NextGP.jl" begin

end


text:
using ProgressMeter, JLD

include("Q1-Ising2D.jl")

cd(dirname(@__FILE__))

datapath = "../../Data/Q1/"

rnum = 30
Llist = hcat(4:20)
βzoom = hcat(0.35:0.002:0.45)
EXPData = zeros(4, length(βzoom), length(Llist), rnum)
PRG = Progress(length(βzoom) * rnum)
for n ∈ 1:rnum
    for i ∈ 1:length(βzoom)
        for l ∈ 1:length(Llist)
            L = Llist[l]
            β = βzoom[i]
            Parameters = Dict(
                :dim => L,
                :MCLSize => 1000,
                :NLSize => 1000,
                :SkipNum => 200,
                :ProgBar => false)
            EXPData[:, i, l, n] = [Ising2D.IsingModel(β; Parameters...)...]
        end
        next!(PRG)
    end
end
save(datapath * "EXPData.jld", "EXPData", EXPData, "βzoom", βzoom, "Llist", Llist)


text:
module MechGlueDiffEqBase
import Base: similar
import Unitfu: AbstractQuantity, Quantity, ustrip, norm, unit, zero
import Unitfu: Dimensions, FreeUnits
import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2
import DiffEqBase: calculate_residuals, @muladd
using RecursiveArrayTools
export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, AbstractQuantity, Quantity
export norm , ArrayPartition, similar, zero

# This is identical to what DiffEqBase defines for Unitful
function value(x::Type{AbstractQuantity{T,D,U}}) where {T,D,U}
    T
end
# This is different from what DiffEqBase defines for Unitful
value(::Type{<:AbstractQuantity{T,D,U}}) where {T,D,U<:Core.TypeofBottom} = Base.undef_ref_str
value(x::Q) where {Q<:AbstractQuantity} = ustrip(x)

# This is identical to what DiffEqBase defines for Unitful
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:AbstractQuantity,N},t) where {N}
    # Support adaptive errors should be errorless for exponentiation
    sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
# This is identical to what DiffEqBase defines for Unitful
@inline function ODE_DEFAULT_NORM(u::Array{<:AbstractQuantity,N},t) where {N}
    sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
# This is identical to what DiffEqBase defines for Unitful
@inline function ODE_DEFAULT_NORM(u::AbstractQuantity, t)
    abs(ustrip(u))
end

@inline function UNITLESS_ABS2(u::AbstractArray{<:AbstractQuantity,N} where N)
    map(UNITLESS_ABS2, u)
end
@inline function UNITLESS_ABS2(u::AbstractArray{Quantity{T},N}) where {N, T}
    map(UNITLESS_ABS2, u)
end

@inline function UNITLESS_ABS2(x::T) where T <: AbstractQuantity
    xul = x / oneunit(T)
    abs2(xul)
end
@inline function UNITLESS_ABS2(x::Quantity{T, D, U}) where {T, D, U}
    xul = x / oneunit(Quantity{T, D, U})
    abs2(xul)::T
end

# Vectors with compatible units, treat as normal
zero(x::Vector{Quantity{T, D, U}}) where {T,D,U} = fill!(similar(x), zero(Quantity{T, D, U}))
# Vectors with incompatible units, special inferreable treatment
function zero(x::Vector{Q}) where {Q<:AbstractQuantity{T, D, U} where {D,U}} where T
    x0 = copy(x)
    for i in eachindex(x0)
        x = x0[i]
        x0[i] = 0 * x * sign(x)
    end
    x0
end

# Vectors with compatible units, treat as normal
similar(x::Vector{Quantity{T, D, U}}) where {T,D,U} = Vector{Quantity{T, D, U}}(undef, size(x,1))
# Vectors with incompatible units, special inferreable treatment
# VERY similar is still similar and (very slightly) different
similar(x::Vector{Q}) where {Q<:AbstractQuantity{T, D, U} where {D,U}} where T = copy(x)

# Vectors with compatible units, treat as normal
similar(x::Vector{Quantity{T, D, U}}, S::Type) where {T,D,U} = Vector{S}(undef, size(x,1))
# Vectors with incompatible units, special inferreable treatment
function similar(x::Vector{Q}, S::Type) where {Q<:AbstractQuantity{T, D, U} where {D, U}} where T
    x0 = Vector{S}(undef, size(x, 1))
end

# KISS pre-compillation to reduce loading times
# This is simply a boiled-down obfuscated test_4.jl

import Unitfu: m, s, kg, N, ∙
let
    r0ul = [1131.340, -2282.343, 6672.423]
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = ArrayPartition(r0ul,v0ul)
    ODE_DEFAULT_NORM(rv0ul, 0.0)
    r0 = [1131.340, -2282.343, 6672.423]∙kg
    v0 = [-5.64305, 4.30333, 2.42879]∙kg/s
    rv0 = ArrayPartition(r0, v0)
    ODE_DEFAULT_NORM(rv0, 0.0s)
    r0 = [1.0kg, 2.0N, 3.0m/s, 4.0m/s]
    v0 = [1.0kg/s, 2.0N/s, 3.0m/s^2, 4m/s^2]
    rv0 = ArrayPartition(r0, v0)
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = ArrayPartition(r0ul, v0ul)
    r0 = [1131.340, -2282.343, 6672.423]∙kg
    r1 = [1kg, 2.0m]
    zero(r0)
    zero(r1)
    rv0 = ArrayPartition(r0)
    zero(rv0)
    rv1 = ArrayPartition(r1)
    zero(rv1)
    r0 = [1.0kg, -2kg, 3m/s, 4m/s]
    zero(r0)
    rv0 = ArrayPartition(r0)
    zer = zero(rv0)
    zer == [0.0kg, 0.0kg, 0.0m/s, 0.0m/s]
    r0 = [1131.340, -2282.343, 6672.423]∙kg
    simi = similar(r0)
    rv0 = ArrayPartition(r0)
    sima = similar(rv0)
    r0 = [1.0kg, -2kg, 3m/s, 4m/s]
    similar(r0)
    rv0 = ArrayPartition(r0)
    sima = similar(rv0)
    r0ul = [1131.340, -2282.343, 6672.423]
    v0ul = [-5.64305, 4.30333, 2.42879]
    rv0ul = ArrayPartition(r0ul,v0ul)
    UNITLESS_ABS2(rv0ul)
    r0 = [1.0kg, 2.0N, 3.0m/s, 4.0m/s]
    v0 = [1.0kg/s, 2.0N/s, 3.0m/s^2, 4m/s^2]
    rv0 = ArrayPartition(r0, v0)
    UNITLESS_ABS2(1.0kg)
    UNITLESS_ABS2(r0)
    UNITLESS_ABS2(rv0)
    nothing
end

end


text:
#!/usr/bin/env julia

using Test, DynamicBoundsBase

const DEqR = DynamicBoundsBase
const TSC = DEqR.TerminationStatusCode

@testset "Integrator Attribute Constructor" begin
    t = TimeIndex(3)
    @test Gradient{Lower}() == Gradient{Lower}(-1, -Inf)
    @test Gradient{Upper}() == Gradient{Upper}(-1, -Inf)
    @test Gradient{Nominal}() == Gradient{Nominal}(-1, -Inf)
    @test Gradient() == Gradient{Undefined}(-1, -Inf)
    @test Gradient{Nominal}(t). index == 3
    @test Gradient{Nominal}(2.2).time == 2.2

    ref_attr = Gradient()
    ref = Ref(ref_attr)
    @test Base.broadcastable(ref_attr).x == ref.x

    @test Subgradient{Lower}() == Subgradient{Lower}(-1, -Inf)
    @test Subgradient{Upper}() == Subgradient{Upper}(-1, -Inf)
    @test Subgradient() == Subgradient{Undefined}(-1, -Inf)
    @test Subgradient{Lower}(t). index == 3
    @test Subgradient{Lower}(2.2).time == 2.2

    @test Bound{Lower}() == Bound{Lower}(-1, -Inf)
    @test Bound{Upper}() == Bound{Upper}(-1, -Inf)
    @test Bound() == Bound{Undefined}(-1, -Inf)
    @test Bound{Lower}(t). index == 3
    @test Bound{Lower}(2.2).time == 2.2

    @test Relaxation{Lower}() == Relaxation{Lower}(-1, -Inf)
    @test Relaxation{Upper}() == Relaxation{Upper}(-1, -Inf)
    @test Relaxation() == Relaxation{Undefined}(-1, -Inf)
    @test Relaxation{Lower}(t). index == 3
    @test Relaxation{Lower}(2.2).time == 2.2

    @test ParameterValue() == ParameterValue(-1)

    @test ParameterBound{Lower}() == ParameterBound{Lower}(-1)
    @test ParameterBound{Upper}() == ParameterBound{Upper}(-1)
    @test ParameterBound() == ParameterBound{Undefined}(-1)
end

struct UndefinedIntegrator <: DEqR.AbstractDERelaxIntegrator end
struct UndefinedProblem <: DEqR.AbstractDERelaxProblem end

mutable struct TestIntegrator <: DEqR.AbstractODERelaxIntegrator
    temp::Float64
end

mutable struct TestProblem <: DEqR.AbstractDERelaxProblem
    temp::Float64
end

struct new_integrator_attribute <: DEqR.AbstractIntegratorAttribute
end

struct new_problem_attribute <: DEqR.AbstractRelaxProblemAttribute
end

@testset "Integrator Attributes Interface" begin

    undefined_problem = UndefinedProblem()
    DEqR.supports(undefined_problem, DEqR.ConstantStateBounds()) == false

    DEqR.supports(::TestIntegrator, ::DEqR.IntegratorName) = true
    DEqR.supports(::TestIntegrator, ::DEqR.Gradient) = true
    DEqR.supports(::TestIntegrator, ::DEqR.Subgradient) = true
    DEqR.supports(::TestIntegrator, ::DEqR.Bound) = true
    DEqR.supports(::TestIntegrator, ::DEqR.Relaxation) = true
    DEqR.supports(::TestIntegrator, ::DEqR.IsNumeric) = true
    DEqR.supports(::TestIntegrator, ::DEqR.IsSolutionSet) = true
    DEqR.supports(::TestIntegrator, ::DEqR.TerminationStatus) = true
    DEqR.supports(::TestIntegrator, ::DEqR.Value) = true

    DEqR.get(t::TestIntegrator, a::DEqR.IntegratorName) = "TestIntegrator"
    DEqR.get(t::TestIntegrator, a::DEqR.Gradient) = t.temp
    DEqR.get(t::TestIntegrator, a::DEqR.Subgradient) = t.temp
    DEqR.get(t::TestIntegrator, a::DEqR.Bound) = t.temp
    DEqR.get(t::TestIntegrator, a::DEqR.Relaxation) = t.temp
    DEqR.get(t::TestIntegrator, a::DEqR.IsNumeric) = true
    DEqR.get(t::TestIntegrator, a::DEqR.IsSolutionSet) = true
    DEqR.get(t::TestIntegrator, a::DEqR.TerminationStatus) = t.temp
    DEqR.get(t::TestIntegrator, a::DEqR.Value) = t.temp

    DEqR.set!(t::TestIntegrator, a::DEqR.Gradient, value) = (t.temp = value)
    DEqR.set!(t::TestIntegrator, a::DEqR.Subgradient, value) = (t.temp = value)
    DEqR.set!(t::TestIntegrator, a::DEqR.Bound, value) = (t.temp = value)
    DEqR.set!(t::TestIntegrator, a::DEqR.Relaxation, value) = (t.temp = value)
    DEqR.set!(t::TestIntegrator, a::DEqR.TerminationStatus, value) = (t.temp = value)
    DEqR.set!(t::TestIntegrator, a::DEqR.Value, value) = (t.temp = value)

    DEqR.supports(::TestProblem, ::DEqR.Value) = true
    DEqR.get(t::TestProblem, a::DEqR.Value) = t.temp

    undefined_integrator = UndefinedIntegrator()
    @test !DEqR.supports(undefined_integrator, DEqR.IntegratorName())
    @test !DEqR.supports(undefined_integrator, DEqR.Gradient{DEqR.Lower}(1.0))
    @test !DEqR.supports(undefined_integrator, DEqR.Subgradient{DEqR.Lower}(1.0))
    @test !DEqR.supports(undefined_integrator, DEqR.Bound{DEqR.Lower}(1.0))
    @test !DEqR.supports(undefined_integrator, DEqR.Relaxation{DEqR.Lower}(1.0))
    @test !DEqR.supports(undefined_integrator, DEqR.IsNumeric())
    @test !DEqR.supports(undefined_integrator, DEqR.IsSolutionSet())
    @test !DEqR.supports(undefined_integrator, DEqR.TerminationStatus())
    @test !DEqR.supports(undefined_integrator, DEqR.Value())

    test_integrator = TestIntegrator(1.0)
    @test @inferred DEqR.supports(test_integrator, DEqR.IntegratorName())
    @test @inferred DEqR.supports(test_integrator, DEqR.Gradient{DEqR.Lower}(1.0))
    @test @inferred DEqR.supports(test_integrator, DEqR.Subgradient{DEqR.Lower}(1.0))
    @test @inferred DEqR.supports(test_integrator, DEqR.Bound{DEqR.Lower}(1.0))
    @test @inferred DEqR.supports(test_integrator, DEqR.Relaxation{DEqR.Lower}(1.0))
    @test @inferred DEqR.supports(test_integrator, DEqR.IsNumeric())
    @test @inferred DEqR.supports(test_integrator, DEqR.IsSolutionSet())
    @test @inferred DEqR.supports(test_integrator, DEqR.TerminationStatus())
    @test @inferred DEqR.supports(test_integrator, DEqR.Value())

    @test_nowarn @inferred DEqR.set!(test_integrator, DEqR.Gradient{DEqR.Lower}(1.1), 1.2)
    val = @inferred DEqR.get(test_integrator, DEqR.Gradient{DEqR.Lower}(1.1))
    @test val === 1.2
    @test_nowarn @inferred DEqR.set!(test_integrator, DEqR.Subgradient{DEqR.Lower}(1.1), 1.3)
    val = @inferred DEqR.get(test_integrator, DEqR.Subgradient{DEqR.Lower}(1.1))
    @test val === 1.3
    @test_nowarn @inferred DEqR.set!(test_integrator, DEqR.Bound{DEqR.Lower}(1.1), 1.4)
    val = @inferred DEqR.get(test_integrator, DEqR.Bound{DEqR.Lower}(1.1))
    @test val === 1.4
    @test_nowarn @inferred DEqR.set!(test_integrator, DEqR.Relaxation{DEqR.Lower}(1.1), 1.5)
    val = @inferred DEqR.get(test_integrator, DEqR.Relaxation{DEqR.Lower}(1.1))
    @test val === 1.5
    @test_nowarn DEqR.set!(test_integrator, DEqR.Value(), 1.53)
    val = @inferred DEqR.get(test_integrator, DEqR.Value())
    @test val === 1.53

    sval = @inferred DEqR.get(test_integrator, DEqR.IntegratorName())
    @test sval === "TestIntegrator"
    @test @inferred DEqR.get(test_integrator, DEqR.IsNumeric())
    @test @inferred DEqR.get(test_integrator, DEqR.IsSolutionSet())

    @test_nowarn @inferred DEqR.set!(test_integrator, DEqR.TerminationStatus(), 1.9)
    val = @inferred DEqR.get(test_integrator, DEqR.TerminationStatus())
    @test val === 1.9

    @test_throws ArgumentError DEqR.get(UndefinedIntegrator(), DEqR.IntegratorName())

    vsbnds1 = DEqR.VariableStateBounds()
    @test vsbnds1.xL == Base.isempty
    @test vsbnds1.xU == Base.isempty

    vsbnds2 = DEqR.VariableStateBounds(x -> 1.0*x, x-> 2.0*x)
    @test vsbnds2.xL(1.0) == 1.0
    @test vsbnds2.xU(1.1) == 2.2

    cbnds1 = DEqR.ConstantStateBounds()
    @test isempty(cbnds1.xL)
    @test isempty(cbnds1.xU)

    cbnds2 = DEqR.ConstantStateBounds([1.0], [2.2])
    @test cbnds2.xL[1] == 1.0
    @test cbnds2.xU[1] == 2.2

    ref_attr = DEqR.VariableStateBounds()
    ref = Ref(ref_attr)
    @test Base.broadcastable(ref_attr).x == ref.x

    pconstr = DEqR.PolyhedralConstraint([1.0 1.0; 2.0 2.1], [3.2; 3.1])
    @test pconstr.A == [1.0 1.0; 2.0 2.1]
    @test pconstr.b == [3.2; 3.1]

    states = DEqR.IntegratorStates()
    @test !states.first_pnt_eval
    @test states.new_decision_box
    @test states.new_decision_pnt
    @test states.termination_status === RELAXATION_NOT_CALLED

    out = Any[]
    new_pa = new_problem_attribute()
    new_ia = new_integrator_attribute()

    @test_throws ArgumentError DEqR.get(undefined_problem, new_pa)
    @test_throws ArgumentError DEqR.get(undefined_integrator, new_ia)

    v = Float64[1.0; 2.0]
    @test_throws ArgumentError DEqR.get(undefined_problem, new_pa)
    @test_throws ArgumentError DEqR.get(undefined_integrator, new_ia)

    @test_throws ArgumentError DEqR.getall!(out, undefined_problem, new_pa)
    @test_throws ArgumentError DEqR.getall!(out, undefined_integrator, new_ia)

    @test_throws ArgumentError DEqR.getall(undefined_problem, new_pa)
    @test_throws ArgumentError DEqR.getall(undefined_integrator, new_ia)

    @test_throws DynamicBoundsBase.UnsupportedRelaxAttribute{new_problem_attribute} DEqR.set!(undefined_problem, new_pa, 1)
    @test_throws DynamicBoundsBase.UnsupportedRelaxAttribute{new_integrator_attribute} DEqR.set!(undefined_integrator, new_ia, 1)

    v = Float64[1.0; 2.0]
    @test_throws DimensionMismatch DEqR.set!(undefined_problem, new_pa, v, [1])
    @test_throws DimensionMismatch DEqR.set!(undefined_integrator, new_ia, v, [1])

    v = Float64[1.0; 2.0]
    @test_throws DEqR.UnsupportedRelaxAttribute{new_problem_attribute}(new_problem_attribute(), "") DEqR.set!(undefined_problem, new_pa, v, [1; 2])
    @test_throws DEqR.UnsupportedRelaxAttribute{new_integrator_attribute}(new_integrator_attribute(), "") DEqR.set!(undefined_integrator, new_ia, v, [1; 2])

    @test DEqR.operation_name(DEqR.SetRelaxAttributeNotAllowed{Gradient{Lower}}(Gradient{Lower}(), "test_me")) == "Setting attribute Gradient{Lower}(-1, -Inf)"
    @test DEqR.message(DEqR.SetRelaxAttributeNotAllowed{Gradient{Lower}}(Gradient{Lower}(), "test_me")) == "test_me"

    @test DEqR.UnsupportedRelaxAttribute(DEqR.UnsupportedRelaxAttribute{Gradient{Lower}}(Gradient{Lower}(),"aaa")) == "Attribute Gradient{Lower}(-1, -Inf)"

    @test_throws DEqR.SetRelaxAttributeNotAllowed{Gradient{Nominal}}(Gradient{Nominal}(-1, -Inf), "") DEqR.setall!(test_integrator, DEqR.Gradient{Nominal}(), [1.0; 2.0; 3.0])
end

@testset "ODE Relax Problem" begin

    x0(p) = [9.0]
    xL = [0.0]
    xU = [10.0]

    function f!(dx, x, p, t)
        dx[1] = p[1] - x[1]^2
        nothing
    end

    function Jx!(dx, x, p, t)
        dx[1] = -2.0*x[1]
        nothing
    end

    function Jp!(dx, x, p, t)
        dx[1] = one(p[1])
        nothing
    end

    tspan = (0.0, 0.1)
    pL = [-1.0]
    pU = [1.0]
    pval = [0.1]
    tsupports = [0.1; 0.2; 0.3]

    prob = DEqR.ODERelaxProb(f!, tspan, x0, pL, pU, xL = xL, xU = xU, Jx! = Jx!,
                             Jp! = Jp!, p = pval, tsupports = tsupports, nx = 1)

    @test DEqR.supports(prob, DEqR.HasStateBounds())
    @test DEqR.supports(prob, DEqR.HasConstantStateBounds())
    @test DEqR.supports(prob, DEqR.HasVariableStateBounds())
    @test DEqR.supports(prob, DEqR.HasUserJacobian())
    @test DEqR.supports(prob, DEqR.ConstantStateBounds())
    @test DEqR.supports(prob, DEqR.PolyhedralConstraint())
    @test DEqR.supports(prob, DEqR.SupportSet())

    @test DEqR.get(prob, DEqR.HasStateBounds())
    @test DEqR.get(prob, DEqR.HasConstantStateBounds())
    @test DEqR.get(prob, DEqR.HasVariableStateBounds())
    @test DEqR.get(prob, DEqR.HasUserJacobian())

    cbnds = DEqR.ConstantStateBounds([1.0], [2.2])
    pconstr = DEqR.PolyhedralConstraint([1.0 1.0; 2.0 2.1], [3.2; 3.1])

    @test_nowarn DEqR.set!(prob, cbnds)
    @test_nowarn DEqR.set!(prob, pconstr)

    DEqR.set!(prob, cbnds)
    DEqR.set!(prob, pconstr)
    @test prob.constant_state_bounds.xL[1] == 1.0
    @test prob.polyhedral_constraint.A[2,2] == 2.1

    csbnds = DEqR.get(prob, DEqR.ConstantStateBounds())
    pconstrs = DEqR.get(prob, DEqR.PolyhedralConstraint())

    @test csbnds.xL[1] == 1.0
    @test csbnds.xU[1] == 2.2

    @test pconstrs.A[2,2] == 2.1
    @test pconstrs.b[1] == 3.2

    prob2 = DEqR.ODERelaxProb(f!, tspan, x0, pL, pU)
    @test prob2.p == 0.5*(pL + pU)

    DEqR.set!(prob2, DEqR.SupportSet([1.0; 2.0]))
    support_set = DEqR.get(prob2, DEqR.SupportSet())
    @test support_set.s[1] == 1.0
    @test support_set.s[2] == 2.0
end


text:
"""
A MappedGrid consists of a grid and a map. Each grid point of the mapped grid
is the map of the corresponding point of the underlying grid.
"""
struct MappedGrid{G,M,T,N} <: AbstractGrid{T,N}
	supergrid	::	G
	map			::	M

	MappedGrid{G,M,T,N}(supergrid::AbstractGrid{T,N}, map) where {G,M,T,N} = new(supergrid, map)
end

const MappedGrid1d{G,M,T<:Number,N} = MappedGrid{G,M,T,N}

MappedGrid(grid::AbstractGrid{T,N}, map::AbstractMap) where {T,N} =
	MappedGrid{typeof(grid),typeof(map),T,N}(grid, map)

name(grid::MappedGrid) = "Mapped grid"

supergrid(g::MappedGrid) = g.supergrid

mapping(g::MappedGrid) = g.map

mapped_grid(grid::AbstractGrid, map::AbstractMap) = MappedGrid(grid, map)

# avoid multiple mappings
mapped_grid(g::MappedGrid, map::AbstractMap) = MappedGrid(supergrid(g), map∘mapping(g))

# Convenience function, similar to apply_map for Dictionary's
apply_map(grid::AbstractGrid, map::AbstractMap) = mapped_grid(grid, map)

for op in (:length, :size, :eachindex, :indextype, :isperiodic)
	@eval $op(g::MappedGrid) = $op(supergrid(g))
end

for op in (:minimum, :maximum, :support)
	@eval $op(g::MappedGrid1d) = mapping(g)*$op(supergrid(g))
end

resize(g::MappedGrid, n::Int) = apply_map(resize(supergrid(g), n), mapping(g))

unsafe_getindex(g::MappedGrid, idx) = applymap(g.map, g.supergrid[idx])

function rescale(g::AbstractGrid1d, a, b)
	m = interval_map(endpoints(support(g))..., a, b)
	mapped_grid(g, m)
end

# Preserve tensor product structure
function rescale(g::ProductGrid, a::SVector{N}, b::SVector{N}) where {N}
	scaled_grids = [ rescale(element(g, i), a[i], b[i]) for i in 1:N]
	ProductGrid(scaled_grids...)
end


text:
export EpsilonGreedyPolicy

"Policy that acts uniformly at random with ϵ chance, but is otherwise greedy."
struct EpsilonGreedyPolicy{P <: PolicySolution, R <: AbstractRNG} <: PolicySolution
    domain::Domain
    policy::P
    epsilon::Float64
    rng::R
end

EpsilonGreedyPolicy(domain, policy, epsilon) =
    EpsilonGreedyPolicy(domain, policy, epsilon, Random.GLOBAL_RNG)

get_action(sol::EpsilonGreedyPolicy, state::State) =
    rand_action(sol, state)
best_action(sol::EpsilonGreedyPolicy, state::State) =
    best_action(sol.policy, state)
get_value(sol::EpsilonGreedyPolicy, state::State) =
    get_value(sol.policy, state)
get_value(sol::EpsilonGreedyPolicy, state::State, action::Term) =
    get_value(sol.policy, state, action)
get_action_values(sol::EpsilonGreedyPolicy, state::State) =
    get_action_values(sol.policy, state)

function rand_action(sol::EpsilonGreedyPolicy, state::State)
    best_act = best_action(sol.policy, state)
    actions = [available(state, sol.domain); best_act]
    probs = ones(length(actions)) * sol.epsilon
    probs[end] = 1 - sol.epsilon
    return sample(sol.rng, actions, Weights(probs, 1.0))
end


text:
using AWSCore.Services: route53

@testset "Endpoints" begin
    @testset "Route53" begin
        try
            route53("GET", "/2013-04-01/accountlimit/MAX_HEALTH_CHECKS_BY_OWNER")
            @test true
        catch
            @test false
        end
    end
end

@testset "rest-json" begin
    @testset "glacier - special case" begin
        vault_names = [string("test-vault-", UUIDs.uuid4()), string("test-vault-", UUIDs.uuid4())]

        @testset "Create Vaults" begin
            for vault in vault_names
                Services.glacier("PUT", "/-/vaults/$vault")
            end
        end

        @testset "Get Vaults" begin
            # If this is an Integer AWS Coral cannot convert it to a String
            # "class com.amazon.coral.value.json.numbers.TruncatingBigNumber can not be converted to an String"
            limit = "1"
            headers = Dict("foo"=>"bar")

            cases = [
                Dict("limit"=>limit, "headers"=>Dict("foo"=>"bar")),
                Dict("limit"=>limit, "headers"=>("foo"=>"bar")),
                Dict("limit"=>limit, "headers"=>[("foo"=>"bar"), ("biz"=>"baz")]),
                Dict("limit"=>limit),
                ("limit"=>limit, "headers"=>Dict("foo"=>"bar")),
                ("limit"=>limit, "headers"=>("foo"=>"bar")),
                ("limit"=>limit, "headers"=>[("foo"=>"bar"), ("biz"=>"baz")]),
                ("limit"=>limit),
                [("limit"=>limit), ("headers"=>Dict("foo"=>"bar"))],
                [("limit"=>limit), ("headers"=>("foo"=>"bar"))],
                [("limit"=>limit), ("headers"=>[("foo"=>"bar"), ("biz"=>"baz")])],
                [("limit"=>limit)]
            ]

            @testset "$case" for case in cases
                result = Services.glacier("GET", "/-/vaults", case)

                @test length(result["VaultList"]) == parse(Int, limit)
            end
        end

        @testset "Delete Vaults" begin
            for vault in vault_names
                Services.glacier("DELETE", "/-/vaults/$vault")
            end

            result = Services.glacier("GET", "/-/vaults")
            vaults = [v["VaultName"] for v in result["VaultList"]]

            for vault in vault_names
                @test !(vault in vaults)
            end
        end
    end
end


text:
w = 500
h = 300
a = Image(w, h)

purple = Color(0xff, 0, 0xff)
green = Color(0, 0xff, 0)
white = Color(0xff, 0xff, 0xff)

fill!(a, green)
for i in 20:220, j in 10:100
    splat!(a, i, j, purple)
end
for i in 180:400, j in 80:200
    splat!(a, i, j, white)
end

fn = "bitmap_write.ppm"
writeppm(fn, a)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using NearestNeighbors, Statistics, Distances

export estimate_dimension, stochastic_indicator
export Euclidean, Chebyshev, Cityblock

#####################################################################################
#                                Estimate Dimension                                 #
#####################################################################################
"""
    estimate_dimension(s::AbstractVector, τ::Int, γs = 1:5, method = "afnn"; kwargs...)

Compute a quantity that can estimate an optimal amount of
temporal neighbors `γ` to be used in [`reconstruct`](@ref) or [`embed`](@ref).

## Description
Given the scalar timeseries `s` and the embedding delay `τ` compute a quantity
for each `γ ∈ γs` based on the "nearest neighbors" in the embedded time series.

The quantity that is calculated depends on the algorithm defined by the string `method`:

* `"afnn"` (default) is Cao's "Averaged False Nearest Neighbors" method [1], which
    gives a ratio of distances between nearest neighbors. This ratio saturates
    around `1.0` near the optimal value of `γ` (see [`afnn`](@ref)).
* `"fnn"` is Kennel's "False Nearest Neighbors" method [2], which gives the
    number of points that cease to be "nearest neighbors" when the dimension
    increases. This number drops down to zero near the optimal value of `γ`.
    This method accepts the keyword arguments `rtol` and `atol`, which stand
    for the "tolerances" required by Kennel's algorithm (see [`fnn`](@ref)).
* `"f1nn"` is Krakovská's "False First Nearest Neighbors" method [3], which
    gives the ratio of pairs of points that cease to be "nearest neighbors"
    when the dimension increases. This number drops down to zero near the
    optimal value of `γ` (see [`f1nn`](@ref)).

`"afnn"` and `"f1nn"` also support the `metric` keyword, which can be any of
`Cityblock(), Euclidean(), Chebyshev()`. This metric is used both
for computing the nearest neighbors (`KDTree`s) as well as the distances necessary for
Cao's method (eqs. (2, 3) of [1]). Defaults to `Euclidean()` (note that [1] used
`Chebyshev`).

Please be aware that in **DynamicalSystems.jl** `γ` stands for the amount of temporal
neighbors and not the embedding dimension (`D = γ + 1`, see also [`embed`](@ref)).

## References

[1] : Liangyue Cao, [Physica D, pp. 43-50 (1997)](https://www.sciencedirect.com/science/article/pii/S0167278997001188?via%3Dihub)

[2] : M. Kennel *et al.*, [Phys. Review A **45**(6), 3403-3411](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.45.3403)
(1992).

[3] : Anna Krakovská *et al.*, [J. Complex Sys. 932750 (2015)](https://doi.org/10.1155/2015/932750)
"""
function estimate_dimension(s::AbstractVector, τ::Int, γs = 1:5, method = "afnn";
    metric = Euclidean(), kwargs...)

    if method == "afnn"
        return afnn(s, τ, γs, metric)
    elseif method == "fnn"
        return fnn(s, τ, γs; kwargs...)
    elseif method == "f1nn"
        return f1nn(s, τ, γs, metric)
    end
end

"""
    afnn(s::AbstractVector, τ:Int, γs = 1:5, metric=Euclidean())

Compute the parameter E₁ of Cao's "averaged false nearest neighbors" method for
determining the minimum embedding dimension of the time series `s`, with
a sequence of `τ`-delayed temporal neighbors [1].

## Description
Given the scalar timeseries `s` and the embedding delay `τ` compute the
values of `E₁` for each `γ ∈ γs`, according to Cao's Method (eq. 3 of [1]).

This quantity is a ratio of the averaged distances between the nearest neighbors
of the reconstructed time series, which quantifies the increment of those
distances when the number of temporal neighbors changes from `γ` to `γ+1`.

Return the vector of all computed `E₁`s. To estimate a good value for `γ` from this,
find `γ` for which the value `E₁` saturates at some value around 1.

*Note: This method does not work for datasets with perfectly periodic signals.*

See also: [`estimate_dimension`](@ref), [`fnn`](@ref), [`f1nn`](@ref).

## References

[1] : Liangyue Cao, [Physica D, pp. 43-50 (1997)](https://www.sciencedirect.com/science/article/pii/S0167278997001188?via%3Dihub)
"""
function afnn(s::AbstractVector{T}, τ::Int, γs = 1:5, metric=Euclidean()) where {T}
    E1s = zeros(length(γs))
    aafter = 0.0
    aprev = _average_a(s, γs[1], τ, metric)
    for (i, γ) ∈ enumerate(γs)
        aafter = _average_a(s, γ+1, τ, metric)
        E1s[i] = aafter/aprev
        aprev = aafter
    end
    return E1s
end
# then use function `saturation_point(γs, E1s)` from ChaosTools

function _average_a(s::AbstractVector{T},γ,τ,metric) where {T}
    #Sum over all a(i,d) of the Ddim Reconstructed space, equation (2)
    Rγ = reconstruct(s[1:end-τ],γ,τ)
    tree2 = KDTree(Rγ, metric)
    nind = (x = knn(tree2, Rγ.data, 2)[1]; [ind[1] for ind in x])
    e = 0.0
    @inbounds for (i,j) ∈ enumerate(nind)
        δ = evaluate(metric, Rγ[i], Rγ[j])
        #If Rγ[i] and Rγ[j] are still identical, choose the next nearest neighbor
        if δ == 0.0
            j = knn(tree2, Rγ[i], 3, true)[1][end]
            δ = evaluate(metric, Rγ[i], Rγ[j])
        end
        e += _increase_distance(δ,s,i,j,γ,τ,metric)/δ
    end
    return e / (length(Rγ)-1)
end

# Function to increase the distance (p-norm) between two points `(i,j)` of
# the embedded time `s`series, by adding one temporal neighbor
_increase_distance(δ, s, i::Int, j::Int, γ::Int, τ::Int, ::Chebyshev) =
    max(δ, abs(s[i+γ*τ+τ] - s[j+γ*τ+τ]))
_increase_distance(δ, s, i::Int, j::Int, γ::Int, τ::Int, ::Euclidean) =
    sqrt(δ*δ + abs2(s[i+γ*τ+τ] - s[j+γ*τ+τ]) )
_increase_distance(δ, s, i::Int, j::Int, γ::Int, τ::Int, ::Cityblock) =
    δ + abs(s[i+γ*τ+τ] - s[j+γ*τ+τ])

"""
    fnn(s::AbstractVector, τ:Int, γs = 1:5; rtol=10.0, atol=2.0)

Calculate the number of "false nearest neighbors" (FNN) of the datasets created
from `s` with a sequence of `τ`-delayed temporal neighbors.

## Description
Given a dataset made by embedding `s` with `γ` temporal neighbors and delay `τ`,
the "false nearest neighbors" (FNN) are the pairs of points that are nearest to
each other at dimension `γ`, but are separated at dimension `γ+1`. Kennel's
criteria for detecting FNN are based on a threshold for the relative increment
of the distance between the nearest neighbors (`rtol`, eq. 4 in [1]), and
another threshold for the ratio between the increased distance and the
"size of the attractor" (`atol`, eq. 5 in [1]). These thresholds are given
as keyword arguments.

The returned value is a vector with the number of FNN for each `γ ∈ γs`. The
optimal value for `γ` is found at the point where the number of FNN approaches
zero.

See also: [`estimate_dimension`](@ref), [`afnn`](@ref), [`f1nn`](@ref).

## References

[1] : M. Kennel *et al.*, "Determining embedding dimension for phase-space
reconstruction using a geometrical construction", *Phys. Review A 45*(6), 3403-3411
(1992).
"""
function fnn(s::AbstractVector, τ::Int, γs = 1:5; rtol=10.0, atol=2.0)
    rtol2 = rtol^2
    Ra = std(s, corrected=false)
    nfnn = zeros(length(γs))
    @inbounds for (k, γ) ∈ enumerate(γs)
        y = reconstruct(s[1:end-τ],γ,τ)
        tree = KDTree(y)
        nind = (x = knn(tree, y.data, 2)[1]; [ind[1] for ind in x])
        for (i,j) ∈ enumerate(nind)
            δ = norm(y[i]-y[j], 2)
            # If y[i] and y[j] are still identical, choose the next nearest neighbor
            # as in Cao's algorithm (not suggested by Kennel, but still advisable)
            if δ == 0.0
                j = knn(tree, y[i], 3, true)[1][end]
                δ = norm(y[i]-y[j])
            end
            δ1 = _increase_distance(δ,s,i,j,γ,τ,Euclidean())
            cond_1 = ((δ1/δ)^2 - 1 > rtol2) # equation (4) of Kennel
            cond_2 = (δ1/Ra > atol)         # equation (5) of Kennel
            if cond_1 | cond_2
                nfnn[k] += 1
            end
        end
    end
    return nfnn
end

"""
    f1nn(s::AbstractVector, τ:Int, γs = 1:5, metric = Euclidean())

Calculate the ratio of "false first nearest neighbors" (FFNN) of the datasets created
from `s` with a sequence of `τ`-delayed temporal neighbors.

## Description
Given a dataset made by embedding `s` with `γ` temporal neighbors and delay `τ`,
the "first nearest neighbors" (FFNN) are the pairs of points that are nearest to
each other at dimension `γ` that cease to be nearest neighbors at dimension
`γ+1` [1].

The returned value is a vector with the ratio between the number of FFNN and
the number of points in the dataset for each `γ ∈ γs`. The optimal value for `γ`
is found at the point where this ratio approaches zero.

See also: [`estimate_dimension`](@ref), [`afnn`](@ref), [`fnn`](@ref).

## References

[1] : Anna Krakovská *et al.*, "Use of false nearest neighbours for selecting
variables and embedding parameters for state space reconstruction", *J Complex
Sys* 932750 (2015), DOI: 10.1155/2015/932750
"""
function f1nn(s::AbstractVector, τ::Int, γs = 1:5, metric = Euclidean())
    f1nn_ratio = zeros(length(γs))
    γ_prev = 0 # to recall what γ has been analyzed before
    Rγ = reconstruct(s[1:end-τ],γs[1],τ) # this is for the first iteration
    for (i, γ) ∈ enumerate(γs)
        if i>1 && γ!=γ_prev+1
            # Re-calculate the series with γ delayed dims if γ does not follow
            # the dimension of the previous iteration
            Rγ = reconstruct(s[1:end-τ],γ,τ)
        end
        (nf1nn, Rγ) = _compare_first_nn(s,γ,τ,Rγ,metric)
        f1nn_ratio[i] = nf1nn/length(Rγ)
        # Trim Rγ for the next iteration
        Rγ = Rγ[1:end-τ,:]
        γ_prev = γ
    end
    return f1nn_ratio
end

function _compare_first_nn(s, γ::Int, τ::Int, Rγ::Dataset{D,T}, metric) where {D,T}
    # This function compares the first nearest neighbors of `s`
    # embedded with Dimensions `γ` and `γ+1` (the former given as input)
    tree = KDTree(Rγ,metric)
    Rγ1 = reconstruct(s,γ+1,τ)
    tree1 = KDTree(Rγ1,metric)
    nf1nn = 0
    # For each point `i`, the fnn of `Rγ` is `j`, and the fnn of `Rγ1` is `k`
    nind = (x = knn(tree, Rγ.data, 2)[1]; [ind[1] for ind in x])
    @inbounds for  (i,j) ∈ enumerate(nind)
        k = knn(tree1, Rγ1.data[i], 2, true)[1][end]
        if j != k
            nf1nn += 1
        end
    end
    # `Rγ1` is returned to re-use it if necessary
    return (nf1nn, Rγ1)
end

"""
    stochastic_indicator(s::AbstractVector, τ:Int, γs = 1:4) -> E₂s

Compute an estimator for apparent randomness in a reconstruction with `γs` temporal
neighbors.

## Description
Given the scalar timeseries `s` and the embedding delay `τ` compute the
values of `E₂` for each `γ ∈ γs`, according to Cao's Method (eq. 5 of [1]).

Use this function to confirm that the
input signal is not random and validate the results of [`estimate_dimension`](@ref).
In the case of random signals, it should be `E₂ ≈ 1 ∀ γ`.

## References

[1] : Liangyue Cao, [Physica D, pp. 43-50 (1997)](https://www.sciencedirect.com/science/article/pii/S0167278997001188?via%3Dihub)
"""
function stochastic_indicator(s::AbstractVector{T},τ, γs=1:4) where T # E2, equation (5)
    #This function tries to tell the difference between deterministic
    #and stochastic signals
    #Calculate E* for Dimension γ+1
    E2s = Float64[]
    for γ ∈ γs
        Rγ1 = reconstruct(s,γ+1,τ)
        tree1 = KDTree(Rγ1[1:end-1-τ])
        method = FixedMassNeighborhood(2)

        Es1 = 0.
        nind = (x = neighborhood(Rγ1[1:end-τ], tree1, method); [ind[1] for ind in x])
        for  (i,j) ∈ enumerate(nind)
            Es1 += abs(Rγ1[i+τ][end] - Rγ1[j+τ][end]) / length(Rγ1)
        end

        #Calculate E* for Dimension γ
        Rγ = reconstruct(s,γ,τ)
        tree2 = KDTree(Rγ[1:end-1-τ])
        Es2 = 0.
        nind = (x = neighborhood(Rγ[1:end-τ], tree2, method); [ind[1] for ind in x])
        for  (i,j) ∈ enumerate(nind)
            Es2 += abs(Rγ[i+τ][end] - Rγ[j+τ][end]) / length(Rγ)
        end
        push!(E2s, Es1/Es2)
    end
    return E2s
end


text:
###################################################################################################
# IO FUNCTIONS
###################################################################################################
"""
    `print_log(MESSAGE)`

    Function that prints MESSAGE to stderr.

    # Examples
    ```julia-repl
    julia> CpelTdm.print_log("Hello")
    [2020-03-30 16:24:18]: Hello
    ```
"""
function print_log(mess::String)::Nothing

    # Right format for date
    date = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")
    println(stderr,"[$(date)]: "  * mess)
    flush(stderr)

    # Return
    return nothing

end
"""
    `sleep_exit(MESSAGE)`

    Function that prints MESSAGE, waits 5 seconds and exits without error.

    # Examples
    ```julia-repl
    julia> CpelTdm.sleep_exit("Bye")
    [2020-03-30 16:24:18]: Bye
    ```
"""
function sleep_exit(mess::String)::Nothing

    # Print message
    print_log(mess)
    sleep(10)
    exit(0)

    # Return
    return nothing

end
"""
    `get_paths(BAMS1,BAMS2,FASTA,BED,OUTDIR,OUTPREFIX)`

    Function that returns output file names.

    # Examples
    ```julia-repl
    julia> out_mml_paths,out_nme_paths,out_diff_paths = CpelTdm.get_paths(bams1,bams2,fasta,bed,outdir,outprefix)
    ```
"""
function get_paths(bams1::Vector{String},bams2::Vector{String},fasta::String,bed::String,outdir::String,
                   outprefix::String)::NTuple{3,Vector{String}}

    # Check BED file exists
    if !isfile(bed)
        sleep_exit("BED file was not found ...")
    end

    # Check FASTA index file exists
    if !all([isfile(fasta),isfile(fasta*".fai")])
        sleep_exit("Fasta reference and/or index were not found ...")
    end

    # Check BAM index files exists
    ind_miss = false
    for bam in vcat(bams1,bams2)
        if !all([isfile(bam),isfile(bam*".bai")])
            ind_miss |= true
            print_log("BAM and/or its index file missing for")
            print_log("$(bam)")
        end
    end
    ind_miss && sleep_exit("Index files missing ...")

    # Individual bedGraph output files
    prefixes = [String(split(basename(bam),".")[1]) for bam in bams1]
    append!(prefixes,[String(split(basename(bam),".")[1]) for bam in bams2])
    out_mml_paths = ["$(outdir)/$(prefix)_mml.bedGraph" for prefix in prefixes]
    out_nme_paths = ["$(outdir)/$(prefix)_nme.bedGraph" for prefix in prefixes]

    # Check for existance of at least an output files
    if any(isfile.(vcat(out_mml_paths,out_nme_paths)))
        sleep_exit("At least an MML or NME file already exists...")
    end

    # Differentiial analysis bedGraph output files
    out_diff_paths = "$(outdir)/$(outprefix)_" .* ["tmml","tnme","tpdm"] .* "_diff_analysis.bedGraph"

    # Check for existance of at least an output files
    if any(isfile.(out_diff_paths))
        sleep_exit("At least a differential output file already exists...")
    end

    # Return paths
    return out_mml_paths,out_nme_paths,out_diff_paths

end
"""
    `get_bed_out_filename(bed)`

    Function that returns name of output BED file.

    # Examples
    ```julia-repl
    julia> CpelTdm.get_bed_out_filename(bed)
    ```
"""
function get_bed_out_filename(bed::String)::String

    # Return name
    return dirname(bed) * "/" * split(basename(bed),".")[1] * "_cpel.bed"

end
"""
    `write_mml_out(OUT_PMAP,MML_PATHS)`

    Function that writes MML in corresponding bedGraph files.

    # Examples
    ```julia-repl
    julia> CpelTdm.write_mml_out(out_pmap,mml_paths)
    ```
"""
function write_mml_out(out_pmap::Vector{RoiData},mml_paths::Vector{String})::Nothing

    # Open streams
    ios = [open(path,"a") for path in mml_paths]

    # Loop over ROIs
    @inbounds for roi in out_pmap

        # Get data from ROI
        mmls = vcat(roi.mmls1,roi.mmls2)

        # Write
        @inbounds for i=1:length(ios)
            write(ios[i],"$(roi.chr)\t$(roi.chrst-1)\t$(roi.chrend)\t$(mmls[i])\n")
        end

    end

    # Close all streams
    close.(ios)

    # Return nothing
    return nothing

end
"""
    `write_nme_out(OUT_PMAP,NME_PATHS)`

    Function that writes NME in corresponding bedGraph files.

    # Examples
    ```julia-repl
    julia> CpelTdm.write_nme_out(out_pmap,nme_paths)
    ```
"""
function write_nme_out(out_pmap::Vector{RoiData},nme_paths::Vector{String})::Nothing

    # Open streams
    ios = [open(path,"a") for path in nme_paths]

    # Loop over ROIs
    @inbounds for roi in out_pmap

        # Get data from ROI
        nmes = vcat(roi.nmes1,roi.nmes2)

        # Write
        @inbounds for i=1:length(ios)
            write(ios[i],"$(roi.chr)\t$(roi.chrst-1)\t$(roi.chrend)\t$(nmes[i])\n")
        end

    end

    # Close all streams
    close.(ios)

    # Return nothing
    return nothing

end
"""
    `write_diff_out(OUT_PMAP,DIFF_PATHS)`

    Function that writes differential output in corresponding bedGraph files.

    # Examples
    ```julia-repl
    julia> CpelTdm.write_diff_out(out_pmap,diff_paths)
    ```
"""
function write_diff_out(out_pmap::Vector{RoiData},diff_paths::Vector{String})::Nothing

    # Open streams
    ios = [open(path,"a") for path in diff_paths]

    # Loop over ROIs
    @inbounds for roi in out_pmap

        # Get data from ROI
        tmml,pmml = roi.mml_test
        tnme,pnme = roi.nme_test
        tpdm,ppdm = roi.pdm_test

        # Write
        write(ios[1],"$(roi.chr)\t$(roi.chrst-1)\t$(roi.chrend)\t$(tmml)\t$(pmml)\n")
        write(ios[2],"$(roi.chr)\t$(roi.chrst-1)\t$(roi.chrend)\t$(tnme)\t$(pnme)\n")
        write(ios[3],"$(roi.chr)\t$(roi.chrst-1)\t$(roi.chrend)\t$(tpdm)\t$(ppdm)\n")

    end

    # Close all streams
    close.(ios)

    # Return nothing
    return nothing

end
"""
    `write_output(OUT_PMAP,OUT_PATHS)`

    Function that write output of pmap into respective files.

    # Examples
    ```julia-repl
    julia> CpelTdm.write_output(out_pmap,out_paths)
    ```
"""
function write_output(out_pmap::Vector{RoiData},out_paths::NTuple{3,Vector{String}})::Nothing

    # Assign paths
    mml_paths = out_paths[1]
    nme_paths = out_paths[2]
    diff_paths = out_paths[3]

    # Write MML
    write_mml_out(out_pmap,mml_paths)

    # Write NMEs
    write_nme_out(out_pmap,nme_paths)

    # Write differential analysis
    write_diff_out(out_pmap,diff_paths)

    # Return nothing
    return nothing

end


text:
module ClimateBase
export NCDataset

include("core/coredefs.jl")
include("core/prettyprint.jl")
include("core/aggregation.jl")
include("core/interpolation.jl")

include("physical_dimensions/spatial.jl")
include("physical_dimensions/spatial_equalarea.jl")
include("physical_dimensions/temporal.jl")

include("io/vector2range.jl")
include("io/netcdf.jl")

# All following will be moved to ClimateTools.jl once its updated
include("climate/solar.jl")
include("tsa/continuation.jl")
include("tsa/decomposition.jl")

include("exports.jl")

using Requires
function __init__()
    @require GeoMakie="db073c08-6b98-4ee5-b6a4-5efafb3259c6" begin
        include("plotting/geomakie.jl")
    end
end

end # module


text:
@testset "Monomial selection for certificate" begin
    @polyvar x y
    @test_throws ErrorException getmonomialsforcertificate([x*y, y^2], :Sparse)
    @test getmonomialsforcertificate([x*y, y^2]) == [y]
end

@testset "Random SOS should be SOS with $(factory.constructor)" for factory in sdp_factories
    @polyvar x y
    x = [1, x, y, x^2, y^2, x*y]
    @test_throws ArgumentError randsos(x, monotype=:Unknown)
    for i in 1:10
        for monotype in [:Classic, :Gram]
            p = randsos(x, monotype=monotype)

            m = SOSModel(factory)

            @constraint m p >= 0

            JuMP.optimize!(m)

            @test JuMP.primal_status(m) == MOI.FeasiblePoint
        end
    end
end


text:
using Distributions:Random
using AdaptiveABC
using SequentialLOB
using Random
using Distributions
using Optim
using NLopt
using JLD2
using Plots
using DelimitedFiles

# Set Simulator Parameters and Seed
parameter_names = ["D", "σ", "nu", "μ", "α"]
plotting_parameter_names = ["D", "σ", "ν", "μ", "α"]
n_replications = 10
Random.seed!(1717617000301)
# Constants
num_paths = 1; M = 400 ; T = 2299 ; p₀ = 238.75 ; L = 200 ; λ = 1.0

# Read in Observed Market Data

observed_price_path = round.(readdlm("$(dirname(@__FILE__))/Original_Price_Bars_2300.csv", ',', Float64, '\n')[:,1], digits=2)
observed_log_returns = diff(log.(observed_price_path[:,1]))
observed_summary_stats = get_summary_stats(observed_log_returns)

# plt = plot(
#     1:T + 1, observed_price_path[:,1], legend=false,
#     xlab="Time", ylab="Mid-Price",
# )

# Define Summary function used in AdaptiveABC. Takes in a vector of parameters and outputs a success boolean and the
# simulated summart statistics
function summary_fn(parameters, n_summary_stats, n_replications)
    try
        D, σ, nu, μ, α_slob = parameters
        model = SLOB(num_paths, T, p₀, M, L, D, σ, nu, α_slob, SourceTerm(λ, μ))
        summary_stats = Array{Float64,2}(undef, n_summary_stats, n_replications)
        for i in 1:n_replications
            sim_price_path = model()
            sim_obs = diff(log.(sim_price_path[:, 1]))
            summary_stats[:, i] = get_summary_stats(observed_log_returns, sim_obs)
        end
        return true, summary_stats
    catch e
        return false, zeros(n_summary_stats, n_replications)
    end
end

# Setup for Calibration of Parameters through AdaptiveABC
# Define Prior distribution using Prior struct
prior = Prior([
    Uniform(0.5, 3.0), # D
    Uniform(0.1, 3.0), # σ
    Uniform(0.05, 0.95), # nu
    Uniform(0.1, 3.0), # μ
    Uniform(50.0, 200.0) # α_slob
])

#####
# Calibration Technique: SLOB ABC-PMC BBWM
#####

# Distance Function to measure distance between observed summary stats and simulated summary stats
boot_weight_matrix = BlockBootstrapWeightMatrix(9201292, vcat([0], observed_log_returns), get_summary_stats, 100, 10_000)
weighted_bootstrap = WeightedBootstrap(vcat([0], observed_log_returns), get_summary_stats, boot_weight_matrix)

abc_input_wb = ABCInput(
    prior,
    parameter_names,
    observed_summary_stats,
    summary_fn,
    17,
    weighted_bootstrap
)
abc_pmc_out_wb = JLD2.load_object("$(dirname(@__FILE__))/abc_pmc_out_wb.jld2")
# abc_pmc_out_wb = ABC_PMC(
#     abc_input_wb,
#     100,
#     n_replications,
#     0.25,
#     10_000,
#     10;
#     parallel=true,
#     batch_size=250,
#     seed=18291002012543
# )
posterior_means = AdaptiveABC.parameter_means(abc_pmc_out_wb)

println()
println("ABC-PMC BBWM: D, σ, nu, μ, α_slob")
print(round.(posterior_means[:,end], digits=3))
println()

# JLD2.save_object("$(dirname(@__FILE__))/abc_pmc_out_wb.jld2", abc_pmc_out_wb)

# ABC-PMC BBWM: D, σ, nu, μ, α_slob
# [2.087, 0.421, 0.603, 1.202, 106.504]

param_inds = [1,2,3,4,5]
plt = plot(abc_pmc_out_wb, iteration_colours=cgrad(:blues, 5, categorical=true),
    iterations=[1, abc_pmc_out_wb.n_iterations ÷ 2, abc_pmc_out_wb.n_iterations],
    prior_dists=abc_input_wb.prior.distribution[param_inds],
    param_inds=param_inds, param_names=plotting_parameter_names[param_inds])
plt
savefig(plt, "$(dirname(@__FILE__))/abc_pmc_wb.pdf")

#####
# Calibration Technique: ABC-PMC MADWE
#####
abc_input_we = ABCInput(
    prior,
    parameter_names,
    observed_summary_stats,
    summary_fn,
    17,
    WeightedEuclidean(observed_log_returns, get_summary_stats, ones(17), "MAD")
)

abc_pmc_out_we = JLD2.load_object("$(dirname(@__FILE__))/abc_pmc_out_we.jld2")
# abc_pmc_out_we = ABC_PMC(
#     abc_input_we,
#     100,
#     n_replications,
#     0.15,
#     10_000,
#     10;
#     parallel=true,
#     batch_size=250,
#     seed=77261511
# )
posterior_means = AdaptiveABC.parameter_means(abc_pmc_out_we)
println("")
println("ABC-PMC MADWE: D, σ, nu, μ, α")
print(round.(posterior_means[:,end], digits=3))
println("")
# ABC-PMC MADWE: D, σ, nu, μ, α
# [1.661, 1.127, 0.667, 1.987, 127.368]

# JLD2.save_object("$(dirname(@__FILE__))/abc_pmc_out_we.jld2", abc_pmc_out_we)

param_inds = [1,2,3,4,5]
plt = plot(abc_pmc_out_we, iteration_colours=cgrad(:blues, 5, categorical=true),
    iterations=[1,abc_pmc_out_we.n_iterations],
    prior_dists=abc_input_we.prior.distribution[param_inds],
    param_inds=param_inds, param_names=plotting_parameter_names[param_inds])
savefig(plt, "$(dirname(@__FILE__))/abc_pmc_we.pdf")
plt

max_stars_repo_name: GantZA/CALOBMTMD

max_issues_repo_name: GantZA/CALOBMTMD

max_forks_repo_name: GantZA/CALOBMTMD


text:
@withkw immutable RateSynapseParameter
    lr::Float = 1e-3
end

@withkw type RateSynapse
    param::RateSynapseParameter = RateSynapseParameter()
    colptr::Vector{Int} # column pointer of sparse W
    I::Vector{Int}      # postsynaptic index of W
    W::Vector{Float}  # synaptic weight
    rI::Vector{Float} # postsynaptic rate
    rJ::Vector{Float} # presynaptic rate
    g::Vector{Float}  # postsynaptic conductance
    records::Dict = Dict()
end

function RateSynapse(pre, post; σ = 0.0, p = 0.0)
    w = σ / √(p * pre.N) * sprandn(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    rI, rJ, g = post.r, pre.r, post.g
    RateSynapse(;@symdict(colptr, I, W, rI, rJ, g)...)
end

@replace function forward!(c::RateSynapse, param::RateSynapseParameter)
    fill!(g, zero(eltype(g)))
    @inbounds for j in 1:(length(colptr) - 1)
        rJj = rJ[j]
        for s = colptr[j]:(colptr[j+1] - 1)
            g[I[s]] += W[s] * rJj
        end
    end
end

@replace function plasticity!(c::RateSynapse, param::RateSynapseParameter, dt::Float, t::Float)
    @inbounds for j in 1:(length(colptr) - 1)
        s_row = colptr[j]:(colptr[j+1] - 1)
        rIW = zero(Float)
        for s in s_row
            rIW += rI[I[s]] * W[s]
        end
        Δ = lr * (rJ[j] - rIW)
        for s in s_row
            W[s] += rI[I[s]] * Δ
        end
    end
end


text:
using QML
using Qt5QuickControls_jll
using Observables

const qml_file = "Observables.qml"
const input = Observables(1.0)
const output = Observables(0.0)

on(output) do x
    println("Output changed to ", x)
end

load(qml_file, input = input, output = output)
exec_async()

max_stars_repo_name: ErikQQY/Julie

max_issues_repo_name: ErikQQY/Julie

max_forks_repo_name: ErikQQY/Julie


text:
losses = Float64[10, 8, 9, 10, 11, 12, 12, 13, 14, 15, 16, 17, 16]

# codecov:
@test EarlyStopping._min(nothing, 5) == 5

@testset "Never" begin
    @test stopping_time(Never(), losses) == 0
end

@testset "NotANumber" begin
    @test stopping_time(NotANumber(), losses) == 0
    N = 5
    losses2 = fill(123.4, N)
    @test all(reverse(eachindex(losses2))) do j
        losses2[j] = NaN
        stopping_time(NotANumber(), losses2) == j
    end
    losses2 = Float64[1, 2, 3, 1, NaN, 3, 1, 2, 3]
    is_training = Bool[1, 1, 0, 1, 1, 0, 1, 1, 0]
    @test stopping_time(NotANumber(), losses2, is_training) == 2
    losses2 = Float64[1, 2, 3, 1, 2, NaN, 1, 2, 3]
    @test stopping_time(NotANumber(), losses2, is_training) == 2
    losses2 = Float64[1, 2, 3, 1, 2, 3, NaN, 2, 3]
    @test stopping_time(NotANumber(), losses2, is_training) == 3
    losses2 = Float64[1, 2, 3, 1, 2, 3, 1, 2, 3]
    @test stopping_time(NotANumber(), losses2, is_training) == 0
    @test_logs((:info, r"loss updates: 0"),
               (:info, r"state: true"),
               (:info, r"loss updates: 1"),
               (:info, r"state: true"),
               stopping_time(NotANumber(),
                             [NaN, 1],
                             [true, false],
                             verbosity=1))
end

struct SleepyIterator{T}
    iter::T
    t::Float64
end
SleepyIterator(iter; t=0.1) = SleepyIterator(iter, t)

Base.iterate(iter::SleepyIterator) = (sleep(iter.t); iterate(iter.iter))
Base.iterate(iter::SleepyIterator, state) =
    (sleep(iter.t); iterate(iter.iter, state))

@testset "TimeLimit" begin
    @test_throws ArgumentError TimeLimit(t=0)
    @test TimeLimit(1).t == Millisecond(3_600_000)
    @test TimeLimit(t=Day(2)).t == Millisecond(48*3_600_000)
    sleepy_losses = SleepyIterator(losses; t=0.1)
    @test stopping_time(TimeLimit(t=Millisecond(600)), sleepy_losses) == 7
    # codecov:
    @test EarlyStopping.update_training(TimeLimit(), 42.0) <= now()
end

@testset "GL" begin
    @test !EarlyStopping.needs_in_and_out_of_sample(GL())

    # constructor:
    @test_throws ArgumentError GL(alpha=0)
    @test GL(alpha=1).alpha === 1.0

    # stopping times:
    n = @test_logs((:info, r"loss updates: 1"),
                   (:info, r"state: \(loss = 10.0, min_loss = 10.0\)"),
                   (:info, r"loss updates: 2"),
                   (:info, r"state: \(loss = 8.0, min_loss = 8.0\)"),
                   (:info, r"loss updates: 3"),
                   (:info, r"state: \(loss = 9.0, min_loss = 8.0\)"),
                   stopping_time(GL(alpha=12), losses, verbosity=1))
    @test n == 3
    @test stopping_time(GL(alpha=20), losses) == 4
    @test stopping_time(GL(alpha=40), losses) == 6
    @test stopping_time(GL(alpha=90), losses) == 11
    @test stopping_time(GL(alpha=110), losses) == 12
    @test stopping_time(GL(alpha=1000), losses) == 0
end

@testset "PQ" begin
    v = [:c, :b, :a]
    @test EarlyStopping.prepend(v, :d, 4) == [:d, :c, :b, :a]
    @test EarlyStopping.prepend(v, :d, 3) == [:d, :c, :b]
    @test EarlyStopping.progress(Float64[2, 3, 6, 7]) ≈ 1000*(4.5 - 2)/2
    @test_throws ArgumentError PQ(alpha=0)
    @test_throws ArgumentError PQ(k=-Inf)
    @test_throws ArgumentError PQ(k=1.2)
    @test_throws ArgumentError PQ(k=1)

    c = PQ(alpha=10, k=2)
    @test EarlyStopping.needs_in_and_out_of_sample(c)

    # first update must be training:
    @test_throws Exception EarlyStopping.update(c, 1.0)

    state = EarlyStopping.update_training(c, 10.0)
    # at least two training updates before out-of-sample update:
    @test_throws Exception EarlyStopping.update(c, state, 10.0)

    state = EarlyStopping.update_training(c, 10.0, state)
    state = EarlyStopping.update(c, 10.0, state)
    @test EarlyStopping.done(c, state) # progress = 0

    state = EarlyStopping.update_training(c, 10.0, state)
    # can't be done if last update was a training update:
    @test !EarlyStopping.done(c, state)

    #                 k=2                progress GL    PQ    t
    losses2 = [9.5, 9.3, 10,            # 10.8     0     0     1
              9.3, 9.1, 8.9, 8,        # 11.2     0     0     2
              8.3, 8.4, 9,             # 6.02     12.5  2.08  3
              9.9, 9.5, 10,            # 21.2     25.0  1.18  4
              10.6, 10.4, 11,          # 9.61     37.5  3.90  5
              8888, 11.8, 11.7, 12,    # 4.27     50.0  11.7  6
              11.6, 11.4, 12,          # 8.77     50.0  5.70  7
              12.2, 12.1, 13,          # 4.1      62.5  15.2  8
              14.5, 14.1, 14,          # 14.2     75    5.28  9
              13.9, 13.7, 15,          # 7.30     87.5  11.9  10
              12.5, 12.3, 16,          # 8.13     100   12.3  11
              11.2, 11.0, 17,          # 9.09     112.5 12.4  12
              10.5, 10.3, 16]          # 9.71     100   10.3  13
    is_training = Bool[1, 1, 0,
                       1, 1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0,
                       1, 1, 0]

    @test stopping_time(PQ(alpha=2.0, k=2), losses2, is_training) == 3
    @test stopping_time(PQ(alpha=3.8, k=2), losses2, is_training) == 5
    @test stopping_time(PQ(alpha=11.6, k=2), losses2, is_training) == 6
    @test stopping_time(PQ(alpha=15.1, k=2), losses2, is_training) == 8
    @test stopping_time(PQ(alpha=15.3, k=2), losses2, is_training) == 0
end

@testset "Patience" begin
    @test_throws ArgumentError Patience(n=0)
    @test stopping_time(Patience(n=6), losses) == 0
    @test stopping_time(Patience(n=5), losses) == 12
    @test stopping_time(Patience(n=4), losses) == 6
    @test stopping_time(Patience(n=3), losses) == 5
    @test stopping_time(Patience(n=2), losses) == 4
    @test stopping_time(Patience(n=1), losses) == 3
end

true


text:
# Main Interface to BAMF

"""
    function bamf_roi()

run BAMF on single ROI and return chain

This is the main entry point to the BAMF algorithm

```
function bamf_roi(data::BAMFData, psf::MicroscopePSFs.PSF, prior_photons::Distributions.UnivariateDistribution, prior_background::Distributions.UnivariateDistribution;
    bndpixels::Real = 0.0,
    σ_xy::Real = 0.05,
    σ_split::Real = 0.5,
    σ_photons::Real = sqrt(mean(prior_photons)) / 10,
    σ_background::Real = sqrt(mean(prior_background)) / 10,
    burnin::Int=5000,
    iterations::Int=5000,
    jumpprobability = [.5, .1, .1, .1, .1, .1]
)
```

# Arguments
- `data::BAMFData`: sub region dataa
- `psf::MicroscopePSFs.PSF`: point spread function
- `prior_photons::Distributions.UnivariateDistribution`: prior probability on photons/emitter/frame
- `prior_background::Distributions.UnivariateDistribution`: prior probability on bg photons/pixel

# Optional Keyword Arguments
- `bndpixels::Real = 0.0`: boundary region. +1 means allow positoins 1 pixel larger than ROI.
- `σ_xy::Real = 0.05`: parameterizes MCMC jumps in x,y
- `σ_split::Real = 0.5`: parameterizes MCMC jumps in seperation in splitting
- `σ_photons::Real = sqrt(mean(prior_photons)) / 10`: parameterizes MCMC jumps in photons/emitter
- `σ_background::Real = sqrt(mean(prior_background)) / 10`: parameterizes MCMC jumps bg
- `burnin::Int=5000`:
- `iterations::Int=5000`:
- `jumpprobability = [.5, .1, .1, .1, .1, .1]`: [move, bg, add,remove,split,merge]. Will be normalized internally

`jumpprobability` can be used to fix number of emitters, bg, or emitter x,y,photons by setting jump probability to zero.

"""
function bamf_roi(data::BAMFData, psf::MicroscopePSFs.PSF, prior_photons::Distributions.UnivariateDistribution, prior_background::Distributions.UnivariateDistribution;
    bndpixels::Real = 0.0,
    σ_xy::Real = 0.05,
    σ_split::Real = 0.5,
    σ_photons::Real = min(sqrt(mean(prior_photons)) / 10 , std(prior_photons) / 10),
    σ_background::Real = sqrt(mean(prior_background)) / 10,
    burnin::Int=5000,
    iterations::Int=5000,
    jumpprobability = [.5, .1, .1, .1, .1, .1]
)
    currentdata=deepcopy(data)
    println(typeof(currentdata))
    testdata=deepcopy(data)

    myRJ = RJStruct(data, psf, prior_photons, prior_background, bndpixels,
        σ_xy, σ_split, σ_photons, σ_background, currentdata, testdata)

    ## setup the RJMCMC.jl model
    # Jumptypes are: move, bg, add, remove, split, merge
    njumptypes = 6
    jumpprobability = jumpprobability / sum(jumpprobability)

    # create an RJMCMC structure with all model info
    acceptfuns = [accept_move,  accept_bg,  accept_add,  accept_remove,  accept_split,  accept_merge] #array of functions
    propfuns = [propose_move,  propose_bg,  propose_add,  propose_remove,  propose_split,  propose_merge] #array of functions
    myRJMCMC = ReversibleJumpMCMC.RJMCMCStruct(burnin, iterations, njumptypes, jumpprobability, propfuns, acceptfuns)

    #create an intial state
    state1 =  calcintialstate(myRJ)

    ## run chain. This is the call to the main algorithm
    return ReversibleJumpMCMC.buildchain(myRJMCMC, myRJ, state1)

end


text:
@static if Base.VERSION >= v"1.7-"
    const write_manifest = Pkg.Types.write_manifest
else
    function write_manifest(io::IO, manifest::Dict)
        print(io, "# This file is machine-generated - editing it directly is not advised\n\n")
        TOML.print(io, manifest, sorted=true) do x
            x isa UUID || x isa SHA1 || x isa VersionNumber || pkgerror("unhandled type `$(typeof(x))`")
            return string(x)
        end
        return nothing
    end
end


text:
module MesocosmMetabolism

using DifferentialEquations
using Interpolations
using Turing
using DataFrames
using Plots
using Distributions
using LinearAlgebra
using OrdinaryDiffEq
using FillArrays

import Distributions.truncated
import Distributions.Normal
import Distributions.Cauchy
import Distributions.LogNormal
import OrdinaryDiffEq.AutoVern7
import OrdinaryDiffEq.Rodas4

include("processModel.jl")
include("initializeProcessModel.jl")
include("statisticalModel.jl")
include("simulateData.jl")

export mesocosm_metabolism!
export initialize_process_model
export statisticalModel
export simulate_data

end

max_stars_repo_name: elizabethjmohr/mesocosm-metabolism

max_issues_repo_name: elizabethjmohr/mesocosm-metabolism

max_forks_repo_name: elizabethjmohr/mesocosm-metabolism


text:
using Test
using SparseIR
using MathLink
sphericalbesselj_M(n, x) = weval(W"SphericalBesselJ"(W"n", W"x"); n, x)

@testset "_specfuncs.jl" begin
    @testset "domain" begin
        @test_throws DomainError SparseIR.sphericalbesselj(-2, 0.3)
        @test !isnan(SparseIR.sphericalbesselj(4, 1e20))
    end

    @testset "accuracy (with pn = $pn, px = $px)" for pn in 0:7, px in 0:25
        n = 2^pn
        x = float(2^px)
        @test SparseIR.sphericalbesselj(n, x) ≈ sphericalbesselj_M(n, x)
        @test SparseIR.sphericalbesselj(n + 1, x) ≈ sphericalbesselj_M(n + 1, x)
    end
end


text:
Rectangular() = Uniform(-0.5, +0.5)

function _string(kernel::Uniform)
    a = kernel.a
    b = kernel.b
    "Rectangular kernel (U[$a,$b])"
end

const LocationScaleDists = Union{Cosine,Epanechnikov,SymTriangularDist,Triweight}

_string(kernel::SymTriangularDist) = "Triangular kernel"
_string(kernel::LocationScaleDists) = "$(Base.typename(typeof(kernel))) kernel"

const SupportedKernels = Union{LocationScaleDists,Uniform}

export Uniform, Epanechnikov, SymTriangularDist

function setbandwidth(kernel::LocationScaleDists, h::Number)
    typeof(kernel)(zero(kernel.μ), h * kernel.σ)
end

function setbandwidth(kernel::Uniform, h::Number)
    a = kernel.a
    b = kernel.b

    if a + b != 0
        error("Use symmetric Uniform/Rectangular kernel")
    end
    Uniform(a * h, b * h)
end

weights(::Uniform, ZsR::RunningVariable) = uweights(length(ZsR))
weights(D::Distribution, ZsR::RunningVariable) = pdf.(D, ZsR)

struct EquivalentKernel{K,T}
    kernel::K #only order 1 for now
    kernel_mult::T
end

function EquivalentKernel(kernel::SupportedKernels)
    ub = support(kernel).ub
    E(f) = quadgk(f, 0, ub)[1]
    cov_x = [x->1 x->x; x->x x->x^2]
    kernel_mult = [1 0] * inv(map(f -> E(u -> f(u) * pdf(kernel, u)), cov_x))
    EquivalentKernel(kernel, kernel_mult)
end

Distributions.support(eq::EquivalentKernel) = Distributions.support(eq.kernel)

function Distributions.pdf(eq::EquivalentKernel, u)
    kernel_mult = eq.kernel_mult
    kernel = eq.kernel
    dot(kernel_mult, [1; u]) * pdf(kernel, u)
end

function kernel_moment(kernel, ::Val{j}) where {j}
    ub = support(kernel).ub
    quadgk(x -> x^j * pdf(kernel, x), 0, ub)[1]
end

kernel_moment(::SupportedKernels, ::Val{0}) = 0.5
kernel_moment(::SupportedKernels, ::Val{1}) = 0
kernel_moment(::SupportedKernels, ::Val{3}) = 0

function kernel_moment(eq::EquivalentKernel, ::Val{j}) where {j}
    ub = support(eq).ub
    E(f) = quadgk(f, 0, ub)[1]
    E(x -> x^j * pdf(eq, x))
end

function squared_kernel_moment(kernel, ::Val{j}) where {j}
    ub = support(kernel).ub
    quadgk(x -> x^j * pdf(kernel, x)^2, 0, ub)[1]
end

squared_kernel_moment(::SupportedKernels, ::Val{1}) = 0
squared_kernel_moment(::SupportedKernels, ::Val{3}) = 0

function squared_kernel_moment(eq::EquivalentKernel, ::Val{j}) where {j}
    ub = support(eq).ub
    E(f) = quadgk(f, 0, ub)[1]
    E(x -> x^j * pdf(eq, x)^2)
end


text:
@testset "TMJets algorithm (TMJets21b)" begin
    prob, tspan = vanderpol()

    # default algorithm for nonlinear systems
    sol = solve(prob, tspan=tspan)
    @test sol.alg isa TMJets

    # pass the algorithm explicitly
    sol = solve(prob, tspan=tspan, TMJets())
    @test sol.alg isa TMJets

    # pass options outside algorithm
    sol = solve(prob, T=0.2, maxsteps=2100)
    @test sol.alg isa TMJets && sol.alg.maxsteps == 2100
    sol = solve(prob, T=0.2, maxsteps=2100) # alias
    @test sol.alg isa TMJets && sol.alg.maxsteps == 2100
    sol = solve(prob, T=0.2, orderT=7, orderQ=1)
    @test sol.alg isa TMJets && sol.alg.orderT == 7 && sol.alg.orderQ == 1
    sol = solve(prob, T=0.2, abstol=1e-11)
    @test sol.alg isa TMJets && sol.alg.abstol == 1e-11
    sol = solve(prob, T=0.2, abstol=1e-11)
    @test sol.alg isa TMJets && sol.alg.abstol == 1e-11

    @test ReachabilityAnalysis.tspan(shift(sol,1.0))==ReachabilityAnalysis.tspan(sol)+1.0

    # split initial conditions
    X0, S = initial_state(prob), system(prob)
    X0s = split(X0, [2, 1]) # split along direction x
    sols = solve(IVP(S, X0s), T=0.1, threading=false) # no threading
    @test flowpipe(sols) isa MixedFlowpipe

    sols = solve(IVP(S, X0s), T=0.1, threading=true) # with threading (default)
    @test flowpipe(sols) isa MixedFlowpipe
end

@testset "TMJets algorithm (TMJets21b): linear IVPs" begin
    prob, dt = exponential_1d()
    sol = solve(prob, tspan=dt, TMJets())
    @test sol.alg isa TMJets

    # getter functions for a taylor model reach-set
    R = sol[1]
    @test domain(R) == tspan(R)
    @test diam(remainder(R)[1]) < 1e-9
    @test get_order(R) == [8]
    @test polynomial(R) isa Vector{Taylor1{TaylorN{Float64}}}
    @test expansion_point(R) ≈ [IA.Interval(0.0)]

    # test intersection with invariant
    prob, dt = exponential_1d(invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
    sol_inv = solve(prob, tspan=dt, TMJets())
    @test [0.3] ∈ overapproximate(sol_inv[end], Zonotope)
    m = length(sol_inv)
    # check that the following reach-set escapes the invariant
    @test [0.3] ∉ overapproximate(sol[m+1], Zonotope)

    # TODO test higher order system
#    prob, tspan = linear5D_homog()
#    sol = solve(prob, tspan=tspan, TMJets())
#    @test sol.alg isa TMJets

    # TODO test linear system with input
#    prob, tspan = linear5D()
#    sol = solve(prob, tspan=tspan, TMJets())
#    @test sol.alg isa TMJets
end

@testset "TMJets algorithm (TMJets20): linear IVPs" begin
    prob, dt = exponential_1d()
    sol = solve(prob, tspan=dt, TMJets20())
    @test sol.alg isa TMJets20

    # getter functions for a taylor model reach-set
    R = sol[1]
    @test domain(R) == tspan(R)
    @test diam(remainder(R)[1]) < 1e-13
    @test get_order(R) == [8]
    @test polynomial(R) isa Vector{Taylor1{TaylorN{Float64}}}
    @test expansion_point(R) ≈ [IA.Interval(0.0)]

    # test intersection with invariant
    prob, dt = exponential_1d(invariant=HalfSpace([-1.0], -0.3)) # x >= 0.3
    sol_inv = solve(prob, tspan=dt, TMJets20())
    @test [0.3] ∈ overapproximate(sol_inv[end], Zonotope)
    m = length(sol_inv)
    # check that the following reach-set escapes the invariant
    @test [0.3] ∉ overapproximate(sol[m+1], Zonotope)
end

@testset "1D Burgers equation (TMJets21b)" begin
    L0 = 1. # domain length
    U0 = 1. # Re = 20.
    x = range(-0.5*L0, 0.5*L0, length=4)
    # Initial velocity
    X0 = Singleton(-U0*sin.(2*π/L0*x))
    # IVP definition
    prob = @ivp(x' = burgers!(x), dim=4, x(0) ∈ X0)
    sol = solve(prob, tspan=(0.0, 1.0), alg=TMJets());
    @test dim(sol) == 4
end

#=
alg = TMJets(abstol=1e-10, orderT=10, orderQ=2)

# reach mode
sol = solve(P, T=7.0, alg)
@test set(sol[1]) isa Hyperrectangle # check default set representation
=#


text:
"""
    calculateSWI(data, options)

Returns the calculated SWI using 'data' and 'options'.
"""
function calculateSWI(data, options=Options())
    if !isnothing(options.writesteps) mkpath(options.writesteps) end
    getswimag(data, options) .* getswiphase(data, options)
end

createMIP(S, d=7) = createIntensityProjection(S, minimum, d)

function createIntensityProjection(S, func, d=7)
    [func(S[x,y,z:z+d-1]) for x in 1:size(S,1), y in 1:size(S,2), z in 1:size(S,3)-d+1]
end


text:
using DelimitedFiles

function file_to_array(infile)
	 array = readdlm(infile, ',', Int, '\n')
	 return array
end
function number_pair_search(report_array::Array{Int})
	 found_values = []
	 for reportnum::Int in report_array
	     complement::Int = 2020 - reportnum
	     found = findall(isequal(complement), report_array)
	     if size(found)[1] == 1
	     	found_values = [reportnum, complement]
		break
	     end
	 end
	 if size(found_values)[1] == 2
	    return found_values
	 else
	    return "Could not find a pair of values that sum to 2020.\n"
	 end
end

function number_triplet_search(report_array::Array{Int})
	 found_values = []
	 for i in 1:size(report_array,1)
	     complement::Int = 2020 - report_array[i]
	     for j in 2:size(report_array,1)
	         found = findall(isequal(complement-report_array[j]), report_array)
	     	 if size(found)[1] == 1
	     	    found_values = [report_array[i], report_array[j], complement-report_array[j]]
		    break
		 end
	     end
	 end
	 if size(found_values)[1] == 3
	    return found_values
	 else
	    return "Could not find a triplet of values that sum to 2020.\n"
	 end
end

function calculate_product(values)
	 product = 1
	 for i in 1:size(values,1)
	     product *= values[i]
	 end

	 return product
end

function main()
    val_array = file_to_array("input/01.txt")
    pair_values = number_pair_search(val_array)
    @assert size(pair_values)[1] == 2
    pair_product = calculate_product(pair_values)
    println("Pair values product: $pair_product")
    triplet_values = number_triplet_search(val_array)
    @assert size(triplet_values)[1] == 3
    triple_product = calculate_product(triplet_values)
    println("Triplet values product: $triple_product")
end

main()


text:
using StructWalk
using Test

@testset "StructWalk.jl" begin
    @test postwalk(x->x isa Expr ? eval(x) : x, :(2*3+5)) == 2*3+5
    @test postwalk(x->x isa Expr ? eval(x) : x==:+ ? :* : x==:* ? :+ : x, :(2*3+5)) == (2+3)*5
    @test postwalk(x-> x isa Int ? Float32(x) : x, Dict(:a=>3, :b=>5)) == Dict(:a=>3f0, :b=>5f0)

    struct C
        x::Int
        y::Int
        z::Int
    end
    x = (a=3, b=5, c=(x=4, y=4, z=6), d=((((5,),),),))
    @test postwalk(x) do x
        x isa NamedTuple{(:x, :y, :z)} && return C(x...)
        x isa Tuple && return first(x)
        x isa Integer && return float(x)
        return x
    end == (a=3.0, b=5.0, c=C(4., 4., 6.), d=5.)

    f(x) = x isa Integer ? StructWalk.LeafNode(x // 2) : x isa Tuple ? =>(x .+ 1...) : x
    @test prewalk(f, (3, 5)) == (2//1 => 3//1)
end


text:
using OpenEEW
using Test

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