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

using Markdown
using InteractiveUtils

# ╔═╡ 03cc6704-70eb-11eb-2f5f-d727783a1f53
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)
end

# ╔═╡ 6d52e124-70eb-11eb-2c56-5dcea746627e
begin
	using MLJ
	using CSV
	using DataFrames
	using PlutoUI
end

# ╔═╡ b738e190-70ef-11eb-04da-8105c8390b73
using Random

# ╔═╡ a9f436f6-70fc-11eb-394e-f7a9755bb779
using Test

# ╔═╡ 28cfd44c-70fc-11eb-0155-f984980e96f6
using Printf # Pkg.add("Printf")

# ╔═╡ 62c32226-70ea-11eb-3914-730e6d2ec0a4
md"""
## C02 - w01: Simple Linear Regression

In this notebook we will use data on house sales in King County to predict house prices using simple (one input) linear regression. we will
  - compute important summary statistics
  - Write a function to compute the Simple Linear Regression weights using the closed form solution
  - Write a function to make predictions of the output given the input feature
  - Turn the regression around to predict the input given the output
  - Compare two different models for predicting house prices
  - *Finally use a predefined MLJ solution*
"""

# ╔═╡ ec3daae4-70eb-11eb-08df-416f041130f6

# ╔═╡ fa6be52e-70eb-11eb-2050-0bd51777d5d0
sales = CSV.File("../../ML_UW_Spec/C02/w01/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ 39ca1334-70ed-11eb-2ee9-5bda94ba96f5
first(sales, 7)

# ╔═╡ 53b101ac-70f1-11eb-3844-3d5f21d9b7aa
size(sales)

# ╔═╡ d2a6fa36-70f7-11eb-3817-b36402af12f2
eltype.(eachcol(sales))

# ╔═╡ 299273f6-70f4-11eb-11c9-035af4af7032
md"""
### Split data into training and testing set
"""

# ╔═╡ 55c029d2-70f4-11eb-00ac-0b7698d7b6e2
function train_test_split(df; split=0.8, seed=42)
	Random.seed!(seed)
	(nr, nc) = size(df)
	nrp = round(Int, nr * split)
	row_ixes = shuffle(1:nr)
	df_train = view(df[row_ixes, :], 1:nrp, 1:nc)
	df_test = view(df[row_ixes, :], nrp+1:nr, 1:nc)
	(df_train, df_test)
end

# ╔═╡ d403506c-70f4-11eb-386f-ad0152978098
sales_train, sales_test = train_test_split(sales);

# ╔═╡ fff80614-70f3-11eb-137d-cbf3fcbaf98f
first(sales_train, 7)

# ╔═╡ 24949900-70f0-11eb-24f4-bd5c33121397
first(sales_test, 3)

# ╔═╡ 11aafb94-70f4-11eb-2438-8320f2e6cb5b
size(sales_train), size(sales_test)

# ╔═╡ 1fb75280-70f8-11eb-1cc3-29a08a026bcc
# mean and sum of a column
mean(sales_train[!, :bedrooms]), sum(sales_train[!, :bedrooms])

# ╔═╡ 799514b8-70f8-11eb-0d03-e5fe5b1c0cb4
sales_train[!, :bedrooms] .* sales_train[!, :bedrooms]

# ╔═╡ 7c9c2c1a-70f9-11eb-3ddc-d90788056518
length(sales_train[!, :bedrooms])

# ╔═╡ cc20e6e0-70f9-11eb-2b0d-7bd7f07d4733
typeof(sales_train[!, :bedrooms])

# ╔═╡ 416a6240-70f9-11eb-1962-e77cb5215607
md"""
### Build a generic simple linear regression function
"""

# ╔═╡ 5b5efd00-70f9-11eb-36f1-75d629b148e6
function slope(X, y)::Float64
	sx, sy = sum(X), sum(y)
	n = length(X)
	num = sum(X .* y) - (1/n) * sx * sy
	den = sum(X .* X) - (1/n) * sx * sx
	num / den
end

# ╔═╡ c40f3ea0-70f9-11eb-2ff0-9b1e9c4f3084
intercept(X, y, slope) = mean(y) - slope * mean(X)

# ╔═╡ 04b3c642-70fa-11eb-3e0b-b9b86108e1af
function simple_linear_regression(X, y)::Tuple{Real, Real}
	_slope = slope(X, y)
	_intercept = intercept(X, y, _slope)
	(_intercept, _slope)
end

# ╔═╡ 4e0f51e4-70fa-11eb-2511-b9813e9c82bb
begin
	test_feature = collect(1:5)
	test_output = 1 .+ test_feature

	(test_intercept, test_slope) = simple_linear_regression(test_feature, test_output)

	@test test_intercept ≈ 1.0
	@test test_slope ≈ 1.0
end

# ╔═╡ 64258a92-70fb-11eb-08d4-b11d5f75639f
begin
	sqft_intercept, sqft_slope = simple_linear_regression(
		sales_train[!, :sqft_living], sales_train[!, :price])

	with_terminal() do
		@printf("Intercept: %.4f\n", sqft_intercept)
		@printf("Slope:  %.4f\n", sqft_slope)
	end
end

# ╔═╡ beee3bd2-70fc-11eb-38cd-e7e8fda2f883
md"""
### Predicting Values

Now that we have the model parameters: intercept and slope we can make predictions.

Create a function to return the predicted output given the input_feature, slope and intercept.

"""

# ╔═╡ deda4470-70fc-11eb-20a2-15820ab34231
function get_regression_predictions(input_feature, intercept::Real,
		slope::Real)
    slope .* input_feature .+ intercept
end

# ╔═╡ 1ed980ae-70fd-11eb-07ed-09b7648af658
md"""
Now that we can calculate a prediction given the slope and intercept let's make a prediction.

**Quiz Question**: Using your Slope and Intercept from above, what is the predicted price for a house with 2650 sqft?
"""

# ╔═╡ 490bda66-70fd-11eb-3cd3-e7ccd064c7ec
with_terminal() do
	house_sqft = 2650
	estimated_price = get_regression_predictions(house_sqft,
		sqft_intercept, sqft_slope)

	@printf("The estimated price for a house with %d squarefeet is %.2f\n",
		house_sqft,
		estimated_price)
end

# ╔═╡ 96f5469a-70fd-11eb-1f17-675eeff21970
md"""
### Residual Sum of Squares

Now that we have a model and can make predictions let's evaluate our model using Residual Sum of Squares (RSS). Recall that RSS is the sum of the squares of the residuals and the residuals is just a fancy word for the difference between the predicted output and the true output.

Create a function (`get_residual_sum_of_squares`) to compute the RSS of a simple linear regression model given the input_feature, output, intercept and slope.

"""

# ╔═╡ b5e12f56-70fd-11eb-1bf7-cbdfa36fc673
function get_residual_sum_of_squares(input_feature, output,
		intercept::Real, slope::Real)
    ## 1. get the predictions
    yhat = get_regression_predictions(input_feature, intercept, slope)
    diff = (yhat .- output)   ## 2. compute the residuals
    rss = sum(diff .* diff)  ## 3. square the residuals and add them up
    return rss
end

# ╔═╡ 25479c96-70ff-11eb-2b86-29c052eb2f87
test_slope

# ╔═╡ 064ec37c-70fe-11eb-3bce-a79f7208c876
md"""
Let's test our `get_residual_sum_of_squares` function by applying it to the test model where the data lie exactly on a line. Since they lie exactly on a line the residual sum of squares should be zero!
"""

# ╔═╡ 1b0643d0-70fe-11eb-250c-53e54279e747
@test get_residual_sum_of_squares(test_feature, test_output,
	test_intercept, test_slope) ≈ 0.0

# ╔═╡ 5e2a432e-70ff-11eb-3e61-bffdc0ebd46e
md"""
Now use your function to calculate the RSS on training data from the squarefeet model calculated above.

**Quiz Question**: According to this function and the slope and intercept from the squarefeet model What is the RSS for the simple linear regression using squarefeet to predict prices on *training* data?

"""

# ╔═╡ 7e8bcf82-70ff-11eb-1c79-a56db6635447
begin
	rss_prices_on_sqft = get_residual_sum_of_squares(sales_train[!, :sqft_living],
		sales_train[!, :price],
        sqft_intercept,
		sqft_slope)
	with_terminal() do
		@printf("The RSS of predicting Prices based on Square Feet is: %1.3e\n", rss_prices_on_sqft)
	end
end

# ╔═╡ 0afc9a02-7100-11eb-0578-8b470fbb1d32
md"""
### Predict the squarefeet given price

What if we want to predict the squarefoot given the price?

Since we have an equation $y = a + b \times x$ we can solve the function for x, so that if we have the intercept (a) and the slope (b) and the price (y) we can solve for the estimated squarefeet (x) with:
$x = \frac{y - a}{b}$.

Create a function (`inverse_regression_predictions`) to compute the inverse regression estimate, i.e. predict the input_feature given the output.
"""

# ╔═╡ 5936bff2-7100-11eb-3673-03eb5174f2b0
function inverse_regression_predictions(output, intercept::Real, slope::Real)
    (output .- intercept) ./ slope
end

# ╔═╡ 7f548db8-7100-11eb-1a71-cfa1b2817f83
md"""
Now that we have a function to compute the squarefeet given the price from our simple regression model let's see how big we might expect a house that costs \$800,000 to be.

**Quiz Question**: According to this function and the regression slope and intercept from (3) what is the estimated square-feet for a house costing \$800,000?

"""

# ╔═╡ 967a4d8e-7100-11eb-369c-a7eef43335a2
begin
	house_price = 800000
	estimated_squarefeet = inverse_regression_predictions(house_price,
		sqft_intercept,
		sqft_slope)

	with_terminal() do
		@printf("The estimated squarefeet for a house worth %.2f is %.3f\n", house_price, estimated_squarefeet)
	end
end

# ╔═╡ f8cb4844-7100-11eb-2036-f1691f9277ba
md"""
### New Model: estimate prices from bedrooms

We have made one model for predicting house prices using squarefeet, but there are many other features in the sales DataFrame.

Use your simple linear regression function to estimate the regression parameters from predicting Prices based on number of bedrooms. Use the *training* data!

"""

# ╔═╡ 1a89af3e-7101-11eb-274d-43bbcf50db13
bedr_intercept, bedr_slope = simple_linear_regression(
	sales_train[!, :bedrooms],
	sales_train[!, :price]
)

# ╔═╡ 4d0ca242-7101-11eb-1001-7bb71d3b8de6
md"""
### Test your Linear Regression Algorithm

Now we have two models for predicting the price of a house. How do we know which one is better? Calculate the RSS on the *test* data (remember this data wasn't involved in learning the model). Compute the RSS from predicting prices using bedrooms and from predicting prices using squarefeet.

**Quiz Question**: Which model (square feet or bedrooms) has lowest RSS on TEST data? Think about why this might be the case.

"""

# ╔═╡ 69a38e9e-7101-11eb-2e57-83829ceb7ab9
begin
	## Compute RSS when using bedrooms on test data:
	rss_prices_on_bedr = get_residual_sum_of_squares(sales_test[!, :bedrooms],
        sales_test[!, :price],
        bedr_intercept,
		bedr_slope)
	with_terminal() do
		@printf("The RSS of predicting Prices based on number of Bedrooms is: %1.3e\n ", rss_prices_on_bedr)
	end
end

# ╔═╡ d043648a-7101-11eb-1b30-8d7ff3b37b7a
begin
	## Compute RSS when using squarefeet on test data:
	rss_prices_on_sqft_t = get_residual_sum_of_squares(sales_test[!, :sqft_living],
        sales_test[!, :price],
        sqft_intercept,
		sqft_slope)
	with_terminal() do
		@printf("The RSS of predicting Prices based on Squarefeet is: %1.3e\n ", rss_prices_on_sqft_t)
	end
end

# ╔═╡ Cell order:
# ╟─62c32226-70ea-11eb-3914-730e6d2ec0a4
# ╠═03cc6704-70eb-11eb-2f5f-d727783a1f53
# ╠═6d52e124-70eb-11eb-2c56-5dcea746627e
# ╟─ec3daae4-70eb-11eb-08df-416f041130f6
# ╠═fa6be52e-70eb-11eb-2050-0bd51777d5d0
# ╠═39ca1334-70ed-11eb-2ee9-5bda94ba96f5
# ╠═53b101ac-70f1-11eb-3844-3d5f21d9b7aa
# ╠═d2a6fa36-70f7-11eb-3817-b36402af12f2
# ╟─299273f6-70f4-11eb-11c9-035af4af7032
# ╠═b738e190-70ef-11eb-04da-8105c8390b73
# ╠═55c029d2-70f4-11eb-00ac-0b7698d7b6e2
# ╠═d403506c-70f4-11eb-386f-ad0152978098
# ╠═fff80614-70f3-11eb-137d-cbf3fcbaf98f
# ╠═24949900-70f0-11eb-24f4-bd5c33121397
# ╠═11aafb94-70f4-11eb-2438-8320f2e6cb5b
# ╠═1fb75280-70f8-11eb-1cc3-29a08a026bcc
# ╠═799514b8-70f8-11eb-0d03-e5fe5b1c0cb4
# ╠═7c9c2c1a-70f9-11eb-3ddc-d90788056518
# ╠═cc20e6e0-70f9-11eb-2b0d-7bd7f07d4733
# ╟─416a6240-70f9-11eb-1962-e77cb5215607
# ╠═5b5efd00-70f9-11eb-36f1-75d629b148e6
# ╠═c40f3ea0-70f9-11eb-2ff0-9b1e9c4f3084
# ╠═04b3c642-70fa-11eb-3e0b-b9b86108e1af
# ╠═a9f436f6-70fc-11eb-394e-f7a9755bb779
# ╠═4e0f51e4-70fa-11eb-2511-b9813e9c82bb
# ╠═28cfd44c-70fc-11eb-0155-f984980e96f6
# ╠═64258a92-70fb-11eb-08d4-b11d5f75639f
# ╟─beee3bd2-70fc-11eb-38cd-e7e8fda2f883
# ╠═deda4470-70fc-11eb-20a2-15820ab34231
# ╟─1ed980ae-70fd-11eb-07ed-09b7648af658
# ╠═490bda66-70fd-11eb-3cd3-e7ccd064c7ec
# ╟─96f5469a-70fd-11eb-1f17-675eeff21970
# ╠═b5e12f56-70fd-11eb-1bf7-cbdfa36fc673
# ╠═25479c96-70ff-11eb-2b86-29c052eb2f87
# ╟─064ec37c-70fe-11eb-3bce-a79f7208c876
# ╠═1b0643d0-70fe-11eb-250c-53e54279e747
# ╟─5e2a432e-70ff-11eb-3e61-bffdc0ebd46e
# ╠═7e8bcf82-70ff-11eb-1c79-a56db6635447
# ╟─0afc9a02-7100-11eb-0578-8b470fbb1d32
# ╠═5936bff2-7100-11eb-3673-03eb5174f2b0
# ╟─7f548db8-7100-11eb-1a71-cfa1b2817f83
# ╠═967a4d8e-7100-11eb-369c-a7eef43335a2
# ╟─f8cb4844-7100-11eb-2036-f1691f9277ba
# ╠═1a89af3e-7101-11eb-274d-43bbcf50db13
# ╟─4d0ca242-7101-11eb-1001-7bb71d3b8de6
# ╠═69a38e9e-7101-11eb-2e57-83829ceb7ab9
# ╠═d043648a-7101-11eb-1b30-8d7ff3b37b7a


text:
const W1SOL = ("(S2.*(S3P.*ws1+(-1).*S1P.*ws3).^2+S2P.^2.*(S3.*ws1.^2+ws3.*((-2).* ...
  S13.*ws1+S1.*ws3))+(-2).*S2P.*(S1P.*S3.*ws1.*ws2+(-1).*S13.*S3P.* ...
  ws1.*ws2+(-1).*S12.*S3P.*ws1.*ws3+(-1).*S13.*S1P.*ws2.*ws3+S1.* ...
  S3P.*ws2.*ws3+S12.*S1P.*ws3.^2+S23.*ws1.*(S3P.*ws1+(-1).*S1P.*ws3) ...
  )+ws2.*(S3P.^2.*((-2).*S12.*ws1+S1.*ws2)+2.*S1P.*S3P.*(S23.*ws1+( ...
  -1).*S13.*ws2+S12.*ws3)+S1P.^2.*(S3.*ws2+(-2).*S23.*ws3))).^(-1).* ...
  (S2P.^2.*S3.*ws1+(-2).*S23.*S2P.*S3P.*ws1+S2.*S3P.^2.*ws1+(-1).* ...
  S1P.*S2P.*S3.*ws2+S1P.*S23.*S3P.*ws2+(-1).*S2P.*S3.*SGP.*ws1.*ws2+ ...
  S23.*S3P.*SGP.*ws1.*ws2+S1P.*S3.*SGP.*ws2.^2+S1P.*S23.*S2P.*ws3+( ...
  -1).*S1P.*S2.*S3P.*ws3+S23.*S2P.*SGP.*ws1.*ws3+(-1).*S2.*S3P.* ...
  SGP.*ws1.*ws3+(-2).*S1P.*S23.*SGP.*ws2.*ws3+S1P.*S2.*SGP.*ws3.^2+ ...
  S13.*S2P.*(S3P.*ws2+(-1).*S2P.*ws3)+(-1).*S12.*S3P.*(S3P.*ws2+(-1) ...
  .*S2P.*ws3)+(-1).*S13.*SGP.*ws2.*(S3P.*ws2+(-1).*S2P.*ws3)+S12.* ...
  SGP.*ws3.*(S3P.*ws2+(-1).*S2P.*ws3)+(-1).*(S3P.*ws2+(-1).*S2P.* ...
  ws3).^2.*((S3P.*ws2+(-1).*S2P.*ws3).^(-2).*((-1).*S1.*S2P.^2.*S3+ ...
  2.*S1.*S23.*S2P.*S3P+S12.^2.*S3P.^2+(-1).*S1.*S2.*S3P.^2+(-2).* ...
  S12.*S2P.*S3.*SGP.*ws1+2.*S12.*S23.*S3P.*SGP.*ws1+S2P.^2.*S3.*SG.* ...
  ws1.^2+(-2).*S23.*S2P.*S3P.*SG.*ws1.^2+S2.*S3P.^2.*SG.*ws1.^2+ ...
  S23.^2.*SGP.^2.*ws1.^2+(-1).*S2.*S3.*SGP.^2.*ws1.^2+2.*S1.*S2P.* ...
  S3.*SGP.*ws2+(-2).*S1.*S23.*S3P.*SGP.*ws2+(-2).*S12.*S3P.^2.*SG.* ...
  ws1.*ws2+2.*S12.*S3.*SGP.^2.*ws1.*ws2+S1.*S3P.^2.*SG.*ws2.^2+(-1) ...
  .*S1.*S3.*SGP.^2.*ws2.^2+S13.^2.*(S2P+(-1).*SGP.*ws2).^2+(-2).* ...
  S1.*S23.*S2P.*SGP.*ws3+(-2).*S12.^2.*S3P.*SGP.*ws3+2.*S1.*S2.* ...
  S3P.*SGP.*ws3+2.*S12.*S2P.*S3P.*SG.*ws1.*ws3+(-2).*S12.*S23.* ...
  SGP.^2.*ws1.*ws3+(-2).*S1.*S2P.*S3P.*SG.*ws2.*ws3+2.*S1.*S23.* ...
  SGP.^2.*ws2.*ws3+S1.*S2P.^2.*SG.*ws3.^2+S12.^2.*SGP.^2.*ws3.^2+( ...
  -1).*S1.*S2.*SGP.^2.*ws3.^2+S1P.^2.*(S23.^2+(-1).*S2.*S3+S3.*SG.* ...
  ws2.^2+(-2).*S23.*SG.*ws2.*ws3+S2.*SG.*ws3.^2)+(-2).*S1P.*(S13.*(( ...
  -1).*S2.*S3P+S3P.*SG.*ws2.^2+S23.*(S2P+(-1).*SGP.*ws2)+S2.*SGP.* ...
  ws3+(-1).*S2P.*SG.*ws2.*ws3)+S12.*((-1).*S2P.*S3+S23.*S3P+S3.* ...
  SGP.*ws2+(-1).*S23.*SGP.*ws3+(-1).*S3P.*SG.*ws2.*ws3+S2P.*SG.* ...
  ws3.^2)+ws1.*(S23.^2.*SGP+(-1).*S2.*S3.*SGP+S2P.*S3.*SG.*ws2+S2.* ...
  S3P.*SG.*ws3+(-1).*S23.*SG.*(S3P.*ws2+S2P.*ws3)))+(-2).*S13.*( ...
  S12.*(S2P+(-1).*SGP.*ws2).*(S3P+(-1).*SGP.*ws3)+ws1.*(S23.*SGP.*(( ...
  -1).*S2P+SGP.*ws2)+S2P.*SG.*((-1).*S3P.*ws2+S2P.*ws3)+S2.*SGP.*( ...
  S3P+(-1).*SGP.*ws3))))).^(1/2))")

