text:
using BioProfiling
using Test
using DataFrames
using Statistics
using StatsBase
using LinearAlgebra: I
using Random
using RCall
using Distributed
using ParallelDataTransfer

@testset "freqtable" begin
	d = DataFrame(Any[0.05131 0.32830 "Exp1"; 0.83296 0.97647 "Exp1"; 0.66463 0.66939 "Exp2";
	                  0.30651 0.58938 "Exp2"; 0.71313 0.18477 "Exp2"; 0.81810 0.16309 "Exp2";
	                  0.05657 0.06012 "Exp1"; 0.02205 0.17055 "Exp2"; 0.49819 0.91871 "Exp1";
	                  0.90857 0.18794 "Exp2"; 0.12327 0.00619 "Exp2"; 0.34146 0.62640 "Exp1"])
	rename!(d, [:Ft1, :Ft2, :Experiment])

	e1 = Experiment(d, description = "Test Experiment")

	d1 = "Reject cells with high Ft2 values"
	f1 = Filter(0.9, :Ft2, compare = isless,
	            description = d1)

	@test freqtable(e1, f1) == [2,10]
	filter_entries!(e1, f1)

	d2 = "Select only a single experiment"
    f2 = Filter("Exp1", :Experiment, description = d2)

    @test freqtable(e1, f2) == [7,3]

    @test freqtable(e1, :Experiment) == [3,7]
end

@testset "decorrelate" begin
    X = DataFrame([[1,2,3],[3,2,1],[0,1,2],[1,0,1]])
    @test decorrelate(X) == [1,4]
    @test decorrelate(X, ordercol = [3,1,2]) == [3]
end

@testset "hellinger" begin
	Random.seed!(777)
	# The covariance matrix s is expected to be
    # a symmetric positive definite matrix
    c = rand(5)
    s = rand(5,5)
    s = s*s' + I
    # Distance of a distribution to itself should be zero, up to
    # machine precision. Note: approximation includes relative term
    # so isapprox(1e-200, 0) = false
	@test hellinger(c, s, c, s) + 1 ≈ 1
	function testpositivedefinite()
		c = rand(5)
	    s = rand(5,5)
	    s = s*s' + I
	    return(hellinger(c, s, c .+1, s) > 0)
	end
	@test all([testpositivedefinite() for x in 1:20])
end

@testset "logtransform" begin
    @test logtransform(1) == 0
    @test_throws MethodError logtransform("str")
    @test round.(logtransform(1:5), digits=2) == [0.0, 0.69, 1.1, 1.39, 1.61]
end

@testset "mahalanobis" begin
	Random.seed!(777)
    c = rand(5)
    s = rand(5,5)
	@test mahalanobis(c, c, s) == 0
	function testpositivedefinite()
		c = rand(5)
		# The covariance matrix s is expected to be
	    # a symmetric positive definite matrix
	    s = rand(5,5)
	    s = s*s' + I
	    return(mahalanobis(c, c .+1, s) > 0)
	end
	@test all([testpositivedefinite() for x in 1:20])
end

@testset "normtransform" begin
    x = 1:5
	@test round.(normtransform(x,x), digits = 2) == [-1.35, -0.67, 0.0, 0.67, 1.35]
    @test round.(normtransform(x,x[1:3]), digits = 2)  == [-0.67, 0.0, 0.67, 1.35, 2.02]
    @test_throws MethodError normtransform("str",x)
    @test_throws MethodError normtransform(x,"str")
end

@testset "Experiment" begin
	# Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])
	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)
	@test e1.description == "No description provided"
	@test e1.selected_entries == 1:12

	# Additional checks that could be performed:
	# Test print format
end

@testset "Filter" begin
	d1 = "Select only a single experiment"
    f1 = Filter("Exp1", :Experiment, description = d1)
	d2 = "Reject cells in high density regions"
	f2 = Filter(0.2, :Intensity_MedianIntensity_NeurDensity, compare = isless,
	            description = d2)
	@test f1.description == d1
	@test f2.description == d2

	# Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])
	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)

	filter_entries!(e1, f1)
	@test e1.selected_entries == [1,2,7,9,12]

	filter_entries!(e1, f2)
	@test e1.selected_entries == [7]

	e2 = Experiment(d)
	filter_entries!(e2, [f1,f2])
	@test e2.selected_entries == [7]

	e3 = Experiment(d)
	f3 = Filter(0.8, :Ft1, compare = >, description = "Large feature 1")
	cf1 = CombinationFilter(f1,f2,intersect)
	cf2 = CombinationFilter(cf1,f3,union)

	@test filter_entries(e3, cf1) == [7]
	@test filter_entries(e3, cf2) == [2,6,7,10]

	# Additional checks that could be performed:
	# Filter.compare::Function -> Make sure it takes 2 arguments and return 1?
	# CombinationFilter.operator::Function -> Make sure it takes 2 lists and return 1?

	filter_entries!(e3, f2)

	# Add missing values
	e3.data.Ft1 = Array{Union{Missing, Float64},1}(e3.data.Ft1)
	e3.data.Experiment = Array{Union{Missing, String},1}(e3.data.Experiment)
	e3.data.Ft1[4:6] .= missing
	e3.data.Experiment[[5,10,12]] .= missing

	mf1 = MissingFilter(:Ft1)
	mf2 = MissingFilter()

	@test filter_entries(e3, mf1) == [7,8,10,11]
	@test filter_entries(e3, mf2) == [7,8,11]
end

@testset "Selector" begin
	# Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])

	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)

	s1 = Selector(x -> eltype(x) <: Number)
	s2 = Selector(x -> mean(x) > 0.5, subset = x -> x.Experiment .== "Exp1",
				  description = "High mean for Exp1")

	@test s1.subset === nothing
	@test s2.description == "High mean for Exp1"

	@test select_features(e1, s1) == [1,2]
	select_features!(e1, [s1, s2])
	@test e1.selected_features == [2]

	strToRemove = ["MedianIntensity", "MorePatterns"]
	s3 = NameSelector(x -> !any(occursin.(strToRemove, String(x))))
	e2 = Experiment(d)
	@test select_features(e2, s3) == [1,3]

	select_features!(e2, [s1, s2, s3])
	@test length(e2.selected_features) == 0

	e3 = Experiment(d)
	s4 = deepcopy(s1)
	# Inverse function: keeps textual features
	s4.summarize = !s4.summarize
	s5 = deepcopy(s3)
	# Inverse function: keeps features including "MedianIntensity"
	s5.summarize = !s5.summarize

	cs1 = CombinationSelector(s4,s5,union)
	@test select_features(e3, cs1) == [2,3]
end

@testset "filter!" begin
    # Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])
	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)

	# Create selectors
	s1 = Selector(x -> !(eltype(x) <: Number),
				  description = "Keep textual features")
	keptPatterns = ["MedianIntensity", "MorePatterns"]
	s2 = NameSelector(x -> any(occursin.(keptPatterns, String(x))),
					  description = "Keep features including chosen patterns")
	cs1 = CombinationSelector(s1,s2,union)

	# Create filters
	d1 = "Select only a single experiment"
    f1 = Filter("Exp1", :Experiment, description = d1)
	d2 = "Reject cells in high density regions"
	f2 = Filter(0.2, :Intensity_MedianIntensity_NeurDensity, compare = isless,
	            description = d2)

	filter!(e1,[cs1,f1,f2])
	@test length(e1.selected_features) == 2
	@test e1.selected_entries == [7]
end

@testset "diagnostic" begin
    # Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])

	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)

	d1 = "Select only a single experiment"
	f1 = Filter("Exp1", :Experiment, description = d1)
	d2 = "Reject cells in high density regions"
	f2 = Filter(0.2, :Intensity_MedianIntensity_NeurDensity, compare = isless,
	            description = d2)

	cf1 = CombinationFilter(f1,f2,intersect)

	@test diagnostic(e1, cf1, features = [:Ft1]) == DataFrame(Ft1 = 0.0565727)

	@test diagnostic_path(e1, cf1, :Ft1) == [0.0565727]
	@test diagnostic_path(e1, cf1, :Experiment, rgx = [r".*" => s"example.png"]) == ["example.png"]

	@test diagnostic_images(e1, cf1, :Experiment, rgx = [r".*" => s"example.png"], saveimages = false)
	# Additional checks that could be performed:
	# Centers of diagnostic_path
	# getColorImage [internal]
	# colimgifrgb [internal]
	# rgb parameters of diagnostic_images
	# output of diagnostic_images
end

@testset "interpret" begin
	# Define example dataset
	d = DataFrame([[0,2,4,6,8,21],[0,1,2,3,4,0],[0,3,6,9,12,3],'A':'F'])
	rename!(d, [Symbol.("Ft".*string.(1:3))..., :Class])

	f1 = Filter('F', :Class, compare = !=, description = "Exclude row")
	s1 = NameSelector(x -> occursin("Ft", String(x)), "Keep numerical features")

	e = Experiment(d)
	select!(e, [f1,s1])

	@test most_variable_features(e) == ["Ft3", "Ft1", "Ft2"]

	# Test equivalence to sorting by mad
	top1ft = most_variable_features(e, top = 1)
	decorrelate_by_mad!(e)
	@test names(e.data)[e.selected_features] == top1ft

    d = DataFrame(rand(120,3))
	rename!(d, Symbol.("Ft".*string.(1:3)))
	d.Class = repeat(["Ref", "A", "B"], 40)

	# Add specific differences for each class
	d[d.Class .== "A",:Ft2] .+= 2;
	d[d.Class .== "B",:Ft3] .+= 2;

	f1 = Filter(0.9, :Ft1, compare = <=, description = "Exclude 10% of entries")
	s1 = NameSelector(x -> occursin("Ft", String(x)), "Keep numerical features")

	e = Experiment(d)
	select!(e, [f1,s1])

	filt_Ref = Filter("Ref", :Class)
	filt_A = Filter("A", :Class)
	filt_B = Filter("B", :Class)

	@test characteristic_features(e, filt_Ref, filt_A, top = 1) == ["Ft2"]
	@test characteristic_features(e, filt_Ref, filt_B, top = 1) == ["Ft3"]
	@test characteristic_features(e, filt_A, filt_B)[3] == "Ft1"
	@test characteristic_features(e, filt_A, filt_B) == characteristic_features(e, filt_B, filt_A)

	correlated_ref = 3 .* getdata(e).Ft2;
	@test most_correlated(e, correlated_ref, top = 1)[1] == "Ft2"
	@test most_correlated(e, :Ft1)[1] == "Ft1"
	e.data.Ft2 .-= 2 .* e.data.Ft3
	@test most_correlated(e, :Ft3) == ["Ft3", "Ft2", "Ft1"]
end

