text:
using Base.Threads
using Distributed

include("util.jl")
include("./dataset.jl")
include("./iterator.jl")

function get_dataset(datasetdir)
    train_paris, val_pairs = make_trainval_pairs(datasetdir)
    train_dataset = Dataset(train_paris)
    val_dataset = Dataset(val_pairs)
    return train_dataset, val_dataset
end

function cachecheck()
    datasetdir = expanduser("~/dataSSD120GB/food-101")
    train_dataset, val_dataset = get_dataset(datasetdir)
    println("Hello")
    batchsize=32
    mag = 10
     val_iter = SerialIterator(val_dataset,batchsize*mag,shuffle=false)
    @time for batch in val_iter
        break
    end

     val_iter = SerialIterator(val_dataset,batchsize*mag,shuffle=false)
    @time for batch in val_iter
        break
    end

     val_iter = SerialIterator(val_dataset,batchsize*mag,shuffle=false)
    @time for batch in val_iter
        break
    end
    println(typeof(val_dataset.data[1]))
    println(typeof(val_iter.get_example))
    #println(sort([k for k in keys(val_dataset.image_cache)]))
end

#cachecheck()
dic=Dict{Int,Any}()
datasetdir = expanduser("~/dataSSD120GB/food-101")
train_dataset, val_dataset = get_dataset(datasetdir)
for (i, (f,l)) in enumerate(val_dataset.data[1:100])
    dic[i]=load(f)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
	fourier_matrix(n::Int; normalized=true)

Returns a ``n``-by-``n`` Fourier matrix with optional normalization (`true` by default).
"""
function fourier_matrix(n::Int; normalized = true)

    U = Array{ComplexF64}(undef, n, n)

    for i in 1:n
        for j in 1:n
            U[i,j] = exp((2im * pi / n)* (i-1) * (j-1))
        end
    end

    if normalized
        return 1/sqrt(n) * U
    else
        return U
    end
end

"""
	hadamard_matrix(n::Int; normalized=true)

Returns a ``n``-by-``n`` Hadamard matrix with optional normalization (`true` by default).
"""
function hadamard_matrix(n::Int, normalized = true )

	H = Array{ComplexF64}(undef, n, n)

	for i in 1:n
		for j in 1:n
			bi = map(x->parse(Int,x), split(bitstring(Int8(i-1)),""))
			bj = map(x->parse(Int,x), split(bitstring(UInt8(j-1)),""))
			H[i,j] = (-1)^(sum(bi.&bj))
		end
	end

	if normalized
		return 1/sqrt(n) * H
	else
		return H
	end

end

"""
	sylvester_matrix(p::Int; normalized=true)

