text:
@testset "ElasticDistortion" begin
    @test (ElasticDistortion <: Augmentor.AffineOperation) == false
    @test (ElasticDistortion <: Augmentor.ImageOperation) == true

    @testset "constructor" begin
        @test_throws MethodError ElasticDistortion()
        let op = @inferred ElasticDistortion(3)
            @test op.gridheight == 3
            @test op.gridwidth == 3
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 3)"
            @test str_showconst(op) == "ElasticDistortion(3, 3)"
            @test str_showcompact(op) == "Distort using a smoothed and normalized 3×3 grid with pinned border"
        end
        let op = @inferred ElasticDistortion(3, 4)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4)"
            @test str_showconst(op) == "ElasticDistortion(3, 4)"
            @test str_showcompact(op) == "Distort using a smoothed and normalized 3×4 grid with pinned border"
        end
        let op = @inferred ElasticDistortion(3, 4, 0.6)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6)"
            @test str_showcompact(op) == "Distort using a smoothed and normalized 3×4 grid with pinned border"
        end
        let op = @inferred ElasticDistortion(3, 4, 0.6, 0.2)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test op.sigma == 0.2
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6, sigma=0.2)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6, sigma=0.2)"
            @test str_showcompact(op) == "Distort using a smoothed and normalized 3×4 grid with pinned border"
        end
        let op = @inferred ElasticDistortion(3, 4, 0.6, 0.2, 2, false, true)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test op.sigma == 0.2
            @test op.iterations == 2
            @test op.border == false
            @test op.normalize == true
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2)"
            @test str_showcompact(op) == "Distort using a 2-times smoothed and normalized 3×4 grid with pinned border"
        end
        let op = @inferred ElasticDistortion(3, 4, 0.6, 0.2, 1, true, false)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test op.sigma == 0.2
            @test op.iterations == 1
            @test op.border == true
            @test op.normalize == false
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6, sigma=0.2, border=true, norm=false)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6, sigma=0.2, border=true, norm=false)"
            @test str_showcompact(op) == "Distort using a smoothed 3×4 grid"
        end
        let op = ElasticDistortion(3, 4, 0.6, sigma=0.2, iter=2, norm=false, border=false)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test op.sigma == 0.2
            @test op.iterations == 2
            @test op.border == false
            @test op.normalize == false
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2, norm=false)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2, norm=false)"
            @test str_showcompact(op) == "Distort using a 2-times smoothed 3×4 grid with pinned border"
        end
        let op = ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2, norm=false, border=false)
            @test op.gridheight == 3
            @test op.gridwidth == 4
            @test op.scale == 0.6
            @test op.sigma == 0.2
            @test op.iterations == 2
            @test op.border == false
            @test op.normalize == false
            @test str_show(op) == "Augmentor.ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2, norm=false)"
            @test str_showconst(op) == "ElasticDistortion(3, 4, scale=0.6, sigma=0.2, iter=2, norm=false)"
            @test str_showcompact(op) == "Distort using a 2-times smoothed 3×4 grid with pinned border"
        end
    end
    imgs = [
        (rect),
        (Augmentor.prepareaffine(rect)),
        (OffsetArray(rect, -2, -1)),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3))),
    ]
    @testset "eager" begin
        @test Augmentor.supports_eager(ElasticDistortion) === false
        # TODO: better tests
        @testset "single image" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(ElasticDistortion(4,4), img_in))
                @test size(res) == size(rect)
                @test typeof(res) <: Array{eltype(img_in),2}
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_in = (img_in1, img_in2)
                res = @inferred(Augmentor.applyeager(ElasticDistortion(4,4), img_in))
                @test res[1] == res[2]
                @test typeof(res) <: NTuple{2,Array{eltype(img_in1),2}}
            end
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(ElasticDistortion) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(ElasticDistortion) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(ElasticDistortion) === true
        @testset "single image" begin
            for img in imgs
                res = @inferred(Augmentor.applylazy(ElasticDistortion(4,4), img))
                @test size(res) == size(rect)
                @test typeof(res) <: Augmentor.DistortedView{eltype(rect)}
                @test parent(res) === img
            end
        end
        @testset "multiple images" begin
            for img_in1 in imgs, img_in2 in imgs
                img_in = (img_in1, img_in2)
                res = @inferred(Augmentor.applylazy(ElasticDistortion(4,4), img_in))
                @test res[1] == res[2]
                @test res[1].field.itp.coefs === res[2].field.itp.coefs
                @test typeof(res) <: NTuple{2,Augmentor.DistortedView{eltype(img_in1)}}
            end
        end
    end
    @testset "view" begin
        @test Augmentor.supports_view(ElasticDistortion) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(ElasticDistortion) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(ElasticDistortion) === false
    end
end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

pArr = [100, 200, 500]
elemArr = [(5,5), (8, 7), (50, 25), (21, 20), (30, 30)]
n = 300
est      = Array{Any}(undef, 9)   # number of methods

rep   = parse(Int,ARGS[1])
ip    = parse(Int,ARGS[2])
iElem = parse(Int,ARGS[3])
dir = ARGS[4]

p = pArr[ip]
Random.seed!(1234)

# generate model
Ωx = Matrix{Float64}(I, p, p)
ρ = 0.7
for k=1:p-1
    for l=1:p-k
        Ωx[l  , l+k] = ρ^k
        Ωx[l+k, l  ] = ρ^k
    end
end
Δ = zeros(p, p)
# generate Delta
for j=1:5
    i = 5 + (j-1)
    Δ[i, i] = rand(Uniform(0.1, 0.2))
end
for j=1:5
    i = 5 + (j-1)
    v = rand(Uniform(0.2, 0.5))
    Δ[i  , i+1] = v
    Δ[i+1, i  ] = v
end
Ωy = Ωx - Δ
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σx = inv(Symmetric(D * Ωx * D))
Σy = inv(Symmetric(D * Ωy * D))

tΔ = D * Δ * D

dist_X = MvNormal(convert(Matrix, Σx))
dist_Y = MvNormal(convert(Matrix, Σy))

# generate data
Random.seed!(1234 + rep)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

@show (ip, iElem)
supportOracle =  LinearIndices(tΔ)[findall(!iszero, tΔ)]
ri, ci = elemArr[iElem]

@time est[1], eΔ, esuppΔ, eω, esuppω, eSupport = DiffPrecTest.estimate(SymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy)
@time est[2], _, _, _, _, _ = DiffPrecTest.estimate(AsymmetricNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy, Δ=eΔ, suppΔ=esuppΔ)
@time est[3] = DiffPrecTest.estimate(SeparateNormal(), X, Y, ri, ci; Sx=Sx, Sy=Sy)
@time est[4] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport)
@time est[5] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, eSupport)

@time est[6] = DiffPrecTest.estimate(SymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[7] = DiffPrecTest.estimate(AsymmetricOracleNormal(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[8] = DiffPrecTest.estimate(SymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle)
@time est[9] = DiffPrecTest.estimate(AsymmetricOracleBoot(), Sx, Sy, X, Y, ri, ci, supportOracle)

@save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est eΔ esuppΔ eω esuppω eSupport

ext: jl


text:
function sbt(x)
  K = length(x)
  offset = (K + 1) .- cumsum(ones(K), dims=1)
  z = sigmoid.(x .- log.(offset))
  z_cumprod = cumprod(1 .- z, dims=1)
  # y = cat(z, pad, dims=dims) .* cat(pad, z_cumprod, dims=1)
  y = vcat(z, 1) .* vcat(1, z_cumprod)

  return y
end

function sbt_inv(p)
  # see: https://pytorch.org/docs/stable/_modules/torch/distributions/transforms.html
  error("unimplemented")
end

function sbt_logabsdetJ(x, p)
  K = length(x)

  offset = (K + 1) .- cumsum(ones(K), dims=1)
  z = sigmoid.(x .- log.(offset))
  detJ = sum(log.(1 .- z) .+ log.(p[1:end-1]))

  return detJ
end

#= TEST
x = param(randn(3))
p = sbt(x)
sbt_logabsdetJ(x, p)
logpdf(Dirichlet(ones(4)), p) + sbt_logabsdetJ(x, p)
=#


text:
struct RS{bits} <: Encoding
    key::MbedTLS.PKContext
    is_private_key::Bool
end

alg(::RS{bits}) where {bits} = "RS$(bits)"

function Base.show(io::IO, encoding::RS)
    print(io, alg(encoding))

    if encoding.is_private_key
        print(io, " Private Key")
    else
        print(io, " Public Key")
    end
end

"RSASSA using SHA-256 hash algorithm"
const RS256 = RS{256}

"RSASSA using SHA-384 hash algorithm"
const RS384 = RS{384}

md_hash_alg(::RS256) = MbedTLS.MD_SHA256
md_hash_alg(::RS384) = MbedTLS.MD_SHA384

function _try_base64decode(str::AbstractString) :: Union{Nothing, String}
    try
        return String(Base64.base64decode(str))
    catch
        return nothing
    end
end

function _try_isfile(str::AbstractString) :: Bool
    try
        return isfile(str)
    catch
        return false
    end
end

@inline function has_key_prefix(str::AbstractString)
    return startswith(str, "-----BEGIN PUBLIC KEY-----") || startswith(str, "-----BEGIN RSA PRIVATE KEY-----")
end

@inline convert_string_to_bytes(str::AbstractString) :: Vector{UInt8} = convert(Vector{UInt8}, codeunits(str))

"""
    RS{bits}(key_or_filepath::AbstractString) where {bits}

`key_or_filepath` can be either the key content as plain text or base64 encoded string,
or the filepath to the key file.
"""
function RS{bits}(key_or_filepath::AbstractString) where {bits}
    local key_as_bytes::Vector{UInt8}
    local key_as_string::String

    error_msg = "$key_or_filepath is not a valid RSA public or private key."

    if has_key_prefix(key_or_filepath)
        # plain key string
        key_as_string = String(key_or_filepath)
        key_as_bytes = convert_string_to_bytes(key_as_string)
    else
        # base64 encoded key string
        decoded_key = _try_base64decode(key_or_filepath)
        if (decoded_key != nothing) && (has_key_prefix(decoded_key))
            key_as_string = String(decoded_key)
            key_as_bytes = convert_string_to_bytes(decoded_key)
        elseif _try_isfile(key_or_filepath)
            # filepath
            key_as_bytes = read(open(key_or_filepath, "r"))
            key_as_string = String(copy(key_as_bytes))
        else
            throw(ArgumentError(error_msg))
        end
    end

    context = MbedTLS.PKContext()

    if startswith(key_as_string, "-----BEGIN PUBLIC KEY-----")
        # public key
        MbedTLS.parse_public_key!(context, key_as_bytes)
        return RS{bits}(context, false)
    elseif startswith(key_as_string, "-----BEGIN RSA PRIVATE KEY-----")
        # private key
        MbedTLS.parse_key!(context, key_as_bytes)
        return RS{bits}(context, true)
    else
        throw(ArgumentError(error_msg))
    end
end

function sign(encoding::RS, data)
    @assert encoding.is_private_key "Must sign using a private key."
    md = md_hash_alg(encoding)
    signature = MbedTLS.sign(encoding.key, md, MbedTLS.digest(md, data), MersenneTwister(0))
    return base64url_encode(signature)
end

function has_valid_signature(encoding::RS, header_and_claims_encoded::AbstractString, signature_encoded::AbstractString) :: Bool
    try
        md = md_hash_alg(encoding)
        _hash = MbedTLS.digest(md, header_and_claims_encoded)
        return MbedTLS.verify(encoding.key, md, _hash, base64url_decode(signature_encoded)) == 0
    catch e
        return false
    end
end


text:
function checkJL(expr, val)
	def = expr
	def = replace(def, r"(\w+)::(\w+):(\d+)" => s"(\1:\3)::\2")

	X = @eval module $(gensym())
		using CBinding
		c`-std=c99`
		c"""
		#include <stdalign.h>
		"""s
		const X = @c_str $(expr)
	end

	x = Base.invokelatest(X.X)
	result = []
	push!(result, "$(string(sizeof(x))): $(bytes2hex(UInt8[getfield(x, :mem)...,]))")
	types = Base.invokelatest(fieldtypes, X.X)
	for (ind, prop) in enumerate(Base.invokelatest(fieldnames, X.X))
		v1 = sizeof(types[ind]) < sizeof(val) ? Core.Intrinsics.trunc_int(types[ind], val) : reinterpret(types[ind], val)
		x = Base.invokelatest(X.X; (prop => v1,)...)
		v2 = Base.invokelatest(getproperty, x, prop)
		x = Base.invokelatest(X.X; (prop => v2,)...)
		v3 = Base.invokelatest(getproperty, x, prop)
		@test v2 === v3
		push!(result, "$(String(prop)): $(bytes2hex(UInt8[getfield(x, :mem)...,]))")
	end
	return result
end

function checkC(expr, val)
	props = []
	for line in split(expr, '\n')
		m = match(r"^\s*(?:_Alignas\S+\s+)?((?:(?:char|short|int|long|float|double|}|}\s+__attribute__\S+)\s+)+)([cexsilfd]+)[:\d]*;\s*$", line)
		isnothing(m) || push!(props, m.captures[2] => strip(m.captures[1]))
	end
	m = match(r"^\s*(struct|union)\s+(\S+)", expr)
	(kind, name) = m.captures[end-1:end]

	use = []
	push!(use, """
		$(kind) $(name) x; memset(&x, 0, sizeof(x));
		printf("%ld: ", sizeof(x));
		for (int i = 0;i < sizeof(x);i++) printf("%02x", (unsigned int)(*(((unsigned char*)&x)+i)));
		printf("\\n");
		"""
	)
	for (prop, typ) in props
		setprop = "x.$(prop)"
		if typ == "float"
			setprop = "*((unsigned int *)&($(setprop))) = $(repr(Core.Intrinsics.trunc_int(UInt64, val)));"
		elseif typ == "double"
			setprop = "*((unsigned long long *)&($(setprop))) = $(repr(Core.Intrinsics.trunc_int(UInt64, val)));"
		elseif typ == "long double"
			setprop = """
				*(((unsigned long long *)&($(setprop)))+0) = $(repr(Core.Intrinsics.trunc_int(UInt64, val)));
				*(((unsigned long long *)&($(setprop)))+1) = $(repr(Core.Intrinsics.trunc_int(UInt64, val>>64)));"""
		else
			setprop = "$(setprop) = $(repr(Core.Intrinsics.trunc_int(UInt64, val)));"
		end
		push!(use, """
			$(kind) $(name) x; memset(&x, 0, sizeof(x));
			$(setprop)
			printf("%s: ", $(repr(prop)));
			for (int i = 0;i < sizeof(x);i++) printf("%02x", (unsigned int)(*(((unsigned char*)&x)+i)));
			printf("\\n");
			"""
		)
	end
	use = "{\n"*join(use, "}\n{\n")*"}\n"

	code = """
		#include <stdio.h>
		#include <stdalign.h>
		#include <string.h>

		$(expr)

		int main(int argc, char **argv) {
			$(use)

			return 0;
		}
		"""

	return CBinding.libclang.Clang_jll.clang() do bin
		tmp = tempname()
		open(f -> write(f, code), tmp*".c", "w+")
		run(`gcc -Wno-overflow -Wno-constant-conversion -Wno-address-of-packed-member -std=c99 -o $(tmp) $(tmp).c`)
		return readlines(`$(tmp)`)
	end
end

const LAYOUT_TESTS = [
"""struct T {
};""",

"""struct T {
	char c;
};""",

"""struct T {
	short s;
};""",

"""struct T {
	int i;
};""",

"""struct T {
	long l;
};""",

"""struct T {
	long long ll;
};""",

"""struct T {
	char x;
	char c;
};""",

"""struct T {
	char x;
	short s;
};""",

"""struct T {
	char x;
	int i;
};""",

"""struct T {
	char x;
	long l;
};""",

"""struct T {
	char x;
	long long ll;
};""",

"""struct T {
	char x;
	float f;
};""",

"""struct T {
	char x;
	double d;
};""",

"""struct T {
	char x;
	long double ld;
};""",

"""struct T {
	char c:3;
};""",

"""struct T {
	short s:3;
};""",

"""struct T {
	int i:3;
};""",

"""struct T {
	long l:3;
};""",

"""struct T {
	long long ll:3;
};""",

"""struct T {
	char x:3;
	char c;
};""",

"""struct T {
	char x:3;
	short s;
};""",

"""struct T {
	char x:3;
	int i;
};""",

"""struct T {
	char x:3;
	long l;
};""",

"""struct T {
	char x:3;
	long long ll;
};""",

"""struct T {
	char x:3;
	float f;
};""",

"""struct T {
	char x:3;
	double d;
};""",

"""struct T {
	char x:3;
	long double ld;
};""",

"""struct T {
	char c:3;
	char x;
};""",

"""struct T {
	short s:3;
	char x;
};""",

"""struct T {
	int i:3;
	char x;
};""",

"""struct T {
	long l:3;
	char x;
};""",

"""struct T {
	long long ll:3;
	char x;
};""",

"""struct T {
	char x:3;
	char c:3;
};""",

"""struct T {
	char x:3;
	short s:3;
};""",

"""struct T {
	char x:3;
	int i:3;
};""",

"""struct T {
	char x:3;
	long l:3;
};""",

"""struct T {
	char x:3;
	long long ll:3;
};""",

"""struct T {
	char x:3;
	char c:5;
};""",

"""struct T {
	char x:3;
	short s:13;
};""",

"""struct T {
	char x:3;
	int i:29;
};""",

"""struct T {
	char x:3;
	long l:29;
};""",

"""struct T {
	char x:3;
	long long ll:61;
};""",

"""struct T {
	char x:3;
	struct {
		char c:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		short s:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		int i:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		long l:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		long long ll:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		char c:3;
	};
};""",

"""struct T {
	short x:3;
	struct {
		short s:3;
	};
};""",

"""struct T {
	int x:3;
	struct {
		int i:3;
	};
};""",

"""struct T {
	long x:3;
	struct {
		long l:3;
	};
};""",

"""struct T {
	long long x:3;
	struct {
		long long ll:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		char c:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		short s:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		int i:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		long l:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		long long ll:3;
	};
};""",

"""struct T {
	char x:3;
	union {
		char c:3;
	};
};""",

"""struct T {
	short x:3;
	union {
		short s:3;
	};
};""",

"""struct T {
	int x:3;
	union {
		int i:3;
	};
};""",

"""struct T {
	long x:3;
	union {
		long l:3;
	};
};""",

"""struct T {
	long long x:3;
	union {
		long long ll:3;
	};
};""",

"""struct T {
	char x:3;
	_Alignas(sizeof(char)*2) struct {
		char c:3;
	};
};""",

"""struct T {
	short x:3;
	_Alignas(sizeof(short)*2) struct {
		short s:3;
	};
};""",

"""struct T {
	int x:3;
	_Alignas(sizeof(int)*2) struct {
		int i:3;
	};
};""",

"""struct T {
	long x:3;
	_Alignas(sizeof(long)*2) struct {
		long l:3;
	};
};""",

"""struct T {
	long long x:3;
	_Alignas(sizeof(long long)*2) struct {
		long long ll:3;
	};
};""",

"""struct T {
	char x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	};
};""",

"""struct T {
	short x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	};
};""",

"""struct T {
	int x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	};
};""",

"""struct T {
	long x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	};
};""",

"""struct T {
	long long x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	};
};""",

"""struct T {
	char x:3;
	struct {
		char c:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		short s:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		int i:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		long l:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		long long ll:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		char c:3;
	} __attribute__((packed));
};""",

"""struct T {
	short x:3;
	struct {
		short s:3;
	} __attribute__((packed));
};""",

"""struct T {
	int x:3;
	struct {
		int i:3;
	} __attribute__((packed));
};""",

"""struct T {
	long x:3;
	struct {
		long l:3;
	} __attribute__((packed));
};""",

"""struct T {
	long long x:3;
	struct {
		long long ll:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		char c:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		short s:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		int i:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		long l:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		long long ll:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	union {
		char c:3;
	} __attribute__((packed));
};""",

"""struct T {
	short x:3;
	union {
		short s:3;
	} __attribute__((packed));
};""",

"""struct T {
	int x:3;
	union {
		int i:3;
	} __attribute__((packed));
};""",

"""struct T {
	long x:3;
	union {
		long l:3;
	} __attribute__((packed));
};""",

"""struct T {
	long long x:3;
	union {
		long long ll:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	_Alignas(sizeof(char)*2) struct {
		char c:3;
	} __attribute__((packed));
};""",

"""struct T {
	short x:3;
	_Alignas(sizeof(short)*2) struct {
		short s:3;
	} __attribute__((packed));
};""",

"""struct T {
	int x:3;
	_Alignas(sizeof(int)*2) struct {
		int i:3;
	} __attribute__((packed));
};""",

"""struct T {
	long x:3;
	_Alignas(sizeof(long)*2) struct {
		long l:3;
	} __attribute__((packed));
};""",

"""struct T {
	long long x:3;
	_Alignas(sizeof(long long)*2) struct {
		long long ll:3;
	} __attribute__((packed));
};""",

"""struct T {
	char x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	} __attribute__((packed));
};""",

"""struct T {
	short x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	} __attribute__((packed));
};""",

"""struct T {
	int x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	} __attribute__((packed));
};""",

"""struct T {
	long x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	} __attribute__((packed));
};""",

"""struct T {
	long long x:3;
	struct {
		_Alignas(sizeof(int)*2) char c;
	} __attribute__((packed));
};""",

"""struct T {
	enum {
		VAL1,
		VAL2,
		VAL3,
	} e;
};""",

"""struct T {
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} e;
};""",

"""struct T {
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} e;
};""",

"""struct T {
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} e;
};""",

"""struct T {
	char x;
	enum {
		VAL1,
		VAL2,
		VAL3,
	} e;
};""",

"""struct T {
	char x;
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} e;
};""",

"""struct T {
	char x;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} e;
};""",

"""struct T {
	char x;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1,
		VAL2,
		VAL3,
	} e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} e;
};""",

"""struct T {
	enum {
		VAL1,
		VAL2,
		VAL3,
	} __attribute__((packed)) e;
};""",

"""struct T {
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __attribute__((packed)) e;
};""",

"""struct T {
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __attribute__((packed)) e;
};""",

"""struct T {
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x;
	enum {
		VAL1,
		VAL2,
		VAL3,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x;
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1,
		VAL2,
		VAL3,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1,
		VAL2 = VAL1,
		VAL3 = 1 << 10,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 1,
	} __attribute__((packed)) e;
};""",

"""struct T {
	char x:3;
	enum {
		VAL1 = -1,
		VAL2 = 0,
		VAL3 = 0xffffffff,
	} __attribute__((packed)) e;
};""",
]


text:
using Test
using Example2018

# Polynomial form a ring
#1. addition is associative
#2. additive identity
#3. additive inverse
#4. addition commutes
#5. multiplication is ...

@testset "Polynomials form a ring" begin
    z = Polynomial([0])
    p = Polynomial([1,2,3])
    for n in 1:100
        deg = rand(0:100)
        p = Polynomial(rand(-100000:100000, deg+1))

        @test p + z == p
        @test z + p == p
        @test p + (-p) == z
   end

end


text:
@info("loading libraries...")
using NBodyIPs, StaticArrays, BenchmarkTools, JuLIP, Test

# using JuLIP.Potentials: evaluate, evaluate_d
# using NBodyIPs: BondLengthDesc, BondAngleDesc, invariants, invariants_d, descriptor,
#                bo2angles, bo2edges
using NBodyIPs.Polys
using NBodyIPs.EnvIPs
using NBodyIPs.EnvIPs: EnvPoly

TRANSFORM = "r -> (2.9/r)^3"
DEGREES = [14, 12, 10]
ENVDEGS = [3, 2, 1]
RCUT = [7.0, 5.80, 4.5]

function random_ip(DT)
   DD = [ DT(TRANSFORM, (:cos, (0.66*rcut), rcut))  for rcut in RCUT ]

   rn = 1.5 * rnn(:W)
   Vn = ("1.0/(1.0 + exp(1.0 / ($rn - r + 1e-2)))", rn)

   # 2-body potential
   BB = [ envpolys(n+1, DD[n], DEGREES[n], Vn, ENVDEGS[n]) for n = 1:3 ]
   B = vcat(BB...)
   c = rand(length(B))
   IP = NBodyIP(B, c)

   return IP, fast(IP)
end

##

@info("Generate random EnvIP...")
IP, IPf = random_ip(BondLengthDesc)

@info("Generate faster EnvIP => EnvPoly ...")
V2 = IPf.components[1:4]
V3 = IPf.components[5:7]
V4 = IPf.components[8:9]

IPff = NBodyIP( [ EnvPoly([v.Vr for v in V2], V2[1].Vn, V2[1].str_Vn),
                  EnvPoly([v.Vr for v in V3], V3[1].Vn, V3[1].str_Vn),
                  EnvPoly([v.Vr for v in V4], V4[1].Vn, V4[1].str_Vn) ] )

##

at = rattle!(bulk(:W, cubic=true) * 6, 0.02)
forces(IPff, at);
forces(IPff, at);
# Profile.clear(); @profile forces(IPff, at);
# Profile.print()

# @code_warntype forces(IPff, at)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
export entropy

import StatsBase

@inline function StatsBase.entropy(m::ConditionalModel, N::Int=DEFAULT_SAMPLE_SIZE)
    return _entropy(getmoduletypencoding(m.model), m.model, argvals(m), Val(N))
end

sourceEntropy(m::Model, N::Int=DEFAULT_SAMPLE_SIZE) = sourceEntropy()(m, Val(N))

export sourceEntropy

function sourceEntropy()

    function(_m::Model, ::Val{_N}) where {_N}
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))

        function proc(_m, st::Sample)
            if isleaf(_m, st.x)
                return quote
                    _H += mean(entropy($(st.rhs)))
                end
            else
                return quote
                    $(st.x) = parts($(st.rhs), $_N)
                    _H += mean(entropy($(st.rhs)))
                end
            end

        end

        proc(_m, st::Return)  = nothing # :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),variables(_m))

        wrap(kernel) = @q begin
            _H = 0.0
            $kernel
            return _H
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

# TODO: Get entropy going again
# StatsBase.entropy(d::iid) = prod(d.size) * entropy(d.dist)
# StatsBase.entropy(d::For) = sum(entropy ∘ d.f, d.θ)

@gg function _entropy(M::Type{<:TypeLevel}, _m::Model, _args, _n::Val{_N}) where {_N}
    body = sourceEntropy()(type2model(_m), _n()) |> loadvals(_args, NamedTuple())
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end

@gg function _entropy(M::Type{<:TypeLevel}, _m::Model, _args::NamedTuple{()}, _n::Val{_N}) where {_N}
    body = sourceEntropy()(type2model(_m), _n)
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
# TODO: Document this.
macro order_by(src, exprs...)
    for expr in exprs
        # TODO: Add proper error messages here.
        @assert !has_alias(expr)
    end

    Expr(
        :call,
        :OrderBy,
        esc(src),
        Expr(
            :tuple,
            expression_expr.(exprs)...,
        ),
    )
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule JasPer_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("JasPer")
JLLWrappers.@generate_main_file("JasPer", UUID("f6ca3df5-1ae7-521a-b96b-cb83f57f3f1f"))
end  # module JasPer_jll


text:
using PyPlot

"""
    dotplot(p, color, marker)

Plot dots with `color` and `marker`
"""
function dotplot(p, color, marker)
    if marker == "o" || marker == "s"
        plot(p[1,:], p[2,:], color=color, linestyle="None", marker=marker, markerfacecolor="none", markeredgewidth=0.5)
    else
        plot(p[1,:], p[2,:], color=color, linestyle="None", marker=marker, ms=2)
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
"""
    EGPpowermix(σ, ξ, κ₁, κ₂, p)

*EGPpowermix* corresponds to the second extended GP model of Naveau et al. (2016), which is a mixture of power laws with G(v) = pv^κ₁ + (1-p)v^κ₂.

It is a five parameters family: κ₁ controls the lower tail behavior, κ₂ modifies the shape of the density in its central part, σ is a scale parameter, and ξ is a shape parameter.

```julia
EGPpowermix(σ, ξ, κ₁, κ₂, p)   # EGP of Naveau et al. (2016) (type 2) with parameters σ, ξ, κ₁, κ₂, p.

params(d)           # Get the parameters, i.e. (σ, ξ, κ₁, κ₂, p)
```

Reference :

* Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016), Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection, Water Resour. Res., 52, 2753– 2769, doi:10.1002/2015WR018552.
"""
struct EGPpowermix{T<:Real} <: ContinuousUnivariateDistribution
    σ::T    # scale parameter
    ξ::T    # rate of upper tail decay
    κ₁::T    # shape of the lower tail
    κ₂::T   # shape of the density in the central part
    p::T    # probability

    function EGPpowermix{T}(σ::T, ξ::T, κ₁::T,κ₂::T, p::T) where {T <: Real}
        new{T}(σ, ξ, κ₁, κ₂, p)
    end
end

function EGPpowermix(σ::T, ξ::T, κ₁::T, κ₂::T, p::T; check_args=true) where {T <: Real}
    check_args && @check_args(EGPpowermix, σ > zero(σ) && κ₁ > zero(κ₁) && κ₂ > zero(κ₂))  # manque test pour p
    return EGPpowermix{T}(σ, ξ, κ₁, κ₂, p)
end

EGPpowermix(σ::Real, ξ::Real, κ₁::Real, κ₂::Real, p::Real) = EGPpowermix(promote(σ, ξ, κ₁, κ₂, p)...)
EGPpowermix(σ::Integer, ξ::Integer, κ₁::Integer, κ₂::Integer, p::Integer) = EGPpowermix(float(σ), float(ξ), float(κ₁), float(κ₂), float(p))

minimum(d::EGPpowermix) = 0.0
maximum(d::EGPpowermix) = Inf * (d.ξ >= 0) - ( d.σ / d.ξ ) * (d.ξ < 0)
insupport(d::EGPpowermix, x::Real) = minimum(d) <= x <= maximum(d)

#### Parameters

scale(d::EGPpowermix) = d.σ
decay(d::EGPpowermix) = d.ξ    # noms à vérifier...
shape_tail(d::EGPpowermix) = d.κ₁
shape_central(d::EGPpowermix) = d.κ₂
prob(d::EGPpowermix) = d.p

params(d::EGPpowermix) = (d.σ, d.ξ, d.κ₁, d.κ₂, d.p)
partype(::EGPpowermix{T}) where {T} = T

#### Evaluation

function logpdf(d::EGPpowermix{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, κ₁, κ₂, p) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    g(v::Real, p::Real, κ₁::Real, κ₂::Real) = p*κ₁*v^(κ₁-1) + (1-p)*κ₂*v^(κ₂-1)
    lg(v::Real, p::Real, κ₁::Real, κ₂::Real) = log(g(v,p,κ₁,κ₂))

    lf = -log(σ) + Distributions.logpdf(pd, x/σ) + lg(Distributions.cdf(pd, x/σ),p,κ₁,κ₂)

    return lf
end

pdf(d::EGPpowermix, x::Real) = exp(logpdf(d, x))

function logcdf(d::EGPpowermix{T}, x::Real) where T<:Real
    μ = 0
    (σ, ξ, κ₁, κ₂, p) = params(d)

    pd = GeneralizedPareto(μ, 1, ξ)

    G(v::Real, p::Real, κ₁::Real, κ₂::Real) = p*v^κ₁ + (1-p)*v^κ₂

    lf = log(G(Distributions.cdf(pd, x/σ), p, κ₁, κ₂))

    return lf
end

cdf(d::EGPpowermix, x::Real) = exp(logcdf(d, x))

function quantile(d::EGPpowermix{T}, p::Real) where T<:Real
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    (σ, ξ, κ₁, κ₂, m) = params(d)

    G(v::Real, m::Real, κ₁::Real, κ₂::Real) = m*v^κ₁ + (1-m)*v^κ₂

    function invG(p, m, κ₁, κ₂)
        a = 0
        b = 1
        δ = 0.00001
        crit = 1
        X = NaN
        while crit > 2*δ
            X = (a + b)/2
            F = G(X, m, κ₁, κ₂)
            if F <= p
                a = X
            else
                b = X
            end
            crit = b - a
        end
        return X
    end

    x = (σ/ξ)*(((1 - invG(p,m,κ₁, κ₂))^(-ξ)) - 1)

    return x
end

#### Sampling

function rand(rng::AbstractRNG, d::EGPpowermix)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    G(v::Real, m::Real, κ₁::Real, κ₂::Real) = m*v^κ₁ + (1-m)*v^κ₂

    function invG(p, m, κ₁, κ₂)
        a = 0
        b = 1
        δ = 0.00001
        crit = 1
        X = NaN
        while crit > 2*δ
            X = (a + b)/2
            F = G(X, m, κ₁, κ₂)
            if F <= p
                a = X
            else
                b = X
            end
            crit = b - a
        end
        return X
    end

    return (d.σ/d.ξ)*(((1 - invG(u,d.p,d.κ₁,d.κ₂))^(-d.ξ)) - 1)
end

sampler(d::EGPpowermix) = d


text:
function mmek(;
  S=301, E=120, SE=0, P=0,
  k1=0.00166, k2=0.0001, k3=0.1)
  # Initialize network
  m = Network("Michaelis-Menten")

  # Species Definitions
  m <= Species("S",  S)
  m <= Species("E",  E)
  m <= Species("SE", SE)
  m <= Species("P",  P)

  # Reaction Definitions
  m <= Reaction("binding",      k1, "S + E --> SE")
  m <= Reaction("dissociation", k2, "SE --> S + E")
  m <= Reaction("conversion",   k3, "SE --> E + P")

  return m
end

model = mmek()


text:
using OpenCV
using Test
using Pkg.Artifacts
import LibGit2

if "OPENCV_TEST_DATA_PATH" in keys(ENV)
    test_dir = joinpath(ENV["OPENCV_TEST_DATA_PATH"], "cv")

else
    print("Skipping tests as test data not found.")
    exit(0)

    # artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
    # testdata_hash = artifact_hash("OPENCV_TEST_DATA_PATH", artifact_toml)

    # if testdata_hash === nothing || !artifact_exists(testdata_hash)
    #     testdata_hash = create_artifact() do artifact_dir
    #         opencv_extra_zip = "https://github.com/opencv/opencv_extra/archive/refs/heads/master.zip"
    #         download(opencv_extra_zip, artifact_dir)

    #         LibGit2.clone(opencv_extra_git, artifact_dir)
    #     end
    #     bind_artifact!(artifact_toml, "OPENCV_TEST_DATA_PATH", testdata_hash)
    # end

    # test_dir = joinpath(artifact_path(testdata_hash), "opencv_extra", "cv")

end

include("test_mat.jl")
include("test_feature2d.jl")
include("test_imgproc.jl")
include("test_objdetect.jl")
include("test_dnn.jl")


text:
using SparseArrays
using LinearAlgebra
using Plots

include("vekit.jl")
include("invvekit.jl")

let
    n=100
    # Matrix definieren
    B=Matrix(n*n*spdiagm(-1 => -ones(n-2),0 => 2*ones(n-1), 1 => -ones(n-2)))
    # y erstellen
    y=[]
    for k in 1:n-1
        y=[y;k]
    end
    m=0
    e=10^-3
    # größter Eigenvektor
    l_1,y_1=vekit(B,y,e)
    # kleinster Eigenvektor
    l_n1,y_n1=invvekit(B,y,m,e) #TODO: Tut nicht, wegen Fehler in invvekit

    #TODO: Plot
end

max_stars_repo_name: SoulKathi/vekit

max_issues_repo_name: SoulKathi/vekit

max_forks_repo_name: SoulKathi/vekit


text:
macro bin_str(x)
    a = BitVector((i for i in x) .== '1')
    return a
end

tobig(x) = parse(BigInt, bytes2hex(reverse(x)), base=16)

function int2bytes(x::Integer)
    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end

    return reverse(hex2bytes(hex))
end

struct StaticBigInt{N} <: Integer
    x::NTuple{N, UInt8}
end

function StaticBigInt(x::Integer; n = bitlength(x))
    bytes = int2bytes(x)

    N = div(n, 8, RoundUp)

    append!(bytes, UInt8[0 for i in 1:N-length(bytes)])

    sv = NTuple{N, UInt8}(bytes)
    return StaticBigInt(sv)
end

Base.convert(::Type{BigInt}, x::StaticBigInt) = tobig(x.x)
Base.convert(::Type{Integer}, x::StaticBigInt) = convert(BigInt, x)

Base.BigInt(x::StaticBigInt) = tobig(x.x)

style(x, n) = "\33[1;$(n)m$x\33[0m"

function Base.show(io::IO, a::StaticBigInt)
    show(io, BigInt(a))
end

Base.display(x::StaticBigInt) = show(x)

#######################

function _hex2bytes(x::String)

    normalized = join(split(x, " "), "")

    N = length(normalized)

    if mod(N, 2) != 0
        normalized = "0" * normalized
    end

    return hex2bytes(normalized)
end

macro hex_str(x)
    return _hex2bytes(x)
end

function bytes2bits(x::Vector{UInt8})
    bv = BitVector(u << -i % Bool for u in x for i in 7:-1:0)
    return bv
end

hex2bits(x::String) = bytes2bits(_hex2bytes(x))

tobits(x::String) = hex2bits(x)
tobits(x::Vector{UInt8}) = bytes2bits(x)
tobits(x) = convert(BitVector, x)

nbits(x::Integer) = sizeof(x) * 8
nbits(x::BigInt) = x.size * 64

function tobits(a::Integer) # I could still associate keyword arguments
    s = BitVector(a << -i % Bool for i in 0:(nbits(a)-1))
    n = findlast(x->x==1, s)
    if n == nothing
        n = 0
    end
    return reverse(s[1:n])
end

function tobits(a::Integer, n::Integer)
    s = BitVector(a << -i % Bool for i in 0:(nbits(a)-1))

    if length(s) < n
        resize!(s, n)
        return reverse(s)
    else
        return reverse(s[1:n])
    end
end

#####################

struct StaticBitVector{N}
    len::Int
    chunks::NTuple{N, UInt64}
end

StaticBitVector(x::BitVector) = StaticBitVector(x.len, Tuple(x.chunks))

function Base.convert(::Type{BitVector}, x::StaticBitVector)

    a = BitVector(undef, x.len)
    a.chunks = UInt64[x.chunks...]

    return a
end

Base.show(io::IO, x::StaticBitVector) = print(io, join(i ? "1" : "0" for i in convert(BitVector, x)))

modinv(s, q) = mod(gcdx(s, q)[2], q)

Base.length(x::StaticBitVector) = x.len

########################

static(x::BitVector) = StaticBitVector(x)
static(x::BigInt) = StaticBigInt(x)
static(x::Integer) = x # BigInt is the only exception
static(::Nothing) = nothing

struct StaticSymbol
    x::UInt128

    function StaticSymbol(s::Symbol)
        x = string(s)
        bytes = Vector{UInt8}(x)
        padded_bytes = append!(bytes, UInt8[0 for i in 1:16-length(bytes)])
        uint = reinterpret(UInt128, padded_bytes)
        return new(uint[1])
    end
end

function Base.convert(::Type{Symbol}, s::StaticSymbol)

    x = s.x

    bytes = reinterpret(UInt8, [x])

    N = findlast(x->x!=0, bytes)

    if N == nothing
        error("Nothing encoded")
    else
        return Symbol(String(bytes[1:N]))
    end
end

static(x::Symbol) = StaticSymbol(x)

# Goeas through the keyword arguments and returns a coresponding namedtuple where on each argument static is being called
static(; kwargs...) = NamedTuple((key, static(value)) for (key, value) in pairs(kwargs))


text:
module TestInterface
    using Test
    using CBSOData3

    include("mock_getodata.jl")

    tbl = CBSOData3.ODataTable("80137ned")
    @test length(tbl) == 675
end


text:
export FourierPreprocessor, PolynomialPreprocessor, TilingPreprocessor

"""
    FourierPreprocessor(order::Int)
Transform a scalar to a vector of `order+1` Fourier bases.
"""
struct FourierPreprocessor <: AbstractPreprocessor
    order::Int
end

(p::FourierPreprocessor)(s::Number) = [cos(i * π * s) for i = 0:p.order]

"""
    PolynomialPreprocessor(order::Int)
Transform a scalar to vector of maximum `order` polynomial.
"""
struct PolynomialPreprocessor <: AbstractPreprocessor
    order::Int
end

(p::PolynomialPreprocessor)(s::Number) = [s^i for i = 0:p.order]

"""
    TilingPreprocessor(tilings::Vector{<:Tiling})
Use each `tilings` to encode the state and return a vector.
"""
struct TilingPreprocessor{Tt<:Tiling} <: AbstractPreprocessor
    tilings::Vector{Tt}
end

(p::TilingPreprocessor)(s::Union{<:Number,<:Array}) = [encode(t, s) for t in p.tilings]


text:
using Test
using NeXLCore

@testset "MatUncertain" begin
    mat = "Au60Ag40"
    inplbls = [
        MassFractionLabel(mat, n"Ag"),
        AtomicWeightLabel(mat, n"Ag"),
        MassFractionLabel(mat, n"Au"),
        AtomicWeightLabel(mat, n"Au"),
    ]
    inpvals = Float64[0.4020, 107.87, 0.5950, 196.97]
    inpcovs = Float64[(0.0090)^2 0 0 0; 0 0 0 0; 0 0 (0.0120)^2 0; 0 0 0 0]
    inputs = uvs(inplbls, inpvals, inpcovs)

    @testset "AuAg - MFtoAF" begin

        afs = propagate(MFtoAF(mat), inputs)

        @test isapprox(value(afs[AtomicFractionLabel(mat, n"Ag")]), 0.5523, atol = 0.0001)
        @test isapprox(value(afs[AtomicFractionLabel(mat, n"Au")]), 0.4477, atol = 0.0001)

        @test isapprox(σ(afs[AtomicFractionLabel(mat, n"Ag")]), 0.0075, atol = 0.0001)
        @test isapprox(σ(afs[AtomicFractionLabel(mat, n"Au")]), 0.0075, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(
                afs,
                AtomicFractionLabel(mat, n"Ag"),
                AtomicFractionLabel(mat, n"Au"),

            ),
            -1.0,
            atol = 0.0001,
        )
    end

    @testset "AuAg - MFtoNMF" begin

        afs = propagate(MFtoNMF(mat), inputs)

        @test isapprox(value(afs[NormMassFractionLabel(mat, n"Ag")]), 0.4032, atol = 0.0001)
        @test isapprox(value(afs[NormMassFractionLabel(mat, n"Au")]), 0.5968, atol = 0.0001)

        @test isapprox(σ(afs[NormMassFractionLabel(mat, n"Ag")]), 0.0073, atol = 0.0001)
        @test isapprox(σ(afs[NormMassFractionLabel(mat, n"Au")]), 0.0073, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(
                afs,
                NormMassFractionLabel(mat, n"Ag"),
                NormMassFractionLabel(mat, n"Au"),
            ),
            -1.0,
            atol = 0.0001,
        )
    end

    @testset "AuAg - MatStats" begin

        afs = propagate(MatStats(mat), inputs)

        @test isapprox(value(afs[MeanZ(mat)]), 65.8990, atol = 0.0001)
        @test isapprox(value(afs[MeanAz(mat)]), 160.56, atol = 0.01)

        @test isapprox(σ(afs[MeanZ(mat)]), 1.0381, atol = 0.0001)
        @test isapprox(σ(afs[MeanAz(mat)]), 2.5552, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(afs, MeanZ(mat), MeanAz(mat)),
            1.0,
            atol = 0.001,
        )
    end

    @testset "AuAg - Comp" begin
        afs = mf2comp(mat, inputs)

        @test isapprox(value(afs[AtomicFractionLabel(mat, n"Ag")]), 0.5523, atol = 0.0001)
        @test isapprox(value(afs[AtomicFractionLabel(mat, n"Au")]), 0.4477, atol = 0.0001)

        @test isapprox(σ(afs[AtomicFractionLabel(mat, n"Ag")]), 0.0075, atol = 0.0001)
        @test isapprox(σ(afs[AtomicFractionLabel(mat, n"Au")]), 0.0075, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(
                afs,
                AtomicFractionLabel(mat, n"Ag"),
                AtomicFractionLabel(mat, n"Au")
            ),
            -1.0,
            atol = 0.0001,
        )

        @test isapprox(value(afs[NormMassFractionLabel(mat, n"Ag")]), 0.4032, atol = 0.0001)
        @test isapprox(value(afs[NormMassFractionLabel(mat, n"Au")]), 0.5968, atol = 0.0001)

        @test isapprox(σ(afs[NormMassFractionLabel(mat, n"Ag")]), 0.0073, atol = 0.0001)
        @test isapprox(σ(afs[NormMassFractionLabel(mat, n"Au")]), 0.0073, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(
                afs,
                NormMassFractionLabel(mat, n"Ag"),
                NormMassFractionLabel(mat, n"Au")
            ),
            -1.0,
            atol = 0.0001,
        )

        @test isapprox(value(afs[MeanZ(mat)]), 65.8990, atol = 0.0001)
        @test isapprox(value(afs[MeanAz(mat)]), 160.56, atol = 0.01)

        @test isapprox(σ(afs[MeanZ(mat)]), 1.0381, atol = 0.0001)
        @test isapprox(σ(afs[MeanAz(mat)]), 2.5552, atol = 0.0001)
        @test isapprox(
            NeXLUncertainties.correlation(afs, MeanZ(mat), MeanAz(mat)),
            1.0,
            atol = 0.001,
        )
    end

    @testset "AuAg - MacU" begin

        mat, cxr = material("Au60Ag40", n"Ag" => 0.4020, n"Au" => 0.5950), n"O K-L3"
        inplbls = [
            MassFractionLabel(mat.name, n"Ag"),
            AtomicWeightLabel(mat.name, n"Ag"),
            MassFractionLabel(mat.name, n"Au"),
            AtomicWeightLabel(mat.name, n"Au"),
            μoρElementLabel(n"Ag", cxr),
            μoρElementLabel(n"Au", cxr),
        ]
        macAg, macAu = macU(n"Ag", cxr), macU(n"Au", cxr)
        @test isapprox(value(macAg), mac(n"Ag", cxr), rtol = 1.0e-6)
        @test isapprox(value(macAu), mac(n"Au", cxr), rtol = 1.0e-6)

        inpvals = Float64[0.4020, 107.87, 0.5950, 196.97, value(macAg), value(macAu)]
        inpcovs = Float64[
            (0.0090)^2 0 0 0 0 0
            0 0 0 0 0 0
            0 0 (0.0120)^2 0 0 0
            0 0 0 0 0 0
            0 0 0 0 variance(macAg) 0
            0 0 0 0 0 variance(macAu)
        ]

        inputs = uvs(inplbls, inpvals, inpcovs)
        model = μoρMaterial(mat.name, cxr)
        res = model(inputs)
        mc_res = mcpropagate(
            model,
            inputs,
            1000,
            parallel = false,
            rng = MersenneTwister(0xBADF00D),
        )
        #println(res)
        #println(mc_res)
        lbl = μoρLabel(mat.name, cxr)
        @test isapprox(value(res[lbl]), mac(mat, cxr), rtol = 1.0e-6)
        @test isapprox(value(res[lbl]), value(mc_res[lbl]), rtol = 0.01)
        @test isapprox(σ(res[lbl]), σ(mc_res[lbl]), rtol = 0.01)
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

t_start = 0.0
t_stop = 8.0

# case 1: ME-FMU with state events

pathToFMU = joinpath(dirname(@__FILE__), "..", "model", ENV["EXPORTINGTOOL"], "SpringFrictionPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateME(fmuStruct, t_start, t_stop)
@test length(solution.u) > 0
@test length(solution.t) > 0

@test solution.t[1] == t_start
@test solution.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test solution.u[1] == [0.5, 0.0]
@test sum(abs.(solution.u[end] - [1.06736, -1.03552e-10])) < 0.01
fmiUnload(myFMU)

# case 2: ME-FMU with state and time events

if ENV["EXPORTINGTOOL"] == "Dymola/2020x"
    pathToFMU = joinpath(dirname(@__FILE__), "..", "model", ENV["EXPORTINGTOOL"], "SpringTimeFrictionPendulum1D.fmu")

    myFMU = fmiLoad(pathToFMU)

    comp = fmiInstantiate!(myFMU; loggingOn=false)
    @test comp != 0

    # choose FMU or FMUComponent
    fmuStruct = nothing
    envFMUSTRUCT = ENV["FMUSTRUCT"]
    if envFMUSTRUCT == "FMU"
        fmuStruct = myFMU
    elseif envFMUSTRUCT == "FMUCOMPONENT"
        fmuStruct = comp
    end
    @assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

    solution = fmiSimulateME(fmuStruct, t_start, t_stop)
    @test length(solution.u) > 0
    @test length(solution.t) > 0

    @test solution.t[1] == t_start
    @test solution.t[end] == t_stop

    # reference values from Simulation in Dymola2020x (Dassl)
    @test solution.u[1] == [0.5, 0.0]
    @test sum(abs.(solution.u[end] - [1.05444, 1e-10])) < 0.01
    fmiUnload(myFMU)
end


text:
using ProgressMeter
using Dates

function main()

    start_t = now()
    out_dir  = "result"
    PARAMDAT = "PARAMDAT.json"
    fwrite   = "write"

    nval = 4
    Rd   = 287.0
    R    = 8.314

    xmax, ymax, nodes, vecAx, vecAy = read_allgrid()
    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau, cfl,
    init_rho, init_u, init_v, init_p, init_T, specific_heat_ratio, Rd, bdcon = input_para(PARAMDAT)

    init_p = init_rho*Rd*init_T

    Qbase, cellxmax, cellymax, restartnum = set_initQbase(xmax, ymax, restart_file, init_rho, init_u, init_v, init_p, init_T,
                                                        specific_heat_ratio, out_file_front, out_ext, out_dir, restartnum, Rd, nval)

    # init Delta_Qcon_hat
    volume = set_volume(nodes, cellxmax, cellymax)
    dx, dy = set_dx_lts(nodes, cellxmax, cellymax)
    reset_write(fwrite)

    # main loop
    print("threads num : ")
    println(Threads.nthreads())

    # check bd
    check_bd(bdcon)

    #throw(UndefVarError(:x))
    # main loop
    prog = Progress(nt,1)
    @time for t in 1:nt
        next!(prog)

        evalnum = t+restartnum
        if time_integ == "1"
            # exlicit scheme

            Qbase    = set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, specific_heat_ratio, nval)
            Qcon     = base_to_conservative(Qbase, cellxmax, cellymax, specific_heat_ratio, nval)
            Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval)

            # initial_setup
            mu     = set_mu(Qbase, cellxmax, cellymax, specific_heat_ratio, Rd)
            lambda = set_lambda(Qbase, cellxmax, cellymax, mu, specific_heat_ratio, Rd)

            # RHS
            # advection_term
            #E_adv_hat, F_adv_hat = AUSM(Qbase,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio)
            E_adv_hat, F_adv_hat = AUSM_plus(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)

            # viscos_term
            E_vis_hat, F_vis_hat = central_diff(Qbase, Qcon, cellxmax, cellymax, mu, lambda,
                                                vecAx, vecAy, specific_heat_ratio, volume, Rd, nval)

            RHS = setup_RHS(cellxmax, cellymax, E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, nval, volume)

            # time integral
            Qcon_hat = time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, nval)

            Qcon  = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval)
            Qbase = conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval)

        elseif time_integ == "2"
            output_physicaltime(fwrite, t, dt)
            # Qcn
            Qbasen    = copy(Qbase)
            Qconn     = base_to_conservative(Qbasen, cellxmax, cellymax, specific_heat_ratio, nval)
            Qconn_hat = setup_Qcon_hat(Qconn, cellxmax, cellymax, volume, nval)

            Qbasem = copy(Qbase)

            for tau in 1:in_nt
                # LHS (A_adv_hat=jacobian)
                Qbasem = set_boundary(Qbasem, cellxmax, cellymax, vecAx, vecAy, bdcon, Rd, specific_heat_ratio, nval)

                Qcon     = base_to_conservative(Qbasem, cellxmax, cellymax, specific_heat_ratio, nval)
                Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval)

                # initial_setup
                mu     = set_mu(Qbasem, cellxmax, cellymax, specific_heat_ratio, Rd)
                lambda = set_lambda(Qbasem, cellxmax, cellymax, mu, specific_heat_ratio, Rd)
                dtau   = set_lts(Qbasem, cellxmax, cellymax, mu, dx, dy, vecAx, vecAy, volume, specific_heat_ratio, cfl)

                # RHS
                #advection_term
                E_adv_hat, F_adv_hat = AUSM_plus(Qbasem, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)

                # viscos_term
                E_vis_hat, F_vis_hat = central_diff(Qbasem, Qcon, cellxmax, cellymax, mu, lambda,
                                                vecAx, vecAy, specific_heat_ratio, volume, Rd, nval)

                RHS = setup_RHS(cellxmax, cellymax, E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, nval, volume)

                # lusgs_advection_term
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_shig, B_beta_shig = one_wave(Qbasem, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)
                # lusgs_viscos_term
                jalphaP, jbetaP = central_diff_jacobian(Qbasem, Qcon, cellxmax, cellymax, mu, lambda, vecAx, vecAy, specific_heat_ratio, volume, nval)

                # LUSGS
                delta_Q = zeros(cellxmax, cellymax, nval)

                ite = 0
                norm2 = zeros(nval)
                while true
                    delta_Q_temp = copy(delta_Q)

                    delta_Q = lusgs(dt, dtau, Qcon_hat, Qconn_hat, delta_Q, A_adv_hat_p,  A_adv_hat_m,  B_adv_hat_p,  B_adv_hat_m,  A_beta_shig,  B_beta_shig, jalphaP,  jbetaP, RHS, cellxmax, cellymax, volume, nval)

                    res   = set_res(delta_Q, delta_Q_temp, cellxmax, cellymax, nval)
                    norm2 = check_converge(res, RHS, cellxmax, cellymax, init_small, nval)

                    ite += 1
                    if norm2[1] < norm_ok && norm2[2] < norm_ok && norm2[3] < norm_ok && norm2[4] < norm_ok
                        break
                    end
                    if ite % 100 ==0
                        println(" now cal norm2 ")
                        println(norm2)
                    end
                end
                output_innertime(fwrite, tau, norm2, nval)

                for l in 1:nval
                    for j in 2:cellymax-1
                        for i in 2:cellxmax-1
                            Qcon_hat[i,j,l] = Qcon_hat[i,j,l] + delta_Q[i,j,l]
                        end
                    end
                end

                Qcon = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval)
                Qbasem = conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval)

            end
            Qbase = copy(Qbasem)
        end

        if round(evalnum) % every_outnum == 0
            println("\n")
            println("nt_______________________________"*string(round(evalnum)))
            output_result(evalnum, Qbase, cellxmax, cellymax, specific_heat_ratio, out_file_front, out_ext, out_dir, Rd, nval)
        end

        check_divrege(Qbase, cellxmax, cellymax, Rd, fwrite)
    end

    end_t = now()
    output_fin(fwrite, start_t, end_t, nt, dt, in_nt, cellxmax, cellymax)