# To move from Mathematica, call ToMatlab on the equation after installing the appropriate package
#also need to find-replace out the backslashes
#the only purpose of this function is to handle the above monstrosity
function matlab2julia(input::String = W1SOL, minLineBreak = 80)
notationOLD::Vector{Pair} =   [
    "..." => "",
    ".*" => " * ",
    ".^" => "^",
    "+" => " + ",
    "S1P" => "SGP[1]",
    "S2P" => "SGP[2]",
    "S3P" => "SGP[3]",
    "S12" => "TG.S[1,2]",
    "S23" => "TG.S[2,3]",
    "S13" => "TG.S[1,3]",
    "S1" => "TG.S[1,1]",
    "S2" => "TG.S[2,2]",
    "S3" => "TG.S[3,3]",
    "ws1" => "TG.wₛ[1]",
    "ws2" => "TG.wₛ[2]",
    "ws3" => "TG.wₛ[3]",
    "SG" => "Θ.SG",
    "SGP" => "Θ.SGP",
    "Θ.Θ" => "Θ",
    "Θ.SGP[" => "SGP[",
    ";" => "",
    "\r\n" => "",
    "\n" => "",
    "  " => " ",
    "  " => " ",
    "  " => " "
    ]

    notation::Vector{Pair} =   [
        "..." => "",
        ".*" => " * ",
        ".^" => "^",
        "+" => " + ",
        "S1P" => "CSGP[1]",
        "S2P" => "CSGP[2]",
        "S3P" => "CSGP[3]",
        "S12" => "CS[1,2]",
        "S23" => "CS[2,3]",
        "S13" => "CS[1,3]",
        "S1" => "CS[1,1]",
        "S2" => "CS[2,2]",
        "S3" => "CS[3,3]",
        "ws1" => "Cwₛ[1]",
        "ws2" => "Cwₛ[2]",
        "ws3" => "Cwₛ[3]",
        "SG" => "Θ.SG",
        "SGP" => "Θ.SGP",
        "Θ.Θ" => "Θ",
        "Θ.SGP[" => "SGP[",
        ";" => "",
        "\r\n" => "",
        "\n" => "",
        "  " => " ",
        "  " => " ",
        "  " => " "
        ]

    out::String = input
    for i ∈ 1:length(notation)
      out = replace(out, notation[i])
    end

    outarray::Vector{Char} = Vector{Char}(out)

    #insert line breaks for readability (although it won't help much)
    widthctr::Int = 0
    for i::Int ∈ 1:length(outarray)
      widthctr += 1
      if widthctr > minLineBreak && outarray[i] == ' '
        outarray[i] = '\n'
        widthctr = 0
      end
    end

    out = String(outarray)

