text:
# Pre-image for Gaussian kernel
# From Kwok and Tsang, "The Pre-Image problem in kernel methods", ICML 2003
# (based on matlab code provided by authors)
# Also:
# Mika, et al. "Kernel PCA and Denoising in Feature Spaces", NIPS 1998
# and
# Teixeira et al. "KPCA Denoising and the pre-image problem revisited", DSP 2008
#
# Danny Perry (dperry@cs.utah.edu)
# May 2015

using Debug

include("Kernels.jl")

# Fixed point method to preimage
# y - projected data (trying to recover into input space)
# U - eigenvectors of K
# X - training data
# z_init - initial guess at z
# sigma - Gaussian kernel parameter
# tolerance - convergence criteria
# maxiters - max number of iters
function GaussianKernelPreImage(y,U, X, z_init, sigma, tolerance, maxiters)

	gamma = U * y # results in n x 1 vector

	z = copy(z_init)
	last_z = copy(z)
	iter = 0
	for iter=1:maxiters
		last_z = copy(z)
		diff = X - ones(size(X,1),1)*z
		diff .^= 2
		kz = gamma .* exp( -sum(diff,2) ./ sigma^2)
		if sum(kz) == 0
			println("Preimage fixed point iteration failed: initial guess too far away - reverting to initial guess.")
			return z_init # can't get any closer, returning initial guess.
		end
		z = sum(X .* (kz*ones(1,size(X,2))),1)/sum(kz)
		if norm(last_z-z)/sqrt(norm(z)*norm(last_z)) < tolerance
			break
		end
	end

	if false
		println("iters: ", iter)
		println("err: ", norm(last_z-z)/sqrt(norm(z)*norm(last_z)) , " < ", tolerance)
	end

	if iter == maxiters
		println("warning, did not converge.")
	end

	return z, iter
end