end

# -- main --
main()


text:
@contract begin
    require(a > 0, b > a, c)
    ensure(b < a, a < 0, !c)
    function foo(a::Int64, b::Int64, c::Bool)
        a *= -1
        b *= -1
        c = !c
        return
    end
end

@contract begin
    require(a > 0, b > a, c)
    ensure(b < a, a > b, !c)
    function flawedEnsureFoo(a::Int64, b::Int64, c::Bool)
        a -= 2
        b += 2
        c = true
        return
    end
end

@testset "Assertion Errors" begin
    @test_throws ContractBreachException foo(-1, 10, true)
    @test_throws ContractBreachException foo(5, 2, true)
    @test_throws ContractBreachException foo(5, 10, false)

    @test_throws ContractBreachException flawedEnsureFoo(10, 100, true)
    @test_throws ContractBreachException flawedEnsureFoo(5, 100, true)
    @test_throws ContractBreachException flawedEnsureFoo(2, 100, true)
end

@testset "Contract Breach Description" begin
    try
        foo(-1, 10, true)
    catch e
        b = IOBuffer()
        showerror(b, e)
        @test String(take!(b)) ==
              "ContractBreachException: Breach on Requirement Expression 'a > 0' in function 'foo'"
    end

    try
        flawedEnsureFoo(10, 100, true)
    catch e
        b = IOBuffer()
        showerror(b, e)
        @test String(take!(b)) ==
              "ContractBreachException: Breach on Ensure Expression 'b < a' in function 'flawedEnsureFoo'"
    end