return "\n---------------------------\n(\n$out\n)"
end

ext: jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function funk1(A)
    (a,b) = size(A)
    s = 0
    for i in 1:a
        for j in 1:b
            if rem(i+j,2) == 0
                s+=A[i,j]
            end
        end
    end
    return s
end

function funk2(A)
    (a,b) = size(A)
    s = 0
    br = 0
    for i in 1:a
        for j in 1:b
            if(rem(j,2)==1)
                if rem(A[i,j],2)==0
                    s+=A[i,j]
                    br+=1
                end
            end
        end
    end
    return s/br
end

function funk3(A)
     range10 = zeros(0)
     range100 = zeros(0)
     range1000 = zeros(0)
    n = size(A)[2]
    for i in 1:n
        if A[i]>=0 && A[i] <10
            push!(range10,A[i])
        elseif  A[i]>=10 && A[i]<100
            push!(range100,A[i])
        elseif A[i]>=100 && A[i]<1000
            push!(range1000,A[i])
        end
    end
    return range10, range100, range1000
end


text:
function latex_files(dir::String, output = "./pdfs")
    for lecture_folder in readdir(dir)
        lecture_path = joinpath(dir,lecture_folder)
        if isdir(lecture_path) && occursin("Lecture_", lecture_folder)
            println("Compiling $lecture_folder...")
            handout = lecture_folder * "-slide_handout.tex"
            annotated = lecture_folder * "-slide_annotated.tex"
            note = lecture_folder * "-notes.tex"
            cd(lecture_path)
            latex_file(output, handout)
            latex_file(output, annotated)
            latex_file(output, note)
            cd("../..")
        end
    end
