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

using Markdown
using InteractiveUtils

# ╔═╡ 1bcf25d8-71c2-11eb-05ca-6fe7c4e1b24f
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)
end

# ╔═╡ 1f7ea604-71c2-11eb-028a-a166f6f4dc5d
begin
	using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	using Random
	using Test
	using Printf
end

# ╔═╡ f377f6d4-71c9-11eb-0249-f3d63e87062b
begin
	using LinearAlgebra

	function feature_derivative(errors::Vector{T},
			feature::Vector{T}) where {T <: Real}
    	return 2 * dot(feature, errors) # ≡ 2 * feature' * errors
	end
end

# ╔═╡ 9dcffb1c-71c1-11eb-2402-39eaa1f1b5e1
md"""
## Week 2: Multiple Regression (Gradient Descent)
"""

# ╔═╡ c39f9348-71c1-11eb-1691-1d28c6b1d916
md"""

In this notebook we will cover estimating multiple regression weights via gradient descent, performing the following:

 - Add a constant column of 1's to a DataFrame to account for the intercept
 - Write a `predict_output()` function
 - Write a function to compute the derivative of the regression weights with respect to a single feature
 - Write gradient descent function to compute the regression weights given an initial weight vector, step size and tolerance.
 - Use the gradient descent function to estimate regression weights for multiple features
"""

# ╔═╡ 3e5d907e-71c2-11eb-2eb3-4375fb21b884
md"""
### Load in house sales data
"""