end


text:
using Turing
using DelimitedFiles

using BenchmarkHelper

if !haskey(BenchmarkSuite, "nuts")
    BenchmarkSuite["nuts"] = BenchmarkGroup(["nuts"])
end

fname = joinpath(dirname(@__FILE__), "sv_nuts.data")
y, header = readdlm(fname, ',', header=true)

# Stochastic volatility (SV)
@model function sv_nuts(y, dy, ::Type{T}=Vector{Float64}) where {T}
    N = size(y,1)

    τ ~ Exponential(1/100)
    ν ~ Exponential(1/100)
    s = T(undef, N)

    s[1] ~ Exponential(1/100)
    for n in 2:N
        s[n] ~ Normal(log(s[n-1]), τ)
        s[n] = exp(s[n])
        dy = log(y[n] / y[n-1]) / s[n]
        dy ~ TDist(ν)
    end
end

# Sampling parameter settings
n_samples = 10_000

# Sampling
BenchmarkSuite["nuts"]["sv"] = @benchmarkable sample(sv_nuts(y, NaN), NUTS(0.65), n_samples)


text:
using Plots

function finish()
    plot!(legend = false, aspect_ratio = 1, axis = false, grid = false)
end

gr()


text:
export TR_step_computation

function TR_step_computation(h  :: Any,
                             g  :: Any,
                             dN :: Any,
                             Δ  :: Any)
  if true in (h .> 0.0)
    if true in (g .> 0.0)
      d = max.(-Δ, dN)
    else
      d = min.(dN, Δ)
    end
  else
    if true in (g .> 0.0)
      d = .- Δ
    else
      d = Δ
    end
  end

  d = d[1]

  return [d]