end

function latex_file(output::String, latex_file::String)
    !isdir(output) && mkdir(output)

    command = `pdflatex -file-line-error -interaction=batchmode -synctex=1 -output-directory=$(output) $(latex_file)`
    run(command)
    run(command)

    # Removes intermediate files from latex compilation
    for file in readdir(output)
        if !occursin(".pdf", file)
            rm(joinpath(output,file))
        end
    end
end


text:
"""
	$(SIGNATURES)

Deploy package docs to my website.

This is to be included in each package's `make.jl`. The benefit over hosting on github is to work around issues with private repos and private registries that would arise when TravisCI is used.

# Example
```
# This gives the package directory, if called from `docs/make.jl`.
pkgDir = rstrip(normpath(@__DIR__, ".."), '/');
@assert endswith(pkgDir, "FilesLH")
deploy_docs(pkgDir);
```

# Arguments
- `pkgDir`: points to package directory. If `pkgDir = /a/b/Foo` the docs end up in `julia/Foo` on my website.
"""
function deploy_docs(pkgDir :: AbstractString;
    trialRun :: Bool = false,
    tgBaseDir :: AbstractString = "/Users/lutz/Documents/data/web/hendri54.github.io/julia/")

    pkgName = basename(pkgDir);
    srcDir = joinpath(pkgDir, "docs", "build");
    @assert isdir(srcDir);
    tgDir = joinpath(tgBaseDir, pkgName);
    cm = rsync_command(srcDir, tgDir; trialRun = trialRun,
        doDelete = true,  inclSubDirs = true);
    run(cm);