@testset "negation" begin
	# Define example dataset
	d = DataFrame(Any[0.0513198 0.328301 "Exp1"; 0.832986 0.976474 "Exp1"; 0.664634 0.669392 "Exp2";
	                  0.306519 0.58938 "Exp2"; 0.71313 0.184778 "Exp2"; 0.818107 0.163095 "Exp2";
	                  0.0565727 0.0601279 "Exp1"; 0.022015 0.170559 "Exp2"; 0.498196 0.918719 "Exp1";
	                  0.908576 0.187947 "Exp2"; 0.123237 0.00619995 "Exp2"; 0.341462 0.626406 "Exp1"])

	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity, :Experiment])
	d.Ft1 = convert.(Float64, d.Ft1)
	d.Intensity_MedianIntensity_NeurDensity = convert.(Float64, d.Intensity_MedianIntensity_NeurDensity)
	d.Experiment = convert.(String, d.Experiment)

	e1 = Experiment(d)

	# Test negation of simple filter
	d1 = "Select only a single experiment"
	f1 = Filter("Exp1", :Experiment, description = d1)
	nf1 = negation(f1)

	@test nf1.description == "Do not "*f1.description
	@test all(e1.data.Experiment[filter_entries(e1, f1)] .== "Exp1")
	@test all(e1.data.Experiment[filter_entries(e1, nf1)] .== "Exp2")

	# Test negation of simple selector
	s1 = Selector(x -> eltype(x) <: Number, description = "Keep numeric features")
	ns1 = negation(s1)

	@test ns1.description == "Do not "*s1.description

	ft_ns1 = select_features(e1, ns1)
	@test ft_ns1 == [3]
	# The union of the columns selected by a selector and its negation should be the set of all columns
	append!(ft_ns1, select_features(e1, s1))
	@test Set(ft_ns1) == Set(1:ncol(e1.data))

    # Test negation of name selector
	s2 = NameSelector(x -> occursin("Ft1", String(x)), "Keep Ft1")
	ns2 = negation(s2)

	@test ns2.description == "Do not "*s2.description

	ft_ns2 = select_features(e1, ns2)
	@test ft_ns2 == [2,3]
	# The union of the columns selected by a selector and its negation should be the set of all columns
	append!(ft_ns2, select_features(e1, s2))
	@test Set(ft_ns2) == Set(1:ncol(e1.data))
end

@testset "getdata" begin
    # Here we test whether we can access the data filtered
    # and use it for helper functions

    # Generate test data
	d = DataFrame(rand(12,2))
	rename!(d, [:Ft1, :Intensity_MedianIntensity_NeurDensity])
	d.Experiment = repeat(["Exp1", "Exp2"], 6)

	e1 = Experiment(d)

	# Filter data
	d1 = "Select only a single experiment"
	f1 = Filter("Exp1", :Experiment, description = d1)
	s1 = Selector(x -> eltype(x) <: Number,
	              description = "Keep numerical features")

	filter!(e1,[f1,s1])

	# Test `getdata`
	@test length(e1.selected_features) == 2
	@test length(e1.selected_entries) == 6
	@test getdata(e1) == e1.data[e1.selected_entries, e1.selected_features]

	# Test `logtransform`

	# Keep track of original values
	d_copy = copy(getdata(e1))
	e2 = deepcopy(e1)

	logtransform!(e1)

	# Check the transformation of one column
	@test getdata(e1).Ft1 == d_copy.Ft1 |> x -> log.(1 .+ x .- minimum(x))
	# Check that the original data is modified but not the copy
	@test getdata(e1).Ft1 == d.Ft1[(1:6)*2 .- 1]
	@test getdata(e2) == d_copy

	# Test `normtransform`

	# Keep track of original values
	e3 = deepcopy(e2)

	thr = sort(getdata(e2).Intensity_MedianIntensity_NeurDensity)[end-2]
	f2 = Filter(thr, :Intensity_MedianIntensity_NeurDensity,
	       compare = >=, description = "Top 3 values")

	normtransform!(e2,f2)

	# Must include at least one 0 per column (centered on median of a subset of entries)
	for c in eachcol(getdata(e2))
		@test 0 in c
	end
	# At most 1 value > 0 (after substracting median of top 3 values)
	@test sum(getdata(e2).Intensity_MedianIntensity_NeurDensity .> 0) <= 1
	# Check that the original data is modified but not the copy
	@test getdata(e3) != getdata(e2)
	@test getdata(e3) == d_copy

	# Test `decorrelate`

	# Add correlated columns
	d.Ft3 = -d.Ft1
	d.Ft4 = copy(d.Intensity_MedianIntensity_NeurDensity)
	d.Ft4[1] = 1
	d.Intensity_MedianIntensity_NeurDensity[1] = 0

	e4 = Experiment(d)
	filter!(e4,[f1,s1])
	decorrelate!(e4)
	# First column always kept, inverse column removed
	@test 1 in e4.selected_features
	@test !(4 in e4.selected_features)

	e4 = Experiment(d)
	filter!(e4,[f1,s1])
	decorrelate!(e4, threshold = 1)
	# Not perfectly correlated by design
	@test 5 in e4.selected_features

	e4 = Experiment(d)
	filter!(e4,[f1,s1])
	decorrelate!(e4, ordercol = [4,3,2,1])
	# Reversed order: first column always kept, inverse column removed
	@test 5 in e4.selected_features
	if cor(e4.data[e4.selected_entries, :Ft4],
		   e4.data[e4.selected_entries, :Ft3]) < 0.8
		@test 4 in e4.selected_features
	else
		@test !(4 in e4.selected_features)
	end

	@test !(1 in e4.selected_features)

	# Test sorting by mad
	e5 = Experiment(DataFrame([[0,2,4,6,8],[0,1,2,3,4],[0,3,6,9,12]]))
	decorrelate_by_mad!(e5)
	@test e5.selected_features == [3]

	# Test type issues
	d.Ft4 = [collect(1:9)..., missing, Inf, NaN]
	e5 = Experiment(d)

	# Includes String column
	@test_throws AssertionError decorrelate!(e5)

	filter!(e5,[f1,s1])
	decorrelate!(e5) # This should now work
	append!(e5.selected_features, 5)

	# Includes Inf values
	@test_throws AssertionError decorrelate!(e5)

	e5.data[11,5] = 3
	e6 = deepcopy(e5)
	decorrelate!(e6) # This should now work
	append!(e5.selected_entries, 10)

	# Includes Missing values
	@test_throws AssertionError decorrelate!(e5)

	e5.data[10,5] = 3
	e6 = deepcopy(e5)
	decorrelate!(e6) # This should now work
	append!(e5.selected_entries, 12)

	# Includes NaN values
	@test_throws AssertionError decorrelate!(e5)

	e5.data[12,5] = 3
	decorrelate!(e5) # This should now work

	# Test handling of non-float variables
	d.Ft2 = rand(1:4, size(d,1))
	e5 = Experiment(d)

	thr = sort(getdata(e5).Intensity_MedianIntensity_NeurDensity)[end-2]
	f3 = Filter(thr, :Intensity_MedianIntensity_NeurDensity,
	       compare = >=, description = "Top 3 values")

	@test_throws AssertionError logtransform!(e5)
	@test_throws AssertionError normtransform!(e5,f3)

	filter!(e5,[s1])

	e6 = deepcopy(e5)
	logtransform!(e5)
	normtransform!(e6,f3)

	@test eltype(getdata(e5).Ft2) == Float64
	@test eltype(getdata(e6).Ft2) == Float64
end

@testset "umap" begin
    # Create dataset with 4 random and 1 index column
    d = DataFrame(rand(40,4))
	rename!(d, Symbol.("Ft" .* string.(1:4)))
	d.Ft5 = 1:40;
	e = Experiment(d)

	# Set filters
	f = Filter(10, :Ft5, compare = >)
	s = negation(NameSelector(x -> occursin("Ft3", string(x))))
	filter!(e,[f,s])

	# Test dimensions and arguments
	@test size(umap(e)) == (2, 30)
	@test size(umap(e, 3, n_neighbors = 10, min_dist = 1, n_epochs = 10)) == (3, 30)
	@test_throws ArgumentError size(umap(e, 4))
end

@testset "biaseddistances" begin
	Random.seed!(777)

	# We want significantly more points than dimensions or the covariance
	# matrix can be singular and the results would not make any sense
	d = rand(100,5)
	@test distance_mahalanobis_center(d, 1:50, 1:50) + 1 ≈ 1

	d = hcat(d, (1:100)./50)
	# If some points are share with reference, the distance should be
	# smaller than to a set of completely different points
	dist1 = distance_mahalanobis_center(d, 51:100, 1:50)
	dist2 = distance_mahalanobis_center(d, 25:74, 1:50)
	@test dist1 > 0
	@test dist2 > 0
	@test dist1 > dist2

	dist0 = distance_mahalanobis_median(d, 1:50, 1:50)
	dist1 = distance_mahalanobis_median(d, 25:74, 1:50)
	dist2 = distance_mahalanobis_median(d, 51:100, 1:50)
	@test 0 < dist0 < dist1 < dist2

	# Testing permutation tests:
	# If both ref and pert are sampled from the same distribution,
	# p-value must be high.
	# If ref and pert are sampled from non-overlapping distributions,
	# p-value must be 0.
	d = rand(100,5)
	params = (d, 51:100, 1:50)
	@test 0 < mean(shuffled_distance_mahalanobis_center(params...) .<
	               distance_mahalanobis_center(params...)) < 1
	@test 0 < mean(shuffled_distance_mahalanobis_median(params...) .<
	               distance_mahalanobis_median(params...)) < 1

	d[51:100, :] .+= 2
	@test mean(shuffled_distance_mahalanobis_center(params...) .<
	           distance_mahalanobis_center(params...)) == 1
	@test mean(shuffled_distance_mahalanobis_median(params...) .<
	           distance_mahalanobis_median(params...)) == 1
end

@testset "robustdistances" begin
	Random.seed!(777)

	# First, RCall must be running correctly
	@test_throws RCall.REvalError R"""
	library(NotALibrary)
	"""

	# Robustbase must be installed
	R"""
	if (!require("robustbase")) install.packages("robustbase",
												  repos = "https://cloud.r-project.org")
	library(robustbase)
	"""
	d = rand(100,5)
	@test distance_robust_hellinger(d, 1:50, 1:50) + 1 ≈ 1

	d = hcat(d, (1:100)./50)
	# If some points are share with reference, the distance should be
	# smaller than to a set of completely different points
	dist0 = distance_robust_mahalanobis_median(d, 1:50, 1:50)
	dist1 = distance_robust_mahalanobis_median(d, 25:74, 1:50)
	dist2 = distance_robust_mahalanobis_median(d, 51:100, 1:50)
	@test 0 < dist0 < dist1 < dist2

	dist0 = distance_robust_hellinger(d, 1:50, 1:50)
	dist1 = distance_robust_hellinger(d, 25:74, 1:50)
	dist2 = distance_robust_hellinger(d, 51:100, 1:50)
	@test dist0 < dist1 < dist2

	# Testing permutation tests:
	# If both ref and pert are sampled from the same distribution,
	# p-value must be high.
	# If ref and pert are sampled from non-overlapping distributions,
	# p-value must be 0.
	d = rand(100,5)
	params = (d, 51:100, 1:50)
	@test 0 < mean(shuffled_distance_robust_mahalanobis_median(params...) .<
	               distance_robust_mahalanobis_median(params...)) < 1
	@test 0 < mean(shuffled_distance_robust_hellinger(params...) .<
	               distance_robust_hellinger(params...)) < 1

	d[51:100, :] .+= 1
	@test mean(shuffled_distance_robust_mahalanobis_median(params...) .<=
	           distance_robust_mahalanobis_median(params...)) > 0.9
	@test mean(shuffled_distance_robust_hellinger(params...) .<=
	           distance_robust_hellinger(params...)) > 0.9
end

@testset "robust_morphological_perturbation_value" begin
	Random.seed!(777)

	d = DataFrame(rand(250,5))
	d.Condition = sample('A':'D', 250);

	# Make one condition stand out
	d[d.Condition .== 'D',1:5] .+= 1;

	e = Experiment(d)

	# Select numerical columns
	slt = NameSelector(x -> x != "Condition")
	select_features!(e, slt)
	# Filter out some rows
	# Useful to check that the RMPV computation correctly works on subset of the data
	row_filter = Filter(0.9, :x1, compare = (x,y) -> round(x, digits = 1) != y)
	filter_entries!(e, row_filter)

	# Define reference condition
	f = Filter('C', :Condition)

	@test_throws DomainError robust_morphological_perturbation_value(e,
																     :Condition,
																     f,
																     dist = :IncorrectValue)

	rmpv = robust_morphological_perturbation_value(e, :Condition, f)

	# 4 conditions, 3 columns
	@test size(rmpv) == (4,3)

	# Distance of control to itself == 0
	@test rmpv[rmpv.Condition .== 'C', :Distance][1] + 1 ≈ 1

	# Shifted distribution should be the most different
	@test maximum(rmpv.Distance) == rmpv.Distance[rmpv.Condition .== 'D'][1]

	# A, B and C have the same distribution but D has a different one
	@test rmpv.RMPV[rmpv.Condition .== 'C'][1] > 0.1
	@test rmpv.RMPV[rmpv.Condition .== 'D'][1] < 0.1

	rmpv2 = robust_morphological_perturbation_value(e,
													:Condition,
													'C',
													nb_rep = 50,
													dist = :RobMedMahalanobis)

	# A, B and C have the same distribution but D has a different one
	@test rmpv2.RMPV[rmpv2.Condition .== 'C'][1] > 0.1
	@test rmpv2.RMPV[rmpv2.Condition .== 'D'][1] < 0.1

	rmpv3 = robust_morphological_perturbation_value(e,
													:Condition,
													'C',
													nb_rep = 50,
													dist = :MedMahalanobis)

	# A, B and C have the same distribution but D has a different one
	@test rmpv3.RMPV[rmpv3.Condition .== 'C'][1] > 0.1
	@test rmpv3.RMPV[rmpv3.Condition .== 'D'][1] < 0.1

	rmpv4 = robust_morphological_perturbation_value(e,
													:Condition,
													'C',
													nb_rep = 50,
													dist = :CenterMahalanobis)

	# A, B and C have the same distribution but D has a different one
	@test rmpv4.RMPV[rmpv4.Condition .== 'C'][1] > 0.1
	@test rmpv4.RMPV[rmpv4.Condition .== 'D'][1] < 0.1

	# Test reproducibility
    Random.seed!(777)
    rmpv_run1 = robust_morphological_perturbation_value(e,
                                                        :Condition,
                                                        'C',
                                                        nb_rep = 4,
                                                        dist = :RobMedMahalanobis)
    Random.seed!(777)
    rmpv_run2 = robust_morphological_perturbation_value(e,
                                                        :Condition,
                                                        'C',
                                                        nb_rep = 4,
                                                        dist = :RobMedMahalanobis)
    Random.seed!(777)
    rmpv_run3 = robust_morphological_perturbation_value(e,
                                                    :Condition,
                                                    'C',
                                                    nb_rep = 4,
                                                    dist = :RobMedMahalanobis,
                                                    r_seed = false)
    @test rmpv_run1 == rmpv_run2
    # NB: without changing the seeds, it is possible that
    # the MCD converges so the following is a bad test:
    # @test rmpv_run1 != rmpv_run3

    Random.seed!(777)
    rmpv_run1 = robust_morphological_perturbation_value(e,
                                                        :Condition,
                                                        'C',
                                                        nb_rep = 4,
                                                        dist = :RobHellinger)
    Random.seed!(777)
    rmpv_run2 = robust_morphological_perturbation_value(e,
                                                        :Condition,
                                                        'C',
                                                        nb_rep = 4,
                                                        dist = :RobHellinger)
    Random.seed!(777)
    rmpv_run3 = robust_morphological_perturbation_value(e,
                                                    :Condition,
                                                    'C',
                                                    nb_rep = 4,
                                                    dist = :RobHellinger,
                                                    r_seed = false)
    @test rmpv_run1 == rmpv_run2
    # @test rmpv_run1 != rmpv_run3
end

@testset "parallel_rmpv" begin
	d = DataFrame(rand(100,5))
	d.Condition = sample('A':'D', 100);

	e = Experiment(d)

    addprocs(4)
    pool = CachingPool(workers())
    @everywhere using BioProfiling

	slt = NameSelector(x -> x != "Condition")
	select_features!(e, slt)
	f = Filter('C', :Condition)

	rmpv = robust_morphological_perturbation_value(e,
												   :Condition,
												   f,
												   process_pool = pool)

	# 4 conditions, 3 columns
	@test size(rmpv) == (4,3)
end


text:
using InformationDistances
using CodecXz: XzCompressor
using TranscodingStreams: Noop

using Test

struct MockCompressor <: InformationDistances.AbstractCompressor end

function InformationDistances.compressed_length(::MockCompressor, data::InformationDistances.ByteData)

    return length(data)
end

s0 = ""
s1 = "x"
s2 = "xy"

@testset "InformationDistances.jl" begin

    @testset "AbstractCompressor interface" begin
        @test compressed_lengths(MockCompressor(), [s0, s1, s2]) ==
            map(s -> compressed_length(MockCompressor(), s), [s0, s1, s2])

        @test compressed_length(MockCompressor(), s2) ==
            compressed_length(MockCompressor(), codeunits(s2))
    end

    @testset "NormalizedCompressionDistance" begin

        @test NormalizedCompressionDistance() isa NormalizedCompressionDistance{CodecCompressor{XzCompressor}}

        @test NormalizedCompressionDistance(MockCompressor()) isa
                NormalizedCompressionDistance{MockCompressor}

        len_s1 = compressed_length(MockCompressor(), s1)
        len_s2 = compressed_length(MockCompressor(), s2)
        len_s1s2 = compressed_length(MockCompressor(), s1 * s2)

        @test NormalizedCompressionDistance(MockCompressor())(s1, s2) ==
            (len_s1s2 - min(len_s1, len_s2)) / max(len_s1, len_s2)

        @test NormalizedCompressionDistance()("", "") == 0.0
        @test NormalizedCompressionDistance()(s2, s2) == 0.0
        # NCD is not symmetric in general, but it should be for the empty string
        @test NormalizedCompressionDistance()("", s2) ==
            NormalizedCompressionDistance()(s2, "")
        @test 0.0 <= NormalizedCompressionDistance()(s1, s2) <= 1.0
    end

    @testset "LibDeflateCompressor" begin

        @test LibDeflateCompressor().compresslevel == 12
        @test LibDeflateCompressor(;compresslevel=5).compresslevel == 5

        compressor = LibDeflateCompressor()

        data_CodeUnits = codeunits.([s0, s1, s2])
        data_Vector = Vector(data_CodeUnits)

        @test compressed_lengths(compressor, data_CodeUnits) ==
              compressed_lengths(compressor, data_Vector) ==
              map(s -> compressed_length(compressor, s), data_CodeUnits) ==
              map(s -> compressed_length(compressor, s), data_Vector)

    end

    @testset "CodecCompressor" begin

        @test CodecCompressor{Noop}() isa CodecCompressor{Noop}
        @test collect(CodecCompressor{Noop}(arg1="1", arg2=2).kwargs) ==
            [:arg1 => "1", :arg2 => 2]

        compressor = CodecCompressor{Noop}()

        data_CodeUnits = codeunits.([s0, s1, s2])
        data_Vector = Vector(data_CodeUnits)

        @test compressed_lengths(compressor, data_CodeUnits) ==
              compressed_lengths(compressor, data_Vector) ==
              map(s -> compressed_length(compressor, s), data_CodeUnits) ==
              map(s -> compressed_length(compressor, s), data_Vector)

    end

end

ext: jl


text:
# split data into train and test
"""
    partitionTrainTest(data, at)
Split a dataset into train and test datasets (array).
"""
function partitionTrainTest(data, at)
    # if size(data) |> length != 1
    #     error("Invalid data type")
    # end
    n = length(data)
    idx = Random.shuffle(1:n)
    train_idx = view(idx, 1:floor(Int, at*n))
    test_idx = view(idx, (floor(Int, at*n)+1):n)
    data[train_idx], data[test_idx]
end


text:
"""
    histrecon_u(dims, C2, S2, philen)

Reconstruct a binary image on the basis of input C2 and S2. The size of the
reconstructed image is determined by ```dims```, and the number of pixels set to
one by ```philen```

This returns a full suite of statistics about the reconstructed object to check
stability.

# Example
    dims, C2, S2, philen = get_C2_S2(fname)
    guess, S2n, C2n, S2_BN1, C2_BN1, SN1 = histrecon((200, 200)), C2, S2, 12000)
"""
function histrecon(dims, C2, S2, philen)
    guess, pix = make_rand_im(philen, dims)
    #display(heatmap(guess, color=:grays, aspect_ratio=1))
    maxrng = Int64(round(sqrt(dims[1]^2 + dims[2]^2)+1))
    wpick = rand(1:length(pix[1]))
    bpick = rand(1:length(pix[2]))
    kern = CartesianIndices((-1:1, -1:1))

    original = label_components(guess, trues(3,3))
    SN = SN_comp(dims, maxrng)
    typeof(SN)
    S2_BN = S2_initialise(pix, maxrng)
    C2_BN = C2_initialise(guess, maxrng)
    S2n = naninf(S2_finalise(SN, S2_BN, guess))
    C2n = naninf(S2_finalise(SN, C2_BN, guess))
    #print(SN, "\n\n", S2_BN, "\n\n", C2_BN)
    Err1 = 2.0
    count = 1
    cont_break=0;
    T0 = 1.0
    cont_surf_opt=0
    pflag = 0
    true_count = 0
    while Err1 > 0.0001
        if (cont_surf_opt !== 1)
            wpick, bpick = pick_rand_coord(pix)
            #println("stoachastic mode")
        else
            surf_res = surf_opt(guess)
            shortlist = surf_rand(surf_res, pix)
            wpick, bpick = pick_rand_coord(shortlist)
            wpick, bpick = find_equivalent(shortlist, pix, (wpick, bpick))
            #println("surface optimisation mode")
            #println("recon shortlist: \twlength = ", length(shortlist[1]), " blength = ", length(shortlist[2]))
            #println("wpick = ", wpick, "\tbpick = ", bpick)

        end

        pix  = swap_pix(guess, pix, wpick, bpick)
        modified = label_components(guess, trues(3,3))
        S2_BN_d = update_S2_BN(S2_BN, pix, wpick, bpick, philen, maxrng)
        C2_BN_d = update_C2_BN(C2_BN, pix, original, modified, wpick, bpick, maxrng)
        #println(length(S2_BN), "\n", length(S2_BN_d))
        S2_BN .+= S2_BN_d
        C2_BN .+= C2_BN_d
        #println("\n", S2_BN_d, "\n")
        #println("C2_BN_d = ", C2_BN_d)
        S2n = naninf(S2_finalise(SN, S2_BN, guess))
        C2n = naninf(S2_finalise(SN, C2_BN, guess))
        #println("C2n = ", C2n)

        true_count +=1

        Err2 =  ediff_C2(C2n, C2) + (ediff_S2(S2n, S2))# +
        #SN_comp(size(guess), length(guess), maxrng)
        delErr = Err2 - Err1
        if count == 1
            T0=-abs(delErr)/log(0.5)
            println("T0 = ", T0)
        end
        T = T0*0.99^count
        if (rand(0:1) < min(exp(-delErr/T), 1))
            count += 1
            cont_break=0
            #print("change commited at count ", count, "\t")
            original = modified
            Err1 = Err2
        else
            swap_pix(guess, pix, wpick, bpick)
            S2_BN .-= S2_BN_d
            C2_BN .-= C2_BN_d
            cont_break += 1

        end
         # we bake in the change here - this must be undone
        if (mod(true_count, 1000) == 0) & (count !=pflag)
            println("\ncount is ", count, "\t  True count is ", true_count, "\t prob of rejection =", min(exp(delErr/T), 1), "\t  S2 Err = ", ediff_S2_fft(S2n,  S2)*1250/length(guess), "\t C2 Err = ",  ediff_C2(C2n, C2))
            println("S2_delerr = ", ediff_S2(S2n, S2), "\tC2_delerr = ", ediff_C2(C2n, C2))
            println("C2[1] = ", C2n[1], "\t\tS2[1] = ", S2n[1])
            if cont_surf_opt==1
                println("Surface optimisation mode")
            else
                println("Stoachastic mode")
            end
            #println(Err1)
            pflag = count
        end
        if count==4000
            cont_surf_opt=1
            #println("Changed to surface optimisation mode")
        end
        if true_count >= 10

            break
        end
    end
    return guess, S2n, C2n, S2_BN, C2_BN, SN
end

"""
    get_C2_S2(fname)

Compute C2 and S2 for an existing binary image ```fname```
"""
function get_C2_S2(fname)
    test = loadim(fname)
    dims = size(test)
    maxrng = Int64(round(sqrt(dims[1]^2 + dims[2]^2)+1))
    philen = sum(test)
    SN = SN_comp(dims, maxrng)

    idx = CartesianIndices(dims)
    wpix = findall(x->x==true, test)
    bpix = setdiff(idx, wpix)
    pix = (wpix, bpix)
    S2N = S2_initialise(pix, maxrng)
    S2 = naninf(S2_finalise(SN, S2N, test))
    #S2[length(S2)÷2:end] .= 0
    C2N = C2_initialise(test, maxrng)
    C2 = naninf(S2_finalise(SN, C2N, test))
    return (dims, C2, S2, philen)
end


text:
using Documenter, FiniteDifferences

makedocs(
    modules=[FiniteDifferences],
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages=["Home" => "index.md", "API" => "pages/api.md"],
    authors="Invenia Labs",
    checkdocs=:exports,
    linkcheck=true,
    repo="https://github.com/JuliaDiff/FiniteDifferences.jl/blob/{commit}{path}#L{line}",
    sitename="FiniteDifferences.jl",
    strict=true,
)

deploydocs(repo="github.com/JuliaDiff/FiniteDifferences.jl.git")


text:
using QuasiArrays, IntervalSets, Test
import QuasiArrays: ApplyQuasiArray

@testset "Continuous" begin
    @testset "Inclusion" begin
        # note behaviour changes with ContinuumArrays
        x = Inclusion(0..1)
        @test eltype(x) == Int
        @test axes(x) ≡ (x,)
        @test_throws InexactError x[0.1]
        @test x[1] ≡ 1
        @test_throws BoundsError x[2]
        @test Base.unsafe_getindex(x,2) ≡ 2

        for x in (Inclusion{Float64}(0..1),Inclusion(0.0..1))
            @test Inclusion(x) ≡ Inclusion{Float64}(x) ≡ convert(Inclusion,x) ≡
                    convert(Inclusion{Float64},x) ≡ convert(AbstractQuasiArray,x) ≡
                    convert(AbstractQuasiArray{Float64},x) ≡ convert(AbstractQuasiVector{Float64},x) ≡
                    convert(AbstractQuasiVector,x) ≡ x
            @test axes(x) ≡ (x,)
            @test x[0.1] ≡ 0.1
            @test x[1] ≡ 1.0
            @test_throws BoundsError x[2]
            @test Base.unsafe_getindex(x,2) ≡ 2.0
            @test x[BigFloat(π)/4] == π/4
        end

        @test Inclusion{Float64}(0..1) == Inclusion(0.0..1)
        @test Inclusion{Float64}(0..1) !== Inclusion(0.0..1)

        @test findfirst(isequal(0.1), Inclusion(0.0..1)) ≡ findlast(isequal(0.1), Inclusion(0..1)) ≡ 0.1
        @test findfirst(isequal(2.3), Inclusion(0.0..1)) ≡ findlast(isequal(2.3), Inclusion(0..1)) ≡ nothing
        @test findall(isequal(0.1), Inclusion(0.0..1)) == [0.1]
        @test findall(isequal(2.3), Inclusion(0.0..1)) == Float64[]
    end

    @testset "QuasiDiagonal" begin
        x = Inclusion(0.0..1)
        D = QuasiDiagonal(x)
        @test D[0.1,0.2] ≡ 0.0
        @test D[0.1,0.1] ≡ 0.1
        @test D[1,1] ≡ 1.0
        @test_throws BoundsError D[1.1,1.2]
    end
end


text:
"""
```
    Observations
```
    A one-dimensional array of 2D points wrapped in a tuple called `data`.

    The reason I wrap the points in a tuple is that I intend to refactor `Observations` into a separate package
    in the future, and different elements of the tuple may represent corresponding points between images, for example.
"""
struct Observations{N,T} <: AbstractObservations
    data::NTuple{N,T}
end

# Convert an N x 2 matrix to a one-dimensional array of 2D points (i.e. each point is a row of the N x 2 matrix) and wrap in a tuple.
Observations(data::AbstractMatrix; colwise = true) = Observations(tuple(svectors(transpose(data), Val{2}())))

"""
```
    UncertainObservations
```
    A one-dimensional array of 2D points wrapped in a tuple and a one-dimensional array of corresponding 2D covariance matrices (also wrapped in a tuple).

    The reason I wrap the points in a tuple is that I intend to refactor `UncertainObservations` into a separate package
    in the future, and different elements of the tuple may represent corresponding points between images, for example.
"""
struct UncertainObservations{N,T₁, T₂} <: AbstractUncertainObservations
    data::NTuple{N,T₁}
    covariance_matrices::NTuple{N,T₂}
end


text:
mutable struct GridworldPlanner <: Planner
    gridworld::Gridworld
    num_expansions::Int64
    actions::Vector{GridworldAction}
    residuals::Matrix{Float64}
end

function GridworldPlanner(gridworld::Gridworld, num_expansions::Int64 = 10)
    residuals = fill(0.0, (gridworld.size, gridworld.size))
    GridworldPlanner(gridworld, num_expansions, getActions(gridworld), residuals)
end

function getHeuristic(planner::GridworldPlanner, state::GridworldState)
    manhattan_dist = manhattanHeuristic(planner, state)
    return manhattan_dist + planner.residuals[state.x, state.y]
end

function manhattanHeuristic(planner::GridworldPlanner, state::GridworldState)
    abs(planner.gridworld.goal_state.x - state.x) +
    abs(planner.gridworld.goal_state.y - state.y)
end

function getSuccessors(
    planner::GridworldPlanner,
    state::GridworldState,
    action::GridworldAction,
)
    next_state, cost = step(planner.gridworld, state, action)
    return next_state, cost
end

function checkGoal(planner::GridworldPlanner, state::GridworldState)
    checkGoal(planner.gridworld, state)
end

function getActions(planner::GridworldPlanner, state::GridworldState)
    planner.actions
end

function act(planner::GridworldPlanner, state::GridworldState)
    best_action, info = astar(planner, planner.num_expansions, state)
    return best_action, info
end

function updateResiduals!(planner::GridworldPlanner, info::Dict)
    best_node_f = info["best_node_f"]
    for node in info["closed"]
        state = node.state
        g = node.g
        planner.residuals[state.x, state.y] =
            best_node_f - g - manhattanHeuristic(planner, state)
    end
end


text:
@testset "with_container()" begin
    configs = [
        DockerConfig(; base_image = "julia:latest", verbose = true),
        DockerConfig(; base_image = "julia:latest", verbose = false),
        DockerConfig(;
            base_image = "julia:latest",
            verbose = true,
            docker_build_stdout = Base.devnull,
            docker_build_stderr = Base.devnull,
        ),
    ]
    for config in configs
        with_container() do container
            code = """
            println("This was a success.")
            """
            @test success(container, config, `julia -e $(code)`)
        end

        with_container() do container
            code = """
            throw(ErrorException("This was a failure."))
            """
            @test !success(container, config, `julia -e $(code)`)
        end

        with_container() do container
            code = """
            println("This was a success.")
            """
            p = run(container, config, `julia -e $(code)`; wait = false)
            wait(p)
            @test success(p)
        end

        with_container() do container
            code = """
            throw(ErrorException("This was a failure."))
            """
            p = run(container, config, `julia -e $(code)`; wait = false)
            wait(p)
            @test !success(p)
        end
    end
end


text:
"""
    save(f, tree)

Wrapper function for `FileTrees.save` to configure whether to use parallel
saving using `FileTrees` or to just use a basic serial implementation.
Typically the simpler serial code will be faster unless the project is very
large.
"""
function save(f, tree)
    if DAGGER[]
        FileTrees.save(f, tree)
    else
        for file in FileTrees.files(tree)
            dir = dirname(file)
            isdir(dir) || mkpath(dir)
            f(file)
        end
    end
end

# # HTML

"""
Convert the given `src` project to a collection of HTML files.
"""
function html(src, dst=nothing; keywords...)
    p = Project(src; keywords...)
    p === nothing && return src
    h = p.env["publish"]["html"]
    h["template"]["string"] = String(exec(p.tree[h["template"]["file"]][]))
    sandbox(dst) do
        default_html_pages(p) # TODO: handle as part of template?
        tree = rename(p.tree, pwd())
        mapping = page_neighbours(p.pages)
        searchmd = p"search.md"
        tree = touch(tree, searchmd; value=load_markdown("# Search\n

"))
        mapping[searchmd] = (prev=searchmd, next=searchmd)
        save(f -> _html(p, tree, f, mapping), tree)
    end
    return src
end

function init(p::Project, ::typeof(html); port=nothing, dir=nothing, kws...)
    html(p, dir; kws...)
    LiveServer.serve(; port=port, dir=dir)
end

_html(p::Project, t::FileTree, f::File, m::Dict) = _html(p, exec(f[]), relative(path(f), basename(t)), m)

function _html(p::Project, node::CommonMark.Node, path::AbstractPath, mapping::Dict)
    if haskey(mapping, path)
        dst = with_extension(path, "html")
        ## Setup.
        let pub = p.env["publish"]
            pub["mapping"] = mapping
            pub["smartlink-engine"] = (_,_,n,_)->toc_link(n, p, p.env["publish"], path)
            ast = exec(p.tree[pub["toc"]][])
            pub["html"]["toc"] = CommonMark.html(ast, p.env["publish"])
            pub["html"]["prev"], pub["html"]["next"] = mapping[path]
        end
        ## Writing.
        relative_paths(p, path) do pub
            pub["html"]["prev"] = with_extension(pub["html"]["prev"], "html")
            pub["html"]["next"] = with_extension(pub["html"]["next"], "html")
            dir, name = splitdir(dst)
            cd(isempty(dir) ? "." : dir) do
                open(name, "w") do io
                    pub["template-engine"] = Mustache.render
                    pub["smartlink-engine"] = (_,_,n,_)->html_link(n, p, pub, path)
                    CommonMark.html(io, node, pub)
                end
            end
        end
        write("search.json", JSON.json(json_search_data(p)))
        ## Cleanup.
        let pub = p.env["publish"]
            delete!(pub, "mapping")
            delete!(pub, "template-engine")
            delete!(pub, "smartlink-engine")
            delete!(pub["html"], "toc")
            delete!(pub["html"], "prev")
            delete!(pub["html"], "next")
        end
    end
    return nothing
end
_html(::Project, data::Vector{UInt8}, path::AbstractPath, ::Dict) = write(path, data)
_html(p::Project, t::FileTrees.Thunk, path::AbstractPath, env::Dict) = _html(p, exec(t), path, env)
_html(::Project, ::Any, ::AbstractPath, ::Dict) = nothing

function toc_link(node, project, pub, path)
    obj = deepcopy(node.t)
    ## Change to toc location based on the current path.
    reltoc = relpath(pub["toc"], dirname(joinpath(".", string(path))))
    ## Adjust the toc link's path based on the new toc root.
    obj.destination = joinpath(dirname(reltoc), obj.destination)
    return html_link(obj, node, project, pub, path)
end

html_link(node, project, pub, path) = html_link(deepcopy(node.t), node, project, pub, path)

function html_link(obj, node, project, pub, path)
    if obj.destination == "#"
        function docs_func!(literal::AbstractString)
            dict = frontmatter(exec(project.tree[path][]))
            module_binding = binding(get(dict, "module", findmodule(project.env)))
            if Docs.defined(module_binding)
                target_binding = binding(Docs.resolve(module_binding), literal)
                if Docs.defined(target_binding)
                    rel = relpath("docstrings/$target_binding.html", string(dirname(path)))
                    obj.destination = rel
                    @goto END
                end
            end
            @warn "cross-reference link '$literal' on page '$path' cound not be found."
            @label END
            return nothing
        end
        function header_func!(literal::AbstractString)
            slug = CommonMark.slugify(literal)
            for each in project.pages, (node, enter) in exec(project.tree[each][])
                if enter && get(node.meta, "id", nothing) == slug
                    name = with_extension(relpath(each, dirname(path)), "html")
                    obj.destination = "$name#$slug"
                    obj.title = ""
                    @goto END
                end
            end
            @warn "cross-reference link '$literal' on page '$path' could not be found."
            @label END
            return nothing
        end
        ## `#` is used for cross-references. The link is determined by either
        ## the provided `.title` field of the link, or the contents of the link.
        if isempty(obj.title)
            ## No title provided so we use the contents of the link.
            (!CommonMark.isnull(node.first_child) && node.first_child.t isa CommonMark.Code) ?
                docs_func!(node.first_child.literal) : header_func!(node.first_child.literal)
        else
            ## The `.title` is available, so use that to determine the link.
            m = match(r"^`(.+)`$", obj.title)
            m === nothing ? header_func!(obj.title) : docs_func!(m[1])
        end
    elseif startswith(obj.destination, "#")
        ## Skip these kind of links, they're just page-local.
    else
        dst = Path(normpath(dirname(string(path)), obj.destination))
        if haskey(pub["mapping"], dst)
            ## If it's in the project's page mapping then we change the extension.
            obj.destination = with_extension(obj.destination, "html")
        end
    end
    return obj
end

function default_html_pages(p::Project)
    if !isempty(p.pages)
        content =
            """
            <!DOCTYPE html>

            <meta http-equiv = "refresh" content = "0; url = $(with_extension(first(p.pages), "html"))" />

            """
        write("index.html", content)
    end
    return nothing
end

"""
Extract JSON search data from a project for use in lunr.js.
"""
function json_search_data(project::Project)
    dict = Dict{String,String}()
    root = dirname(joinpath(".", project.env["publish"]["toc"]))
    for page in project.pages
        path = relpath(string(page), root)
        path = with_extension(path, "html")
        id = path
        if hasfile(project.tree, page)
            for (node, enter) in exec(project.tree[page][])
                if enter
                    if haskey(node.meta, "id")
                        id = "$path#$(node.meta["id"])"
                    end
                    if (node.t isa CommonMark.Text || node.t isa CommonMark.Code)
                        if haskey(dict, id)
                            dict[id] = "$(dict[id]) $(node.literal)"
                        else
                            dict[id] = node.literal
                        end
                    end
                end
            end
        end
    end
    json = Dict{String,String}[]
    for (id, body) in dict
        push!(json, Dict("id" => id, "body" => body))
    end
    return json
end

# # PDF

"""
Convert the given `src` project to a PDF file.
"""
function pdf(src, dst=nothing; keywords...)
    p = Project(src; keywords...)
    p === nothing && return nothing
    sandbox(dst) do
        tree = rename(p.tree, pwd())
        save(f -> _pdf(p, tree, f), tree)
        tocroot = joinpath(".", dirname(p.env["publish"]["toc"]))
        io = IOBuffer()
        println(io, "```{=latex}")
        for page in p.pages
            rpath = relpath(string(page), tocroot)
            name, _ = splitext(rpath)
            name = unix_style_path(name) # Path adjustments for Windows.
            println(io, "\\include{$name}")
        end
        println(io, "```")
        ast = load_markdown(io)
        project_file = p.env["name"] * ".tex"
        t = p.env["publish"]["latex"]
        t["template"]["string"] = String(exec(p.tree[t["template"]["file"]][]))
        p.env["publish"]["template-engine"] = Mustache.render
        open(project_file, "w") do handle
            CommonMark.latex(handle, ast, p.env["publish"])
        end
        ## Build the final PDF document using tectonic.
        Tectonic.tectonic() do path
            run(`$path $project_file`)
        end
    end
    return src
end

function init(p::Project, ::typeof(pdf); dir=nothing, kws...)
    pdf(p, dir; kws...)
    pdf_viewer = Sys.iswindows() ? "start" : Sys.isapple() ? "open" : "xdg-open"
    run(`$pdf_viewer $(joinpath(dir, p.env["name"] * ".pdf"))`)
end

_pdf(p::Project, t::FileTree, f::File) = _pdf(p, exec(f[]), relative(path(f), basename(t)))

function _pdf(p::Project, node::CommonMark.Node, path::AbstractPath)
    pub = p.env["publish"]
    pub["smartlink-engine"] = (_, _, n, _) -> tex_link(n)
    dst = with_extension(path, "tex")
    dir, name = splitdir(dst)
    cd(isempty(dir) ? "." : dir) do
        open(name, "w") do io
            CommonMark.latex(io, node, pub)
        end
    end
end
_pdf(::Project, data::Vector{UInt8}, path::AbstractPath) = write(path, data)
_pdf(::Project, ::Any, ::AbstractPath) = nothing

function tex_link(n::CommonMark.Node)
    ## TODO: make links work.
    obj = deepcopy(n.t)
    obj.destination = ""
    return obj
end


text:
# EUtils
# ======
#
# APIs for E-Utilities.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Entrez Programming Utilities (or E-Utilities) module.

The `EUtils` module provides a programming interface to the Entrez databases at
NCBI. Nine functions are exported from this module:
1. `einfo`: database statistics
2. `esearch`: text search
3. `epost`: UID upload
4. `esummary`: document summary download
5. `efetch`: data record download
6. `elink`: Entrez link
7. `egquery`: global query
8. `espell`: spelling suggestion
9. `ecitmatch`: batch citation searching in PubMed

See "Entrez Programming Utilities Help"
(https://www.ncbi.nlm.nih.gov/books/NBK25501/) for more details. Especially,
["E-utilities Quick Start"](https://www.ncbi.nlm.nih.gov/books/NBK25500/) is a
good starting point and ["A General Introduction to the
E-utilities"](https://www.ncbi.nlm.nih.gov/books/NBK25497/) will serve useful
information about its concepts and functions. The implemented APIs are based on
the manual on January 23, 2015.
"""
module EUtils

export
    einfo,
    esearch,
    epost,
    esummary,
    efetch,
    elink,
    egquery,
    espell,
    ecitmatch,
    set_context!

import XMLDict
import JSON
import HTTP

const baseURL = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/"

# APIs of E-utilities
# -------------------

"""
    einfo(ctx=Dict(); params...)

Retrieve a list of databases or statistics for a database.

Parameters: db, version, retmode.
"""
function einfo(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    return HTTP.request("GET", string(baseURL, "einfo.fcgi"), query=params)
end

"""
    esearch(ctx=Dict(); params...)

Retrieve a list of UIDs matching a text query.

Parameters: db, term, usehistory, WebEnv, query_key, retstart, retmax, rettype,
retmode, sort, field, datetype, reldate, mindate, maxdate.
"""
function esearch(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    res = HTTP.request("GET", string(baseURL, "esearch.fcgi"), query=params, retry_non_idempotent=true)
    if get(params, :usehistory, "") == "y"
        set_context!(ctx, res)
    end
    return res
end

"""
    epost(ctx=Dict(); params...)

Upload or append a list of UIDs to the Entrez History server.

Parameters: db, id, WebEnv.
"""
function epost(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    body = HTTP.escapeuri(params)
    # added user-agent header as workaround for EOFError - HTTP.jl issue #342
    headers = Dict("user-agent"=>"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_14_1) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/70.0.3538.102 Safari/537.36")
    res = HTTP.request("POST", string(baseURL, "epost.fcgi"), headers, body=body)
    set_context!(ctx, res)
    return res
end

"""
    esummary(ctx=Dict(); params...)

Retrieve document summaries for a list of UIDs.

Parameters: db, id, query_key, WebEnv, retstart, retmax, retmode, version.
"""
function esummary(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    body = HTTP.escapeuri(params)
    return HTTP.request("POST", string(baseURL, "esummary.fcgi"), body=body)
end

"""
    efetch(ctx=Dict(); params...)

Retrieve formatted data records for a list of UIDs.

Parameters: db, id, query_key, WebEnv, retmode, rettype, retstart, retmax,
strand, seq_start, seq_stop, complexity.
"""
function efetch(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    body = HTTP.escapeuri(params)
    return HTTP.request("POST", string(baseURL, "efetch.fcgi"), body=body, retry_non_idempotent=true)
end

"""
    elink(ctx=Dict(); params...)

Retrieve UIDs linked to an input set of UIDs.

Parameters: db, dbfrom, cmd, id, query_key, WebEnv, linkname, term, holding,
datetype, reldate, mindate, maxdate.
"""
function elink(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    body = HTTP.escapeuri(params)
    return HTTP.request("POST", string(baseURL, "elink.fcgi"), body=body)
end

"""
    egquery(ctx=Dict(); params...)

Retrieve the number of available records in all databases by a text query.

Parameters: term.
"""
function egquery(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    return HTTP.request("GET", string(baseURL, "egquery.fcgi"), query=params)
end

"""
    espell(ctx=Dict(); params...)

Retrieve spelling suggestions.

Parameters: db, term.
"""
function espell(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    return HTTP.request("GET", string(baseURL, "espell.fcgi"), query=params)
end

"""
    ecitmatch(ctx=Dict(); params...)

Retrieve PubMed IDs that correspond to a set of input citation strings.

Parameters: db, rettype, bdata.
"""
function ecitmatch(ctx::AbstractDict=empty_context(); params...)
    params = process_parameters(params, ctx)
    return HTTP.request("GET", string(baseURL, "ecitmatch.cgi"), query=params)
end

# Create an empty context.
function empty_context()
    return Dict{Symbol,Any}()
end

# Set :WebEnv and :query_key values to the context `ctx` from `res`.
function set_context!(ctx, res)
    if res.status != 200
        return ctx
    end

    # extract WebEnv and query_key from the response
    contenttype = Dict(res.headers)["Content-Type"]
    body = deepcopy(res.body)
    data = String(res.body)

    if startswith(contenttype, "text/xml")
        doc = XMLDict.parse_xml(data)
        ctx[:WebEnv] = doc["WebEnv"]
        ctx[:query_key] = doc["QueryKey"]
    elseif startswith(contenttype, "application/json")
        dict = JSON.parse(data)
        ctx[:WebEnv] = dict["esearchresult"]["webenv"]
        ctx[:query_key] = dict["esearchresult"]["querykey"]
    end
    res.body = body

    return ctx
end

# Process query parameters.
function process_parameters(params, ctx)
    # merge context `ctx` into `params`
    params = merge(ctx, Dict(params))

    # flatten a set of IDs into a comma-separated string
    if haskey(params, :id)
        ids = params[:id]
        if isa(ids, AbstractString)
            ids = [ids]
        end
        params[:id] = join([string(id) for id in ids], ',')
    end

    # normalize the usehistory parameter
    if haskey(params, :usehistory) && isa(params[:usehistory], Bool)
        if params[:usehistory]
            params[:usehistory] = "y"
        else
            delete!(params, :usehistory)
        end
    end

    # stringify all values
    for (key, val) in params
        params[key] = string(val)
    end

    return params
end

end


text:
# dense
function reach_blocks_wrapping_effect!(
        ϕ::AbstractMatrix{NUM},
        Xhat0::Vector{<:LazySet{NUM}},
        U::Union{ConstantInput, Nothing},
        overapproximate::Function,
        overapproximate_inputs::Function,
        n::Int,
        N::Int,
        output_function::Nothing, # ignored
        blocks::AbstractVector{Int},
        partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
        δ::NUM,
        termination::Function,
        res::Vector{<:AbstractReachSet}
       )::Tuple{Int, Bool} where {NUM}
    X_store = CartesianProductArray(Xhat0)
    t0 = zero(δ)
    t1 = δ
    store!(res, 1, X_store, t0, t1, NUM)
    terminate, skip = termination(1, X_store, t0)
    if terminate
        return 1, skip
    end

    b = length(partition)
    Xhatk = Vector{LazySet{NUM}}(undef, b)

    if U != nothing
        Whatk = Vector{LazySet{NUM}}(undef, b)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
    end

    arr_length = (U == nothing) ? length(partition) : length(partition) + 1
    arr = Vector{LazySet{NUM}}(undef, arr_length)
    k = 2
    p = Progress(N, 1, "Computing successors ")
    @inbounds while true
        ProgressMeter.update!(p, k)
        for i in 1:b
            bi = partition[blocks[i]]
            for (j, bj) in enumerate(partition)
                arr[j] = ϕ[bi, bj] * array(set(res[k-1]))[j]
            end
            if U != nothing
                arr[arr_length] = Whatk[i]
            end
            Xhatk[i] = overapproximate(blocks[i], MinkowskiSumArray(arr))
        end
        X_store = CartesianProductArray(copy(Xhatk))
        t0 = t1
        t1 += δ
        store!(res, k, X_store, t0, t1, NUM)

        terminate, skip = termination(k, X_store, t0)
        if terminate
            break
        end

        if U != nothing
            for i in 1:b
                bi = partition[blocks[i]]
                Whatk[i] =
                    overapproximate_inputs(k, blocks[i], row(ϕ, bi) * inputs)
            end
        end

        k += 1
    end

    return k, skip
end


text:
# Plotting the data structure for the state of an SRA analysis
#   See also: sra_model_object.jl for the saved states
#   See also: sra_model_plotting.jl for plotting methods
# include("../utils/main_algorithm_utils.jl")
include("../paper/scripts/paper_settings.jl");
using BSON: @save

#####
##### Model fitting
#####
"""
Step 1a: Fit a naive sindy model to the data
    Note: only called first time around

fit_first_model(m::sra_stateful_object, initial_noise)
Input:
    initial_noise - the initial guess for the noise std

Sensitive to parameters:
noise_factor
initial_subsampling
use_turing

See methods:
print_current_equations
print_true_equations

Expectation: this model will not be very good!
"""
function fit_first_model(m::sra_stateful_object, initial_noise)
    # TODO: move noise level to parameter object
    p = m.parameters
    # Change noise estimate based on previous guess
    ensemble_p = p.sindyc_ensemble_parameters
    use_control = is_using_control(m)
    if use_control
        ensemble_p[:selection_dist] = Normal(0.0, initial_noise)
    end

    # Note: both of the branches call sindyc_ensemble
    if p.initial_subsampling
        println("Initial subsampling; control=$use_control")
        (out) = calc_best_random_subsample(p.model_template,
                                        m.dat,
                                        m.numerical_grad,
                                        num_pts=1000,
                                        num_subsamples=10,
                                        val_list=p.sindy_terms_list;
                                        sindyc_ensemble_params=ensemble_p,
                                        use_control=use_control,
                                        max_error=initial_noise)
        sindy_model = out[:best_model]
        all_criteria = out[:all_errL2]
        best_criterion = minimum(all_criteria)
    else
        (sindy_model,best_criterion,all_criteria,all_models) =
            sindyc_ensemble(model_template,
                    m.dat,
                    m.numerical_grad,
                    p.sindy_library,
                    p.sindy_terms_list;
                    var_names=p.variable_names,
                    ensemble_p...)
    end
    m.sindy_model = sindy_model
    m.ensemble_sindy_criteria = all_criteria
    m.best_sindy_criteria = best_criterion
    m.i = 1 # Resets counter
end

"""
Step 1b: Fit a model based on subsampled data
    Note: Mostly same as above, for non-initial iterations

Sensitive to parameters:
TODO: noise_factor

See methods:
print_current_equations
print_true_equations
plot_subsampled_derivatives
plot_subsampled_simulation
"""
function fit_model(m::sra_stateful_object)
    p = m.parameters
    # New: change noise estimate based on previous guess
    ensemble_p = p.sindyc_ensemble_parameters
    # ensemble_p[:selection_dist] = Normal(0.0, m.noise_guess)
    # Same
    if is_using_control(m)
            f = sindyc_ensemble
        else
            f = sindy_ensemble
        end
    (sindy_model,best_criterion,all_criteria,all_models) =
         f(m.parameters.model_template,
         m.dat[:, m.subsample_ind],
         m.numerical_grad[:, m.subsample_ind],
         m.parameters.sindy_terms_list; ensemble_p...)
    # (sindy_model,best_criterion,all_criteria,all_models) =
    #     sindyc_ensemble(
    #             m.dat[:, m.subsample_ind],
    #             m.numerical_grad[:, m.subsample_ind],
    #             p.sindy_library,
    #             p.sindy_terms_list;
    #             var_names=p.variable_names,
    #             ensemble_p...)
    # Only save if this is an improvement!
    if best_criterion < m.best_sindy_criteria
        if !m.is_saved
            save_model_variables(m)
        end
        println("Iteration $(m.i+1) improved the fit (err=$best_criterion)")
        # These overwrite the previous step, which has been saved
        m.sindy_model = sindy_model
        m.sindy_dat = nothing # Reset
        m.ensemble_sindy_criteria = all_criteria
        m.best_sindy_criteria = best_criterion
        m.is_saved = false
        m.i += 1
        is_improved = true
    else
        println("Iteration $(m.i+1) did not improve the fit")
        # println("(old=$(m.best_sindy_criteria);new=$best_criterion)")
        println("Found equation:")
        print_equations(sindy_model, digits=4)
        is_improved = false
    end
    return all_models, is_improved
end

#####
##### Residual calculation
#####
"""
Step 2: Subsample based on previous model residuals

 calculate_subsampled_ind(m::sra_stateful_object,
                        not_enough_pts_mode="error")

Input:
    not_enough_pts_mode
        - if "error", then raises an error if not enough points
        - if "warning", then raises a warning and uses all of the points

Sensitive to parameters:
noise_factor
use_turing

See methods:
plot_subsampled_points
"""
function calculate_subsampled_ind(m::sra_stateful_object,
                                not_enough_pts_mode="warning")
    p = m.parameters
    if !m.parameters.use_turing
        sindy_grad = m.sindy_model(m.dat, 0)
        residual = m.numerical_grad .- sindy_grad
        # TODO: don't hardcode the index here
        noise_guess = p.noise_factor*abs(median(residual[1,:])) # Really should use Turing here
    else
        error("Turing not implemented in stateful version")
    end
    # TODO: don't hardcode length (and other parameters)
    # TODO: Should this really be a loop?
    accepted_ind = []
    tmp_noise_factor = 1.0
    i = 1
    while length(accepted_ind) < p.num_pts
        accepted_ind = subsample_using_residual(residual,
                    tmp_noise_factor*noise_guess, min_length=4);
        i += 1;
        i > 20 && break
        tmp_noise_factor *= 3.0;
    end

    if length(accepted_ind) < p.num_pts+p.start_ind-1
        if not_enough_pts_mode == "error"
            error("DataError: Not enough accepted points; increase noise_factor or decrease num_pts")
        else
            @warn("DataError: Not enough accepted points with noise level $(tmp_noise_factor*noise_guess)")
        end
        m.subsample_ind = accepted_ind[p.start_ind:end]
    else
        m.subsample_ind = accepted_ind[p.start_ind:p.num_pts+p.start_ind-1]
    end
    m.noise_guess = noise_guess
end

#####
##### Saving
#####
"""
Save current step for continuation
"""
function save_model_variables(m::sra_stateful_object)
    pushfirst!(m.all_sindy_models, m.sindy_model)
    pushfirst!(m.all_ensemble_sindy_criteria, m.ensemble_sindy_criteria)
    pushfirst!(m.all_best_sindy_criteria, m.best_sindy_criteria)
    pushfirst!(m.all_sindy_dat, m.sindy_dat)
    pushfirst!(m.all_residuals, m.residual)
    pushfirst!(m.all_subsample_ind, m.subsample_ind)
    pushfirst!(m.all_noise_guesses, m.noise_guess)

    m.is_saved = true;
end

#####
##### Evaluation
#####

"""
Calculate the L2 error between the found and true models

function calc_coefficient_error(m::sra_stateful_object,
                                t::sra_truth_object)

return err
"""
function calc_coefficient_error(m::sra_stateful_object,
                                t::sra_truth_object)
    err = calc_coefficient_error(m.sindy_model, t.core_dyn_true)
    return err
end

#####
##### Other steps
#####
"""
Helper function: Integrate the current model

Note: not necessary for simulation, just for plotting

See methods:
plot_sindy_model
"""
function simulate_model(m::sra_stateful_object)
    # p = m.parameters
    # Callback aborts if it blows up
    # condition(u,t,integrator) = any(abs.(u).>1e4)
    # cb = DiscreteCallback(condition, terminate!)
    # prob = ODEProblem(m.sindy_model, m.u0, m.tspan, [0], callback=cb)
    # m.sindy_dat = Array(solve(prob, Tsit5(), saveat=m.ts));
    m.sindy_dat = simulate_model(m.sindy_model, m.u0)
end

"""
Calculates the residual of the data and the current model
"""
function calc_residual(m::sra_stateful_object)
    sindy_grad = m.sindy_model(m.dat, 0)
    residual = m.numerical_grad .- sindy_grad
    return residual
end

"""
Saves the data in a standardized format in four different files

Note: requires the paper_settings.jl file to specify the DAT_FOLDERNAME variable

Usage:
save_for_plotting(m::sra_stateful_object,
                t::sra_truth_object,
                this_dat_name="test")
"""
function save_for_plotting(m::sra_stateful_object,
                            t::sra_truth_object,
                            this_dat_name="test")
    this_dat_name = DAT_FOLDERNAME*this_dat_name;

    dat_raw = t.dat_raw
    fname = this_dat_name*"uncontrolled.bson"
    @save fname dat_raw

    dat, U_true = m.dat, t.U_true
    fname = this_dat_name*"ode_vars.bson";
    @save fname dat U_true

    noise_guess, accepted_ind = m.noise_guess, m.subsample_ind
    residual = calc_residual(m)
    fname = this_dat_name*"naive_vars_bayes.bson";
    @save fname noise_guess residual accepted_ind #ctr_guess

    dat_ctr = simulate_model(m)
    ctr_guess2 = process_residual(residual, noise_guess) # TODO: should be Bayesian
    fname = this_dat_name*"ctr_vars_sindy.bson";
    @save fname dat_ctr ctr_guess2 #sindy_grad_ctr
end

# Note: just saving matrices in a standard way
# function save_for_plotting(noise_vals,
#                 vec_naive, std_naive, vec_err, std_err,
#                 vec_deriv, std_deriv, vec_naive_deriv, std_naive_deriv,
#                 this_dat_name="test")
#     this_dat_name = DAT_FOLDERNAME*this_dat_name;
#
#     fname = this_dat_name*"coefficients.bson"
#     @save fname noise_vals vec_naive std_naive vec_err std_err
#
#     fname = this_dat_name*"derivatives.bson";
#     @save fname vec_deriv std_deriv vec_naive_deriv std_naive_deriv
# end

function save_for_plotting(noise_vals, control_signal_vals, coef_norm,
                all_naive_err, all_err,
                all_err_deriv_subsample, all_naive_err_deriv,
                this_dat_name="test")
    this_dat_name = DAT_FOLDERNAME*this_dat_name;

    fname = this_dat_name*"metadata.bson"
    @save fname noise_vals control_signal_vals coef_norm

    fname = this_dat_name*"coefficients.bson"
    @save fname all_naive_err all_err

    fname = this_dat_name*"derivatives.bson";
    @save fname all_err_deriv_subsample all_naive_err_deriv
end

###
### Utilities
###

is_using_control(m::sra_stateful_object) =
    isa(m.parameters.model_template, sindycModel)

##### Export
export fit_first_model, fit_model, calculate_subsampled_ind,
    save_model_variables, simulate_model, save_for_plotting,
    calc_residual, calc_coefficient_error, is_using_control


text:
# ---
# title: 39. Combination Sum
# id: problem39
# author: zhwang
# date: 2022-02-17
# difficulty: Medium
# categories: Array, Backtracking
# link: <https://leetcode.com/problems/combination-sum/description/>
# hidden: true
# ---
#
# Given an array of **distinct** integers `candidates` and a target integer
# `target`, return _a list of all **unique combinations** of _`candidates`
# _where the chosen numbers sum to_`target` _._ You may return the combinations
# in **any order**.
#
# The **same** number may be chosen from `candidates` an **unlimited number of
# times**. Two combinations are unique if the frequency of at least one of the
# chosen numbers is different.
#
# It is **guaranteed** that the number of unique combinations that sum up to
# `target` is less than `150` combinations for the given input.
#
#
#
# **Example 1:**
#
#
#
#     Input: candidates = [2,3,6,7], target = 7
#     Output: [[2,2,3],[7]]
#     Explanation:
#     2 and 3 are candidates, and 2 + 2 + 3 = 7. Note that 2 can be used multiple times.
#     7 is a candidate, and 7 = 7.
#     These are the only two combinations.
#
#
# **Example 2:**
#
#
#
#     Input: candidates = [2,3,5], target = 8
#     Output: [[2,2,2,2],[2,3,3],[3,5]]
#
#
# **Example 3:**
#
#
#
#     Input: candidates = [2], target = 1
#     Output: []
#
#
# **Example 4:**
#
#
#
#     Input: candidates = [1], target = 1
#     Output: [[1]]
#
#
# **Example 5:**
#
#
#
#     Input: candidates = [1], target = 2
#     Output: [[1,1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= candidates.length <= 30`
#   * `1 <= candidates[i] <= 200`
#   * All elements of `candidates` are **distinct**.
#   * `1 <= target <= 500`
#
#
## @lc code=start
using LeetCode

function combination_sum_39(
    nums::Union{Vector{Int},SubArray}, target::Int
)::Vector{Vector{Int}}
    target == 0 && return [Int[]]
    isempty(nums) && return Vector{Int}[]
    res, num = Vector{Int}[], nums[end]
    for i in 0:(target ÷ num)
        for vec in combination_sum_39(@view(nums[1:(end - 1)]), target - i * num)
            append!(res, [append!(vec, fill(num, i))])
        end
    end
    return res
end
## @lc code=end


text:
package java.lang;

public class NumberFormatException extends IllegalArgumentException {
    public NumberFormatException() { }
    public NumberFormatException(String s) { super(s); }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using mbRK
using mbPred
using mbUtil

# procedure usually underestimates length of period
"""
    findCycle(H, t0, y0, transientIterations, transientStepSize,
	    steadyStateIterations, steadyStateStepSize)
Returns all points from the `transientIterations`-times numerical integration of `H` with
fixed step-size `transientStepSize`, and all points from the `steadyStateIterations`-times
integration with fixed step-size `steadyStateStepSize` (assuming it hit the steady state
part), and the periodicity of the found cycle.
"""
function findCycle(
	H::Function, t₀::Float64, y₀::Vector,
	TransientIterations::Integer, TransientStepSize::Float64,
	SteadyStateIterations::Integer, SteadyStateStepSize::Float64
	)

	local dataTransient, dataSteadyState
	dataTransient, dataSteadyState = Array{Float64}(0,length(y₀)), Array{Float64}(0,length(y₀))
	function callbackTransient(t,y,ft) dataTransient = [dataTransient; y'] end
	function callbackSteadyState(t,y,ft) dataSteadyState = [dataSteadyState; y'] end

	# transient part
	if TransientIterations <= 0
		dataTransient = y₀'
	else
		# rk1(H, t₀, y₀, TransientStepSize/10, predCount(TransientIterations); callback=callbackTransient)
		rk4(H, t₀, y₀, TransientStepSize, predCount(TransientIterations); callback=callbackTransient)
	end

	# (hopefully) steady state part
	local t₁, y₁
	t₁, y₁ = TransientIterations*TransientStepSize, squeeze(dataTransient[end,:]', (2))
	# rk1(H, t₁, y₁, SteadyStateStepSize/10, predCount(SteadyStateIterations); callback=callbackSteadyState)
	rk4(H, t₁, y₁, SteadyStateStepSize, predCount(SteadyStateIterations); callback=callbackSteadyState)

	# autocorrelate trajectory, standardize, take minimum
	# 	(and-connect them, only interested in periodicity of ALL components)
	local ac, peaks, P, cyc

	ac = autocorrelation(dataSteadyState)
	ac = mapslices(minimum, standardize(ac), [2])

	# find peaks in AC, take second highest (not interested in 0-shift)
	peaks = findPeaks(ac)[2:end-1]
	P = peaks[findmax(ac[peaks])[2]]

	return dataTransient, dataSteadyState, P
end

"""
    prepareCycle(data, h, P[, fac])
cut single cycle of length `P*fac` from data, resample, shift s.t. `X(0)≈0`, Fourier transform.
"""
function prepareCycle(data::Array{Float64}, h::Float64, P::Integer; fac::Integer=1)
	local ω, cyc, y

	# cut out singe period
	ω = 2pi / (P*fac * h)
	cyc = data[end-P*fac+2:end,:]

	# resample
	# y = mapslices(V->mbUtil.interpolate(V, precision, linspace(1.0, size(cyc,1), samples)), cyc, [1])

	# local I = linspace(1.0, size(cyc,1), Samples)
	# y = mapslices(V->(V[floor(Integer,I)] + V[ceil(Integer,I)])/2, cyc, [1])

	# rotate s.t. x(0)≈0
	# this (↓) is ultimately important! optimization is super sensitive to it...
	y = circshift(cyc,[-findmin(abs(cyc[:,1]))[2]+1])

	return y, ω
end

"""
    findCyclePoincare(F, y[, plane, clusterRating, nIntersections,
	    maxCycles, sampleSize, transientIterations, transientStepSize,
		steadyStateStepSize])
extracts a single cycle of the steady state of ode `F` using poincare cuts through
the `plane`.
"""
function findCyclePoincare(
		F :: Function,
		y₀ :: Vector;
		plane :: Function = first,	# x ∈ plane ⇔ plane(x) = 0
		clusterRating :: Function = x -> var(x) / length(x),	# use penalty for more points, i.e., normalize
		nIntersections :: Int = 120,
		maxCycles :: Int = 30,
		sampleSize :: Int = 512,
		transientIterations :: Int = 2000,
		transientStepSize :: Float64 = 0.1,
		steadyStateStepSize :: Float64 = 0.1)

	dim = length(y₀)
	flagIntersected = false
	y⁻ = y₀
	y₁ = ode(F, y₀, transientStepSize, function (y, i, t, V)
		plane(y⁻) ≤ 0 ≤ plane(y) && (flagIntersected = true)
		y⁻ = y
		return i < transientIterations
	end)

	flagIntersected || error("No intersection in transient phase.")

	# generate intersections with plane
	h = steadyStateStepSize
	y⁻ = y₁
	intersections = zeros(dim + 1, 0)
	ode(F, y₁, h, (y, i, t, V) -> begin
		if (sign(plane(y)) > 0 &&  sign(plane(y⁻)) <= 0)
			# find h′ ∈ [0,h] s.t. plane(y⁻ + h′⋅V(h′,y⁻)) = 0
			step(h′) = y⁻ + h′*V(h′,y⁻)
			h′ = mbUtil.bisection(h′ -> plane(step(h′)), 0, h, ϵ=1e-4)
			y′ = step(h′)
			intersections = [intersections [y′; t+h′]]
		end
		y⁻ = y
		i > 100*transientStepSize/steadyStateStepSize*transientIterations && error("No intersection in steady state phase.")
		nIntersections > size(intersections, 2)
	end)

	#scatter(intersections[2,:], intersections[3,:])

	# find number of cycles s.t. rating is minimum
	n = size(intersections, 2)
	ratings = map(1:maxCycles) do i
		mean( [ sum(var(intersections[1:end-1, j:i:n], [2])) for j in 1:i ] )
	end

	mn₁,m₁ = minimum(ratings), mean(ratings)
	period = findfirst(1:maxCycles) do i
		X = ratings[i:i:end]
		# mn₂,m₂ = minimum(ratings), mean(ratings)
		std(X)/mean(X) < .5 && all(X) do x
			abs(x-mn₁)/m₁ < .5
		end
	end

	if period == 0
		period = 1
		warn("A valid period could not be found.")
	end

	# cut out steady state trajectory
	y₂ = intersections[1:dim,1]
	T = intersections[dim+1,1+period] - intersections[dim+1,1]
	dataSteadyState = zeros(0,dim)
	ode(F, y₂, steadyStateStepSize, (y, i, t, V) -> begin
		dataSteadyState = [dataSteadyState; y']
		return t < T
	end)

	dataSteadyState, period, 2π/T
end

# externalize or merge with rk4
function ode(f, y, h, pred)
	RK  = [1/2 0 0 0; 0 1/2 0 0; 0 0 1 0; 1/6 2/6 2/6 1/6]; # classic
	s   = size(RK, 1);
	dim = length(y);
	t   = 0

	function V(h, y)
		k = zeros(dim, s);

		k[:,1] = f(t, y);
		for j=1:s-1
			k[:,j+1] = f( t + h * sum(RK[j,:]), y + h * k * RK[j,:]' );
		end

		k * RK[s,:]';
	end

	i = 0
	while pred(y, i, t, V) !== false
		y = y + h * V(h, y);
		t += h
		i += 1
	end

	y
end

# y = rand(3)
# T = zeros(0,3)
# y = ode((t,v)->f(t,[v;350]), y, .01, (y, i, t, V) -> begin
# 	global T = [T; y']
# 	return i < 1000
# end); println(y)
#
# f = figure()
# gca(projection="3d")
# plot(T[:,1], T[:,2], T[:,3])
# f[:show]()


text:
module DerivativeTest

import Calculus

using Test
using Random
using ForwardDiff
using DiffTests

include(joinpath(dirname(@__FILE__), "utils.jl"))

Random.seed!(1)

########################
# test vs. Calculus.jl #
########################

const x = 1

for f in DiffTests.NUMBER_TO_NUMBER_FUNCS
    println("  ...testing $f")
    v = f(x)
    d = ForwardDiff.derivative(f, x)
    @test isapprox(d, Calculus.derivative(f, x), atol=FINITEDIFF_ERROR)

    out = DiffResults.DiffResult(zero(v), zero(v))
    out = ForwardDiff.derivative!(out, f, x)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.derivative(out), d)
end

for f in DiffTests.NUMBER_TO_ARRAY_FUNCS
    println("  ...testing $f")
    v = f(x)
    d = ForwardDiff.derivative(f, x)

    @test !(eltype(d) <: ForwardDiff.Dual)
    @test isapprox(d, Calculus.derivative(f, x), atol=FINITEDIFF_ERROR)

    out = similar(v)
    out = ForwardDiff.derivative!(out, f, x)
    @test isapprox(out, d)

    out = DiffResults.DiffResult(similar(v), similar(d))
    out = ForwardDiff.derivative!(out, f, x)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.derivative(out), d)
end

for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS
    println("  ...testing $f!")
    m, n = 3, 2
    y = fill(0.0, m, n)
    f = x -> (tmp = similar(y, promote_type(eltype(y), typeof(x)), m, n); f!(tmp, x); tmp)
    v = f(x)
    cfg = ForwardDiff.DerivativeConfig(f!, y, x)
    d = ForwardDiff.derivative(f, x)

    fill!(y, 0.0)
    @test isapprox(ForwardDiff.derivative(f!, y, x), d)
    @test isapprox(v, y)

    fill!(y, 0.0)
    @test isapprox(ForwardDiff.derivative(f!, y, x, cfg), d)
    @test isapprox(v, y)

    out = similar(v)
    fill!(y, 0.0)
    ForwardDiff.derivative!(out, f!, y, x)
    @test isapprox(out, d)
    @test isapprox(v, y)

    out = similar(v)
    fill!(y, 0.0)
    ForwardDiff.derivative!(out, f!, y, x, cfg)
    @test isapprox(out, d)
    @test isapprox(v, y)

    out = DiffResults.DiffResult(similar(v), similar(d))
    out = ForwardDiff.derivative!(out, f!, y, x)
    @test isapprox(v, y)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.derivative(out), d)

    out = DiffResults.DiffResult(similar(v), similar(d))
    out = ForwardDiff.derivative!(out, f!, y, x, cfg)
    @test isapprox(v, y)
    @test isapprox(DiffResults.value(out), v)
    @test isapprox(DiffResults.derivative(out), d)
end

@testset "exponential function at base zero" begin
    @test (x -> ForwardDiff.derivative(y -> x^y, -0.5))(0.0) === -Inf
    @test (x -> ForwardDiff.derivative(y -> x^y,  0.0))(0.0) === -Inf
    @test (x -> ForwardDiff.derivative(y -> x^y,  0.5))(0.0) === 0.0
    @test (x -> ForwardDiff.derivative(y -> x^y,  1.5))(0.0) === 0.0
end

@testset "dimension error for derivative" begin
    @test_throws DimensionMismatch ForwardDiff.derivative(sum, fill(2pi, 3))
end

end # module


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: codestar_notifications

using AWS.Compat
using AWS.UUIDs
"""
    CreateNotificationRule()

Creates a notification rule for a resource. The rule specifies the events you want notifications about and the targets (such as SNS topics) where you want to receive them.

# Required Parameters
- `DetailType`: The level of detail to include in the notifications for this resource. BASIC will include only the contents of the event as it would appear in AWS CloudWatch. FULL will include any supplemental information provided by AWS CodeStar Notifications and/or the service for the resource for which the notification is created.
- `EventTypeIds`: A list of event types associated with this notification rule. For a list of allowed events, see EventTypeSummary.
- `Name`: The name for the notification rule. Notifictaion rule names must be unique in your AWS account.
- `Resource`: The Amazon Resource Name (ARN) of the resource to associate with the notification rule. Supported resources include pipelines in AWS CodePipeline, repositories in AWS CodeCommit, and build projects in AWS CodeBuild.
- `Targets`: A list of Amazon Resource Names (ARNs) of SNS topics to associate with the notification rule.

# Optional Parameters
- `ClientRequestToken`: A unique, client-generated idempotency token that, when provided in a request, ensures the request cannot be repeated with a changed parameter. If a request with the same parameters is received and a token is included, the request returns information about the initial request that used that token.  The AWS SDKs prepopulate client request tokens. If you are using an AWS SDK, an idempotency token is created for you.
- `Status`: The status of the notification rule. The default value is ENABLED. If the status is set to DISABLED, notifications aren't sent for the notification rule.
- `Tags`: A list of tags to apply to this notification rule. Key names cannot start with \"aws\".
"""

create_notification_rule(DetailType, EventTypeIds, Name, Resource, Targets; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/createNotificationRule", Dict{String, Any}("DetailType"=>DetailType, "EventTypeIds"=>EventTypeIds, "Name"=>Name, "Resource"=>Resource, "Targets"=>Targets, "ClientRequestToken"=>string(uuid4())); aws_config=aws_config)
create_notification_rule(DetailType, EventTypeIds, Name, Resource, Targets, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/createNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DetailType"=>DetailType, "EventTypeIds"=>EventTypeIds, "Name"=>Name, "Resource"=>Resource, "Targets"=>Targets, "ClientRequestToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    DeleteNotificationRule()

Deletes a notification rule for a resource.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule you want to delete.

"""

delete_notification_rule(Arn; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
delete_notification_rule(Arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), args)); aws_config=aws_config)

"""
    DeleteTarget()

Deletes a specified target for notifications.

# Required Parameters
- `TargetAddress`: The Amazon Resource Name (ARN) of the SNS topic to delete.

# Optional Parameters
- `ForceUnsubscribeAll`: A Boolean value that can be used to delete all associations with this SNS topic. The default value is FALSE. If set to TRUE, all associations between that target and every notification rule in your AWS account are deleted.
"""

delete_target(TargetAddress; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteTarget", Dict{String, Any}("TargetAddress"=>TargetAddress); aws_config=aws_config)
delete_target(TargetAddress, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/deleteTarget", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("TargetAddress"=>TargetAddress), args)); aws_config=aws_config)

"""
    DescribeNotificationRule()

Returns information about a specified notification rule.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule.

"""

describe_notification_rule(Arn; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/describeNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
describe_notification_rule(Arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/describeNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), args)); aws_config=aws_config)

"""
    ListEventTypes()

Returns information about the event types available for configuring notifications.

# Optional Parameters
- `Filters`: The filters to use to return information by service or resource type.
- `MaxResults`: A non-negative integer used to limit the number of returned results. The default number is 50. The maximum number of results that can be returned is 100.
- `NextToken`: An enumeration token that, when provided in a request, returns the next batch of the results.
"""

list_event_types(; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listEventTypes"; aws_config=aws_config)
list_event_types(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listEventTypes", args; aws_config=aws_config)

"""
    ListNotificationRules()

Returns a list of the notification rules for an AWS account.

# Optional Parameters
- `Filters`: The filters to use to return information by service or resource type. For valid values, see ListNotificationRulesFilter.  A filter with the same name can appear more than once when used with OR statements. Filters with different names should be applied with AND statements.
- `MaxResults`: A non-negative integer used to limit the number of returned results. The maximum number of results that can be returned is 100.
- `NextToken`: An enumeration token that, when provided in a request, returns the next batch of the results.
"""

list_notification_rules(; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listNotificationRules"; aws_config=aws_config)
list_notification_rules(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listNotificationRules", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns a list of the tags associated with a notification rule.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) for the notification rule.

"""

list_tags_for_resource(Arn; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTagsForResource", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
list_tags_for_resource(Arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), args)); aws_config=aws_config)

"""
    ListTargets()

Returns a list of the notification rule targets for an AWS account.

# Optional Parameters
- `Filters`: The filters to use to return information by service or resource type. Valid filters include target type, target address, and target status.  A filter with the same name can appear more than once when used with OR statements. Filters with different names should be applied with AND statements.
- `MaxResults`: A non-negative integer used to limit the number of returned results. The maximum number of results that can be returned is 100.
- `NextToken`: An enumeration token that, when provided in a request, returns the next batch of the results.
"""

list_targets(; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTargets"; aws_config=aws_config)
list_targets(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/listTargets", args; aws_config=aws_config)

"""
    Subscribe()

Creates an association between a notification rule and an SNS topic so that the associated target can receive notifications when the events described in the rule are triggered.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule for which you want to create the association.
- `Target`:

# Optional Parameters
- `ClientRequestToken`: An enumeration token that, when provided in a request, returns the next batch of the results.
"""

subscribe(Arn, Target; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/subscribe", Dict{String, Any}("Arn"=>Arn, "Target"=>Target); aws_config=aws_config)
subscribe(Arn, Target, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/subscribe", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "Target"=>Target), args)); aws_config=aws_config)

"""
    TagResource()

Associates a set of provided tags with a notification rule.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule to tag.
- `Tags`: The list of tags to associate with the resource. Tag key names cannot start with \"aws\".

"""

tag_resource(Arn, Tags; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/tagResource", Dict{String, Any}("Arn"=>Arn, "Tags"=>Tags); aws_config=aws_config)
tag_resource(Arn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/tagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    Unsubscribe()

Removes an association between a notification rule and an Amazon SNS topic so that subscribers to that topic stop receiving notifications when the events described in the rule are triggered.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule.
- `TargetAddress`: The ARN of the SNS topic to unsubscribe from the notification rule.

"""

unsubscribe(Arn, TargetAddress; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/unsubscribe", Dict{String, Any}("Arn"=>Arn, "TargetAddress"=>TargetAddress); aws_config=aws_config)
unsubscribe(Arn, TargetAddress, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/unsubscribe", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "TargetAddress"=>TargetAddress), args)); aws_config=aws_config)

"""
    UntagResource()

Removes the association between one or more provided tags and a notification rule.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule from which to remove the tags.
- `TagKeys`: The key names of the tags to remove.

"""

untag_resource(Arn, TagKeys; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/untagResource", Dict{String, Any}("Arn"=>Arn, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(Arn, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/untagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn, "TagKeys"=>TagKeys), args)); aws_config=aws_config)

"""
    UpdateNotificationRule()

Updates a notification rule for a resource. You can change the events that trigger the notification rule, the status of the rule, and the targets that receive the notifications.  To add or remove tags for a notification rule, you must use TagResource and UntagResource.

# Required Parameters
- `Arn`: The Amazon Resource Name (ARN) of the notification rule.

# Optional Parameters
- `DetailType`: The level of detail to include in the notifications for this resource. BASIC will include only the contents of the event as it would appear in AWS CloudWatch. FULL will include any supplemental information provided by AWS CodeStar Notifications and/or the service for the resource for which the notification is created.
- `EventTypeIds`: A list of event types associated with this notification rule.
- `Name`: The name of the notification rule.
- `Status`: The status of the notification rule. Valid statuses include enabled (sending notifications) or disabled (not sending notifications).
- `Targets`: The address and type of the targets to receive notifications from this notification rule.
"""

update_notification_rule(Arn; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/updateNotificationRule", Dict{String, Any}("Arn"=>Arn); aws_config=aws_config)
update_notification_rule(Arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = codestar_notifications("POST", "/updateNotificationRule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Arn"=>Arn), args)); aws_config=aws_config)


text:
using Test
using Smoothers

@testset "stl" begin

    @test_throws AssertionError stl(rand(100),10; ns=5)

    x = stl(rand(100),10)
    @test x isa Matrix
    @test size(x) == (100,3)

end


text:
println("Advent of Code 2017, day 6")

f = open("input6")
input = read(f,String)
close(f)

banks = [parse(Int,n) for n in split(input)]
n = length(banks)

# Part 1

function next(j)
  if j<n return j+1 else return 1 end
end

function redistribute(bk)
  b = copy(bk)
  (m,j) = findmax(b)
  b[j] = 0
  while m > 0
    j = next(j)
    b[j] += 1
    m -= 1
  end
  return b
end

cycles = [ banks ]
b = redistribute(banks)
i = 1
while !(b in cycles)
  append!(cycles,[b])
  global b = redistribute(b)
  global i += 1
end

println("Part 1: ",i)

# Part 2

index = findfirst(x -> x==b,cycles)

println("Part 2: ",i-index+1)

max_stars_repo_name: venanzio/AoC

max_issues_repo_name: venanzio/AoC

max_forks_repo_name: venanzio/AoC


text:
# Autogenerated wrapper script for wannier90_jll for x86_64-linux-musl-libgfortran3
export libwannier, wannier90

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("wannier90")
JLLWrappers.@declare_library_product(libwannier, "libwannier.so")
JLLWrappers.@declare_executable_product(wannier90)
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libwannier,
        "lib/libwannier.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        wannier90,
        "bin/wannier90",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()