end


text:
# Lagrangian for kuka iiwa via RigidBodyDynamics

using LinearAlgebra
using ForwardDiff
using RigidBodyDynamics
using MeshCatMechanisms

nq = 7 + 3 # number of states - kuka arm joint positions + particle in 3-space
nu = 7 # number of controls - kuka arm joint torques
nγ = 1 # contact points/ number of contact forces
ns = nγ

h = 0.05 # s

m_p = 0.01 # kg - mass of particle

urdf = "/home/taylor/Research/contact_implicit_trajectory_optimization/models/kuka/temp/kuka.urdf"
mechanismkuka = parse_urdf(urdf)
visualskuka = URDFVisuals(urdf)

function get_kuka_ee(kuka)
    ee_body = findbody(kuka, "iiwa_link_ee")
    ee_point = Point3D(default_frame(ee_body),0.,0.,0.)
    return ee_body, ee_point
end

function get_kuka_ee_position_fun(kuka::Mechanism,statecache=StateCache(kuka)) where {O}
    ee_body, ee_point = get_kuka_ee(kuka)
    world = root_frame(kuka)
    nn = num_positions(kuka)

    function ee_position(x::AbstractVector{T}) where T
        state = statecache[T]
        set_configuration!(state, x[1:nn])
        RigidBodyDynamics.transform(state, ee_point, world).v
    end
end

end_effector_function = get_kuka_ee_position_fun(parse_urdf(urdf,remove_fixed_tree_joints=false))

# ee_pos = end_effector_function(zeros(7))
# state = MechanismState(mechanism)

const statecachekuka = StateCache(mechanismkuka)

function 𝓛(q::AbstractVector{T},q̇::AbstractVector{T}) where T
    q_kuka = q[1:7]
    q̇_kuka = q̇[1:7]

    q_p = q[8:10]
    q̇_p = q̇[8:10]

    state = statecachekuka[T]

    set_configuration!(state,q_kuka)
    set_velocity!(state,q̇_kuka)

    kinetic_energy(state)  + 0.5*m_p*q̇_p'*q̇_p - m_p*9.81*q_p[3] - gravitational_potential_energy(state)
end

𝓛d(q1,q2) = let 𝓛=𝓛, h=h
    h*𝓛(0.5*(q1+q2),(q2-q1)/h)
end
𝓛d(z) = 𝓛d(z[1:nq],z[nq .+ (1:nq)])

D𝓛d(z) = ForwardDiff.gradient(𝓛d,z)
D1𝓛d(z) = D𝓛d(z)[1:nq]
D2𝓛d(z) = D𝓛d(z)[nq .+ (1:nq)]
D1𝓛d(q1,q2) = D1𝓛d([q1;q2])
D2𝓛d(q1,q2) = D2𝓛d([q1;q2])

δD1𝓛dδq2(z) = ForwardDiff.jacobian(D1𝓛d,z)[1:nq,nq .+ (1:nq)]
δD1𝓛dδq2(q1,q2) = δD1𝓛dδq2([q1;q2])

d = 0.5

ϕ(q) = let d=d
    q_kuka = q[1:7]
    ee = end_effector_function(q_kuka)
    q_p = q[8:10]
    # println(ee)
    # println(q_p)
    d - norm(ee - q_p)
end
∇ϕ(q) = ForwardDiff.gradient(ϕ,q)

B = zeros(nq,nu)
B[1:nu,1:nu] = Diagonal(ones(nu))
B

F(q,q̇) = let m=mass,g=9.81
    zero(q)
end

Fd(q1,q2) = let F=F, h=h
    h*F(0.5*(q1+q2), (q2-q1)/h)
end

var_int(q1,q2,u,γ,s,q⁺) = let B=B
    D2𝓛d(q1,q2) + D1𝓛d(q2,q⁺) + 0.5*(Fd(q1,q2) + Fd(q2,q⁺)) + B[:,:]*u[:] + ∇ϕ(q⁺)[:,:]*γ[:]
end
var_int(rand(nq),rand(nq),rand(nu),rand(nγ),rand(ns),rand(nq))

function var_int(z)
    q1 = z[1:nq]
    q2 = z[nq .+ (1:nq)]
    u = z[2nq .+ (1:nu)]
    γ = z[(2nq + nu) .+ (1:nγ)]
    s = z[(2nq + nu + nγ) .+ (1:ns)]
    q⁺ = z[(2nq + nu + nγ + ns) .+ (1:nq)]
    var_int(q1,q2,u,γ,s,q⁺)
end
∇var_int(z) = ForwardDiff.jacobian(var_int,z)
∇var_int(q1,q2,u,γ,s,q⁺) = ∇var_int([q1;q2;u;γ;s;q⁺])
∇var_int(rand(nq),rand(nq),rand(nu),rand(nγ),rand(ns),rand(nq))


text:
#==============================================================================#
# SNS/test/runtests.jl
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

using AWSCore
using AWSLambda
using AWSSNS
using Retry
using SymDict
using JSON
using Base.Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(
                         region = "ap-northeast-1",
                         lambda_bucket = "ocaws.jl.lambdatest.tokyo",
                         #region = "us-east-1",
                         #lambda_bucket = "ocaws.jl.lambdatest",
                         lambda_packages = ["Requests",
                                            "Nettle",
                                            "LightXML",
                                            "JSON",
                                            "DataStructures",
                                            "StatsBase",
                                            "DataFrames",
                                            "DSP",
                                            "GZip",
                                            "ZipFile",
                                            "IniFile",
                                            "SymDict",
                                            "XMLDict",
                                            "Retry"
                                           ])