end

"""
	$(SIGNATURES)

Given a path and several base directories, return that path hanging off each of the base directories.

# Arguments
- srcPath :: String
    Absolute or relative path. If absolute, it must include one of the base directories.
"""
function common_base_dir(srcPath :: String, baseDirV :: Vector{String})
    @assert length(baseDirV) >= 1
    @assert all(isabspath.(baseDirV)) "All base directories must be absolute paths"

    # Find which baseDir we match.
    relPath = "";
    if isabspath(srcPath)
        found = false;
        for baseDir in baseDirV
            if startswith(srcPath, baseDir)
                relPath = relpath(srcPath, baseDir);
                found = true;
                break;
            end
        end
        if !found
			error("srcPath must be inside a base directory:  $srcPath")
        end
    else
        relPath = srcPath;
    end

    outDirV = similar(baseDirV);
    for (j, baseDir) in enumerate(baseDirV)
        outDirV[j] = joinpath(baseDir, relPath);
    end

	return outDirV
end

"""
	$(SIGNATURES)

Make a pair of local and remote paths, both in the same location relative to homedir().

# Arguments
- srcPath
	can be absolute path on either machine or relative to `homedir()`
"""
function remote_and_local_path(srcPath :: String;
	srcComp = defaultLocal,  tgComp = defaultRemote)

	srcComp = get_computer(srcComp);
    tgComp = get_computer(tgComp);

    outDirV = common_base_dir(srcPath,  [srcComp.homeDir, tgComp.homeDir]);
    localDir = outDirV[1];
    remoteDir = outDirV[2];

	# # Make path relative to home on whatever machine it matches.
    # if isabspath(srcPath)
    # 	if startswith(srcPath, srcComp.homeDir)
	# 		srcPath = relpath(srcPath, srcComp.homeDir);
	# 	elseif startswith(srcPath, tgComp.homeDir)
	# 		srcPath = relpath(srcPath, tgComp.homeDir);
	# 	else
	# 		error("srcPath must be inside home directory:  $srcPath")
	# 	end
	# end

    # localDir = joinpath(srcComp.homeDir, srcPath);
    # remoteDir = joinpath(tgComp.homeDir, srcPath);
	return localDir, remoteDir
end

"""
	$(SIGNATURES)

Make a directory on the remote, optionally creating all parents as well.
Creates a command of the form
`ssh lhendri@longleaf.unc.edu 'mkdir -p ~/home/abc'`.
Shows on-screen error if the directory already exists.

Generates on-screen output from ssh login.

Note that other shell commands can be run with the same syntax.

See https://unix.stackexchange.com/questions/8612/programmatically-creating-a-remote-directory-using-ssh

# Arguments
- `sshStr`: e.g. "lhendri@longleaf.unc.edu
"""
function make_remote_dir(remoteDir :: String,  sshStr :: String;
    trialRun :: Bool = false, createParents :: Bool = true)

    if createParents
        pStr = " -p";
    else
        pStr = " ";
    end
    rcCmd = `ssh $sshStr mkdir $pStr $remoteDir`;

    if trialRun
        println(rcCmd);
    else
        run(rcCmd);
    end
    return rcCmd
end

"""
	$(SIGNATURES)

Copy a file between local and remote computer using `scp`.
`scp` replaces existing files without warning.

Syntax

	`scp /path/local lhendri@longleaf.unc.edu:/path/remote`

# Arguments
- sshStr :: String
	e.g. "lhendri@longleaf.unc.edu"
"""
function remote_copy(localPath :: String, remotePath :: String, sshStr :: String,
	upDown :: Symbol; trialRun :: Bool = false)

	if upDown == :up
		rcCmd = `scp $localPath $(sshStr):$(remotePath)`
	elseif upDown == :down
		rcCmd = `scp $(sshStr):$(remotePath) $localPath`
	else
		error("Invalid direction: $upDown")
	end
	if trialRun
		println(rcCmd)
	else
		run(rcCmd)
	end
	return rcCmd