# Linear algebra approach to preimage
function GaussianKernelPreImage(distsq, X, neibsize)
	sidx = sortperm(distsq)
	sdistsq = distsq[sidx]
	XH = X[sidx[1:neibsize],:]
	Xmean = mean(XH,1)
	XH = XH - ones(neibsize,1) * Xmean # centered around neighborhood mean

	UM,SM,VM = svd(XH')
	rankM = rank(diagm(SM),1e-5)
	UM = UM[:,1:rankM]
	SM = SM[1:rankM]
	VM = VM[:,1:rankM]
	transf = UM*diagm(1./SM)*(VM'./2)

	sd0 = zeros(neibsize,1)
	ZM = diagm(SM)*VM'
	for i=1:neibsize
		sd0[i] = (ZM[:,i]'*ZM[:,i])[1]
	end
	result = transf * (vec(sd0) - vec(sdistsq[1:neibsize])) + vec(Xmean)
end

# X - training data (n x d)
# K - uncentered training data Gram (n x n)
# Ktest - uncentered test data Gram with training (n x nt)
function GaussianKernelPreImage(X,K,Ktest, neibsize, sigma)

	n  = size(K,1)
	nt = size(Ktest,2)
	d = size(X,2)

	spectrum_pct = 0
	target_dim = n
	centering = "additive"
	Kc,P,V,S,Y = KernelPCA(K, spectrum_pct, target_dim, centering)

	H = eye(n)-ones(n,n)/n # centering matrix
	HMH = H*P*P'*H
	cK = mean(K,2)
	meanK = mean(cK)

	result = zeros(nt,d)
	neibs = zeros(Int64, nt, neibsize)
	for j=1:nt
		# calculate the distance between the testing point and training points
		k_x = Ktest[:,j]
		gammaC = HMH*(k_x-cK);
		PphiNormC = ((k_x+cK)'*gammaC + meanK)[1]
		d2 = zeros(n);
		for i = 1:n
			PphiProjC = (K[i,:]*gammaC)[1]+cK[i]
			d2[i]  = -log(abs((1-PphiNormC+2*PphiProjC)/2))*(sigma*2)
		end
		result[j,:] = GaussianKernelPreImage(d2,X,neibsize);

		closestind = sortperm(d2)
		neibs[j,:] = closestind[1:neibsize]
	end

	return result,neibs
end

function GaussianKernelDenoise(X,Xtest, iters, neibsize, sigma)

	n  = size(X,1)
	nt = size(Xtest,1)
	d = size(X,2)

	K = GaussianKernel(X,X,sigma)

	spectrum_pct = 0
	target_dim = n
	centering = "additive"
	Kc,P,V,S,Y = KernelPCA(K, spectrum_pct, target_dim, centering)

	H = eye(n)-ones(n,n)/n # centering matrix
	HMH = H*P*P'*H
	cK = mean(K,2)
	meanK = mean(cK)

	result = copy(Xtest)
	neibs = zeros(Int64, nt, neibsize)
	for iter=1:iters
		for j=1:nt
			# calculate the distance between the testing point and training points
			k_x = GaussianKernel(X,Xtest[j,:],sigma)
			gammaC = HMH*(k_x-cK);
			PphiNormC = ((k_x+cK)'*gammaC + meanK)[1]
			d2 = zeros(n);
			for i = 1:n
				PphiProjC = (vec(K[i,:])'*vec(gammaC)+cK[i])[1]
				d2[i]  = -log((1-PphiNormC+2*PphiProjC)/2)*sigma
			end
			result[j,:] = GaussianKernelPreImage(d2,X,neibsize);

			closestind = sortperm(d2)
			neibs[j,:] = closestind[1:neibsize]
		end
	end

	return result,neibs
end


text:
function add(obj::AbstractFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "add", FieldMatrix, (FieldMatrix,), arg0)
end

function add_to_entry(obj::AbstractFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "addToEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function copy(obj::AbstractFieldMatrix)
    return jcall(obj, "copy", FieldMatrix, ())
end

function copy_sub_matrix(obj::AbstractFieldMatrix, arg0::Vector{jint}, arg1::Vector{jint}, arg2::Vector{Vector{FieldElement}})
    return jcall(obj, "copySubMatrix", void, (Vector{jint}, Vector{jint}, Vector{Vector{FieldElement}}), arg0, arg1, arg2)
end

function copy_sub_matrix(obj::AbstractFieldMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::Vector{Vector{FieldElement}})
    return jcall(obj, "copySubMatrix", void, (jint, jint, jint, jint, Vector{Vector{FieldElement}}), arg0, arg1, arg2, arg3, arg4)
end

function create_matrix(obj::AbstractFieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", FieldMatrix, (jint, jint), arg0, arg1)
end

function equals(obj::AbstractFieldMatrix, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_column(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getColumn", Vector{FieldElement}, (jint,), arg0)
end

function get_column_dimension(obj::AbstractFieldMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_column_matrix(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getColumnMatrix", FieldMatrix, (jint,), arg0)
end

function get_column_vector(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getColumnVector", FieldVector, (jint,), arg0)
end

function get_data(obj::AbstractFieldMatrix)
    return jcall(obj, "getData", Vector{Vector{FieldElement}}, ())
end

function get_entry(obj::AbstractFieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", FieldElement, (jint, jint), arg0, arg1)
end

function get_field(obj::AbstractFieldMatrix)
    return jcall(obj, "getField", Field, ())
end

function get_row(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getRow", Vector{FieldElement}, (jint,), arg0)
end

function get_row_dimension(obj::AbstractFieldMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function get_row_matrix(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getRowMatrix", FieldMatrix, (jint,), arg0)
end

function get_row_vector(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "getRowVector", FieldVector, (jint,), arg0)
end

function get_sub_matrix(obj::AbstractFieldMatrix, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(obj, "getSubMatrix", FieldMatrix, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function get_sub_matrix(obj::AbstractFieldMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(obj, "getSubMatrix", FieldMatrix, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function get_trace(obj::AbstractFieldMatrix)
    return jcall(obj, "getTrace", FieldElement, ())
end

function hash_code(obj::AbstractFieldMatrix)
    return jcall(obj, "hashCode", jint, ())
end

function is_square(obj::AbstractFieldMatrix)
    return jcall(obj, "isSquare", jboolean, ())
end

function map(obj::FieldMatrix, arg0::Function)
    return jcall(obj, "map", FieldMatrix, (Function,), arg0)
end

function map_to_self(obj::FieldMatrix, arg0::Function)
    return jcall(obj, "mapToSelf", FieldMatrix, (Function,), arg0)
end

function multiply(obj::AbstractFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "multiply", FieldMatrix, (FieldMatrix,), arg0)
end

function multiply_entry(obj::AbstractFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "multiplyEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function multiply_transposed(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "multiplyTransposed", FieldMatrix, (FieldMatrix,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function operate(obj::AbstractFieldMatrix, arg0::FieldVector)
    return jcall(obj, "operate", FieldVector, (FieldVector,), arg0)
end

function operate(obj::AbstractFieldMatrix, arg0::Vector{FieldElement})
    return jcall(obj, "operate", Vector{FieldElement}, (Vector{FieldElement},), arg0)
end

function power(obj::AbstractFieldMatrix, arg0::jint)
    return jcall(obj, "power", FieldMatrix, (jint,), arg0)
end

function pre_multiply(obj::AbstractFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "preMultiply", FieldMatrix, (FieldMatrix,), arg0)
end

function pre_multiply(obj::AbstractFieldMatrix, arg0::FieldVector)
    return jcall(obj, "preMultiply", FieldVector, (FieldVector,), arg0)
end

function pre_multiply(obj::AbstractFieldMatrix, arg0::Vector{FieldElement})
    return jcall(obj, "preMultiply", Vector{FieldElement}, (Vector{FieldElement},), arg0)
end

function scalar_add(obj::AbstractFieldMatrix, arg0::FieldElement)
    return jcall(obj, "scalarAdd", FieldMatrix, (FieldElement,), arg0)
end

function scalar_multiply(obj::AbstractFieldMatrix, arg0::FieldElement)
    return jcall(obj, "scalarMultiply", FieldMatrix, (FieldElement,), arg0)
end

function set_column(obj::AbstractFieldMatrix, arg0::jint, arg1::Vector{FieldElement})
    return jcall(obj, "setColumn", void, (jint, Vector{FieldElement}), arg0, arg1)
end

function set_column_matrix(obj::AbstractFieldMatrix, arg0::jint, arg1::FieldMatrix)
    return jcall(obj, "setColumnMatrix", void, (jint, FieldMatrix), arg0, arg1)
end

function set_column_vector(obj::AbstractFieldMatrix, arg0::jint, arg1::FieldVector)
    return jcall(obj, "setColumnVector", void, (jint, FieldVector), arg0, arg1)
end

function set_entry(obj::AbstractFieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "setEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function set_row(obj::AbstractFieldMatrix, arg0::jint, arg1::Vector{FieldElement})
    return jcall(obj, "setRow", void, (jint, Vector{FieldElement}), arg0, arg1)
end

function set_row_matrix(obj::AbstractFieldMatrix, arg0::jint, arg1::FieldMatrix)
    return jcall(obj, "setRowMatrix", void, (jint, FieldMatrix), arg0, arg1)
end

function set_row_vector(obj::AbstractFieldMatrix, arg0::jint, arg1::FieldVector)
    return jcall(obj, "setRowVector", void, (jint, FieldVector), arg0, arg1)
end

function set_sub_matrix(obj::AbstractFieldMatrix, arg0::Vector{Vector{FieldElement}}, arg1::jint, arg2::jint)
    return jcall(obj, "setSubMatrix", void, (Vector{Vector{FieldElement}}, jint, jint), arg0, arg1, arg2)
end

function subtract(obj::AbstractFieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "subtract", FieldMatrix, (FieldMatrix,), arg0)
end

function to_string(obj::AbstractFieldMatrix)
    return jcall(obj, "toString", JString, ())
end

function transpose(obj::AbstractFieldMatrix)
    return jcall(obj, "transpose", FieldMatrix, ())
end

function transpose_multiply(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "transposeMultiply", FieldMatrix, (FieldMatrix,), arg0)
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end

function walk_in_column_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_column_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_column_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_column_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_row_order(obj::AbstractFieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_row_order(obj::AbstractFieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

ext: jl


text:
using Test
using LanguageServer

@testset "URI conversion" begin

    if Sys.iswindows()
        @test LanguageServer.URIs2.uri2filepath(uri"file://SERVER/foo/bar") == "\\\\SERVER\\foo\\bar"
        @test LanguageServer.URIs2.uri2filepath(uri"file://wsl%24/foo/bar") == "\\\\wsl\$\\foo\\bar"
        @test LanguageServer.URIs2.uri2filepath(uri"file:///D:/foo/bar") == "d:\\foo\\bar"
        @test LanguageServer.URIs2.uri2filepath(uri"file:///foo/bar") == "\\foo\\bar"

        @test LanguageServer.filepath2uri("\\\\SERVER\\foo\\bar") == uri"file://SERVER/foo/bar"
        @test LanguageServer.filepath2uri("\\\\wsl\$\\foo\\bar") == uri"file://wsl%24/foo/bar"
        @test LanguageServer.filepath2uri("d:\\foo\\bar") == uri"file:///d%3A/foo/bar"
    else
        @test LanguageServer.URIs2.uri2filepath(uri"file:///foo/bar") == "/foo/bar"

        @test LanguageServer.filepath2uri("/foo/bar") == uri"file:///foo/bar"
    end

end

@testset "URI comparison" begin

    if Sys.iswindows()
        @test string(filepath2uri("C:\\foo\\bar")) == "file:///c%3A/foo/bar"
        @test string(filepath2uri("\\\\wsl\$\\foo\\bar")) == "file://wsl\$/foo/bar"

        @test hash(uri"file:///d:/FOO/bar") == hash(uri"file:///d%3A/FOO/bar")
        @test hash(uri"file:///c:/foo/space bar") == hash(uri"file:///c%3A/foo/space%20bar")
        @test hash(uri"file://wsl$/foo/bar") == hash(uri"file://wsl%24/foo/bar")

        @test uri"file:///d:/FOO/bar" == uri"file:///d%3A/FOO/bar"
        @test uri"file:///c:/foo/space bar" == uri"file:///c%3A/foo/space%20bar"
        @test uri"file://wsl$/foo/bar" == uri"file://wsl%24/foo/bar"
    else
        @test string(uri"file:///foo/bar") == "file:///foo/bar"

        @test hash(uri"file:///foo/bar") == hash(uri"file:///foo/bar")
        @test hash(uri"file:///foo/space bar") == hash(uri"file:///foo/space%20bar")

        @test uri"file:///foo/bar" == uri"file:///foo/bar"
        @test uri"file:///foo/space bar" == uri"file:///foo/space%20bar"
    end

end

@testset "is_in_target_dir_of_package" begin
    @test LanguageServer.is_in_target_dir_of_package(@__DIR__, "test")
    @test !LanguageServer.is_in_target_dir_of_package(pathof(LanguageServer), "test")
end


text:
using FireFlower
using Test

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


text:
@testset "1004.max-consecutive-ones-iii.jl" begin
    @test longest_ones([1,1,1,0,0,0,1,1,1,1,0], 2) == 6
    @test longest_ones([0,0,1,1,0,0,1,1,1,0,1,1,0,0,0,1,1,1,1], 3) == 10
    @test longest_ones([0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0], 6) == 17
end


text:
module Tokens

import Base.eof

export Token

include("token_kinds.jl")

iskeyword(k::Kind) = begin_keywords < k < end_keywords
isliteral(k::Kind) = begin_literal < k < end_literal
isoperator(k::Kind) = begin_ops < k < end_ops
iserror(k::Kind) = begin_errors < k < end_errors
iscontextualkeyword(k::Kind) = begin_contextual_keywords < k < end_contextual_keywords

function iswordoperator(k::Kind)
    # Keyword-like operators
    k == Tokens.IN ||
    k == Tokens.ISA ||
    k == Tokens.WHERE
end

# Create string => keyword kind
const KEYWORDS = Dict{String, Kind}()

function _add_kws()
    for k in instances(Kind)
        if iskeyword(k)
            KEYWORDS[lowercase(string(k))] = k
        end
    end
end
_add_kws()

# Error kind => description
TOKEN_ERROR_DESCRIPTION = Dict{Kind, String}(
    EOF_MULTICOMMENT => "unterminated multi-line comment #= ... =#",
    EOF_CHAR => "unterminated character literal",
    INVALID_NUMERIC_CONSTANT => "invalid numeric constant",
    INVALID_OPERATOR => "invalid operator",
    INVALID_INTERPOLATION_TERMINATOR => "interpolated variable ends with invalid character; use `\$(...)` instead",
    ERROR => "unknown error",
)

struct Token
    kind::Kind
    # Offsets into a string or buffer
    startbyte::Int # The byte where the token start in the buffer
    endbyte::Int # The byte where the token ended in the buffer
    dotop::Bool
    suffix::Bool
end
function Token(kind::Kind, startbyte::Int, endbyte::Int)
    Token(kind, startbyte, endbyte, false, false)
end
Token() = Token(ERROR, 0, 0, false, false)

const EMPTY_TOKEN = Token()

function kind(t::Token)
    isoperator(t.kind) && return OP
    iskeyword(t.kind) && return KEYWORD
    iserror(t.kind) && return ERROR
    return t.kind
end
exactkind(t::Token) = t.kind

startbyte(t::Token) = t.startbyte
endbyte(t::Token) = t.endbyte

function untokenize(t::Token, str::String)
    String(codeunits(str)[1 .+ (t.startbyte:t.endbyte)])
end

function Base.show(io::IO, t::Token)
    print(io, rpad(string(startbyte(t), "-", endbyte(t)), 11, " "))
    print(io, rpad(kind(t), 15, " "))
end

end # module


text:
# Install Julia 0.6.4 and then run
# include("install.jl") (Only once)
# to install all required packages.

# It is not advised to run this file directly. Rather, start a fresh
# Julia session and run the code between the "================================"
# delimeters. Afterwards, start a fresh Julia session and run the next block.
# This is to ensure that the number of worker processes is correct. Also,
# each block will require a couple of hours to finish.

# For customization, the parameter arguments can be changed/supplied in the
# of the L-shaped solvers in dayahead_benchmark.jl. Use `?LShapedSolver`
# in Julia prompt after `using LShapedSolvers` for a docstring on its usage.
# The parameters and default value of each algorithm is accessed by for example
# `?LShaped`.

# Serial benchmarks. Used to generate Fig. 2
# ================================ #
include("dayahead_benchmark.jl")
serial_benchmark = prepare_sbenchmark()
warmup_benchmarks!(serial_benchmark)
run_benchmarks!(serial_benchmark) # Takes a couple of hours
save_results!(serial_benchmark,"serial_benchmarks.json")
# Plotting
using Plots
pyplot()
p = plot(TimePlot(serial_benchmark)) # Creates Fig. 2
savefig(p,"serial_benchmarks.pdf")
# ================================ #

# Distributed benchmarks. Used to generate Fig. 3 and Fig. 4
# ================================ #
# No worker cores (serial)
include("dayahead_benchmark.jl")
distributed_benchmark = prepare_dbenchmark()
warmup_benchmarks!(distributed_benchmark)
run_benchmarks!(distributed_benchmark) # Takes a couple of hours
save_results!(distributed_benchmark,"scaling_1.json")
# ================================ #
# 2 worker cores
addprocs(2)
include("dayahead_benchmark.jl")
distributed_benchmark = prepare_dbenchmark()
warmup_benchmarks!(distributed_benchmark)
run_benchmarks!(distributed_benchmark) # Takes a couple of hours
save_results!(distributed_benchmark,"scaling_2.json")
# ================================ #
# 4 worker cores
addprocs(4)
include("dayahead_benchmark.jl")
distributed_benchmark = prepare_dbenchmark()
warmup_benchmarks!(distributed_benchmark)
run_benchmarks!(distributed_benchmark) # Takes a couple of hours
save_results!(distributed_benchmark,"scaling_4.json")
# ================================ #
# 8 worker cores
addprocs(8)
include("dayahead_benchmark.jl")
distributed_benchmark = prepare_dbenchmark()
warmup_benchmarks!(distributed_benchmark)
run_benchmarks!(distributed_benchmark) # Takes a couple of hours
save_results!(distributed_benchmark,"scaling_8.json")
# ================================ #
# 16 worker cores
addprocs(16)
include("dayahead_benchmark.jl")
distributed_benchmark = prepare_dbenchmark()
warmup_benchmarks!(distributed_benchmark)
run_benchmarks!(distributed_benchmark) # Takes a couple of hours
save_results!(distributed_benchmark,"scaling_16.json")
# ================================ #
# Collect all distributed benchmarks and plot
distributed_benchmark = load_results!("scaling_1.json")
distributed_benchmark.results["1"] = copy(distributed_benchmark.results["1000"])
delete!(distributed_benchmark.results, "1000")
scaling_2 = load_results!("scaling_2.json")
distributed_benchmark.results["2"] = copy(scaling_2.results["1000"])
scaling_4 = load_results!("scaling_4.json")
distributed_benchmark.results["4"] = copy(scaling_4.results["1000"])
scaling_8 = load_results!("scaling_8.json")
distributed_benchmark.results["8"] = copy(scaling_8.results["1000"])
scaling_16 = load_results!("scaling_16.json")
distributed_benchmark.results["16"] = copy(scaling_16.results["1000"])
save_results!(distributed_benchmark,"distributed_benchmarks.json")
using Plots
pyplot()
p = plot(TimePlot(distributed_benchmark)) # Creates Fig. 3
savefig(p,"strong_scaling_comptime.pdf")
p = plot(ScalingPlot(distributed_benchmark)) # Creates Fig. 4
savefig(p,"strong_scaling_efficiency.pdf")
# ================================ #

max_stars_repo_name: martinbiel/paw-atm-reproducibility

max_issues_repo_name: martinbiel/paw-atm-reproducibility

max_forks_repo_name: martinbiel/paw-atm-reproducibility


text:
module SimplePool

# simple scan into a list of free buffers

import ..@pool_timeit, ..actual_alloc, ..actual_free

using CUDAdrv

## tunables

# how much larger a buf can be to fullfil an allocation request.
# lower values improve efficiency, but increase pressure on the underlying GC.
function max_oversize(sz)
    if sz <= 2^20       # 1 MiB
        # small buffers are fine no matter
        return typemax(Int)
    elseif sz <= 2^20   # 32 MiB
        return 2^20
    else
        return 2^22
    end
end

## pooling

const available = Set{Mem.Buffer}()
const allocated = Set{Mem.Buffer}()

function scan(sz)
    for buf in available
        if sz <= sizeof(buf) <= max_oversize(sz)
            delete!(available, buf)
            return buf
        end
    end
    return
end

function reclaim(sz)
    freed = 0
    while freed < sz && !isempty(available)
        buf = pop!(available)
        actual_free(buf)
        freed += sizeof(buf)
    end

    return freed
end

function pool_alloc(sz)
    buf = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc(false)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc(true)" GC.gc(true)
        end

        @pool_timeit "$phase.1 scan" begin
            buf = scan(sz)
        end
        buf === nothing || break

        @pool_timeit "$phase.2 alloc" begin
            buf = actual_alloc(sz)
        end
        buf === nothing || break

        @pool_timeit "$phase.3 reclaim + alloc" begin
            reclaim(sz)
            buf = actual_alloc(sz)
        end
        buf === nothing || break
    end

    return buf
end

function pool_free(buf)
    push!(available, buf)
end

## interface

init() = return

function deinit()
    @assert isempty(allocated) "Cannot deinitialize memory pool with outstanding allocations"

    for buf in available
        actual_free(buf)
    end
    empty!(available)

    return
end

function alloc(sz)
    buf = pool_alloc(sz)
    if buf !== nothing
        push!(allocated, buf)
    end
    return buf
end

function free(buf)
    delete!(allocated, buf)
    pool_free(buf)
    return
end

used_memory() = isempty(allocated) ? 0 : sum(sizeof, allocated)

cached_memory() = isempty(available) ? 0 : sum(sizeof, available)

dump() = return

end


text:
using FlowWorkspace
using Test

######################################################### loading
load("./data/den.wsp"; files="./data/den/*.fcs")

######################################################### quadgates
load("./data/compensation.wsp"; files="./data/compensation/*fcs")


text:
using Gadfly, Compose

set_default_plot_size(6inch, 3inch)

p1 = plot(x=[1,2,3], y=[4,5,6], Geom.point);
p2 = plot(x=[1,2,3], y=[4,5,6], Geom.line);
title(vstack(p1,p2), "foo", font("helvetica"))


text:
export pol_scalar, pol_scalar_spiral, pol_x, pol_y, pol_circ, pol_circ_spiral
# define a bunch of standard pupil polarizations

"""
    pol_scalar(T, xypos)

ignores polarization aspects in the calculation but calculates only based on (high-NA) scalar theory.
This is a lot faster but not as accurate.

"""
function pol_scalar(T, xypos)  # returns only one "polarization" indicating that the whole calculation is to be performed scalar
    (one(T),)
end

"""
    pol_scalar_spiral(T, xypos)

ignores polarization aspects in the calculation but calculates only based on (high-NA) scalar theory.
This version still includes a (scalar) phase spiral.
This is a lot faster but not as accurate.
"""
function pol_scalar_spiral(T, xypos) # e.g. for STED microscopy
    cis.(atan(xypos...))
end

"""
    pol_x(T, xypos)

assumes x-polarization in illumination or an x-oriented polarizer in detection.
In a high-NA objective this is converted into XYZ electric fields at the focus.
"""
function pol_x(T, xypos)
    (one(T),zero(T))
end

"""
    pol_y(T, xypos)

assumes y-polarization in illumination or an x-oriented polarizer in detection.
In a high-NA objective this is converted into XYZ electric fields at the focus.
"""
function pol_y(T, xypos)
    (zero(T), one(T))
end

"""
    pol_circ(T, xypos)

assumes circular polarization in illumination or an x-oriented polarizer in detection.
In a high-NA objective this is converted into XYZ electric fields at the focus.
"""
function pol_circ(T, xypos)
    (one(T)/sqrt(2)+0im, one(T)/sqrt(2)*1im)
end

"""
    pol_circ_spiral(T, xypos)

assumes circular polarization in illumination or an x-oriented polarizer in detection.
This version includes phase spiral defining the local (xypos-dependent) phase of both x and y polarization.
"""
function pol_circ_spiral(T, xypos) # e.g. for STED microscopy
    pol_circ(T, xypos) .* cis.(atan(xypos...))
end


text:
# indexing

export global_size, synchronize_threads, linear_index

# thread indexing functions
for f in (:blockidx, :blockdim, :threadidx, :griddim)
    @eval $f(ctx::AbstractKernelContext)::Int = error("Not implemented") # COV_EXCL_LINE
    @eval export $f
end

"""
    global_size(ctx::AbstractKernelContext)

Global size == blockdim * griddim == total number of kernel execution
"""
@inline function global_size(ctx::AbstractKernelContext)
    griddim(ctx) * blockdim(ctx)
end

"""
    linear_index(ctx::AbstractKernelContext)

linear index corresponding to each kernel launch (in OpenCL equal to get_global_id).

"""
@inline function linear_index(ctx::AbstractKernelContext)
    (blockidx(ctx) - 1) * blockdim(ctx) + threadidx(ctx)
end

"""
    linearidx(A, ctxsym = :ctx)

Macro form of `linear_index`, which calls return when out of bounds.
So it can be used like this:

    ```julia
    function kernel(ctx::AbstractKernelContext, A)
        idx = @linearidx A ctx
        # from here on it's save to index into A with idx
        @inbounds begin
            A[idx] = ...
        end
    end
    ```
"""
macro linearidx(A, ctxsym = :ctx)
    quote
        x1 = $(esc(A))
        i1 = linear_index($(esc(ctxsym)))
        i1 > length(x1) && return
        i1
    end
end

"""
    cartesianidx(A, ctxsym = :ctx)

Like [`@linearidx(A, ctxsym = :ctx)`](@ref), but returns a N-dimensional `CartesianIndex`.
"""
macro cartesianidx(A, ctxsym = :ctx)
    quote
        x = $(esc(A))
        i = @linearidx(x, $(esc(ctxsym)))
        CartesianIndices(x)[i]
    end
end


text:
function Base.fft!{T<:BigFloat}(data::Vector{T})
	nn=int(length(data)/2)
	Tpi=convert(T,pi)
	n=nn << 1
	j=1
	for i=1:2:n-1
		if j>i
			data[j],data[i]=data[i],data[j]
			data[j+1],data[i+1]=data[i+1],data[j+1]
		end
		m=nn
		while m >= 2 && j > m
			j -= m
			m >>= 1
		end
		j += m
	end
	mmax=2
	while n > mmax
		istep=mmax << 1
		theta=-2Tpi/mmax
		wtemp=sin(theta/2)
		wpr = -2wtemp*wtemp
		wpi=sin(theta)
		wr=one(T)
		wi=zero(T)
		for m=1:2:mmax-1
			for i=m:istep:n
				j=i+mmax
				tempr=wr*data[j]-wi*data[j+1]
				tempi=wr*data[j+1]+wi*data[j]
				data[j]=data[i]-tempr
				data[j+1]=data[i+1]-tempi
				data[i] += tempr
				data[i+1] += tempi
			end
			wr=(wtemp=wr)*wpr-wi*wpi+wr
			wi=wi*wpr+wtemp*wpi+wi
		end
        mmax=istep
	end
	return data
end

function Base.fft{T<:BigFloat}(x::Vector{Complex{T}})
	n=length(x)
	y = Array(T,2n)
	for i = 1:n
		y[2i-1] = x[i].re
		y[2i] = x[i].im
	end
	Base.fft!(y)
	return complex(y[1:2:2n-1],y[2:2:2n])
end
Base.fft{T<:BigFloat}(x::Vector{T}) = Base.fft(complex(x))

function Base.ifft{T<:BigFloat}(x::Vector{Complex{T}})
	n=length(x)
	y = Array(T,2n)
	for i = 1:n
		y[2i-1] = x[i].re
		y[2i] = -x[i].im
	end
	Base.fft!(y)
	x = complex(y[1:2:2n-1],-y[2:2:2n])/n
end


text:
#=const CSVC = Int32(0)
const NuSVC = Int32(1)
const OneClassSVM = Int32(2)
const EpsilonSVR = Int32(3)
const NuSVR = Int32(4)

const Linear = Int32(0)
const Polynomial = Int32(1)
const RBF = Int32(2)
const Sigmoid = Int32(3)
const Precomputed = Int32(4)=#

verbosity = false

immutable SVMNode
    index::Int32
    value::Float64
end

immutable SVMProblem
    l::Int32
    y::Ptr{Float64}
    x::Ptr{Ptr{SVMNode}}
end

immutable SVMParameter
    svm_type::Int32
    kernel_type::Int32
    degree::Int32
    gamma::Float64
    coef0::Float64

    cache_size::Float64
    eps::Float64
    C::Float64
    nr_weight::Int32
    weight_label::Ptr{Int32}
    weight::Ptr{Float64}
    nu::Float64
    p::Float64
    shrinking::Int32
    probability::Int32
end

# immutable SVMModel
#   param::SVMParameter
#   nr_class::Int32
#   l::Int32
#   SV::Ptr{Ptr{SVMNode}}
#   sv_coef::Ptr{Ptr{Float64}}
#   rho::Ptr{Float64}
#   probA::Ptr{Float64}
#   probB::Ptr{Float64}
#   sv_indices::Ptr{Int32}

#   label::Ptr{Int32}
#   nSV::Ptr{Int32}

#   free_sv::Int32
# end

type SVMModel{T}
    ptr::Ptr{Void}
    param::Vector{SVMParameter}
    # Prevent these from being garbage collected
    problem::Vector{SVMProblem}
    nodes::Array{SVMNode}
    nodeptr::Vector{Ptr{SVMNode}}

    labels::Vector{T}
    weight_labels::Vector{Int32}
    weights::Vector{Float64}
    nfeatures::Int
    verbose::Bool
end

let libsvm = C_NULL
    global get_libsvm
    function get_libsvm()
        if libsvm == C_NULL
          @windows_only begin
              libsvm = Libdl.dlopen(joinpath(Pkg.dir(), "Gabi", "deps","windows","libsvm-3.21","libsvm.dll"))
          end
          @unix_only begin
              libsvm = Libdl.dlopen(joinpath(Pkg.dir(), "Gabi", "deps","unix","libsvm-3.21","libsvm.so.2"))
          end
            # libsvm = Libdl.dlopen("/usr/local/Cellar/libsvm/3.21/lib/libsvm.2.dylib")
            ccall(Libdl.dlsym(libsvm, :svm_set_print_string_function), Void,
                (Ptr{Void},), cfunction(svmprint, Void, (Ptr{UInt8},)))
        end
        libsvm
    end
end

macro cachedsym(symname)
    cached = gensym()
    quote
        let $cached = C_NULL
            global ($symname)
            ($symname)() = ($cached) == C_NULL ?
                ($cached = Libdl.dlsym(get_libsvm(), $(string(symname)))) : $cached
        end
    end
end
@cachedsym svm_train
@cachedsym svm_predict_values
@cachedsym svm_predict_probability
@cachedsym svm_free_model_content

function grp2idx{T, S <: Real}(::Type{S}, labels::AbstractVector,
    label_dict::Dict{T, Int32}, reverse_labels::Vector{T})

    idx = Array(S, length(labels))
    nextkey = length(reverse_labels) + 1
    for i = 1:length(labels)
        key = labels[i]
        if (idx[i] = get(label_dict, key, nextkey)) == nextkey
            label_dict[key] = nextkey
            push!(reverse_labels, key)
            nextkey += 1
        end
    end
    idx
end

function instances2nodes{U<:Real}(instances::AbstractMatrix{U})
    nfeatures = size(instances, 1)
    ninstances = size(instances, 2)
    nodeptrs = Array(Ptr{SVMNode}, ninstances)
    nodes = Array(SVMNode, nfeatures + 1, ninstances)

    for i=1:ninstances
        k = 1
        for j=1:nfeatures
            nodes[k, i] = SVMNode(Int32(j), Float64(instances[j, i]))
            k += 1
        end
        nodes[k, i] = SVMNode(Int32(-1), NaN)
        nodeptrs[i] = pointer(nodes, (i-1)*(nfeatures+1)+1)
    end

    (nodes, nodeptrs)
end

function instances2nodes{U<:Real}(instances::SparseMatrixCSC{U})
    ninstances = size(instances, 2)
    nodeptrs = Array(Ptr{SVMNode}, ninstances)
    nodes = Array(SVMNode, nnz(instances)+ninstances)

    j = 1
    k = 1
    for i=1:ninstances
        nodeptrs[i] = pointer(nodes, k)
        while j < instances.colptr[i+1]
            val = instances.nzval[j]
            nodes[k] = SVMNode(Int32(instances.rowval[j]), Float64(val))
            k += 1
            j += 1
        end
        nodes[k] = SVMNode(Int32(-1), NaN)
        k += 1
    end

    (nodes, nodeptrs)
end

function svmprint(str::Ptr{UInt8})
    if verbosity::Bool
        print(bytestring(str))
    end
    nothing
end

function indices_and_weights{T, U<:Real}(labels::AbstractVector{T},
        instances::AbstractMatrix{U},
        weights::Union{Dict{T, Float64}, Void}=nothing)
    label_dict = Dict{T, Int32}()
    reverse_labels = Array(T, 0)
    idx = grp2idx(Float64, labels, label_dict, reverse_labels)

    if length(labels) != size(instances, 2)
        error("""Size of second dimension of training instance matrix
        ($(size(instances, 2))) does not match length of labels
        ($(length(labels)))""")
    end

    # Construct SVMParameter
    if weights == nothing || length(weights) == 0
        weight_labels = Int32[]
        weights = Float64[]
    else
        weight_labels = grp2idx(Int32, keys(weights), label_dict,
            reverse_labels)
        weights = Float64(values(weights))
    end

    (idx, reverse_labels, weights, weight_labels)
end

function svmtrain{U<:Real,T}(labels::Vector{T},
        instances::Matrix{U}; svm_type::Int32=CSVC,
        kernel_type::Int32=RBF, degree::Integer=3,
        gamma::Float64=1.0/size(instances, 2), coef0::Float64=0.0,
        C::Float64=1.0, nu::Float64=0.5, p::Float64=0.1,
        cache_size::Float64=100.0, eps::Float64=0.001, shrinking::Bool=true,
        probability_estimates::Bool=false,
        weights::Union{Dict{T, Float64}, Void}=nothing,
        verbose::Bool=false)
    global verbosity

    (idx, reverse_labels, weights, weight_labels) = indices_and_weights(labels,
        instances, weights)

    param = Array(SVMParameter, 1)
    param[1] = SVMParameter(svm_type, kernel_type, Int32(degree), Float64(gamma),
        coef0, cache_size, eps, C, Int32(length(weights)),
        pointer(weight_labels), pointer(weights), nu, p, Int32(shrinking),
        Int32(probability_estimates))

    # Construct SVMProblem
    (nodes, nodeptrs) = instances2nodes(instances)
    problem = SVMProblem[SVMProblem(Int32(size(instances, 2)), pointer(idx),
        pointer(nodeptrs))]

    verbosity = verbose
    ptr = ccall(svm_train(), Ptr{Void}, (Ptr{SVMProblem},
        Ptr{SVMParameter}), problem, param)

    model = SVMModel(ptr, param, problem, nodes, nodeptrs, reverse_labels,
        weight_labels, weights, size(instances, 1), verbose)
    finalizer(model, svmfree)
    model
end

function svmcv{T, U<:Real, V<:Real, X<:Real}(labels::AbstractVector{T},
        instances::AbstractMatrix{U}, nfolds::Int=5,
        C::Union{V, AbstractArray{V}}=2.0.^(-5:2:15),
        gamma::Union{X, AbstractArray{X}}=2.0.^(3:-2:-15);
        svm_type::Int32=CSVC, kernel_type::Int32=RBF, degree::Integer=3,
        coef0::Float64=0.0, nu::Float64=0.5, p::Float64=0.1,
        cache_size::Float64=100.0, eps::Float64=0.001, shrinking::Bool=true,
        weights::Union{Dict{T, Float64}, Void}=nothing,
        verbose::Bool=false)
    global verbosity
    verbosity = verbose
    (idx, reverse_labels, weights, weight_labels) = indices_and_weights(labels,
        instances, weights)

    # Construct SVMParameters
    params = Array(SVMParameter, length(C), length(gamma))
    degree = Int32(degree)
    nweights = Int32(length(weights))
    shrinking = Int32(shrinking)
    for i = 1:length(C), j = 1:length(gamma)
        params[i, j] = SVMParameter(svm_type, kernel_type, Int32(degree),
            Float64(gamma[j]), coef0, cache_size, eps, Float64(C[i]), nweights,
            pointer(weight_labels), pointer(weights), nu, p, shrinking,
            Int32(0))
    end

    # Get information about classes
    (nodes, nodeptrs) = instances2nodes(instances)
    n_classes = length(reverse_labels)
    n_class = zeros(Int, n_classes)
    for id in idx
        n_class[id] += 1
    end
    by_class = [Array(Int, n) for n in n_class]
    idx_class = zeros(Int, n_classes)
    for i = 1:length(idx)
        cl = idx[i]
        by_class[cl][(idx_class[cl] += 1)] = i
    end
    for i = 1:n_classes
        shuffle!(by_class[i])
    end
    prop_class = n_class/length(idx)

    # Perform cross-validation
    decvalues = Array(Float64, length(reverse_labels))
    fold_classes = Array(Range1{Int}, n_classes)
    perf = zeros(Float64, length(C), length(gamma))
    for i = 1:nfolds
        # Get range for test for each class
        fold_ntest = 0
        for j = 1:n_classes
            fold_classes[j] =
                div((i-1)*n_class[j], nfolds)+1:div(i*n_class[j], nfolds)
            fold_ntest += length(fold_classes[j])
        end

        # Get indices of test and training instances
        fold_ntrain = length(idx) - fold_ntest
        fold_train = Array(Int, fold_ntrain)
        fold_test = Array(Int, fold_ntest)
        itrain = 0
        itest = 0
        for j = 1:n_classes
            train_range = fold_classes[j]
            class_idx = by_class[j]

            n = first(train_range) - 1
            if n > 0
                fold_train[itest+1:itest+n] = class_idx[1:n]
                itest += n
            end
            n = n_class[j] - last(train_range)
            if n > 0
                fold_train[itest+1:itest+n] = class_idx[last(train_range)+1:end]
                itest += n
            end

            n = length(train_range)
            fold_test[itrain+1:itrain+n] = class_idx[train_range]
            itrain += n
        end

        fold_train_nodeptrs = nodeptrs[fold_train]
        fold_train_labels = idx[fold_train]
        problem = SVMProblem[SVMProblem(Int32(fold_ntrain),
            pointer(fold_train_labels),
            pointer(fold_train_nodeptrs))]

        for j = 1:length(params)
            ptr = ccall(svm_train(), Ptr{Void}, (Ptr{SVMProblem},
                Ptr{SVMParameter}), problem, pointer(params, j))
            correct = 0
            for k in fold_test
                correct += ccall(svm_predict_values(), Float64, (Ptr{Void},
                    Ptr{SVMNode}, Ptr{Float64}), ptr, nodeptrs[k], decvalues) ==
                    idx[k]
            end
            ccall(svm_free_model_content(), Void, (Ptr{Void},), ptr)
            perf[j] += correct/fold_ntest
        end
    end

    best = ind2sub(size(perf), indmax(perf))
    (C[best[1]], gamma[best[2]], perf/nfolds)
end

svmfree(model::SVMModel) = ccall(svm_free_model_content(), Void, (Ptr{Void},),
    model.ptr)

function svmpredict{U<:Real,T}(model::SVMModel{T},
        instances::Matrix{U})
    global verbosity
    ninstances = size(instances, 2)

    if size(instances, 1) != model.nfeatures
        error("Model has $(model.nfeatures) but $(size(instances, 1)) provided")
    end

    (nodes, nodeptrs) = instances2nodes(instances)
    class = Array(T, ninstances)
    nlabels = nothing
    decvalues = nothing
    verbosity = model.verbose
    svm_tp::Int =  model.param[1].svm_type
    fn = model.param[1].probability == 1 ? svm_predict_probability() :
        svm_predict_values()
     if (svm_tp <= 2)
        #Classification
        nlabels = length(model.labels)
        decvalues = Array(Float64, nlabels, ninstances)

        for i = 1:ninstances
            output = ccall(fn, Float64, (Ptr{Void}, Ptr{SVMNode}, Ptr{Float64}),
                model.ptr, nodeptrs[i], pointer(decvalues, nlabels*(i-1)+1))
            class[i] = model.labels[round(Int,output)]
        end
    else
      #Regression
      nlabels = length(model.labels)
      decvalues = Array(Float64, nlabels, ninstances)
      for i = 1:ninstances
          output = ccall(fn, Float64, (Ptr{Void}, Ptr{SVMNode}, Ptr{Float64}),
              model.ptr, nodeptrs[i],pointer(decvalues, nlabels*(i-1)+1))
          class[i] = output
      end
    end
    (class, decvalues)
end


text:
#=

  Curious numbers in ConstraintSolver.jl

  """
  Curious Numbers from "Amusements in Mathematics, Dudeney", number 114.

  The number 48 has this peculiarity, that if you add 1 to it the result
  is a square number, and if you add 1 to its half, you also get a
  square number. Now, there is no limit to the numbers that have this
  peculiarity, and it is an interesting puzzle to find three more of
  them---the smallest possible numbers. What are they?
  """

  The least such numbers are:
  [
   [48,49,7,24,25,5],
   [1680,1681,41,840,841,29],
   [57120,57121,239,28560,28561,169],
   [1940448,1940449,1393,970224,970225,985]
  ]

  Note: This model use the experimental (and very slow)
  non-linear mult constraint (mult(model,A,B,C) -> C = A*B)
  It thus only handle the first two results shows above.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function cur_num(max_val=100,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 6
    @variable(model, 1 <= x[1:n] <= max_val, Int)
    X,A,B,C,D,E = x

    @constraint(model,X + 1 == A) # if you add 1 to it
    # @constraint(model,A == B * B) # the result is a square number
    mult(model,B,B,A)

    # @cons)raint(model,X == 2 * C) # if you to its half
    mult(model,2,C,X)
    @constraint(model,C + 1 == D) # add 1
    # @constraint(model,D == E * E) # you also get a square number
    mult(model,E,E,D)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

# for max_val in [100,2000]
for max_val in [100,2000]
    println("\nmax_val:$max_val")
    @time cur_num(max_val)
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function recursive_backtracker!(grid, start_pos=(1, 1))
    start_at = grid[start_pos...]
    stack = [start_at]
    while !isempty(stack)
        current = stack[end]
        nbs = filter(x -> isempty(x.connections), neighbors(grid, current))
        if isempty(nbs)
            pop!(stack)
        else
            neighbor = rand(nbs)
            connect!(current, neighbor)
            push!(stack, neighbor)
        end
    end
    return grid
end


text:
using Compat
import Compat.String

@testset "Types" begin
    solver = NelderMead()
    T = typeof(solver)
    trace = OptimizationTrace{Float64, T}()
    push!(trace,OptimizationState{Float64, T}(1,1.0,1.0,Dict()))
    push!(trace,OptimizationState{Float64, T}(2,1.0,1.0,Dict()))
    @test length(trace) == 2
    @test trace[end].iteration == 2

    prob = MultivariateProblems.UnconstrainedProblems.examples["Rosenbrock"]
    f_prob = MVP.objective(prob)
    for g_free in (NelderMead(), SimulatedAnnealing())
        res = Optim.optimize(f_prob, prob.initial_x, g_free)
        @test typeof(f_prob(prob.initial_x)) == typeof(Optim.minimum(res))
        @test eltype(prob.initial_x) == eltype(Optim.minimizer(res))

        io = IOBuffer()
        show(io, res)
        s = String(take!(io))

        lines = split(s, '\n')
        @test lines[1] == "Results of Optimization Algorithm"
        @test startswith(lines[2], " * Algorithm: ")
        @test startswith(lines[3], " * Starting Point: ")
        @test startswith(lines[4], " * Minimizer: [")
        @test startswith(lines[5], " * Minimum: ")
        @test startswith(lines[6], " * Iterations: ")
        @test startswith(lines[7], " * Convergence: ")
        if res.method == "Nelder-Mead"
            @test startswith(lines[8], "   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: ")
            @test startswith(lines[9], "   * Reached Maximum Number of Iterations: ")
            @test startswith(lines[10], " * Objective Calls: ")
        elseif res.method == "Simulated Annealing"
            @test startswith(lines[8], "   * |x - x'| ≤ ")
            @test startswith(lines[9], "   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|")
            @test startswith(lines[10], "   * |g(x)| ≤ ")
            @test startswith(lines[11], "   * Stopped by an increasing objective:")
            @test startswith(lines[12], "   * Reached Maximum Number of Iterations: ")
            @test startswith(lines[13], " * Objective Calls: ")
        end
    end

    res = Optim.optimize(MVP.objective(prob), MVP.gradient(prob), prob.initial_x, LBFGS())
    @test typeof(f_prob(prob.initial_x)) == typeof(Optim.minimum(res))
    @test eltype(prob.initial_x) == eltype(Optim.minimizer(res))

    io = IOBuffer()
    show(io, res)
    s = String(take!(io))

    lines = split(s, '\n')
    @test lines[1] == "Results of Optimization Algorithm"
    @test startswith(lines[2], " * Algorithm: ")
    @test startswith(lines[3], " * Starting Point: ")
    @test startswith(lines[4], " * Minimizer: [")
    @test startswith(lines[5], " * Minimum: ")
    @test startswith(lines[6], " * Iterations: ")
    @test startswith(lines[7], " * Convergence: ")
    @test startswith(lines[8], "   * |x - x'| ≤ ")
    @test startswith(lines[9],  "     |x - x'| = ")
    @test startswith(lines[10], "   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|:") # TODO: regexp
    @test startswith(lines[11], "     |f(x) - f(x')| = 7.43e+10 |f(x)|") # TODO: regexp
    @test startswith(lines[12], "   * |g(x)| ≤ ")
    @test startswith(lines[13], "     |g(x)| = ")
    @test startswith(lines[14], "   * Stopped by an increasing objective:")
    @test startswith(lines[15], "   * Reached Maximum Number of Iterations: ")
    @test startswith(lines[16], " * Objective Calls: ")
    @test startswith(lines[17], " * Gradient Calls: ")
    if res.method in ("Newton's Method", "Newton's Method (Trust Region)")
        @test startswith(lines[18], " * Hessian Calls: ")
    end
end


text:
# ---------------------------------------------------------------------------
# density

@recipe function f(::Type{Val{:density}}, x, y, z; trim=false)
    newx, newy = violin_coords(y, trim=trim)
    if Plots.isvertical(plotattributes)
        newx, newy = newy, newx
    end
    x := newx
    y := newy
    seriestype := :path
    ()
end
Plots.@deps density path

# ---------------------------------------------------------------------------
# cumulative density

@recipe function f(::Type{Val{:cdensity}}, x, y, z; trim=false,
                   npoints = 200)
    newx, newy = violin_coords(y, trim=trim)

    if Plots.isvertical(plotattributes)
        newx, newy = newy, newx
    end

    newy = [sum(newy[1:i]) for i = 1:length(newy)] / sum(newy)

    x := newx
    y := newy
    seriestype := :path
    ()
end
Plots.@deps cdensity path

ea_binnumber(y, bin::AbstractVector) = error("You cannot specify edge locations for equal area histogram")
ea_binnumber(y, bin::Real) = (floor(bin) == bin || error("Only integer or symbol values accepted by bins"); Int(bin))
ea_binnumber(y, bin::Int) = bin
ea_binnumber(y, bin::Symbol) = Plots._auto_binning_nbins((y,), 1, mode = bin)

@recipe function f(::Type{Val{:ea_histogram}}, x, y, z)
    bin = ea_binnumber(y, plotattributes[:bins])
    bins := quantile(y, range(0, stop = 1, length = bin+1))
    normalize := :density
    seriestype := :barhist
    ()
end
Plots.@deps histogram barhist

push!(Plots._histogram_like, :ea_histogram)

@shorthands ea_histogram

@recipe function f(::Type{Val{:testhist}}, x, y, z)
    markercolor --> :red
    seriestype := :scatter
    ()
end
@shorthands testhist
# ---------------------------------------------------------------------------
# Compute binsizes using Wand (1997)'s criterion
# Ported from R code located here https://github.com/cran/KernSmooth/tree/master/R

"Returns optimal histogram edge positions in accordance to Wand (1995)'s criterion'"
Plots.wand_edges(x::AbstractVector, args...) = (binwidth = wand_bins(x, args...); minimum(x)-binwidth:binwidth:maximum(x)+binwidth)

"Returns optimal histogram bin widths in accordance to Wand (1995)'s criterion'"
function wand_bins(x, scalest = :minim, gridsize = 401, range_x = extrema(x), trun = true)

    n = length(x)
    minx, maxx = range_x
    gpoints = range(minx, stop = maxx, length = gridsize)
    gcounts = linbin(x, gpoints, trun = trun)

    scalest = if scalest == :stdev
        sqrt(var(x))
    elseif scalest == :iqr
        (quantile(x, 3//4) - quantile(x, 1//4))/1.349
    elseif scalest == :minim
        min((quantile(x, 3//4) - quantile(x, 1//4))/1.349, sqrt(var(x)))
    else
        error("scalest must be one of :stdev, :iqr or :minim (default)")
    end

    scalest == 0 && error("scale estimate is zero for input data")
    sx = (x .- mean(x))./scalest
    sa = (minx - mean(x))/scalest
    sb = (maxx - mean(x))/scalest

    gpoints = range(sa, stop = sb, length = gridsize)
    gcounts = linbin(sx, gpoints, trun = trun)

    hpi = begin
        alpha = ((2/(11 * n))^(1/13)) * sqrt(2)
        psi10hat = bkfe(gcounts, 10, alpha, [sa, sb])
        alpha = (-105 * sqrt(2/pi)/(psi10hat * n))^(1//11)
        psi8hat = bkfe(gcounts, 8, alpha, [sa, sb])
        alpha = (15 * sqrt(2/pi)/(psi8hat * n))^(1/9)
        psi6hat = bkfe(gcounts, 6, alpha, [sa, sb])
        alpha = (-3 * sqrt(2/pi)/(psi6hat * n))^(1/7)
        psi4hat = bkfe(gcounts, 4, alpha, [sa, sb])
        alpha = (sqrt(2/pi)/(psi4hat * n))^(1/5)
        psi2hat = bkfe(gcounts, 2, alpha, [sa, sb])
        (6/(-psi2hat * n))^(1/3)
    end

    scalest * hpi
end

function linbin(X, gpoints; trun = true)
    n, M = length(X), length(gpoints)

    a, b = gpoints[1], gpoints[M]
    gcnts = zeros(M)
    delta = (b-a)/(M-1)

    for i in 1:n
        lxi = ((X[i]-a)/delta) + 1
        li = floor(Int, lxi)
        rem = lxi - li

        if 1 <= li < M
            gcnts[li] += 1-rem
            gcnts[li+1] += rem
        end

        if !trun
            if lt < 1
                gcnts[1] += 1
            end

            if li >= M
                gcnts[M] += 1
            end
        end
    end
    gcnts
end

"binned kernel function estimator"
function bkfe(gcounts, drv, bandwidth, range_x)
    bandwidth <= 0 && error("'bandwidth' must be strictly positive")

    a, b = range_x
    h = bandwidth
    M = length(gcounts)
    gpoints = range(a, stop = b, length = M)

    ## Set the sample size and bin width

    n = sum(gcounts)
    delta = (b-a)/(M-1)

    ## Obtain kernel weights

    tau = 4 + drv
    L = min(floor(tau*h/delta), M)

    lvec = 0:L
    arg = lvec .* delta/h

    kappam = pdf(Normal(),arg)/(h^(drv+1))
    hmold0, hmnew = ones(length(arg)), ones(length(arg))
    hmold1 = arg

    if drv >= 2
        for i in (2:drv)
            hmnew = arg .* hmold1 .- (i-1) .* hmold0
            hmold0 = hmold1       # Compute mth degree Hermite polynomial
            hmold1 = hmnew        # by recurrence.
        end
    end
    kappam = hmnew .* kappam

    ## Now combine weights and counts to obtain estimate
    ## we need P >= 2L+1L, M: L <= M.
    P = 2^(ceil(log(M+L+1)/log(2)))
    kappam = [kappam; zeros(P-2*L-1); reverse(kappam[2:end])]
    Gcounts = [gcounts; zeros(P-M)]
    kappam = fft(kappam)
    Gcounts = fft(Gcounts)

    sum(gcounts .* (real(ifft(kappam .* Gcounts)))[1:M] )/(n^2)
end


text:
g_u(x,t) = VectorValue(0.0,0.0,0.0)
g_u(t::Real) = x -> g_u(x,t)
∂tu(t) = x -> VectorValue(0.0,0.0,0.0)
∂tu(x,t) = ∂tu(t)(x)
∂t(::typeof(g_u)) = ∂tu

g_j(x,t) = VectorValue(0.0,0.0,0.0)
g_j(t::Real) = x -> g_j(x,t)
∂tj(t) = x -> VectorValue(0.0,0.0,0.0)
∂tj(x,t) = ∂tj(t)(x)
∂t(::typeof(g_j)) = ∂tj

function transient_duct_flow(;nx::Int=3, ny::Int=3, Re::Float64 = 10.0,
    Ha::Float64 = 10.0, Δt::Float64 = 0.1, t0::Float64 = 0.0, tF::Float64 = 1.0,
    θ::Float64 = 1.0, f_u = nothing, resultsfile = nothing)

  L = 1.0
  N = Ha^2/Re
  K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
  ∂p∂z = -Re * K / L^3

  B = VectorValue(0.0,1.0,0.0)

  if f_u === nothing
    _f_u = VectorValue(0.0,0.0, -∂p∂z) * L/Re^2 # Assumed ν=L=1 by default
  else
    _f_u = f_u
  end

  # Discretization
  order = 2
  domain = (-1.0,1.0,-1.0,1.0,0.0,0.1)
  map(x) = VectorValue(sign(x[1])*(abs(x[1])*0.5)^0.5,
                       sign(x[2])*(abs(x[2])*0.5)^0.5, x[3]*sqrt(2)/2)*2/sqrt(2)

  dirichlet_tags_u = append!(collect(1:20),[23,24,25,26])
  dirichlet_tags_j = append!(collect(1:20),[23,24,25,26])

  partition=(nx,ny,3)
  model = CartesianDiscreteModel(domain,partition;
    isperiodic=(false,false,true),map=map)

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet_u",dirichlet_tags_u)
  add_tag_from_tags!(labels,"dirichlet_j",dirichlet_tags_j)

  Vu = FESpace(model, ReferenceFE(lagrangian,VectorValue{3,Float64},order);
      conformity=:H1, dirichlet_tags="dirichlet_u")

  Vp = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:P);
      conformity=:L2, constraint=:zeromean)

  Vj = FESpace(model, ReferenceFE(raviart_thomas,Float64,order-1);
      conformity=:Hdiv, dirichlet_tags="dirichlet_j")

  Vφ = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:Q);
      conformity=:L2, constraint=:zeromean)

  U = TransientTrialFESpace(Vu,g_u)
  P = TransientTrialFESpace(Vp)
  J = TransientTrialFESpace(Vj,g_j)
  Φ = TransientTrialFESpace(Vφ)

  Y = MultiFieldFESpace([Vu, Vp, Vj, Vφ])
  X = TransientMultiFieldFESpace([U, P, J, Φ])

  trian = Triangulation(model)
  degree = 2*(order)
  dΩ = Measure(trian,degree)

  res(t,x,xt,y) = ∫(xt[1]⋅y[1])*dΩ +
    InductionlessMHD.dimensionless_residual(x, y, Re, N, B, _f_u,dΩ)
  jac(t,x,xt,dx,y) = InductionlessMHD.dimensionless_jacobian(x, dx, y, Re, N, B,dΩ)
  jac_t(t,x,xt,dxt,y) = ∫(dxt[1]⋅y[1])*dΩ

  op  = TransientFEOperator(res,jac,jac_t,X,Y)

  uh0 = interpolate(VectorValue(0.0,0.0,0.0),U(0.0))
  ph0 = interpolate(0.0,P(0.0))
  jh0 = interpolate(VectorValue(0.0,0.0,0.0),J(0.0))
  φh0 = interpolate(0.0,Φ(0.0))
  xh0 = interpolate([uh0,ph0,jh0,φh0],X(0.0))

  ls = LUSolver()
  nls = NLSolver(ls;
    show_trace=true, method=:newton, linesearch=BackTracking(), iterations=2)
  odes = ThetaMethod(nls,Δt,θ)
  solver = TransientFESolver(odes)

  xh_t = solve(solver,op,xh0,t0,tF)

  result = Base.iterate(xh_t)

  if resultsfile != nothing
    startPVD(resultsfile)
    write_timestep(resultsfile,t0,trian,
      cellfields=["uh"=>uh0, "ph"=>ph0, "jh"=>jh0, "φh"=>φh0])
    for (xh, t) in xh_t
      uh, ph, jh, φh = xh
      write_timestep(resultsfile,t,trian,
        cellfields=["uh"=>uh, "ph"=>ph, "jh"=>jh, "φh"=>φh])
    end
    closePVD(resultsfile)
  end

  (xh_t, trian, dΩ)
end


text:
#####
##### Weighted Essentially Non-Oscillatory (WENO) fifth-order advection scheme
#####

using OffsetArrays
using Oceananigans.Grids: with_halo
using Oceananigans.Architectures: arch_array
using Adapt
import Base: show

const two_32 = Int32(2)

const C3₀ = 3/10
const C3₁ = 3/5
const C3₂ = 1/10

const ƞ = Int32(2) # WENO exponent
const ε = 1e-6

"""
    struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, W} <: AbstractUpwindBiasedAdvectionScheme{2}

Weighted Essentially Non-Oscillatory (WENO) fifth-order advection scheme.

$(TYPEDFIELDS)
"""
struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, W} <: AbstractUpwindBiasedAdvectionScheme{2}

    "coefficient for ENO reconstruction on x-faces"
    coeff_xᶠᵃᵃ::XT
    "coefficient for ENO reconstruction on x-centers"
    coeff_xᶜᵃᵃ::XT
    "coefficient for ENO reconstruction on y-faces"
    coeff_yᵃᶠᵃ::YT
    "coefficient for ENO reconstruction on y-centers"
    coeff_yᵃᶜᵃ::YT
    "coefficient for ENO reconstruction on z-faces"
    coeff_zᵃᵃᶠ::ZT
    "coefficient for ENO reconstruction on z-centers"
    coeff_zᵃᵃᶜ::ZT

    "coefficient for WENO smoothness indicators on x-faces"
    smooth_xᶠᵃᵃ::XS
    "coefficient for WENO smoothness indicators on x-centers"
    smooth_xᶜᵃᵃ::XS
    "coefficient for WENO smoothness indicators on y-faces"
    smooth_yᵃᶠᵃ::YS
    "coefficient for WENO smoothness indicators on y-centers"
    smooth_yᵃᶜᵃ::YS
    "coefficient for WENO smoothness indicators on z-faces"
    smooth_zᵃᵃᶠ::ZS
    "coefficient for WENO smoothness indicators on z-centers"
    smooth_zᵃᵃᶜ::ZS
end

"""
    WENO5([FT = Float64;] grid = nothing, stretched_smoothness = false, zweno = false)

Construct a fifth-order weigthed essentially non-oscillatory advection scheme. The constructor allows
construction of WENO schemes on either uniform or stretched grids.

Keyword arguments
=================

  - `grid`: (defaults to `nothing`)
  - `stretched_smoothness`: When `true` it results in computing the coefficients for the smoothness
    indicators β₀, β₁ and β₂ so that they account for the stretched `grid`. (defaults to `false`)
  - `zweno`: When `true` implement a Z-WENO formulation for the WENO weights calculation. (defaults to
    `false`)

!!! warn "No support for WENO5 on curvilinear grids"
    Currently, WENO 5th-order advection schemes don't work for for curvilinear grids.
    Providing `WENO5(::AbstractCurvilinearGrid)` defaults to uniform setting, i.e.
    `WENO5(::AbstractCurvilinearGrid) = WENO5()`.

Not providing any keyword argument, `WENO5()` defaults to the uniform 5th-order coefficients ("uniform
setting) in all directions, using a JS-WENO formulation.

```jldoctest; filter = [r".*┌ Warning.*", r".*└ @ Oceananigans.*"]
julia> using Oceananigans

julia> WENO5()
┌ Warning: defaulting to uniform WENO scheme with Float64 precision, use WENO5(grid = grid) if this was not intended
└ @ Oceananigans.Advection .../src/Advection/weno_fifth_order.jl:90
WENO5 advection scheme with:
    ├── X regular
    ├── Y regular
    └── Z regular
```

`WENO5(grid = grid)` defaults to uniform interpolation coefficient for each of the grid directions that
is uniform (`typeof(Δc) <: Number`) while it precomputes the ENO coefficients for reconstruction for all
grid directions that are stretched. (After testing "on-the-fly" calculation of coefficients for stretched
directions ended up being way too expensive and, therefore, is not supported.)

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size = (3, 4, 5), x = (0, 1), y = (0, 1), z = [-10, -9, -7, -4, -1.5, 0]);

julia> WENO5(grid = grid)
WENO5 advection scheme with:
    ├── X regular
    ├── Y regular
    └── Z stretched
```

`WENO5(grid = grid, stretched_smoothness = true)` behaves similarly to `WENO5(grid = grid)` but,
additionally, it also computes the smoothness indicators coefficients, ``β₀``, ``β₁``, and ``β₂``,
taking into account the stretched dimensions.

`WENO5(zweno = true)` implements a Z-WENO formulation for the WENO weights calculation

Comments
========

All methods have the roughly the same execution speed except for `stretched_smoothness = true` that
requires more memory and is less computationally efficient, especially on GPUs. In addition, it has
not been found to be much impactful on the tested cases. As such, most of the times we urge users
to use `WENO5(grid = grid)`, as this increases accuracy on a stretched mesh  but does decreases
memory utilization (and also results in a slight speed-up).

(The above claims were made after some preliminary tests. Thus, we still users to perform some
benchmarks/checks before performing, e.g., a large simulation on a "weirdly" stretched grid.)

On the other hand, a Z-WENO formulation is *most of the times* beneficial (also in case of a uniform
mesh) with roughly the same performances (just a slight slowdown). The same can be said for the
stretched `WENO5(grid = grid)` formulation in case of stretched grids.

References
==========

Shu, Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic
    Conservation Laws, 1997, NASA/CR-97-206253, ICASE Report No. 97-65

Castro et al, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation
    laws, 2011, Journal of Computational Physics, 230(5), 1766-1792
"""
function WENO5(FT = Float64; grid = nothing, stretched_smoothness = false, zweno = false)

    metrics   = (:xᶠᵃᵃ, :xᶜᵃᵃ, :yᵃᶠᵃ, :yᵃᶜᵃ, :zᵃᵃᶠ, :zᵃᵃᶜ)
    dirsize   = (:Nx, :Nx, :Ny, :Ny, :Nz, :Nz)

    if grid isa Nothing
        @warn "defaulting to uniform WENO scheme with $(FT) precision, use WENO5(grid = grid) if this was not intended"
        for metric in metrics
            @eval $(Symbol(:coeff_ , metric)) = nothing
            @eval $(Symbol(:smooth_, metric)) = nothing
        end
    elseif !(grid isa RectilinearGrid)
        FT = Float32
        @warn "Stretched WENO is not supported with grids other than Rectilinear, defaulting to Uniform WENO"
        for metric in metrics
            @eval $(Symbol(:coeff_ , metric)) = nothing
            @eval $(Symbol(:smooth_, metric)) = nothing
        end
    else
        FT       = eltype(grid)
        arch     = grid.architecture
        new_grid = with_halo((4, 4, 4), grid)

        for (dir, metric) in zip(dirsize, metrics)
            @eval $(Symbol(:coeff_ , metric)) = calc_interpolating_coefficients($FT, $new_grid.$metric, $arch, $new_grid.$dir)
            @eval $(Symbol(:smooth_, metric)) = calc_smoothness_coefficients($FT, $Val($stretched_smoothness), $new_grid.$metric, $arch, $new_grid.$dir)
        end
    end

    XT = typeof(coeff_xᶠᵃᵃ)
    YT = typeof(coeff_yᵃᶠᵃ)
    ZT = typeof(coeff_zᵃᵃᶠ)
    XS = typeof(smooth_xᶠᵃᵃ)
    YS = typeof(smooth_yᵃᶠᵃ)
    ZS = typeof(smooth_zᵃᵃᶠ)
    zweno ? W = Number : W = Nothing

    return WENO5{FT, XT, YT, ZT, XS, YS, ZS, W}(coeff_xᶠᵃᵃ , coeff_xᶜᵃᵃ , coeff_yᵃᶠᵃ , coeff_yᵃᶜᵃ , coeff_zᵃᵃᶠ , coeff_zᵃᵃᶜ ,
                                                smooth_xᶠᵃᵃ, smooth_xᶜᵃᵃ, smooth_yᵃᶠᵃ, smooth_yᵃᶜᵃ, smooth_zᵃᵃᶠ, smooth_zᵃᵃᶜ)
end

const JSWENO = WENO5{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Nothing}
const ZWENO  = WENO5{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any}

function Base.show(io::IO, a::WENO5{FT, RX, RY, RZ}) where {FT, RX, RY, RZ}
    print(io, "WENO5 advection scheme with: \n",
              "    ├── X $(RX == Nothing ? "regular" : "stretched") \n",
              "    ├── Y $(RY == Nothing ? "regular" : "stretched") \n",
              "    └── Z $(RZ == Nothing ? "regular" : "stretched")" )
end

Adapt.adapt_structure(to, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS, W}) where {FT, XT, YT, ZT, XS, YS, ZS, W} =
     WENO5{FT, typeof(Adapt.adapt(to, scheme.coeff_xᶠᵃᵃ)),
               typeof(Adapt.adapt(to, scheme.coeff_yᵃᶠᵃ)),
               typeof(Adapt.adapt(to, scheme.coeff_zᵃᵃᶠ)),
               typeof(Adapt.adapt(to, scheme.smooth_xᶠᵃᵃ)),
               typeof(Adapt.adapt(to, scheme.smooth_yᵃᶠᵃ)),
               typeof(Adapt.adapt(to, scheme.smooth_zᵃᵃᶠ)), W}(
        Adapt.adapt(to, scheme.coeff_xᶠᵃᵃ),
        Adapt.adapt(to, scheme.coeff_xᶜᵃᵃ),
        Adapt.adapt(to, scheme.coeff_yᵃᶠᵃ),
        Adapt.adapt(to, scheme.coeff_yᵃᶜᵃ),
        Adapt.adapt(to, scheme.coeff_zᵃᵃᶠ),
        Adapt.adapt(to, scheme.coeff_zᵃᵃᶜ),
        Adapt.adapt(to, scheme.smooth_xᶠᵃᵃ),
        Adapt.adapt(to, scheme.smooth_xᶜᵃᵃ),
        Adapt.adapt(to, scheme.smooth_yᵃᶠᵃ),
        Adapt.adapt(to, scheme.smooth_yᵃᶜᵃ),
        Adapt.adapt(to, scheme.smooth_zᵃᵃᶠ),
        Adapt.adapt(to, scheme.smooth_zᵃᵃᶜ))

@inline boundary_buffer(::WENO5) = 2

@inline symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, centered_fourth_order, c)

@inline symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, ::WENO5, u) = symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, centered_fourth_order, u)
@inline symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, ::WENO5, v) = symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, centered_fourth_order, v)
@inline symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, ::WENO5, w) = symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, centered_fourth_order, w)

# Unroll the functions to pass the coordinates in case of a stretched grid
@inline left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_xᶠᵃᵃ(i, j, k, scheme, ψ, i, Face)
@inline left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_yᵃᶠᵃ(i, j, k, scheme, ψ, j, Face)
@inline left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_zᵃᵃᶠ(i, j, k, scheme, ψ, k, Face)

@inline right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_xᶠᵃᵃ(i, j, k, scheme, ψ, i, Face)
@inline right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_yᵃᶠᵃ(i, j, k, scheme, ψ, j, Face)
@inline right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_zᵃᵃᶠ(i, j, k, scheme, ψ, k, Face)

@inline left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_xᶠᵃᵃ(i+1, j, k, scheme, ψ, i, Center)
@inline left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_yᵃᶠᵃ(i, j+1, k, scheme, ψ, j, Center)
@inline left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ)  = weno_left_biased_interpolate_zᵃᵃᶠ(i, j, k+1, scheme, ψ, k, Center)

@inline right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_xᶠᵃᵃ(i+1, j, k, scheme, ψ, i, Center)
@inline right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_yᵃᶠᵃ(i, j+1, k, scheme, ψ, j, Center)
@inline right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ) = weno_right_biased_interpolate_zᵃᵃᶠ(i, j, k+1, scheme, ψ, k, Center)

# Stencil to calculate the stretched WENO weights and smoothness indicators
@inline left_stencil_x(i, j, k, ψ, args...) = @inbounds ( (ψ[i-3, j, k], ψ[i-2, j, k], ψ[i-1, j, k]), (ψ[i-2, j, k], ψ[i-1, j, k], ψ[i, j, k]), (ψ[i-1, j, k], ψ[i, j, k], ψ[i+1, j, k]) )
@inline left_stencil_y(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j-3, k], ψ[i, j-2, k], ψ[i, j-1, k]), (ψ[i, j-2, k], ψ[i, j-1, k], ψ[i, j, k]), (ψ[i, j-1, k], ψ[i, j, k], ψ[i, j+1, k]) )
@inline left_stencil_z(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j, k-3], ψ[i, j, k-2], ψ[i, j, k-1]), (ψ[i, j, k-2], ψ[i, j, k-1], ψ[i, j, k]), (ψ[i, j, k-1], ψ[i, j, k], ψ[i, j, k+1]) )

@inline right_stencil_x(i, j, k, ψ, args...) = @inbounds ( (ψ[i-2, j, k], ψ[i-1, j, k], ψ[i, j, k]), (ψ[i-1, j, k], ψ[i, j, k], ψ[i+1, j, k]), (ψ[i, j, k], ψ[i+1, j, k], ψ[i+2, j, k]) )
@inline right_stencil_y(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j-2, k], ψ[i, j-1, k], ψ[i, j, k]), (ψ[i, j-1, k], ψ[i, j, k], ψ[i, j+1, k]), (ψ[i, j, k], ψ[i, j+1, k], ψ[i, j+2, k]) )
@inline right_stencil_z(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j, k-2], ψ[i, j, k-1], ψ[i, j, k]), (ψ[i, j, k-1], ψ[i, j, k], ψ[i, j, k+1]), (ψ[i, j, k], ψ[i, j, k+1], ψ[i, j, k+2]) )

#####
##### biased pₖ for û calculation
#####

@inline left_biased_p₀(scheme, ψ, args...) = @inbounds sum(coeff_left_p₀(scheme, args...) .* ψ)
@inline left_biased_p₁(scheme, ψ, args...) = @inbounds sum(coeff_left_p₁(scheme, args...) .* ψ)
@inline left_biased_p₂(scheme, ψ, args...) = @inbounds sum(coeff_left_p₂(scheme, args...) .* ψ)

@inline right_biased_p₀(scheme, ψ, args...) = @inbounds sum(coeff_right_p₀(scheme, args...) .* ψ)
@inline right_biased_p₁(scheme, ψ, args...) = @inbounds sum(coeff_right_p₁(scheme, args...) .* ψ)
@inline right_biased_p₂(scheme, ψ, args...) = @inbounds sum(coeff_right_p₂(scheme, args...) .* ψ)

##₁###
##₂### Jiang & Shu (1996) WENO smoothness indicators. See also Equation 2.63 in Shu (1998)
#####

@inline left_biased_β₀(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * (3ψ[1] - 4ψ[2] +  ψ[3])^two_32
@inline left_biased_β₁(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1]         -  ψ[3])^two_32
@inline left_biased_β₂(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1] - 4ψ[2] + 3ψ[3])^two_32

@inline right_biased_β₀(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1] - 4ψ[2] + 3ψ[3])^two_32
@inline right_biased_β₁(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1]         -  ψ[3])^two_32
@inline right_biased_β₂(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * (3ψ[1] - 4ψ[2] +  ψ[3])^two_32

#####
##### Stretched smoothness indicators gathered from precomputed values.
##### The stretched values for β coefficients are calculated from
##### Shu, NASA/CR-97-206253, ICASE Report No. 97-65
##### by hardcoding that p(x) is a 2nd order polynomial
#####

@inline function biased_left_β(ψ, scheme, r, dir, i, location)
    @inbounds begin
        stencil = retrieve_left_smooth(scheme, r, dir, i, location)
        wᵢᵢ = stencil[1]
        wᵢⱼ = stencil[2]
        # horrible but have to do this for GPU execution (broadcast doesn't work apparently)
        result = 0
        for j = 1:3
            result += ψ[j] * ( wᵢᵢ[j] * ψ[j] + wᵢⱼ[j] * dagger(ψ)[j] )
        end
    end
    return result
end

@inline function biased_right_β(ψ, scheme, r, dir, i, location)
    @inbounds begin
        stencil = retrieve_right_smooth(scheme, r, dir, i, location)
        wᵢᵢ = stencil[1]
        wᵢⱼ = stencil[2]
        # horrible but have to do this for GPU execution (broadcast doesn't work apparently sum(ψ.*(wᵢᵢ.*ψ.+wᵢⱼ.*dagger(ψ))) )
        result = 0
        for j = 1:3
            result += ψ[j] * ( wᵢᵢ[j] * ψ[j] + wᵢⱼ[j] * dagger(ψ)[j] )
        end
    end
    return result
end

@inline left_biased_β₀(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 0, args...)
@inline left_biased_β₁(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 1, args...)
@inline left_biased_β₂(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 2, args...)

@inline right_biased_β₀(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 2, args...)
@inline right_biased_β₁(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 1, args...)
@inline right_biased_β₂(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 0, args...)

# Right-biased smoothness indicators are a reflection or "symmetric modification" of the left-biased smoothness
# indicators around grid point `i-1/2`.

@inline left_biased_α₀(FT, ψ, T, scheme, args...) = FT(C3₀) / (left_biased_β₀(FT, ψ, T, scheme, args...) + FT(ε))^ƞ
@inline left_biased_α₁(FT, ψ, T, scheme, args...) = FT(C3₁) / (left_biased_β₁(FT, ψ, T, scheme, args...) + FT(ε))^ƞ
@inline left_biased_α₂(FT, ψ, T, scheme, args...) = FT(C3₂) / (left_biased_β₂(FT, ψ, T, scheme, args...) + FT(ε))^ƞ

@inline right_biased_α₀(FT, ψ, T, scheme, args...) = FT(C3₂) / (right_biased_β₀(FT, ψ, T, scheme, args...) + FT(ε))^ƞ
@inline right_biased_α₁(FT, ψ, T, scheme, args...) = FT(C3₁) / (right_biased_β₁(FT, ψ, T, scheme, args...) + FT(ε))^ƞ
@inline right_biased_α₂(FT, ψ, T, scheme, args...) = FT(C3₀) / (right_biased_β₂(FT, ψ, T, scheme, args...) + FT(ε))^ƞ

#####
##### Z-WENO-5 reconstruction (Castro et al: High order weighted essentially non-oscillatory WENO-Z schemesfor hyperbolic conservation laws)
#####

@inline function left_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, T, scheme::ZWENO, args...)
    β₀ = left_biased_β₀(FT, ψ₀, T, scheme, args...)
    β₁ = left_biased_β₁(FT, ψ₁, T, scheme, args...)
    β₂ = left_biased_β₂(FT, ψ₂, T, scheme, args...)

    τ₅ = abs(β₂ - β₀)
    α₀ = FT(C3₀) * (1 + (τ₅ / (β₀ + FT(ε)))^ƞ)
    α₁ = FT(C3₁) * (1 + (τ₅ / (β₁ + FT(ε)))^ƞ)
    α₂ = FT(C3₂) * (1 + (τ₅ / (β₂ + FT(ε)))^ƞ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function right_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, T, scheme::ZWENO, args...)
    β₀ = right_biased_β₀(FT, ψ₀, T, scheme, args...)
    β₁ = right_biased_β₁(FT, ψ₁, T, scheme, args...)
    β₂ = right_biased_β₂(FT, ψ₂, T, scheme, args...)

    τ₅ = abs(β₂ - β₀)
    α₀ = FT(C3₂) * (1 + (τ₅ / (β₀ + FT(ε)))^ƞ)
    α₁ = FT(C3₁) * (1 + (τ₅ / (β₁ + FT(ε)))^ƞ)
    α₂ = FT(C3₀) * (1 + (τ₅ / (β₂ + FT(ε)))^ƞ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

#####
##### JS-WENO-5 reconstruction
#####

@inline function left_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, T, scheme::JSWENO, args...)
    α₀ = left_biased_α₀(FT, ψ₀, T, scheme, args...)
    α₁ = left_biased_α₁(FT, ψ₁, T, scheme, args...)
    α₂ = left_biased_α₂(FT, ψ₂, T, scheme, args...)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function right_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, T, scheme::JSWENO, args...)
    α₀ = right_biased_α₀(FT, ψ₀, T, scheme, args...)
    α₁ = right_biased_α₁(FT, ψ₁, T, scheme, args...)
    α₂ = right_biased_α₂(FT, ψ₂, T, scheme, args...)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

#####
##### Biased interpolation functions
#####

@inline function weno_left_biased_interpolate_xᶠᵃᵃ(i, j, k, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = left_stencil_x(i, j, k, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, XS, scheme, Val(1), args...)
    return w₀ * left_biased_p₀(scheme, ψ₀, XT, Val(1), args...) +
           w₁ * left_biased_p₁(scheme, ψ₁, XT, Val(1), args...) +
           w₂ * left_biased_p₂(scheme, ψ₂, XT, Val(1), args...)
end

@inline function weno_left_biased_interpolate_yᵃᶠᵃ(i, j, k, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = left_stencil_y(i, j, k, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, YS, scheme, Val(2), args...)
    return w₀ * left_biased_p₀(scheme, ψ₀, YT, Val(2), args...) +
           w₁ * left_biased_p₁(scheme, ψ₁, YT, Val(2), args...) +
           w₂ * left_biased_p₂(scheme, ψ₂, YT, Val(2), args...)
end

@inline function weno_left_biased_interpolate_zᵃᵃᶠ(i, j, k, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = left_stencil_z(i, j, k, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, ZS, scheme, Val(3), args...)
    return w₀ * left_biased_p₀(scheme, ψ₀, ZT, Val(3), args...) +
           w₁ * left_biased_p₁(scheme, ψ₁, ZT, Val(3), args...) +
           w₂ * left_biased_p₂(scheme, ψ₂, ZT, Val(3), args...)
end

@inline function weno_right_biased_interpolate_xᶠᵃᵃ(i, j, k, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = right_stencil_x(i, j, k, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, XS, scheme, Val(1), args...)
    return w₀ * right_biased_p₀(scheme, ψ₀, XT, Val(1), args...) +
           w₁ * right_biased_p₁(scheme, ψ₁, XT, Val(1), args...) +
           w₂ * right_biased_p₂(scheme, ψ₂, XT, Val(1), args...)
end

@inline function weno_right_biased_interpolate_yᵃᶠᵃ(i, j, k, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = right_stencil_y(i, j, k, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, YS, scheme, Val(2), args...)
    return w₀ * right_biased_p₀(scheme, ψ₀, YT, Val(2), args...) +
           w₁ * right_biased_p₁(scheme, ψ₁, YT, Val(2), args...) +
           w₂ * right_biased_p₂(scheme, ψ₂, YT, Val(2), args...)
end

@inline function weno_right_biased_interpolate_zᵃᵃᶠ(i, j, k,scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS}, ψ, args...) where {FT, XT, YT, ZT, XS, YS, ZS}
    ψ₂, ψ₁, ψ₀ = right_stencil_z(i, j, k, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights(FT, ψ₂, ψ₁, ψ₀, ZS, scheme, Val(3), args...)
    return w₀ * right_biased_p₀(scheme, ψ₀, ZT, Val(3), args...) +
           w₁ * right_biased_p₁(scheme, ψ₁, ZT, Val(3), args...) +
           w₂ * right_biased_p₂(scheme, ψ₂, ZT, Val(3), args...)
end

#####
##### Coefficients for stretched (and uniform) ENO schemes (see Shu NASA/CR-97-206253, ICASE Report No. 97-65)
#####

@inline coeff_left_p₀(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (  FT(1/3),    FT(5/6), - FT(1/6))
@inline coeff_left_p₁(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (- FT(1/6),    FT(5/6),   FT(1/3))
@inline coeff_left_p₂(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (  FT(1/3),  - FT(7/6),  FT(11/6))

@inline coeff_right_p₀(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₂(scheme, Nothing, args...))
@inline coeff_right_p₁(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₁(scheme, Nothing, args...))
@inline coeff_right_p₂(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₀(scheme, Nothing, args...))

@inline coeff_left_p₀(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 0, dir, i ,loc)
@inline coeff_left_p₁(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 1, dir, i ,loc)
@inline coeff_left_p₂(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 2, dir, i ,loc)

@inline coeff_right_p₀(scheme, T, dir, i, loc) = retrieve_coeff(scheme, -1, dir, i ,loc)
@inline coeff_right_p₁(scheme, T, dir, i, loc) = retrieve_coeff(scheme,  0, dir, i ,loc)
@inline coeff_right_p₂(scheme, T, dir, i, loc) = retrieve_coeff(scheme,  1, dir, i ,loc)

@inline retrieve_coeff(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.coeff_xᶠᵃᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.coeff_xᶜᵃᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.coeff_yᵃᶠᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.coeff_yᵃᶜᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.coeff_zᵃᵃᶠ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.coeff_zᵃᵃᶜ[r+2][i]

@inline retrieve_left_smooth(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.smooth_xᶠᵃᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.smooth_xᶜᵃᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.smooth_yᵃᶠᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.smooth_yᵃᶜᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.smooth_zᵃᵃᶠ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.smooth_zᵃᵃᶜ[r+1][i]

@inline retrieve_right_smooth(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.smooth_xᶠᵃᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.smooth_xᶜᵃᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.smooth_yᵃᶠᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.smooth_yᵃᶜᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.smooth_zᵃᵃᶠ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.smooth_zᵃᵃᶜ[r+4][i]

@inline calc_interpolating_coefficients(FT, coord::OffsetArray{<:Any, <:Any, <:AbstractRange}, arch, N) = nothing
@inline calc_interpolating_coefficients(FT, coord::AbstractRange, arch, N)                              = nothing

@inline calc_smoothness_coefficients(FT, ::Val{false}, args...) = nothing
@inline calc_smoothness_coefficients(FT, ::Val{true}, coord::OffsetArray{<:Any, <:Any, <:AbstractRange}, arch, N) = nothing
@inline calc_smoothness_coefficients(FT, ::Val{true}, coord::AbstractRange, arch, N) = nothing

function calc_interpolating_coefficients(FT, coord, arch, N)

    cpu_coord = Array(parent(coord))
    cpu_coord = OffsetArray(cpu_coord, coord.offsets[1])

    s1 = create_interp_coefficients(FT,-1, cpu_coord, arch, N)
    s2 = create_interp_coefficients(FT, 0, cpu_coord, arch, N)
    s3 = create_interp_coefficients(FT, 1, cpu_coord, arch, N)
    s4 = create_interp_coefficients(FT, 2, cpu_coord, arch, N)

    return (s1, s2, s3, s4)
end

function create_interp_coefficients(FT, r, cpu_coord, arch, N)

    stencil = NTuple{3, FT}[]
    @inbounds begin
        for i = 0:N+1
            push!(stencil, interp_weights(r, cpu_coord, i, 0, -))
        end
    end
    return OffsetArray(arch_array(arch, stencil), -1)
end

function calc_smoothness_coefficients(FT, beta, coord, arch, N)

    cpu_coord = Array(parent(coord))
    cpu_coord = OffsetArray(cpu_coord, coord.offsets[1])

    s1 = create_smoothness_coefficients(FT, 0, -, cpu_coord, arch, N)
    s2 = create_smoothness_coefficients(FT, 1, -, cpu_coord, arch, N)
    s3 = create_smoothness_coefficients(FT, 2, -, cpu_coord, arch, N)
    s4 = create_smoothness_coefficients(FT, 0, +, cpu_coord, arch, N)
    s5 = create_smoothness_coefficients(FT, 1, +, cpu_coord, arch, N)
    s6 = create_smoothness_coefficients(FT, 2, +, cpu_coord, arch, N)

    return (s1, s2, s3, s4, s5, s6)
end

function create_smoothness_coefficients(FT, r, op, cpu_coord, arch, N)

    # derivation written on overleaf

    stencil = NTuple{2, NTuple{3, FT}}[]
    @inbounds begin
        for i = 0:N+1

            bias1 = Int(op == +)
            bias2 = bias1 - 1

            Δcᵢ = cpu_coord[i + bias1] - cpu_coord[i + bias2]

            Bᵢ  = prim_interp_weights(r, cpu_coord, i, bias1, op)
            bᵢ  =      interp_weights(r, cpu_coord, i, bias1, op)
            bₓᵢ = der1_interp_weights(r, cpu_coord, i, bias1, op)
            Aᵢ  = prim_interp_weights(r, cpu_coord, i, bias2, op)
            aᵢ  =      interp_weights(r, cpu_coord, i, bias2, op)
            aₓᵢ = der1_interp_weights(r, cpu_coord, i, bias2, op)

            pₓₓ = der2_interp_weights(r, cpu_coord, i, op)
            Pᵢ  =  (Bᵢ .- Aᵢ)

            wᵢᵢ = Δcᵢ  .* (bᵢ .* bₓᵢ .- aᵢ .* aₓᵢ .- pₓₓ .* Pᵢ)  .+ Δcᵢ^4 .* (pₓₓ .* pₓₓ)
            wᵢⱼ = Δcᵢ  .* (star(bᵢ, bₓᵢ)  .- star(aᵢ, aₓᵢ) .- star(pₓₓ, Pᵢ)) .+
                                                 Δcᵢ^4 .* star(pₓₓ, pₓₓ)

            push!(stencil, (wᵢᵢ, wᵢⱼ))
        end
    end

    return OffsetArray(arch_array(arch, stencil), -1)
end

@inline dagger(ψ)    = (ψ[2], ψ[3], ψ[1])
@inline star(ψ₁, ψ₂) = (ψ₁ .* dagger(ψ₂) .+ dagger(ψ₁) .* ψ₂)

# Integral of ENO coefficients for 2nd order polynomial reconstruction at the face
function prim_interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        prod = 1
                        sum  = 0
                        for q = 0:3
                            if q != m && q != l
                                prod *= coord[op(i, r-q+1)]
                                sum  += coord[op(i, r-q+1)]
                            end
                        end
                        num += coord[i+bias]^3 / 3 - sum * coord[i+bias]^2 / 2 + prod * coord[i+bias]
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# Second derivative of ENO coefficients for 2nd order polynomial reconstruction at the face
function der2_interp_weights(r, coord, i, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        num += 2
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# first derivative of ENO coefficients for 2nd order polynomial reconstruction at the face
function der1_interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        sum = 0
                        for q = 0:3
                            if q != m && q != l
                                sum += coord[op(i, r-q+1)]
                            end
                        end
                        num += 2 * coord[i+bias] - sum
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# ENO coefficients for 2nd order polynomial reconstruction at the face
function interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        prod = 1
                        for q = 0:3
                            if q != m && q != l
                                prod *= (coord[i+bias] - coord[op(i, r-q+1)])
                            end
                        end
                        num += prod
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end


text:
# general zlib constants, definitions

@static if Sys.iswindows() && VERSION < v"1.3.0-rc2.9"
    # The name of the binary for Windows changed with
    # https://github.com/JuliaLang/julia/pull/33125
    const _zlib = "zlib1"
else
    const _zlib = "libz"
end

# Constants

zlib_version =  unsafe_string(ccall((:zlibVersion, _zlib), Ptr{UInt8}, ()))
ZLIB_VERSION = tuple([parse(Int, c) for c in split(zlib_version, '.')]...)

# Flush values
const Z_NO_FLUSH       =  Int32(0)
const Z_PARTIAL_FLUSH  =  Int32(1)
const Z_SYNC_FLUSH     =  Int32(2)
const Z_FULL_FLUSH     =  Int32(3)
const Z_FINISH         =  Int32(4)
const Z_BLOCK          =  Int32(5)
const Z_TREES          =  Int32(6)

# Return codes
const Z_OK             =  Int32(0)
const Z_STREAM_END     =  Int32(1)
const Z_NEED_DICT      =  Int32(2)
const Z_ERRNO          =  Int32(-1)
const Z_STREAM_ERROR   =  Int32(-2)
const Z_DATA_ERROR     =  Int32(-3)
const Z_MEM_ERROR      =  Int32(-4)
const Z_BUF_ERROR      =  Int32(-5)
const Z_VERSION_ERROR  =  Int32(-6)

# Zlib errors as Exceptions
zerror(e::Integer) =  unsafe_string(ccall((:zError, _zlib), Ptr{UInt8}, (Int32,), e))
mutable struct ZError <: Exception
    err::Int32
    err_str::AbstractString

    ZError(e::Integer) = (e == Z_ERRNO ? new(e, strerror()) : new(e, zerror(e)))
end

# Compression Levels
const Z_NO_COMPRESSION      =  Int32(0)
const Z_BEST_SPEED          =  Int32(1)
const Z_BEST_COMPRESSION    =  Int32(9)
const Z_DEFAULT_COMPRESSION =  Int32(-1)

# Compression Strategy
const Z_FILTERED             =  Int32(1)
const Z_HUFFMAN_ONLY         =  Int32(2)
const Z_RLE                  =  Int32(3)
const Z_FIXED                =  Int32(4)
const Z_DEFAULT_STRATEGY     =  Int32(0)

# data_type field
const Z_BINARY    =  Int32(0)
const Z_TEXT      =  Int32(1)
const Z_ASCII     = Z_TEXT
const Z_UNKNOWN   =  Int32(2)

# deflate compression method
const Z_DEFLATED    =  Int32(8)

# misc
const Z_NULL   = C_NULL

# Constants for use with gzbuffer
const Z_DEFAULT_BUFSIZE = 8192
const Z_BIG_BUFSIZE = 131072

# Constants for use with gzseek
const SEEK_SET =  Int32(0)
const SEEK_CUR =  Int32(1)

# Create ZFileOffset alias
# Use 64bit if the *64 functions are available or zlib is compiles with 64bit
# file offset.

# Get compile-time option flags
const zlib_compile_flags = ccall((:zlibCompileFlags, _zlib), UInt, ())
const z_off_t_sz = 2 << ((zlib_compile_flags >> 6) & UInt(3))
if z_off_t_sz == 8 || (!Sys.iswindows() && Libdl.dlsym_e(Libdl.dlopen(_zlib), :gzopen64) != C_NULL)
    const ZFileOffset = Int64
elseif z_off_t_sz == 4
    const ZFileOffset = Int32
else
    error("Can't figure out what to do with ZFileOffset. sizeof(z_off_t) = ", z_off_t_sz)
end


text:
using ForwardDiff

abstract type JacobianParameters{T} end

struct JacobianParametersUser{T, JT <: Function} <: JacobianParameters{T}
    J::JT
end

struct JacobianParametersAD{T, FT <: Function, JT <: ForwardDiff.JacobianConfig} <: JacobianParameters{T}
    F!::FT
    Jconfig::JT
    tx::Vector{T}
    ty::Vector{T}
end

function JacobianParametersAD(F!::FT, Jconfig::JT, tx::Vector{T}, ty::Vector{T}) where {T, FT, JT}
    @assert length(tx) == length(ty)
    JacobianParametersAD{T, FT, JT}(F!, Jconfig, tx, ty)
end

function JacobianParametersAD(F!::FT, T, n::Int) where {FT <: Function}
    F!rev = (y,x) -> F!(x,y)
    tx = zeros(T, n)
    ty = zeros(T, n)
    Jconfig = ForwardDiff.JacobianConfig(F!rev, ty, tx)
    Jparams = JacobianParametersAD(F!rev, Jconfig, tx, ty)
end

struct JacobianParametersFD{T, FT} <: JacobianParameters{T}
    ϵ::T
    F!::FT
    f1::Vector{T}
    f2::Vector{T}
    e::Vector{T}
    tx::Vector{T}
end

function JacobianParametersFD(F!::FT, ϵ, T, n::Int) where {FT <: Function}
    f1 = zeros(T, n)
    f2 = zeros(T, n)
    e  = zeros(T, n)
    tx = zeros(T, n)
    Jparams = JacobianParametersFD{T, typeof(F!)}(ϵ, F!, f1, f2, e, tx)
end

function computeJacobian(x::Vector{T}, J::Matrix{T}, params::JacobianParametersUser{T}) where {T}
    params.J(x, J)
end

function computeJacobian(x::Vector{T}, J::Matrix{T}, params::JacobianParametersFD{T}) where {T}
    @assert length(x) == size(J, 1) == size(J, 2)

    local ϵⱼ::T

    for j in eachindex(x)
        ϵⱼ = params.ϵ * x[j] + params.ϵ
        fill!(params.e, 0)
        params.e[j] = 1
        params.tx .= x .- ϵⱼ .* params.e
        params.F!(params.tx, params.f1)
        params.tx .= x .+ ϵⱼ .* params.e
        params.F!(params.tx, params.f2)
        for i in eachindex(x)
            J[i,j] = (params.f2[i]-params.f1[i])/(2ϵⱼ)
        end
    end
end

function computeJacobian(x::Vector{T}, J::Matrix{T}, Jparams::JacobianParametersAD{T}) where {T}
    ForwardDiff.jacobian!(J, Jparams.F!, Jparams.ty, x, Jparams.Jconfig)
end

function computeJacobianFD(x::Vector{T}, J::Matrix{T}, F!::FT, ϵ::T) where{T, FT <: Function}
    params = JacobianParametersFD(F!, ϵ, T, length(x))
    computeJacobian(x, J, params)
end

function computeJacobianAD(x::Vector{T}, J::Matrix{T}, F!::FT) where{T, FT <: Function}
    params = JacobianParametersAD(F!, T, length(x))
    computeJacobian(x, J, params)
end

function getJacobianParameters(J, F!, T, n)
    if J == nothing
        if get_config(:jacobian_autodiff)
            Jparams = JacobianParametersAD(F!, T, n)
        else
            Jparams = JacobianParametersFD(F!, get_config(:jacobian_fd_ϵ), T, n)
        end
    else
        Jparams = JacobianParametersUser{T, typeof(J)}(J)
    end
    return Jparams
end