#create_jl_lambda_base(aws)

#using AWSS3
#s3_copy(aws, "ocaws.jl.lambdatest", "jl_lambda_base.zip",
#             to_bucket="ocaws.jl.lambdatest.tokyo", to_path= "jl_lambda_base.zip")

#-------------------------------------------------------------------------------
# Lambda tests
#-------------------------------------------------------------------------------

# Count primes in the cloud...

λ = @λ aws function count_primes(low::Int, high::Int)
    count = length(primes(low, high))
    println("$count primes between $low and $high.")
    return count
end

@test invoke_lambda(aws, "count_primes", low = 10, high = 100)[:jl_data] == "21"

# Count primes in parallel...

function count_primes(low::Int, high::Int)
    w = 500000000
    counts = amap(λ, [(i, min(high,i + w)) for i = low:w:high])
    count = sum(counts)
    println("$count primes between $low and $high.")
    return count
end

@test count_primes(10, 10000000000) == 455052507

mktempdir() do tmp
    cd(tmp) do

        # Create a test module under "tmp"...

        mkpath("TestModule")
        open(io->write(io, """
            module TestModule

            export test_function

            __precompile__()

            test_function(x) = x * x

            end
        """), "TestModule/TestModule.jl", "w")

        push!(LOAD_PATH, "TestModule")

        # Create a lambda that uses the TestModule...
        λ = @λ aws function lambda_test(x)

            # Check that precompile cache is being used...
            @assert !Base.stale_cachefile("/var/task/TestModule/TestModule.jl",
                                          "/var/task/TestModule.ji")
            using TestModule
            return test_function(x)
        end

        @test λ(4) == 16
    end
end

#-------------------------------------------------------------------------------
# API Gateway tests
#-------------------------------------------------------------------------------

if false

for api in apigateway_restapis(aws)
    apigateway(aws, "DELETE", "/restapis/$(api["id"])")
end

apigateway_create(aws, "count_primes", (:low, :high))

end

#==============================================================================#
# End of file.
#==============================================================================#


text:
"Problem 5 in the Hock-Schittkowski suite"
function hs5()

  nlp = Model()

  l = [-1.5, -3]
  u = [4, 3]
  @variable(nlp, l[i] ≤ x[i=1:2] ≤ u[i], start=0.0)

  @NLobjective(
    nlp,
    Min,
    sin(x[1] + x[2]) + (x[1] - x[2])^2 - 1.5 * x[1] + 2.5 * x[2] + 1
  )

  return nlp
end


text:
# Low Rank Singular Value Decomposition
"""
    get_approximate_basis(A, l::Int; niter::Int = 2, M = nothing) -> Q
Return Matrix ``Q`` with ``l`` orthonormal columns such that ``Q Q^H A`` approximates ``A``. If ``M`` is specified, then ``Q`` is such that ``Q Q^H (A - M)`` approximates ``A - M``.
"""

struct Normal_QR{QT, RT}
    Q::QT
    R::RT
    function Normal_QR(Q::QT, R::RT) where {QT<:AbstractMatrix, RT<:AbstractMatrix}
        size(Q, 2) == size(R, 1) || throw(DimensionMismatch("$(size(Q)), $(length(R)) not compatible"))
        new{QT, RT}(Q, R)
    end
end

function private_qr(A::AbstractMatrix{T}) where T
	res = LinearAlgebra.qr(A)
	Normal_QR(Matrix(res.Q), res.R)
end

function get_approximate_basis(
    A::AbstractSparseMatrix{T}, l::Int, niter::Int = 2, M::Union{AbstractMatrix{T}, Nothing} = nothing) where T
    m, n = size(A)
    Ω = rand(T, (n, l))
    if M === nothing
        F_j = private_qr(A * Ω)
        for j = 1:niter
            F_H_j = private_qr(A' * Matrix(F_j.Q))
            F_j = qr(A * Matrix(F_H_j.Q))
        end
    else
        F_j = qr(A * Ω .- M * Ω)
        for j = 1:niter
            F_H_j = private_qr(A' * Matrix(F_j.Q) - M' * sum(Matrix(F_j.Q), dims=1))
            F_j = private_qr(A * Matrix(F_H_j.Q) .- M * Matrix(F_H_j.Q))
        end
    end
    Matrix(F_j.Q)
end

"""
    low_rank_svd(A, l::Int; niter::Int = 2, M = nothing) -> U, S, Vt
Return the singular value decomposition of a sparse matrix ``A`` with estimated rank ``l`` such that ``A ≈ U diag(S) Vt``. In case row vector ``M`` is given, then SVD is computed for the matrix ``A - M``.
"""

function low_rank_svd(A::AbstractSparseMatrix{T}, l::Int, niter::Int = 2, M::Union{AbstractMatrix{T}, Nothing} = nothing) where T

    if M === nothing
        Q = get_approximate_basis(A, l, niter, M)
        B = Q' * A
    else
        size(M, 2) == size(A, 2) || throw(DimensionMismatch())
        Q = get_approximate_basis(A, l, niter, M)
        B = Q' * A - sum(Q', dims=2) * M
    end

    dense_svd = svd(B)
    U = Q * dense_svd.U
    return U, dense_svd.S, dense_svd.Vt
end


text:
function append_depthfirst!(g, lookup, node)
    push!(lookup, node)
    pv = nv(g)
    for c in node["children"]
        add_vertex!(g)
        add_edge!(g, pv, nv(g))
        append_depthfirst!(g, lookup, c)
    end
    g, lookup
end

"""
    g, vertexlist = ontology(filename)
    g, vertexlist = ontology("mouse")

Return the structure graph `g` and list `vertexlist` of metadata for nodes.
"""
function ontology(filename)
    str = readstring(dataset(filename, "ontology"))
    hierarchy = JSON.parse(str)
    g = DiGraph()
    add_vertex!(g)  # for the root node
    vertexlist = Dict{String,Any}[]
    append_depthfirst!(g, vertexlist, hierarchy)
    g, vertexlist
end

"""
    matches = findvertices(pat, vertexlist, field="acronym")

Return the list of vertices that Regex-match to `pat` for the
specified field ( `"acronym"` or `"name"`) of each node listed in `vertexlist`.
"""
findvertices(pat, vertexlist, field="acronym") = find(x->ismatch(pat, x[field]), vertexlist)

"""
    ids = structureids(v, vertexlist)

Return the structure_id associated with each vertex in `v`.
"""
structureids(v, vertexlist) = map(i->vertexlist[i]["id"], v)
structureids(pat::Regex, vertexlist, field="acronym") =
    structureids(findvertices(pat, vertexlist, field), vertexlist)

function annotation(resolution=50; species="mouse")
    anno = load(dataset(species, "annotation_"*string(resolution)))
    # The NRRD header is incorrect, we should fix the orientation
    r = resolution*μm
    return AxisArray(anno.data, (:P, :I, :R), (r, r, r))
end

function download_annotation(resolution=50)
    dest = joinpath(mousedir, "annotation_$resolution.nrrd")
    download("http://download.alleninstitute.org/informatics-archive/current-release/mouse_ccf/annotation/ccf_2016/annotation_$resolution.nrrd", dest)
end

"""
    bbleft, bbright = boundingbox(f, annotated)

Return the BoundingBoxes encapsulating structure(s) for which
`f(annotated[i])` returns `true`, separately for the two hemispheres
of the brain. `annotated` is a labeled volume. Units of `bb` are
microns.
"""
function boundingbox(f, annotated)
    scalebb(bx, px) = BoundingBox(map((iv, s) -> minimum(iv)*s .. maximum(iv)*s, bx.intervals, ps))
    local bb
    isfirst = true
    # Search only the left half of the volume
    inds = mapfilter(i->first(i):(first(i)+last(i))÷2, indices(annotated),
                     Axis{:R}, axes(annotated))
    for i in CartesianRange(inds)
        if f(annotated[i])
            if isfirst
                bb = BoundingBox(i, i)
                isfirst = false
            else
                bb |= i
            end
        end
    end
    dim = find(x->x==:R, axisnames(annotated))[1]
    bbr = BoundingBox(ntuple(ndims(annotated)) do d
        iv = bb.intervals[d]
        if d == dim
            id = indices(annotated, d)
            l = first(id) + last(id)
            l-maximum(iv) .. l-minimum(iv)
        else
            iv
        end
    end)
    ps = pixelspacing(annotated)
    scalebb(bb, ps), scalebb(bbr, ps)
end


text:
@testset "sparse" begin
    @testset "conversion" begin
        nrow = 4
        ncol = 5
        nnz = 7
        mat = [0 1 4 0 0;
            0 0 0 0 0;
            0 2 0 6 0;
            0 3 5 7 0]

        csr_rowptr = [1, 3, 3, 5, 8]
        csr_col = [2, 3, 2, 4, 2, 3, 4]
        csr_nzv = [1, 4, 2, 6, 3, 5, 7]

        csc_colptr = [1, 1, 4, 6, 8, 8]
        csc_row = [1, 3, 4, 1, 4, 3, 4]
        csc_nzv = collect(1:nnz)

        nzval, colptr, row = Pickle.csr_to_csc(nrow, ncol, csr_nzv, csr_rowptr, csr_col)
        @test colptr == csc_colptr
        @test row == csc_row
        @test nzval == csc_nzv
    end

    @testset "load" begin
        spA = SparseMatrixCSC(Array(to_nparray(pyload("./test_pkl/csr_mat.pkl"))))
        @test Pickle.npyload("./test_pkl/csr_mat.pkl") == spA
    end
end


text:
@testset "Principal Value (MRP)" begin
    for i = 1:1000
        p = MRP(5.0 * randn(), 5.0 * randn(), 5.0 * randn())
        p_prin = principal_value(p)
        @test p_prin ≈ p
        @test (p_prin.x^2 + p_prin.y^2 + p_prin.z^2) < (1.0 + 1.0e-13)
        # end
    end
end

@testset "Principal Value (UnitQuaternion)" begin
    for i = 1:1000
        qq = rand(UnitQuaternion)
        qq_prin = principal_value(qq)
        @test 0.0 <= qq_prin.w
        @test qq_prin ≈ qq
    end
end

@testset "Principal Value (Angle Axis)" begin
    for i = 1:1000
        aa = AngleAxis(100.0 * randn(), randn(), randn(), randn())
        aa_prin = principal_value(aa)
        @test 0.0 <= aa_prin.theta
        @test aa_prin ≈ aa
    end
end

@testset "Principal Value (Rotation Vector)" begin
    for i = 1:1000
        rv = RotationVec(100.0 * randn(), 100.0 * randn(), 100.0 * randn())
        rv_prin = principal_value(rv)
        @test rotation_angle(rv_prin) <= pi
        @test rv_prin ≈ rv
    end
    rv = RotationVec(0.0, 0.0, 0.0)
    rv_prin = principal_value(rv)
    @test rotation_angle(rv_prin) <= pi
    @test rv_prin ≈ rv
end

@testset "Principal Value ($(rot_type))" for rot_type in [:RotX, :RotY, :RotZ] begin
        @eval begin
            for i = 1:1000
                r = $(rot_type)(100.0 * randn())
                r_prin = principal_value(r)
                @test -pi <= r_prin.theta <= pi
                @test r_prin ≈ r
            end
        end
    end
end

@testset "Principal Value ($(rot_type))" for rot_type in [:RotXY, :RotYX, :RotZX, :RotXZ, :RotYZ, :RotZY] begin
        @eval begin
            for i = 1:1000
                r = $rot_type(100.0 * randn(), 100.0 * randn())
                r_prin = principal_value(r)
                @test -pi <= r_prin.theta1 <= pi
                @test -pi <= r_prin.theta2 <= pi
                @test r_prin ≈ r
            end
        end
    end
end

@testset "Principal Value ($(rot_type))" for rot_type in [:RotXYX, :RotYXY, :RotZXZ, :RotXZX, :RotYZY, :RotZYZ, :RotXYZ, :RotYXZ, :RotZXY, :RotXZY, :RotYZX, :RotZYX] begin
        @eval begin
            for i = 1:1000
                r = $(rot_type)(100.0 * randn(), 100.0 * randn(), 100.0 * randn())
                r_prin = principal_value(r)
                @test -pi <= r_prin.theta1 <= pi
                @test -pi <= r_prin.theta2 <= pi
                @test -pi <= r_prin.theta3 <= pi
                @test r_prin ≈ r
            end
        end
    end
end