# ╔═╡ 4c736078-71c2-11eb-225b-816745199225
sales = CSV.File("../../ML_UW_Spec/C02/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ 58e6255a-71c2-11eb-2094-9dbe7fb64151
first(sales, 3)

# ╔═╡ 64c97662-71c2-11eb-2894-2d544642b136
typeof(sales.bedrooms)

# ╔═╡ 80fbf0f8-71c2-11eb-265a-fb539d764180
typeof(sales.bathrooms)

# ╔═╡ 35b5d05a-71cc-11eb-15a9-852c5a534fab
md"""
### Convert to Julia Matrix/Vector
"""

# ╔═╡ a2928dc6-71c2-11eb-3724-29dc1448a0f3
md"""
Now we will write a function that will accept a DataFrame, a list of feature names (e.g. [:sqft_living, :bedrooms') and a target feature e.g. (:price) and will return two things:

 - A matrix whose columns are the desired features plus a constant column (this is how we create an 'intercept')
 - An array containing the values of the output

With this in mind, let's write the `get_data` function:

"""

# ╔═╡ e4da992e-71c2-11eb-36b8-93dd614051c1
function get_data(df, features, output)
	df[:, :constant] .= 1.0 # df.constant = fill(1.0, size(df, 1))
	features = [:constant, features...]
	X_matrix = convert(Matrix, select(df, features)) # to get a matrix
	y = df[!, output]                                # => to get a vector
	return (X_matrix, y)
end

# ╔═╡ d62e0368-71c2-11eb-16f7-5fda56951d21
begin
	(ex_features, ex_output) = get_data(sales, [:sqft_living], :price)
	with_terminal() do
		println(ex_features[1:5, :], " / ", typeof(ex_features))
		println(ex_output[1], " / ", typeof(ex_output))
	end
end

# ╔═╡ 498b3458-71cc-11eb-3b55-8381d58d04bb
md"""
### Predicting output given regression weights
"""

# ╔═╡ daf8aefc-71c7-11eb-04fa-e71d7b2cde21
md"""
The predictions from all the observations are just the dot product between the features matrix (on the left) and the weights vector (on the right).

With this in mind write the following `predict_output` function to compute the predictions given the feature matrix (X) and the weights:
"""

# ╔═╡ 02f0c99e-71c5-11eb-1fe1-b3f820998d54
function predict_output(X::Matrix{T}, weights::Vector{T}) where {T <: Real}
    # assume feature_matrix is a matrix containing the features as columns
	# and weights is a corresponding array
    X * weights
end

# ╔═╡ 9cd636f2-71c5-11eb-2957-afbaf437aaa5
begin
	my_weights = Float64[1., 1.]
	test_predictions = predict_output(ex_features, my_weights)

	@test test_predictions[1] == 1431.0  # should be 1431.0
	@test test_predictions[2] == 2951.0  # should be 2951.0
end

# ╔═╡ 21635b6e-71c8-11eb-2272-edec98c58339
md"""
### Computing the derivative

We are now going to compute the derivative of the regression cost function. Recall that the cost function is the sum over the data points of the squared difference between an observed output and a predicted output.

Since the derivative of a sum is the sum of the derivatives we can compute the derivative for a single data point and then sum over data points. We can write the squared difference between the observed output and predicted output for a single point as follows:

$$(w[0]\times[CONST] + w[1]\times[feature_1] + ... + w[i]\times[feature_i] + ... +  w[k]\times[feature_k] - output)^2$$

Where we have $k$ features and a constant. So the derivative with respect to weight $w[i]$ by the chain rule is:

$$2\times(w[0]\times[CONST] + w[1]\times[feature_1] + ... + w[i]\times[feature_i] + ... +  w[k]\times[feature_k] - output)\times[feature_i]$$

The term inside the parenthesis is just the error (difference between prediction and output). So we can re-write this as:

$$2\times error\times[feature_i]$$

That is, the derivative for the weight for feature $i$ is:
  - the sum (over data points) of 2 × the product of the error and the feature itself.
    In the case of the constant then this is just twice the sum of the errors!

Recall that twice the sum of the product of two vectors is just:
  - twice the dot product of the two vectors.
    Therefore the derivative for the weight for $feature_i$ is just two times the dot product between the values of $feature_i$ and the current errors.

With this in mind let's write the following derivative function which computes the derivative of the weight given the value of the feature (over all data points) and the errors (over all data points).

`doc ref. https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/`
"""

# ╔═╡ 1301a4d2-71ca-11eb-339c-d1f747a3ffe7
begin
	(nex_features, nex_output) = get_data(sales, [:sqft_living], :price)
	ya_weights = Float64[0., 0.] # this makes all the predictions 0
	test_preds = predict_output(nex_features, ya_weights)

	# just like SFrames 2 numpy arrays can be elementwise subtracted with '-':
	errors = test_preds .- nex_output # prediction errors in this case is just the -example_output
	ya_feature = nex_features[:, 1]   # let's compute the derivative with respect to 'constant', the ":" indicates "all rows"
	der = feature_derivative(errors, ya_feature)
	alt_der = -sum(nex_output) * 2

	@test der == alt_der  # should be the same as derivative
end

# ╔═╡ 49ce0c84-71cb-11eb-3947-6141a9498ffb
with_terminal() do
	@printf("derivative: %2.5e / alt_derivative: %2.5e\n", der, alt_der)
end

# ╔═╡ 239ff6ca-71cc-11eb-341a-df1453d60a15
md"""
### Gradient Descent [GD]

Now we will write a function that performs a gradient descent. The basic premise is simple. Given a starting point we update the current weights by moving in the negative gradient direction. Recall that the gradient is the direction of increase and therefore the negative gradient is the direction of decrease and we're trying to minimize a cost function.

The amount by which we move in the negative gradient direction is called the *step size* denoted by η. We stop when we are 'sufficiently close' to the optimum. We define this by requiring that the magnitude (length) of the gradient vector to be smaller than a fixed *tolerance* denoted by ϵ.

With this in mind, wriet the following gradient descent function below using the derivative function above. For each step in the gradient descent we update the weight for each feature befofe computing our stopping criteria
"""

# ╔═╡ 958f81a6-71cc-11eb-13e8-61dddeb014a7
function regression_gradient_descent(f_matrix::Matrix{T},
		output::Vector{T}, init_weights::Vector{T}, η::T, ϵ::T) where {T <: Real}
    weights = init_weights # make sure it's a vector

    while true
		preds = predict_output(f_matrix, weights) # compute the predictions
        errors = preds .- output                  # compute the errors
        ∇_sum_squares = 0.0                       # init the gradient sum of squares

		## Update the weights
        for ix ∈ 1:length(weights) # loop over each weight
            ## Recall that feature_matrix[:, ix] is the feature column
		    ## associated with weights[i]
            deriv_ix = feature_derivative(errors, f_matrix[:, ix])
			## add squared value of derivative to ∇_sum_squares (convergence)
            ∇_sum_squares += deriv_ix * deriv_ix

			## subtract the step size times the derivative from the CURRENT weight
            weights[ix] -= η * deriv_ix
		end
        ##  compute the square-root of the ∇_sum_squares to get the ∇ magnitude:
        ∇_magnitude = √(∇_sum_squares)
        ∇_magnitude < ϵ && break
	end
    weights
end

# ╔═╡ 4454239e-71ce-11eb-295e-1755bc0549e6
md"""
### Running the Gradient Descent as Simple Regression
"""

# ╔═╡ 0623d920-71ce-11eb-0936-870685f2e788
begin
function train_test_split(df; split=0.8, seed=42, shuffled=true)
	Random.seed!(seed)
	(nr, nc) = size(df)
	nrp = round(Int, nr * split)
	row_ixes = shuffled ? shuffle(1:nr) : collect(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

sales_train, sales_test = train_test_split(sales);
end

# ╔═╡ 6a1ee460-71ce-11eb-1b49-4de942aefe9d
begin
	# let's test out the gradient descent
	s_features = [:sqft_living]
	s_out = :price
	(simple_f_matrix, out) = get_data(sales_train, s_features, s_out)

	init_weights = [-47000., 1.]
	η = 7e-12
	ϵ = 2.5e7
end

# ╔═╡ 9883c1ee-71cf-11eb-3b48-a14df232588f
gd_weights = regression_gradient_descent(simple_f_matrix, Vector{Float64}(out), init_weights, η, ϵ)

# ╔═╡ 06cf79e0-71d0-11eb-3eab-ad040213a64b
md"""
Use your newly estimated weights and your `predict_output()` function to compute the predictions on all the test data.

you will need to create a julia array of the test feature_matrix and test output first, as follows:
"""

# ╔═╡ 26199ccc-71d0-11eb-0b4e-5f2f8ebc7b76
(test_s_feature_matrix, test_out) = get_data(sales_test, s_features, s_out);

# ╔═╡ 4be11b60-71d0-11eb-1550-6d2854691b05
md"""
Now compute your predictions using `test_s_feature_matrix` and the `gd_weights` from above.
"""

# ╔═╡ 65c26930-71d0-11eb-248a-3751d927b2f6
ya_preds = predict_output(test_s_feature_matrix, gd_weights)

# ╔═╡ 749979f6-71d0-11eb-0b3a-459257d95ca9
md"""
**Quiz Question: What is the predicted price for the 1st house in the TEST data set for model 1 (round to nearest dollar)?**

"""

# ╔═╡ 898fa5b2-71d0-11eb-18c3-ad6d0dbf9bf8
begin
	pred_price_first_house = ya_preds[1]

	with_terminal() do
		@printf("price: %6.2f / price rounded to nearest dollar: %6.0f\n", pred_price_first_house, round(pred_price_first_house))
	end
end

# ╔═╡ 5f0bdea2-71d3-11eb-150c-8722eee3f23a
md"""
Now that you have the predictions on test data, compute the RSS on the test data set. Save this value for comparison later. Recall that RSS is the sum of the squared errors (difference between prediction and output).
"""

# ╔═╡ 6d6eb4b0-71d3-11eb-24f0-f9f796040473
begin
	rss_test = sum((ya_preds - test_out).^2)

	with_terminal() do
		@printf("rss on test set: %15.2f / in scientific notation: %2.4e\n", rss_test, rss_test)
	end
end

# ╔═╡ e13d4ca6-71d0-11eb-1ad7-119b4bf3c830
md"""
### Running a multiple regression

Now we will use more than one actual feature. Use the following code to produce the weights for a second model with the following parameters:

"""

# ╔═╡ ec505642-71d0-11eb-2619-0f410fee1f94
begin
	mr_model_features = [:sqft_living, :sqft_living15]
	mr_target = :price
	(mr_feature_matrix, mr_output) = get_data(sales_train, mr_model_features,
			mr_target)
	mr_init_weights = [-100000., 1., 1.]
	ηᵣ = 4e-12
	ϵᵣ = 1e9
end

# ╔═╡ 59a6c7bc-71d1-11eb-2836-71ec957bd538
md"""
Use the above parameters to estimate the model weights. Record these values for your quiz.
"""

# ╔═╡ 6a907c30-71d1-11eb-3adb-1355e084f741
weights_mr = regression_gradient_descent(mr_feature_matrix,
	Vector{Float64}(mr_output), mr_init_weights, ηᵣ, ϵᵣ)

# ╔═╡ c47ea9ec-71d1-11eb-0491-9ba85adf1e0a
md"""
Use your newly estimated weights and the predict_output function to compute the predictions on the test data.

*Don't forget to create a Julia array for these features from the test set first!*
"""

# ╔═╡ da6f7c6a-71d1-11eb-381f-5f75054e14e7
begin
	(test_mr_features_matrix, test_mr_output) = get_data(sales_test,
		mr_model_features, mr_target)

	preds_mr = predict_output(test_mr_features_matrix, weights_mr)
end

# ╔═╡ 05450824-71d2-11eb-0a05-2b895313421f
md"""
**Quiz Question: What is the predicted price for the 1st house in the TEST data set for model 2 (round to nearest dollar)?**

"""

# ╔═╡ 1706d680-71d2-11eb-3cd6-4d156d964fd7
begin
	pred_mr_price_first_house = preds_mr[1]

	with_terminal() do
		@printf("price: %6.2f / price rounded to nearest dollar: %6.0f\n", pred_mr_price_first_house, round(pred_mr_price_first_house))
	end
end

# ╔═╡ 44c15f32-71d2-11eb-2232-5fd5a07b2bfc
md"""
What is the actual price for the 1st house in the test data set?
"""

# ╔═╡ 4eeb7218-71d2-11eb-1fa7-21fec4e8c846
sales_test[!, :price][1]

# ╔═╡ 84462372-71d2-11eb-1627-1fc0ccb0d650
md"""
**Quiz Question: Which estimate was closer to the true price for the 1st house on the *test* data set, model 1 or model 2?**

  - [x] Model 1
  - [ ] Model 2

"""

# ╔═╡ 8ef08466-71d2-11eb-0cbd-17d088b840e0
md"""
Now use your predictions and the output to compute the RSS for model 2 on *test* data.
"""

# ╔═╡ 9db31f5e-71d2-11eb-1527-29c0a844ce47
begin
	rss_test_mr = sum((preds_mr - test_mr_output).^2)

	with_terminal() do
		@printf("rss on test set: %15.2f / in scientific notation: %2.4e\n", rss_test_mr, rss_test_mr)
	end
end

# ╔═╡ 0d254b32-71d3-11eb-2c92-b1a2d1125c5b
md"""
**Quiz Question: Which model (1 or 2) has lowest RSS on all of the test data?**

  - [ ] model 1
  - [x] model 2
"""

# ╔═╡ Cell order:
# ╟─9dcffb1c-71c1-11eb-2402-39eaa1f1b5e1
# ╟─c39f9348-71c1-11eb-1691-1d28c6b1d916
# ╠═1bcf25d8-71c2-11eb-05ca-6fe7c4e1b24f
# ╠═1f7ea604-71c2-11eb-028a-a166f6f4dc5d
# ╟─3e5d907e-71c2-11eb-2eb3-4375fb21b884
# ╠═4c736078-71c2-11eb-225b-816745199225
# ╠═58e6255a-71c2-11eb-2094-9dbe7fb64151
# ╠═64c97662-71c2-11eb-2894-2d544642b136
# ╠═80fbf0f8-71c2-11eb-265a-fb539d764180
# ╟─35b5d05a-71cc-11eb-15a9-852c5a534fab
# ╟─a2928dc6-71c2-11eb-3724-29dc1448a0f3
# ╠═e4da992e-71c2-11eb-36b8-93dd614051c1
# ╠═d62e0368-71c2-11eb-16f7-5fda56951d21
# ╟─498b3458-71cc-11eb-3b55-8381d58d04bb
# ╟─daf8aefc-71c7-11eb-04fa-e71d7b2cde21
# ╠═02f0c99e-71c5-11eb-1fe1-b3f820998d54
# ╠═9cd636f2-71c5-11eb-2957-afbaf437aaa5
# ╟─21635b6e-71c8-11eb-2272-edec98c58339
# ╠═f377f6d4-71c9-11eb-0249-f3d63e87062b
# ╠═1301a4d2-71ca-11eb-339c-d1f747a3ffe7
# ╠═49ce0c84-71cb-11eb-3947-6141a9498ffb
# ╟─239ff6ca-71cc-11eb-341a-df1453d60a15
# ╠═958f81a6-71cc-11eb-13e8-61dddeb014a7
# ╟─4454239e-71ce-11eb-295e-1755bc0549e6
# ╠═0623d920-71ce-11eb-0936-870685f2e788
# ╠═6a1ee460-71ce-11eb-1b49-4de942aefe9d
# ╠═9883c1ee-71cf-11eb-3b48-a14df232588f
# ╟─06cf79e0-71d0-11eb-3eab-ad040213a64b
# ╠═26199ccc-71d0-11eb-0b4e-5f2f8ebc7b76
# ╟─4be11b60-71d0-11eb-1550-6d2854691b05
# ╠═65c26930-71d0-11eb-248a-3751d927b2f6
# ╟─749979f6-71d0-11eb-0b3a-459257d95ca9
# ╠═898fa5b2-71d0-11eb-18c3-ad6d0dbf9bf8
# ╠═5f0bdea2-71d3-11eb-150c-8722eee3f23a
# ╠═6d6eb4b0-71d3-11eb-24f0-f9f796040473
# ╟─e13d4ca6-71d0-11eb-1ad7-119b4bf3c830
# ╠═ec505642-71d0-11eb-2619-0f410fee1f94
# ╟─59a6c7bc-71d1-11eb-2836-71ec957bd538
# ╠═6a907c30-71d1-11eb-3adb-1355e084f741
# ╟─c47ea9ec-71d1-11eb-0491-9ba85adf1e0a
# ╠═da6f7c6a-71d1-11eb-381f-5f75054e14e7
# ╟─05450824-71d2-11eb-0a05-2b895313421f
# ╠═1706d680-71d2-11eb-3cd6-4d156d964fd7
# ╟─44c15f32-71d2-11eb-2232-5fd5a07b2bfc
# ╠═4eeb7218-71d2-11eb-1fa7-21fec4e8c846
# ╟─84462372-71d2-11eb-1627-1fc0ccb0d650
# ╟─8ef08466-71d2-11eb-0cbd-17d088b840e0
# ╠═9db31f5e-71d2-11eb-1527-29c0a844ce47
# ╟─0d254b32-71d3-11eb-2c92-b1a2d1125c5b


text:
using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays, ArrayLayouts
import LazyArrays: CachedArray, colsupport, rowsupport, LazyArrayStyle, broadcasted,
            PaddedLayout, ApplyLayout, BroadcastLayout, AddArray, LazyLayout

@testset "Lazy MemoryLayout" begin
    @testset "ApplyArray" begin
        A = [1.0 2; 3 4]
        @test eltype(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))) == Float64
        @test @inferred(MemoryLayout(typeof(AddArray(A, Fill(0, (2, 2)), Zeros(2, 2))))) ==
            ApplyLayout{typeof(+)}()
    end

    @testset "BroadcastArray" begin
        A = [1.0 2; 3 4]

        @test @inferred(MemoryLayout(typeof(BroadcastArray(+, A, Fill(0, (2, 2)), Zeros(2, 2))))) ==
            BroadcastLayout{typeof(+)}()

        @test MemoryLayout(typeof(Diagonal(BroadcastArray(exp,randn(5))))) == DiagonalLayout{LazyLayout}()
    end

    @testset "Vcat" begin
        @test @inferred(MemoryLayout(typeof(Vcat(Ones(10),Zeros(10))))) == PaddedLayout{FillLayout}()
        @test @inferred(MemoryLayout(typeof(Vcat([1.],Zeros(10))))) == PaddedLayout{DenseColumnMajor}()
    end