end

"""
    $(SIGNATURES)

Rsync to or from remote machine. Only for directories.

Fails if local or remote paths contain spaces.

# Arguments
- srcDirIn, tgDirIn :: String
    Source and target directories.
    Remote directory must include part required for `ssh` (e.g. `lhendri@longleaf.unc.edu:`)
- `doDelete` :: Bool
    Should orphan files be deleted on target? `False` by default.
- `trialRun` :: Bool
- `progress` :: Bool
    Show progress indicator

#ToDo: How to make the interpolation work when there are quotes around the objects?
"""
function rsync_command(srcDirIn :: String, tgDirIn :: String;
    trialRun :: Bool = false, doDelete :: Bool = false,
    doCompress :: Bool = true, verbose :: Bool = true,
    inclSubDirs :: Bool = true, progress :: Bool = false,
    excludeV = nothing)

	# This ensures that `srcDir` ends in `/`
    srcDir = joinpath(srcDirIn, "");
    @assert srcDir[end] == filesep();

    # Ensure that remote dir does NOT end in "/"
    tgDir = rstrip(tgDirIn, filesep());

    # Need to concatenate commands, not strings for interpolation to work.
    if doDelete
        deleteStr = ` --delete`;
    else
        deleteStr = ``;
    end

    # Switches
    switchV = Vector{Char}();
    trialRun  &&  push!(switchV, 'n');
    doCompress &&  push!(switchV, 'z');
    verbose &&  push!(switchV, 'v');
    inclSubDirs  &&  push!(switchV, 'a');
    progress  &&  push!(switchV, 'P');
    vecStr = *(switchV...);
    switchStr = ` -tu$vecStr `;

    exclCmd = exclude_commands(excludeV);

    cmdStr = `rsync $switchStr $deleteStr $exclCmd --exclude .git $srcDir $tgDir`;
	return cmdStr
end

exclude_commands(exclV :: Nothing) = ``;
exclude_commands(exclV :: AbstractString) = `--exclude $exclV`;
function exclude_commands(exclV)
    if isempty(exclV)
        return ``;
    end
    cmdV = [`--exclude $x` for x in exclV];
    return cat_commands(cmdV);
end

# function exclude_command(x)
#     return `--exclude $x`;
# end

# Concatenate several commands
function cat_commands(cmdV :: Vector{Cmd})
    if isempty(cmdV)
        return ``
    end
    c = cmdV[1];
    if length(cmdV) > 1
        for j = 2 : length(cmdV)
            c = `$c $(cmdV[j])`;
        end
    end
    return c
end

"""
	$(SIGNATURES)

Transfer a directory to remote using rsync.
`delete` is by default off.
"""
function rsync_dir(srcDir :: String, tgDir :: String; kwargs...)
    # trialRun :: Bool = false, doDelete :: Bool = false)

	rCmd = rsync_command(srcDir,  tgDir; kwargs...);
	println(rCmd);
	# `trialRun` is built into the command
    run(rCmd);
    return rCmd
end

"""
	$(SIGNATURES)

Upload a repo to git.
"""
function git_upload_dir(repoDir :: String; trialRun :: Bool = false)
    @assert isdir(repoDir);
    addCmd = Cmd(`git add .`, dir = repoDir);
    commitCmd = Cmd(`git commit -am "update"`, dir = repoDir);
    pushCmd = Cmd(`git push origin master`, dir = repoDir);
    success = true;

    if trialRun
        println(addCmd);
        println(commitCmd);
        println(pushCmd);
    else
        run(addCmd)
        # Commit fails if there is nothing to commit
        try
            run(commitCmd);
            run(pushCmd);
        catch
            @warn "git commit failed"
            success = false;
        end
    end
    return success
end

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