text:
@testset "arraymath.jl" begin
    @testset "inv(::Matrix{$T})" for T in (Float64, ComplexF64)
        B = generate_well_conditioned_matrix(T, 3)
        test_frule(inv, B)
        test_rrule(inv, B)
    end

    @testset "*: $T" for T in (Float64, ComplexF64)
        ⋆(a) = round.(5*randn(T, a))  # Helper to generate nice random values
        ⋆(a, b) = ⋆((a, b))  # matrix
        ⋆() = only(⋆(()))  # scalar

        @testset "Scalar-Array $dims" for dims in ((3,), (5,4), (2, 3, 4, 5))
            test_rrule(*, ⋆(), ⋆(dims))
            test_rrule(*, ⋆(dims), ⋆())
        end

        @testset "AbstractMatrix-AbstractVector n=$n, m=$m" for n in (2, 3), m in (4, 5)
            @testset "Array" begin
                test_rrule(*, n ⋆ m, ⋆(m))
            end
        end

        @testset "AbstractVector-AbstractMatrix n=$n, m=$m" for n in (2, 3), m in (4, 5)
            @testset "Array" begin
                test_rrule(*, ⋆(n), 1 ⋆ m)
            end
        end

        @testset "dense-matrix n=$n, m=$m, p=$p" for n in (2, 5), m in (2, 4), p in (2, 3)
            @testset "Array" begin
                test_rrule(*, (n⋆m), (m⋆p))
            end

            @testset "SubArray - $indexname" for (indexname, m_index) in (
                ("fast", :), ("slow", m:-1:1)
            )
                test_rrule(*, view(n⋆m, :, m_index), view(m⋆p, m_index, :))
                test_rrule(*, n⋆m, view(m⋆p, m_index, :))
                test_rrule(*, view(n⋆m, :, m_index), m⋆p)
            end

            @testset "Adjoints and Transposes" begin
                test_rrule(*, Transpose(m⋆n) ⊢ Transpose(m⋆n), Transpose(p⋆m) ⊢ Transpose(p⋆m))
                test_rrule(*, Adjoint(m⋆n) ⊢ Adjoint(m⋆n), Adjoint(p⋆m) ⊢ Adjoint(p⋆m))

                test_rrule(*, Transpose(m⋆n) ⊢ Transpose(m⋆n), (m⋆p))
                test_rrule(*, Adjoint(m⋆n) ⊢ Adjoint(m⋆n), (m⋆p))

                test_rrule(*, (n⋆m), Transpose(p⋆m) ⊢ Transpose(p⋆m))
                test_rrule(*, (n⋆m), Adjoint(p⋆m) ⊢ Adjoint(p⋆m))
            end
        end

        @testset "Diagonal" begin
            test_rrule(*, Diagonal([1.0, 2.0, 3.0]), Diagonal([4.0, 5.0, 6.0]))
            test_rrule(*, Diagonal([1.0, 2.0, 3.0]), rand(3))

            # Needs to not try and inplace, as `mul!` will do wrong.
            # see https://github.com/JuliaDiff/ChainRulesCore.jl/issues/411
            test_rrule(*, Diagonal([1.0, 2.0, 3.0]), rand(3,3))
        end

        @testset "Covector * Vector n=$n" for n in (3, 5)
            @testset "$f" for f in (adjoint, transpose)
                # This should be same as dot product and give a scalar
                test_rrule(*, f(⋆(n)) ⊢ f(⋆(n)), ⋆(n))
            end
        end
    end

    @testset "muladd: $T" for T in (Float64, ComplexF64)
        @testset "add $(typeof(z))" for z in [rand(T), rand(T, 3), rand(T, 3, 3), false]
            @testset "matrix * matrix" begin
                A = rand(T, 3, 3)
                B = rand(T, 3, 3)
                test_rrule(muladd, A, B, z)
                test_rrule(muladd, A', B, z)
                test_rrule(muladd, A , B', z)

                C = rand(T, 3, 5)
                D = rand(T, 5, 3)
                test_rrule(muladd, C, D, z)
            end
            if ndims(z) <= 1
                @testset "matrix * vector" begin
                    A, B = rand(T, 3, 3), rand(T, 3)
                    test_rrule(muladd, A, B, z)
                    test_rrule(muladd, A, B ⊢ rand(T, 3,1), z)
                end
                @testset "adjoint * matrix" begin
                    At, B = rand(T, 3)', rand(T, 3, 3)
                    test_rrule(muladd, At, B, z')
                    test_rrule(muladd, At ⊢ rand(T,1,3), B, z')
                end
            end
            if ndims(z) == 0
                @testset "adjoint * vector" begin # like dot
                    A, B = rand(T, 3)', rand(T, 3)
                    test_rrule(muladd, A, B, z)
                    test_rrule(muladd, A ⊢ rand(T,1,3), B, z')
                end
            end
            if ndims(z) == 2 # other dims lead to e.g. muladd(ones(4), ones(1,4), 1)
                @testset "vector * adjoint" begin # outer product
                    A, B = rand(T, 3), rand(T, 3)'
                    test_rrule(muladd, A, B, z)
                    test_rrule(muladd, A, B ⊢ rand(T,1,3), z)
                end
            end
        end
    end

    @testset "$f" for f in (/, \)
        @testset "Matrix" begin
            for n in 3:5, m in 3:5
                A = randn(m, n)
                B = randn(m, n)
                test_rrule(f, A, B; check_inferred=false) # ChainRulesCore #407
            end
        end
        @testset "Vector" begin
            x = randn(10)
            y = randn(10)
            test_rrule(f, x, y; check_inferred=false) # ChainRulesCore #407
        end
        if f == (\)
            @testset "Matrix $f Vector" begin
                X = randn(10, 4)
                y = randn(10)
                test_rrule(f, X, y)
            end
            @testset "Vector $f Matrix" begin
                x = randn(10)
                Y = randn(10, 4)
                test_rrule(f, x, Y; output_tangent=Transpose(rand(4)))
            end
        else
            A = rand(2, 4)
            B = rand(4, 4)
            test_rrule(f, A, B; check_inferred=false) # ChainRulesCore #407
        end
    end

    @testset "/ and \\ Scalar-AbstractArray" begin
        A = round.(10 .* randn(3, 4, 5), digits=1)
        test_rrule(/, A, 7.2)
        test_rrule(\, 7.2, A)

        C = round.(10 .* randn(6) .+ im .* 10 .* randn(6), digits=1)
        test_rrule(/, C, 7.2+8.3im)
        test_rrule(\, 7.2+8.3im, C)
    end

    @testset "negation" begin
        A = randn(4, 4)
        Ā = randn(4, 4)

        test_rrule(-, A)
        test_rrule(-, Diagonal(A); output_tangent=Diagonal(Ā))
    end

    @testset "addition" begin
        test_rrule(+, randn(4, 4), randn(4, 4), randn(4, 4))
        test_rrule(+, randn(3), randn(3,1), randn(3,1,1))
    end
end


text:
using Test, LinearAlgebra
using DiffEqSensitivity, StochasticDiffEq
using Random

@info "SDE Checkpointing"

seed = 100
Random.seed!(seed)

u₀ = [0.5]
tstart = 0.0
tend = 0.1
dt = 0.005
trange = (tstart, tend)
t = tstart:dt:tend
tarray = collect(t)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

p2 = [1.01,0.87]

f_oop_linear(u,p,t) = p[1]*u
σ_oop_linear(u,p,t) = p[2]*u

dt1 = tend/1e3

Random.seed!(seed)
prob_oop = SDEProblem(f_oop_linear,σ_oop_linear,u₀,trange,p2)
sol_oop = solve(prob_oop,EulerHeun(),dt=dt1,adaptive=false,save_noise=true)

@show length(sol_oop)

res_u0, res_p = adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

res_u0a, res_pa = adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()),
      checkpoints=sol_oop.t[1:2:end])

@test isapprox(res_u0, res_u0a, rtol = 1e-5)
@test isapprox(res_p, res_pa, rtol = 1e-2)

res_u0a, res_pa =  adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()),
      checkpoints=sol_oop.t[1:10:end])

@test isapprox(res_u0, res_u0a, rtol = 1e-5)
@test isapprox(res_p, res_pa, rtol = 1e-1)

dt1 = tend/1e4

Random.seed!(seed)
prob_oop = SDEProblem(f_oop_linear,σ_oop_linear,u₀,trange,p2)
sol_oop = solve(prob_oop,EulerHeun(),dt=dt1,adaptive=false,save_noise=true)

@show length(sol_oop)

res_u0, res_p = adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,
      sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

res_u0a, res_pa = adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()),
      checkpoints=sol_oop.t[1:2:end])

@test isapprox(res_u0, res_u0a, rtol = 1e-6)
@test isapprox(res_p, res_pa, rtol = 1e-3)

res_u0a, res_pa =  adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()),
      checkpoints=sol_oop.t[1:10:end])

@test isapprox(res_u0, res_u0a, rtol = 1e-6)
@test isapprox(res_p, res_pa, rtol = 1e-2)

res_u0a, res_pa =  adjoint_sensitivities(sol_oop,EulerHeun(),dg!,tarray
      ,dt=dt1,adaptive=false,sensealg=InterpolatingAdjoint(noise=DiffEqSensitivity.ZygoteNoise()),
      checkpoints=sol_oop.t[1:500:end])

@test isapprox(res_u0, res_u0a, rtol = 1e-3)
@test isapprox(res_p, res_pa, rtol = 1e-2)


text:
using Test, LinearFractionalTransformations

f = LFT(1,0,0,1)  # identity
@test f(1+im) == 1+im
@test f[1+im] == 1+im

f = LFT(2,0,0,1)  # f(x) = 2x
@test f(-im) == -2im

g = inv(f)
@test g(im) == 0.5im

g = f*f
@test g(2-im) == 4*(2-im)

f = LFT(1,2+im, 3,0, 1-im,5)
@test f(1)==2+im
@test f(3)==0
@test f(1-im)==5

g = inv(f)
@test g(0) == 3

@test LFT() == g*f

f = LFT(0,1,1,0)
@test f(Inf) == 0
@test isinf(f(0))

f = LFT(3-im, 2+im, 6)
@test f(3-im) == 0
@test f(2+im) == 1
@test isinf(f(6))


text:
@testset "sampler.jl" begin
    @testset "SampleFromPrior and SampleUniform" begin
        @model function gdemo(x, y)
            s ~ InverseGamma(2, 3)
            m ~ Normal(2.0, sqrt(s))
            x ~ Normal(m, sqrt(s))
            y ~ Normal(m, sqrt(s))
        end

        model = gdemo(1.0, 2.0)
        N = 1_000

        chains = sample(model, SampleFromPrior(), N; progress = false)
        @test chains isa Vector{<:VarInfo}
        @test length(chains) == N

        # Expected value of ``X`` where ``X ~ N(2, ...)`` is 2.
        @test mean(vi[@varname(m)] for vi in chains) ≈ 2 atol = 0.1

        # Expected value of ``X`` where ``X ~ IG(2, 3)`` is 3.
        @test mean(vi[@varname(s)] for vi in chains) ≈ 3 atol = 0.1

        chains = sample(model, SampleFromUniform(), N; progress = false)
        @test chains isa Vector{<:VarInfo}
        @test length(chains) == N

        # Expected value of ``X`` where ``X ~ U[-2, 2]`` is ≈ 0.
        @test mean(vi[@varname(m)] for vi in chains) ≈ 0 atol = 0.1

        # Expected value of ``exp(X)`` where ``X ~ U[-2, 2]`` is ≈ 1.8.
        @test mean(vi[@varname(s)] for vi in chains) ≈ 1.8 atol = 0.1
    end
    @testset "Initial parameters" begin
        # dummy algorithm that just returns initial value and does not perform any sampling
        struct OnlyInitAlg end
        function DynamicPPL.initialstep(
            rng::Random.AbstractRNG,
            model::Model,
            ::Sampler{OnlyInitAlg},
            vi::AbstractVarInfo;
            kwargs...,
            )
            return vi, nothing
        end
        DynamicPPL.getspace(::Sampler{OnlyInitAlg}) = ()

        # model with one variable: initialization p = 0.2
        @model function coinflip()
            p ~ Beta(1, 1)
            10 ~ Binomial(25, p)
        end
        model = coinflip()
        sampler = Sampler(OnlyInitAlg())
        lptrue = logpdf(Binomial(25, 0.2), 10)
        chain = sample(model, sampler, 1; init_params = 0.2, progress = false)
        @test chain[1].metadata.p.vals == [0.2]
        @test getlogp(chain[1]) == lptrue

        # parallel sampling
        chains = sample(
            model, sampler, MCMCThreads(), 1, 10;
            init_params = 0.2, progress = false,
        )
        for c in chains
            @test c[1].metadata.p.vals == [0.2]
            @test getlogp(c[1]) == lptrue
        end

        # model with two variables: initialization s = 4, m = -1
        @model function twovars()
            s ~ InverseGamma(2, 3)
            m ~ Normal(0, sqrt(s))
        end
        model = twovars()
        lptrue = logpdf(InverseGamma(2, 3), 4) + logpdf(Normal(0, 2), -1)
        chain = sample(model, sampler, 1; init_params = [4, -1], progress = false)
        @test chain[1].metadata.s.vals == [4]
        @test chain[1].metadata.m.vals == [-1]
        @test getlogp(chain[1]) == lptrue

        # parallel sampling
        chains = sample(
            model, sampler, MCMCThreads(), 1, 10;
            init_params = [4, -1], progress = false,
        )
        for c in chains
            @test c[1].metadata.s.vals == [4]
            @test c[1].metadata.m.vals == [-1]
            @test getlogp(c[1]) == lptrue
        end

        # set only m = -1
        chain = sample(model, sampler, 1; init_params = [missing, -1], progress = false)
        @test !ismissing(chain[1].metadata.s.vals[1])
        @test chain[1].metadata.m.vals == [-1]

        # parallel sampling
        chains = sample(
            model, sampler, MCMCThreads(), 1, 10;
            init_params = [missing, -1], progress = false,
        )
        for c in chains
            @test !ismissing(c[1].metadata.s.vals[1])
            @test c[1].metadata.m.vals == [-1]
        end
    end