end
include("applytests.jl")
include("multests.jl")
include("ldivtests.jl")
include("addtests.jl")
include("setoptests.jl")
include("macrotests.jl")
include("lazymultests.jl")
include("concattests.jl")
include("broadcasttests.jl")

@testset "Kron"  begin
    A = [1,2,3]
    B = [4,5,6,7]

    @test_throws MethodError Array(@inferred(Kron(A))) == A
    K = @inferred(Kron(A,B))
    @test size(K) == (12,)
    @test size(K,1) == 12
    @test size(K,2) == 1
    @test axes(K) == (Base.OneTo(12),)
    @test [K[k] for k=1:length(K)] == Array(K) == kron(A,B)

    A = randn(3)
    K = @inferred(Kron(A,B))
    @test K isa Kron{Float64}
    @test all(K.args .=== (A,B))
    @test [K[k] for k=1:length(K)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B)

    # C = [7,8,9,10,11]
    # K = Kron(A,B,C)
    # @time [K[k] for k=1:length(K)] == Array(Kron(A,B)) == kron(A,B)

    A = randn(3,2)
    B = randn(4,6)
    K, k = Kron(A,B), kron(A,B)
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == k
    @test det(K) == 0  # kronecker of rectangular factors
    @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
    @test tr(K) ≈ tr(k)

    K, k = Kron(A,B'), kron(A,B')
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B')) == k
    @test_throws DimensionMismatch det(K)
    @test_throws DimensionMismatch tr(K)

    K, k = Kron(A',B), kron(A',B)
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B)) == k
    @test_throws DimensionMismatch det(K)
    @test_throws DimensionMismatch tr(K)

    K, k = Kron(A',B'), kron(A',B')
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B')) == k
    @test det(K) == 0  # kronecker of rectangular factors
    @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
    @test tr(K) ≈ tr(k)

    A = randn(3,3)
    B = randn(6,6)
    C = randn(2,2)
    K, k = Kron(A,B,C), kron(A,B,C)
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B,C)) == k
    @test det(K) ≈ det(k)
    @test tr(K) ≈ tr(k)

    A = randn(3,2)
    B = randn(4,6)
    K = Kron(A,B)
    C = similar(K)
    copyto!(C, K)
    @test C == kron(A,B)
    @test_throws DimensionMismatch randn(2,2) .= K

    A = rand(Int,3,2)
    K = Kron(A,B)
    @test K isa Kron{Float64}
    @test all(K.args .=== (A,B))
    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(K) == Array(Kron{Float64}(A,B)) == kron(A,B)

    K = @inferred(Kron{Float64}(Eye{Float64}(1), zeros(4)))
    @test Array(K) == zeros(4,1)

    @testset "Applied bug" begin
        A = randn(5,5)
        K = Kron(A,A)
        k = applied(kron,A,A)
        @test K[1,1] == k[1,1] == A[1,1]^2
        x = randn(5)
        K = Kron(A,x)
        k = applied(kron,A,x)
        @test K[1,1] == k[1,1] == A[1,1]*x[1]
        K = Kron(x,x)
        k = applied(kron,x,x)
        @test K[1] == k[1] == x[1]^2
        K = Kron(A)
        k = applied(kron,A)
        @test K[1,1] == k[1,1] == A[1,1]
        K = Kron(x)
        k = applied(kron,x)
        @test K[1] == k[1] == x[1]
    end

    @testset "triple vector" begin
        x = randn(5)
        y = randn(6)
        z = randn(4)
        K = Kron(x,y,z)
        @test K[1] == K[1,1] == x[1]y[1]z[1]
        @test K == kron(x,y,z)
    end
end

@testset "Cache" begin
    A = 1:10
    C = cache(A)
    @test size(C) == (10,)
    @test axes(C) == (Base.OneTo(10),)
    @test all(Vector(C) .=== Vector(A))
    @test cache(C) isa CachedArray{Int,1,Vector{Int},UnitRange{Int}}
    C2 = cache(C)
    @test C2.data !== C.data
    @test C[:] == C[Base.Slice(Base.OneTo(10))] == C

    @test cache(A)[2,1] == 2
    @test_throws BoundsError cache(A)[2,2]

    A = reshape(1:10^2, 10,10)
    C = cache(A)
    @test size(C) == (10,10)
    @test axes(C) == (Base.OneTo(10),Base.OneTo(10))
    @test all(Array(C) .=== Array(A))
    @test cache(C) isa CachedArray{Int,2,Matrix{Int},typeof(A)}
    C2 = cache(C)
    @test C2.data !== C.data

    A = reshape(1:10^3, 10,10,10)
    C = cache(A)
    @test size(C) == (10,10,10)
    @test axes(C) == (Base.OneTo(10),Base.OneTo(10),Base.OneTo(10))
    @test all(Array(C) .=== Array(A))

    A = reshape(1:10^3, 10,10,10)
    C = cache(A)
    LazyArrays.resizedata!(C,5,5,5)
    LazyArrays.resizedata!(C,8,8,8)
    @test all(C.data .=== Array(A)[1:8,1:8,1:8])

    @test C[1:3,1,:] == A[1:3,1,:]

    @testset "Matrix cache" begin
        A = collect(1:5)
        C = cache(A)
        @test C isa Vector{Int}
        C[1] = 2
        @test A[1] ≠ 2

        A = cache(Matrix(reshape(1:6,2,3)))
        C = cache(A)
        @test C isa Matrix{Int}
        C[1,1] = 2
        @test A[1,1] ≠ 2
        C[1] = 3
        @test A[1,1] ≠ 3
    end

    @testset "setindex!" begin
        A = (1:5)
        C = cache(A)
        C[1] = 2
        @test C[1] == 2

        A = cache(reshape(1:6,2,3))
        C = cache(A)
        C[1,1] = 2
        @test A[1,1] ≠ 2
        @test C[1,1] == 2
        C[1] = 3
        @test C[1,1] == 3

        @test_throws BoundsError C[3,1]
        @test_throws BoundsError C[7]
    end

    @testset "colsupport past size" begin
        C = cache(Zeros(5,5)); C[5,1];
        @test colsupport(C,1) == Base.OneTo(5)
        @test colsupport(C,3) == 1:0
        @test rowsupport(C,1) == Base.OneTo(1)
    end

    @testset "broadcast" begin
        x = CachedArray([1,2,3],1:8);
        y = 1:8;
        f = (x,y) -> cos(x*y)
        @test f.(x,y) isa CachedArray
        @test @inferred(broadcast(f,x,y)) == f.(Vector(x), Vector(y))

        @test (x + y) isa CachedArray
        @test (x + y).array isa AbstractRange
        @test (x + y) == Vector(x) + Vector(y)

        z = CachedArray([1,4],Zeros{Int}(8));
        @test (x .+ z) isa CachedArray
        @test (x + z) isa CachedArray
        @test Vector( x .+ z) == Vector( x + z) == Vector(x) + Vector(z)

        # Lazy mixed with Static treats as Lazy
        s = SVector(1,2,3,4,5,6,7,8)
        @test f.(x , s) isa CachedArray
        @test f.(x , s) == f.(Vector(x), Vector(s))
    end

    @testset "padded CachedVector getindex" begin
        v = CachedArray([1,2,3],Zeros{Int}(1000))
        @test v[3:100] == [3; zeros(Int,97)]
        @test v[3:100] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test v[4:end] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test all(iszero,v[4:end])
        @test isempty(v[4:0])
        v = CachedArray([1,2,3],Fill{Int}(1,1000))
        @test v[3:100] == [3; ones(97)]
    end