text:
"""
    function ll_twogroup_newvariant_infboost_var_trans(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate;σ = 0.16,ω = 1/180,ct_min2 = 0.445)

Log-likelihood of the case and serology data in `model` assuming that new variant is more transmissible rather than immune evading.
"""
function ll_twogroup_newvariant_infboost_var_trans(θ,model::KenyaCoVSD.CoVAreaModel,seroreversionrate;σ = 0.16,ω = 1/180,ct_min2 = 0.445,ι=1.0)
    @unpack PCR_cases,sero_cases,baseline_sero_array,PCR_array,sero_sensitivity,sero_specificity,N,M_BB,prob,α,γ,relative_testing_rate = model
    @unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    #Set variance scalers
    clustering_factor_PCR = 0.5
    M_PCR = 30.
    T = eltype(R₀)

    #Set transmission parameters
	N₁ = P_eff*N
	N₂ = N - N₁
    p = convert.(T,[R₀,α,γ,ϵ,σ,N₁,N₂,ω,schooleffect,ct_min1,ct_min2,ι])
    u0 = convert.(T,[N₁,E₀,0.,0.,0.,0.,
                     N₂,E₀,0.,0.,0.,0.,
                     0.0,0.0,
		     0.0,0.0,0.0,0.0])

    #Create new variant introduction
    function new_variant_effect!(integrator)
        integrator.p[1] *= extra_transmissibility
        integrator.u[2] += influx_exposed_new_variant #Lower SES
        integrator.u[8] += influx_exposed_new_variant #Higher SES
    end

    janendpoint = (Date(2021,1,30) - Date(2020,2,20)).value
	aprilendpoint = (Date(2021,4,30) - Date(2020,2,20)).value
    variant_cb = PresetTimeCallback([janendpoint],new_variant_effect!)

    #Sero-waning
    sero_array = vcat(baseline_sero_array[1:30],[(1-seroreversionrate)^k for k in 1:length(baseline_sero_array[31:end])])

    #Solve for daily incidence by age (not currently broken down by severity of disease)
    LL = T(0.)
    try
        sol = solve(prob, BS3();tspan = (0,aprilendpoint),reltol = 1e-3,
				callback = variant_cb,
				u0=u0,
				p=p,
				saveat = 1,
				verbose = true,
			        isoutofdomain=(u,p,t) -> any(x -> x < 0, u))

        ι₁ = diff(sol[:C₁])
        ι₂ = diff(sol[:C₂])
        ι_sero₁ = diff(sol[:C_sero₁])
        ι_sero₂ = diff(sol[:C_sero₂])
        PCR₁ = KenyaCoVSD.simple_conv(ι₁,PCR_array)
        PCR₂ = KenyaCoVSD.simple_conv(ι₂,PCR_array)
    	sero₁ = sero_sensitivity.*KenyaCoVSD.simple_conv(ι_sero₁,sero_array)
	    sero₂ = sero_sensitivity.*KenyaCoVSD.simple_conv(ι_sero₂,sero_array)

        #Calculate log-likelihood for PCR testing
        for t in 55:janendpoint
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = relative_testing_rate[t]*(p_test₁*PCR₁[t] + p_test₂*PCR₂[t]) + 0.001
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = p_choose1*(χ₁*PCR₁[t]/((χ₁-1)*PCR₁[t] + sum(N₁))) + (1-p_choose1)*(χ₂*PCR₂[t]/((χ₂-1)*PCR₂[t] + sum(N₂))) + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing --- proportion postive
            end
        end

        #Calculate log-likelihood for PCR testing after new variant introduction
        for t in (janendpoint+1):min(size(PCR_cases,1),size(PCR₁,1))
            #Convert from μ,α parameterisation to p,r parameterisation for using Neg. binomial model
            μ = p_test_boost*relative_testing_rate[t]*(p_test₁*PCR₁[t] + p_test₂*PCR₂[t]) + 0.001
            σ² = μ + clustering_factor_PCR*μ^2
            p_negbin = 1 - (clustering_factor_PCR*μ^2/σ²)
            r_negbin  = 1/clustering_factor_PCR
            LL += logpdf(NegativeBinomial(r_negbin,p_negbin),PCR_cases[t,1])#likelihood contribution from PCR testing --- positive case detection

            if PCR_cases[t,2] >= 0 #Negative tests available
                p_PCR_pred = p_choose1*(χ_boost*χ₁*PCR₁[t]/((χ_boost*χ₁-1)*PCR₁[t] + sum(N₁))) + (1-p_choose1)*(χ_boost*χ₂*PCR₂[t]/((χ_boost*χ₂-1)*PCR₂[t] + sum(N₂))) + 0.001
                #Convert from p_PCR,inv_M_PCR parameterisation of Beta-Binomial to standard α,β parameterisation
                LL += logpdf(BetaBinomial(PCR_cases[t,2],p_PCR_pred*M_PCR,(1-p_PCR_pred)*M_PCR),PCR_cases[t,1])#likelihood contribution from PCR testing that day in that age group
            end
        end
        #Calculate log-likelihood contribution from serological testing
        for t in 1:min(size(sero_cases,1),size(sero₁,1))
            #Convert from the p_hat,M_BB parameterisation of the Betabinomial distribution to the α_BB,β_BB parameterisation

            #Overall serological testing
            if sero_cases[t,1,1] + sero_cases[t,2,1] > 0
                p_sero_pred  = (sero₁[t] + sero₂[t])./(N₁ + N₂)
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,1] + sero_cases[t,2,1],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,1])#Likelihood contribution from sero testing
            end
            #Lower SES serological testing
            if sero_cases[t,1,2] + sero_cases[t,2,2] > 0
                p_sero_pred  = sero₁[t]./N₁
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,2] + sero_cases[t,2,2],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,2])#Likelihood contribution from sero testing
            end
            #Higher SES serological testing
            if sero_cases[t,1,3] + sero_cases[t,2,3] > 0
                p_sero_pred  = sero₂[t]./N₂
                p_hat  = p_sero_pred + (1-p_sero_pred)*(1-sero_specificity)
                LL += logpdf(BetaBinomial(sero_cases[t,1,3] + sero_cases[t,2,3],M_BB*p_hat,M_BB*(1-p_hat)),sero_cases[t,1,3])#Likelihood contribution from sero testing
            end

        end
    catch
        LL = T(-Inf)
    end
    return LL::T
end

"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/180,num_chains = 1)

Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search starts at `q₀` in the transformed domain defined by `trans`.
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.

    Serowaning rate can be added, as well as additional chains for cross-comparison and validation.

"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/365,num_chains = 1,σ = 0.16,ω = 1/180,ι = 1.0)
        println("Starting MCMC parameter inference for $(model.areaname)")
        n = length(q₀)
        l_area(x) = model.log_likelihood(x,model,serowaningrate;σ = σ,ω = ω,ι=ι) + model.log_priors(x)
        ℓ = TransformedLogDensity(trans, l_area)#transformed log-likelihood
        ∇ℓ = LogDensityProblems.ADgradient(:ForwardDiff, ℓ)#transformed log-likelihood gradient wrt the parameters

        results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric))#,
                                # reporter = NoProgressReport())
        # results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
        #                         initialization = (q = q₀,κ=GaussianKineticEnergy(D)))

        transformed_results = TransformVariables.transform.(trans,results.chain)
        val = zeros(length(transformed_results),length(transformed_results[1]),1)

        for i = 1:size(val,1),j = 1:size(val,2)
                val[i,j,1] = transformed_results[i][j]
        end

        for chain_reps = 2:num_chains
                println("Running additional chain $(chain_reps) for $(model.areaname)")
                results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                        initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                        warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric))#,
                                        # reporter = NoProgressReport())
                transformed_results = TransformVariables.transform.(trans,results.chain)
                for i = 1:size(val,1),j = 1:size(val,2)
                        val[i,j,chain_reps] = transformed_results[i][j]
                end
        end

        chn = Chains(val,[String(k) for k in keys(transformed_results[1])])

        model.MCMC_results = KenyaCoVSD.MCMCResults(chn,
                                        [l_area(transformed_results[i]) - model.log_priors(transformed_results[i]) for i = 1:length(transformed_results)],
                                        results.tree_statistics)

        return nothing
end

"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,searchrange;serowaningrate = 1/180,num_chains = 1)

Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search has a pre-optimisation using the default adaptive differential evolution algorithm implemented by the `BlackBoxOptim` package,
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.

    serowaning rate can be added, as well as additional chains for cross-comparison and validation.