end


text:
function dir_to_idx(dir::Symbol)
    if dir === :x
        1
    elseif dir === :y
        2
    elseif dir === :z
        3
    else
        0
    end
end

dir_to_idx(i::Int) = i

getdomain(g::AbstractGrid) = getfield(g, :grid)
getdomain(f::AbstractPICDataStructure) = getfield(f, :grid)
unwrapdata(f::AbstractPICDataStructure) = getfield(f, :data)

function newstruct(f::AbstractPICDataStructure, data, grid; name=:same)
    name = name == :same ? nameof(f) : name
    parameterless_type(f)(data, grid, name)
end

Base.nameof(f::AbstractPICDataStructure) = getfield(f, :name)

"""
    axisnames(f::AbstractPICDataStructure; include_units=true)

Get the names of the axis of the grid corresponding to the input data structure `f`.
In the case of Unitful quantities, they can be excluded by setting `include_units` to `false`.
"""
axisnames(f::AbstractPICDataStructure; include_units=true) = axisnames(getdomain(f); include_units)

function unwrap(f)
    _f = hasunit(f) ? ustrip(f) : f

    grid = getdomain(_f)
    data = unwrapdata(_f)

    return grid, data
end

expandgrid(f::Observable{T}) where T = expandgrid(domain_discretization(T), f)

function expandgrid(::LatticeGrid, f)
    dirs = propertynames(getdomain(f[]))
    map(dirs) do dir
        lift(f) do val_f
            g = getdomain(val_f)
            getproperty(g, dir)
        end
    end
end

function expandgrid(::ParticleGrid, f)
    @lift getdomain($f)
end

@inline function apply_reinterpret(f, data)
    T = recursive_unitless_bottom_eltype(data)
    data_units = recursive_bottom_unit(data)
    raw_data = reinterpret(T, data)
    resized_data = f(raw_data)
    U = typeof(one(T) * data_units)
    reinterpret(U, resized_data)
end

function unwrap(f::Observable)
    _f = @lift hasunit($f) ? ustrip($f) : $f

    grid = expandgrid(_f)
    data = @lift unwrapdata($_f)

    return grid, data
end

dimensionality(::AbstractGrid{N}) where N = N
dimensionality(::Type{<:AbstractGrid{N}}) where N = N
dimensionality(::AbstractPICDataStructure{T,N,G}) where {T,N,G} = dimensionality(G)
dimensionality(::Type{<:AbstractPICDataStructure{T,N,G}}) where {T,N,G} = dimensionality(G)
# use Tuple{T, Vararg{T,N}} since NTuple{N,T} can have N==0
dimensionality(::Type{<:NamedTuple{Names, Tuple{T, Vararg{T,N}}}}) where {Names, N, T} = N + 1

function mapgrid(f, grid::AbstractAxisGrid)
    map(f, Iterators.product(grid...))
end

function mapgrid(f, grid::ParticlePositions)
    map(f, zip(grid...))
end

mapgrid(f, field::AbstractPICDataStructure) = mapgrid(f, getdomain(field))

function scalarfield(f, grid; name="")
    data = mapgrid(f, grid)
    ScalarField(data, grid, name)
end

function scalarvariable(f, grid; name="")
    data = mapgrid(f, grid)
    ScalarVariable(data, grid, name)
end

function vectorfield(f, grid; name="")
    data = mapgrid(f, grid)
    VectorField(data, grid; name)
end


text:
# Mountain Car problem for continuous reinforcement learning
# As described in XXX

@with_kw struct MountainCar <: MDP{Tuple{Float64,Float64},Float64}
  discount::Float64 = 0.99
  cost::Float64 = -1.# reward at each state not at the goal (should be a negative number)
  jackpot::Float64 = 0.0 # reward at the top
end

actions(::MountainCar) = [-1., 0., 1.]

reward(mc::MountainCar,
       s::Tuple{Float64,Float64},
       a::Float64,
       sp::Tuple{Float64,Float64}) = isterminal(mc,sp) ? mc.jackpot : mc.cost

initialstate(mc::MountainCar, ::AbstractRNG) = (-0.5,0.,)

isterminal(::MountainCar,s::Tuple{Float64,Float64}) = s[1] >= 0.5
discount(mc::MountainCar) = mc.discount

function gen(::DDNNode{:sp},
             mc::MountainCar,
             s::Tuple{Float64,Float64},
             a::Float64,
             ::AbstractRNG)
  x,v = s
  v_ = v + a*0.001+cos(3*x)*-0.0025
  v_ = max(min(0.07,v_),-0.07)
  x_ = x+v_
  #inelastic boundary
  if x_ < -1.2
      x_ = -1.2
      v_ = 0.
  end
  sp = (x_,v_,)
  return sp
end

function convert_s(::Type{A}, s::Tuple{Float64,Float64}, mc::MountainCar) where A<:AbstractArray
    v = copyto!(A(undef, 2), s)
    return v
end
convert_s(::Type{Tuple{Float64,Float64}}, s::A, mc::MountainCar) where A<:AbstractArray = (s[1], s[2])

# Example policy -- works pretty well
mutable struct Energize <: Policy end
action(::Energize,s::Tuple{Float64,Float64}) = sign(s[2])


text:
@symbol_func function a_P(cur_reactor::AbstractReactor)
  cur_a = 1.04

  cur_a *= a(cur_reactor)

  cur_a
end


text:
using EchelleCCFs
using RvSpectMLBase
using Test

@testset "Utilities" begin
    @testset "Physics" begin
        @test EchelleCCFs.λ_vac_to_air(EchelleCCFs.λ_air_to_vac(5000.0)) ≈ 5000.0
    end

    @testset "Make mask utils" begin
        Δv = 7e3
        Δz = Δv/EchelleCCFs.speed_of_light_mps
        @test_nowarn EchelleCCFs.ChunkWidthFixedΔlnλ(Δz)
    end

end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# project.jl

using SPECTrecon: SPECTplan
using SPECTrecon: project!
using MAT
using LinearAlgebra: norm

function cal_nrmse(x, xtrue)
    len = size(xtrue, 3)
    err = zeros(eltype(xtrue), len)
    for l = 1:len
        err[l] = norm(vec(x[:,:,l]) - vec(xtrue[:,:,l])) / norm(vec(xtrue[:,:,l]))
    end
    nrmse = sum(err) / len
    return nrmse
end

function project_error()
    T = Float32
    path = "../data/"
    file = matopen(path*"mumap208.mat")
    mumap = read(file, "mumap208")
    close(file)

    file = matopen(path*"psf_208.mat")
    psfs = read(file, "psf_208")
    close(file)

    file = matopen(path*"xtrue.mat")
    xtrue = convert(Array{Float32, 3}, read(file, "xtrue"))
    close(file)

    file = matopen(path*"proj_jeff_newmumap.mat")
    proj_jeff = read(file, "proj_jeff")
    close(file)

    dy = T(4.7952)
    nview = size(psfs, 4)

    (nx, ny, nz) = size(xtrue)

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            println(string(interpmeth)*", "*string(mode))
            plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
            views = zeros(T, nx, nz, nview)
            project!(views, xtrue, plan)
            nrmse = cal_nrmse(views, proj_jeff)
            @show nrmse
        end
    end

    # 1d interp: 2e-7 nrmse
    # 2d interp: 4e-3 nrmse
end

# shows accuracy
project_error()
nothing
# plot(nrmse1d)


text:
function iliev03()
#
#  generalization of Iliev example:
#
#    (x-1) (x-2)^2 (x-3)^3
#
    iliev(8)
end


text:
include("problems.jl")

function Basic_GRAPE(prob)
        ctrl_amps = prob.init_guess
        init_guess_amps = copy(ctrl_amps)
        n_ctrl_Hams = length(prob.ctrl_Hams)
        Δt = prob.transfer_time/prob.n_steps

        for ctrl_Ham in 1:n_ctrl_Hams
                for timeslice in 1:prob.n_steps
                        λ = back_target_op(prob, ctrl_amps, timeslice)
                        ρ = density_op(prob, ctrl_amps, timeslice)
                        thre = 1 - performance_func(prob, ctrl_amps, timeslice)
                        if thre < 1 - prob.min_fid
                                println("solution converged at step = ", ctrl_Ham*timeslice)
                                @goto esc
                        end
                        A = λ
                        B = -1im * Δt * commutator(prob.ctrl_Hams[ctrl_Ham], ρ)
                        ctrl_amps[timeslice, ctrl_Ham] = ctrl_amps[timeslice, ctrl_Ham] - expval(A,B)
                end
        end
        println("solution couldn't be found try another parameters")
        @label esc
        return init_guess_amps, ctrl_amps
end

function GRAPE(prob; noise=false)
        if noise
                # to be added ....
        else
                init_guess_amps, ctrl_amps = Basic_GRAPE(prob)
                return init_guess_amps, ctrl_amps
        end
end

function C1_grad(prob, ctrl_amps)
        grad = Array{Float64}(undef, prob.n_steps, length(prob.ctrl_Hams))
        P = tr(prob.target_density_op'*prob.init_density_op)*prob.target_density_op'
        X = prob.init_density_op
        Δt = prob.transfer_time/prob.n_steps
        for j in reverse(1:prob.n_steps)
                for k in 1:length(prob.ctrl_Hams)
                        grad[j,k] = 2*Δt*tr(P*prob.ctrl_Hams[k]*X)
                end
                X = evolve(prob, ctrl_amps, j)'*X
                P = P*evolve(prob, ctrl_amps, j)
        end
        return grad
end

function GRAPE_auto_diff(prob)
        grad = Array{Float64}(undef, prob.n_steps, length(prob.ctrl_Hams))
        P = tr(prob.target_density_op'*prob.init_density_op)*prob.target_density_op'
        X = prob.init_density_op
        Δt = prob.transfer_time/prob.n_steps

        ctrl_amps = prob.init_guess
        init_guess_amps = copy(ctrl_amps)

        infid = C1_cost(X, prob.target_density_op, 2)

        while infid > 1
                println(infid)
                for j in reverse(1:prob.n_steps)
                        for k in 1:length(prob.ctrl_Hams)
                                grad[j,k] = 2*Δt*tr(P*prob.ctrl_Hams[k]*X)
                                ctrl_amps[j,k] = ctrl_amps[j,k] + 3*grad[j,k]
                        end
                        X = evolve(prob, ctrl_amps, j)'*X
                        P = P*evolve(prob, ctrl_amps, j)
                end
                infid = C1_cost(X, prob.target_density_op, 2)
        end
        #return ctrl_amps
end

#=
function AUTO_DIFF(prob)
        cost(u) = 1 - density_op_fid(prob, u, prob.init_density_op, prob.target_density_op)
        init_ctrl_arr = prob.init_guess

        res = optimize(cost, init_ctrl_arr)
        final_ctrl_arr = Optim.minimizer(res)

        return init_ctrl_arr, final_ctrl_arr
end
=#

function OPTIMIZE(prob::GateTransfer)
    cost(u) = gate_infid(prob.drift_Ham, prob.ctrl_Hams, prob.n_steps, prob.transfer_time, u, prob.tar_gate)
    init_ctrl_arr = prob.init_guess

    res = optimize(cost, init_ctrl_arr)
    final_ctrl_arr = Optim.minimizer(res)
    fid = 1 - Optim.minimum(res)
    println("Fidelity = ", fid)
    return init_ctrl_arr, final_ctrl_arr
end

function OPTIMIZE(prob::StateTransfer)
    cost(u) = state_infid(prob.drift_Ham, prob.ctrl_Hams, prob.n_steps, prob.transfer_time, u, prob.init_state, prob.tar_state)
    init_ctrl_arr = prob.init_guess

    res = optimize(cost, init_ctrl_arr)
    final_ctrl_arr = Optim.minimizer(res)
    fid = 1 - Optim.minimum(res)
    println("Fidelity = ", fid)
    return init_ctrl_arr, final_ctrl_arr
end