end

@testset "Diff and Cumsum" begin
    x = Vcat([3,4], [1,1,1,1,1], 3)
    y = @inferred(cumsum(x))
    @test y isa Vcat
    @test y == cumsum(Vector(x)) == Cumsum(x)
    @test @inferred(diff(x)) == diff(Vector(x)) == Diff(x)
    @test diff(x) isa Vcat

    @test sum(x) == sum(Vector(x)) == last(y)
    @test cumsum(Vcat(4)) === Vcat(4)
    @test diff(Vcat(4)) === Vcat{Int}()

    A = randn(3,4)
    @test Cumsum(A; dims=1) == cumsum(A; dims=1)
    @test Cumsum(A; dims=2) == cumsum(A; dims=2)
    @test Diff(A; dims=1) == diff(A; dims=1)
    @test Diff(A; dims=2) == diff(A; dims=2)

    @test_broken cumsum(Vcat(Int[], 1:5)) == cumsum(1:5)

    @test cumsum(BroadcastArray(exp, 1:10)) === Cumsum(BroadcastArray(exp, 1:10))
    @test cumsum(ApplyArray(+, 1:10)) === Cumsum(ApplyArray(+, 1:10))
end

@testset "col/rowsupport" begin
    A = randn(5,6)
    @test rowsupport(A,1) === Base.OneTo(6)
    @test colsupport(A,1) === Base.OneTo(5)
    D = Diagonal(randn(5))
    @test rowsupport(D,3) === colsupport(D,3) === 3:3
    Z = Zeros(5)
    @test rowsupport(Z,1) === colsupport(Z,1) === 1:0
    @test_broken cache(D)
    C = cache(Array,D);
    @test colsupport(C,2) === 2:2
    @test colsupport(C,1) === 1:1
    @test colsupport(cache(Zeros(5,5)),1) == 1:0
    C = cache(Zeros(5));
    @test colsupport(C,1) == 1:0
    C[3] = 1
    @test colsupport(C,1) == 1:3

    LazyArrays.zero!(C)
    @test colsupport(C,1) == 1:3
    @test C == zeros(5)

    # bug from BandedMartrices.jl
    @test LazyArrays.convexunion(7:10,9:8) == LazyArrays.convexunion(9:8,7:10) == 7:10
end

@testset "triu/tril" begin
    A = ApplyArray(triu,randn(2,2))
    @test A isa ApplyArray{Float64}
    @test A == triu(A.args[1])
    A = ApplyArray(tril,randn(2,2))
    @test A isa ApplyArray{Float64}
    @test A == tril(A.args[1])
    A = ApplyArray(triu,randn(2,2),1)
    @test A isa ApplyArray{Float64}
    @test A == triu(A.args[1],1)
    A = ApplyArray(tril,randn(2,2),-1)
    @test A isa ApplyArray{Float64}
    @test A == tril(A.args[1],-1)
end