Returns a ``2^p``-by-``2^p`` Sylvester matrix with optional normalization
(`true` by default) following [https://arxiv.org/abs/1502.06372](https://arxiv.org/abs/1502.06372).
"""
function sylvester_matrix(p; normalized = true)

    """returns the sylvester matrix of size 2^p"""

    function sylvester_element(i,j,p)

        """returns the element A(i,j) for the sylvester matrix of size 2^p"""
        # following https://arxiv.org/abs/1502.06372
        if p == 0
            return 1
        else
            i_b = digits(i-1, base = 2, pad = p+1)
            j_b = digits(j-1, base = 2, pad = p+1)
            return (-1)^(sum(i_b .* j_b))
        end
    end

    A = Array{ComplexF64}(undef, 2^p, 2^p)
    for i in 1:2^p
        for j in 1:2^p
            A[i,j] = sylvester_element(i,j,p)
        end
    end

    normalized == true ? 1/sqrt(2^p) * A : A

end

"""
	rand_haar(n::Int)

Returns a ``n``-by-``n`` Haar distributed unitary matrix following [https://case.edu/artsci/math/esmeckes/Meckes_SAMSI_Lecture2.pdf](https://case.edu/artsci/math/esmeckes/Meckes_SAMSI_Lecture2.pdf).
"""
function rand_haar(n::Int)

    """generates a Haar distributed unitary matrix of size n*n"""
    # follows https://case.edu/artsci/math/esmeckes/Meckes_SAMSI_Lecture2.pdf

    qr(randn(ComplexF64, n,n)).Q
end

function matrix_test(n::Int)

    """matrix of 1,2,3... for testing your code"""

    U = Array{Int64}(undef, n, n)

    for i in 1:n
        for j in 1:n
            U[i,j] = (i-1) * n + j
        end
    end

    U

end

function antihermitian_test_matrix(n::Int)

    h = randn(ComplexF64, n, n)

    for i = 1:n
        h[i,i] = 1im*imag(h[i,i])
    end

    for i = 1:n
        for j = i+1:n
            h[j,i] = -conj(h[i,j])
        end
    end

    h

end

function hermitian_test_matrix(n::Int)

    h = randn(ComplexF64, n,n)

    for i = 1:n
        h[i,i] = real(h[i,i])
    end

    for i = 1:n
        for j = i+1:n
            h[j,i] = conj(h[i,j])
        end
    end

    h

end

function permutation_matrix_special_partition_fourier(n::Int)

    """P * [1,1,1,1,1,0,0,0,0,0] = [1, 0, 1, 0, 1, 0, 1, 0, 1, 0]"""

    P = zeros(Int, n,n)

    for i in 1:2:n
        P[i, Int((i+1)/2)] = 1
    end

    for i in 2:2:n
        P[i, Int(n/2 + i/2)] = 1
    end

    P

end

"""
	rand_gram_matrix_from_orthonormal_basis(n::Int, r::Int)

Returns a ``n``-by-``n`` random Gram matrix that generates orthonormal basis of `r` vectors.
"""
function rand_gram_matrix_from_orthonormal_basis(n,r)

	"""generates a random gram matrix with the property that the
	generating vector basis of r vectors is orthonormal (as is
	strangely the case in Drury's work)"""

	function normalized_random_vector(n)
		v = rand(ComplexF64, n)
		1/norm(v) .* v
	end

	if r >= n
		throw(ArgumentError("need rank < dim to have a non trivial result"))
	end

	generating_vectors = 0.
	#while det(generating_vectors) ≈ 0. #check it's a basis
	    generating_vectors = hcat([normalized_random_vector(n) for i = 1:r]...)
	#### TODO caveat : we don't check for linear independance so there is a chance it doesn't work though highly unlikely
	#end

	generating_vectors = modified_gram_schmidt(generating_vectors)

	is_orthonormal(generating_vectors)

	generating_vectors * generating_vectors'

end

function modified_gram_schmidt(input_vectors)

	"""does the modified gram schmidt (numerically stable) on the columns of the matrix input_vectors"""

	function projector(a, u)
		"""projector of a on u"""
		dot(u,a)/dot(u,u) .* u
	end

	final_vectors = copy(input_vectors)

	for i in 1:size(input_vectors)[2]

		if i == 1
			#normalize first vector
			final_vectors[:,1] /= norm(final_vectors[:, 1])
		else
			for j in 1:i-1
				final_vectors[:, i] -= projector(final_vectors[:, i], final_vectors[:, j])
			end
			final_vectors[:, i] /= norm(final_vectors[:, i])

		end

	end

	is_orthonormal(final_vectors)

	final_vectors

end

"""
	rand_gram_matrix(n::Int)

Returns a full rank ``n``-by-``n`` random Gram matrix.
"""
function rand_gram_matrix(n::Int)

	"""random gram matrix with full rank """

	function normalized_random_vector(n)
		v = rand(ComplexF64, n)
		1/norm(v) .* v
	end

	generating_vectors = 0.
	while det(generating_vectors) ≈ 0. #check it's a basis
	    generating_vectors = hcat([normalized_random_vector(n) for i = 1:n]...)

	end

	generating_vectors' * generating_vectors

end

"""
	rand_gram_matrix_rank(n::Int, r::Int)

Returns a ``n``-by-``n`` random Gram matrix of maximum rank and great likelihood `r`.
"""
function rand_gram_matrix_rank(n, r)

	"""random gram matrix with rank at most r, and with great likelihood r
	"""

	function normalized_random_vector(r)
		v = rand(ComplexF64, r)
		1/norm(v) .* v
	end

	generating_vectors = 0.
	generating_vectors = hcat([normalized_random_vector(r) for i = 1:n]...)

	generating_vectors' * generating_vectors

end

"""
	rand_gram_matrix_real(n::Int)

Returns a real ``n``-by-``n`` random Gram matrix of full rank.
"""
function rand_gram_matrix_real(n)

	"""random gram matrix with full rank and real el"""

	function normalized_random_vector(n)
		v = rand(Float64, n)
		1/norm(v) .* v
	end

	generating_vectors = 0.
	while det(generating_vectors) ≈ 0. #check it's a basis
	    generating_vectors = hcat([normalized_random_vector(n) for i = 1:n]...)

	end

	generating_vectors' * generating_vectors

end

"""
	rand_gram_matrix_positive(n::Int)

Returns a positive elements ``n``-by-``n`` random Gram matrix of full rank.
"""
function rand_gram_matrix_positive(n::Int)

	"""random gram matrix with full rank and positive el"""

	function normalized_random_vector(n)
		v = abs.(rand(Float64, n))
		1/norm(v) .* v
	end

	generating_vectors = 0.
	while det(generating_vectors) ≈ 0. #check it's a basis
	    generating_vectors = hcat([normalized_random_vector(n) for i = 1:n]...)

	end

	generating_vectors' * generating_vectors

end

# functions to clear matrices for better readability of a matrix by a human

function set_small_values_to_zero(x, eps)
	abs(x) < eps ? 0.0 : x
end

function set_small_imag_real_to_zero(x, eps = 1e-7)
	if abs(real(x)) < eps
		x = 0.0 + 1im * imag(x)
	end
	if abs(imag(x)) < eps
		x = real(x) + 0.0im
	end
	x
end

function make_matrix_more_readable(M)
	set_small_imag_real_to_zero.(M)
end

function display_matrix_more_readable(M)
	pretty_table(make_matrix_more_readable(M))
end

####### sub matrices ########

function remove_row_col(A, rows_to_remove, cols_to_remove)

	"""hands back the matrix A without the rows labelled in rows_to_remove, columns in cols_to_remove

	ex :
	rows_to_remove = [1]
	cols_to_remove = [2]

	to have A(1,2) in the notations of Minc"""

	n = size(A)[1]
	m = size(A)[2]

	range_rows = collect(1:n)
	range_col = collect(1:m)

	setdiff!(range_rows, rows_to_remove)
	setdiff!(range_col, cols_to_remove)

	A[range_rows, range_col]

end

### permanent related quantities ###

function sub_permanent_matrix(A)

	n = size(A)[1]

	sub_permanents = similar(A)

	for i in 1 : n
		for j in 1 : n

			sub_permanents[i,j] = ryser(remove_row_col(A, [i], [j]))

		end
	end

	sub_permanents

end

function gram_from_n_r_vectors(M)

	# generate a random matrix of n r-vectors M
	# normalize it to have a gram matrix

	n = size(M)[1]

	for i in  1 : n
		M[:,i] ./= norm(M[:,i])
	end

	M' * M # our rank r n*n gram matrix

end

"""
	gram_matrix_one_param(n::Int, x::Real)

Returns a ``n``-by-``n`` Gram matrix parametrized by the real ``0 ≤ x ≦ 1``.
"""
function gram_matrix_one_param(n::Int, x::Real)

	@argcheck x>=0 && x<= 1

	S = 1.0 * Matrix(I, n, n)
	for i in 1:n
		for j in 1:n
			i!=j ? S[i,j] = x : continue
		end
	end

	S

end

function column_normalize(M)
   for col in 1:size(M,2) #column normalize
	   M[:, col] = M[:, col]./norm(M[:, col])
   end
   M
end

"""
	perturbed_gram_matrix(M, epsilon)

Returns a Gram matrix generated by the columns of `M` which are perturbed by
a Gaussian quantity of variance epsilon once normalized.
"""
function perturbed_gram_matrix(M, epsilon)

    """M defines the set of vectors generating the gram matrix
    each column is a generating vector for the gram matrix
    we perturb them by some random gaussian amount with set variance epsilon once normalized """

    M = column_normalize(M)

    d = Normal(0.0, epsilon)

    perturbation_vector = rand(d,size(M))

    M += perturbation_vector

    M = column_normalize(M)

    M' * M

end

"""
	perturbed_unitary(U, epsilon)

Returns a unitary matrix whose columns are generating vector perturbed by a
random Gaussian quantity with variance epsilon once normalized.
"""
function perturbed_unitary(U, epsilon)

    """U a unitary matrix each column is a generating vector
    we perturb them by some random gaussian amount with set variance epsilon once normalized """

    d = Normal(0.0, epsilon)

    perturbation_vector = rand(d,size(U))

    U += perturbation_vector

    U = modified_gram_schmidt(U)

    U

end

ext: jl


text:
struct WasmImport
    wasm_importtype_ptr::Ptr{wasm_importtype_t}
    extern_kind::wasm_externkind_enum
    import_module::String
    name::String

    function WasmImport(wasm_importtype_ptr::Ptr{wasm_importtype_t})
        name_vec_ptr = wasm_importtype_name(wasm_importtype_ptr)
        name = name_vec_ptr_to_str(name_vec_ptr)

        import_module_ptr = wasm_importtype_module(wasm_importtype_ptr)
        import_module = name_vec_ptr_to_str(import_module_ptr)

        externtype_ptr = wasm_importtype_type(wasm_importtype_ptr)
        extern_kind = wasm_externkind_enum(wasm_externtype_kind(externtype_ptr))

        new(wasm_importtype_ptr, extern_kind, import_module, name)
    end
end

Base.unsafe_convert(::Type{Ptr{wasm_importtype_t}}, wasm_import::WasmImport) = wasm_import.wasm_importtype_ptr
Base.show(io::IO, wasm_import::WasmImport) = print(
    io,
    "WasmImport($(wasm_import.extern_kind), \"$(wasm_import.import_module)\", \"$(wasm_import.name)\")",
)

struct WasmImports{M<:AbstractWasmModule}
    wasm_module::M
    wasm_imports::Vector{WasmImport}

    function WasmImports(wasm_module::M, wasm_imports_vec) where {M<:AbstractWasmModule}
        wasm_imports = map(imp -> wasm_importtype_copy(imp) |> WasmImport, wasm_imports_vec)
        new{M}(wasm_module, wasm_imports)
    end
end

function Base.show(io::IO, wasm_imports::WasmImports)
    print(io::IO, "WasmImports(")
    show(io, wasm_imports.wasm_imports)
    print(")")
end


text:
# SPDX-License-Identifier: MIT
# Copyright (c) 2020: Pablo Zubieta

module ReactionCoordinates

using LinearAlgebra
using StaticArrays

export Acylindricity, Angle, Asphericity, Barycenter, DihedralAngle, DistanceFrom,
    GyrationTensor, PairwiseKernel, PrincipalMoments, RadiusOfGyration,
    RouseMode, Separation, ShapeAnysotropy, TorsionAngle, WeightedBarycenter,
    WeightedGyrationTensor, WeightedRadiusOfGyration

### Type definitions

abstract type ReactionCoordinate <: Function end
abstract type AbstractGyrationTensor <: ReactionCoordinate end
abstract type AbstractBarycenter <: ReactionCoordinate end

struct Angle <: ReactionCoordinate end
struct DihedralAngle <: ReactionCoordinate end
struct GyrationTensor <: AbstractGyrationTensor end
struct RadiusOfGyration{T} <: ReactionCoordinate end
struct Separation{T} <: ReactionCoordinate end

const TorsionAngle = DihedralAngle

struct WeightedGyrationTensor{T} <: AbstractGyrationTensor
    ws::Vector{T}
end

struct WeightedRadiusOfGyration{T} <: ReactionCoordinate
    ws::Vector{T}
end

struct PrincipalMoments{T <: AbstractGyrationTensor} <: ReactionCoordinate
    S::T
end

struct Asphericity{T <: AbstractGyrationTensor} <: ReactionCoordinate
    S::T
end

struct Acylindricity{T <: AbstractGyrationTensor} <: ReactionCoordinate
    S::T
end

struct ShapeAnysotropy{T <: AbstractGyrationTensor} <: ReactionCoordinate
    S::T
end

struct Barycenter{F} <: AbstractBarycenter
    f::F
end

struct WeightedBarycenter{F, T} <: AbstractBarycenter
    f::F
    ws::Vector{T}
end

struct PairwiseKernel{F} <: ReactionCoordinate
    f::F
end

struct DistanceFrom{T} <: ReactionCoordinate
    r::SVector{T}
end

struct RouseMode{N, P, T} <: ReactionCoordinate
    c::T
    ws::Vector{T}

    function RouseMode{N, 0}(::Type{T} = Float64) where {N, T}
        @assert (N > 0) "Parameters must satisfy N > P ≥ 0. Got N = $N and P = 0"
        c = sqrt(one(T) / N)
        return new{N, 0, T}(c, T[])
    end

    function RouseMode{N, P}(::Type{T} = Float64) where {N, P, T}
        @assert (N > P ≥ 0) "Parameters must satisfy N > P ≥ 0. Got N = $N and P = $P"
        c = sqrt(2 * one(T) / N)
        ws = T[cospi(P * (2i - 1) / 2N) for i = 1:N]
        return new{N, P, T}(c, ws)
    end
end

struct TransformMatrix{T} <: FieldMatrix{3, 3, T}
    xx::T
    yx::T
    zx::T
    xy::T
    yy::T
    zy::T
    xz::T
    yz::T
    zz::T
end

#=

Explore using FunctionalWrappers.jl for this

struct MultiReactionCoordinate{T <: Tuple} <: ReactionCoordinate
    ξs::T
end

=#

### Outer Constructors

function WeightedGyrationTensor(ws::Vector)
    w̃s = ws / sum(ws)  # weights normalized at construction
    T = eltype(w̃s)
    return WeightedGyrationTensor{T}(w̃s)
end

function WeightedRadiusOfGyration(ws::Vector)
    w̃s = ws / sum(ws)
    T = eltype(w̃s)
    return WeightedRadiusOfGyration{T}(w̃s)
end

function WeightedBarycenter(ws::Vector; f::F = identity) where {F}
    w̃s = ws / sum(ws)
    T = eltype(w̃s)
    return WeightedBarycenter{F, T}(f, w̃s)
end

PrincipalMoments() = PrincipalMoments(GyrationTensor())
PrincipalMoments(ws::Vector) = PrincipalMoments(WeightedGyrationTensor(ws))

Asphericity() = Asphericity(GyrationTensor())
Asphericity(ws::Vector) = Asphericity(WeightedGyrationTensor(ws))

Acylindricity() = Acylindricity(GyrationTensor())
Acylindricity(ws::Vector) = Acylindricity(WeightedGyrationTensor(ws))

ShapeAnysotropy() = ShapeAnysotropy(GyrationTensor())
ShapeAnysotropy(ws::Vector) = ShapeAnysotropy(WeightedGyrationTensor(ws))

Barycenter() = Barycenter(identity)

### Getters

gyration_tensor(ξ::PrincipalMoments) = ξ.S
gyration_tensor(ξ::Asphericity) = ξ.S
gyration_tensor(ξ::Acylindricity) = ξ.S
gyration_tensor(ξ::ShapeAnysotropy) = ξ.S

weights(ξ::WeightedGyrationTensor) = ξ.ws
weights(ξ::WeightedRadiusOfGyration) = ξ.ws
weights(ξ::WeightedBarycenter) = ξ.ws
weights(ξ::RouseMode) = ξ.ws

op(ξ::AbstractBarycenter) = ξ.f
op(ξ::PairwiseKernel) = ξ.f

reference(ξ::DistanceFrom) = ξ.r

coeff(ξ::RouseMode) = ξ.c

### Functors methods

(ξ::Angle)(rs) = @inbounds angle(rs[1], rs[2], rs[3])
(ξ::DihedralAngle)(rs) = @inbounds dihedral_angle(rs[1], rs[2], rs[3], rs[4])
(ξ::GyrationTensor)(rs) = gyration_tensor(rs)
(ξ::WeightedGyrationTensor)(rs) = gyration_tensor(rs, weights(ξ))
(ξ::RadiusOfGyration)(rs) = radius_of_gyration(rs)
(ξ::WeightedRadiusOfGyration)(rs) = radius_of_gyration(rs, weights(ξ))
(ξ::PrincipalMoments)(rs) = principal_moments(gyration_tensor(ξ)(rs))
(ξ::Asphericity)(rs) = asphericity(gyration_tensor(ξ)(rs))
(ξ::Acylindricity)(rs) = acylindricity(gyration_tensor(ξ)(rs))
(ξ::ShapeAnysotropy)(rs) = shape_anysotropy(gyration_tensor(ξ)(rs))
(ξ::Barycenter)(rs) = barycenter(op(ξ), rs)
(ξ::WeightedBarycenter)(rs) = barycenter(op(ξ), rs, weights(ξ))
(ξ::PairwiseKernel)(rs₁, rs₂) = pairwise(op(ξ), rs₁, rs₂)
(ξ::DistanceFrom)(rs) = distance(rs, reference(ξ))
(ξ::Separation)(rs) = @inbounds distance(rs[1], rs[2])
(ξ::RouseMode)(rs) = rouse_mode(rs, coeff(ξ), weights(ξ))
(ξ::RouseMode{N, 0})(rs) where {N} = rouse_mode₀(rs, coeff(ξ))

#==========#
#  Angles  #
#==========#

"""    angle(r₁, r₂, r₃)

Computes the angle between the two vectors defined by three points in space (around the
point in the middle).
"""
angle(r₁, r₂, r₃) = angle(r₁ - r₂, r₃ - r₂)
#
@inline angle(a, b) = atan(norm(a × b), a ⋅ b)

"""    dihedral_angle(r₁, r₂, r₃, r₄)

Computes the dihedral (or torsion) angle defined by four points in space (around the line
defined by the two central points).
"""
dihedral_angle(r₁, r₂, r₃, r₄) = dihedral_angle(r₂ - r₁, r₃ - r₂, r₄ - r₃)
#
@inline function dihedral_angle(a, b, c)
    p = a × b
    q = b × c
    return atan((p × q) ⋅ b, (p ⋅ q) * norm(b))
end

#==============#
#  Box Volume  #
#==============#

volume(H::TransformMatrix) = det(H)

#=====================#
#  Shape Descriptors  #
#=====================#

function gyration_tensor(rs)
    # Alternative implementation:
    # f = r -> r .* r'
    # S = sum(f, rs)
    S = accumulator(GyrationTensor, rs)
    @simd for r in rs
        S .= muladd.(r, r', S)
    end
    return Symmetric(SMatrix(S ./= length(rs)))
end

function gyration_tensor(rs, ws)
    # Alternative implementation:
    # f = ((w, r),) -> w .* r .* r'
    # S = sum(f, rs)
    S = accumulator(WeightedGyrationTensor, rs, ws)
    @inbounds @simd for i in eachindex(ws, rs)
        w, r = ws[i], rs[i]
        S .= muladd.(w, r .* r', S)
    end
    return Symmetric(SMatrix(S))
end

radius_of_gyration(rs) = sum(r -> r ⋅ r, rs) / length(rs)

function radius_of_gyration(rs, ws)
    # Alternative implementation:
    # f = ((w, r),) -> w * (r ⋅ r)
    # R² = sum(f, rs)
    R² = accumulator(WeightedRadiusOfGyration, rs, ws)
    @inbounds @simd for i in eachindex(ws, rs)
        w, r = ws[i], rs[i]
        R² = muladd(w, r ⋅ r, R²)
    end
    return R²
end

const principal_moments = LinearAlgebra.eigvals

function asphericity(S)
    λ₁², λ₂², λ₃² = principal_moments(S)
    return λ₃² - (λ₁² + λ₂²) / 2
end

function acylindricity(S)
    λ₁², λ₂², λ₃² = principal_moments(S)
    return (λ₂² - λ₁²)
end

function shape_anysotropy(S)
    λ₁², λ₂², λ₃² = principal_moments(S)
    λ₁⁴ = λ₁²^2
    λ₂⁴ = λ₂²^2
    λ₃⁴ = λ₃²^2
    return (3 * (λ₁⁴ + λ₂⁴ + λ₃⁴) / (λ₁² + λ₂² + λ₃²)^2 - 1) / 2
end

### Accumulators

@inline function accumulator(::Type{GyrationTensor}, rs)
    R = eltype(eltype(rs))
    T = typeof(zero(R) / 1)
    return zeros(MMatrix{3, 3, T})
end

@inline function accumulator(::Type{WeightedGyrationTensor}, rs, ws)
    W = eltype(ws)
    R = eltype(eltype(rs))
    T = typeof(zero(W) * zero(R))
    return zeros(MMatrix{3, 3, T})
end

@inline function accumulator(::Type{WeightedRadiusOfGyration}, rs, ws)
    W = eltype(ws)
    R = eltype(eltype(rs))
    T = typeof(zero(W) * zero(R))
    return zero(T)
end

#========================#
#  Particle Coordinates  #
#========================#

const Identity = typeof(identity)

getx(r) = @inbounds(r[1])
gety(r) = @inbounds(r[2])
getz(r) = @inbounds(r[3])

barycenter(f::F, rs) where {F} = sum(f, rs) / length(rs)

function barycenter(::Identity, rs, ws)
    # Alternative implementation:
    # f = ((w, r),) -> w * r
    # R = sum(f, rs)
    R = accumulator(WeightedBarycenter{Identity}, rs, ws)
    @inbounds @simd for i in eachindex(ws, rs)
        w, r = ws[i], rs[i]
        R .= muladd.(w, r, R)
    end
    return SVector(R)
end

function barycenter(f::F, rs, ws) where {F}
    # Alternative implementation:
    # g = ((w, r),) -> w * f(r)
    # R = sum(g, rs)
    R = accumulator(WeightedBarycenter, rs, ws)
    @inbounds @simd for i in eachindex(ws, rs)
        w, r = ws[i], rs[i]
        R = muladd(w, f(r), R)
    end
    return R
end

### Accumulators

@inline function accumulator(::Type{CV}, rs, ws) where {CV <: WeightedBarycenter}
    W = eltype(ws)
    R = eltype(eltype(rs))
    T = typeof(zero(W) * zero(R) / 1)
    return _accumulator(CV, T)
end

@inline _accumulator(::Type{<:WeightedBarycenter}, ::Type{T}) where {T} = zero(T)
@inline function _accumulator(::Type{<:WeightedBarycenter{Identity}}, ::Type{T}) where {T}
    return zeros(MVector{3, T})
end

#====================#
#  Pairwise Kernels  #
#====================#

function pairwise(f::F, r₁, r₂) where {F}
    ξ = accumulator(PairwiseKernel, f, r₁, r₂)
    @inbounds for i in eachindex(r₁)
        @simd for j in eachindex(r₂)
            ξ += f(r₁[i], r₂[j])
        end
    end
    return ξ
end

@inline function accumulator(::Type{PairwiseKernel}, f::F, r₁, r₂) where {F}
    T₁ = eltype(eltype(r₁))
    T₂ = eltype(eltype(r₂))
    T = typeof(f(zero(T₁), zero(T₂)))
    return zero(T)
end

### Common kernels

struct Gaussian{T} <: Function
    μ::T
    σ²::T
end

function Gaussian(μ, σ)
    μ̃, σ̃² = promote(μ, σ^2)
    T = typeof(μ̃)
    return Gaussian{T}(μ̃, σ̃²)
end

mean(f::Gaussian) = f.μ
var(f::Gaussian) = f.σ²

(f::Gaussian)(r) = gaussian(mean(f), var(f), r)

gaussian(μ, σ², r) = exp(-(r - μ)^2 / 2σ²)

#=============#
#  Distances  #
#=============#

distance(r, s) = norm(r - s)

#=============#
#  Utilities  #
#=============#

abstract type DecayingFunction <: Function end

struct RationalDecay{M, N, T} <: DecayingFunction
    d::T
    r::T

    function RationalDecay{M, N}(d, r) where {M, N}
        @assert N > M > 0
        d, r = promote(d₀, r₀)
        return new{M, N, typeof(d)}(d, r)
    end
end

numerator_exponent(::RationalDecay{M}) where {M} = M
denominator_exponent(::RationalDecay{M, N}) where {M, N} = N

origin(f::RationalDecay) = f.d
width(f::RationalDecay) = f.r

RationalDecay(; d₀ = 0.0, r₀ = 1.0) = RationalDecay{6, 12}(d₀, r₀)
RationalDecay{M, N}(; d₀ = 0.0, r₀ = 1.0) where {M, N} = RationalDecay{M, N}(d₀, r₀)

(f::RationalDecay{6, 12})(r) = rational_decay⁶₁₂(r; d₀ = origin(f), r₀ = width(f))
(f::RationalDecay{8, 12})(r) = rational_decay⁸₁₂(r; d₀ = origin(f), r₀ = width(f))
function (f::RationalDecay)(r)
    m = numerator_exponent(f)
    n = denominator_exponent(f)
    return rational_decay(m, n, r; d₀ = origin(f), r₀ = width(f))
end

### More accurate implementation for M = 6, N = 12
@inline function rational_decay⁶₁₂(r; d₀ = 0.0, r₀ = 1.0)
    ρ = ((r - d₀) / r₀)
    if ρ < 0
        return one(ρ)
    end
    return 1 / (1 + ρ^6)
end

### More accurate implementation for M = 8, N = 12
@inline function rational_decay⁸₁₂(r; d₀ = 0.0, r₀ = 1.0)
    ρ = ((r - d₀) / r₀)
    if ρ < 0
        return one(ρ)
    end
    ρ⁴ = ρ^4
    return 1 / (ρ⁴ + 1 / (1 + ρ⁴))
end

@inline function rational_decay(m, n, r; d₀ = 0.0, r₀ = 1.0)
    @assert n > m > 0
    ρ = (r - d₀) / r₀
    if ρ < 0
        return one(ρ)
    elseif isone(ρ)
        return one(ρ) * m / n
    end
    ρᵐ = ρ^m
    if isinf(ρᵐ)
        return zero(ρ)
    end
    return (1 - ρᵐ) / (1 - ρ^n)
end

#===============#
#  Rouse Modes  #
#===============#

function rouse_modeₒ(rs, c)
    x₀ = sum(rs)
    return c * norm(x₀)
end
#
function rouse_mode(rs, c, ws)
    xₚ = sum(((w, r),) -> w * r, zip(ws, rs))
    return c * norm(xₚ)
end

end  # module ReactionCoordinates


text:
"""
This file holds all the includes commands for easier test and dev purpose
"""

include("TypeDef.jl")
include("lossFuns.jl")
include("actFuns.jl")
include("initializations.jl")
include("cnn/includes.jl")
# include("backForProp.jl")
include("parallelBackForProp.jl")
include("assistance.jl")
include("chain.jl")


text:
"""
standard Armijo linesearch to determine step length
default parameters c = τ = 0.5 as in Armijo's paper (1966)
"""
function backtracking_linesearch(X, U, Z, dZ, dU, μ0 = 1, c = 0.5, τ = 0.5)
    μ = μ0
    val = norm(X - U*Z')^2
    m = -(norm(dU)^2 + norm(dZ)^2)
    A = X - U*Z'
    a = [dU*Z', U*dZ', dU*dZ']
    while norm(A + μ*(a[1] + a[2]) - μ^2*a[3])^2 > val + μ*c*m
        μ *= τ
    end
    return μ
end

"""
augments low rank representation with relevant directions from normal cone
"""
function augment_rank(N, U0, Z0, σ)
    V, Σ, P = svd(N)
    r = sum(Σ .>= σ)
    N_U = V[:, 1:r]
    N_Z = P[:, 1:r]*Diagonal(Σ[1:r])
    U = cat(U0, N_U, dims=2)
    Z = cat(Z0, N_Z, dims=2)
    project_to_low_rank!(U, Z, U0, Z0)
    return U, Z
end

"""
finds normal to low rank matrix manifold
"""
normal_direction(U, Z, L_U, L_Z) = (I - U*U')*L_U*(L_Z*(I - Z*pinv(Z'*Z)*Z'))
normal_direction(U, Z, L) = (I - U*U')*L*(I - Z*pinv(Z'*Z)*Z')

"""
finds orthonormal ̄U such that || (U-̄U)Z || is minimal (via gradient descent)
"""
function orthonormalize!(U, Z; ϵ = 1e-8, μ = 0.1, iter_max = 1000)
    iter = 0
    K = U'*U
    A = Matrix{Float64}(I, size(K))
    Ainv = A
    dA = A'*K*A - I
    while norm(dA) > ϵ && iter < iter_max
        A += μ*dA
        Ainv -= μ*Ainv*A*Ainv
        dA = - K * A * (A'*K*A - I)
        iter += 1
    end

    if iter == iter_max
        @warn "Maximum number of iterations exceeded in orthonormalization step"
    end

    U .= U*A
    Z .= Z*Ainv'
end

"""
finds orthonormal U and Z such that || UZ' - U0 Z0'|| is minimal (via gradient
descent)
"""
function project_to_low_rank!(U, Z, U0, Z0; μ = 0.1, ϵ = 1e-8, iter_max = 1000)
    iter = 0
    while iter < iter_max
        dZ = (Z - Z0*(U0'*U))
        dU = -(I - U*U')*U0*(Z0'*Z*pinv(Z'*Z)) # probably better to use a backsolve
        Z .-= μ*dZ
        U .-= μ*dU
        orthonormalize!(U, Z; μ = 0.1)
        if √(norm(dZ)^2 + norm(dU)^2) < ϵ
            break
        end
        iter += 1
    end
    if iter == iter_max
        @warn "Maximum number of iterations exceeded in projection step"
    end
end

"""
finds orthonormal U and Z such that || X - U Z'|| is minimal (via gradient
descent)
"""
function project_to_low_rank!(U, Z, X; μ = 0.1, ϵ = 1e-8, iter_max = 1000)
    iter = 0
    while iter < iter_max
        dZ = (Z - X*U)
        dU = -(I - U*U')*X*Z*pinv(Z'*Z) # probably better to use a backsolve
        Z .-= μ*dZ
        U .-= μ*dU
        orthonormalize!(U, Z; μ = 0.1)
        if √(norm(dZ)^2 + norm(dU)^2) < ϵ
            break
        end
        iter += 1
    end

    if iter == iter_max
        @warn "Maximum number of iterations exceeded in projection step"
    end
end


text:
function change(k)
    if k==1
        return 0
    else
        return 1
    end
end

function mark(r, k)
    if k==1
        putmarker!(r)
    end
end

function draw(r)
    k = 1
    right = 0
    bot = 0
    putmarker!(r)
    while !isborder(r, West)
        move!(r, West)
        k = change(k)
        right+=1
    end
    while !isborder(r, Nord)
        move!(r, Nord)
        k = change(k)
        bot+=1
    end

    mark(r, k)
    while !isborder(r, Ost)
        move!(r, Ost)
        k = change(k)
        mark(r, k)
    end

    while !isborder(r, Sud)
        move!(r, Sud)
        k = change(k)
        mark(r, k)

        while !isborder(r, West)
            move!(r, West)
            k = change(k)
            mark(r, k)
        end

        if !isborder(r, Sud)
            move!(r, Sud)
            k = change(k)
            mark(r, k)
            while !isborder(r, Ost)
                move!(r, Ost)
                k = change(k)
                mark(r, k)
            end
        end
    end

    while !isborder(r, West)
        move!(r, West)
    end
    while !isborder(r, Nord)
        move!(r, Nord)
    end

    for _ = 1:right
        move!(r, Ost)
    end
    for _ = 1:bot
        move!(r, Sud)
    end
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using BinDeps

@BinDeps.setup

libf477 = library_dependency("libf477")

const cpcsrcdir = joinpath(BinDeps.depsdir(libf477),"src","f477")
provides(SimpleBuild,
    (@build_steps begin
        ChangeDirectory(cpcsrcdir)
        `make install`
    end), libf477)

@BinDeps.install Dict([(:libf477, :libf477)])


text:
crawl(::Val{:braces}, args...) = _crawl_dict(args...)

function _crawl_dict(args...)
    return :(JSAST(
        :object,
        $(_crawl_objectpair.(args)...)
    ))
end

# Enable `@js { foo }` syntax
function _crawl_objectpair(sym::Symbol)
    return :(
        JSAST(:objectpair, JSTerminal($(string(sym))), $(crawl(sym)))
    )
end

# Enable `@js { foo=foo }` and `@js { :foo => foo }` syntax
function _crawl_objectpair(expr::Expr)
    if expr.head == :call && expr.args[1] == :(=>)
        lhs, rhs = expr.args[2:3]
        # Ensure that the LHS of the pair is a Symbol literal (which will be
        # wrapped in a QuoteNode in the Expr) or a string.
        if lhs isa QuoteNode
            lhs = string(_unquote(lhs))
        elseif lhs isa AbstractString
            lhs = string(lhs)
        else
            lhs_repr = lhs isa Symbol ? "identifier $(lhs)" : string(lhs)
            error(
                "Invalid key expression in object literal syntax " *
                "(LHS of pair expression must be a Symbol or String, " *
                "got $(lhs_repr))."
            )
        end
    elseif expr.head == :(=)
        lhs, rhs = expr.args
        if !(lhs isa Symbol)
            error(
                "Invalid key expression in object literal syntax " *
                "(LHS of equals expression must be a literal identifier, " *
                "got $(lhs))."
            )
        end
        lhs = string(lhs)
    end
    return :(
        JSAST(:objectpair, JSTerminal($(lhs::String)), $(crawl(rhs)))
    )
end

# Disallow things like `@js { "foo" }`
function _crawl_objectpair(expr::Any)
    error(
        "Invalid object literal item " *
        "(expected pair expression, equals expression, or literal identifier)."
    )
end

function deparse(::Val{:object}, pairs...)
    return jsstring(
        "{",
        join(deparse.(pairs), ", "),
        "}",
    )
end#

function deparse(::Val{:objectpair}, key, value)
    return jsstring(
        deparse(key), ": ",
        deparse(value),
    )
end

"""
Determine whether the tuple expression with args `items` is a namedtuple.

This is used when crawling tuples, which is defined in `arrays.jl`.
"""
function _is_namedtuple(items...)
    if isempty(items)
        return false
    end
    firstelt = first(items)
    return firstelt isa Expr && firstelt.head == :(=)
end

function _crawl_namedtuple(items...)
    return :(JSAST(:object, $(_crawl_namedtuple_arg.(items)...)))
end

function _crawl_namedtuple_arg(item)
    if !(item isa Expr && item.head == :(=))
        error("Invalid namedtuple item (expected = expression, got $(item)).")
    end
    lhs, rhs = item.args
    if !(lhs isa Symbol)
        error("Invalid namedtuple key (expected Symbol, got $(lhs)).")
    end
    return :(JSAST(:objectpair, JSTerminal($(string(lhs))), $(crawl(rhs))))
end


text:
# ____________________________________________________________________
# ironopt
#
function ironopt_read(file)
    a = readdlm(file, '\t', comments=true)
    df = OrderedDict()
    df[:line] = string.(a[:,1])
    df[:transition] = string.(a[:,2])
    df[:wavelength] = fill(NaN, length(df[:line]));
    df[:ul]  = fill(NaN, length(df[:line]));
    df[:aat] = fill(NaN, length(df[:line]));
    df[:wht] = fill(NaN, length(df[:line]));

    ii = findall(typeof.(a[ :,3]) .== Float64)
    df[:ul][ii] = float.(a[ii,3])

    ii = findall(typeof.(        a[ :,4]) .== Float64)
    df[:wavelength][ii] = float.(a[ii,4])

    ii = findall(typeof.( a[ :,5]) .== Float64)
    df[:aat][ii] = float.(a[ii,5])

    ii = findall(typeof.( a[ :,6]) .== Float64)
    df[:wht][ii] = float.(a[ii,6])

    ii = findall(isfinite.(df[:wavelength])  .&  isfinite.(df[:wht]))
    for (key, val) in df
        df[key] = val[ii]
    end

    return df
end

mutable struct ironopt <: AbstractComponent
    file::String
    fwhm::Float64
    L::Vector{Float64}
    norm::Parameter

    function ironopt(file::String, fwhm::Number)
        out = new(file, float(fwhm), Vector{Float64}(), Parameter(1))
        out.norm.low = 0
        return out
    end
end

function prepare!(comp::ironopt, domain::Domain{1})
    df = ironopt_read(comp.file)
    # Drop Balmer lines (they are accounted for in the main QSFit code)
    ii = findall(getindex.(df[:line], Ref(1:2)) .!= "H\$")
    for (key, val) in df
        df[key] = val[ii]
    end

    λ, L = delta_conv_gauss(df[:wavelength], df[:wht], comp.fwhm / 2.355)
    L ./= int_tabulated(λ, L)[1]
    comp.L = Spline1D(λ, L, k=1, bc="zero")(domain[:])
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::ironopt, domain::Domain{1},
                   norm)
    buffer .= norm .* comp.L
end

ironopt_broad( fwhm) = ironopt(qsfit_data() * "/VC2004/TabA1", fwhm)
ironopt_narrow(fwhm) = ironopt(qsfit_data() * "/VC2004/TabA2", fwhm)


text:
import .Tracker
using .Tracker: TrackedReal, TrackedVector, TrackedMatrix, TrackedArray, @grad

_eps(::Type{<:TrackedReal{T}}) where {T} = eps(T)

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inversed{<:ADBijector{<:TrackerAD}}},
    x::Real
)
    return Tracker.data(Tracker.gradient(b, x)[1])
end
function jacobian(
    b::Union{<:ADBijector{<:TrackerAD}, Inversed{<:ADBijector{<:TrackerAD}}},
    x::AbstractVector{<:Real}
)
    # We extract `data` so that we don't return a `Tracked` type
    return Tracker.data(Tracker.jacobian(b, x))
end

# implementations for Shift bijector
function _logabsdetjac_shift(a::TrackedReal, x::Real, ::Val{0})
    return Tracker.param(_logabsdetjac_shift(Tracker.data(a), Tracker.data(x), Val(0)))
end
function _logabsdetjac_shift(a::TrackedReal, x::AbstractVector{<:Real}, ::Val{0})
    return Tracker.param(_logabsdetjac_shift(Tracker.data(a), Tracker.data(x), Val(0)))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractVector{<:Real},
    ::Val{1}
)
    return Tracker.param(_logabsdetjac_shift(Tracker.data(a), Tracker.data(x), Val(1)))
end
function _logabsdetjac_shift(
    a::Union{TrackedReal, TrackedVector{<:Real}},
    x::AbstractMatrix{<:Real},
    ::Val{1}
)
    return Tracker.param(_logabsdetjac_shift(Tracker.data(a), Tracker.data(x), Val(1)))
end

# implementations for Scale bijector
# Adjoints for 0-dim and 1-dim `Scale` using `Real`
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Val{0})
    return Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::Real, ::Val{0})
    return _logabsdetjac_scale(Tracker.data(a), Tracker.data(x), Val(0)), Δ -> (inv(Tracker.data(a)) .* Δ, nothing, nothing)
end
# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Val{0})
    return Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Val{0})
    da = Tracker.data(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, Tracker.data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Val{0})
    return Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val(0))
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Val{0})
    da = Tracker.data(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, Tracker.data(x), Val(0)), Δ -> (transpose(J) * Δ, nothing, nothing)
end
# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    return Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Val{1})
    # ∂ᵢ (∑ⱼ log|aⱼ|) = ∑ⱼ δᵢⱼ ∂ᵢ log|aⱼ|
    #                 = ∂ᵢ log |aᵢ|
    #                 = (1 / aᵢ) ∂ᵢ aᵢ
    #                 = (1 / aᵢ)
    da = Tracker.data(a)
    J = inv.(da)
    return _logabsdetjac_scale(da, Tracker.data(x), Val(1)), Δ -> (J .* Δ, nothing, nothing)
end
function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    return Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val(1))
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Val{1})
    da = Tracker.data(a)
    Jᵀ = repeat(inv.(da), 1, size(x, 2))
    return _logabsdetjac_scale(da, Tracker.data(x), Val(1)), Δ -> (Jᵀ * Δ, nothing, nothing)
end
# TODO: implement analytical gradient for scaling a vector using a matrix
# function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     Tracker.track(_logabsdetjac_scale, a, Tracker.data(x), Val{1})
# end
# @grad function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Val{1})
#     throw
# end

# implementations for Stacked bijector
function logabsdetjac(b::Stacked, x::TrackedMatrix{<:Real})
    return Tracker.collect([logabsdetjac(b, x[:, i]) for i = 1:size(x, 2)])
end
# TODO: implement custom adjoint since we can exploit block-diagonal nature of `Stacked`
function (sb::Stacked)(x::TrackedMatrix{<:Real})
    init = reshape(sb(x[:, 1]), :, 1)
    return Tracker.collect(mapreduce(i -> sb(x[:, i]), hcat, 2:size(x, 2); init = init))
end


text:
module ReactiveMPCorrectionTest

using Test
using ReactiveMP
using Random

using LinearAlgebra

@testset "Correction" begin

    rng = MersenneTwister(1234)

    for n in [ 3, 5, 10 ]
        a = rand(rng)
        A = rand(rng, n, n)

        b = ReactiveMP.correction!(NoCorrection(), a)
        B = ReactiveMP.correction!(NoCorrection(), A)

        @test a == b
        @test a === b
        @test A == B
        @test A === B

        c = ReactiveMP.correction!(TinyCorrection(), ReactiveMP.tiny / 2)
        @test c >= ReactiveMP.tiny

        C = ReactiveMP.correction!(TinyCorrection(), copy(A))

        @test A ≈ C
        @test mapreduce((d) -> d[1] + ReactiveMP.tiny === d[2], &, zip(diag(A), diag(C)))

        D = rand(rng, n, n)
        E = ReactiveMP.correction!(TinyCorrection(), D)
        @test D === E

        v = 1e-10 * rand(rng)

        f = ReactiveMP.correction!(FixedCorrection(v), v / 2)
        @test f >= v

        F = ReactiveMP.correction!(FixedCorrection(v), copy(A))
        @test A ≈ F
        @test mapreduce((d) -> d[1] + v === d[2], &, zip(diag(A), diag(F)))

        G = rand(rng, n, n)
        H = ReactiveMP.correction!(FixedCorrection(v), G)
        @test G === H

        j = ReactiveMP.correction!(ClampEigenValuesCorrection(10.0), 5.0)
        @test j >= 10.0

        J   = ReactiveMP.correction!(ClampEigenValuesCorrection(10.0), copy(A))
        S_J = svd(J)

        @test mapreduce((d) -> d >= 10.0 || d ≈ 10.0, &, S_J.S)

        K = ReactiveMP.correction!(ClampEigenValuesCorrection(1e-12), copy(A))
        @test K ≈ A
    end

end

end


text:
using HttpServer

"""
A `FujiServer` represents the webserver itself. There currently should only ever
be one of these, but this may be changed later. It currently has three fields:

* `after`: An optional block of code that may be run after requests are
    evaluated
* `before`: An optional block of code that may be run before requests are
    evaluated
* `routes`: An array of the routes that the server is currently serving
"""
type FujiServer
    after::Nullable{Function}
    before::Nullable{Function}
    routes::Array{Route,1}
end

"""
    after(action)

Adds a block of code that will be evaluated after each request. Can read the
request and read/modify the response.

# Examples

```julia-repl
julia> after() do req, res
           println(req.resource, " was hit")
       end
(::#1) (generic function with 1 method)
```
"""
after(action::Function) = server.after = action

"""
    before(action)

Adds a block of code that will be evaluated before each request. Can read the
request and read/modify the response.

# Examples

```julia-repl
julia> before() do req, res
           println(req.resource, " will be hit")
       end
(::#1) (generic function with 1 method)
```
"""
before(action::Function) = server.before = action

"""
    route(action, endpoint[, methods])

Adds a route to the Fuji webserver. Assumes that the route is only for GET
requests by default, but passing `methods` will allow you to specify multiple
HTTP methods.

# Examples

```julia-repl
julia> route("/hello") do req, res
           "hi there!"
       end
Fuji.Route(#1, "/hello", String["GET"])
```

```julia-repl
julia> route("/hello", ["GET", "POST"]) do req, res
           "hi there!"
       end
Fuji.Route(#1, "/hello", String["GET", "POST"])
```
"""
function route(action::Function, endpoint::String, methods::Array{String,1}=["GET"])
    route = Route(action, endpoint, methods)

    # ensure that there is only one route with a given endpoint
    before = size(server.routes, 1)
    unroute(route)
    after = size(server.routes, 1)

    if before > after
        warn("There was already a route with the endpoint ", endpoint, ". Removing the previous one now...")
    end

    # add the route to the routes list, effectively turning it on
    push!(server.routes, route)

    route
end

"""
    post(action, endpoint)

Adds a POST route to the Fuji webserver.

# Examples

```julia-repl
julia> post("/hello") do req, res
           "hello"
       end
Fuji.Route(#1, "/hello", String["POST"])
```
"""
post(action::Function, endpoint::String) = route(action, endpoint, ["POST"])

"""
    put(action, endpoint)

Adds a PUT route to the Fuji webserver.

# Examples

```julia-repl
julia> put("/hello") do req, res
           "hello"
       end
Fuji.Route(#1, "/hello", String["PUT"])
```
"""
put(action::Function, endpoint::String) = route(action, endpoint, ["PUT"])

"""
    patch(action, endpoint)

Adds a PATCH route to the Fuji webserver.

# Examples

```julia-repl
julia> patch("/hello") do req, res
           "hello"
       end
Fuji.Route(#1, "/hello", String["PATCH"])
```
"""
patch(action::Function, endpoint::String) = route(action, endpoint, ["PATCH"])

"""
    delete(action, endpoint)

Adds a DELETE route to the Fuji webserver.

# Examples

```julia-repl
julia> delete("/hello") do req, res
           "hello"
       end
Fuji.Route(#1, "/hello", String["DELETE"])
```
"""
delete(action::Function, endpoint::String) = route(action, endpoint, ["DELETE"])

"""
    unroute(route)

Removes a route from the Fuji webserver. If other routes exist with the same
endpoint but with different HTTP methods, those will be preserved.

# Examples

```julia-repl
julia> r = route("/hello") do req, res
           "hello"
       end
Fuji.Route(#1, "/hello", String["GET"])

julia> unroute(r)
0-element Array{Fuji.Route,1}
```
"""
function unroute(route::Route)
    # remove all routes with the same endpoint as the one that was just passed
    filter!(server.routes) do r
        if r.endpoint == route.endpoint
            return length(intersect(r.methods, route.methods)) == 0
        end

        true
    end
end

"""
    start([host, port])

Starts the Fuji webserver so it can begin receiving and interpreting requests.

# Examples

```julia-repl
julia> Fuji.start()
 * Running on http://127.0.0.1:8000/ (Press CTRL+C to quit)
```

```julia-repl
jjulia> Fuji.start(IPv4(127, 0, 0, 1), 8001)
 * Running on http://127.0.0.1:8001/ (Press CTRL+C to quit)
```
"""
function start(host=IPv4(127, 0, 0, 1), port=8000)
    http = HttpHandler() do req, res
        timestamp = Dates.format(now(), "u d, yyyy HH:MM:SS")
        request = FujiRequest(req)
        route_found = false

        for route in server.routes
            if ismatch(req, route)
                request = FujiRequest(req, route)

                if !isnull(server.before)
                    server.before.value(request, res)
                end

                value = route.action(request, res)

                if isa(value, Int64)
                    res.status = value
                elseif isa(value, AbstractString)
                    res.data = value
                end

                route_found = true

                break
            end
        end

        if !route_found
            if !isnull(server.before)
                server.before.value(request, res)
            end

            res.status = 404
        end

        log("[", timestamp, "] ", request.method, " ", request.resource, " -> ", res.status)

        if !isnull(server.after)
            server.after.value(request, res)
        end

        res
    end

    http.events["listen"] = (saddr) -> log(" * Running on http://$saddr/ (Press CTRL+C to quit)")
    web_server = Server(http)

    try
        run(web_server, host=host, port=port)
    catch e
        if isa(e, InterruptException)
            log("Stopping server...")
        end
    end
end


text:
function phi_galerk_sgpe(phi, V, t, dt, lambda, r, g, k, xstep, T, D)
    noise = sqrt(T * dt / xstep) * sampler_complex(size(phi))
    d = length(size(phi))
    R = CartesianIndices(phi)
    out = similar(phi)
    @. out = (phi - 1im * dt * (g * abs.(phi)^2 + r + lambda * V[1]) * phi + noise)
    out = fft(out)
    for I in R
        out[I] = 1.0 / (1.0 + 1im * dt * sum([k[I[i]]^2 for i = 1:d])) * out[I]
    end
    out = ifft(out)
    @. phi = out
end

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
using Test
using LogFixPoint16s

# test for 7-11 fraction bits
LogFixPoint16s.set_nfrac(7)
include("logfixpoint16.jl")

LogFixPoint16s.set_nfrac(8)
include("logfixpoint16.jl")

LogFixPoint16s.set_nfrac(9)
include("logfixpoint16.jl")

LogFixPoint16s.set_nfrac(10)
include("logfixpoint16.jl")

LogFixPoint16s.set_nfrac(11)
include("logfixpoint16.jl")

# check changing rounding mode
LogFixPoint16s.set_nfrac(9)
LogFixPoint16s.set_rounding_mode(:log)
LogFixPoint16s.set_rounding_mode(:lin)
include("logfixpoint16.jl")


text:
## Function so the user does not need to type it every time
# --------------------------------------> For now, resolution of simulation is every 1. Need to add modification in case user wants to go lower
function defModStruct()
    model_def = Dict();

    model_def["NameF"] = []; # Name for the model (scripts will be stored with this name) --> String
    model_def["nStat"] = []; # Nuber of states in the model --> Integer
    model_def["nPar"] = []; # Number of parameters of the mode --> Integer
    model_def["nInp"] = []; # Nuber of stimuly (inducer) in the model --> Integer
    model_def["stName"] = []; # Name for all the states --> Vector of strings
    model_def["parName"] = []; # Name for all the parameters --> Vector of strings
    model_def["inpName"] = []; # Name for all the inputs --> Vector of strings
    model_def["eqns"] = []; # Equations defining the model (left and right hand side) --> Vector of strings
    model_def["Y0eqs"] = []; # Steady state equations --> Vector of strings or empty vector if not used. If some element of the equation requires an experimental value for the calculation, please add exp at the beginning of the state(exemple: Cmrna -> expCmrna). Please do not name anything alp.
    model_def["Y0Sim"] = []; # If the ON simulation is required for the steady state --> true, false or a string saying Yes/yes, No/no. If empty the default would be No. Time-scale for the simulation is assumed in minutes. If this wants to be changed, what should be introduced here is a number for the time conversion (e.g. 1/60 if to convert to days or 60 if to convert to seconds)
    model_def["tols"] = []; # Relative and absolute tolerances (in that order) --> vector of 2 floats (can be left empty, where 1e-6 will be taken as default for both)
    model_def["solver"] = []; # IVP solver to solve the ODEs. If nothing specified, the default will be Tsit5(). For more info check https://diffeq.sciml.ai/v2.0/tutorials/ode_example.html

    return(model_def)
end

## Function with all the checks
# NOTE: some checks have been commented since they considered the option that the user would only introduce ODE and steady state equations and nothing else. Code has been changed to allow definition of variables, prints, ifs, fors, etc.

function checkStruct(model_def)

    # Check taht all the dictionary entries are correct
    entries = ["stName","inpName","eqns","nPar","NameF","nInp","parName","nStat","Y0eqs","Y0Sim","tols","solver"]
    if symdiff(entries,keys(model_def))!=[] && symdiff(entries,keys(model_def))!=["modelpath"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the entries of the dictionary, there is something wrong...")
        println(symdiff(entries,keys(model_def)))
        return
    end

    # Check one to be sure there is no empty entry
    if model_def["NameF"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce a name for the model!")
        return
    elseif model_def["nStat"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the number of states!")
        return
    elseif model_def["nPar"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the number of parameters!")
        return
    elseif model_def["nInp"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the number of stimuli!")
        return
    elseif model_def["stName"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the state names!")
        return
    elseif model_def["parName"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the parameter names!")
        return
    elseif model_def["inpName"] == [] && (model_def["nInp"] != 0 && model_def["nInp"] != [0])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the stimuli names!")
        return
    elseif model_def["eqns"] == []
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, introduce the equations!")
        return
    end

    # Checks to see if the content of the fields is the expected (also consier if single value entries have been given as a number or a vector)
    if ((typeof(model_def["NameF"])!=Array{String,1}) && (typeof(model_def["NameF"])!=String))
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field NameF!")
        return
    elseif (typeof(model_def["nStat"][1])!=Int) && (typeof(model_def["nStat"])!=Int)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field nStat!")
        return
    elseif (typeof(model_def["nPar"][1])!=Int) && (typeof(model_def["nPar"])!=Int)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field nPar!")
        return
    elseif (typeof(model_def["nInp"][1])!=Int) && (typeof(model_def["nInp"])!=Int)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field nInp!")
        return
    elseif (typeof(model_def["stName"])!=Array{String,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field stName!")
        return
    elseif (typeof(model_def["parName"])!=Array{String,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field parName!")
        return
    elseif (typeof(model_def["inpName"])!=Array{String,1}) && (typeof(model_def["inpName"])!=Array{Any,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field inpName!")
        return
    elseif (typeof(model_def["eqns"])!=Array{String,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field eqns!")
        return
    elseif (typeof(model_def["Y0eqs"])!=Array{Any,1}) && (typeof(model_def["Y0eqs"])!=Array{String,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Y0eqs!")
        return
    elseif ((typeof(model_def["solver"])!=Array{String,1}) && (typeof(model_def["solver"])!=String)) && ((model_def["solver"])!=[])
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field solver!")
        return
    elseif model_def["Y0Sim"]!=true && model_def["Y0Sim"]!="Yes" && model_def["Y0Sim"]!="yes" &&
    model_def["Y0Sim"]!= false && model_def["Y0Sim"] != "No" && model_def["Y0Sim"] != "no"&&
    model_def["Y0Sim"]!=[true] && model_def["Y0Sim"]!=["Yes"] && model_def["Y0Sim"]!=["yes"] &&
    model_def["Y0Sim"]!= [false] && model_def["Y0Sim"] != ["No"] && model_def["Y0Sim"] != ["no"]&&
    model_def["Y0Sim"]!= [] && model_def["Y0Sim"] != "" && (typeof(model_def["Y0Sim"])!=Array{Int,1}) &&
    (typeof(model_def["Y0Sim"])!=Array{Float64,1}) && (typeof(model_def["Y0Sim"])!=Array{Float32,1}) &&
    (typeof(model_def["Y0Sim"])!=Int) && (typeof(model_def["Y0Sim"])!=Float64) && (typeof(model_def["Y0Sim"])!=Float32)

        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field Y0Sim!")
        return
    elseif (typeof(model_def["tols"])!=Array{Any,1}) && (typeof(model_def["tols"])!=Array{Float64,1}) && (typeof(model_def["tols"])!=Array{Float32,1})
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field tols!")
        return
    elseif ((typeof(model_def["tols"])==Array{Float64,1}) || (typeof(model_def["tols"])==Array{Float32,1})) && (length(model_def["tols"]) != 2)
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the field tols!")
        return
    end

    if "alp" in model_def["parName"] || "alp" in model_def["inpName"] || "alp" in model_def["stName"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but the name alp is reserved for something else in the model...")
        return
    end

    # Take single elements outside vector (if so) to generalise things after
    if typeof(model_def["NameF"]) == Array{String,1}
        model_def["NameF"] = model_def["NameF"][1];
    end

    if typeof(model_def["nStat"]) == Array{Int,1}
        model_def["nStat"] = model_def["nStat"][1];
    end

    if typeof(model_def["nPar"]) == Array{Int,1}
        model_def["nPar"] = model_def["nPar"][1];
    end

    if typeof(model_def["nInp"]) == Array{Int,1}
        model_def["nInp"] = model_def["nInp"][1];
    end

    if model_def["Y0Sim"]==true || model_def["Y0Sim"]=="Yes" || model_def["Y0Sim"]=="yes" ||
    model_def["Y0Sim"]==[true] || model_def["Y0Sim"]==["Yes"] || model_def["Y0Sim"]==["yes"]
        model_def["Y0Sim"]=true;
    elseif model_def["Y0Sim"]== false || model_def["Y0Sim"] == "No" || model_def["Y0Sim"] == "no"||
    model_def["Y0Sim"]== [false] || model_def["Y0Sim"] == ["No"] || model_def["Y0Sim"] == ["no"]||
    model_def["Y0Sim"]== [] || model_def["Y0Sim"] == ""
        model_def["Y0Sim"]=false;
    elseif typeof(model_def["Y0Sim"]) <: Array
        model_def["Y0Sim"]=model_def["Y0Sim"][1];
    end

    if model_def["tols"] == []
        model_def["tols"] = [1e-6, 1e-6];
    end

    if model_def["solver"] == []
        model_def["solver"] = "Tsit5";
    elseif typeof(model_def["solver"]) == Array{String,1}
        model_def["solver"] = model_def["solver"][1];
    end

    # Check that the user has introduced a correct solver
    if (model_def["solver"] != "AutoTsit5(Rosenbrock23())") && (model_def["solver"] != "AutoVern7(Rodas5())") &&
        (model_def["solver"] != "Tsit5") && (model_def["solver"] != "BS3") && (model_def["solver"] != "Vern7") &&
        (model_def["solver"] != "Rodas4") && (model_def["solver"] != "Rodas5") && (model_def["solver"] != "KenCarp4") &&
        (model_def["solver"] != "TRBDF2") && (model_def["solver"] != "RadauIIA") && (model_def["solver"] != "CVODE_BDF")
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Sorry, but you have introduced a wrong solver!")
        println("The options are AutoTsit5(Rosenbrock23()), AutoVern7(Rodas5()), Tsit5, BS3, Vern7, Rodas4, Rodas5, ")
        println("          KenCarp4, TRBDF2, RadauIIA, CVODE_BDF.")
        println("For more info, please check https://diffeq.sciml.ai/stable/tutorials/ode_example/#ode_example-1")
        return
    end

    # Check tha there are no spaces in determined fields
    if findfirst(" ", join(model_def["stName"])) != nothing
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, do NOT include spaces in the entries of the field stName!")
        return
    end

    if findfirst(" ", join(model_def["inpName"])) != nothing
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, do NOT include spaces in the entries of the field inpName!")
        return
    end

    if findfirst(" ", join(model_def["NameF"])) != nothing
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, do NOT include spaces in the entries of the field NameF!")
        return
    end

    if findfirst(" ", join(model_def["parName"])) != nothing
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, do NOT include spaces in the entries of the field parName!")
        return
    end

    # Check that contents make sense
    if length(model_def["stName"]) != model_def["nStat"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check stName and nStat, they do not make sense")
        return
    end
    if length(model_def["inpName"]) != model_def["nInp"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check inpName and nInp, they do not make sense")
        return
    end

    if length(model_def["parName"]) != model_def["nPar"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check parName and nPar, they do not make sense")
        return
    end

    #     if length(model_def["eqns"]) != model_def["nStat"]
    #         println("-------------------------- Process STOPPED!!! --------------------------")
    #         println("Please, chec eqns and nStat, the number of equations is not the same as the number of states")
    #         return
    #     end

    # ----------------------------------------> Might need to double check this part
    # Check if in the ODEs all the states are in the equations and with a d in front
    # First check that the equations contain an = sign and only 1

    # Check that in the ODEs there is only 1 = sign
    if (findfirst.("=", model_def["eqns"])!=findlast.("=", model_def["eqns"]))
        println("-------------------------- WARNING!!! --------------------------")
        println("Please, check the number of = signs in your OED equations!")
        println("Make sure it is something desired (ifs, fors, prints, etc.) and not that you put 2 equations in one same string or forgot a right-hand side. ")
    end
    if model_def["Y0eqs"] != [] # Avoind the check if no Steady state equations are given
        # Check that in the Steady State Equations there is only 1 equal sign
        if sum((findfirst.("=", model_def["Y0eqs"]).==findlast.("=", model_def["Y0eqs"]))) < model_def["nStat"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, check the number of = signs in your Steady State equations! You might have put 2 equations in one same string!")
            return
        end
        # Check that in the Steady state equations all entries have left and right hand side
        # if sum(occursin.("=", model_def["Y0eqs"])) != length(model_def["Y0eqs"])
        #     println("-------------------------- Process STOPPED!!! --------------------------")
        #     println("Please, check the equations! You might have forgoten a left hand side or add 2 equations in one same string in the Steady State Equation.")
        #     return
        # end
        # Check that in the Steady State equations all the states are given as equations
        ds2 = [model_def["Y0eqs"][i][1:findfirst.("=", model_def["Y0eqs"][i])[1]] for i in 1:length(model_def["Y0eqs"]) if !isnothing(findfirst.("=", model_def["Y0eqs"][i]))]
        if sum(occursin.(model_def["stName"], join(ds2))) !=  model_def["nStat"]
            println("-------------------------- Process STOPPED!!! --------------------------")
            println("Please, check the equations! You might have forgotten a Steady State Equation!.")
            return
        end
    end

    # Check that in the ODEs all entries have left and right hand side
    # if sum(occursin.("=", model_def["eqns"])) != length(model_def["eqns"])
    #     println("-------------------------- Process STOPPED!!! --------------------------")
    #     println("Please, check the equations! You might have forgoten a left hand side of the ODEs or add 2 equations in one same string.")
    #     return
    # end

    # Check that in the ODEs all the states are given as equations
    ds = [model_def["eqns"][i][1:findfirst.("=", model_def["eqns"][i])[1]] for i in 1:length(model_def["eqns"]) if !isnothing(findfirst.("=", model_def["eqns"][i]))];
    if sum(occursin.(string.("d", model_def["stName"]), join(ds))) != model_def["nStat"]
        println("-------------------------- Process STOPPED!!! --------------------------")
        println("Please, check the equations! You might have forgotten an ODE!.")
        println("Remember, ODEs should start with a d and no other equation can.")
        return
    end

    return(model_def)
end

## Function to generate all the necesary Julia scripts for the model simulations

function GenerateModel(model_def)

    model_def = checkStruct(model_def);

    # Packages needed for the simulations
    Head = string("
    using DifferentialEquations
    using OrdinaryDiffEq
    using DiffEqBase
    using Sundials
    using ODEInterfaceDiffEq
        ");

    # Generate ODEs function
    # String containing the equtions, ODEs and others
    tet = Array{Any,1}(undef, length(model_def["eqns"])); # Vector of strings that will contain the eauqtions for the function
    # odec = 1; # Counter for the ODEs to take into account equations that are not ODEs
    for i in 1:length(model_def["eqns"]) # Loop over number of equations
        if replace(model_def["eqns"][i], " "=>"")[1]=='d' # If the equation starts with a d (differential)
            for j in 1:model_def["nStat"]
                eqsg = findfirst("=", replace(model_def["eqns"][i], " "=>""))[1];
                if replace(model_def["eqns"][i], " "=>"")[1:eqsg-1] == string("d", model_def["stName"][j])
                    tet[i] = string("        du[", j, "] = ", model_def["eqns"][i],";\n");
                    break;
                else
                    tet[i] = string("        ", model_def["eqns"][i],";\n");
                end
            end
        else # Other equations
            tet[i] = string("            ", model_def["eqns"][i],";\n");
        end
    end;

    fun1 = string("
        # du -> Derivative
        # u -> State at time t
        # p-> paramter vector
        # t-> time as tuple (init, end)

        function ",join(model_def["NameF"]),"ODE!(du,u,p,t)
            ",join(model_def["stName"], ", ")," = u;

            ",join(model_def["parName"], ", "),", ",join(model_def["inpName"], ", ")," = p;

",join(tet),"

        end
        ");

    # Generate Steady State function
    y0eq = ""; # Vector that will contain the steady state equations (SSE)
    if model_def["Y0eqs"]==[] # If there is no equations, make the function return the same Y0 that the user gives
        y0eq = (string(join(model_def["stName"], ", ")," = I; \n        ",

                "       return(I) \n "));
    else
        tet = Array{Any,1}(undef, length(model_def["Y0eqs"])); # Vector of strings that will contain the eauqtions for the function
        odec = 1 # Counter for the SSEs to take into account equations that are not SSEs
        df = Array{Any,1}(undef, length(model_def["Y0eqs"]));
        for i in 1:length(model_def["Y0eqs"]) # Check for the equations that start with the name of the states (to tell apart from other equations)
            if !isnothing(findfirst.("=", model_def["Y0eqs"][i]))
                df[i] = model_def["Y0eqs"][i][1:findfirst.("=", model_def["Y0eqs"][i])[1]]
            else
                df[i] = "";
            end
        end;
        for i in 1:length(model_def["Y0eqs"]) # Loop over equations
            if sum(occursin.(model_def["stName"], df[i]))!=0 # If the equation looked begins with the name of a state
                tet[i] = string("        alp[", odec, "] = ", model_def["Y0eqs"][i],";\n");
                odec +=1;
            else
                tet[i] = string("        ", model_def["Y0eqs"][i],";\n");
            end
        end

        y0eq = (string(join([string("exp",model_def["stName"][k]) for k in 1:model_def["nStat"]], ", ")," = I; \n        ",

        join(model_def["parName"], ", "),", ",join(model_def["inpName"], ", ")," = p;\n \n \n ",

                "        alp = zeros(",length(model_def["stName"]),"); \n \n",

                join(tet),
                "        return(alp) \n "));

    end

    fun2 = string("
        # p -> parameter vector (plus inducer vector at the end)
        # I -> initial Y0 vector

        function ",join(model_def["NameF"]),"SteadyState(p,I)

            ",y0eq,"

        end
        ");

    # Generate Step Wise simulation function
    ONSim = "";
    if model_def["Y0Sim"]==false
        ONSim = "y0 = y_al;";
    elseif model_def["Y0Sim"]==true
        if model_def["solver"] == "CVODE_BDF"
            ONSim = string("
                prob = ODEProblem(",join(model_def["NameF"]),"ODE!,y_al,(0.0,Float64(24*60-1)),pSS);
                ssv = Sundials.solve(prob, CVODE_BDF(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],");
                y0 = ssv[:,end]; ");
        else
            ONSim = string("
                prob = ODEProblem(",join(model_def["NameF"]),"ODE!,y_al,(0.0,Float64(24*60-1)),pSS);
                ssv = DifferentialEquations.solve(prob, ",join(model_def["solver"]),"(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],");
                y0 = ssv[:,end]; ");
        end
    elseif typeof(model_def["Y0Sim"]) == Float64 || typeof(model_def["Y0Sim"]) == Float32 || typeof(model_def["Y0Sim"]) == Int
        if model_def["solver"] == "CVODE_BDF"
            ONSim = string("
                prob = ODEProblem(",join(model_def["NameF"]),"ODE!,y_al,(0.0,Float64(",model_def["Y0Sim"],"*24*60-1)),pSS);
                ssv = Sundials.solve(prob, CVODE_BDF(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],");
                y0 = ssv[:,end]; ");
        else
            ONSim = string("
                prob = ODEProblem(",join(model_def["NameF"]),"ODE!,y_al,(0.0,Float64(",model_def["Y0Sim"],"*24*60-1)),pSS);
                ssv = DifferentialEquations.solve(prob, ",join(model_def["solver"]),"(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],");
                y0 = ssv[:,end]; ");
        end
    end

    if model_def["solver"] == "CVOED_BDF"
        soso = "Sundials.";
    else
        soso = "DifferentialEquations.";
    end

    fun3 = string("
        # ts -> time vector going from t=0 to t=end every 1
        # p -> parameter vector
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # pre -> Vector with the inducers in the ON. It can be empty

        function ",join(model_def["NameF"]),"_solvecoupledODE(ts, p, sp, inputs, ivss, pre=[])

            maxtime = length(ts);
            Nsp = length(sp);
            Nevents = length(sp)-1;
            Neq = ",join(model_def["nStat"]),";
            if pre != []
                pSS = vcat(p, pre);
            else
                pSS = p;
            end

            final = zeros(maxtime, Neq);

            y_al = ",join(model_def["NameF"]),"SteadyState(pSS,ivss) # Calculation of initial guesses for steady state

            ",join(ONSim),"
            initialV = y0;
            i = 1;

            for q in collect(1:Nevents)
                lts = length(ts[(sp[q]+1):sp[q+1]+1]);  # General way to define the number of elements in each event series
                Tevent = ts[(sp[q]+1):sp[q+1]+1];  # General way to extract the times of each event
                I = inputs[i:(i+(",model_def["nInp"]-1,"))];
                pSte = vcat(p, I);

                if q == 1
                    prob = ODEProblem(",join(model_def["NameF"]),"ODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                    part1 = ",join(soso),"solve(prob, ",join(model_def["solver"]),"(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],",saveat=1);
                else
                    prob = ODEProblem(",join(model_def["NameF"]),"ODE!,initialV,(ts[(sp[q]+1)],ts[sp[q+1]+1]),pSte);
                    part1 = ",join(soso),"solve(prob, ",join(model_def["solver"]),"(),reltol=",model_def["tols"][1],",abstol=",model_def["tols"][2],",saveat=1);
                end

                initialV = part1[lts];
                i+=",model_def["nInp"],";

                for d in collect((sp[q]+1):(sp[q]+lts))
                    final[d,:] = part1[d-sp[q]];
                end
            end
            return(final)

        end
        ");

    # Function to simulate the system as many times as parameter draws are given
    fun4 = string("
        # ts -> time vector going from t=0 to t=end every 1
        # pD -> Matrix containing all the draws for the parameters (can be a single vector if nly 1 draw of parameters is considered)
        # sp -> vector with switching times (t=0 and t=end needs to be included). Length of vector will be length(inducer steps)+1
        # inputs -> Vecor with the inducers (if more than 1 inducer in the system, the shape will be a1,b1,c1,..., a2,b2,c2,... aend,bend,cend,... where each leter represents a different inducer)
        # ivss -> Initial Y0 vector
        # samps -> Vector of Time points that want to be extracted from the simulations (it considers the initial time point as 0)
        # pre -> Vector with the inducers in the ON. It can be empty

        function ",join(model_def["NameF"]),"_SolveAll(ts, pD, sp, inputs, ivss, samps, pre=[])

            if length(size(pD)) == 1
                pD = reshape(pD,size(pD)[1],1);
            end

            if size(pD)[2] != ",model_def["nPar"],"
                pD = pD';
            end

            if length(ivss)/",model_def["nStat"]," > 1
                if size(ivss)[2] != ",model_def["nStat"],"
                    ivss = ivss';
                end
            end

            AllSolTest = zeros(length(samps), ",model_def["nStat"],", length(pD[:,1]));

            for drawInd in collect(1:length(pD[:,1]))
                p = pD[drawInd,:];
                if length(ivss)/",model_def["nStat"]," > 1
                    ivss2 = ivss[drawInd,:];
                else
                    ivss2 = ivss;
                end
                temp = ",join(model_def["NameF"]),"_solvecoupledODE(ts, p, sp, inputs, ivss2, pre);
                AllSolTest[:,:,drawInd] = temp[convert.(Int,samps.+1),:];
            end

            return(AllSolTest);
        end
        ");

    finscr = string(Head, fun1, fun2, fun3, fun4);

    fidi = pwd(); #@__DIR__;

    if !isdir(string(fidi,"\\ModelsFunctions"))
        mkdir(string(fidi,"\\ModelsFunctions"));
    end

    open(string(fidi,"\\ModelsFunctions\\", model_def["NameF"], "_Model.jl"), "w") do io
       write(io, finscr);
    end;

    include(string(fidi,"\\ModelsFunctions\\", model_def["NameF"], "_Model.jl"))

    model_def["modelpath"] = string(fidi,"\\ModelsFunctions\\", model_def["NameF"], "_Model.jl");

    println("")
    println("----------------------------------------- MODEL GENERATION -----------------------------------------")
    println("The model has been generated in the directory: ")
    println(string("                 ", model_def["modelpath"]))
    println("--------------------------------------------------------------------------------------")
    println("")

    return(model_def)

end


text:
@enum Statistics fermion boson

statistics_sign(stat::Statistics)::Int64 = stat == fermion ? -1 : 1

PyObject(stat::Statistics)::PyObject = (stat == fermion ? "F" : "B")


text:
using SciMLTutorials
SciMLTutorials.weave_file("Testing","test.jmd")


text:
"""
    AMHCopula{d,T}

Fields:
  - θ::Real - parameter

Constructor

    AMHCopula(d, θ)

The [AMH](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Most_important_Archimedean_copulas) copula in dimension ``d`` is parameterized by ``\\theta \\in [0,1)``. It is an Archimedean copula with generator :

```math
\\phi(t) = 1 - \\frac{1-\\theta}{e^{-t}-\\theta}
```
"""
struct AMHCopula{d,T} <: ArchimedeanCopula{d}
    θ::T
end
AMHCopula(d,θ) = (0 <= θ <= 1) ? AMHCopula{d,typeof(θ)}(θ) : @error "Theta must be between 0 and 1"
ϕ(  C::AMHCopula,t) = (1-C.θ)/(exp(t)-C.θ)
ϕ⁻¹(  C::AMHCopula,t) = log(C.θ + (1-C.θ)/t)

τ(C::AMHCopula) = 1 - 2(C.θ+(1-C.θ)^2*log(1-C.θ))/(3C.θ^2) # no closed form inverse...
function τ⁻¹(::Type{AMHCopula},τ)
    if τ == zero(τ)
        return τ
    end
    if τ > 1/3
        @warn "AMHCopula cannot handle kendall tau's greater than 1/3. We capped it to 1/3."
        return 1
    end
    return Roots.fzero(θ -> 1 - 2(θ+(1-θ)^2*log(1-θ))/(3θ^2) - τ,0.5)
end

radial_dist(C::AMHCopula) = 1 + Distributions.Geometric(1-C.θ)


text:
using Test
using JLD2
using Statistics
using LinearAlgebra
using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation.EnsembleKalmanInversions: iterate!, FullEnsembleDistribution, NaNResampler, resample!
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize

data_path = "convective_adjustment_test.jld2"
Nensemble = 3

@testset "EnsembleKalmanInversions tests" begin
    #####
    ##### Build InverseProblem
    #####

    observation = SyntheticObservations(data_path, field_names=(:u, :v, :b))

    Nz = observation.grid.Nz
    Hz = observation.grid.Hz
    Lz = observation.grid.Lz
    Δt = observation.metadata.parameters.Δt

    Qᵘ, Qᵇ, N², f = [zeros(Nensemble, 1) for i = 1:4]
    Qᵘ[:, 1] .= observation.metadata.parameters.Qᵘ
    Qᵇ[:, 1] .= observation.metadata.parameters.Qᵇ
    N²[:, 1] .= observation.metadata.parameters.N²
    f[:, 1] .= observation.metadata.coriolis.f

    file = jldopen(observation.path)
    closure = file["serialized/closure"]
    close(file)

    column_ensemble_size = ColumnEnsembleSize(Nz=Nz, ensemble=(Nensemble, 1), Hz=Hz)
    ensemble_grid = RectilinearGrid(size = column_ensemble_size, topology = (Flat, Flat, Bounded), z = (-Lz, 0))

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
                                                 tracers = (:b,),
                                                 buoyancy = BuoyancyTracer(),
                                                 boundary_conditions = (; u=u_bcs, b=b_bcs),
                                                 coriolis = FPlane.(f),
                                                 closure = [deepcopy(closure) for i = 1:Nensemble, j=1:1])

    ensemble_simulation = Simulation(ensemble_model; Δt=Δt, stop_time=observation.times[end])

    lb, ub = [0.9, 1.1]
    priors = (
        convective_κz = ConstrainedNormal(0.0, 1.0, lb, ub),
        background_κz = ConstrainedNormal(0.0, 1.0, 5e-5, 2e-4)
    )

    Nparams = length(priors)

    free_parameters = FreeParameters(priors)

    calibration = InverseProblem(observation, ensemble_simulation, free_parameters)

    #####
    ##### Test EKI
    #####

    eki = EnsembleKalmanInversion(calibration; noise_covariance = 0.01)

    iterations = 5
    iterate!(eki; iterations = iterations, show_progress = false)

    @test length(eki.iteration_summaries) == iterations + 1
    @test eki.iteration == iterations

    θθ_std_arr = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...))
    @test size(θθ_std_arr) == (Nparams, iterations + 1)

    # Test that parameters stay within bounds
    parameters = eki.iteration_summaries[0].parameters
    convective_κzs = getproperty.(parameters, :convective_κz)
    @test all(convective_κzs .> lb)
    @test all(convective_κzs .< ub)

    # Test that parameters change
    @test convective_κzs[1] != convective_κzs[2]

    iterate!(eki; iterations = 1, show_progress = false)

    @test length(eki.iteration_summaries) == iterations + 2
    @test eki.iteration == iterations + 1

    #####
    ##### Test Resampler
    #####

    resampler = NaNResampler(; abort_fraction=1.0, distribution=FullEnsembleDistribution())

    θ = rand(Nparams, Nensemble)
    p1 = deepcopy(θ[:, 1])
    p2 = deepcopy(θ[:, 2])

    # Fake a forward map output with NaNs
    G = eki.inverting_forward_map(θ)
    G[:, 2] .= NaN

    resample!(resampler, G, θ, eki)

    @test any(isnan.(G)) == false
    @test θ[:, 1] == p1
    @test θ[:, 2] != p2

end


text:
using Logging
using Test
import TrajectoryOptimization: InnerLoop, print_header, print_row

solver = Solver(Dynamics.pendulum...)
solver.opts.verbose = true
logger = TrajectoryOptimization.default_logger(solver)
@test_logs (:info, "hi") @info "hi"

print_header(logger,InnerLoop)
with_logger(logger) do
    @logmsg InnerLoop :iter value=1
    @logmsg InnerLoop "stuff"
    @logmsg InnerLoop :newcol value=10.2
end
print_row(logger,InnerLoop)

@info "hi"

logger.leveldata[InnerLoop]


text:
function write_netcdf2d_reg(fname,var,lon,lat,varname,filenc)

(nsx,)=size(lon)
(nsy,)=size(lat)
(ns,ns,nt)=size(var)

mode=NC_NETCDF4

ncin= NetCDF.open(filenc)
lonatt= ncin.vars["lon"].atts
latatt= ncin.vars["lat"].atts
timeatt= ncin.vars["time"].atts
varatt= ncin.vars[varname].atts
tim=NetCDF.readvar(ncin,"time");
NetCDF.close(ncin)

#ncgetatt(filenc, "global", gatts)

varattr=Dict()
for k in collect(keys(varatt))
   if((k!="_FillValue")&&(k!="add_offset")&&(k!="scale_factor"))
       varattr[k]=varatt[k]
   end
end

isfile(fname) ? rm(fname) : nothing
nccreate( fname, varname , "lon" , nsx, "lat", nsy, "time", tim[1:nt] ,timeatt, atts=varattr,mode=mode,t=NC_FLOAT);
nccreate( fname, "lon" , "lon" , nsx,  atts=lonatt,mode=mode,t=NC_FLOAT);
nccreate( fname, "lat" , "lat" , nsy,  atts=latatt,mode=mode,t=NC_FLOAT);

ncwrite(var,fname,varname)
ncwrite(lon,fname,"lon")
ncwrite(lat,fname,"lat")

ncclose(fname)
end


text:
@testset "SlidingSum2LP: $(fct_type), $(array_size) items, $(slide_length) length, $(T)" for fct_type in ["vector of variables", "vector affine function"], array_size in [3, 4], slide_length in [1, 2], T in [Int, Float64]
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.SlidingSum2LP{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.SlidingSum{T},
    )

    if T == Int
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
    elseif T == Float64
        x_array = MOI.add_variables(model, array_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x_array)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x_array)
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.SlidingSum(T(2), T(5), slide_length, array_size))

    for i in 1:array_size
        @test MOI.is_valid(model, x_array[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.SlidingSum{T}}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Count{MOI.EqualTo{T}}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.Interval{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.Interval{T}}()) == array_size - slide_length

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.Interval{T}}()) == bridge.cons
    end

    @testset "Constraints" begin
        @test length(bridge.cons) == array_size - slide_length

        for i in 1:(array_size - slide_length)
            @test MOI.is_valid(model, bridge.cons[i])
            @test MOI.get(model, MOI.ConstraintSet(), bridge.cons[i]) == MOI.Interval(T(2), T(5))

            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons[i])
            @test length(f.terms) == slide_length
            @test f.constant === zero(T)

            for j in 1:slide_length
                t = f.terms[j]
                @test t.coefficient === one(T)
                @test t.variable == x_array[i + j - 1]
            end
        end
    end
end


text:
function source(settings, data, ac, shielding, n_t, idx_src, input_src)

    # Extract inputs
    TS = input_src[idx_src["TS"][1]:idx_src["TS"][2]]

    # Get type of input vector
    T = eltype(input_src)

    # Initialize source mean-square acoustic pressure
    msap_source = zeros(T, (n_t, settings.N_f))

    if settings.fan_inlet
        msap_fan_inlet = fan(settings, data, ac, n_t, shielding, idx_src, input_src, "fan_inlet")
        msap_source = msap_source .+ msap_fan_inlet
    end
    if settings.fan_discharge
        msap_fan_discharge = fan(settings, data, ac, n_t, shielding, idx_src, input_src, "fan_discharge")
        msap_source = msap_source .+ msap_fan_discharge
    end

    if settings.core
        msap_core = core(settings, data, ac, n_t, idx_src, input_src)
        if settings.suppression && settings.case_name in ["nasa_stca_standard", "stca_enginedesign_standard"]
            msap_core[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_core)[findall(TS.*ones(1, settings.N_f).>0.8)]
        end
        msap_source = msap_source .+ msap_core
    end

    if settings.jet_mixing
        msap_jet_mixing = jet_mixing(settings, data, ac, n_t, idx_src, input_src)
        if settings.suppression && settings.case_name in ["nasa_stca_standard", "stca_enginedesign_standard"]
            msap_jet_mixing[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_jet_mixing)[findall(TS.*ones(1, settings.N_f).>0.8)]
        end
        msap_source = msap_source .+ msap_jet_mixing
    end
    if settings.jet_shock
        msap_jet_shock = jet_shock(settings, data, ac, n_t, idx_src, input_src)
        if settings.suppression && settings.case_name in ["nasa_stca_standard", "stca_enginedesign_standard"]
            msap_jet_shock[findall(TS.*ones(1, settings.N_f).>0.8)] = (10. ^(-2.3 / 10.) * msap_jet_shock)[findall(TS.*ones(1, settings.N_f).>0.8)]
        end
        msap_source = msap_source .+ msap_jet_shock
    end

    if settings.airframe
        msap_af = airframe(settings, data, ac, n_t, idx_src, input_src)
        msap_source = msap_source .+ msap_af
    end

    return msap_source
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
import Base.conj
export conj
export sign, curvature, monotonicity, evaluate, conic_form!
struct ConjugateAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function ConjugateAtom(x::AbstractExpr)
        children = (x,)
        return new(:conj, hash(children), children, (x.size[1], x.size[2]))
    end
end

function sign(x::ConjugateAtom)
    return sign(x.children[1])
end

function monotonicity(x::ConjugateAtom)
    return (Nondecreasing(),)
end

function curvature(x::ConjugateAtom)
    return ConstVexity()
end

function evaluate(x::ConjugateAtom)
    return conj(evaluate(x.children[1]))
end

function conic_form!(x::ConjugateAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        objective = conic_form!(x.children[1], unique_conic_forms)
        for var in keys(objective)
            x1 = conj(objective[var][1])
            x2 = conj(objective[var][2])
            objective[var] = (x1,x2)
        end
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

conj(x::AbstractExpr) = ConjugateAtom(x)
conj(x::Constant) = Constant(conj(x.value))


text:
@doc raw"""
# Lorenz Attractor

"""
module LorenzAttractor

    using GeometricEquations

    export lorenz_attractor_ode, plot_lorenz_attractor

    Δt = 0.01
    nt = 1000

    const q₀ = [1., 1., 1.]
    const default_params = (σ = 10., ρ = 28., β = 8/3)

    function lorenz_attractor_v(t, x, v, params=default_params)
        σ, ρ, β = params
        v[1] = σ * (x[2] - x[1])
        v[2] = x[1] * (ρ - x[3]) - x[2]
        v[3] = x[1] * x[2] - β * x[3]
        nothing
    end

    function lorenz_attractor_ode(q₀=q₀)
        ODE(lorenz_attractor_v, q₀)
    end

end


text:
using Test
using TestSetExtensions

@testset ExtendedTestSet "ThorinDistributions tests" begin
    @includetests ARGS
end

# @testset verbose=true "ThorinDistributions.jl" begin
#     include("kStats.jl")
#     include("ThorinMoments.jl")
#     include("UnivariateGammaConvolution.jl")
#     include("HankelMatrices.jl")
# end


text:
__precompile__()

module SolarSystemLib

include("DateUtils.jl")

export eccfaccor,
    declination,
    sunsethourlyangle,
    middayhour,
    sunsethour,
    sunrisehour

# ARTICLES
#  - http://old.sevierriver.org/et/t110.html
#
#

@doc """
Correction of eccentricity factor
"""
eccfaccor(julia_year::Signed) = 1+(0.033*cos(180/pi*julia_year/365))
eccfaccor(date::Date) = eccfaccor(DateUtils.julian_year(date))

@doc """
Solar Declination
"""
declination(julian_year::Signed) = 23.45 * pi / 180 * (sin((2 * pi / 365) * (julian_year - 80)))
declination(date::Date) = declination(DateUtils.julian_year(date))

@doc """

"""
sunsethourlyangle(lat_rad::AbstractFloat, declination::AbstractFloat) = acos(-tan(lat_rad) * tan(declination))

@doc """
Middle Day Hour
"""
middayhour(hourly_angle::AbstractFloat) = hourly_angle * (180 / pi / 15)

@doc """
Sunset hour
"""
sunsethour(hourly_angle::AbstractFloat) = middayhour(hourly_angle) + 12

@doc """
Sunrise hour
"""
sunrisehour(hourly_angle::AbstractFloat) = 12 - middayhour(hourly_angle)

#@doc """
#
#References:
#    - http://www.ufvjm.edu.br/disciplinas/agr049/files/2014/04/RADIA%C3%87%C3%83O-SOLAR-E-TERRESTRE-1-2-2014-.pdf
#
#"""

end # module


text:
module KNITRO

const _DEPS_FILE = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(_DEPS_FILE)
    include(_DEPS_FILE)
else
    error("KNITRO.jl not properly installed. Please run `] build KNITRO`")
end

using Libdl, SparseArrays
import Base: show

const KNITRO_VERSION = if libknitro == "julia_registryci_automerge"
    VersionNumber(11, 0, 0) # Fake a valid version for AutoMerge
else
    len = 15
    out = zeros(Cchar,len)
    ccall((:KTR_get_release, libknitro), Any, (Cint, Ptr{Cchar}), len, out)
    res = String(strip(String(convert(Vector{UInt8},out)), '\0'))
    VersionNumber(split(res, " ")[2])
end

if KNITRO_VERSION < v"11.0"
    error("You have installed version $KNITRO_VERSION of Artelys Knitro, which is not supported
    by KNITRO.jl. We require a Knitro version greater than 11.0.
    ")
end

include("kn_common.jl")
include("kn_env.jl")

include("kn_model.jl")
include("kn_defines.jl")
include("kn_params.jl")
include("kn_variables.jl")
include("kn_attributes.jl")
include("kn_constraints.jl")
include("kn_residuals.jl")
include("kn_solve.jl")
include("kn_callbacks.jl")

# the MathOptInterface wrapper works only with the new API
include("MOI_wrapper.jl")

end


text:
# # Wealth distribution model

# This model is a simple agent-based economy that is modelled according
# to the work of [Dragulescu *et al.*](https://arxiv.org/abs/cond-mat/0211175).
# This work introduces statistical mechanics concepts to study wealth distributions.
# What we show here is also referred to as "Boltzmann wealth distribution" model.

# This model has a version with and without space.
# The rules of the space-less game are quite simple:
# 1. There is a pre-determined number of agents.
# 2. All agents start with one unit of wealth.
# 3. At every step an agent gives 1 unit of wealth (if they have it) to some other agent.

# Even though this rule-set is simple, it can still recreate the basic
# properties of wealth distributions, e.g. power-laws distributions.

# ## Core structures: space-less
# We start by defining the Agent type and initializing the model.
using Agents
mutable struct WealthAgent <: AbstractAgent
    id::Int
    wealth::Int
end

# Notice that this agent does not have a `pos` field. That is okay, because
# there is no space structure to this example.
# We can also make a very simple [`AgentBasedModel`](@ref) for our model.

function wealth_model(; numagents = 100, initwealth = 1)
    model = ABM(WealthAgent, scheduler = random_activation)
    for i in 1:numagents
        add_agent!(model, initwealth)
    end
    return model
end

model = wealth_model()

# The next step is to define the agent step function
function agent_step!(agent, model)
    agent.wealth == 0 && return # do nothing
    ragent = random_agent(model)
    agent.wealth -= 1
    ragent.wealth += 1
end
nothing # hide

# We use `random_agent` as a convenient way to just grab a second agent.
# (this may return the same agent as `agent`, but we don't care in the long run)

# ## Running the space-less model
# Let's do some data collection, running a large model for a lot of time
N = 5
M = 2000
adata = [:wealth]
model = wealth_model(numagents = M)
data, _ = run!(model, agent_step!, N; adata)
data[(end - 20):end, :]

# What we mostly care about is the distribution of wealth,
# which we can obtain for example by doing the following query:

wealths = filter(x -> x.step == N - 1, data)[!, :wealth]

# and then we can make a histogram of the result.
# With a simple visualization we immediately see the power-law distribution:

using UnicodePlots
UnicodePlots.histogram(wealths)

# ## Core structures: with space
# We now expand this model to (in this case) a 2D grid. The rules are the same
# but agents exchange wealth only with their neighbors.

# It is also available from the `Models` module as [`Models.wealth_distribution`](@ref).

# We therefore have to add a `pos` field as the second field of the agents:

mutable struct WealthInSpace <: AbstractAgent
    id::Int
    pos::NTuple{2,Int}
    wealth::Int
end

function wealth_model_2D(; dims = (25, 25), wealth = 1, M = 1000)
    space = GridSpace(dims, periodic = true)
    model = ABM(WealthInSpace, space; scheduler = random_activation)
    for i in 1:M # add agents in random positions
        add_agent!(model, wealth)
    end
    return model
end

model2D = wealth_model_2D()

# The agent actions are a just a bit more complicated in this example.
# Now the agents can only give wealth to agents that exist on the same or
# neighboring positions (their "neighbors").

function agent_step_2d!(agent, model)
    agent.wealth == 0 && return # do nothing
    neighboring_positions = collect(nearby_positions(agent.pos, model))
    push!(neighboring_positions, agent.pos) # also consider current position
    rpos = rand(neighboring_positions) # the position that we will exchange with
    available_ids = ids_in_position(rpos, model)
    if length(available_ids) > 0
        random_neighbor_agent = model[rand(available_ids)]
        agent.wealth -= 1
        random_neighbor_agent.wealth += 1
    end
end
nothing # hide

# ## Running the model with space
using Random # hide
Random.seed!(5) # hide
init_wealth = 4
model = wealth_model_2D(; wealth = init_wealth)
adata = [:wealth, :pos]
data, _ = run!(model, agent_step!, 10; adata = adata, when = [1, 5, 9])
data[(end - 20):end, :]

# Okay, now we want to get the 2D spatial wealth distribution of the model.
# That is actually straightforward:
using Plots
gr() # hide

function wealth_distr(data, model, n)
    W = zeros(Int, size(model.space))
    for row in eachrow(filter(r -> r.step == n, data)) # iterate over rows at a specific step
        W[row.pos...] += row.wealth
    end
    return W
end

W1 = wealth_distr(data, model2D, 1)
Plots.heatmap(W1)

#

W5 = wealth_distr(data, model2D, 5)
Plots.heatmap(W5)

#

W10 = wealth_distr(data, model2D, 9)
Plots.heatmap(W10)

# What we see is that wealth gets more and more localized.


text:
"""
    AbstractModel{S}

Abstract type for any model over states of type `S`.
"""
abstract type AbstractModel{S} end

"""
    state_type(arg)

Returns the type of the hidden state in `arg`.
Supported argument types:
* Any subtype of AbstractModel
* Any subtype of AbstractFilteringProblem
"""
function state_type(arg) end

"""
    obs_type(arg)

Returns the type of the observed state in `arg`.
Supported argument types:
* Any subtype of ObservationModel
* Any subtype of AbstractFilteringProblem
"""
function obs_type(arg) end

"""
    time_type(model::AbstractModel)

Returns the time type of `model`.
"""
function time_type(model::AbstractModel) end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters