text:
@testset "Helpers" begin
    ch = DiscordChannel(; id="255", type=0, guild_id="1")
    r = Role(; id=0xff, name="foo")
    u = User(; id="255", username="foo")
    e1 = Emoji(; id="255", name="foo")
    e2 = Emoji(; id=nothing, name="foo")

    @testset "String representations" begin
        @test string(ch) == "<#255>"
        @test string(r) == "<@&255>"
        @test string(u) == "<@255>"
        @test string(e1) == "<:foo:255>"
        @test string(e2) == "foo"
        m = Member(u, "foo", [], now(), true, true)
        @test string(m) == "<@!255>"
        m = Member(u, nothing, [], now(), true, true)
        @test string(m) == string(u)
        m = Member(u, missing, [], now(), true, true)
        @test string(m) == string(u)
        # Make sure that string interpolation works normally (julia#21982).
        @test "$u" == string(u)
        @test "foo $u bar" == "foo " * string(u) * " bar"
    end

    @testset "split_message" begin
        # A vector is returned even if no splitting happens.
        chunks = split_message("foo")
        @test chunks == String["foo"]

        # Formatting is preserved between chunks.
        chunks = split_message(string(repeat('.', 1995), "**foobar**"))
        @test length(chunks) == 2
        @test chunks[2] == "**foobar**"
        chunks = split_message(string(repeat('.', 1995), "```julia\n1 + 1\n```"))
        @test chunks[2] == "```julia\n1 + 1\n```"

        # Chunks are separated by spaces.
        chunks = split_message(string(repeat('.', 1995), " foo bar baz"))
        @test length(chunks) == 2
        @test chunks[1] == string(repeat('.', 1995), " foo")
        @test chunks[2] == "bar baz"
    end

    @testset "plaintext" begin
        msg = Message(;
            id="1",
            channel_id="1",
            content="<@255> <@!255>",
            mentions=[JSON.lower(u)],
        )
        @test plaintext(msg) == "@foo @foo"
    end

    @testset "@fetch/@fetchval" begin
        @testset "wrapfn!" begin
            ex = :(foo(1, 2, 3))

            # Calls to the designated functions get wrapped in another function.
            @test wrapfn!(ex, (:foo,), :bar) == :($(esc(:bar))($(esc(:(foo(1, 2, 3))))))
            # But only if we specify that function to be wrapped.
            @test wrapfn!(ex, (:baz,), :bar) == ex

            ex = :(x = foo(1))
            expected = :($(esc(:x)) = $(esc(:baz))($(esc(:(foo(1))))))

            # The function operates recursively.
            @test wrapfn!(ex, (:foo,), :baz) == expected
        end

        @testset "validate_fetch" begin
            @test_throws ArgumentError validate_fetch(:foo, :(begin end))
            @test_throws ArgumentError validate_fetch(:create, :foo, :(begin end))
            @test_throws ArgumentError validate_fetch(:create, :(x = 1))
            @test_throws ArgumentError validate_fetch(:create, :retrieve)
            validate_fetch(:create, :(begin end))
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ cbcca572-7614-11eb-2fe5-699acc8e2746
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

	# using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	# using Random
	using Test
	using Printf
	using Plots
end

# ╔═╡ f1eefcb4-7614-11eb-0173-d5112c9036d7
include("./utils.jl")

# ╔═╡ 84a95dde-7614-11eb-20da-8191d2c5103e
md"""
## Predicting house prices using k-nearest neighbors regression

In this notebook, we will implement k-nearest neighbors regression. We will:
  - Find the k-nearest neighbors of a given query input
  - Predict the output for the query input using the k-nearest neighbors
  - Choose the best value of k using a validation set
"""

# ╔═╡ db51098e-7614-11eb-2dcd-57f4526747f6
md"""
### Load in house sales data (train, valid and test sets)
"""

# ╔═╡ f20a9c80-7614-11eb-1400-e55b390c3437
begin
  	train = CSV.File("../../ML_UW_Spec/C02/data/kc_house_data_small_train.csv";
		header=true) |> DataFrame;
  	valid = CSV.File("../../ML_UW_Spec/C02/data/kc_house_data_small_validation.csv";
		header=true) |> DataFrame;
	test = CSV.File("../../ML_UW_Spec/C02/data/kc_house_data_small_test.csv";
		header=true) |> DataFrame;

	size(train), size(valid), size(test)
end

# ╔═╡ e0ad9bee-7615-11eb-0226-617fef643e14
md"""
We need `get_data()` and `normalize_features()` from previous week. Let's import them from our utils module.
"""

# ╔═╡ f1d25e60-7614-11eb-2db3-69a5c37b2856
md"""
### Extract features and normalize

Using all of the numerical inputs listed in feature_list, transform the training, test, and validation DataFrames into Julia arrays:
"""

# ╔═╡ f1a1aa7c-7614-11eb-24b8-3183535eb4c1
md"""
In computing distances, it is crucial to normalize features. Otherwise, for example, the sqft_living feature (typically on the order of thousands) would exert a much larger influence on distance than the bedrooms feature (typically on the order of ones). We divide each column of the training feature matrix by its 2-norm, so that the transformed column has unit norm.

*IMPORTANT*: Make sure to store the norms of the features in the training set. The features in the test and validation sets must be divided by these same norms, so that the training, test, and validation sets are normalized consistently.

"""

# ╔═╡ f1b90668-7614-11eb-1313-0b2a3b91b243
begin
	feature_list = [:bedrooms, :bathrooms,
                :sqft_living, :sqft_lot,
                :floors, :waterfront,
                :view, :condition,
                :grade, :sqft_above,
                :sqft_basement, :yr_built,
                :yr_renovated, :lat, :long,
                :sqft_living15, :sqft_lot15]

	features_train, output_train = get_data(train, feature_list, :price)
	features_test, output_test = get_data(test, feature_list, :price)
	features_valid, output_valid = get_data(valid, feature_list, :price)

	## Normalize
	features_train, norms = normalize_features(features_train) # norm. train set features (columns)
	features_test = features_test ./ norms    # normalize test set ...
	features_valid = features_valid ./ norms  # ... and valid set with training set norms

	size(features_train), size(features_test), size(features_valid)
end

# ╔═╡ 0249baea-761a-11eb-09fa-efa1de647bb8
md"""
### Compute a single distance

To start, let's just explore computing the "distance" between two given houses. We will take our query house to be the first house of the test set and look at the distance between this house and the 10th house of the training set.

To see the features associated with the query house, print the first row (index 1) of the test feature matrix. You should get an 18-dimensional vector whose components (in absolute value) are between 0 and 1.

"""

# ╔═╡ 022e64ac-761a-11eb-1daf-d9ad5926c387
begin
	qry_house = features_test[1, :]

	@test size(qry_house)[1] == 18
	@test all(x -> 0. ≤ abs(x) ≤ 1., qry_house)
end

# ╔═╡ 020058f8-761a-11eb-0b03-1be1e9780d1c
md"""
**Quiz Question**

What is the Euclidean distance between the query house and the 10th house of the training set?

"""