"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal;serowaningrate = 1/365,num_chains = 1,σ = 0.16,ω = 1/180,ι = 1.0)
        println("Searching for a good initial condition for $(model.areaname)")
        l_area(x) = model.log_likelihood(x,model,serowaningrate;σ = σ,ω = ω,ι=ι) + model.log_priors(x)
        f(x) = -transform_logdensity(trans, l_area,x)
        searchrange = fill((-3.,3.),TransformVariables.dimension(trans))
        res = bboptimize(f; SearchRange = searchrange,PopulationSize=1000,MaxSteps=30000)#,TraceMode = :silent)
        q₀ = best_candidate(res)
        inferparameters!(model,samples,trans,stepsize,D,q₀;serowaningrate=serowaningrate,num_chains=num_chains,σ = σ,ω = ω,ι=ι)
        return nothing
end


text:
struct FieldDataset{F, M, P}
      fields :: F
    metadata :: M
    filepath :: P
end

"""
    FieldDataset(filepath; architecture=CPU(), grid=nothing, ArrayType=array_type(architecture), backend=InMemory(), metadata_paths=["metadata"])

Returns a `Dict` containing a `FieldTimeSeries` for each field in the JLD2 file located at `filepath`.
Note that model output must have been saved with halos. The `InMemory` backend will store the data
fully in memory as a 4D multi-dimensional array while the `OnDisk` backend will lazily load field time
snapshots when the `FieldTimeSeries` is indexed linearly.

`metadata_paths` is a list of JLD2 paths to look for metadata. By default it looks in `file["metadata"]`.

A `grid` and `ArrayType` may be specified to override the grid and array type used in the JLD file.
"""
function FieldDataset(filepath; architecture=CPU(), grid=nothing, ArrayType=array_type(architecture), backend=InMemory(), metadata_paths=["metadata"])
    file = jldopen(filepath)

    field_names = keys(file["timeseries"])
    filter!(k -> k != "t", field_names)  # Time is not a field.

    ds = Dict{String, FieldTimeSeries}(
        name => FieldTimeSeries(filepath, name; architecture, grid, ArrayType, backend)
        for name in field_names
    )

    metadata = Dict(
        k => file["$mp/$k"]
        for mp in metadata_paths if haskey(file, mp)
        for k in keys(file["$mp"])
    )

    close(file)

    return FieldDataset(ds, metadata, abspath(filepath))
end

Base.getindex(fds::FieldDataset, inds...) = Base.getindex(fds.fields, inds...)

Base.show(io::IO, fds::FieldDataset) where {X, Y, Z, K, A} =
    print(io, "FieldDataset with $(length(fds.fields)) fields and $(length(fds.metadata)) metadata entries.")


text:
@testset "FNT 1D" begin
    for t in BigInt.(1:13)
        @show t
        x = [1:2^(t+1);] .|> BigInt
        @show length(x)
        g = 2 |> BigInt
        q = 2^2^t + 1 |> BigInt
        @time  @test ifnt(fnt(x, g, q), g, q) == x
    end
end

@testset "FNT and NTT coherence" begin
    # FNT is special case with faster implementation of general NTT and should
    # give the same results
    t = 2
    F_t = 2^2^t + 1 #a Fermat number
    x = [1:8;]

    @test fnt(x, 2, F_t) == ntt(x, 2, F_t)
end

@testset "FNT 1D convolution" begin
    using DSP

    t = 3
    q = 2^2^t + 1
    g = 2
    # padding inputs with zeros to get non circular convolution
    x = [1:8; zeros(Int, 8);]
    h = [1:8; zeros(Int, 8);]
    X = fnt(x, g, q)
    H = fnt(h, g, q)
    Y = mod.(X .* H, q)
    y = ifnt(Y, g, q)

    @test y[1:15] == DSP.conv(1:8, 1:8)
end

@testset "FNT 2D" begin
    t = 2
    q = 2^2^t + 1
    g = 2
    x = mod.(rand(Int, 8, 8), q)
    @test ifnt(fnt(x, g, q), g, q) == x

    t = 4
    q = 2^2^t + 1
    g = 314
    x = mod.(rand(Int, 512, 512), q)
    @test ifnt(fnt(x, g, q), g, q) == x
end

@testset "FNT 2D convolution" begin
    using DSP

    t = 3
    q = 2^2^t + 1
    g = 2
    x = mod.(rand(Int, 8, 8), 5)
    x_padded = zeros(Int64, 16, 16)
    x_padded[1:8, 1:8] = x
    h = mod.(rand(Int, 8, 8), 3)
    h_padded = zeros(Int64, 16, 16)
    h_padded[1:8, 1:8] = h
    X = fnt(x_padded, g, q)
    H = fnt(h_padded, g, q)
    Y = mod.(X .* H, q)
    y = ifnt(Y, g, q)

    @test y[1:15, 1:15] == DSP.conv(x, h)
end

@testset "FNT in-place" begin
    t = 2
    q = 2^2^t + 1
    g = 2

    x = mod.(rand(Int, 8), q)
    y = fnt(x, g, q)
    @test x != y
    fnt!(x, g, q)
    @test x == y

    x = mod.(rand(Int, 8, 8), q)
    y = fnt(x, g, q)
    @test x != y
    fnt!(x, g, q)
    @test x == y

end

@testset "FNT fails for non fermat modulus" begin
    (g, q, n) = (2, 31, 5)
    x = [1:n;]
    @test_throws AssertionError fnt(x, g, q)

    (g, q, n) = (2, 33, 10)
    x = [1:n;]
    @test_throws AssertionError fnt(x, g, q)
end

@testset "isfermat() tests" begin
    known = map(n->2^2^n+1, [0:4;])
    expected = map(v -> v in known, [1:maximum(known);])
    actual = isfermat.([1:maximum(known);])

    @test expected == actual
end

@testset "modfermat() tests" begin

    for t in 0:3
        q = 2^2^t+1
        x = 0:(q-1)^2
        @test mod.(x, q) == modfermat.(x, q)
    end

    for t in 0:10
        q = BigInt(2)^2^t+1
        limit = (q-1)^2
        x = mod.(rand(0:limit, 1000), (q-1)^2)
        @test mod.(x, q) == modfermat.(x, q)
    end
end