@testset "BroadcastArray" begin
    bc = broadcasted(exp,[1,2,3])
    v = BroadcastArray(exp, [1,2,3])
    @test BroadcastArray(bc) == BroadcastVector(bc) == BroadcastVector{Float64,typeof(exp),typeof(bc.args)}(bc) ==
        v == BroadcastVector(exp, [1,2,3]) == exp.([1,2,3])

    Base.IndexStyle(typeof(BroadcastVector(exp, [1,2,3]))) == IndexLinear()

    bc = broadcasted(exp,[1 2; 3 4])
    M = BroadcastArray(exp, [1 2; 3 4])
    @test BroadcastArray(bc) == BroadcastMatrix(bc) == BroadcastMatrix{Float64,typeof(exp),typeof(bc.args)}(bc) ==
        M == BroadcastMatrix(BroadcastMatrix(bc)) == BroadcastMatrix(exp,[1 2; 3 4]) == exp.([1 2; 3 4])

    @test exp.(v') isa BroadcastMatrix
    @test exp.(transpose(v)) isa BroadcastMatrix
    @test exp.(M') isa BroadcastMatrix
    @test exp.(transpose(M)) isa BroadcastMatrix
end

ext: jl


text:
##### Beginning of file

import LibGit2
import Pkg

function maketempdir()::String
    path::String = mktempdir()
    atexit(() -> _force_remove_path(path))
    return path
end

function _force_remove_path(path)::Nothing
    try
        rm(path; force = true, recursive = true)
    catch
    end
    return nothing
end

function pairwise_equality(
        x::AbstractVector,
        y::AbstractVector,
        )::Matrix{Bool}
    length_x::Int = length(x)
    length_y::Int = length(y)
    result_matrix::Matrix{Bool} = Matrix{Bool}(undef, length_x, length_y)
    for i = 1:length(x)
        for j = 1:length(y)
            result_matrix[i,j] = x[i] == y[j]
        end
    end
    return result_matrix
end

any_pairwise_equality(x, y)::Bool = any(pairwise_equality(x, y))

function version_string_equality(x::AbstractString,y::AbstractString,)::Bool
    x_converted::String = convert(String, x)
    y_converted::String = convert(String, y)
    result = version_string_equality(x_converted, y_converted,)
    return result
end

function version_string_equality(x::String, y::String,)::Bool
    x_stripped::String = strip(x)
    x_stripped_lowercase::String = lowercase(x_stripped)
    all_x = [
        x_stripped,
        x_stripped_lowercase,
        VersionNumber(x_stripped),
        VersionNumber(x_stripped_lowercase),
        strip(string(VersionNumber(x_stripped))),
        strip(string(VersionNumber(x_stripped_lowercase))),
        VersionNumber(strip(string(VersionNumber(x_stripped)))),
        VersionNumber(strip(string(VersionNumber(x_stripped_lowercase)))),
        ]
    y_stripped::String = strip(y)
    y_stripped_lowercase::String = lowercase(y_stripped)
    all_y = [
        y_stripped,
        y_stripped_lowercase,
        VersionNumber(y_stripped),
        VersionNumber(y_stripped_lowercase),
        strip(string(VersionNumber(y_stripped))),
        strip(string(VersionNumber(y_stripped_lowercase))),
        VersionNumber(strip(string(VersionNumber(y_stripped)))),
        VersionNumber(strip(string(VersionNumber(y_stripped_lowercase)))),
        ]
    result::Bool = any_pairwise_equality(
        all_x,
        all_y,
        )
    return result
end

function is_valid_version_string(x::AbstractString)::Bool
    x_converted::String = convert(String, x,)
    result = is_valid_version_string(x_converted)
    return result
end

function is_valid_version_string(x::String)::Bool
    x_stripped::String = strip(x)
    result_stripped::Bool = try
        isa(VersionNumber(x_stripped), VersionNumber)
    catch
        false
    end
    x_stripped_lowercase::String = lowercase(x_stripped)
    result_stripped_lowercase::Bool = try
        isa(VersionNumber(x_stripped_lowercase), VersionNumber)
    catch
        false
    end
    result::Bool = result_stripped && result_stripped_lowercase
    return result
end

function compare_external_registry(
        registry_path::AbstractString,
        external_registry::Pkg.RegistrySpec;
        job::AbstractString = ENV["JOB"],
        config_file::AbstractString = joinpath(
            registry_path,
            "ci",
            "ci.toml",
            ),
        )::Nothing
    original_directory::String = pwd()
    this_job_interval::AbstractInterval = _construct_interval(
        convert(String, strip(job))
        )
    configuration::Dict{String,Any} = Pkg.TOML.parsefile(config_file)
    registry_configuration_myregistry = Pkg.TOML.parsefile(
        joinpath(registry_path,"Registry.toml",)
        )
    name_to_path_myregistry = Dict{String, String}()
    for pair in registry_configuration_myregistry["packages"]
        name = pair[2]["name"]
        path = pair[2]["path"]
        name_to_path_myregistry[name] = path
    end
    all_packages::Vector{String} = collect(keys(name_to_path_myregistry))
    sort!(all_packages)
    unique!(all_packages)
    n = length(all_packages)
    @debug("all_packages ($(n)):")
    for i = 1:n
        @debug("$(i). $(all_packages[i])")
    end
    _this_job_interval_contains_x(x) = _interval_contains_x(
        this_job_interval,x,
        )
    packages_in_this_job_interval = all_packages[
        _this_job_interval_contains_x.(all_packages)
        ]
    unique!(packages_in_this_job_interval)
    sort!(packages_in_this_job_interval)
    n = length(packages_in_this_job_interval)
    @debug("packages_in_this_job_interval ($(n)):")
    for i = 1:n
        @debug("$(i). $(packages_in_this_job_interval[i])")
    end
    my_depot::String = joinpath(maketempdir(), "depot",)
    my_environment::String = joinpath(maketempdir(), "depot",)
    rm(my_depot; force = true, recursive = true,)
    rm(my_environment;force = true,recursive = true,)
    mkpath(my_depot)
    mkpath(my_environment)
    original_depot_path = [x for x in Base.DEPOT_PATH]
    empty!(Base.DEPOT_PATH)
    pushfirst!(Base.DEPOT_PATH, my_depot,)
    Pkg.activate(my_environment)
    Pkg.Registry.add(external_registry)
    external_registry_afteradding::Pkg.Types.RegistrySpec = first(
        Pkg.Types.collect_registries()
        )
    registry_configuration_externalregistry = Pkg.TOML.parsefile(
        joinpath(external_registry_afteradding.path,"Registry.toml",)
        )
    name_to_path_externalregistry = Dict{String, String}()
    for pair in registry_configuration_externalregistry["packages"]
        name = pair[2]["name"]
        path = pair[2]["path"]
        name_to_path_externalregistry[name] = path
    end
    packages_in_external_registry = Set(
        collect(keys(name_to_path_externalregistry))
        )
    n = length(packages_in_this_job_interval)
    for i = 1:n
        ### @debug(string("Checking git-tree-sha1 values for \"$(name)\" ","",))
        name = packages_in_this_job_interval[i]
        @debug(
            string(
                "Comparing git-tree-sha1 values for \"$(name)\" ",
                "(package $(i) of $(n))",
                )
            )
        my_path = name_to_path_myregistry[name]
        if name in packages_in_external_registry
            their_path = name_to_path_externalregistry[name]
            my_versions_toml = Pkg.TOML.parsefile(
                joinpath(
                    registry_path,
                    my_path,
                    "Versions.toml",
                    )
                )
            their_versions_toml = Pkg.TOML.parsefile(
                joinpath(
                    external_registry_afteradding.path,
                    their_path,
                    "Versions.toml",
                    )
                )
            my_version_strings = collect(keys(my_versions_toml))
            their_version_strings = collect(keys(their_versions_toml))
            for j = 1:length(my_version_strings)
                my_version = my_version_strings[j]
                my_git_tree_sha1 = strip(
                    my_versions_toml[
                        my_version]["git-tree-sha1"]
                    )
                for k = 1:length(their_version_strings)
                    their_version = their_version_strings[k]
                    if version_string_equality(my_version, their_version)
                        their_git_tree_sha1 = strip(
                            their_versions_toml[
                                their_version]["git-tree-sha1"]
                            )
                        if my_git_tree_sha1 == their_git_tree_sha1
                            @debug(
                                "good news: git-tree-sha1 values match",
                                my_version,
                                their_version,
                                my_git_tree_sha1,
                                their_git_tree_sha1,
                                )
                        else
                            @error(
                                "git-tree-sha1 mismatch",
                                my_version,
                                their_version,
                                my_git_tree_sha1,
                                their_git_tree_sha1,
                                )
                            error("git-tree-sha1 mismatch")
                        end

                    end
                end
            end
        end
    end
    rm(my_depot; force = true, recursive = true,)
    rm(my_environment;force = true,recursive = true,)
    empty!(Base.DEPOT_PATH)
    for x in original_depot_path
        push!(Base.DEPOT_PATH, x,)
    end
    unique!(Base.DEPOT_PATH)
    cd(original_directory)
    return nothing
end

function test_registry(
        registry_path::AbstractString;
        additional_registries::Vector{Pkg.RegistrySpec} = Pkg.RegistrySpec[
            Pkg.RegistrySpec(name="General"),
            ],
        job::AbstractString = ENV["JOB"],
        config_file::AbstractString = joinpath(
            registry_path,
            "ci",
            "ci.toml",
            ),
        )::Nothing
    original_directory::String = pwd()
    this_job_interval::AbstractInterval = _construct_interval(
        convert(String, strip(job))
        )
    configuration::Dict{String,Any} = Pkg.TOML.parsefile(config_file)
    all_packages::Vector{String} = String[]
    for (root, dirs, files) in walkdir(registry_path)
        for file in files
            if endswith(lowercase(strip(file)), ".toml")
                file_contents = Pkg.TOML.parsefile(
                    joinpath(root, file,)
                    )
                if lowercase(strip(file)) == "registry.toml"
                    append!(
                        all_packages,
                        [x["name"] for x in
                            collect(values(file_contents["packages"]))],
                        )
                end
                if lowercase(strip(file)) == "package.toml"
                    push!(
                        all_packages,
                        file_contents["name"],
                        )
                end
                if lowercase(strip(file)) == "versions.toml"
                    for version_number in keys(file_contents)
                        if !is_valid_version_string(version_number)
                            @error(
                                "Not a valid version number",
                                version_number,
                                root,
                                file,
                                )
                            error("Not a valid version number")
                        end
                    end
                end
            end
        end
    end
    unique!(all_packages)
    sort!(all_packages)
    n = length(all_packages)
    @debug("all_packages ($(n)):")
    for i = 1:n
        @debug("$(i). $(all_packages[i])")
    end
    clone_ignore::Vector{String} = configuration["clone"]["ignore"]
    unique!(clone_ignore)
    sort!(clone_ignore)
    n = length(clone_ignore)
    @debug("clone_ignore ($(n)):")
    for i = 1:n
        @debug("$(i). $(clone_ignore[i])")
    end
    packages_to_clone::Vector{String} = strip.(
        setdiff(all_packages, clone_ignore,)
        )
    unique!(packages_to_clone)
    sort!(packages_to_clone)
    n = length(packages_to_clone)
    @debug("packages_to_clone ($(n)):")
    for i = 1:n
        @debug("$(i). $(packages_to_clone[i])")
    end
    _this_job_interval_contains_x(x) = _interval_contains_x(
        this_job_interval,
        x,
        )
    packages_to_clone_in_this_job_interval = packages_to_clone[
        _this_job_interval_contains_x.(
            packages_to_clone
            )
        ]
    unique!(packages_to_clone_in_this_job_interval)
    sort!(packages_to_clone_in_this_job_interval)
    n = length(packages_to_clone_in_this_job_interval)
    @debug("packages_to_clone_in_this_job_interval ($(n)):")
    for i = 1:n
        @debug("$(i). $(packages_to_clone_in_this_job_interval[i])")
    end
    my_depot::String = joinpath(maketempdir(), "depot",)
    my_environment::String = joinpath(maketempdir(), "depot",)
    rm(my_depot; force = true, recursive = true,)
    rm(my_environment;force = true,recursive = true,)
    mkpath(my_depot)
    mkpath(my_environment)
    original_depot_path = [x for x in Base.DEPOT_PATH]
    empty!(Base.DEPOT_PATH)
    pushfirst!(Base.DEPOT_PATH, my_depot,)
    Pkg.activate(my_environment)
    Pkg.Registry.add(Pkg.RegistrySpec(path=registry_path,))
    for additional_registry in additional_registries
        Pkg.Registry.add(additional_registry)
        Pkg.Registry.update(additional_registry)
    end
    n = length(packages_to_clone_in_this_job_interval)
    for i = 1:n
        name = packages_to_clone_in_this_job_interval[i]
        @debug("Adding \"$(name)\" (package $(i) of $(n))")
        rm(my_environment;force = true,recursive = true,)
        mkpath(my_environment)
        Pkg.add(name)
        Pkg.build(name)
    end
    rm(my_depot; force = true, recursive = true,)
    rm(my_environment;force = true,recursive = true,)
    empty!(Base.DEPOT_PATH)
    for x in original_depot_path
        push!(Base.DEPOT_PATH, x,)
    end
    unique!(Base.DEPOT_PATH)
    registry_configuration = Pkg.TOML.parsefile(
        joinpath(registry_path,"Registry.toml",)
        )
    name_to_path = Dict{String, String}()
    for pair in registry_configuration["packages"]
        name = pair[2]["name"]
        path = pair[2]["path"]
        name_to_path[name] = path
    end
    n = length(packages_to_clone_in_this_job_interval)
    for i = 1:n
        name = packages_to_clone_in_this_job_interval[i]
        @debug(
            string(
                "Checking git-tree-sha1 values for \"$(name)\" ",
                "(package $(i) of $(n))",
                )
            )
        path = name_to_path[name]
        previous_directory = pwd()
        package_configuration = Pkg.TOML.parsefile(
            joinpath(path,"Package.toml")
            )
        versions_configuration = Pkg.TOML.parsefile(
            joinpath(path,"Versions.toml")
            )
        git_tree_sha1_list = String[]
        for version in keys(versions_configuration)
            push!(
                git_tree_sha1_list,
                versions_configuration[version]["git-tree-sha1"],
                )
        end
        repo_url = package_configuration["repo"]
        tmp_repo_clone_path = maketempdir()
        Base.shred!(LibGit2.CachedCredentials()) do creds
            LibGit2.with(
                Pkg.GitTools.clone(
                    repo_url,
                    tmp_repo_clone_path;
                    header = "git-repo from $(repr(repo_url))",
                    credentials = creds,
                    )
                ) do repo
            end
        end
        cd(tmp_repo_clone_path)
        for git_tree_sha1_value in git_tree_sha1_list
            cat_file_type = lowercase(
                strip(
                    read(
                        `git cat-file -t $(git_tree_sha1_value)`,
                        String,
                        )
                    )
                )
            if cat_file_type != "tree"
                @debug(
                    "git_tree_sha1 does not correspond to a tree",
                    name,
                    path,
                    repo_url,
                    git_tree_sha1_value,
                    cat_file_type,
                    )
                error("git_tree_sha1 does not correspond to a tree")
            end
        end
        cd(previous_directory)
        rm(
            tmp_repo_clone_path;
            force = true,
            recursive = true,
            )
    end
    cd(original_directory)
    return nothing
end

##### End of file


text:
function central_diff(Qbase, Qcon, cellxmax, cellymax, cellzmax, mu, lambda, vecAx, vecAy, vecAz, specific_heat_ratio, volume, Rd)
    E_vis_hat = zeros(cellxmax+1,   cellymax,   cellzmax, 5)
    F_vis_hat = zeros(  cellxmax, cellymax+1,   cellzmax, 5)
    G_vis_hat = zeros(  cellxmax,   cellymax, cellzmax+1, 5)

    for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            for k in 2:cellzmax-1

                u_av = 0.5*(Qbase[i,j,k,2]+Qbase[i-1,j,k,2])
                v_av = 0.5*(Qbase[i,j,k,3]+Qbase[i-1,j,k,3])
                w_av = 0.5*(Qbase[i,j,k,4]+Qbase[i-1,j,k,4])

                mu_av     = 0.5*(mu[i-1,j,k]+mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i-1,j,k]+1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k]+volume[i-1,j,k])

                dudxi = Qbase[i,j,k,2] - Qbase[i-1,j,k,2]
                dvdxi = Qbase[i,j,k,3] - Qbase[i-1,j,k,3]
                dwdxi = Qbase[i,j,k,4] - Qbase[i-1,j,k,4]
                dTdxi = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd) - Qbase[i-1,j,k,5]/(Qbase[i-1,j,k,1]*Rd)

                dudeta = 0.25*(Qbase[i,j+1,k,2] - Qbase[i,j-1,k,2] + Qbase[i-1,j+1,k,2] - Qbase[i-1,j-1,k,2])
                dvdeta = 0.25*(Qbase[i,j+1,k,3] - Qbase[i,j-1,k,3] + Qbase[i-1,j+1,k,3] - Qbase[i-1,j-1,k,3])
                dwdeta = 0.25*(Qbase[i,j+1,k,4] - Qbase[i,j-1,k,4] + Qbase[i-1,j+1,k,4] - Qbase[i-1,j-1,k,4])
                dTdeta = 0.25*( Qbase[i,j+1,k,5] /   (Qbase[i,j+1,k,1]*Rd) -
                                Qbase[i,j-1,k,5] /   (Qbase[i,j-1,k,1]*Rd) +
                                Qbase[i-1,j+1,k,5] / (Qbase[i-1,j+1,k,1]*Rd) -
                                Qbase[i-1,j-1,k,5] / (Qbase[i-1,j-1,k,1]*Rd) )

                dudzeta = 0.25*(Qbase[i,j,k+1,2] - Qbase[i,j,k-1,2] + Qbase[i-1,j,k+1,2] - Qbase[i-1,j,k-1,2])
                dvdzeta = 0.25*(Qbase[i,j,k+1,3] - Qbase[i,j,k-1,3] + Qbase[i-1,j,k+1,3] - Qbase[i-1,j,k-1,3])
                dwdzeta = 0.25*(Qbase[i,j,k+1,4] - Qbase[i,j,k-1,4] + Qbase[i-1,j,k+1,4] - Qbase[i-1,j,k-1,4])
                dTdzeta = 0.25*(Qbase[i,j,k+1,5] /   (Qbase[i,j,k+1,1]*Rd) -
                                Qbase[i,j,k-1,5] /   (Qbase[i,j,k-1,1]*Rd) +
                                Qbase[i-1,j,k+1,5] / (Qbase[i-1,j,k+1,1]*Rd) -
                                Qbase[i-1,j,k-1,5] / (Qbase[i-1,j,k-1,1]*Rd) )

                vecAx_xav = vecAx[i,j,k,1]
                vecAx_yav = vecAx[i,j,k,2]
                vecAx_zav = vecAx[i,j,k,3]

                vecAy_xav = 0.25*( vecAy[i,j,k,1] + vecAy[i,j+1,k,1] + vecAy[i-1,j,k,1] + vecAy[i-1,j+1,k,1] )
                vecAy_yav = 0.25*( vecAy[i,j,k,2] + vecAy[i,j+1,k,2] + vecAy[i-1,j,k,2] + vecAy[i-1,j+1,k,2] )
                vecAy_zav = 0.25*( vecAy[i,j,k,3] + vecAy[i,j+1,k,3] + vecAy[i-1,j,k,3] + vecAy[i-1,j+1,k,3] )

                vecAz_xav = 0.25*( vecAz[i,j,k,1] + vecAz[i,j,k+1,1] + vecAz[i-1,j,k,1] + vecAz[i-1,j,k+1,1] )
                vecAz_yav = 0.25*( vecAz[i,j,k,2] + vecAz[i,j,k+1,2] + vecAz[i-1,j,k,2] + vecAz[i-1,j,k+1,2] )
                vecAz_zav = 0.25*( vecAz[i,j,k,3] + vecAz[i,j,k+1,3] + vecAz[i-1,j,k,3] + vecAz[i-1,j,k+1,3] )

                dudx = vecAx_xav*dudxi + vecAy_xav*dudeta +vecAz_xav*dudzeta
                dudy = vecAx_yav*dudxi + vecAy_yav*dudeta +vecAz_yav*dudzeta
                dudz = vecAx_zav*dudxi + vecAy_zav*dudeta +vecAz_zav*dudzeta

                dvdx = vecAx_xav*dvdxi + vecAy_xav*dvdeta +vecAz_xav*dvdzeta
                dvdy = vecAx_yav*dvdxi + vecAy_yav*dvdeta +vecAz_yav*dvdzeta
                dvdz = vecAx_zav*dvdxi + vecAy_zav*dvdeta +vecAz_zav*dvdzeta

                dwdx = vecAx_xav*dwdxi + vecAy_xav*dwdeta +vecAz_xav*dwdzeta
                dwdy = vecAx_yav*dwdxi + vecAy_yav*dwdeta +vecAz_yav*dwdzeta
                dwdz = vecAx_zav*dwdxi + vecAy_zav*dwdeta +vecAz_zav*dwdzeta

                tau_xx = 2/3*mu_av*(2*dudx - dvdy - dwdz)
                tau_yy = 2/3*mu_av*(2*dvdy - dwdz - dudx)
                tau_zz = 2/3*mu_av*(2*dwdz - dudx - dvdy)

                tau_xy = mu_av*(dudy + dvdx)
                tau_yz = mu_av*(dvdz + dwdy)
                tau_zx = mu_av*(dwdx + dudz)

                dTdx = vecAx_xav*dTdxi + vecAy_xav*dTdeta + vecAz_xav*dTdzeta
                dTdy = vecAx_yav*dTdxi + vecAy_yav*dTdeta + vecAz_yav*dTdzeta
                dTdz = vecAx_zav*dTdxi + vecAy_zav*dTdeta + vecAz_zav*dTdzeta

                betax = u_av*tau_xx + v_av*tau_xy + w_av*tau_zx + lambda_av * dTdx
                betay = u_av*tau_xy + v_av*tau_yy + w_av*tau_yz + lambda_av * dTdy
                betaz = u_av*tau_zx + v_av*tau_yz + w_av*tau_zz + lambda_av * dTdz

                E_vis_hat[i,j,k,1] = 0.0
                E_vis_hat[i,j,k,2] = (vecAx_xav*tau_xx + vecAx_yav*tau_xy + vecAx_zav*tau_zx) / volume_av
                E_vis_hat[i,j,k,3] = (vecAx_xav*tau_xy + vecAx_yav*tau_yy + vecAx_zav*tau_yz) / volume_av
                E_vis_hat[i,j,k,4] = (vecAx_xav*tau_zx + vecAx_yav*tau_yz + vecAx_zav*tau_zz) / volume_av
                E_vis_hat[i,j,k,5] = (vecAx_xav*betax  + vecAx_yav*betay  + vecAx_zav*betaz )  / volume_av
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1
            for k in 2:cellzmax -1
                u_av = 0.5*(Qbase[i,j,k,2]+Qbase[i,j-1,k,2])
                v_av = 0.5*(Qbase[i,j,k,3]+Qbase[i,j-1,k,3])
                w_av = 0.5*(Qbase[i,j,k,4]+Qbase[i,j-1,k,4])

                mu_av     = 0.5*(mu[i,j-1,k]+mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i,j-1,k]+1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k]+volume[i,j-1,k])

                dudxi = 0.25*(Qbase[i+1,j,k,2] - Qbase[i-1,j,k,2] + Qbase[i+1,j-1,k,2] - Qbase[i-1,j-1,k,2])
                dvdxi = 0.25*(Qbase[i+1,j,k,3] - Qbase[i-1,j,k,3] + Qbase[i+1,j-1,k,3] - Qbase[i-1,j-1,k,3])
                dwdxi = 0.25*(Qbase[i+1,j,k,4] - Qbase[i-1,j,k,4] + Qbase[i+1,j-1,k,4] - Qbase[i-1,j-1,k,4])
                dTdxi = 0.25*(  Qbase[i+1,j,k,4]/  (Qbase[i+1,j,k,1]*Rd)   -
                                Qbase[i-1,j,k,4]/  (Qbase[i-1,j,k,1]*Rd)   +
                                Qbase[i+1,j-1,k,4]/(Qbase[i+1,j-1,k,1]*Rd) -
                                Qbase[i-1,j-1,k,4]/(Qbase[i-1,j-1,k,1]*Rd) )

                dudeta = Qbase[i,j,k,2]-Qbase[i,j-1,k,2]
                dvdeta = Qbase[i,j,k,3]-Qbase[i,j-1,k,3]
                dwdeta = Qbase[i,j,k,4]-Qbase[i,j-1,k,4]
                dTdeta = (Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd) - Qbase[i,j-1,k,5]/(Qbase[i,j-1,k,1]*Rd))

                dudzeta = 0.25*(Qbase[i,j,k+1,2] - Qbase[i,j,k-1,2] + Qbase[i,j-1,k+1,2] - Qbase[i,j-1,k-1,2])
                dvdzeta = 0.25*(Qbase[i,j,k+1,3] - Qbase[i,j,k-1,3] + Qbase[i,j-1,k+1,3] - Qbase[i,j-1,k-1,3])
                dwdzeta = 0.25*(Qbase[i,j,k+1,4] - Qbase[i,j,k-1,4] + Qbase[i,j-1,k+1,4] - Qbase[i,j-1,k-1,4])
                dTdzeta = 0.25*(Qbase[i,j,k+1,4]/  (Qbase[i,j,k+1,1]*Rd)   -
                                Qbase[i,j,k-1,4]/  (Qbase[i,j,k-1,1]*Rd)   +
                                Qbase[i,j-1,k+1,4]/(Qbase[i,j-1,k+1,1]*Rd) -
                                Qbase[i,j-1,k-1,4]/(Qbase[i,j-1,k-1,1]*Rd) )

                vecAx_xav = 0.25*( vecAx[i,j,k,1] + vecAx[i+1,j,k,1] + vecAx[i,j-1,k,1] + vecAx[i+1,j-1,k,1] )
                vecAx_yav = 0.25*( vecAx[i,j,k,2] + vecAx[i+1,j,k,2] + vecAx[i,j-1,k,2] + vecAx[i+1,j-1,k,2] )
                vecAx_zav = 0.25*( vecAx[i,j,k,3] + vecAx[i+1,j,k,3] + vecAx[i,j-1,k,3] + vecAx[i+1,j-1,k,3] )

                vecAy_xav = vecAy[i,j,k,1]
                vecAy_yav = vecAy[i,j,k,2]
                vecAy_zav = vecAy[i,j,k,3]

                vecAz_xav = 0.25*( vecAz[i,j,k,1] + vecAz[i,j,k+1,1] + vecAz[i,j-1,k,1] + vecAz[i,j-1,k+1,1] )
                vecAz_yav = 0.25*( vecAz[i,j,k,2] + vecAz[i,j,k+1,2] + vecAz[i,j-1,k,2] + vecAz[i,j-1,k+1,2] )
                vecAz_zav = 0.25*( vecAz[i,j,k,3] + vecAz[i,j,k+1,3] + vecAz[i,j-1,k,3] + vecAz[i,j-1,k+1,3] )

                dudx = vecAx_xav*dudxi + vecAy_xav*dudeta +vecAz_xav*dudzeta
                dudy = vecAx_yav*dudxi + vecAy_yav*dudeta +vecAz_yav*dudzeta
                dudz = vecAx_zav*dudxi + vecAy_zav*dudeta +vecAz_zav*dudzeta

                dvdx = vecAx_xav*dvdxi + vecAy_xav*dvdeta +vecAz_xav*dvdzeta
                dvdy = vecAx_yav*dvdxi + vecAy_yav*dvdeta +vecAz_yav*dvdzeta
                dvdz = vecAx_zav*dvdxi + vecAy_zav*dvdeta +vecAz_zav*dvdzeta

                dwdx = vecAx_xav*dwdxi + vecAy_xav*dwdeta +vecAz_xav*dwdzeta
                dwdy = vecAx_yav*dwdxi + vecAy_yav*dwdeta +vecAz_yav*dwdzeta
                dwdz = vecAx_zav*dwdxi + vecAy_zav*dwdeta +vecAz_zav*dwdzeta

                tau_xx = 2/3*mu_av*(2*dudx - dvdy - dwdz)
                tau_yy = 2/3*mu_av*(2*dvdy - dwdz - dudx)
                tau_zz = 2/3*mu_av*(2*dwdz - dudx - dvdy)

                tau_xy = mu_av*(dudy + dvdx)
                tau_yz = mu_av*(dvdz + dwdy)
                tau_zx = mu_av*(dwdx + dudz)

                dTdx = vecAx_xav*dTdxi + vecAy_xav*dTdeta + vecAz_xav*dTdzeta
                dTdy = vecAx_yav*dTdxi + vecAy_yav*dTdeta + vecAz_yav*dTdzeta
                dTdz = vecAx_zav*dTdxi + vecAy_zav*dTdeta + vecAz_zav*dTdzeta

                betax = u_av*tau_xx + v_av*tau_xy + w_av*tau_zx + lambda_av * dTdx
                betay = u_av*tau_xy + v_av*tau_yy + w_av*tau_yz + lambda_av * dTdy
                betaz = u_av*tau_zx + v_av*tau_yz + w_av*tau_zz + lambda_av * dTdz

                F_vis_hat[i,j,k,1] = 0.0
                F_vis_hat[i,j,k,2] = (vecAy_xav*tau_xx + vecAy_yav*tau_xy + vecAy_zav*tau_zx) / volume_av
                F_vis_hat[i,j,k,3] = (vecAy_xav*tau_xy + vecAy_yav*tau_yy + vecAy_zav*tau_yz) / volume_av
                F_vis_hat[i,j,k,4] = (vecAy_xav*tau_zx + vecAy_yav*tau_yz + vecAy_zav*tau_zz) / volume_av
                F_vis_hat[i,j,k,5] = (vecAy_xav*betax  + vecAy_yav*betay  + vecAy_zav*betaz ) / volume_av
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax -1
            for k in 2:cellzmax+1 -1
                u_av = 0.5*(Qbase[i,j,k,2]+Qbase[i,j,k-1,2])
                v_av = 0.5*(Qbase[i,j,k,3]+Qbase[i,j,k-1,3])
                w_av = 0.5*(Qbase[i,j,k,4]+Qbase[i,j,k-1,4])

                mu_av     = 0.5*(mu[i,j,k-1]+mu[i,j,k])
                lambda_av = (0.5*(1/lambda[i,j,k-1]+1/lambda[i,j,k]))^(-1)

                volume_av = 0.5*(volume[i,j,k]+volume[i,j,k-1])

                dudxi = 0.25*(Qbase[i+1,j,k,2] - Qbase[i-1,j,k,2] + Qbase[i+1,j,k-1,2] - Qbase[i-1,j,k-1,2])
                dvdxi = 0.25*(Qbase[i+1,j,k,3] - Qbase[i-1,j,k,3] + Qbase[i+1,j,k-1,3] - Qbase[i-1,j,k-1,3])
                dwdxi = 0.25*(Qbase[i+1,j,k,4] - Qbase[i-1,j,k,4] + Qbase[i+1,j,k-1,4] - Qbase[i-1,j,k-1,4])
                dTdxi = 0.25*(  Qbase[i+1,j,k,4]/  (Qbase[i+1,j,k,1]*Rd)   -
                                Qbase[i-1,j,k,4]/  (Qbase[i-1,j,k,1]*Rd)   +
                                Qbase[i+1,j,k-1,4]/(Qbase[i+1,j,k-1,1]*Rd) -
                                Qbase[i-1,j,k-1,4]/(Qbase[i-1,j,k-1,1]*Rd) )

                dudeta = 0.25*(Qbase[i,j+1,k,2] - Qbase[i,j-1,k,2] + Qbase[i,j+1,k-1,2] - Qbase[i,j-1,k-1,2])
                dvdeta = 0.25*(Qbase[i,j+1,k,3] - Qbase[i,j-1,k,3] + Qbase[i,j+1,k-1,3] - Qbase[i,j-1,k-1,3])
                dwdeta = 0.25*(Qbase[i,j+1,k,4] - Qbase[i,j-1,k,4] + Qbase[i,j+1,k-1,4] - Qbase[i,j-1,k-1,4])
                dTdeta = 0.25*( Qbase[i,j+1,k,4]/  (Qbase[i,j+1,k,1]*Rd)   -
                                Qbase[i,j-1,k,4]/  (Qbase[i,j-1,k,1]*Rd)   +
                                Qbase[i,j+1,k-1,4]/(Qbase[i,j+1,k-1,1]*Rd) -
                                Qbase[i,j-1,k-1,4]/(Qbase[i,j-1,k-1,1]*Rd) )

                dudzeta = Qbase[i,j,k,2]-Qbase[i,j,k-1,2]
                dvdzeta = Qbase[i,j,k,3]-Qbase[i,j,k-1,3]
                dwdzeta = Qbase[i,j,k,4]-Qbase[i,j,k-1,4]
                dTdzeta = (Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd) - Qbase[i,j,k-1,5]/(Qbase[i,j,k-1,1]*Rd))

                vecAx_xav = 0.25*( vecAx[i,j,k,1] + vecAx[i+1,j,k,1] + vecAx[i,j,k-1,1] + vecAx[i+1,j,k-1,1] )
                vecAx_yav = 0.25*( vecAx[i,j,k,2] + vecAx[i+1,j,k,2] + vecAx[i,j,k-1,2] + vecAx[i+1,j,k-1,2] )
                vecAx_zav = 0.25*( vecAx[i,j,k,3] + vecAx[i+1,j,k,3] + vecAx[i,j,k-1,3] + vecAx[i+1,j,k-1,3] )

                vecAy_xav = 0.25*( vecAy[i,j,k,1] + vecAy[i,j+1,k,1] + vecAy[i,j,k-1,1] + vecAy[i,j+1,k-1,1] )
                vecAy_yav = 0.25*( vecAy[i,j,k,2] + vecAy[i,j+1,k,2] + vecAy[i,j,k-1,2] + vecAy[i,j+1,k-1,2] )
                vecAy_zav = 0.25*( vecAy[i,j,k,3] + vecAy[i,j+1,k,3] + vecAy[i,j,k-1,3] + vecAy[i,j+1,k-1,3] )

                vecAz_xav = vecAz[i,j,k,1]
                vecAz_yav = vecAz[i,j,k,2]
                vecAz_zav = vecAz[i,j,k,3]

                dudx = vecAx_xav*dudxi + vecAy_xav*dudeta +vecAz_xav*dudzeta
                dudy = vecAx_yav*dudxi + vecAy_yav*dudeta +vecAz_yav*dudzeta
                dudz = vecAx_zav*dudxi + vecAy_zav*dudeta +vecAz_zav*dudzeta

                dvdx = vecAx_xav*dvdxi + vecAy_xav*dvdeta +vecAz_xav*dvdzeta
                dvdy = vecAx_yav*dvdxi + vecAy_yav*dvdeta +vecAz_yav*dvdzeta
                dvdz = vecAx_zav*dvdxi + vecAy_zav*dvdeta +vecAz_zav*dvdzeta

                dwdx = vecAx_xav*dwdxi + vecAy_xav*dwdeta +vecAz_xav*dwdzeta
                dwdy = vecAx_yav*dwdxi + vecAy_yav*dwdeta +vecAz_yav*dwdzeta
                dwdz = vecAx_zav*dwdxi + vecAy_zav*dwdeta +vecAz_zav*dwdzeta

                tau_xx = 2/3*mu_av*(2*dudx - dvdy - dwdz)
                tau_yy = 2/3*mu_av*(2*dvdy - dwdz - dudx)
                tau_zz = 2/3*mu_av*(2*dwdz - dudx - dvdy)

                tau_xy = mu_av*(dudy + dvdx)
                tau_yz = mu_av*(dvdz + dwdy)
                tau_zx = mu_av*(dwdx + dudz)

                dTdx = vecAx_xav*dTdxi + vecAy_xav*dTdeta + vecAz_xav*dTdzeta
                dTdy = vecAx_yav*dTdxi + vecAy_yav*dTdeta + vecAz_yav*dTdzeta
                dTdz = vecAx_zav*dTdxi + vecAy_zav*dTdeta + vecAz_zav*dTdzeta

                betax = u_av*tau_xx + v_av*tau_xy + w_av*tau_zx + lambda_av * dTdx
                betay = u_av*tau_xy + v_av*tau_yy + w_av*tau_yz + lambda_av * dTdy
                betaz = u_av*tau_zx + v_av*tau_yz + w_av*tau_zz + lambda_av * dTdz

                G_vis_hat[i,j,k,1] = 0.0
                G_vis_hat[i,j,k,2] = (vecAz_xav*tau_xx + vecAz_yav*tau_xy + vecAz_zav*tau_zx) / volume_av
                G_vis_hat[i,j,k,3] = (vecAz_xav*tau_xy + vecAz_yav*tau_yy + vecAz_zav*tau_yz) / volume_av
                G_vis_hat[i,j,k,4] = (vecAz_xav*tau_zx + vecAz_yav*tau_yz + vecAz_zav*tau_zz) / volume_av
                G_vis_hat[i,j,k,5] = (vecAz_xav*betax  + vecAz_yav*betay  + vecAz_zav*betaz ) / volume_av
            end
        end
    end

    return E_vis_hat, F_vis_hat, G_vis_hat
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
module HTTPUtils

using HTTP

function Base.Dict(req::HTTP.Request) :: Dict{String,String}
  result = Dict{String,String}()
  for (k,v) in Dict(req.headers)
    result[lowercase(string(k))] = lowercase(string(v))
  end

  result
end

end