# ╔═╡ f15192e4-7614-11eb-1fca-3b6fc4634159
begin
	function calc_euclidean_distance(row::Vector, row_ref::Vector)
  		√(sum((row - row_ref) .* (row - row_ref)))
	end

	function calc_euclidean_distance(rows::Matrix, row_ref::Vector)
		sqrt.(sum(((rows .- row_ref').^2), dims=2))
	end
end

# ╔═╡ f137c3e6-7614-11eb-190e-4d331f41b02e
begin
	d10th = calc_euclidean_distance(features_train[10, :], qry_house)

	with_terminal() do
		@printf("distance %1.8f / rounded: %1.3f\n", d10th, d10th)
	end
end

# ╔═╡ 25250e4c-761b-11eb-2c45-f54d7410184a
md"""
### Compute multiple distances

Of course, to do nearest neighbor regression, we need to compute the distance between our query house and all houses in the training set.

To visualize this nearest-neighbor search, let's first compute the distance from our query house (features_test[1]) to the first 10 houses of the training set (features_train[1:10]) and then search for the nearest neighbor within this small set of houses. Through restricting ourselves to a small set of houses to begin with, we can visually scan the list of 10 distances to verify that our code for finding the nearest neighbor is working.

Compute the Euclidean distance from the query house to each of the first 10 houses in the training set.
"""

# ╔═╡ 67931b0c-761b-11eb-3f3c-cfa948d07ebf
dist_10 = calc_euclidean_distance(features_train[1:10, :], qry_house)

# ╔═╡ e504f550-761f-11eb-1471-0b3ae40edef8
sort(collect(enumerate(dist_10[:,])),
	by=((_ix, x)=t) -> x, rev=false)[1:2]

# ╔═╡ 67441dd6-761b-11eb-1008-9f7fe69785b8
md"""
**Quiz Question**

Among the first 10 training houses, which house is the closest to the query house?
"""

# ╔═╡ 6699cfcc-761b-11eb-0c6a-89637d586437
## Resp.
begin
	ix_10 = argmin(dist_10)
	d_10 = dist_10[ix_10]
	(index=ix_10[1], distance=d_10)
end

# ╔═╡ 667cb82c-761b-11eb-1925-153041dc9a6a
md"""
### Perform k-nearest neighbor regression

For k-nearest neighbors, we need to find a set of k houses in the training set closest to a given query house. We then make predictions based on these k nearest neighbors.
Fetch k-nearest neighbors

Using the functions above, implement a function that takes in
  - the value of k;
  - the feature matrix for the training houses; and
  - the feature vector of the query house

and returns the indices of the k closest training houses. For instance, with 2-nearest neighbor, a return value of [5, 10] would indicate that the 5th and 10th training houses are closest to the query house.
"""

# ╔═╡ 6661ff46-761b-11eb-2817-0bf3e5c48275
function k_nearest(k::Int, feature_matrix::Matrix, feature_vect::Vector)
	@assert k ≥ 1
  	dist = calc_euclidean_distance(feature_matrix, feature_vect)
	sort(collect(enumerate(dist[:])),
		by=((_ix, x)=t) -> x, rev=false)[1:k] |>
		a -> map(t -> t[1], a)
end

# ╔═╡ 6645eb26-761b-11eb-2de2-212f9a229bef
md"""
**Quiz Question**

Take the query house to be third house of the test set (features_test[3]). What are the indices of the 4 training houses closest to the query house?
"""

# ╔═╡ 6628a162-761b-11eb-386d-5d6c4b9d8646
## for quiz, index answers are expected to be 0-based
## Julia however is 1-based index, thus removed 1 to this answers
k_nearest(4, features_train, features_test[3, :])

# ╔═╡ 660da98c-761b-11eb-34af-b3b35c66288c
md"""
### Make a single prediction by averaging k nearest neighbor outputs

Now that we know how to find the k-nearest neighbors, write a function that predicts the value of a given query house. For simplicity, take the average of the prices of the k nearest neighbors in the training set. The function should have the following parameters:

  - the value of k;
  - the feature matrix for the training houses;
  - the output values (prices) of the training houses; and
  - the feature vector of the query house, whose price we are predicting.

The function should return a predicted value of the query house.

"""

# ╔═╡ f7a44f48-7622-11eb-1d99-e33a5cbe3f29
function predict(k::Int, feature_matrix::Matrix, feature_vect::Vector,
		output::Vector)
  ixes = k_nearest(k, feature_matrix, feature_vect)
  sum(output[ixes]) / k
end

# ╔═╡ f7809ec2-7622-11eb-2d16-615ec1f7fb49
md"""
### Quiz Question

Again taking the query house to be third house of the test set (features_test[3]), predict the value of the query house using k-nearest neighbors with k=4 and the simple averaging method described and implemented above.
"""

# ╔═╡ f7637df6-7622-11eb-076f-0f6ab652bfd2
predict(4, features_train, features_test[3, :], output_train)
# cmp
# 413987.5 (k-nearest) vs 249000 (1-nearest)

# ╔═╡ f74b6982-7622-11eb-37d2-21ed0ec6551d
md"""
### Make multiple predictions

Write a function to predict the value of each and every house in a query set. (The query set can be any subset of the dataset, be it the test set or validation set.) The idea is to have a loop where we take each house in the query set as the query house and make a prediction for that specific house. The new function should take the following parameters:

  - the value of k;
  - the feature matrix for the training houses;
  - the output values (prices) of the training houses; and
  - the feature matrix for the query set.

The function should return a set of predicted values, one for each house in the query set.

"""

# ╔═╡ f728dac0-7622-11eb-094f-6782145366d8
function predict_all(k::Int, feature_matrix::Matrix, feature_matrix_qs::Matrix,
		output::Vector)
	n = size(feature_matrix_qs)[1]
	vresp = zeros(Float64, n)
	for ix ∈ 1:n
		vresp[ix] = predict(k, feature_matrix, feature_matrix_qs[ix, :], output)
	end
	vresp
end

# ╔═╡ 5b3578e0-7626-11eb-1c95-27a497721fd1
size(features_test), typeof(features_test)

# ╔═╡ c1c949fc-7625-11eb-212e-7316c164b585
size(features_test[1:10, :]), typeof(features_test[1:10, :])

# ╔═╡ a307a67c-7624-11eb-0021-09f621585464
size(features_test[1:10, :])[1]

# ╔═╡ c51c48a4-7628-11eb-07cd-97f29b9fb5fa
output_train

# ╔═╡ 65f46ecc-761b-11eb-1617-0fe428d5a1bd
md"""
**Quiz Question**

Make predictions for the first 10 houses in the test set using k-nearest neighbors with k=10.

  1. What is the index of the house in this query set that has the lowest predicted value?
  2. What is the predicted value of this house?

"""

# ╔═╡ d4f46ee6-7625-11eb-3b88-3d68b2e1c58b
features_test[2, :], size(features_test[2, :])

# ╔═╡ 65dba81a-761b-11eb-3bb8-47bf75fafa40
begin
	sol = predict_all(10, features_train, features_test[1:10, :], output_train)
	m_ix = argmin(sol)
	with_terminal() do
		println((solution=sol, index=m_ix[1], predicted_price=sol[m_ix]))
	end
end

# ╔═╡ 65bfee54-761b-11eb-16e8-4b4ca6decc9d
md"""
### Choosing the best value of k using a validation set

There remains a question of choosing the value of k to use in making predictions. Here, we use a validation set to choose this value. Write a loop that does the following:

```
  For k ∈ [1, 2, ..., 15]
    Makes predictions for each house in the validation set using the k-nearest neighbors from the training set.
    Computes the RSS for these predictions on the validation set
    Stores the RSS computed above in rss_all
  End
  Report which k produced the lowest RSS on validation set.
```

"""

# ╔═╡ 9cb7ace4-762a-11eb-3e80-f16cbf712a62
function rss_fn(preds, output)
	sum((preds .- output) .^ 2)
end

# ╔═╡ 726f0674-7629-11eb-2c27-7d634b6d1d4b
function find_best_k(features_matrix::Matrix, features_matrix_qs::Matrix,
		output::Vector, output_val::Vector; range=1:15)
 	rss_all = []
  	best_k, best_rss = nothing, nothing
  	for k ∈ range
    	preds = predict_all(k, features_matrix, features_matrix_qs, output)
    	rss = rss_fn(preds, output_val)
    	push!(rss_all, rss)
    	if isnothing(best_rss) || rss < best_rss
			best_rss, best_k = rss, k
		end
	end
  	(best_k, best_rss, rss_all)
end

# ╔═╡ 65a2c806-761b-11eb-2f1c-3773d2584857
(best_k, best_rss, rss_all) = find_best_k(features_train, features_valid, output_train, output_valid)

# ╔═╡ 6577f8ec-761b-11eb-01f5-67cb002a7070
md"""
To visualize the performance as a function of k, plot the RSS on the *validation* set for each considered k value.
"""

# ╔═╡ 377da4be-762a-11eb-31be-c1041f7c306b
plot(collect(1:15), rss_all, legend=false, color=[:lightblue], marker=".")

# ╔═╡ 754a565c-762a-11eb-0b81-9d30ddf27406
md"""
**Quiz Question**

What is the RSS on the test data using the value of k found above? To be clear, sum over all houses in the test set.

"""

# ╔═╡ 9210589a-762a-11eb-3439-f3130e570162
begin
	test_preds = predict_all(best_k, features_train, features_test, output_train)
	test_rss = rss_fn(test_preds, output_test)

	with_terminal() do
		@printf("test rss: %1.7e\n", test_rss)
	end
end

# ╔═╡ Cell order:
# ╟─84a95dde-7614-11eb-20da-8191d2c5103e
# ╠═cbcca572-7614-11eb-2fe5-699acc8e2746
# ╟─db51098e-7614-11eb-2dcd-57f4526747f6
# ╠═f20a9c80-7614-11eb-1400-e55b390c3437
# ╟─e0ad9bee-7615-11eb-0226-617fef643e14
# ╠═f1eefcb4-7614-11eb-0173-d5112c9036d7
# ╟─f1d25e60-7614-11eb-2db3-69a5c37b2856
# ╟─f1a1aa7c-7614-11eb-24b8-3183535eb4c1
# ╠═f1b90668-7614-11eb-1313-0b2a3b91b243
# ╟─0249baea-761a-11eb-09fa-efa1de647bb8
# ╠═022e64ac-761a-11eb-1daf-d9ad5926c387
# ╟─020058f8-761a-11eb-0b03-1be1e9780d1c
# ╠═f15192e4-7614-11eb-1fca-3b6fc4634159
# ╠═f137c3e6-7614-11eb-190e-4d331f41b02e
# ╟─25250e4c-761b-11eb-2c45-f54d7410184a
# ╠═67931b0c-761b-11eb-3f3c-cfa948d07ebf
# ╠═e504f550-761f-11eb-1471-0b3ae40edef8
# ╟─67441dd6-761b-11eb-1008-9f7fe69785b8
# ╠═6699cfcc-761b-11eb-0c6a-89637d586437
# ╟─667cb82c-761b-11eb-1925-153041dc9a6a
# ╠═6661ff46-761b-11eb-2817-0bf3e5c48275
# ╟─6645eb26-761b-11eb-2de2-212f9a229bef
# ╠═6628a162-761b-11eb-386d-5d6c4b9d8646
# ╟─660da98c-761b-11eb-34af-b3b35c66288c
# ╠═f7a44f48-7622-11eb-1d99-e33a5cbe3f29
# ╟─f7809ec2-7622-11eb-2d16-615ec1f7fb49
# ╠═f7637df6-7622-11eb-076f-0f6ab652bfd2
# ╟─f74b6982-7622-11eb-37d2-21ed0ec6551d
# ╠═f728dac0-7622-11eb-094f-6782145366d8
# ╠═5b3578e0-7626-11eb-1c95-27a497721fd1
# ╠═c1c949fc-7625-11eb-212e-7316c164b585
# ╠═a307a67c-7624-11eb-0021-09f621585464
# ╠═c51c48a4-7628-11eb-07cd-97f29b9fb5fa
# ╟─65f46ecc-761b-11eb-1617-0fe428d5a1bd
# ╠═d4f46ee6-7625-11eb-3b88-3d68b2e1c58b
# ╠═65dba81a-761b-11eb-3bb8-47bf75fafa40
# ╟─65bfee54-761b-11eb-16e8-4b4ca6decc9d
# ╠═9cb7ace4-762a-11eb-3e80-f16cbf712a62
# ╠═726f0674-7629-11eb-2c27-7d634b6d1d4b
# ╠═65a2c806-761b-11eb-2f1c-3773d2584857
# ╟─6577f8ec-761b-11eb-01f5-67cb002a7070
# ╠═377da4be-762a-11eb-31be-c1041f7c306b
# ╟─754a565c-762a-11eb-0b81-9d30ddf27406
# ╠═9210589a-762a-11eb-3439-f3130e570162

ext: jl


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, SpecialFunctions, Test, LazyArrays
import ApproxFunBase: Multiplication, testraggedbelowoperator, testbandedoperator, interlace, ∞, ℵ₀

@testset "ODE" begin
    @testset "Airy" begin
        d=Interval(-10.,5.);
        S=Chebyshev(d)

        Bm=Evaluation(d,leftendpoint(d));
        Bp=Evaluation(d,rightendpoint(d));
        B=[Bm;Bp];
        D2=Derivative(d,2);
        X=Multiplication(Fun(x->x,d));

        testbandedoperator(D2-X)
        testraggedbelowoperator([B;D2-X])

        @time u = [B;D2-X] \ [airyai.(endpoints(d))...,0.];
        @test Number.(Array(B*u)) ≈ [airyai.(endpoints(d))...]

        @test ≈(u(0.),airyai(0.);atol=10ncoefficients(u)*eps())

        @time u=[Bm;D2-X;Bp]\[airyai(leftendpoint(d)),0.,airyai(rightendpoint(d))];
        @test ≈(u(0.),airyai(0.);atol=10ncoefficients(u)*eps())

        @time u=[D2-X;Bm;Bp]\[0.,airyai(leftendpoint(d)),airyai(rightendpoint(d))];
        @test ≈(u(0.),airyai(0.);atol=10ncoefficients(u)*eps())

        d=Interval(-1000.,5.);
        Bm=Evaluation(d,leftendpoint(d));
        Bp=Evaluation(d,rightendpoint(d));
        B=[Bm;Bp];
        D2=Derivative(d,2);
        X=Multiplication(Fun(x->x,d));

        u=[B;D2-X]\[airyai(leftendpoint(d)),airyai(rightendpoint(d)),0.];
        @test ≈(u(0.),airyai(0.);atol=10ncoefficients(u)*eps())

        B=Neumann(d);
        A=[B;D2-X];
        b=[[airyaiprime(leftendpoint(d)),airyaiprime(rightendpoint(d))],0.];

        @time u=A\b;

        @test ≈(u(0.),airyai(0.);atol=10ncoefficients(u)*eps())

        ##) Neumann condition
    end

    @testset "Exp" begin
        f=Fun(x->-x^2)
        g=Fun(t->exp(-t^2))

        fp=f';
        Bm=Evaluation(domain(f),leftendpoint(domain(f)));
        u=[Bm,Derivative(domain(f)) - fp]\[exp(f(leftendpoint(domain(f)))),0.];
        @test norm(u-g)<100eps()
    end

    @testset "Oscillatory integral" begin
        f=Fun(exp);
        D=Derivative(domain(f));
        w=10.;
        B=ApproxFunBase.SpaceOperator(BasisFunctional(floor(w)),Chebyshev(),ApproxFunBase.ConstantSpace(Float64));
        A=[B; D+1im*w*I];

        @time u = A\[0.,f];
        @test (u(1.)exp(1im*w)-u(-1.)exp(-1im*w)) ≈ (-0.18575766879136255 + 0.17863980562549928im)
    end

    @testset "Bessel" begin
        d=ChebyshevInterval()
        D=Derivative(d)
        x=Fun(identity,d)
        A=x^2*D^2+x*D+x^2
        testbandedoperator(x^2*D^2)
        testbandedoperator(ToeplitzOperator([0.5],[0.0,0.5]))
        testbandedoperator(HankelOperator(Float64[]))
        testbandedoperator(A)
        u=[ldirichlet(d);A]\[besselj(0,leftendpoint(d)),0.];

        @test u(0.1) ≈ besselj(0.,0.1)
        @test norm(A*u)<10eps()
        @test norm(Fun(A.ops[1]*u,d)-x.^2 .* differentiate(u,2))<eps()
        @test norm(Fun(A.ops[2]*u,d)-x.*u') < eps()
        @test norm(Fun(A.ops[end]*u,d)-x.^2 .* u) < eps()
        @test norm(x.^2 .*u'' + x.*u' + x.^2 .* u)<10eps()
    end

    @testset "QR" begin
        S=Chebyshev()
        B=Dirichlet(S)
        D=Derivative(S)

        Q,R=qr([B;D^2+I])
        @test Q[1,1] == -0.5773502691896257
        @test size(Q') == (ℵ₀,ℵ₀)
        u=R\(Q'*[[cos(-1.0),cos(1.0)],0.0])

        @test u(0.) ≈ cos(0.0)

        S=Chebyshev()
        A=[Dirichlet(S);Derivative(S)^2 - I]
        QRf = qr(A)
        @test (QRf\[[1.,0],0])(0.0) ≈ 0.3240271368319427
        Q,R = qr(A)
        u=(R\(Q'*[[1.,0.0],0.0]))
        @test u(0.0)  ≈ 0.3240271368319427
        # matrix RHS
        U = QRf \ [[[1.,0.],0] [[0.,1.0],0.]]
        @test U[1,1] ≈ u

        x=Fun(S)
        A=[Dirichlet(S); Derivative(S)^2 - exp(im*x)]
        QRf = qr(A)

        u=(QRf\[[1.,0.0],0.0])
        @test u(0.0) ≈ (0.3329522068795961 + 0.024616008954634165im)
    end

    @testset "Union of intervals" begin
        x=Fun(identity,(-2..15) \ [-1,0])
        sp=space(x)

        B = [Dirichlet(sp);continuity(sp,0:1)]

        # We don't want to concat piecewise space
        @test !(continuity(sp,0) isa ApproxFunBase.VectorInterlaceOperator)
        @test B isa ApproxFunBase.VectorInterlaceOperator

        D=Derivative(sp)
        A=[B;D^2-x]

        ApproxFunBase.testraggedbelowoperator(A)
        QR=qr(A)

        @time u=QR\[[airyai(-2.),0.0],zeros(4),0.0]

        @test u(0.0) ≈ airyai(0.)
    end

    @testset "Vector" begin
        d=ChebyshevInterval()
        D=Derivative(d);
        B=ldirichlet();
        Bn=lneumann();

        f=Fun(x->[exp(x),cos(x)],d)

        A=[B 0;
           Bn 0;
           0 B;
           D^2-I 2.0I;
           0 D+I]

        # makes sure ops are in right order
        @test A.ops[4,1] isa ApproxFunBase.PlusOperator
        QR=qr(A)
        v=Any[0.,0.,0.,f...]
        @test (QR\v)(0.0) ≈ [0.0826967758420519,0.5553968826533497]

        Q,R=qr(A)
        v=Any[0.,0.,0.,f...]
        @test (QR\v)(0.0) ≈ [0.0826967758420519,0.5553968826533497]
    end

    @testset "Auto-space" begin
        t=Fun(identity,0..1000)
        D = Derivative()
        L=D^2+2I  # our differential operator, 𝒟 is equivalent to Derivative()

        u=[ivp();L]\[0.;0.;cos(100t)]
        @test ≈(u(1000.0),0.00018788162639452911;atol=1000eps())

        x=Fun(identity,1..2000)
        d=domain(x)
        B=Dirichlet()
        ν=100.
        L=(x^2*D^2) + x*D + (x^2 - ν^2)   # our differential operator

        @time u=[B;L]\[[besselj(ν,leftendpoint(d)),besselj(ν,rightendpoint(d))],0.]
        @test ≈(u(1900.),besselj(ν,1900.);atol=1000eps())
    end

    @testset "Complex RHS for real operator" begin
        D=Derivative(Chebyshev())
        B=ldirichlet()
        u1=[B;D]\[0.;Fun(exp)+0im]
        u2=[B;D]\[0.;Fun(exp)]
        @test u1(0.1) ≈ u2(0.1)
    end

    @testset "Delay" begin
        S = Chebyshev(0..1)
        D = Derivative(S)
        L = [D            0I      0I;
             I             D      0I;
             0I           I        D]

        B = [ldirichlet(S)   0                 0;
             rdirichlet(S)  -ldirichlet(S)     0;
             0              rdirichlet(S)      -ldirichlet(S)]

        v = [B; L] \ [2.3; zeros(2); zeros(3)]

        u = v[1] + setdomain(v[2], Domain(1..2)) + setdomain(v[3], Domain(2..3))

        @test abs(u'(1.3) + u(1.3-1) ) ≤ 10eps()
    end

    @testset "Tenth order" begin
        x = Fun()
        D = Derivative()

        L=D^10 + cosh(x)*D^8 + x^3*D^6 + x^4*D^4 + cos(x)*D^2 + x^2
        d = ChebyshevInterval()
        B = [Dirichlet(d) ;
             Neumann(d)   ;
             [Evaluation(d,leftendpoint,k) for k=2:4]... ;
             [Evaluation(d,rightendpoint,k) for k=2:4]...]

        rs = rangespace(B)
        @test ApproxFunBase.blocklengths(rs) == [10]

        A = [B; L]
        rs = rangespace(A)
        # @test ApproxFunBase.blocklengths(rs) isa Vcat

        u = [B; L] \ [ [0.,0.], [1.,1.], zeros(6)..., exp(x)]
        @test u(0.5) ≈ -0.4024723414410859 # Empirical
    end

    @testset "constant forcing" begin
        d = Interval(0.,50.)
        D = Derivative(d)
        t = Fun(identity,d)

        F = D^2 +.5D + I

        A= [ 0    ldirichlet(d);
            0    lneumann(d);
            0    rdirichlet(d);
            -1    F; ]

        @time u,x=A\[1.,0.,2.,0.]

        @test norm(F*x-u)<1000eps()
    end
end


text:
struct LightSource
    n::Integer # number of atoms
    Em::Vector # line magnitudes
    νm::Vector # central frequencies of lines in GHz
    σ::Number # Doppler broadening in GHz
    νMin::Number # minimum bandpass frequency in GHz
    νMax::Number # maximum bandpass frequency in GHz
    γRate::Number # photon count rate in GHz
    function LightSource(n::Integer,Em::Vector,νm::Vector,σ::Number,νMin::Number,νMax::Number,γRate::Number)
        @assert length(Em) == length(νm) "Vectors for line magnitudes and frequencies must have the same length"
        @assert νMin < νMax "Bandpass minimum must be less than maximum"
        new(
            n,
            convert(Vector{Complex},Em),
            convert(Vector{Real},νm),
            convert(Real,σ),
            convert(Real,νMin),
            convert(Real,νMax),
            convert(Real,γRate)
        )
    end
end

function LightSource(n::Integer,Em::Vector,νm::Vector,σ::Number, γRate::Number)
    νMin = min(νm...) - 5*σ # automatically set the bandpass minimum to 5σ below the lowest line frequency
    νMax = max(νm...) + 5*σ # automatically set the bandpass maximum to 5σ above the highest line frequency
    return LightSource(n,Em,νm,σ,νMin,νMax,γRate)
end

function LightSource(params::Dict)
    return LightSource(
        params[:n],
        params[:Em],
        params[:νm],
        params[:σ],
        params[:fγ]
    )
end

export LightSource

"""
    ν0(source::LightSource)

Returns the average frequency of spectral lines in source
"""
function ν0(source::LightSource)
    return mean(source.νm)
end

export ν0

"""
    Δm(source::LightSource)

Returns vector of the sources line separation from the average frequency of all lines
"""
function Δm(source::LightSource)
    return source.νm .- ν0(source)
end

export Δm

"""
    lineShifts(source::LightSource)

Returns a vector with the frequency differences between all spectral lines
"""
function lineShifts(source::LightSource)
    return map(x->abs(x[2]-x[1]),subsets(source.νm,2))
end

export lineShifts

"""
    nbar(t::Real,source::LightSource)

Returns the average photon counts received over a duration t from source
"""
function nbar(t::Real,source::LightSource)
    return t*source.γRate
end

export nbar

################################################################################
# Electric field
################################################################################
"""
    eField(νn::Vector,ϕmn::Matrix,source::LightSource)

Container holding frequencies and phases for one realization of the electric field.
"""
struct eField
    νn::Vector
    ϕmn::Matrix
    source::LightSource

    function eField(νn::Vector,ϕmn::Matrix,source::LightSource)
        @assert (size(ϕmn)[1] == length(source.νm) && size(ϕmn)[2] == length(νn)) "Phase array must have shape (m,n)"
        new(νn,ϕmn,source)
    end
end

"""
    eField(source::LightSource, seed::Integer = -1)

Generate frequencies and phases for a single instance of the electric field.
"""
function eField(source::LightSource, seed::Integer = -1)
    if seed != -1
       Random.seed!(seed)
    end
    νDist = Normal(ν0(source),source.σ)
    νn    = rand(νDist,source.n)
    ϕmn   = 2*π*rand(Float64,(length(source.νm),source.n))
    return eField(νn,ϕmn,source)
end

export eField

"""
	function eField(t::Number,field::eField)

Returns the electric field value at time t(ns)
"""
function eFieldT(t::Number,field::eField)
	# generate frequencies
	ωmn = transpose(field.νn) .+ Δm(field.source)
    ωmn .*= 2*π
	# add the phase
	exponentmn = -im*(t*ωmn+field.ϕmn)
	# put them in the exponent
	enm = exp.(exponentmn)
	# multiply by the field magnitude
	fieldnm = field.source.Em .* enm
	# add it all together
	return sum(ivec(fieldnm))
end

export eFieldT

################################################################################
# Calculate intensity
################################################################################
"""
    intensity(et::Number)

Returns the intensity given the value of the electric field at a particular point in time and space
"""
function intensity(et::Number)
    return real(et*conj(et))
end

"""
    intensity(t::Real,field::eField)

Returns the intensity given an instance of the EM field and a time t(ns).
"""
function intensity(t::Number,field::eField)
    efield = eFieldT(t,field)
    return intensity(efield)
end

export intensity

"""
    γIntensity(nbar::T,γvec::Vector{T}) where T<:Real

Returns γIntensity object which contains the average photon count rate and the renormalized intensity time series
"""
struct γIntensity{T<:Real}
    nbar::T
    γvec::Vector{T}
    function γIntensity(nbar::T,γvec::Vector{T}) where T<:Real
        total = sum(γvec)
        coeff = nbar/total
        return new{T}(nbar,coeff*γvec)
    end
end

function γIntensity(nbar::Real,t::Real,dt::Real,source::LightSource)
    field = eField(source)
    return γIntensity(nbar,t,dt,field)
end

function γIntensity(nbar::Real,t::Real,dt::Real,field::eField)
    intensitySeries = map(time->intensity(time,field),0:dt:t)
    return γIntensity(nbar,intensitySeries)
end

export γIntensity

# array indexing and iteration interface for γIntensity object
function Base.length(γint::γIntensity)
    return length(γint.γvec)
end

function Base.getindex(γint::γIntensity,i::Int)
    return γint.γvec[i]
end

function Base.firstindex(γint::γIntensity)
    return γint.γvec[1]
end

function Base.lastindex(γint::γIntensity)
    return γint.γvec[end]
end

function Base.iterate(γint::γIntensity)
    return (Base.firstindex(γint),1)
end

function Base.iterate(γint::γIntensity,state::Int)
    return state >= length(γint) ? nothing : (γint[state+1],state+1)
end

"""
    Beamsplitter(r::Number,t::Number)

Struct to hold beam splitter coefficients.
"""
struct Beamsplitter
    r::Number
    t::Number
    function Beamsplitter(r::Number,t::Number)
        beamNorm = sqrt(r^2+t^2)
        new(r/beamNorm,t/beamNorm)
    end
end

export Beamsplitter


text:
using Test

using CfGRIB
using Conda

if !(haskey(Conda._installed_packages_dict(), "cfgrib"))
    Conda.add_channel("conda-forge")
    Conda.add("cfgrib")
end

const dir_tests = abspath(joinpath(dirname(pathof(CfGRIB)), "..", "test"))
const dir_testfiles = abspath(joinpath(dir_tests, "sample-data"))

@testset "CfGRIB.jl" begin
    @testset "Indexing" begin
        include("indexing.jl")
    end

    @testset "CFMessage" begin
        include("cfmessage.jl")
    end

    @testset "Dataset" begin
        include("dataset.jl")
    end

    @testset "Parity" begin
        include("parity.jl")
    end

    @testset "Backends" begin
        include("backends.jl")
    end
end


text:
```
Declination(args...)

Estimate decliantion of a given data set. By default it is applied an Anomaly Dectection
algorithm to filter the representative data.

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Date|Mandatory|--|Array{Date,1}| Array of dates|
|Rate|Mandatory|--|Array{Number,1}| Array of rates|
|AD|Optional|AD=true|AD=Bool| Activation of Anomaly detection algorithm|
```
function Declination(time,rate; AD=true)

        #Anomally Detection
if AD==true
    epsilon=0.05339
    mu=mean(rate)
    sigma=std(rate)
    NormRate=map(x->(x-mu)/sigma,rate)
    pval=map(x->pdf(Normal(),x),NormRate)
    ev=pval.>epsilon
    Time=time[ev]
    Rate=rate[ev]
else
    Time=time
    Rate=rate
end

        #
        Ti=Dates.days.(Time[1])
        Cum=accumulate(+,Rate)

        D=Reg(Dates.days.(Time).-Ti,Rate, TypeReg="Exponential")
        Q=Reg(Cum,Rate, TypeReg="Linear")

        S=ExponentialDec(D[1],Q[2],Ti)
        DecLine=map(X->S.Qi.*exp(S.Di.*X),Dates.days.(Time).-S.Ti)

            return S, Time, DecLine

end #End Function

function DeclinationForecast(S,DateRange)

    ForecastLine=map(X->S.Qi.*exp(S.Di.*X),Dates.days.(DateRange).-S.Ti)

end # End function


text:
using LogisticRegression
using Test

include("simple_models.jl")


text:
"""
    grm(s, method=:GRM, minmaf=0.01, colinds=nothing)

Compute empirical kinship matrix from a SnpArray `s`. Missing genotypes are imputed
on the fly by column mean.

# Arguments
- `method::Symbol=:GRM`: `:GRM` (default), `:MoM`, or `:Robust`.
- `minmaf::Real=0.01`: columns with MAF less than `minmaf` are excluded; default 0.01.
- `cinds=nothing`: indices or mask of columns to be used for calculating GRM.
- `t::Type{T}=Float64`: Float type for calculating GRM; default is `Float64`.
"""
function grm(
    s::SnpArray;
    method::Symbol = :GRM,
    minmaf::Real = 0.01,
    cinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    t::Type{T} = Float64
    ) where T <: AbstractFloat
    mf = maf(s)
    colinds = something(cinds, mf .≥ minmaf)
    @views grm(s[:, colinds], mf[colinds], method, t)
end # function grm

function grm(
    s::AbstractSnpArray,
    maf::AbstractVector,
    method::Symbol = :GRM,
    ::Type{T} = Float64
    ) where T <: AbstractFloat
    m, n = size(s)
    G = Mmap.mmap(Matrix{T}, m, n) # slightly faster than G = Matrix{T}(undef, m, n)
    if method == :GRM
        Base.copyto!(G, s, model=ADDITIVE_MODEL, impute=true, center=true, scale=true)
        Φ = G * transpose(G)
        Φ ./= 2n
    elseif method == :MoM
        Base.copyto!(G, s, model=ADDITIVE_MODEL, impute=true)
        G .-= 1
        Φ = G * transpose(G)
        c = sum(x -> abs2(x) + abs2(1 - x), maf)
        shft, scal = n / 2 - c, 1 / (n - c)
        @inbounds @simd for i in eachindex(Φ)
            Φ[i] = (Φ[i] / 2 + shft) * scal
        end
    elseif method == :Robust
        Base.copyto!(G, s, model=ADDITIVE_MODEL, center=true, impute=true)
        scal = sum(x -> 4x * (1 - x), maf)
        Φ = G * transpose(G)
        Φ ./= scal
    else
        throw(ArgumentError("method should be :GRM, :MoM, or :Robust"))
    end
    Φ
end # function grm


text:
#Add structure/type that defined actuators and their effects

#Define AbstractActuator type
abstract type AbstractActuator end

#= Define a Realstic Thruster type which includes:
    - Specific impulse (only used when infinitely throttable, to determine fuel used)
                        (set to inf if no fuel is to be used)
    - [x,y,z] position w.r.t. centre of mass of vehicle
    - Transformation matrix between body-axis of vehicle and thruster (thrust in x-direction)
    - NozzleType;
    - TypeFlag; indicates whether or not the thruster is throttable (0 = on/off, 1 = throttable)
                (if this flag is zero, a non-zero thrust setting is interpreted as on)
    - ThrottableFlag; indicates if actuator is finely throttable (0) or has limits (1)
                     if it has limits, the nozzle type will be used to determine thrust
                     as function of mass flow and pressure. If its infinitely throttable,
                     the force generated is the control variable, else its the mass flow.
    - ConstantFlag; indicates whether or not the thruster generates a constant thrust
                    (0 = constant, 1 = variable)
                    (if the flag is zero, the specific impulse is constant
                    (if the flag is one, the altitude function will be used to compute the
                    atmospheric pressure that determines the variation in generated thrust)
    -ControlGroupFlag; indicates if used for translational (0) or rotational control (1)
    - Force generated; Computed
    - Moment generated; Computed

=#
struct RealisticThruster <: AbstractActuator

end #RealisticThruster

#allow for simplistic thruster models
struct SimpleThruster <: AbstractActuator
    ForceMoment::Vector{Float64} #Vector in body-axis frame
    TypeFlag::Bool #0 means that it generates a force, 1 means it generates a moment
end #SimpleThruster

#Define an AerodynamicEffector type
struct AerodynamicEffector <: AbstractActuator

end #AerodynamicEffector


text:
using FlightMechanics.Models

@testset "Dynamic systems conversions" begin
    pos = EarthPosition(1000.0, 500.0, -1000.0)
    att = Attitude(π/4, π/32, π/8)
    vel = [100.0, 0.0, 0.0]
    ang_vel = [π/90, π/18, π/9]
    accel = [1.5, 0.3, 0.5]
    ang_accel = [π/45, π/60, π/30]
    state = State(pos, att, vel, ang_vel, accel, ang_accel)

    # SixDOFEulerFixedMass
    x = convert(SixDOFEulerFixedMass, state)
    state_ = convert(state, x)
    @test isapprox(state, state_)

    # SixDOFQuaternionFixedMass
    x = convert(SixDOFQuaternionFixedMass, state)
    state_ = convert(state, x)
    @test isapprox(state, state_)

    # SixDOFAeroEulerFixedMass
    x = convert(SixDOFAeroEulerFixedMass, state)
    state_ = convert(state, x)
    @test isapprox(state, state_)

    # SixDOFECEFQuaternionFixedMass
    pos_ecef = convert(ECEFPosition, pos)
    state = State(pos_ecef, att, vel, ang_vel, accel, ang_accel)
    x = convert(SixDOFECEFQuaternionFixedMass, state)
    state_ = convert(state, x)
    @test isapprox(state, state_)
end

@testset "Dynamic systems evaluation" begin
    ac = F16()

    att = Attitude(250*DEG2RAD, 20*DEG2RAD, -30*DEG2RAD)
    pos = EarthPosition(0.0, 0.0, 0.0)
    u, v, w = wind2body(500*FT2M, 0, 0, 0.5, -0.2)

    state = State(
        pos,
        att,
        [u, v, w],
        [0.7, -0.8, 0.9],
        [0., 0., 0.],
        [0., 0., 0.]
    )

    # Modify directly the deflexions of the surfaces on the FCS
    fcs = get_fcs(ac)
    set_value!(fcs.de, 20.0*DEG2RAD)
    set_value!(fcs.da, -15.0*DEG2RAD)
    set_value!(fcs.dr, -20*DEG2RAD)
    set_value!(fcs.cpl, 90.0)

    env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")
    aerostate = AeroState(state, env)
    grav = get_gravity(env)

    ac = calculate_aircraft(ac, controls, aerostate, state, grav; consume_fuel=false)
    pfm = ac.pfm
    mass_props = get_mass_props(ac)
    mass = mass_props.mass
    inertia = mass_props.inertia
    h = [160.0*SLUGFT2_2_KGM2, 0.0, 0.0]

    ds_1 = convert(SixDOFEulerFixedMass, state)
    x_1 = get_x(ds_1)
    f = get_state_equation(ds_1)
    x_dot_1 = f(x_1, mass, inertia, pfm.forces, pfm.moments, h)
    ds = SixDOFEulerFixedMass(x_1, x_dot_1)
    st1 = convert(state, ds)

    ds_2 = convert(SixDOFAeroEulerFixedMass, state)
    x_2 = get_x(ds_2)
    f = get_state_equation(ds_2)
    x_dot_2 = f(x_2, mass, inertia, pfm.forces, pfm.moments, h)
    ds = SixDOFAeroEulerFixedMass(x_2, x_dot_2)
    st2 = convert(state, ds)

    ds_3 = convert(SixDOFQuaternionFixedMass, state)
    x_3 = get_x(ds_3)
    f = get_state_equation(ds_3)
    x_dot_3 = f(x_3, mass, inertia, pfm.forces, pfm.moments, h)
    ds = SixDOFQuaternionFixedMass(x_3, x_dot_3)
    st3 = convert(state, ds)

    ds_4 = convert(SixDOFECEFQuaternionFixedMass, state)
    x_4 = get_x(ds_4)
    f = get_state_equation(ds_4)
    x_dot_4 = f(x_4, mass, inertia, pfm.forces, pfm.moments, h)
    ds = SixDOFECEFQuaternionFixedMass(x_4, x_dot_4)
    st4 = convert(state, ds)
    # Note: Transform position to Earth position to be comparable with other states
    st4 = State(
        convert(st1.position, st4.position),
        st4.attitude,
        st4.velocity,
        st4.angular_velocity,
        st4.acceleration,
        st4.angular_acceleration
    )

    @testset "$(pair[1])" for pair in (
            ("SixDOFEulerFixedMass-SixDOFAeroEulerFixedMass", st1, st2),
            ("SixDOFEulerFixedMass-SixDOFQuaternionFixedMass", st1, st3),
            ("SixDOFAeroEulerFixedMass-SixDOFQuaternionFixedMass", st2, st3),
        )
        x, y = pair[2:3]
        @test isapprox(x.position, y.position)
        @test isapprox(x.attitude, y.attitude)
        @test isapprox(x.velocity, y.velocity)
        @test isapprox(x.angular_velocity, y.angular_velocity)
        @test isapprox(x.acceleration, y.acceleration)
        @test isapprox(x.angular_acceleration, y.angular_acceleration)
    end

    @testset "$(pair[1])" for pair in (
            ("SixDOFEulerFixedMass-SixDOFECEFQuaternionFixedMass", st1, st4),
            ("SixDOFAeroEulerFixedMass-SixDOFECEFQuaternionFixedMass", st2, st4),
            ("SixDOFQuaternionFixedMass-SixDOFECEFQuaternionFixedMass", st3, st4),
        )
        x, y = pair[2:3]
        @test isapprox(x.position, y.position)
        @test isapprox(x.attitude, y.attitude)
        @test isapprox(x.velocity, y.velocity)
        @test isapprox(x.angular_velocity, y.angular_velocity)
        # TODO: If ROT_VELOCITY in constants is set to zero this tests passes with atol=0.0
        # Otherwise the systems are not comparable because of the effect of rotating Earth.
        @test isapprox(x.acceleration, y.acceleration, atol=0.01)
        @test_broken isapprox(x.acceleration, y.acceleration)

        @test isapprox(x.angular_acceleration, y.angular_acceleration)
    end

end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct SAAData{T <: AbstractFloat}
    sample_size::Integer = 16
    interval::ConfidenceInterval{T} = ConfidenceInterval(-Inf, Inf, 1.0)
end

@with_kw mutable struct SAAParameters{T <: AbstractFloat}
    confidence::T = 0.95
    num_lower_trials::Integer = 10
    num_upper_trials::Integer = 10
    num_samples::Integer = 100
    num_eval_samples::Integer = 1000
    num_ews_samples::Integer = 1000
    num_eev_samples::Integer = 1000
    tolerance::T = 1e-2
    init_num_samples::Int = 16
    max_num_samples::Int = 5000
    optimize_config::Function = (optimizer, N) -> nothing
    time_limit::T = Inf
    log::Bool = true
    keep::Bool = false
    offset::Int = 5
    indent::Int = 4
end

"""
    SampleAverageApproximation

Functor object for the sample average approximation algorithm.
"""
mutable struct SampleAverageApproximation{T <: AbstractFloat, A <: AbstractVector, SM <: StochasticModel, S <: AbstractSampler}
    model::SM
    sampler::S
    optimizer
    x₀::A
    interval_history::Vector{ConfidenceInterval{T}}
    data::SAAData{T}
    parameters::SAAParameters{T}

    function SampleAverageApproximation(model::StochasticModel, sampler::AbstractSampler, x₀::AbstractVector; optimizer = nothing, kw...)
        T = eltype(x₀)
        A = typeof(x₀)
        SM = typeof(model)
        S = typeof(sampler)
        # Algorithm parameters
        params = SAAParameters{T}(; kw...)
        return new{T,A,SM,S}(model,
                             sampler,
                             optimizer,
                             x₀,
                             Vector{ConfidenceInterval{T}}(),
                             SAAData{T}(),
                             params)
    end
end

function (saa::SampleAverageApproximation)()
    @unpack log, tolerance, confidence = saa.parameters
    @unpack init_num_samples, max_num_samples = saa.parameters
    @unpack num_lower_trials, num_upper_trials, num_eval_samples = saa.parameters
    N = init_num_samples
    α = 1 - confidence
    progress = ProgressThresh(tolerance, 0.0, "SAA gap")
    log && ProgressMeter.update!(progress, Inf,
                                 showvalues = [
                                     ("Confidence interval", NaN),
                                     ("Relative error", Inf),
                                     ("Sample size", NaN),
                                 ])
    status = MOI.OPTIMIZE_NOT_CALLED
    while true
        saa.data.sample_size = N
        CI = confidence_interval(saa.model, saa.sampler)
        saa.data.interval = CI
        push!(saa.interval_history, CI)
        Q = (upper(CI) + lower(CI)) / 2
        gap = length(CI) / abs(Q + 1e-10)
        log && ProgressMeter.update!(progress, gap,
                                     showvalues = [
                                         ("Confidence interval", CI),
                                         ("Relative error", gap),
                                         ("Sample size", N)
                                     ])
        if gap <= tolerance
            status = MOI.OPTIMAL
            break
        end
        N = N * 2
        if N > max_num_samples
            status = MOI.ITERATION_LIMIT
            log && ProgressMeter.update!(progress, 0.0,
                                         showvalues = [
                                             ("Confidence interval", CI),
                                             ("Early termination", status),
                                             ("Sample size", N)
                                         ])
            break
        end
        #optimizer_config(optimizer(stochasticmodel), N)
    end
    # Move cursor down
    for _ in 1:saa.parameters.offset + 2
        print("\r\u1b[B")
    end
    # Return status
    return status
end

function instance(saa::SampleAverageApproximation)
    sp = instantiate(saa.model,
                     saa.sampler,
                     saa.data.sample_size;
                     optimizer = saa.optimizer)
    optimize!(sp)
    Q = objective_value(sp)
    while !(Q ∈ saa.data.interval)
        sp = instantiate(saa.model,
                         saa.sampler,
                         saa.data.sample_size;
                         optimizer = saa.optimizer)
        optimize!(sp)
        Q = objective_value(sp)
    end
    return sp
end


text:
## gradient methods

"""
    gradient(inFLD::MeshArray,Γ::NamedTuple)

Compute spatial derivatives. Other methods:
```
gradient(inFLD::MeshArray,Γ::NamedTuple,doDIV::Bool)
gradient(inFLD::MeshArray,iDXC::MeshArray,iDYC::MeshArray)
```
"""
function gradient(inFLD::MeshArray,Γ::NamedTuple)
(dFLDdx, dFLDdy)=gradient(inFLD,Γ,true)
return dFLDdx, dFLDdy
end

function gradient(inFLD::MeshArray,Γ::NamedTuple,doDIV::Bool)

exFLD=exchange(inFLD,1)
dFLDdx=similar(inFLD)
dFLDdy=similar(inFLD)

for a=1:inFLD.grid.nFaces
  (s1,s2)=size(exFLD.f[a])
  tmpA=view(exFLD.f[a],2:s1-1,2:s2-1)
  tmpB=tmpA-view(exFLD.f[a],1:s1-2,2:s2-1)
  tmpC=tmpA-view(exFLD.f[a],2:s1-1,1:s2-2)
  if doDIV
    dFLDdx.f[a]=tmpB./Γ.DXC.f[a]
    dFLDdy.f[a]=tmpC./Γ.DYC.f[a]
  else
    dFLDdx.f[a]=tmpB
    dFLDdy.f[a]=tmpC
  end
end

return dFLDdx, dFLDdy
end

function gradient(inFLD::MeshArray,iDXC::MeshArray,iDYC::MeshArray)

exFLD=exchange(inFLD,1)
dFLDdx=similar(inFLD)
dFLDdy=similar(inFLD)

for a=1:inFLD.grid.nFaces
  (s1,s2)=size(exFLD.f[a])
  tmpA=view(exFLD.f[a],2:s1-1,2:s2-1)
  tmpB=tmpA-view(exFLD.f[a],1:s1-2,2:s2-1)
  tmpC=tmpA-view(exFLD.f[a],2:s1-1,1:s2-2)
  dFLDdx.f[a]=tmpB.*iDXC.f[a]
  dFLDdy.f[a]=tmpC.*iDYC.f[a]
end

return dFLDdx, dFLDdy
end

## mask methods

function mask(fld::MeshArray)
fldmsk=mask(fld,NaN)
return fldmsk
end

"""
    mask(fld::MeshArray, val::Number)

Replace non finite values with val. Other methods:
```
mask(fld::MeshArray)
mask(fld::MeshArray, val::Number, noval::Number)
```
"""
function mask(fld::MeshArray, val::Number)
  fldmsk=similar(fld)
  for a=1:fld.grid.nFaces
    tmp1=copy(fld.f[a])
    replace!(x -> !isfinite(x) ? val : x, tmp1 )
    fldmsk.f[a]=tmp1
  end
  return fldmsk
end

function mask(fld::MeshArray, val::Number, noval::Number)
  fldmsk=similar(fld)
  for a=1:fld.grid.nFaces
    tmp1=copy(fld.f[a])
    replace!(x -> x==noval ? val : x, tmp1  )
    fldmsk.f[a]=tmp1
  end
  return fldmsk
end

## convergence methods

"""
    convergence(uFLD::MeshArray,vFLD::MeshArray)

Compute convergence of a vector field
"""
function convergence(uFLD::MeshArray,vFLD::MeshArray)

#important note:
#  Normally uFLD, vFLD should not contain any NaN;
#  if otherwise then something this may be needed:
#  uFLD=mask(uFLD,0.0); vFLD=mask(vFLD,0.0);

CONV=similar(uFLD)

(tmpU,tmpV)=exch_UV(uFLD,vFLD)
for a=1:tmpU.grid.nFaces
  (s1,s2)=size(uFLD.f[a])
  tmpU1=view(tmpU.f[a],1:s1,1:s2)
  tmpU2=view(tmpU.f[a],2:s1+1,1:s2)
  tmpV1=view(tmpV.f[a],1:s1,1:s2)
  tmpV2=view(tmpV.f[a],1:s1,2:s2+1)
  CONV.f[a]=tmpU1-tmpU2+tmpV1-tmpV2
end

return CONV
end

## smooth function

"""
    smooth(FLD::MeshArray,DXCsm::MeshArray,DYCsm::MeshArray,Γ::NamedTuple)

Smooth out scales below DXCsm / DYCsm via diffusion
"""
function smooth(FLD::MeshArray,DXCsm::MeshArray,DYCsm::MeshArray,Γ::NamedTuple)

#important note:
#input FLD should be land masked (NaN/1) by caller if needed

#get land masks (NaN/1):
mskC=fill(1.0,FLD) + 0.0 * mask(FLD)
(mskW,mskS)=gradient(FLD,Γ,false)
mskW=fill(1.0,FLD) + 0.0 * mask(mskW)
mskS=fill(1.0,FLD) + 0.0 * mask(mskS)

#replace NaN with 0. in FLD and land masks:
FLD=mask(FLD,0.0)
mskC=mask(mskC,0.0)
mskW=mask(mskW,0.0)
mskS=mask(mskS,0.0)

#get inverse grid spacing:
iDXC=similar(FLD)
iDYC=similar(FLD)
for a=1:FLD.grid.nFaces
  iDXC.f[a]=1.0./Γ.DXC.f[a]
  iDYC.f[a]=1.0./Γ.DYC.f[a]
end

#Before scaling the diffusive operator ...
tmp0=DXCsm*iDXC*mskW;
tmp00=maximum(tmp0);
tmp0=DYCsm*iDYC*mskS;
tmp00=max(tmp00,maximum(tmp0));

#... determine a suitable time period:
nbt=ceil(1.1*2*tmp00^2);
dt=1.;
T=nbt*dt;
#println("nbt="*"$nbt")

#diffusion operator times DYG / DXG
KuxFac=mskW*DXCsm*DXCsm/T/2.0*Γ.DYG
KvyFac=mskS*DYCsm*DYCsm/T/2.0*Γ.DXG

#time steping factor:
dtFac=dt*mskC/Γ.RAC

#loop:
for it=1:nbt
  (dTdxAtU,dTdyAtV)=gradient(FLD,iDXC,iDYC);
  tmpU=similar(FLD)
  tmpV=similar(FLD)
  for a=1:FLD.grid.nFaces
      tmpU.f[a]=dTdxAtU.f[a].*KuxFac.f[a];
      tmpV.f[a]=dTdyAtV.f[a].*KvyFac.f[a];
  end
  tmpC=convergence(tmpU,tmpV);
  for a=1:FLD.grid.nFaces
      FLD.f[a]=FLD.f[a]-dtFac.f[a].*tmpC.f[a];
  end
end

#Apply land mask (NaN/1) to end result:
mskC=mask(mskC,NaN,0.0)
FLD=mskC*FLD

return FLD

end

## ThroughFlow function

"""
    ThroughFlow(VectorField,IntegralPath,Γ::NamedTuple)

Compute transport through an integration path
"""
function ThroughFlow(VectorField,IntegralPath,Γ::NamedTuple)

    #Note: vertical intergration is not always wanted; left for user to do outside

    U=VectorField["U"]
    V=VectorField["V"]

    nd=ndims(U)
    #println("nd=$nd and d=$d")

    n=fill(1,4)
    tmp=size(U)
    n[1:nd].=tmp[1:nd]

    haskey(VectorField,"factors") ? f=VectorField["factors"] : f=Array{String,1}(undef,0)
    haskey(VectorField,"dimensions") ? d=VectorField["dimensions"] : d=Array{String,1}(undef,nd)

    #a bit of a hack to distinguish gcmfaces v gcmarray indexing:
    isdefined(U,:fIndex) ? ndoffset=1 : ndoffset=0
    length(d)!=nd+ndoffset ? error("inconsistent specification of dims") : nothing

    trsp=Array{Float64}(undef,1,n[3],n[4])
    do_dz=sum(f.=="dz")
    do_dxory=sum(f.=="dxory")

    for i3=1:n[3]
        #method 1: quite slow
        #mskW=IntegralPath.mskW
        #do_dxory==1 ? mskW=mskW*Γ.DYG : nothing
        #do_dz==1 ? mskW=Γ.DRF[i3]*mskW : nothing
        #mskS=IntegralPath.mskS
        #do_dxory==1 ? mskS=mskS*Γ.DXG : nothing
        #do_dz==1 ? mskS=Γ.DRF[i3]*mskS : nothing
        #
        #method 2: less slow
        tabW=IntegralPath.tabW
        tabS=IntegralPath.tabS
        for i4=1:n[4]
            #method 1: quite slow
            #trsp[1,i3,i4]=sum(mskW*U[:,:,i3,i4])+sum(mskS*V[:,:,i3,i4])
            #
            #method 2: less slow
            trsp[1,i3,i4]=0.0
            for k=1:size(tabW,1)
                (a,i1,i2,w)=tabW[k,:]
                do_dxory==1 ? w=w*Γ.DYG.f[a][i1,i2] : nothing
                do_dz==1 ? w=w*Γ.DRF[i3] : nothing
                isdefined(U,:fIndex) ? u=U.f[a,i3,i4][i1,i2] : u=U.f[a][i1,i2,i3,i4]
                trsp[1,i3,i4]=trsp[1,i3,i4]+w*u
            end
            for k=1:size(tabS,1)
                (a,i1,i2,w)=tabS[k,:]
                do_dxory==1 ? w=w*Γ.DXG.f[a][i1,i2] : nothing
                do_dz==1 ? w=w*Γ.DRF[i3] : nothing
                isdefined(V,:fIndex) ? v=V.f[a,i3,i4][i1,i2] : v=V.f[a][i1,i2,i3,i4]
                trsp[1,i3,i4]=trsp[1,i3,i4]+w*v
            end
        end
    end

    nd+ndoffset<4 ? trsp=dropdims(trsp,dims=3) : nothing
    nd+ndoffset<3 ? trsp=dropdims(trsp,dims=2) : nothing
    nd+ndoffset==2 ? trsp=trsp[1] : nothing

    return trsp
end

## LatitudeCircles function

"""
    LatitudeCircles(LatValues,Γ::NamedTuple)

Compute integration paths that follow latitude circles
"""
function LatitudeCircles(LatValues,Γ::NamedTuple)

    LatitudeCircles=Array{NamedTuple}(undef,length(LatValues))

    for j=1:length(LatValues)
        mskCint=1*(Γ.YC .>= LatValues[j])
        mskC=similar(mskCint)
        mskW=similar(mskCint)
        mskS=similar(mskCint)

        mskCint=exchange(mskCint,1)

        for i=1:mskCint.grid.nFaces
            tmp1=mskCint.f[i]
            # tracer mask:
            tmp2=tmp1[2:end-1,1:end-2]+tmp1[2:end-1,3:end]+
            tmp1[1:end-2,2:end-1]+tmp1[3:end,2:end-1]
            mskC.f[i]=1((tmp2.>0).&(tmp1[2:end-1,2:end-1].==0))
            # velocity masks:
            mskW.f[i]=tmp1[2:end-1,2:end-1] - tmp1[1:end-2,2:end-1]
            mskS.f[i]=tmp1[2:end-1,2:end-1] - tmp1[2:end-1,1:end-2]
        end

        function MskToTab(msk::MeshArray)
          n=Int(sum(msk .!= 0)); k=0
          tab=Array{Int,2}(undef,n,4)
          for i=1:msk.grid.nFaces
            a=msk.f[i]
            b=findall( a .!= 0)
            for ii in eachindex(b)
              k += 1
              tab[k,:]=[i,b[ii][1],b[ii][2],a[b[ii]]]
            end
          end
          return tab
        end

        LatitudeCircles[j]=(lat=LatValues[j],
        tabC=MskToTab(mskC),tabW=MskToTab(mskW),tabS=MskToTab(mskS))
    end

    return LatitudeCircles

end

##


text:
ac = ExpectedImprovement()
x_premade = rand(2, 10)*15.0 .- [5.0; 0.0]
y_premade = -1 * [branin(x_premade[:, i], noiselevel=0) for i in 1:size(x_premade, 2)]

model_premade = ElasticGPE(2, mean = MeanConst(-10.0),
                           kernel = SEArd([0., 0.], 5.), logNoise = -2.0)
append!(model_premade, x_premade, y_premade)

# (#7)
@testset "Initial Sampling Tracking" begin
    opt = BOpt(x -> branin(x, noiselevel = 0),
               ElasticGPE(2, mean = MeanConst(-10.),
                          kernel = SEArd([0., 0.], 5.),
                          logNoise = -2., capacity = 3000),
               ac,
               MAPGPOptimizer(every = 50, noisebounds = [-4, 3],
                             kernbounds = [[-1, -1, 0], [4, 4, 10]],
                             maxeval = 40),
               [-5., 0.], [10., 15.],
               maxiterations = 10,
               sense = Min,
               verbosity = Silent,
               initializer_iterations = 10)
    boptimize!(opt)

    @test opt.observed_optimum == Int(opt.sense) * maximum(opt.model.y)
    @test length(opt.model.y) == 10
end

# (#7)
@testset "Pre-made model" begin
    opt = BOpt(x -> branin(x, noiselevel = 0),
               model_premade,
               ac,
               MAPGPOptimizer(every = 50, noisebounds = [-4, 3],
                             kernbounds = [[-1, -1, 0], [4, 4, 10]],
                             maxeval = 40),
               [-5., 0.], [10., 15.],
               maxiterations = 10,
               sense = Min,
               verbosity = Silent,
               initializer_iterations = 5)

    @test opt.observed_optimum   == Int(opt.sense) * maximum(y_premade)
    @test opt.observed_optimizer == x_premade[:, argmax(y_premade)]
end

# (#7)
@testset "Initial iterations to 0 on pre-made model" begin
    ac = ExpectedImprovement()
    opt = BOpt(x -> branin(x, noiselevel = 0),
               model_premade,
               ac,
               MAPGPOptimizer(every = 50, noisebounds = [-4, 3],
                             kernbounds = [[-1, -1, 0], [4, 4, 10]],
                             maxeval = 40),
               [-5., 0.], [10., 15.],
               maxiterations = 0,
               sense = Min,
               verbosity = Silent,
               initializer_iterations = 0)
    boptimize!(opt)

    @test opt.acquisition.τ == maximum(y_premade)
    @test length(opt.model.x) == length(x_premade)
    @test length(opt.model.y) == length(y_premade)

    opt.iterations.N = 5
    boptimize!(opt)
    @test length(opt.model.y) == length(y_premade) + 5
end


text:
ev = CPEG.CPEGWorkspace()

ev.scale.dscale = 124.92
ev.scale.tscale = 0.324
ev.scale.uscale = 1.45

# scaled r,v,a
rs = SA[23.45,434.35,-53.23]
vs = SA[-.043,34.2,-34.22]
as = SA[-2.3,4.344,5.32]

# unscale rs vs
r,v = CPEG.unscale_rv(ev.scale,rs,vs)

# rescale and compare
rs2,vs2 = CPEG.scale_rv(ev.scale,r,v)
@test rs ≈ rs2 atol = 1e-12
@test vs ≈ vs2 atol = 1e-12

# unscale vs as
v2,a2 = CPEG.unscale_va(ev.scale,vs,as)

# compare to the first one
@test v2 ≈ v rtol = 1e-6

v2_2s, a2_2s = CPEG.scale_va(ev.scale,v2,a2)

@test v2_2s ≈ vs atol = 1e-12
@test a2_2s ≈ as atol = 1e-12


text:
##
using MRIgeneralizedBloch
using DifferentialEquations
using StaticArrays
using Test
R2slT = precompute_R2sl()

## choose random parameters
Npulse = 500
α = π/2 * rand() * sin.(8π * rand() * (1:Npulse) / Npulse)
α[1] = π
# TRF = similar(α)
# TRF .= 1e-3 # long TRF to make Graham and gBloch match better
TRF = 100e-6 .+ 400e-6 * rand(Npulse)
TRF[1] = 500e-6

TR = 3.5e-3

B1 = 0.7 + 0.6 * rand()
ω0 = 1000 * rand()
m0s = 0.4 * rand()
R1f = 0.15 + 0.2 * rand()
R1s = 2 + rand()
R2f = 1 / (40e-3 + 100e-3 * rand())
T2s = 5e-6 + 10e-6 * rand()
Rx = 40 * rand()

## gBloch model with IDE vs Graham vs linear approximation: complex signal
s_gBloch_IDE    = calculatesignal_gbloch_ide(  α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, Ncyc=4)
s_Graham        = calculatesignal_graham_ode(  α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, Ncyc=4)
s_gBloch_linear = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT)
@test s_gBloch_IDE ≈ s_Graham        rtol = 5e-2
@test s_gBloch_IDE ≈ s_gBloch_linear rtol = 5e-2

## gBloch model with IDE vs Graham vs linear approximation: magnetixation
s_gBloch_IDE = calculatesignal_gbloch_ide(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s; Ncyc=4, output=:realmagnetization)
s_Graham     = calculatesignal_graham_ode(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s; Ncyc=4, output=:realmagnetization)
s_gBloch_linear = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; output=:realmagnetization)
s_gBloch_linear_m = similar(s_gBloch_IDE)
for i ∈ eachindex(s_gBloch_linear)
    s_gBloch_linear_m[i,:] = s_gBloch_linear[i][[1;2;3;5;6]]
end
@test s_gBloch_IDE ≈ s_Graham          rtol = 5e-2
@test s_gBloch_IDE ≈ s_gBloch_linear_m rtol = 5e-2


text:
using FactCheck

import LLVM
import LLVM: Ast, IntPred

include("util.jl")

test_ast(typ, val, str) = begin
    ast = Ast.Mod("<string>", nothing, nothing, [
            Ast.GlobalDef(
                Ast.GlobalVar(name=Ast.UnName(0), typ=typ, init=val)),
            Ast.GlobalDef(
                Ast.GlobalVar(name=Ast.UnName(1), typ=Uint32, init=nothing)),
            Ast.GlobalDef(
                Ast.GlobalVar(name=Ast.UnName(2), typ=Uint32, init=nothing)),
        ])
    astr ="""; ModuleID = '<string>'

@0 = $str
@1 = external global i32
@2 = external global i32
"""
    return (ast, astr)
end

LLVM.with_ctx() do ctx

    facts("test constants") do

        context("integer") do
            ast, asm = test_ast(Uint32,
                                one(Uint32),
                                "global i32 1")
            check_result(ctx, ast, asm)
        end

        context("wide integer") do
            ast, asm  = test_ast(Ast.IntType(65),
                                 Ast.ConstInt(65,1),
                                 "global i65 1")
            check_result(ctx, ast, asm)
        end

        context("big wide integer") do
            ast, asm = test_ast(Ast.IntType(66),
                                Ast.ConstInt(66, 20000000000000000000),
                                "global i66 20000000000000000000")
            check_result(ctx, ast, asm)
        end

        context("negative wide integer") do
            ast, asm = test_ast(Ast.IntType(65),
                                Ast.ConstInt(65, 36893488147419103231),
                                "global i65 -1")
            check_result(ctx, ast, asm)
        end

        context("half") do
            test_float_roundtrip(Float16, ctx, 1e3)
            ast, asm = test_ast(Float16,
                                Ast.ConstFloat(float16(3.5645)),
                                "global half 0xH4321")
            check_result(ctx, ast, asm)
        end

        context("single") do
            test_float_roundtrip(Float32, ctx, 1e3)
            ast, asm = test_ast(Float32,
                                Ast.ConstFloat(1.0f0),
                                "global float 1.000000e+00")
            check_result(ctx, ast, asm)
        end

        context("double") do
            test_float_roundtrip(Float64, ctx, 1e3)
            ast, asm = test_ast(Float64,
                                Ast.ConstFloat(1.0),
                                "global double 1.000000e+00")
            check_result(ctx, ast, asm)
        end

        # crazy float types are missing...

        context("struct") do
            ast, asm = test_ast(Ast.StructType(false, [Uint32, Uint32]),
                                Ast.ConstStruct(nothing,
                                                false,
                                                [Ast.ConstInt(32, 1),
                                                 Ast.ConstInt(32, 1)]),
                                "global { i32, i32 } { i32 1, i32 1 }")
            check_result(ctx, ast, asm)
        end

        context("dataarray") do
            ast, asm = test_ast(Ast.ArrayType(Int32, 3),
                                Ast.ConstArray(Int32,
                                    [Ast.ConstInt(32, uint32(i)) for i in (1,2,1)]),
                                "global [3 x i32] [i32 1, i32 2, i32 1]")
            check_result(ctx, ast, asm)
        end

        context("array") do
            ast, asm = test_ast(
                Ast.ArrayType(Ast.StructType(false, [Uint32]), 3),
                Ast.ConstArray(Ast.StructType(false, [Uint32]),
                    [Ast.ConstStruct(nothing, false, [Ast.ConstInt(32, i)]) for i in (1,2,3)]),
            "global [3 x { i32 }] [{ i32 } { i32 1 }, { i32 } { i32 2 }, { i32 } { i32 3 }]")
            check_result(ctx, ast, asm)
        end

        context("undef") do
            ast, asm = test_ast(Int32, Ast.ConstUndef(Int32), "global i32 undef")
            check_result(ctx, ast, asm)
        end

        context("binop / cast") do
            ast, asm = test_ast(Int64,
                        Ast.ConstAdd(false, false,
                            Ast.ConstPtrToInt(
                                Ast.ConstGlobalRef(Ptr{Int32}, Ast.UnName(1)),
                                Int64),
                            Ast.ConstInt(64, 2)),
                         "global i64 add (i64 ptrtoint (i32* @1 to i64), i64 2)")
            check_result(ctx, ast, asm)
        end

        context("binop / cast signed/unsigned wrap") do
            ast, asm = test_ast(Int64,
                        Ast.ConstAdd(true, false,
                            Ast.ConstPtrToInt(
                                Ast.ConstGlobalRef(Ptr{Int32}, Ast.UnName(1)),
                                Int64),
                            int64(2)),
                        "global i64 add nsw (i64 ptrtoint (i32* @1 to i64), i64 2)")
            check_result(ctx, ast, asm)
        end

        context("icmp") do
            ast, asm = test_ast(Ast.IntType(1),
                                Ast.ConstICmp(
                                    Ast.SGE(), # todo this sucks
                                    Ast.ConstGlobalRef(Ptr{Int32}, Ast.UnName(1)),
                                    Ast.ConstGlobalRef(Ptr{Int32}, Ast.UnName(2))),
                                "global i1 icmp sge (i32* @1, i32* @2)")
            check_result(ctx, ast, asm)
        end

        context("getelementptr") do
            ast, asm = test_ast(Ptr{Int32},
                                Ast.ConstGetElementPtr(true,
                                    Ast.ConstGlobalRef(Ptr{Int32}, Ast.UnName(1)),
                                    [Ast.ConstInt(64, 27)]),
                            "global i32* getelementptr inbounds (i32* @1, i64 27)")
            check_result(ctx, ast, asm)
        end

        #LLVM.dispose!(ctx)
    end
end


text:
using AnyMOD, Cbc, Test

@testset "run_rest" begin

    # create model
    anyM = anyModel("testModel","testModel", objName = "test", shortExp = 10, checkRng = 1e8)

    createOptModel!(anyM)
    setObjective!(:costs,anyM)

    # solve model
    set_optimizer(anyM.optModel,Cbc.Optimizer)
    optimize!(anyM.optModel)

    # write data reports
    reportResults(:summary,anyM, rtnOpt = (:raw,:rawDf,:csv,:csvDf))
    reportResults(:exchange,anyM, rtnOpt = (:raw,:rawDf,:csv,:csvDf))
    reportResults(:costs,anyM, rtnOpt = (:raw,:rawDf,:csv,:csvDf))
    reportTimeSeries(:electricity, anyM, rtnOpt = (:raw,:rawDf,:csv,:csvDf))
    reportTimeSeries(:electricity, anyM, rtnOpt = (:raw,:rawDf,:csv,:csvDf), mergeVar = false)

    # create plots
    plotEnergyFlow(:sankey,anyM, rmvNode = ("electricity; export","import"))
    plotTree(:region,anyM)
    plotEnergyFlow(:graph,anyM)
    moveNode!(anyM,("coal",[0.1,0.1]))

    @test length(anyM.report) == 33
    @test round(objective_value(anyM.optModel),digits = 1) == 137566.6

    # create additional models with several errors
    anyM = anyModel(["testModel","errorTest"],"testModel", objName = "test", shortExp = 10, checkRng = 1e8)

    err = false
    try
        createOptModel!(anyM)
    catch
        err = true
    end
    @test err

    @test length(anyM.report) == 46

end


text:
export ExtensionFrame,
    extensionframe,
    ⇐, →

using BasisFunctions: default_in_support, unsafe_eval_element,
    default_dict_innerproduct, _default_unsafe_eval_element_in_grid, Measure

import BasisFunctions: superdict, support, similardictionary, isbasis,
    isbiorthogonal, isorthogonal, isorthonormal, hasinterpolationgrid, hastransform,
    hasantiderivative, name, string, strings, dict_in_support, iscompatible, *,
    unsafe_eval_element1, interpolation_grid, measure, restrict, innerproduct,
    innerproduct_native, gramdual, modifiersymbol, GridSampling, in_support,
    points, support

include("submeasure.jl")

export ExtensionFrame
"""
    struct ExtensionFrame{S,T} <: DerivedDict{S,T}

An ExtensionFrame is the restriction of a basis to a subset of its domain. This
results in a frame that implicitly represents extensions of functions on the
smaller set to the larger set.

Created by the function `extensionframe`.

See also [`extensionframe`](@ref)
```jldocs
julia> extensionframe(Fourier(10),0.0..0.5)
```
"""
struct ExtensionFrame{S,T,DICT<:Dictionary{S,T},DOM<:Domain} <: DerivedDict{S,T}
    domain      ::  DOM
    basis       ::  DICT

    function ExtensionFrame{S,T,DICT,DOM}(domain::DOM, basis::DICT) where {S,T,DICT<:Dictionary{S,T},DOM<:Domain}
        # @assert isbasis(basis)
        new(domain, basis)
    end
end

ExtensionFrame(domain::Domain, basis::Dictionary{S,T}) where {S,T} =
    ExtensionFrame{S,T,typeof(basis),typeof(domain)}(domain, basis)
ExtensionFrame(basis::Dictionary, domain::Domain) =
    ExtensionFrame(domain, basis)

superdict(f::ExtensionFrame) = f.basis

export basis
"""
    basis(f::ExtensionFrame)

Return the dictionary on the bouding box.
"""
basis(f::ExtensionFrame) = superdict(f)
support(f::ExtensionFrame) = f.domain
support(f::ExtensionFrame, i) = f.domain

similardictionary(f::ExtensionFrame, dict::Dictionary) = ExtensionFrame(support(f), dict)

isbasis(f::ExtensionFrame) = false
isbiorthogonal(f::ExtensionFrame, measure::Measure) = false
isorthogonal(f::ExtensionFrame, measure::Measure) = false
isorthonormal(f::ExtensionFrame, measure::Measure) = false

# The following properties do not hold for extension frames
# - there is no interpolation grid
hasinterpolationgrid(f::ExtensionFrame) = false
# - there is no unitary transform
hastransform(f::ExtensionFrame) = false
hastransform(f::ExtensionFrame, dgs) = false
# - there is no antiderivative (in general)
hasantiderivative(f::ExtensionFrame) = false

iscompatible(d1::ExtensionFrame, d2::ExtensionFrame) = iscompatible(basis(d1),basis(d2))

function (*)(d1::ExtensionFrame, d2::ExtensionFrame, args...)
    @assert iscompatible(d1,d2)
    (mset, mcoef) = (*)(basis(d1),basis(d2),args...)
    df = ExtensionFrame(support(d1) ∩ support(d2),mset)
    (df, mcoef)
end

unsafe_eval_element1(dict::ExtensionFrame, idx::Int, x) =
    in_support(dict, idx, x) ? unsafe_eval_element(basis(dict), idx, x) : zero(codomaintype(dict))
unsafe_eval_element1(dict::ExtensionFrame, idx, x) =
    in_support(dict, idx, x) ? unsafe_eval_element(basis(dict), idx, x) : zero(codomaintype(dict))
unsafe_eval_element1(dict::ExtensionFrame, idx, grid::AbstractGrid) =
    _default_unsafe_eval_element_in_grid(dict, idx, grid)
in_support(dict::ExtensionFrame, idx::Int, x) =
    default_in_support(dict, idx, x)
in_support(dict::ExtensionFrame, idx, x) =
    default_in_support(dict, idx, x)

function interpolation_grid(f::ExtensionFrame)
    @warn "interpolation_grid called on extensionframe $f"
    subgrid(interpolation_grid(basis(f)),support(f))
end

export extensionframe
"""
    extensionframe(basis::Dictionary, domain::Domain)

    extensionframe(domain::Domain, basis::Dictionary)

Make an ExtensionFrame, but match tensor product domains with tensor product sets
in a suitable way.

# examples
For example: an interval ⊗ a disk (= a cylinder) combined with a 3D Fourier series, leads to a
tensor product of a Fourier series on the interval ⊗ a 2D Fourier series on the disk.
```jldocs
julia> extensionframe(Fourier(10)^2,Disk(.4, SVector(0.5,0.5)))
Dictionary 𝔼(F ⊗ F)

𝔼   :   Extension frame, from A mapped domain based on the 2-dimensional unit ball to 0.0..1.0 (Unit) x 0.0..1.0 (Unit)
F   :   Fourier series
            ↳ length = 10
            ↳ Float64 -> Complex{Float64}
            ↳ support = 0.0..1.0 (Unit)

julia> extensionframe(Fourier(10)^2,(0.0..0.5)^2)
Dictionary 𝔼(F) ⊗ 𝔼(F)

𝔼   :   Extension frame, from 0.0..0.5 to 0.0..1.0 (Unit)
F   :   Fourier series
            ↳ length = 10
            ↳ Float64 -> Complex{Float64}
            ↳ support = 0.0..1.0 (Unit)

```
"""
extensionframe(domain::Domain, basis::Dictionary) = ExtensionFrame(domain, basis)
extensionframe(basis::Dictionary, domain::Domain) = extensionframe(domain, basis)

"Make an extension frame (symbol ⇐ is \\Leftarrow)"
⇐(basis::Dictionary, domain::Domain) = extensionframe(basis, domain)

import BasisFunctions: →
→(domain::Domain, basis::Dictionary) = extensionframe(basis, domain)

function extensionframe(domain::ProductDomain, basis::TensorProductDict)
    ExtensionFrames = Dictionary[]
    dc = 1
    for i = 1:ncomponents(domain)
        el = component(domain, i)
        range = dc:dc+dimension(el)-1
        push!(ExtensionFrames, ExtensionFrame(el, tensorproduct(component(basis, j) for j in range)...))
        dc += dimension(el)
    end
    tensorproduct(ExtensionFrames...)
end

using BasisFunctions: TupleProductDict

const ExtensionFrameTensor = Union{TupleProductDict{D,NTuple{D,DICT}} where D where {DICT<:ExtensionFrame},
    TupleProductDict{2,Tuple{DICT1,DICT2}} where {DICT1<:ExtensionFrame,DICT2<:ExtensionFrame},
    TupleProductDict{3,Tuple{DICT1,DICT2,DICT3}} where {DICT1<:ExtensionFrame,DICT2<:ExtensionFrame,DICT3<:ExtensionFrame},
    TupleProductDict{4,Tuple{DICT1,DICT2,DICT3,DICT4}} where {DICT1<:ExtensionFrame,DICT2<:ExtensionFrame,DICT3<:ExtensionFrame,DICT4<:ExtensionFrame}}
const ExtensionFrameSuper = Union{<:ExtensionFrame,<:ExtensionFrameTensor}
export ExtensionFrameTensor, ExtensionFrameSuper

hasmeasure(::ExtensionFrame) = true
measure(f::ExtensionFrame) = submeasure(measure(basis(f)), support(f))

innerproduct_native(f1::ExtensionFrame, i, f2::ExtensionFrame, j, measure; options...) =
    innerproduct(superdict(f1), i, superdict(f2), j, measure; options...)

# This routine will be called if we have two mapped dictionaries, where the measure is a subset of a
# mapped measure. We check for this case and undo the mapping.
innerproduct_native(dict1::MappedDict, i, dict2::MappedDict, j, measure; options...) =
    _innerproduct_native(support(measure), supermeasure(measure), dict1, i, dict2, j, measure; options...)

function _innerproduct_native(domain::AbstractInterval, superμ, dict1::MappedDict, i, dict2::MappedDict, j, measure; options...)
    if BasisFunctions.ismappedmeasure(superμ) && iscompatible(dict1, dict2) && iscompatible(forward_map(dict1), forward_map(superμ))
        supermap = forward_map(superμ)
        newdomain = Interval(applymap(inv(supermap), infimum(domain)), applymap(inv(supermap), supremum(domain)))
        innerproduct_native(superdict(dict1), i, superdict(dict2), j, submeasure(supermeasure(superμ), newdomain))
    else
        dict_default_innerproduct(dict1, i, dict2, j, measure; options...)
    end
end

function _innerproduct_native(domain::UnionDomain, superμ, dict1::MappedDict, i, dict2::MappedDict, j, measure; options...)
    z = zero(coefficienttype(dict1))
    for d in components(domain)
        z += _innerproduct_native(d, superμ, dict1, i, dict2, j, measure; options...)
    end
    z
end

_innerproduct_native(domain, superμ, dict1::MappedDict, i, dict2::MappedDict, j, measure; options...) =
    BasisFunctions.default_dict_innerproduct(dict1, i, dict2, j, measure; options...)

export extensiondual
"""
    extensiondual(dict::ExtensionFrame, measure::Measure; options...)

Return the extension frame of the dual of the dictionary on the boundingbox.
"""
extensiondual(dict::ExtensionFrame, measure::Measure; options...) =
    extensionframe(support(dict), gramdual(superdict(dict), supermeasure(measure); options...),)
extensiondual(dict::ExtensionFrameTensor, measure::Measure; options...) =
    TensorProductDict(map((dicti,measurei)->extensionframe(support(dicti), gramdual(superdict(dicti), supermeasure(measurei); options...)),
        components(dict), components(measure))...)

function gramdual(dict::ExtensionFrame, measure::Measure; options...)
    @debug "Are you sure you want `dualtype=gramdual` and not `extensiondual`"
    default_gramdual(dict, measure; options...)
end

superdict(dict::ExtensionFrameTensor) = TensorProductDict(map(superdict, components(dict))...)
basis(dict::ExtensionFrameTensor) = TensorProductDict(map(basis, components(dict))...)
support(dict::ExtensionFrameTensor) = ProductDomain((map(support, components(dict)))...)

## Printing

Display.object_parentheses(d::ExtensionFrame) = true
Display.stencil_parentheses(d::ExtensionFrame) = true
Display.displaystencil(d::ExtensionFrame) = _stencil(d, superdict(d), support(d))
_stencil(d::ExtensionFrame, dict, domain) = [domain, " → ", dict]

GridSampling(dgs::GridBasis, grid::AbstractGrid, domain::Domain, scaling) =
    GridSampling(GridBasis{coefficienttype(dgs)}(subgrid(grid, domain)), scaling=scaling)


text:
using HDF5
using JLD
using DataFrames

fname = joinpath(tempdir(), "mydata-$(randstring()).jld")

df =  DataFrames.DataFrame(Any[[2:6;], pi*[1:5;]])
df2 = DataFrames.DataFrame(a = [1:5;], b = pi * [1:5;])

file = jldopen(fname, "w")
write(file, "df", df)
write(file, "df2", df2)
close(file)

file = jldopen(fname, "r")
x = read(file, "df")
y = read(file, "df2")
close(file)

using Compat.Test
@test isequal(df, x)
@test isequal(df2, y)

# Testing issue #236
fname = joinpath(tempdir(), "int_str_data-$(randstring()).jld")
df3 = DataFrames.DataFrame(A = [1:4;], B = ["M", "F", "F", "M"])
file = jldopen(fname, "w")
write(file, "df3", df3)
close(file)
file = jldopen(fname, "r")
iob = IOBuffer()
dump(iob, file)
x = read(file, "df3")
@test isequal(df3, x)


text:
ENV["JULIA_CUDA_SILENT"] = true
# Change the working directory to the file location
cd(@__DIR__)

using DrMZ
using Parameters: @with_kw
using Flux.NNlib
using LaTeXStrings
using ColorSchemes
using Flux
using Plots; pyplot()
using Plots.PlotMeasures
using Random
using Printf

default(fontfamily="serif",frame=:box,grid=:hide,palette=:viridis,markeralpha=0.4,dpi=200,legendfontsize=6);
PyPlot.matplotlib.rc("mathtext",fontset="cm")

@with_kw mutable struct Args
    num_sensors::Int = 128;
    num_train_functions::Int = 500;
    num_test_functions::Int = Int(2*num_train_functions);
    num_sol_points::Int = 100;
    L1::Float64 = 0.0;
    L2::Float64 = 2*pi;
    tspan::Tuple = (0.0,1.0);
    n_epoch::Int = 25000;
    N::Int = num_sensors;
    branch_layers::Int = 2; # Branch depth
    branch_neurons::Array{Tuple} = [(Int(num_sensors),num_sensors),(num_sensors,N)]; # Branch layers input and output dimensions - quantity of tuples must match branch depth
    branch_activations::Array = [tanh,identity]; # Activation functions for each branch layer
    trunk_layers::Int = 3; # Trunk depth
    trunk_neurons::Array{Tuple} = [(2,num_sensors),(num_sensors,num_sensors),(num_sensors,N)]; # Trunk layers input and output dimensions - quantity of tuples must match trunk depth
    trunk_activations::Array = [tanh,tanh,tanh]; # Activation functions for each trunk layer
end

function generate_train(pde_function,pde_function_handle;kws...)

    args = Args(;);

    # Generate training and testing data
    train_data, test_data = generate_periodic_train_test(args.tspan,args.num_sensors,args.num_train_functions,args.num_test_functions,args.num_sol_points,pde_function_handle);
    save_data(train_data,test_data,args.num_train_functions,args.num_test_functions,args.num_sol_points,pde_function)

    branch = build_dense_model(args.branch_layers,args.branch_neurons,args.branch_activations);
    trunk = build_dense_model(args.trunk_layers,args.trunk_neurons,args.trunk_activations)

    # Train the operator neural network
    branch, trunk = train_model(branch,trunk,args.n_epoch,train_data,test_data,pde_function);
    save_model(branch,trunk,args.n_epoch,pde_function)

end

function generate_opnn_results(pde_function,pde_function_handle;random_integer = "none",kws...)

    args = Args(;);

    # Load all the data
    branch, trunk, train_ic, train_loc, train_sol, test_ic, test_loc, test_sol = load_data(args.n_epoch,args.num_train_functions,args.num_test_functions,pde_function)

    train_MSE = loss_all(branch,trunk,train_ic,train_loc,train_sol);
    test_MSE = loss_all(branch,trunk,test_ic,test_loc,test_sol);
    println("Training set MSE: $train_MSE, test set MSE: $test_MSE")

    dL = abs(args.L2 - args.L1);
    j = reduce(vcat,[0:1:args.num_sensors-1]);
    x = (dL.*j)./args.num_sensors;
    dt = 1e-3;
    t = reduce(vcat,[0:dt:args.tspan[2]]);

    if random_integer == "none"
        rand_int = rand(1:args.num_test_functions);
        ic = test_ic[:,rand_int*args.num_sol_points];
    elseif random_integer == "exact"
        ic = sin.(x/2).^2;
        rand_int = 0;
    else
        rand_int = random_integer;
        ic = test_ic[:,rand_int*args.num_sol_points];
    end

    # Comparisons vs random test IC
    pltrand_test = plot(x,ic,label = "IC",linewidth=:2,xlims=(x[1],x[end]),xlabel=L"x",ylabel=L"u(x)")
    plot!(pltrand_test,x,ic,label = "Sensor locations",seriestype=:scatter)
    display(pltrand_test)
    savefig(pltrand_test, @sprintf("random_ic_%s_%i.png",pde_function,rand_int))

    # M size Fourier solution for comparisons
    u_fourier_M = generate_fourier_solution(args.L1,args.L2,args.tspan,args.num_sensors,ic,pde_function_handle)[1];

    pltexact = plot(heatmap(periodic_fill_domain(x),t,periodic_fill_solution(u_fourier_M),aspect_ratio=dL/args.tspan[2]),fillcolor=cgrad(ColorSchemes.viridis.colors),label=false,xlabel = L"x", ylabel = L"t",xlims=(x[1],x[end]),ylims=(0,args.tspan[2]))
    display(pltexact)
    savefig(pltexact, @sprintf("exact_fourier_solution_random_ic_%s_%i.png",pde_function,rand_int))

    # Operator neural network prediction
    u_predict = predict(branch,trunk,ic,x,t);

    pltpredict = plot(heatmap(periodic_fill_domain(x),t,periodic_fill_solution(u_predict),aspect_ratio=dL/args.tspan[2]),fillcolor=cgrad(ColorSchemes.viridis.colors),label=false,xlabel = L"x", ylabel = L"t",xlims=(x[1],x[end]),ylims=(0,args.tspan[2]))
    display(pltpredict)
    savefig(pltpredict, @sprintf("opnn_solution_random_ic_epochs_%i_%s_%i.png",args.n_epoch,pde_function,rand_int))

    opnn_error = norm_rel_error(u_fourier_M,u_predict);
    opnn_error = plot(t,opnn_error,legend=:outerbottom,foreground_color_legend = nothing,label="Operator network average error $(@sprintf("%.4e",average_error(t,opnn_error)))",xlabel = L"t", ylabel = "Rel. Error",xlims=(t[1],t[end]))
    display(opnn_error)
    savefig(opnn_error, @sprintf("opnn_vs_exact_solution_random_ic_epochs_%i_%s_%i.png",args.n_epoch,pde_function,rand_int))

end

# Generate data and train the model
generate_train("advection_equation",advection_pde!)

# Generate results - enter "none", an integer value (e.g. 975), or "exact" to specify an initial condition
# If testing ability of network to extrapolate beyond the training interval, comment out the generation functions above and adjust the Arg tspan
generate_opnn_results("advection_equation",advection_pde!;random_integer = "exact")


text:
# using Pkg; haskey(Pkg.installed(),"ZipFile") || Pkg.add("ZipFile")
using ZipFile

"Where to download dataset from"
const TREEBANK_URL = "https://nlp.stanford.edu/sentiment/trainDevTestTrees_PTB.zip"

"Where to download dataset to"
const TREEBANK_DIR = joinpath(@__DIR__,"treebank")

const TREEBANK_ZIPNAME = "trainDevTestTrees_PTB.zip"
const TREEBANK_SPLITS = ("train", "dev", "test")
const UNK = "_UNK_"
"""

This utility loads [Stanford Sentiment Treebank](https://nlp.stanford.edu/sentiment/index.html)
sentiment classification dataset. There are 8544, 1101 and 2210 trees in train,
dev and test splits respectively. Each tree represented as a s-expression in raw
data. There are 5 different scales to represent sentiment negativity/positivity.

```
# Usage:
include(Pkg.dir("Knet/data/treebank.jl"))
trn = load_treebank_data(splits=["train"]) # possible splits: train,dev,test
trn = load_treebank_data("train")

# each instance is a SentimentTree
# typeof(trn[1]) == SentimentTree

# build word/tag vocabularies
l2i, w2i, i2l, i2w = build_treebank_vocabs(trn)

# make_data! overwrites the data - replaces words/tags with assigned indices
make_data!(trn, w2i, l2i)

# see the data loader source code and the example for further documentation
```

"""
mutable struct SentimentTree
    label
    children
    data

    SentimentTree(x) = new(x,nothing,nothing)
    SentimentTree(x,y) = new(x,y,nothing)
end

const SentimentTrees = Array{SentimentTree,1}

function isleaf(t::SentimentTree)
    return t.children == nothing
end

function pretty(t::SentimentTree)
    t.children == nothing && return t.label
    join([t.label; map(pretty, t.children)], " ")
end

let
    items = []

    global leaves
    function leaves(t::SentimentTree)
        empty!(items)
        helper(t)
        return items
    end

    function helper(subtree)
        if isleaf(subtree)
            push!(items, subtree)
        else
            for child in subtree.children
                helper(child)
            end
        end
    end
end

let
    nodes = []
    global nonterms
    function nonterms(t::SentimentTree)
        empty!(nodes)
        helper(t)
        return nodes
    end

    function helper(subtree)
        if !isleaf(subtree)
            push!(nodes, subtree)
            map(helper, subtree.children)
        end
    end
end

# data load function(s)
let
    global load_treebank_data
    function load_treebank_data(datadir=TREEBANK_DIR, splits=["train","dev"])
        datadir = abspath(datadir)
        if !isdatadir(datadir)
            mkpath(datadir)
        end

        if !isdata(datadir)
            !iszip(datadir) && download_zip(datadir)
            extract_files(datadir)
        end

        data = map(s->load_treebank_data(datadir, s), splits)
    end

    function load_treebank_data(datadir, split::T) where T <: String
        !in(split, TREEBANK_SPLITS) && error("no such split as $split")
        filename = split*".txt"
        filepath = joinpath(datadir, filename)
        data = read_file(filepath)
    end

    function extract_files(datadir)
        r = ZipFile.Reader(joinpath(datadir, TREEBANK_ZIPNAME))
        for f in r.files
            _, this_file = splitdir(f.name)
            split, _ = splitext(this_file)
            if split in TREEBANK_SPLITS
                lines = readlines(f)
                text = join(lines, "\n")
                file = joinpath(datadir, split*".txt")
                open(file, "w") do f
                    write(f, text)
                end
            end
        end
        close(r)
    end

    function isdata(datadir)
        for split in TREEBANK_SPLITS
            splitfile = joinpath(datadir, split*".txt")
            !isfile(splitfile) && return false
        end
        return true
    end

    function iszip(datadir)
        fullpath = joinpath(datadir, TREEBANK_ZIPNAME)
        return isfile(fullpath)
    end

    function download_zip(datadir)
        dest = joinpath(datadir, TREEBANK_ZIPNAME)
        download(TREEBANK_URL, dest)
    end

    function isdatadir(datadir)
        return isdir(datadir)
    end

    function read_file(file)
        data = open(file, "r") do f
            map(parse_line, readlines(f))
        end
    end

    function parse_line(line)
        ln = replace(line, "\n"=>"")
        tokens = tokenize_sexpr(ln)
        popfirst!(tokens)
        return within_bracket(tokens)[1]
    end

    function tokenize_sexpr(sexpr)
        tokker = r" +|[()]|[^ ()]+"
        filter(t -> t != " ", collect(m.match for m = eachmatch(tokker, sexpr)))
    end

    function within_bracket(tokens, state=1)
        (label, state) = iterate(tokens, state)
        children = []
        while iterate(tokens, state) != nothing
            (token, state) = iterate(tokens, state)
            if token == "("
                (child, state) = within_bracket(tokens, state)
                push!(children, child)
            elseif token == ")"
                return SentimentTree(label, children), state
            else
                push!(children, SentimentTree(token))
            end
        end
    end
end

let
    global build_treebank_vocabs
    function build_treebank_vocabs(trees::SentimentTrees)
        words = Set()
        labels = Set()
        for tree in trees
            push!(words, map(t->t.label, leaves(tree))...)
            push!(labels, map(t->t.label, nonterms(tree))...)
        end
        push!(words, UNK)
        w2i, i2w = build_vocab(words)
        l2i, i2l = build_vocab(labels)
        return l2i, w2i, i2l, i2w
    end

    function build_vocab(xs)
        x2i = Dict(); i2x = Dict()
        for (i,x) in enumerate(xs)
            x2i[x] = i
            i2x[i] = x
        end
        return x2i, i2x
    end
end

function make_data!(trees::SentimentTrees, w2i, l2i)
    for tree in trees
        for leaf in leaves(tree)
            ind = get(w2i, leaf.label, w2i[UNK])
            leaf.data = ind
        end
        for nonterm in nonterms(tree)
            nonterm.data = l2i[nonterm.label]
        end
    end
end


text:
module TestPreamble

include("preamble.jl")

@testset "==ₜ" begin
    @test 1 ==ₜ 1
    @test !(1.0 ==ₜ 1)
end

end  # module


text:
module TerminalColor

export color,black,red,green,green,yellow,blue,magenta,cyan,white

COLORS = Dict(
"black"         => 0,
"red"           => 1,
"green"         => 2,
"yellow"        => 3,
"blue"          => 4,
"magenta"       => 5,
"cyan"          => 6,
"white"         => 7,
"default"       => 9
)

MODES = Dict(
"default"       => 0,
"bold"          => 1,
"underline"     => 4,
"blink"         => 5,
"swap"          => 7,
"hide"          => 8
)

function color(color::String, str::String; background::String="default", mode::String="default")
    FColor = "\e[$(30 + COLORS[color]);"
    BColor =    "$(40 + COLORS[background])m"
    Action = "\e[$(MODES[mode]);m"
    return FColor * BColor * str * Action
end

black(str)   = color("black",   string(str))
red(str)     = color("red",     string(str))
green(str)   = color("green",   string(str))
yellow(str)  = color("yellow",  string(str))
blue(str)    = color("blue",    string(str))
magenta(str) = color("magenta", string(str))
cyan(str)    = color("cyan",    string(str))
white(str)   = color("white",   string(str))

end

max_stars_repo_name: sonosole/TerminalColor

max_issues_repo_name: sonosole/TerminalColor

max_forks_repo_name: sonosole/TerminalColor