text:
export ent_spec, ent_S, entropy
"""
SchmidtMatrix{Tm <: Number, Ta <: SubArray, Tb <: SubArray, Ti <: Integer}

Struct that is used to construct the Schmidt matrix:

`|Ψ⟩ = ∑ᵢⱼ Mᵢᵢ|ψᵢ⟩ ⊗ |ψⱼ⟩`

Properties:
-----------
- `M`   : Schidt matrix.
- `A`   : View of region A.
- `B`   : View of region B.
- `base`: Base.
"""
struct SchmidtMatrix{Tm <: Number, Ta <: SubArray, Tb <: SubArray, Ti <: Integer}
    M::Matrix{Tm}
    A::Ta
    B::Tb
    base::Ti
end

Matrix(S::SchmidtMatrix) = S.M

"""
    SchmidtMatrix(T::DataType, b::AbstractBasis, Ainds::AbstractVector{Ta}) where Ta <: Integer

Construction of `SchmidtMatrix`.
"""
function SchmidtMatrix(T::DataType, b::AbstractBasis, Ainds::AbstractVector{Ta}) where Ta <: Integer
    L = length(b)
    B = base(b)
    Binds = Vector{Ta}(undef, L-length(Ainds))
    P = 1
    for i in range(one(Ta), stop=convert(Ta, L))
        if !in(i, Ainds)
            Binds[P] = i
            P += 1
        end
    end
    M = zeros(T, B^length(Ainds), B^length(Binds))
    SchmidtMatrix(M, view(digits(b), Ainds), view(digits(b), Binds), B)
end

function addto!(S::SchmidtMatrix, val::Number)
    ia = index(S.A, base=S.base)
    ib = index(S.B, base=S.base)
    S.M[ia, ib] += val
end

"""
    Vector(v::AbstractVector, b::AbstractBasis)

Convert a state in the symmetry sector to the original state in tensor-product Hilbert space.
"""
Vector(v::AbstractVector, b::AbstractBasis) = reshape(schmidt(v, 1:L, b), :)

"""
    ent_spec(v::AbstractVector, Aind::AbstractVector{<:Integer}, b::AbstractBasis)

Conpute the entanglement spectrum.
"""
ent_spec(v::AbstractVector, Aind::AbstractVector{<:Integer}, b::AbstractBasis) = svdvals(schmidt(v, Aind, b))

"""
    entropy(s::AbstractVector{<:Real}; α::Real=1, cutoff::Real=1e-20)

Compute the entropy of Schmidt values.

Inputs:
-------
- `s`     : Schmidt values.
- `α`     : Renyi index.
- `cutoff`: Cutoff of the Schmidt values.

Outputs:
--------
- `S`: Entanglement entropy.
"""
function entropy(s::AbstractVector{<:Real}; α::Real=1, cutoff::Real=1e-20)
    if isone(α)
        shannon_entropy(s, cutoff=cutoff)
    elseif iszero(α)
        renyi_zero_entropy(s, cutoff=cutoff)
    else
        renyi_entropy(s, α)
    end
end

"""
Compute Shannon entropy
"""
function shannon_entropy(s::AbstractVector{<:Real}; cutoff::Real=1e-20)
    ent = 0.0
    for si in s
        if si > cutoff
            ent -= si * log(si)
        else
            break
        end
    end
    ent
end

"""
Compute Renyi-0 entropy, which is thenumber of non-zero Schmidt values.
"""
function renyi_zero_entropy(s::AbstractVector{<:Real}; cutoff::Real=1e-20)
    N = 0
    for si in s
        if si > cutoff
            N += 1
        else
            break
        end
    end
    N
end

"""
Compute general Renyi entropy
"""
renyi_entropy(s::AbstractVector{<:Real}, α::Real) = log(sum(s.^α)) / (1-α)

"""
    ent_S(v::AbstractVector, Aind::AbstractVector{<:Integer}, b::AbstractBasis; α::Real=1, cutoff::Real=1e-20)

Conpute the entanglement entropy of a state.
"""
function ent_S(v::AbstractVector, Aind::AbstractVector{<:Integer}, b::AbstractBasis; α::Real=1, cutoff::Real=1e-20)
    s = ent_spec(v, Aind, b) .^ 2
    entropy(s, α=α, cutoff=cutoff)
end


text:
promote_type(::Type{Float64}, ::Type{Union{}}) = Float64
promote_type(::Type{Float32}, ::Type{Union{}}) = Float32
promote_type(::Type{Float16}, ::Type{Union{}}) = Float16

promote_type(::Type{T}, ::Type{Union{}}) where {T<:MarkableFloat64} = T
promote_type(::Type{T}, ::Type{Union{}}) where {T<:MarkableFloat32} = T
promote_type(::Type{T}, ::Type{Union{}}) where {T<:MarkableFloat16} = T

promote_type(::Type{Float64}, ::Type{T}) where {T<:MarkableFloat64} = T
promote_type(::Type{MarkableF640}, ::Type{MarkableF641}) = MarkableF641
promote_type(::Type{MarkableF640}, ::Type{MarkableF642}) = MarkableF642
promote_type(::Type{MarkableF640}, ::Type{MarkableF643}) = MarkableF643
promote_type(::Type{MarkableF641}, ::Type{MarkableF642}) = MarkableF642
promote_type(::Type{MarkableF641}, ::Type{MarkableF643}) = MarkableF643
promote_type(::Type{MarkableF642}, ::Type{MarkableF643}) = MarkableF643

promote_type(::Type{Float32}, ::Type{T}) where {T<:MarkableFloat32} = T
promote_type(::Type{MarkableF320}, ::Type{MarkableF321}) = MarkableF321
promote_type(::Type{MarkableF320}, ::Type{MarkableF322}) = MarkableF322
promote_type(::Type{MarkableF320}, ::Type{MarkableF323}) = MarkableF323
promote_type(::Type{MarkableF321}, ::Type{MarkableF322}) = MarkableF322
promote_type(::Type{MarkableF321}, ::Type{MarkableF323}) = MarkableF323
promote_type(::Type{MarkableF322}, ::Type{MarkableF323}) = MarkableF323

promote_type(::Type{Float16}, ::Type{T}) where {T<:MarkableFloat16} = T
promote_type(::Type{MarkableF160}, ::Type{MarkableF161}) = MarkableF161
promote_type(::Type{MarkableF160}, ::Type{MarkableF162}) = MarkableF162
promote_type(::Type{MarkableF160}, ::Type{MarkableF163}) = MarkableF163
promote_type(::Type{MarkableF161}, ::Type{MarkableF162}) = MarkableF162
promote_type(::Type{MarkableF161}, ::Type{MarkableF163}) = MarkableF163
promote_type(::Type{MarkableF162}, ::Type{MarkableF163}) = MarkableF163

ext: jl


text:
function Base.show(io::IO, results::Results{Fl}) where Fl
    if !isempty(results)
        println(io, "                             Results                           ")
        println(io, "===============================================================" * adjust_charachters("=", results.coef_table))
        println(io, "Model:                        ", results.model_name)
        println(io, "Number of observations:       ", results.num_observations)
        println(io, "Number of unknown parameters: ", results.num_hyperparameters)
        println(io, "Log-likelihood:               ", pretty_number(results.llk))
        println(io, "AIC:                          ", pretty_number(results.aic))
        println(io, "AICc:                         ", pretty_number(results.aicc))
        println(io, "BIC:                          ", pretty_number(results.bic))
        print_coef_table(io, results.coef_table)
    else
        # Up to discussion of what should be the behaviour
        error("You must call fit! to get results.")
    end
    return nothing
end

function Base.show(io::IO, model::StateSpaceModel)
    print(io, model_name(model))
    return nothing
end

function print_coef_table(io::IO, coef_table::CoefficientTable{Fl}) where Fl
    # The repeats adjust the space for large names
    println(io, "---------------------------------------------------------------" * adjust_charachters("-", coef_table))
    println(io, "Parameter" * adjust_charachters(" ", coef_table) * "      Estimate      Std.Error      z stat      p-value")
    offset = 1
    for i in 1:length(coef_table)
        p_c, p_std, p_z_stat, p_p_val = print_coef(coef_table, offset)
        println(io, build_print(p_c, p_std, p_z_stat, p_p_val, offset, coef_table))
        offset += 1
    end
    return nothing
end

function print_coef(coef_table::CoefficientTable{Fl}, offset::Int) where Fl
    p_c = pretty_number(coef_table.coef[offset])
    p_std = pretty_number(coef_table.std_err[offset])
    p_z_stat = pretty_number(coef_table.z[offset])
    p_p_val = pretty_number(coef_table.p_value[offset])
    return p_c, p_std, p_z_stat, p_p_val
end

function build_print(p_c, p_std, p_z_stat, p_p_val, offset, coef_table)
    max_len_name = maximum(length.(coef_table.names))
    p = ""
    p *= coef_table.names[offset]
    p *= " "^max(0, 23 + max(0, max_len_name - 9) - length(p_c) - length(coef_table.names[offset]))
    p *= p_c
    p *= " "^max(0, 15 - length(p_std))
    p *= p_std
    p *= " "^max(0, 12 - length(p_z_stat))
    p *= p_z_stat
    p *= " "^max(0, 13 - length(p_p_val))
    p *= p_p_val
    return p
end

function pretty_number(num::Number)
    if isnan(num)
        return " - "
    elseif abs(num) > 1e4
        return @sprintf("%.3E", num)
    else
        return @sprintf("%.4f", num)
    end
end

function adjust_charachters(str::String, coef_table::CoefficientTable)
    max_len_name = maximum(length.(coef_table.names))
    return repeat(str, max(0, max_len_name - 9))
end


text:
# approximating tensor grids
# computing the approximating coefficients on
# a tensor product of basis functions
# NOTE: the fastest varying index in v is the one with highest
# index in ibm
function getTensorCoef(ibm::Dict{Integer, Array{T, 2}}, v::Vector{T}) where T <: Real

	# ibm are usually inverse basis matrices

	nall = length(v)	# length value vector
	nbm  = length(ibm)	# number of basis functions / dimensions of v

	# dim(matrices) must be same as length v
	prod = 1
	for (k,va) in ibm
		prod *= size(va,1)
		if !(size(va,1) == size(va,2))
			throw(ArgumentError("all matrices must be square"))
		end
	end
	if !(prod == nall)
		throw(ArgumentError("sizes of matrices not consistent with v"))
	end

	vtmp = 0.0

	# compute product for first matrix
	ks = sort(collect(keys(ibm)))
	v0    = copy(v)
	v1    = zeros(nall)
	stemp = ibm[ks[1]]
	n     = size(stemp,1)
	m     = round(Int,nall / n)
	for i in 1:m, j=1:n
		vtmp = 0.0
		for ji=1:n
			@inbounds vtmp += stemp[j,ji] * v0[n*(i-1) + ji]
		end
		v1[m*(j-1) + i] = vtmp
		# v[m*(j-1) + i] = vtmp
	end

	# compute for all other matrics
	if nbm > 1
		for imat in 2:nbm
			v0    = copy(v1)
			# v0    = copy(v)
			stemp = ibm[ks[imat]]
			n     = size(stemp,1)
			m     = round(Int,nall / n)
			# fill!(v1,0.0)
			# for i in 1:m, j=1:n
			for i in 1:m, j=1:n
				vtmp = 0.0
				for ji=1:n
					@inbounds vtmp += stemp[j,ji] * v0[n*(i-1) + ji]
				end
				# v[m*(j-1) + i] = vtmp
				v1[m*(j-1) + i] = vtmp
			end
		end
	end
	return v1
end

function evalTensor2(mat1::Vector{T}, mat2::Vector{T}, c::Vector) where T

	# TODO
	# sparse: if you had a type BSpline with field "nonzero" that
	# gives the index of the nonzero basis functions, this could
	# be sped up dramatically

	n1 = size(mat1,1)
	n2 = size(mat2,1)
	m1 = size(mat1,2)
	m2 = size(mat2,2)
	row_offset = 0
	col_offset = 0
	factor = 0.0

	r = zeros(n1 * n2)

	# loop over rows of 1
	for row_idx1 in 1:n1
		row_offset = (row_idx1-1) * n2

		# loop over rows of 2
		for row_idx2 in 1:n2

			# loop over cols of 1
			for col_idx1 in 1:m1
				col_offset = (col_idx1-1) * m2
				factor = mat1[row_idx1,col_idx1]

				# loop over cols of 2
				for col_idx2 in 1:m2

					# println("row_offset = $row_offset")
					# println("row_idx1 = $row_idx1")
					# println("row_idx2 = $row_idx2")
					# println("row_offset + row_idx1 = $(row_offset + row_idx1)")
					# println("factor = $factor")

					r[row_offset + row_idx2] += factor * mat2[row_idx2,col_idx2] * c[col_offset + col_idx2]
				end
			end
		end
	end
	return r
end

function evalTensor2(mat1::SparseMatrixCSC{T, Int64}, mat2::SparseMatrixCSC{T, Int64}, c::Vector{T}) where T

	# assume that the basis matrix is stored colwise
	# i.e. the first column are the basis functions
	# evaluated at the first point.

	# FIXME
	# julia 0.4 will introduce sparse row compression
	# so this can be changed in BSplines.jl

	if size(mat1)[2] > 1
		error("only doing column vector so far, sorry")
	end

	m1 = length(mat1)
	m2 = length(mat2)
	offset = 0
	factor = 0.0

	r = 0.0

	# full idiom is
	# for c in 1:size(mat1)[2]
	#    non_zero_indices[c] = mat1[mat1.colptr[c] : mat1.colptr[c+1]-1]
	#

	# for now (one column only) this is safe:
	# non_zero_indices = mat.rowval

	# loop over non-zero entries of mat1
	for idx1 in mat1.rowval
		offset = (idx1-1) * m2
		factor = mat1[idx1]

		# loop over non-zeros of mat2
		for idx2 in mat2.rowval

			r += factor * mat2[idx2] * c[offset + idx2]
		end
	end
	return r
end

function evalTensor3(mat1::SparseMatrixCSC{T, Int64}, mat2::SparseMatrixCSC{T, Int64}, mat3::SparseMatrixCSC{T, Int64}, c::Vector{T}) where T

	n1 = size(mat1,1)
	n2 = size(mat2,1)
	n3 = size(mat3,1)
	m1 = size(mat1,2)
	m2 = size(mat2,2)
	m3 = size(mat3,2)
	row_offset1= 0
	col_offset1= 0
	factor1= 0.0
	row_offset2= 0
	col_offset2= 0
	factor2= 0.0

	r = zeros(n1 * n2 * n3)

	# loop over rows of 1
	for row_idx1 in 1:n1
		row_offset1 = (row_idx1-1) * n2

		# loop over rows of 2
		for row_idx2 in 1:n2
			row_offset2 = (row_offset1 + (row_idx2-1)) * n3

			# loop rows 3
			for row_idx3 in 1:n3
				# println(row_offset2 + row_idx3)

				# loop over cols of 1
				for col_idx1 in 1:m1
					col_offset1 = (col_idx1-1) * m2
					factor1 = mat1[row_idx1,col_idx1]

					# loop over cols of 2
					for col_idx2 in 1:m2
						col_offset2 = (col_offset1 + (col_idx2-1)) * m3
						factor2 = factor1 * mat2[row_idx2,col_idx2]

						# cols mat3
						for col_idx3 in 1:m3

							r[row_offset2 + row_idx3] += factor2 * mat3[row_idx3,col_idx3] * c[col_offset2 + col_idx3]
						end
					end
				end
			end
		end
	end
	return r
end

function evalTensor3(mat1::Array{Float64,1},mat2::Array{Float64,1},mat3::Array{Float64,1},c::Vector)

	m1 = length(mat1)
	m2 = length(mat2)
	m3 = length(mat3)
	col_offset1= 0
	factor1= 0.0
	col_offset2= 0
	factor2= 0.0

	r = 0.0

	# loop over cols of 1
	for col_idx1 in 1:m1
		col_offset1 = (col_idx1-1) * m2
		factor1 = mat1[col_idx1]

		# loop over cols of 2
		for col_idx2 in 1:m2
			col_offset2 = (col_offset1 + (col_idx2-1)) * m3
			factor2 = factor1 * mat2[col_idx2]

			# cols mat3
			for col_idx3 in 1:m3
				r += factor2 * mat3[col_idx3] * c[col_offset2 + col_idx3]
			end
		end
	end
	return r
end

# loops over matrix rows
function evalTensor4(mat1::Array{Float64,2},mat2::Array{Float64,2},mat3::Array{Float64,2},mat4::Array{Float64,2},c::Vector)

	n1,m1 = size(mat1)
	n2,m2 = size(mat2)
	n3,m3 = size(mat3)
	n4,m4 = size(mat4)

	row_offset1 = 0
	col_offset1 = 0
	factor1     = 0.0
	row_offset2 = 0
	col_offset2 = 0
	factor2     = 0.0
	row_offset3 = 0
	col_offset3 = 0
	factor3     = 0.0

	r = zeros(n1 * n2 * n3 * n4)

	# loop over rows of 1
	for row_idx1 in 1:n1
		row_offset1 = (row_idx1-1) * n2

		# loop over rows of 2
		for row_idx2 in 1:n2
			row_offset2 = (row_offset1 + (row_idx2-1)) * n3

			# loop rows 3
			for row_idx3 in 1:n3
				row_offset3 = (row_offset2 + (row_idx3-1)) * n4
				# println(row_offset2 + row_idx3)

				for row_idx4 in 1:n4

					# loop over cols of 1
					for col_idx1 in 1:m1
						col_offset1 = (col_idx1-1) * m2
						factor1 = mat1[row_idx1,col_idx1]

						# loop over cols of 2
						for col_idx2 in 1:m2
							col_offset2 = (col_offset1 + (col_idx2-1)) * m3
							factor2 = factor1 * mat2[row_idx2,col_idx2]

							# cols mat3
							for col_idx3 in 1:m3
								col_offset3 = (col_offset2 + (col_idx3-1)) * m4
								factor3 = factor2 * mat3[row_idx3,col_idx3]

								for col_idx4 in 1:m4

									r[row_offset3 + row_idx4] += factor3 * mat4[row_idx4,col_idx4] * c[col_offset3 + col_idx4]
								end
							end
						end
					end
				end
			end
		end
	end
	return r
end

function evalTensor4(mat1::SparseMatrixCSC{T, Int64}, mat2::SparseMatrixCSC{T, Int64}, mat3::SparseMatrixCSC{T, Int64}, mat4::SparseMatrixCSC{T, Int64}, c::Vector{T}) where T

	if size(mat1)[2] > 1
		error("only doing column vector so far, sorry")
	end

	m1 = length(mat1)
	m2 = length(mat2)
	m3 = length(mat3)
	m4 = length(mat4)
	offset1= 0
	factor1= 0.0
	offset2= 0
	factor2= 0.0
	offset3= 0
	factor3= 0.0

	r = 0.0

	# loop over non-zeros in mat1
	for idx1 in mat1.rowval
		offset1 = (idx1-1) * m2
		factor1 = mat1[idx1]

		# loop over non-zeros in mat2
		for idx2 in mat2.rowval
			offset2 = (offset1 + (idx2-1)) * m3
			factor2 = factor1 * mat2[idx2]

			# cols mat3
			for idx3 in mat3.rowval
				offset3 = (offset2 + (idx3-1)) * m4
				factor3 = factor2 * mat3[idx3]

				for idx4 in mat4.rowval
					r += factor3 * mat4[idx4] * c[offset3 + idx4]
				end
			end
		end
	end
	return r
end

function evalTensor(bs::Dict{Int64, SparseMatrixCSC{T, Int64}}, c::Array{T}) where T

	n = length(bs)
	if n == 2
		evalTensor2(bs[2],bs[1],c)
	elseif n==3
		evalTensor3(bs[3],bs[2],bs[1],c)
	elseif n==4
		evalTensor4(bs[4],bs[3],bs[2],bs[1],c)
	else
		throw(ErrorException("only dims 2,3,4 implemented so far"))
	end
end


text:
module Tracks
using OffsetArrays, Images, Unitful
export Track, SIMTrack

using ..Acquisitions

abstract type AbstractTrack end

struct Track <: AbstractTrack # {{{
    x::AbstractVector{<:Real}
    y::AbstractVector{<:Real}
    view::AxisArray{A,1} where {A<:OffsetArray{<:AbstractGray}}
end # }}}

# TODO: Should be generic over the SIMacquisition type <20-05-21, kunzaatko> #
struct SIMTrack <: AbstractTrack # {{{
    x::AbstractVector{<:Real}
    y::AbstractVector{<:Real}
    view::AxisArray{A,1} where {A<:OffsetArray{<:AbstractGray}}
    view_reconstructed::AxisArray{A,1} where {A<:OffsetArray{<:AbstractGray}}
end # }}}

function Track( # {{{
    x, y, acquisition::Acquisition, times::AbstractVector{<:Unitful.Time}
)
    # TODO: Add check for times, x and y having the same lengths <20-05-21, kunzaatko> #
    x_offsets, y_offsets = round.(Integer, x), round.(Integer, y)
    view = AxisArray(
        [
            OffsetArray(
                acquisition.images[Axis{:x}(:), Axis{:y}(:), Axis{:time}(t)],
                OffsetArrays.Origin(-x_offset, -y_offset),
            ) for (t, x_offset, y_offset) in zip(times, x_offsets, y_offsets)
        ],
        Axis{:time}(times),
    )
    return Track(x, y, view)
end # }}}

function Track(x, y, acquisition::Acquisition, frames::AbstractVector{<:Integer})
    return Track(x, y, acquisition, times(acquisition)[frames])
end

function SIMTrack( # {{{
    x, y, acquisition::SIMAcquisition, times::AbstractVector{<:Unitful.Time}
)
    x_offsets, y_offsets = round.(Integer, x), round.(Integer, y)
    view = AxisArray(
        [
            OffsetArray(
                acquisition.unreconstructedimages[Axis{:x}(:), Axis{:y}(:), Axis{:time}(t)],
                OffsetArrays.Origin(-x_offset, -y_offset),
            ) for (t, x_offset, y_offset) in zip(times, x_offsets, y_offsets)
        ],
        Axis{:time}(times),
    )
    # FIX: This should be meaningfull times as they will be defined in the constructor <20-05-21, kunzaatko> #
    view_reconstructed = AxisArray(
        [OffsetArray(Gray.(zeros(3, 3)), OffsetArrays.Origin(-1, -1)) for _ in [1, 2]],
        Axis{:time}([1u"s", 2u"s"]),
    ) # DUMMY instance
    return SIMTrack(x, y, view, view_reconstructed)
end # }}}

function SIMTrack(x, y, acquisition::SIMAcquisition, frames::AbstractVector{<:Integer})
    return SIMTrack(x, y, acquisition, times(acquisition)[frames])
end

end


text:
"""
    Pearson Correlation

Find the pearson correlation between two variables.

# Example:

julia> PearsonCorrelation([12,11,16,17,19,21],[11,51,66,72,12,15])
-0.2092706263573845

Contribution by: [Aru Bhardwaj](https://github.com/arubhardwaj)

"""

function PearsonCorrelation(x, y)
    mean_x = sum(x)/length(x)
    mean_y = sum(y)/length(y)
    XY = (x.-mean_x).*(y.-mean_y)
    XXs = sum((x.-mean_x).*(x.-mean_x))
    YYs = sum((y.-mean_y).*(y.-mean_y))
    return(sum(XY)/(sqrt(XXs.*YYs)))
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
# ---
# title: 647. Palindromic Substrings
# id: problem647
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/palindromic-substrings/description/>
# hidden: true
# ---
#
# Given a string, your task is to count how many palindromic substrings in this
# string.
#
# The substrings with different start indexes or end indexes are counted as
# different substrings even they consist of same characters.
#
# **Example 1:**
#
#
#
#     Input: "abc"
#     Output: 3
#     Explanation: Three palindromic strings: "a", "b", "c".
#
#
#
#
# **Example 2:**
#
#
#
#     Input: "aaa"
#     Output: 6
#     Explanation: Six palindromic strings: "a", "a", "a", "aa", "aa", "aaa".
#
#
#
#
# **Note:**
#
#   1. The input string length won't exceed 1000.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export get_agent_view, AbstractGridWorld
export get_actions
abstract type AbstractGridWorld <: AbstractEnv end

function get_agent_view end
function get_agent end

Base.convert(::Type{GridWorldBase}, env::AbstractGridWorld) = env.world
get_object(env::AbstractGridWorld) = get_object(convert(GridWorldBase, env))
get_object(env::AbstractGridWorld, x::Type{<:AbstractObject}) = filter(o -> o isa x, get_object(env))
get_object(env::AbstractGridWorld, x::Type{Agent}) = env.agent
get_pos(env::AbstractGridWorld, ::Type{Agent}) = env.agent_pos

get_agent(env::AbstractGridWorld) = get_object(env, Agent)
get_agent_pos(env::AbstractGridWorld) = get_pos(env, Agent)
get_agent_dir(env::AbstractGridWorld) = env |> get_agent |> get_dir

function get_agent_view(env::AbstractGridWorld, agent_view_size=(7,7))
    w = convert(GridWorldBase, env)
    v = BitArray{3}(undef, size(w, 1), agent_view_size...)
    fill!(v, false)
    get_agent_view!(v, env)
end

get_agent_view_inds(env::AbstractGridWorld, s=(7,7)) = get_agent_view_inds(get_agent_pos(env).I, s, get_agent_dir(env))

get_agent_view!(v::BitArray{3}, env::AbstractGridWorld) = get_agent_view!(v, convert(GridWorldBase, env), get_agent_pos(env), get_agent_dir(env))

#####
# RLBase defaults
#####

RLBase.DefaultStateStyle(env::AbstractGridWorld) = RLBase.PartialObservation{Array}()

RLBase.get_state(env::AbstractGridWorld, ::RLBase.PartialObservation{Array}, args...) = get_agent_view(env)

RLBase.get_actions(env::AbstractGridWorld) = (MOVE_FORWARD, TURN_LEFT, TURN_RIGHT)

RLBase.get_reward(env::AbstractGridWorld) = env.reward

function (env::AbstractGridWorld)(action::Union{TurnRight, TurnLeft})
    env.reward = 0.0
    agent = get_agent(env)
    set_dir!(agent, action(get_dir(agent)))
    env
end


text:
using Test
using MPI
using CLIMA.Topologies

function main()
  MPI.Init()
  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  @assert csize == 2

  DFloat = Float64
  Nx = 3
  Ny = 2
  x = range(DFloat(0); length=Nx+1, stop=1)
  y = range(DFloat(0); length=Ny+1, stop=1)

  topology = BrickTopology(comm, (x,y); boundary=[1 3; 2 4],
                           periodicity=(true,true))

  elems       = topology.elems
  realelems   = topology.realelems
  ghostelems  = topology.ghostelems
  sendelems   = topology.sendelems
  elemtocoord = topology.elemtocoord
  elemtoelem  = topology.elemtoelem
  elemtoface  = topology.elemtoface
  elemtoordr  = topology.elemtoordr
  elemtobndy  = topology.elemtobndy
  nabrtorank  = topology.nabrtorank
  nabrtorecv  = topology.nabrtorecv
  nabrtosend  = topology.nabrtosend

  globalelemtoface = [2  2  2  2  2  2
                      1  1  1  1  1  1
                      4  4  4  4  4  4
                      3  3  3  3  3  3]

  globalelemtoordr = ones(Int, size(globalelemtoface))
  globalelemtobndy = zeros(Int, size(globalelemtoface))

  if crank == 0
    nrealelem = 3
    globalelems = [1,2,3,4,5,6]
    elemtoelem_expect = [6 4 2 4 5 6
                         5 3 4 4 5 6
                         2 1 5 4 5 6
                         2 1 5 4 5 6]
    nabrtorank_expect = [1]
    nabrtorecv_expect = UnitRange{Int}[1:3]
    nabrtosend_expect = UnitRange{Int}[1:3]
  elseif crank == 1
    nrealelem = 3
    globalelems = [4,5,6,1,2,3]
    elemtoelem_expect = [6 4 2 4 5 6
                         5 3 4 4 5 6
                         3 6 1 4 5 6
                         3 6 1 4 5 6]
    nabrtorank_expect = [0]
    nabrtorecv_expect = UnitRange{Int}[1:3]
    nabrtosend_expect = UnitRange{Int}[1:3]
  end

  @test elems == 1:length(globalelems)
  @test realelems == 1:nrealelem
  @test ghostelems == nrealelem+1:length(globalelems)
  @test elemtoface[:,realelems] == globalelemtoface[:,globalelems[realelems]]
  @test elemtoelem == elemtoelem_expect
  @test elemtobndy == globalelemtobndy[:,globalelems]
  @test elemtoordr == ones(eltype(elemtoordr), size(elemtoordr))
  @test nabrtorank == nabrtorank_expect
  @test nabrtorecv == nabrtorecv_expect
  @test nabrtosend == nabrtosend_expect

  MPI.Finalize()
end

main()

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
export Renderable
export draw, image_surface, ColorLine, ColorRect, ColorCirc, ColorTri, Sprite, Text

draw(e::Entity) = nothing

abstract type Renderable <: Entity end

# NOTE all points are considered as offsets from the entity's position
# NOTE rotation is ignored for these simple 2d shapes, you
# either can't rotate them or can transform the points

mutable struct ColorLine <: Renderable
  function ColorLine(p1, p2; color=colorant"white", kw...)
    instantiate!(new(); p1=p1, p2=p2, color=color, kw...)
  end
end

function draw(l::ColorLine)
  SDL_SetRenderDrawColor(
      sdl.rnd,
      sdl_colors(l.color)...,
  )
  SDL_RenderDrawLine(sdl.rnd, Cint.((
    l.p1[1]+l.abs_pos.x, l.p1[2]+l.abs_pos.y,
    l.p2[1]+l.abs_pos.x, l.p2[2]+l.abs_pos.y))...)
end

mutable struct ColorRect <: Renderable
  function ColorRect(p, w, h; fill=true, color=colorant"white", kw...)
    instantiate!(new(); p=p, w=w, h=h, color=color, fill=fill, kw...)
  end
end

Base.convert(T::Type{SDL_Rect}, r::ColorRect) = SDL_Rect(
  Cint.((r.p[1]+r.abs_pos.x-Int(round(r.w/2)),
        r.p[2]+r.abs_pos.y-Int(round(r.h/2)),
        r.w, r.h))...)

function draw(r::ColorRect)
  SDL_SetRenderDrawColor(
      sdl.rnd,
      sdl_colors(r.color)...,
  )
  sr = convert(SDL_Rect, r)
  if !r.fill
      SDL_RenderDrawRect(sdl.rnd, Ref(sr))
  else
      SDL_RenderFillRect(sdl.rnd, Ref(sr))
  end
end

mutable struct ColorCirc <: Renderable
  function ColorCirc(p, r; fill=true, color=colorant"white", kw...)
    instantiate!(new(); p=p, r=r, fill=fill, color=color, kw...)
  end
end

# improved circle drawing algorithm. slower but fills completely. needs optimization
function draw(circle::ColorCirc)
  # define the center and needed sides of circle
  centerX = Cint(circle.p[1]+circle.abs_pos.x)
  centerY = Cint(circle.p[2]+circle.abs_pos.y)
  int_rad = Cint(circle.r)
  left = centerX - int_rad
  top = centerY - int_rad

  SDL_SetRenderDrawColor(
      sdl.rnd,
      sdl_colors(circle.color)...,
  )

  # we consider a grid with sides equal to the circle's diameter
  for x in left:centerX
    for y in top:centerY

      # for each pixel in the top left quadrant of the grid we measure the distance from the center.
      dist = sqrt( (centerX - x)^2 + (centerY - y)^2 )

      # if it is close to the circle's radius it and all associated points in the other quadrants are colored in.
      if (dist <= circle.r + 0.5 && dist >= circle.r - 0.5)
        rel_x = centerX - x
        rel_y = centerY - y

        quad1 = (x              , y              )
        quad2 = (centerX + rel_x, y              )
        quad3 = (x              , centerY + rel_y)
        quad4 = (quad2[1]       , quad3[2]       )

        SDL_RenderDrawPoint(sdl.rnd, quad1[1], quad1[2])
        SDL_RenderDrawPoint(sdl.rnd, quad2[1], quad2[2])
        SDL_RenderDrawPoint(sdl.rnd, quad3[1], quad3[2])
        SDL_RenderDrawPoint(sdl.rnd, quad4[1], quad4[2])

        # if we are told to fill in the circle we draw lines between all of the quadrants to completely fill the circle
        if (circle.fill)
          SDL_RenderDrawLine(sdl.rnd, quad1[1], quad1[2], quad2[1], quad2[2])
          SDL_RenderDrawLine(sdl.rnd, quad2[1], quad2[2], quad4[1], quad4[2])
          SDL_RenderDrawLine(sdl.rnd, quad4[1], quad4[2], quad3[1], quad3[2])
          SDL_RenderDrawLine(sdl.rnd, quad3[1], quad3[2], quad1[1], quad1[2])
        end
      end

    end
  end

end

mutable struct ColorTri <: Renderable
  function ColorTri(p1, p2, p3; fill=true, color=colorant"white", kw...)
    instantiate!(new(); p1=p1, p2=p2, p3=p3, color=color, fill=fill, kw...)
  end
end

function draw(tr::ColorTri)
  p1 = Cint.([tr.p1[1]+tr.abs_pos.x, tr.p1[2]+tr.abs_pos.y])
  p2 = Cint.([tr.p2[1]+tr.abs_pos.x, tr.p2[2]+tr.abs_pos.y])
  p3 = Cint.([tr.p3[1]+tr.abs_pos.x, tr.p3[2]+tr.abs_pos.y])
  SDL_SetRenderDrawColor(sdl.rnd, sdl_colors(tr.color)...)
  SDL_RenderDrawLines(sdl.rnd, [p1; p2; p3; p1], Cint(4))

  ymax = max(p1[2], p2[2], p3[2])
  ymin = min(p1[2], p2[2], p3[2])
  if tr.fill && ymin != ymax
    # Set q1, q2 and q3 in descending order of y-value
    q1 = (p1[2] != ymax != p2[2]) * p3 +
        (p2[2] != ymax != p3[2]) * p1 +
        (p3[2] != ymax != p1[2]) * p2 +
        (p1[2] == p2[2] == ymax) * p2 +
        (p2[2] == p3[2] == ymax) * p3 +
        (p3[2] == p1[2] == ymax) * p1
    q3 = (p1[2] != ymin != p2[2]) * p3 +
        (p2[2] != ymin != p3[2]) * p1 +
        (p3[2] != ymin != p1[2]) * p2 +
        (p1[2] == p2[2] == ymin) * p2 +
        (p2[2] == p3[2] == ymin) * p3 +
        (p3[2] == p1[2] == ymin) * p1
    q2 = ((q1 == p1 && q3 == p3) || (q1 == p3 && q3 == p1)) * p2 +
        ((q1 == p1 && q3 == p2) || (q1 == p2 && q3 == p1)) * p3 +
        ((q1 == p2 && q3 == p3) || (q1 == p3 && q3 == p2)) * p1

    n = q1[2] - q2[2]
    x0 = q1[1] + (q2[2] - q1[2]) / (q3[2] - q1[2]) * (q3[1] - q1[1])
    for j = Cint(0):n-Cint(1)
      r1 = [round(Cint, q2[1] + j / n * (q1[1] - q2[1])); q2[2] + j]
      r2 = [round(Cint, x0 + j / n * (q1[1] - x0)); q2[2] + j]
      SDL_RenderDrawLines(sdl.rnd, [r1; r2], Cint(2))
    end
    n = q2[2] - q3[2]
    for j = Cint(1):n-Cint(1)
      r1 = [round(Cint, q2[1] + j / n * (q3[1] - q2[1])); q2[2] - j]
      r2 = [round(Cint, x0 + j / n * (q3[1] - x0)); q2[2] - j]
      SDL_RenderDrawLines(sdl.rnd, [r1; r2], Cint(2))
    end
  end
end

function Base.size(s::Ptr{SDL_Surface})
  ss = unsafe_load(s)
  return (ss.w, ss.h)
end

function image_surface(img)
  sf = IMG_Load(img)
  if sf == C_NULL
    throw("Error loading $image_file")
  end
  return sf
end

mutable struct Sprite <: Renderable
  # allow textures larger than a single sprite
  # by supporting cell_size, cell_ind, and region
  # fields. cell_size turns img into a matrix of
  # adjacent regions of the same size starting in
  # the top left of the texture. cell_ind is a 0-indexed
  # (row,col) index into a particular cell. region
  # overrides both and denotes the top-left corner
  # of a rectangle, its width, and its height, and
  # uses only that region of the texture.
  function Sprite(img; cell_size=[0, 0], region=[0, 0, 0, 0], cell_ind=[0, 0],
    alpha=UInt8(255), scale=XYZ(1,1,1), kw...)
    instantiate!(new(); img=img, cell_size=cell_size,
      region=region, cell_ind=cell_ind,
      alpha=alpha, scale=scale, kw...)
  end
end

function draw(s::Sprite)
  srf = image_surface(s.img)
  w, h = size(srf)
  txt = SDL_CreateTextureFromSurface(sdl.rnd, srf)

  if (s.alpha < 255)
    SDL_SetTextureBlendMode(txt, SDL_BLENDMODE_BLEND)
    SDL_SetTextureAlphaMod(txt, s.alpha)
  end

  # calculations for source
  if s.region != [0, 0, 0, 0]
    srctl = [s.region[1], s.region[2]]
    srcw, srch = [s.region[3], s.region[4]]
  elseif s.cell_size != [0, 0]
    srctl = [s.cell_ind[2] * s.cell_size[1], s.cell_ind[1] * s.cell_size[2]]
    srcw, srch = s.cell_size
  else
    srctl = [0, 0]
    srcw, srch = [w, h]
  end

  # calculations for destination
  if s.region != [0, 0, 0, 0]
    dstw, dsth = s.region[3], s.region[4]
  elseif s.cell_size != [0, 0]
    dstw, dsth = s.cell_size
  else
    dstw, dsth = [w, h]
  end

  dstw = floor(dstw * s.scale.x)
  dsth = floor(dsth * s.scale.y)

  SDL_RenderCopyEx(
    sdl.rnd,
    txt,
    Ref(SDL_Rect(srctl[1], srctl[2], srcw, srch)),
    Ref(SDL_Rect(s.abs_pos.x-Int(round(dstw/2)), s.abs_pos.y-Int(round(dsth/2)), dstw, dsth)),
    s.abs_rot.z, # NOTE angles are in degrees, positive is clockwise
    C_NULL,
    SDL_FLIP_NONE
  )

  SDL_DestroyTexture(txt)
  SDL_FreeSurface(srf)
end

mutable struct Text <: Renderable
  function Text(text, font_name; font_size=12, color=colorant"black", scale=XYZ(1,1,1), alpha=UInt8(255), kw...)
    instantiate!(new(),
      text=text, font_name=font_name,
      font_size=font_size, scale=scale,
      color=color, alpha=alpha, kw...)
  end
end

function draw(t::Text)
  font = TTF_OpenFont(t.font_name, t.font_size)
  srf = TTF_RenderText_Blended(font, t.text, SDL_Color(sdl_colors(t.color)...))
  w, h = size(srf)

  w = floor(w * t.scale.x)
  h = floor(h * t.scale.y)
  txt = SDL_CreateTextureFromSurface(sdl.rnd, srf)

  if t.alpha < 255
    SDL_SetTextureBlendMode(txt, SDL_BLENDMODE_BLEND)
    SDL_SetTextureAlphaMod(txt, t.alpha)
  end

  SDL_RenderCopyEx(
    sdl.rnd,
    txt,
    C_NULL,
    Ref(SDL_Rect(t.abs_pos.x-Int(round(w/2)), t.abs_pos.y-Int(round(h/2)), w, h)),
    t.abs_rot.z,
    C_NULL,
    SDL_FLIP_NONE
  )

  SDL_DestroyTexture(txt)
  SDL_FreeSurface(srf)
  TTF_CloseFont(font)
end


text:
println("\n\n\nPerformance of FinanceLib.Derivatives.Forwards functions")

import FinanceLib.Derivatives.Forwards as Fd
using Dates: Date
using FinanceLib: yearFrac

# fw = Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50, benf = -1)
# @code_warntype Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50, benf = -1)
# @code_warntype fw.FT
# @code_warntype Fd.getS0(fw)
# @code_warntype Fd.valuation(fw, 52, 1/12)

# fw = Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50)
# @code_warntype Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50, benf = -1)
# @code_warntype fw.FT
# @code_warntype Fd.getS0(fw)

# T0 = Date(2015,10,5)
# T1 = Date(2016,1,8)
# TX = yearFrac(T0, T1)

# fw = Fd.Forward_S0(r = 0.03, T = TX, S0 = 50, benf = -1)
# # println(fw)

# dw = Fd.XForward_S0(r = 0.03, T0 = T0, T = T1, S0 = 50, benf = -1.0)
# # println(dw)

# @code_warntype Fd.getS0(dw)

# TM = Date(2015, 12, 18)
# TZ = yearFrac(T0, TM)
# # println(TZ)

# @code_warntype Fd.valuation(fw, 55, TZ)
# @code_warntype Fd.valuation(dw, 55, TM)

# println("Benchmarks\n")

# println("\n\n\nPerformance of FinanceLib.Derivatives.Options functions")

# import FinanceLib.Derivatives.Options as Op

# ct = Op.Call(2_000, 81.75)

# @code_warntype Op.valuation(ct, ST = 1_900)
# @code_warntype Op.profit(ct, ST = 1_900)

# @code_warntype Op.valuation(ct, ST = 2_100)
# @code_warntype Op.profit(ct, ST = 2_100)

# pt = Op.Put(2_000, 79.25)

# @code_warntype Op.valuation(pt, ST = 1_900)
# @code_warntype Op.profit(pt, ST = 1_900)

# @code_warntype Op.valuation(pt, ST = 2_100)
# @code_warntype Op.profit(pt, ST = 2_100)

# cv = Op.CoveredCall(98, Op.Call(105, 8))

# @code_warntype Op.valuation(cv, ST = 110)
# @code_warntype Op.profit(cv, ST = 110)

# @code_warntype Op.valuation(cv, ST = 88)
# @code_warntype Op.profit(cv, ST = 88)

# pt = Op.ProtectivePut(0.875, Op.Put(0.90, 0.075))

# @code_warntype Op.valuation(pt, ST = 0.96)
# @code_warntype Op.profit(pt, ST = 0.96)

# @code_warntype Op.valuation(pt, ST = 0.75)
# @code_warntype Op.profit(pt, ST = 0.75)

println("Benchmarks\n")


text:
using FASTX
using Random
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--sequences", "-s"
        arg_type = String
        required = true
        help = "Path to the FASTA file containing the protein sequences."
    "--interactions", "-i"
        arg_type = String
        required = true
        help = "Path to the tsv file containing the PPIs."
    "--output", "-o"
        arg_type = String
        help = "Path to the FASTA output file containing the reduced set of sequences."
        required = true
        arg_type = String
    "--threshold", "-t"
        default = 0.70
        arg_type = Float64
        help = "The identity threshold to be used for clustering"
        required = false
end

args = parse_args(ARGS, s)

function compute_word_size(threshold::Float64)
    if threshold >= 0.7
        return 5
    elseif threshold > 0.6
        return 4
    elseif threshold > 0.5
        return 3
    else
        return 2
    end
end

SIMILARITY_THRESHOLD = (args["threshold"], compute_word_size(args["threshold"])) # threshold, word size

# Cluster the sequences
println("Clustering the sequences...")
cmd = `cd-hit -i $(args["sequences"]) -o /tmp/temp_clusters.fasta -n $(SIMILARITY_THRESHOLD[2]) -c $(SIMILARITY_THRESHOLD[1]) -T 8`
run(cmd)

# Interaction counting
println("Counting the number of interactions proteins are involved in...")

interaction_counts = Dict()
for pair in readlines(args["interactions"])
    split_pair = split(pair, "\t")
    for protein in split_pair
        if !haskey(interaction_counts, protein)
            interaction_counts[protein] = 1
        else
            interaction_counts[protein] += 1
        end
    end
end

# Add a value of 0 to proteins which are involed in 0 interactions
for protein in FASTA.Reader(open(args["sequences"], "r"))
	id = FASTA.identifier(protein)
	if !haskey(interaction_counts, id)
		interaction_counts[id] = 0
	end
end

println("Identifying the cluster representatives based on interaction counts...")
representatives = []
cluster_file = read("/tmp/temp_clusters.fasta.clstr", String)
clusters = split(cluster_file, ">C")[2:end]

# For every cluster, identify the proteins in the cluster
for cluster in clusters
    lines = split(cluster, "\n")
    if length(lines) == 4
        representative = lines[2]
        push!(representatives, representative)
    else
        best_candidate = ""
        best_score = -1
        for candidate in lines
            m = match(r"[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}", candidate)
            if m == nothing
                continue
            else
                if interaction_counts[m.match] > best_score
                    best_candidate = m.match
                    best_score = interaction_counts[m.match]
                end
            end
        end
        push!(representatives, best_candidate)
    end
end

println("Saving the cluster representatives...")
representative_set = Set(representatives)
selected_proteins = []
for protein in FASTA.Reader(open(args["sequences"], "r"))
    id = FASTA.identifier(protein)
    if in(id, representative_set)
        push!(selected_proteins, protein)
    end
end

w = open(FASTA.Writer, args["output"])
for p in selected_proteins
    write(w, p)
end

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
"""
Adapted from sandreza/Learning/sandbox/oceananigans_converter.jl
https://github.com/sandreza/Learning/blob/master/sandbox/oceananigans_converter.jl
"""

using JLD2

struct LESbraryData{𝒮, 𝒯, 𝒰, 𝒱}
    # Initial conditions
    T⁰ :: 𝒮
    U⁰ :: 𝒮
    V⁰ :: 𝒮

    # Fields at each moment in time
    T :: 𝒯
    U :: 𝒯
    V :: 𝒯

    # Some second order statistics at each moment in time
    wT :: 𝒯
    wu :: 𝒯
    wv :: 𝒯
    # uu :: 𝒯
    # vv :: 𝒯
    # ww :: 𝒯

    # Subfilter fluxes
    # νₑ_∂z_u :: 𝒯
    # νₑ_∂z_v :: 𝒯
    # νₑ_∂z_w :: 𝒯
    # κₑ_∂z_T :: 𝒯

    # Simulation constants
    α  :: 𝒰
    β  :: 𝒰
    f⁰ :: 𝒰
     g :: 𝒰
     L :: 𝒰

    # Time and grid
    t  :: 𝒮
    zC :: 𝒮
    zF :: 𝒮

    # Boundary conditions
    θ_top    :: 𝒰
    u_top    :: 𝒰
    θ_bottom :: 𝒰
    u_bottom :: 𝒰

    # Info about the simulation
    info :: 𝒱
end

function ReadJLD2_LESbraryData(filename)
    les_data = jldopen(filename, "r")
    les_keys = keys(les_data)
    timeseries_keys = keys(les_data["timeseries"]["t"])
    t₀ = timeseries_keys[1]

    # hold the entries for easy constructor creation
    container = []

    # size of arrays
    Nz = length(collect(les_data["grid"]["zC"])) - 6
    Nt = length(timeseries_keys)

    # Initial Conditions
    T⁰ = zeros(typeof(les_data["timeseries/T/$t₀"][1]), Nz)
    U⁰ = zeros(typeof(les_data["timeseries/u/$t₀"][1]), Nz)
    V⁰ = zeros(typeof(les_data["timeseries/v/$t₀"][1]), Nz)

    # Timeseries
    T = zeros(typeof(les_data["timeseries/T/$t₀"][1]), Nz, Nt)
    U = zeros(typeof(les_data["timeseries/u/$t₀"][1]), Nz, Nt)
    V = zeros(typeof(les_data["timeseries/v/$t₀"][1]), Nz, Nt)
    t = zeros(typeof(les_data["timeseries/t/$t₀"]), Nt)

    # Second Order Statistics
    wT  = zeros(typeof(les_data["timeseries/wT/$t₀"][1]), Nz+1, Nt)
    # uu  = zeros(typeof(les_data["timeseries/uu/$t₀"][1]), Nz,   Nt)
    # vv  = zeros(typeof(les_data["timeseries/vv/$t₀"][1]), Nz,   Nt)
    # ww  = zeros(typeof(les_data["timeseries/ww/$t₀"][1]), Nz+1, Nt)
    wu  = zeros(typeof(les_data["timeseries/wu/$t₀"][1]), Nz+1, Nt)
    wv  = zeros(typeof(les_data["timeseries/wv/$t₀"][1]), Nz+1, Nt)

    # Subfilter momentum fluxes
    # νₑ_∂z_u  = zeros(typeof(les_data["timeseries/νₑ_∂z_u/$t₀"][1]), Nz+1, Nt)
    # νₑ_∂z_v  = zeros(typeof(les_data["timeseries/νₑ_∂z_v/$t₀"][1]), Nz+1, Nt)
    # νₑ_∂z_w  = zeros(typeof(les_data["timeseries/νₑ_∂z_w/$t₀"][1]), Nz, Nt)
    # κₑ_∂z_T  = zeros(typeof(les_data["timeseries/κₑ_∂z_T/$t₀"][1]), Nz+1, Nt)

    # grab arrays
    for j in 1:Nt
        key = timeseries_keys[j]

        # Fields
        @. T[:, j] = les_data["timeseries/T/$key"][1, 1, :]
        @. U[:, j] = les_data["timeseries/u/$key"][1, 1, :]
        @. V[:, j] = les_data["timeseries/v/$key"][1, 1, :]

        # Second Order Statistics
        @. wT[:, j] = les_data["timeseries/wT/$key"][1, 1, :]
        @. wu[:, j] = les_data["timeseries/wu/$key"][1, 1, :]
        @. wv[:, j] = les_data["timeseries/wv/$key"][1, 1, :]
        # @. uu[:, j] = les_data["timeseries/uu/$key"][1, 1, :]
        # @. vv[:, j] = les_data["timeseries/vv/$key"][1, 1, :]
        # @. ww[:, j] = les_data["timeseries/ww/$key"][1, 1, :]

        # Subfilter fluxes
        # @. νₑ_∂z_u[:, j] = les_data["timeseries/νₑ_∂z_u/$key"][1, 1, :]
        # @. νₑ_∂z_v[:, j] = les_data["timeseries/νₑ_∂z_v/$key"][1, 1, :]
        # @. νₑ_∂z_w[:, j] = les_data["timeseries/νₑ_∂z_w/$key"][1, 1, :]
        # @. κₑ_∂z_T[:, j] = les_data["timeseries/κₑ_∂z_T/$key"][1, 1, :]

        t[j] = les_data["timeseries"]["t"][key]
    end

    # Set initial Conditions
    @. T⁰ = T[:,1]
    @. U⁰ = U[:,1]
    @. V⁰ = V[:,1]

    # Push initial conditions current stuff into container
    push!(container, T⁰, V⁰, U⁰)

    # Push fields into container
    push!(container, T, U, V)

    # Push second order statistics into container
    # push!(container, wT, wu, wv, uu, vv, ww)
    push!(container, wT, wu, wv)

    # Push subfilter fluxes into container
    # push!(container, νₑ_∂z_u, νₑ_∂z_v, νₑ_∂z_w, κₑ_∂z_T)

    # Now grab parameter
    if haskey(les_data["buoyancy"], "model")
        α = les_data["buoyancy/model/equation_of_state/α"]
        β = les_data["buoyancy/model/equation_of_state/β"]
        g = les_data["buoyancy/model/gravitational_acceleration"]
    else
        α = les_data["buoyancy/equation_of_state/α"]
        β = les_data["buoyancy/equation_of_state/β"]
        g = les_data["buoyancy/gravitational_acceleration"]
    end

    f⁰ = les_data["coriolis"]["f"]
    L = les_data["grid"]["Lz"]

    # Push parameters to container
    push!(container, α, β, f⁰, g, L)

    # grab domain data
    zC = collect(les_data["grid"]["zC"])[4:end-3] # padding of 3 on each side
    zF = collect(les_data["grid"]["zF"])[4:end-3] # padding of 3 on each side

    # push
    push!(container, t, zC, zF)

    # now grab boundary condition data
    # θ_top = les_data["parameters/boundary_condition_θ_top"]
    buoyancy_flux = les_data["parameters/buoyancy_flux"]
    θ_top = buoyancy_flux / α / g
    u_top = les_data["parameters/boundary_condition_u_top"]
    θ_bottom = les_data["parameters/boundary_condition_θ_bottom"]
    u_bottom = les_data["parameters/boundary_condition_u_bottom"]

    # push to container
    push!(container, θ_top, u_top, θ_bottom, u_bottom)

    # Now construct types
    𝒮 = typeof(U⁰)
    𝒯 = typeof(T)
    𝒰 = typeof(α)
    𝒱 = typeof("string")

    # now create data string
    info_string = "The top boundary conditions are flux boundary conditions \n"
    info_string *= "The  bottom boundary condition for temperature is a gradient boundary condition \n"
    info_string *= "The grid data is assumed to be evenly spaced and a power of two \n"

    # push to container
    push!(container, info_string)

    close(les_data)

    return LESbraryData{𝒮, 𝒯, 𝒰, 𝒱}(container...)
end

# TESTING
# file = "2DaySuite/three_layer_constant_fluxes_hr48_Qu0.0e+00_Qb1.2e-07_f1.0e-04_Nh256_Nz128_free_convection_statistics.jld2"
# les = jldopen(file, "r")
# keys(les["timeseries"])
# les["timeseries"]["κₑ_∂z_T"]
# les["timeseries"]["v"]["0"][1,1,:]
# les["grid"]["Δx"]

# using Plots
# file = "2DaySuite/three_layer_constant_fluxes_hr48_Qu0.0e+00_Qb1.2e-07_f1.0e-04_Nh256_Nz128_free_convection_statistics.jld2"
# myles = ReadJLD2_LESbraryData(file)
#
# vbl = myles.κₑ_∂z_T
# xlims=(-1.2e-5, maximum(myles.κₑ_∂z_T))
# vbl = myles.wT
# xlims=(minimum(myles.wT), maximum(myles.wT))
#
# anim = @animate for i=1:288
#     plot(vbl[:,i], myles.zF, legend=false, xlims=xlims, size=(400,600))
# end
# gif(anim, "wT.gif")
#
# Δt = diff(myles.t, dims=1)'
# Nz,Nt = size(myles.T)
# dudt = (myles.u[:,2:Nt] .- myles.u[:,1:Nt-1]) ./ Δt # Nz x (Nt-1) array of approximate dUdt values


text:
############################################################################################
# static entities
############################################################################################
abstract type Entity end
radius(entity::Entity) = entity.radius
id(entity::Entity) = entity.id
distance_between(e1::Entity, e2::Entity) = √((e1.x - e2.x)^2 + (e1.y - e2.y)^2)
angle_between(e1::Entity, e2::Entity) = rad2deg(atan2(e2.y - e1.y, e2.x - e1.x)) % 360

struct Planet <: Entity
    id::Int
    x::Float64
    y::Float64
    health::Int
    radius::Float64
    num_docking_spots::Int
    current_production::Int
    owned::Bool
    owner_id::Int
    docked_ships_ids::Vector{Int}
end
isowned(planet::Planet) = planet.owned
isfull(planet::Planet) = length(planet.docked_ships_ids) >= planet.num_docking_spots

struct Position <: Entity
    x::Float64
    y::Float64
end
radius(p::Position) = 0
id(p::Position) = -1

############################################################################################
# mobile entities
############################################################################################
abstract type MobileEntity <: Entity end

struct Ship <: MobileEntity
    owner_id::Int
    id::Int
    x::Float64
    y::Float64
    hp::Int
    docked::DockedStatus
    docked_planet_id::Int
    progress::Int
    cooldown::Int
end
radius(::Ship) = SHIP_RADIUS
isdocked(ship::Ship) = ship.docked != UNDOCKED
can_dock(s::Ship, p::Planet) = distance_between(s, p) <= radius(p) + DOCK_RADIUS + radius(s)

# struct Group <: MobileEntity
#     x::Float64
#     y::Float64
#     ships::Vector{Ship}
# end
# radius(::Group) = 3*SHIP_RADIUS+0.1

############################################################################################
# bookkeeping structs
############################################################################################
struct Player
    id::Int
    ships::Dict{Int, Ship}
end
Player(id::Int) = Player(id, Dict())
get_ship(player::Player, ship_id::Int) = player.ships[ship_id]
all_ships(player::Player) = values(player.ships)

struct GameMap
    id::Int
    width::Int
    height::Int
    players::Dict{Int, Player}
    planets::Dict{Int, Planet}
end
get_player(game_map::GameMap, player_id::Int) = game_map.players[player_id]
get_player(game_map::GameMap, entity::Entity) = get_player(game_map, entity.owner_id)
get_me(game_map::GameMap) = get_player(game_map, game_map.id)
all_players(game_map::GameMap) = values(game_map.players)
all_planets(game_map::GameMap) = values(game_map.planets)

get_planet(game_map::GameMap, planet_id::Int) = game_map.planets[planet_id]
get_planet(game_map::GameMap, ship::Ship) = get_planet(game_map, ship.docked_planet_id)
function all_docked_ships(game_map::GameMap, planet::Planet)
    planet.owned && return get_ship.(get_player(game_map, planet), planet.docked_ships_ids)
    return Vector{Ship}()
end
get_docked_ship(game_map::GameMap, planet::Planet, ship_id::Int) = get_ship(get_player(game_map, planet), ship_id)
function all_ships(game_map::GameMap)
    ships = Vector{Ship}()
    for player in all_players(game_map)
        append!(ships, all_ships(player))
    end
    return ships
end

struct Game
    id::Int
    botname::String
    width::Int
    height::Int
    initial_game_map::GameMap
end
function Game(botname::String)
    map = parse_map!()
    Game(map.id, "$botname-$(map.id)", map.width, map.height, map)
end

max_stars_repo_name: Cerbes/Halite-II

max_issues_repo_name: Cerbes/Halite-II

max_forks_repo_name: Cerbes/Halite-II


text:
using IPNets, Test, Sockets

@testset "IPNets" begin
    #############
    ## IPv4Net ##
    #############

    ## Constructors
    @test IPv4Net("1.2.3.4") == parse(IPv4Net, "1.2.3.4") ==
          IPv4Net("1.2.3.4/32") == parse(IPv4Net, "1.2.3.4/32") ==
          IPv4Net(ip"1.2.3.4") == IPv4Net(ip"1.2.3.4", 32) ==
          IPv4Net(ip"1.2.3.4", ip"255.255.255.255") ==
          IPNet("1.2.3.4") == IPNet("1.2.3.4/32") ==
          parse(IPv4Net, SubString("1.2.3.4")) ==
          parse(IPNet, "1.2.3.4") == parse(IPNet, "1.2.3.4/32") ==
          IPv4Net(0x01020304, typemax(UInt32))
    @test IPv4Net("1.2.3.0/24") == parse(IPv4Net, "1.2.3.0/24") ==
          IPv4Net(ip"1.2.3.0", ip"255.255.255.0") ==
          IPv4Net(ip"1.2.3.0", 24) == parse(IPNet, "1.2.3.0/24") ==
          IPv4Net(0x01020300, typemax(UInt32) << 8)

    err = ArgumentError("network mask bits must be between 0 and 32, got 33")
    @test_throws err IPv4Net("1.2.3.4/33")
    @test_throws err IPv4Net(ip"1.2.3.4", 33)
    err = ArgumentError("non-contiguous IPv4 subnets not supported, got 255.240.255.0")
    @test_throws err IPv4Net(ip"1.2.3.0", ip"255.240.255.0")
    err = ArgumentError("input 1.2.3.4/24 has host bits set")
    @test_throws err IPv4Net("1.2.3.4/24")
    @test_throws err parse(IPv4Net, "1.2.3.4/24")
    @test_throws err IPv4Net(ip"1.2.3.4", 24)
    @test_throws err parse(IPNet, "1.2.3.4/24")
    err = ArgumentError("malformed IPNet input: 1.2.3.4/32/32")
    @test_throws err IPv4Net("1.2.3.4/32/32")

    ## Print
    ipnet = IPv4Net("1.2.3.0/24")
    @test sprint(print, ipnet) == "1.2.3.0/24"
    @test sprint(show, ipnet) == "IPv4Net(\"1.2.3.0/24\")"

    ## IPNet as collection
    ipnet = IPv4Net("1.2.3.0/24")

    @test ip"1.2.3.4" in ipnet
    @test ip"1.2.3.0" in ipnet
    @test ip"1.2.3.255" in ipnet
    @test !(ip"1.2.4.0" in ipnet)

    @test IPv4Net("1.2.3.4/32") == IPv4Net("1.2.3.4/32")
    @test IPv4Net("1.2.3.0/24") == IPv4Net("1.2.3.0/24")
    @test IPv4Net("1.2.3.4/32") != IPv4Net("1.2.3.4/31")
    @test IPv4Net("1.2.3.4/31") < IPv4Net("1.2.3.4/32")
    @test IPv4Net("1.2.3.4/32") > IPv4Net("1.2.3.4/31")
    @test IPv4Net("1.2.3.0/24") < IPv4Net("1.2.4.0/24")
    @test IPv4Net("1.2.4.0/24") > IPv4Net("1.2.3.0/24")
    nets = map(IPv4Net, ["1.2.3.0/24", "1.2.3.4/31", "1.2.3.4/32", "1.2.4.0/24"])
    @test sort(nets) == nets

    @test length(ipnet) == length(collect(ipnet)) == 256
    @test collect(ipnet) == [x for x in ipnet]
    @test length(IPv4Net("0.0.0.0/0"))::Int64 == Int64(1) << 32
    @test ipnet[0] == #= ipnet[begin] == =# ip"1.2.3.0" # TODO: Requires Julia 1.4
    @test ipnet[1] == ip"1.2.3.1"
    @test ipnet[255] == ipnet[end] == ip"1.2.3.255"
    @test ipnet[0:1] == ipnet[[0, 1]] == [ip"1.2.3.0", ip"1.2.3.1"]
    @test_throws BoundsError ipnet[-1]
    @test_throws BoundsError ipnet[256]
    @test_throws BoundsError ipnet[-1:2]

    @test is_private(ip"127.0.0.1")
    @test !is_global(ip"127.0.0.1")
    @test !is_private(ip"1.1.1.1")
    @test is_global(ip"1.1.1.1")

    #############
    ## IPv6Net ##
    #############

    ## Constructors
    @test IPv6Net("1:2::3:4") == parse(IPv6Net, "1:2::3:4") ==
          IPv6Net("1:2::3:4/128") == parse(IPv6Net, "1:2::3:4/128") ==
          IPv6Net(ip"1:2::3:4") == IPv6Net(ip"1:2::3:4", 128) ==
          parse(IPv6Net, SubString("1:2::3:4")) ==
          parse(IPNet, "1:2::3:4") == parse(IPNet, "1:2::3:4/128") ==
          IPv6Net(0x00010002000000000000000000030004, typemax(UInt128))
    @test IPv6Net("1:2::3:0/112") == parse(IPv6Net, "1:2::3:0/112") ==
          IPv6Net(ip"1:2::3:0", 112) == parse(IPNet, "1:2::3:0/112") ==
          IPv6Net(0x00010002000000000000000000030000, typemax(UInt128) << 16)

    err = ArgumentError("network mask bits must be between 0 and 128, got 129")
    @test_throws err IPv6Net("1:2::3:4/129")
    @test_throws err IPv6Net(ip"1:2::3:4", 129)
    err = ArgumentError("input 1:2::3:4/112 has host bits set")
    @test_throws err IPv6Net("1:2::3:4/112")
    @test_throws err parse(IPv6Net, "1:2::3:4/112")
    @test_throws err IPv6Net(ip"1:2::3:4", 112)
    @test_throws err parse(IPNet, "1:2::3:4/112")
    err = ArgumentError("malformed IPNet input: 1:2::3:4/32/32")
    @test_throws err IPv6Net("1:2::3:4/32/32")

    ## Print
    ipnet = IPv6Net("1:2::3:0/112")
    @test sprint(print, ipnet) == "1:2::3:0/112"
    @test sprint(show, ipnet) == "IPv6Net(\"1:2::3:0/112\")"

    ## IPNet as collection
    ipnet = IPv6Net("1:2::3:0/112")

    @test ip"1:2::3:4" in ipnet
    @test ip"1:2::3:0" in ipnet
    @test ip"1:2::3:ffff" in ipnet
    @test !(ip"1:2::4:0" in ipnet)

    @test IPv6Net("1:2::3:4/128") == IPv6Net("1:2::3:4/128")
    @test IPv6Net("1:2::3:0/112") == IPv6Net("1:2::3:0/112")
    @test IPv6Net("1:2::3:4/128") != IPv6Net("1:2::3:4/127")
    @test IPv6Net("1:2::3:4/127") < IPv6Net("1:2::3:4/128")
    @test IPv6Net("1:2::3:4/128") > IPv6Net("1:2::3:4/127")
    @test IPv6Net("1:2::3:0/112") < IPv6Net("1:2::4:0/112")
    @test IPv6Net("1:2::4:0/112") > IPv6Net("1:2::3:0/112")
    nets = map(IPv6Net, ["1:2::3:0/112", "1:2::3:4/127", "1:2::3:4/128", "1:2::4:0/112"])
    @test sort(nets) == nets

    @test length(ipnet) == length(collect(ipnet)) == 65536
    @test collect(ipnet)::Vector{IPv6} == [x for x in ipnet]
    @test length(IPv6Net("::/0"))::BigInt == BigInt(1) << 128
    @test ipnet[0] == #= ipnet[begin] == =# ip"1:2::3:0" # TODO: Requires Julia 1.4
    @test ipnet[1] == ip"1:2::3:1"
    @test ipnet[65535] == ipnet[end] == ip"1:2::3:ffff"
    @test ipnet[0:1] == ipnet[[0, 1]] == [ip"1:2::3:0", ip"1:2::3:1"]
    @test_throws BoundsError ipnet[-1]
    @test_throws BoundsError ipnet[65536]
    @test_throws BoundsError ipnet[-1:2]

    @test is_private(ip"::1")
    @test !is_global(ip"::1")
    @test !is_private(ip"2606:4700:4700::1111")
    @test is_global(ip"2606:4700:4700::1111")
end


text:
# Direct, naive translation, works only for small (s, l, m), does not work at poles
function sYlm0(::Val{s}, ::Val{l}, ::Val{m}, θ::Real, ϕ::Real) where {s,l,m}
    T = typeof(zero(θ) * zero(ϕ))
    return bitsign(m) *
           sqrt(factorial(l + m) * factorial(l - m) * (2l + 1) / (factorial(l + s) * factorial(l - s) * 4 * T(π))) *
           sin(θ / 2)^2l *
           sum(
               binomial(l - s, r) * binomial(l + s, r + s - m) * bitsign(l - r - s) * cis(m * ϕ) * cot(θ / 2)^(2r + s - m) for
               r in 0:(l - s)
           )
end

export sYlm
@generated function sYlm(::Val{s}, ::Val{l}, ::Val{m}, θ::Real, ϕ::Real) where {s,l,m}
    stmts = []
    push!(stmts, :((sinθ2, cosθ2) = sincos(θ / 2)))
    T = typeof(one(θ) / one(θ))
    quot = (factorial(big(l + m)) * factorial(big(l - m)))//(factorial(big(l + s)) * factorial(big(l - s)))
    scale = bitsign(m) * sqrt(T(quot * (2l + 1)) / (4 * T(π)))
    push!(stmts, :(res = zero($T)))
    # TODO: Use Horner scheme. Each summand differs by cos(θ/2)^2 / sin(θ/2)^2 from the next.
    for r in 0:(l - s)
        α = bitsign(l - r - s) * binomial(big(l - s), r) * binomial(big(l + s), r + s - m)
        if α ≠ 0
            cospower = 2r + s - m
            sinpower = 2l - cospower
            push!(stmts, :(res += $(T(α)) * sinθ2^$sinpower * cosθ2^$cospower))
        end
    end
    push!(stmts, :(return $scale * res * cis(m * ϕ)))
    return quote
        $((stmts)...)
    end
end

sYlm(s::Int, l::Int, m::Int, θ::Real, ϕ::Real) = sYlm(Val(s), Val(l), Val(m), θ, ϕ)
sYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real) = sYlm(Int(s), Int(l), Int(m), θ, ϕ)


text:
module DLM

export influence_matrix, influence_matrix!
export pressure_coefficients
export AbstractKernel, ParabolicKernel, QuarticKernel

# coefficients for 12-term approximation of kernel integrals
const a = (0.000319759140, -0.000055461471, 0.002726074362, 0.005749551566,
    0.031455895072, 0.106031126212, 0.406838011567, 0.798112357155,
    -0.417749229098, 0.077480713894, -0.012677284771, 0.001787032960)
const b = 0.009054814793

# kernel approximations
"""
    AbstractKernel()

Abstract type for approximation of the planar and nonplanar influence terms in
the numerator of the doublet lattice method kernel
"""
abstract type AbstractKernel end

"""
    ParabolicKernel()

Parabolic approximation for the planar and nonplanar influence terms in the
numerator of the doublet lattice method kernel.
"""
struct ParabolicKernel <: AbstractKernel end

"""
    QuarticKernel()

Quartic approximation for the planar and nonplanar influence terms in the
numerator of the doublet lattice method kernel.
"""
struct QuarticKernel <: AbstractKernel end

"""
    approximate_kernel_integrals(u, k)

Solves for the integrals I0 and J0 using a 12-term exponential approximation
given by Desmarais in "AN ACCURATE METHOD FOR EVALUATING THE KERNEL OF THE
INTEGRAL EQUATION RELATING LIFT TO DOWNWASH IN UNSTEADY POTENTIAL FLOW".

# Arguments
- u = (M*R - x0)/(β^2*x0)
- k = ω*r/U
"""
function approximate_kernel_integrals(u, k)

    TF = promote_type(typeof(u), typeof(k))

    I0 = zero(complex(TF))
    J0 = zero(complex(TF))
    @inbounds for j = 1:12
        p = b*2^j
        pk = complex(p, k)
        tmp = a[j]*exp(-pk*u)/pk
        I0 += tmp
        J0 += tmp*(u+1/pk)
    end

    return I0, J0
end

"""
    kernel_coefficients(x0, y0, z0, ω, U, M)

Solves for the kernel coefficients K1 and K2 and their steady state counterparts
K10 and K20.

See "A Doublet-Lattice Method for Calculating Lift Distributions
on Oscillating Surfaces in Subsonic Flows" by Edward Albano and William P Rodden.

# Arguments
- x0 = x-ξ
- y0 = y-η
- z0 = z-ζ
- ω = frequency
- U = Velocity
- M = Mach Number
"""
function kernel_coefficients(x0, y0, z0, ω, U, M)

    r = sqrt(y0^2 + z0^2)
    β = sqrt(1 - M^2)
    R = sqrt(x0^2+β^2*r^2)

    u = (M*R-x0)/(β^2*max(r, eps(typeof(r))))
    k = ω*r/U

    # use the symmetry properties of the integrand if u < 0
    if u < 0
        I0, J0 = approximate_kernel_integrals(0, k)
        I10 = 1 - 1im*k*I0
        I20 = 2 - 1im*k*I0 + k^2*J0
        u1pos = -u
        I0, J0 = approximate_kernel_integrals(u1pos, k)
        invR = 1/R
        invsqrt = 1/sqrt(1+u1pos^2)
        expk = exp(-1im*k*u1pos)
        I11 = (1-u1pos*invsqrt)*expk - 1im*k*I0
        I21 = ((2 + 1im*k*u1pos)*(1 - u1pos*invsqrt) -
            u1pos*invsqrt^3)*expk - 1im*k*I0 + J0*k^2
        I1 = complex(2*real(I10) - real(I11), imag(I11))
        I2 = complex(2*real(I20) - real(I21), imag(I21))
        expk = exp(-1im*k*u)
    else
        I0, J0 = approximate_kernel_integrals(u, k)
        invR = 1/R
        invsqrt = 1/sqrt(1+u^2)
        expk = exp(-1im*k*u)
        I1 = (1-u*invsqrt)*expk - 1im*k*I0
        I2 = ((2 + 1im*k*u)*(1 - u*invsqrt) -
            u*invsqrt^3)*expk - 1im*k*I0 + J0*k^2
    end

    K1 = I1 + M*r*invR*invsqrt*expk
    K2 = -I2 - 1im*k*invsqrt*expk*(M*r*invR)^2 -
       (M*r*invR)*((1 + u^2)*(β*r*invR)^2 +
       2 + M*r*u*invR)*expk*invsqrt^3

    K10 =  1 + x0/R
    K20 = -2 - x0/R*(2+β^2*r^2/R^2)

    return K1, K2, K10, K20
end

"""
    kernel_numerator(ω, U, M, x0, y0, z0, γr, γs)

Solves for the numerator terms of the kernel and subtracts out the steady-state
components of these terms.

See "Further Refinement of the Subsonic Doublet-Lattice Method" by William P. Rodden,
Paul F. Taylor and Samuel C. McIntosh Jr.
"""
function kernel_numerator(ω, U, M, x0, y0, z0, cγr, sγr, cγs, sγs)
    T1, T2 = cγr*cγs + sγr*sγs, (z0*cγr-y0*sγr)*(z0*cγs-y0*sγs)
    K1, K2, K10, K20 = kernel_coefficients(x0, y0, z0, ω, U, M)
    expn = exp(complex(0, -ω*x0/U))
    return (K1*expn-K10)*T1, (K2*expn-K20)*T2
end

"""
influence_coefficient(ω, U, M, (xr, yr, zr), (xi, yi, zi), (xo, yo, zo),
    chord, λ, γr, γs, kernel)

Solves for the influence coefficient of the sending panel on the receiving panel.

The planar and nonplanar influence terms in the numerator are approximated using
the approximation specified by `kernel`.

See "Further Refinement of the Subsonic Doublet-Lattice Method" by William P. Rodden,
Paul F. Taylor and Samuel C. McIntosh Jr.
"""
influence_coefficient

function influence_coefficient(ω, U, M, (xr, yr, zr),
    (xi, yi, zi), (xo, yo, zo), chord, cλ, cγr, sγr, cγs, sγs, ::ParabolicKernel)

    dx = xo-xi
    dy = yo-yi
    dz = zo-zi
    l = sqrt(dx^2 + dy^2 + dz^2)

    xm = (xo+xi)/2
    ym = (yo+yi)/2
    zm = (zo+zi)/2

    # Approximate numerator coefficients as parabolas
    x0i = xr - xi
    y0i = yr - yi
    z0i = zr - zi
    Ki1, Ki2 = kernel_numerator(ω, U, M, x0i, y0i, z0i, cγr, sγr, cγs, sγs)

    x0m = xr - xm
    y0m = yr - ym
    z0m = zr - zm

    Km1, Km2 = kernel_numerator(ω, U, M, x0m, y0m, z0m, cγr, sγr, cγs, sγs)

    x0o = xr - xo
    y0o = yr - yo
    z0o = zr - zo
    Ko1, Ko2 = kernel_numerator(ω, U, M, x0o, y0o, z0o, cγr, sγr, cγs, sγs)

    e = 1/2*l*cλ
    η = y0m*cγs + z0m*sγs
    ζ = -y0m*sγs + z0m*cγs

    if abs(ζ) == zero(ζ)
        F = 2*e/(η^2-e^2)
    else
        F =atan(2*e*abs(ζ),(η^2+ζ^2-e^2))/abs(ζ)
    end

    # steady term: from VLM
    d0 = steady_influence_coefficient(M, x0i, y0i, z0i, x0o, y0o, z0o, sγr, cγr)

    # planar term: ∫(A1*y^2 + B1*y + C1)/((η-y)^2+ζ^2)
    A1 = (Ki1 - 2*Km1 + Ko1)/(2*e^2)
    B1 = (Ko1 - Ki1)/(2*e)
    C1 = Km1

    d1 = ((η^2-ζ^2)*A1 + η*B1 + C1) * F +
        (1/2*B1 + η*A1) * log(((η-e)^2+ζ^2)/((η+e)^2 + ζ^2)) + 2*e*A1

    # nonplanar term: ∫(A1*y^2 + B1*y + C1)/((η-y)^2+ζ^2)^2
    if abs(ζ) == zero(ζ)
        D2 = zero(complex(typeof(d1)))
    else
        A2 = (Ki2 - 2*Km2 + Ko2)/(2*e^2)
        B2 = (Ko2 - Ki2)/(2*e)
        C2 = Km2

        α = (e/ζ)^2 * (1 - (η^2 + ζ^2 - e^2)/(2*e)*F)

        D2 = e/(η^2 + ζ^2 - e^2)*((
            (2*(η^2 + ζ^2 + e^2)*(e^2*A2 + C2) + 4*η*e^2*B2))/
            (((η + e)^2 + ζ^2)*((η - e)^2 + ζ^2)) -
             (α/e^2)*((η^2 + ζ^2)*A2 + η*B2 + C2))
    end

    return chord/(8*pi)*(d0 + d1 + D2)
end

function influence_coefficient(ω, U, M, (xr, yr, zr),
    (xi, yi, zi), (xo, yo, zo), chord, cλ, cγr, sγr, cγs, sγs, ::QuarticKernel)

    dx = xo-xi
    dy = yo-yi
    dz = zo-zi
    l = sqrt(dx^2 + dy^2 + dz^2)

    xm = (xi+xo)/2
    ym = (yi+yo)/2
    zm = (zi+zo)/2

    xmi = (xi+xm)/2
    ymi = (yi+ym)/2
    zmi = (zi+zm)/2

    xmo = (xm+xo)/2
    ymo = (ym+yo)/2
    zmo = (zm+zo)/2

    # Approximate numerator coefficients as quartic
    x0i = xr - xi
    y0i = yr - yi
    z0i = zr - zi
    Ki1, Ki2 = kernel_numerator(ω, U, M, x0i, y0i, z0i, cγr, sγr, cγs, sγs)

    x0mi = xr - xmi
    y0mi = yr - ymi
    z0mi = zr - zmi

    Kmi1, Kmi2 = kernel_numerator(ω, U, M, x0mi, y0mi, z0mi, cγr, sγr, cγs, sγs)

    x0m = xr - xm
    y0m = yr - ym
    z0m = zr - zm

    Km1, Km2 = kernel_numerator(ω, U, M, x0m, y0m, z0m, cγr, sγr, cγs, sγs)

    x0mo = xr - xmo
    y0mo = yr - ymo
    z0mo = zr - zmo

    Kmo1, Kmo2 = kernel_numerator(ω, U, M, x0mo, y0mo, z0mo, cγr, sγr, cγs, sγs)

    x0o = xr - xo
    y0o = yr - yo
    z0o = zr - zo
    Ko1, Ko2 = kernel_numerator(ω, U, M, x0o, y0o, z0o, cγr, sγr, cγs, sγs)

    e = 1/2*l*cλ
    η = y0m*cγs + z0m*sγs
    ζ = -y0m*sγs + z0m*cγs

    if abs(ζ) == zero(ζ)
        F = 2*e/(η^2-e^2)
    else
        F =atan(2*e*abs(ζ), (η^2+ζ^2-e^2))/abs(ζ)
    end

    # steady term: from VLM
    d0 = steady_influence_coefficient(M, x0i, y0i, z0i, x0o, y0o, z0o, sγr, cγr)

    # planar term: ∫(A1*y^4 + B1*y^3 + C1*y^2 + D1*y + E1)/((η-y)^2+ζ^2)
    A1 = -(1/(6*e^2))*(Ki1 - 16*Kmi1 + 30*Km1 - 16*Kmo1 + Ko1)
    B1 =  (1/(6*e  ))*(Ki1 -  8*Kmi1          +  8*Kmo1 - Ko1)
    C1 = Km1
    D1 = -(2/(3*e^3))*(Ki1 -  2*Kmi1          +  2*Kmo1 - Ko1)
    E1 =  (2/(3*e^4))*(Ki1 -  4*Kmi1 +  6*Km1 -  4*Kmo1 + Ko1)

    d1 = ((η^2-ζ^2)*A1 + η*B1 + C1 + η*(η^2 - 3*ζ^2)*D1 +
        (η^4 - 6*η^2*ζ^2 + ζ^4)*E1)*F + (η*A1 + 1/2*B1 +
        1/2*(3*η^2 - ζ^2)*D1 + 2*η*(η^2-ζ^2)*E1)*log(((η - e)^2 + ζ^2)/((η + e)^2 + ζ^2)) +
        2*e*(A1 + 2*η*D1 + (3*η^2 - ζ^2 + 1/3*e^2)*E1)

    # nonplanar term: ∫(A2*y^4 + B2*y^3 + C2*y^2 + D2*y + E2)/((η-y)^2+ζ^2)^2
    if abs(ζ) == zero(ζ)
        d2 = zero(complex(typeof(d1)))
    else
        A2 = -(1/(6*e^2))*(Ki2 - 16*Kmi2 + 30*Km2 - 16*Kmo2 + Ko2)
        B2 =  (1/(6*e  ))*(Ki2 -  8*Kmi2          +  8*Kmo2 - Ko2)
        C2 = Km2
        D2 = -(2/(3*e^3))*(Ki2 -  2*Kmi2          +  2*Kmo2 - Ko2)
        E2 =  (2/(3*e^4))*(Ki2 -  4*Kmi2 +  6*Km2 -  4*Kmo2 + Ko2)

        d2 = 1/(2*ζ^2)*(((η^2 + ζ^2)*A2 + η*B2 + C2 + η*(η^2 + 3*ζ^2)*D2 +
            (η^4 + 6*η^2*ζ^2 - 3*ζ^4)*E2)*F + 1/((η + e)^2 + ζ^2)*(((η^2 + ζ^2)*η +
            (η^2 - ζ^2)*e)*A2 + (η^2 + ζ^2 + η*e)*B2 + (η + e)*C2 +
            (η^4 - ζ^4 + (η^2 - 3*ζ^2)*η*e)*D2 + ((η^4 - 2*η^2*ζ^2 - 3*ζ^4)*η +
            (η^4 - 6*η^2*ζ^2 + ζ^4)*e)*E2) - 1/((η - e)^2 + ζ^2)*(((η^2 + ζ^2)*η -
            (η^2 - ζ^2)*e)*A2 + (η^2 + ζ^2 - η*e)*B2 + (η - e)*C2 +
            (η^4 - ζ^4 - (η^2 - 3*ζ^2)*η*e)*D2 + ((η^4 - 2*η^2*ζ^2 - 3*ζ^4)*η -
            (η^4 - 6*η^2*ζ^2 + ζ^4)*e)*E2) + (ζ^2*log(((η - e)^2 + ζ^2)/((η+e)^2 + ζ^2)))*D2 +
            4*ζ^2*(e + η*log(((η - e)^2 + ζ^2)/((η + e)^2 + ζ^2)))*E2)
    end

    return chord/(8*pi)*(d0 + d1 + d2)
end

"""
    steady_influence_coefficient(M, x0i, y0i, z0i, x0o, y0o, z0o, sγr, cγr)

Solves for the steady component of the influence of the sending panel on the
receiving panel using the vortex lattice method.
"""
function steady_influence_coefficient(M, x0i, y0i, z0i, x0o, y0o, z0o, sγr, cγr)

    # steady term (from VLM)
    ax = x0i/sqrt(1-M^2)
    ay = y0i
    az = z0i
    anrm = sqrt(ax^2 + ay^2 + az^2)
    ainv = 1/(anrm*(anrm - ax))

    bx = x0o/sqrt(1-M^2)
    by = y0o
    bz = z0o
    bnrm = sqrt(bx^2 + by^2 + bz^2)
    binv = 1/(bnrm*(bnrm - bx))

    vy =  az*ainv - bz*binv
    vz = -ay*ainv + by*binv

    ainv = 1/(anrm*bnrm*(anrm*bnrm + ax*bx + ay*by + az*bz))
    vy = vy + (az*bx - ax*bz)*(anrm + bnrm)*ainv
    vz = vz + (ax*by - ay*bx)*(anrm + bnrm)*ainv

    return -(sγr*vy - cγr*vz)
end

"""
    influence_matrix(ω, U, M, xyz::AbstractVector{<:AbstractArray{<:Number, 3}}, symmetric::AbstractVector; kernel=QuarticKernel())
    influence_matrix(ω, U, M, xyz::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())
    influence_matrix(ω, U, M, xyzr::AbstractArray{<:Number, 3}, xyzs::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())

Fills in aerodynamic influence coefficient matrix. Returns aerodynamic influence
coefficient matrix AIC with shape (nir, njr, nis, njs). For a non-allocating version,
see `influence_matrix!`.

Receiving and sending panel points must be aligned with the x-direction. The
freestream velocity must also be aligned in the x-direction.

# Arguments
- `ω`: oscillation frequency
- `U`: freestream velocity
- `M`: Mach number (for Prandtl-Glauert compressibility correction)
- `xyzr`: definition of receiving panel points with shape: (3, nir, njr)
- `xyzs`: definition of sending panel points with shape: (3, nis, njs)
- `symmetric`: flag indicating symmetry of each sending panel
- `kernel`: indicates which approximation of the kernel should be used
"""
influence_matrix

function influence_matrix(ω, U, M, xyz::AbstractVector{<:AbstractArray{<:Number, 3}},
        symmetric::AbstractVector; kernel=QuarticKernel())

    nx = sum(size.(xyz, 2))
    ny = sum(size.(xyz, 3))

    TF = promote_type(typeof(ω), typeof(U), typeof(M), eltype(eltype(xyz)))

    AIC = zeros(complex(TF), nx, ny, nx, ny)

    return influence_matrix!(AIC, ω, U, M, xyz, symmetric, kernel=kernel)
end

function influence_matrix(ω, U, M, xyzr::AbstractArray{<:Number, 3},
    xyzs::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())

    nir = size(xyzr, 2)-1
    njr = size(xyzr, 3)-1
    nis = size(xyzs, 2)-1
    njs = size(xyzs, 3)-1

    TF = promote_type(typeof(ω), typeof(U), typeof(M), eltype(xyzr), eltype(xyzs))

    AIC = zeros(complex(TF), nir, njr, nis, njs)

    return influence_matrix!(AIC, ω, U, M, xyzr, xyzs, symmetric, kernel=kernel)
end

influence_matrix(ω, U, M, xyz::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel()) = influence_matrix(ω, U, M, xyz, xyz, symmetric, kernel=kernel)

"""
    influence_matrix!(AIC, ω, U, M, xyz::AbstractVector{<:AbstractArray{<:Number, 3}}, symmetric::AbstractVector; kernel=QuarticKernel())
    influence_matrix!(AIC, ω, U, M, xyz::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())
    influence_matrix!(AIC, ω, U, M, xyzr::AbstractArray{<:Number, 3}, xyzs::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())

Non-allocating version of [`influence_matrix`](@ref)
"""
influence_matrix!

function influence_matrix!(AIC, ω, U, M, xyz::AbstractVector{<:AbstractArray{<:Number, 3}}, symmetric::AbstractVector; kernel=QuarticKernel())

    nx = size(AIC, 1)
    ny = size(AIC, 2)
    nsurf = length(xyz)

    ir = 0
    jr = 0
    is = 0
    js = 0
    @inbounds for kr = 1:nsurf, ks = 1:nsurf
        nir = size(xyz[ir], 2)
        njr = size(xyz[ir], 3)
        nis = size(xyz[is], 2)
        njs = size(xyz[is], 3)

        ixr = ir+1:ir+nir
        iyr = jr+1:jr+njr
        ixs = is+1:is+nis
        iys = js+1:js+njs

        influence_matrix!(view(AIC, ixr, iyr, ixs, iys), ω, U, M, xyz[kr], xyz[ks], symmetric[ks], kernel=kernel)
    end

    return AIC
end

function influence_matrix!(AIC, ω, U, M, xyzr::AbstractArray{<:Number, 3}, xyzs::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel())
    nir = size(xyzr, 2)-1
    njr = size(xyzr, 3)-1
    nis = size(xyzs, 2)-1
    njs = size(xyzs, 3)-1

    @assert size(AIC) == (nir, njr, nis, njs)

    @inbounds for jr=1:njr, ir = 1:nir, js=1:njs, is=1:nis
        # receiving point at 3/4 chord
        xri = (1/4)*xyzr[1, ir, jr] + (3/4)*xyzr[1, ir+1, jr]
        yri = (1/4)*xyzr[2, ir, jr] + (3/4)*xyzr[2, ir+1, jr]
        zri = (1/4)*xyzr[3, ir, jr] + (3/4)*xyzr[3, ir+1, jr]
        xro = (1/4)*xyzr[1, ir, jr+1] + (3/4)*xyzr[1, ir+1, jr+1]
        yro = (1/4)*xyzr[2, ir, jr+1] + (3/4)*xyzr[2, ir+1, jr+1]
        zro = (1/4)*xyzr[3, ir, jr+1] + (3/4)*xyzr[3, ir+1, jr+1]
        xr = (xri+xro)/2
        yr = (yri+yro)/2
        zr = (zri+zro)/2
        # sending line at 1/4 chord
        xsi = (3/4)*xyzs[1, is, js] + (1/4)*xyzs[1, is+1, js]
        ysi = (3/4)*xyzs[2, is, js] + (1/4)*xyzs[2, is+1, js]
        zsi = (3/4)*xyzs[3, is, js] + (1/4)*xyzs[3, is+1, js]
        xso = (3/4)*xyzs[1, is, js+1] + (1/4)*xyzs[1, is+1, js+1]
        yso = (3/4)*xyzs[2, is, js+1] + (1/4)*xyzs[2, is+1, js+1]
        zso = (3/4)*xyzs[3, is, js+1] + (1/4)*xyzs[3, is+1, js+1]
        # average chord length of sending panel
        chord = (xyzs[1, is+1, js] - xyzs[1, is, js])/2 + (xyzs[1, is+1, js+1] - xyzs[1, is, js+1])/2
        # cos(sweep) of sending doublet line
        cλs = (yso-ysi)/(sqrt((xso-xsi)^2 + (yso-ysi)^2))
        # cos(dihedral) and sin(dihedral) of receiving point
        tmp = sqrt((yro-yri)^2+(zro-zri)^2)
        cγr = (yro-yri)/tmp
        sγr = (zro-zri)/tmp
        # cos(dihedral) and sin(dihedral) of sending line
        tmp = sqrt((yso-ysi)^2+(zso-zsi)^2)
        cγs = (yso-ysi)/tmp
        sγs = (zso-zsi)/tmp
        # compute influence of the sending panel on the receiving panel
        AIC[ir, jr, is, js] = influence_coefficient(ω, U, M, (xr, yr, zr),
            (xsi, ysi, zsi), (xso, yso, zso), chord, cλs, cγr, sγr, cγs, sγs, kernel)
        # add influence of reflected sending panel, if applicable
        if symmetric
            AIC[ir, jr, is, js] += influence_coefficient(ω, U, M, (xr, yr, zr),
                (xso, -yso, zso), (xsi, -ysi, zsi), chord, cλs, cγr, sγr, cγs, -sγs, kernel)
        end
    end

    return AIC
end

influence_matrix!(AIC, ω, U, M, xyz::AbstractArray{<:Number, 3}, symmetric; kernel=QuarticKernel()) = influence_matrix!(AIC, ω, U, M, xyz, xyz, symmetric, kernel=kernel)

"""
    pressure_coefficients(AIC, w)

Compute pressure coefficients `Cp` on each panel, given the downwash w.

# Arguments:
- `AIC`: aerodynamic influence coefficient matrix with shape (nir, njr, nis, njs)
- `w`: downwash on each panel, with shape (nir, njr)
"""
function pressure_coefficients(AIC, w)

    nir, njr, nis, njs = size(AIC)

    AIC_r = reshape(AIC, nir*njr, nis*njs)
    w_r = reshape(w, nir*njr)

    Cp = AIC_r\w_r

    return reshape(Cp, nir, njr)
end

# """
#     downwash!(w, xyz, Vinf, α, β)
#
# Compute the downwash on each of the panels in `xyz` given the freestream
# velocity `Vinf`, angle of attack `α`, sideslip angle `β`, and panel twist `θ`.
# Panel twist has the shape (ni, nj). Non-allocating.
# """
# function downwash!(w, xyz, Vinf, α, β, θ)
#
#     ca, sa = cos(α), sin(α)
#     cb, sb = cos(β), sin(β)
#
#     V = (Vinf*ca*cb, -Vinf*sb, Vinf*sa*cb)
#
#     ni = size(xyz, 2)-1
#     nj = size(xyz, 3)-1
#
#     @assert size(w) == (ni, nj)
#
#     @inbounds for j = 1:nj, i = 1:ni
#         # compute cos(dihedral) and sin(dihedral)
#         yri = (1/4)*xyz[2, i, j] + (3/4)*xyz[2, i+1, j]
#         zri = (1/4)*xyz[3, i, j] + (3/4)*xyz[3, i+1, j]
#         yro = (1/4)*xyz[2, i, j+1] + (3/4)*xyz[2, i+1, j+1]
#         zro = (1/4)*xyz[3, i, j+1] + (3/4)*xyz[3, i+1, j+1]
#         tmp = sqrt((yro-yri)^2+(zro-zri)^2)
#         cd = (yro-yri)/tmp
#         sd = (zro-zri)/tmp
#         # compute cos(theta) and sin(theta)
#         ct = cos(θ[i,j])
#         st = sin(θ[i,j])
#         # compute downwash on each panel
#         w[i, j] = V[1]*st + V[2]*(-ct*sd) + V[3]*(ct*cd)
#     end
#     return w
# end

# """
#     downwash(xyz, Vinf, α, β)
#
# Compute the downwash on each of the panels in `xyz` given the freestream
# velocity `Vinf`, angle of attack `α`, sideslip angle `β`, and panel twist `θ`.
# Panel twist has the shape (ni, nj). For a non-allocating version see `downwash!`.
# """
# function downwash(xyz, Vinf, α, β, θ)
#
#     TF = promote_type(eltype(xyz), typeof(Vinf), typeof(α), typeof(β), eltype(θ))
#
#     w = zeros(TF, size(xyz,2)-1, size(xyz,3)-1)
#
#     return downwash!(w, xyz, Vinf, α, β, θ)
# end

end # module


text:
struct TestData
    x_raw
    x_table
    y_raw
    y_cat
    n_train
    n_test
    train_idx
    test_idx

    function TestData(;rng=MersenneTwister(0), fraction_train=0.5, n_samples=100, n_dim=10)
        X_raw = rand(rng, n_dim, n_samples)
        X_table = MLJ.table(collect(X_raw'))
        y_raw = rand(rng, (CLASS_NORMAL, CLASS_OUTLIER), n_samples)
        y_cat = to_categorical(y_raw)
        train, test  = MLJ.partition(eachindex(y_cat), fraction_train, rng=rng);
        n_train = Int(n_samples * fraction_train)
        n_test = n_samples - n_train

        new(X_raw, X_table, y_raw, y_cat, n_train, n_test, train, test)
    end
end


text:
using ImageQuilting
using Base.Test

# the output of a homogeneous image is also homogeneous
TI = ones(20,20,20)
reals = iqsim(TI, 10, 10, 10, size(TI)...)
@test reals[1] == ones(TI)

# categories are obtained from training image only
ncateg = 3; TI = rand(RandomDevice(), 0:ncateg, 20, 20, 20)
reals = iqsim(TI, 10, 10, 10, size(TI)..., categorical=true)
@test Set(reals[1]) ⊆ Set(TI)

# trends with soft data
TI = [zeros(10,20,1); ones(10,20,1)]
trend = [zeros(20,10,1) ones(20,10,1)]
reals = iqsim(TI, 10, 10, 1, size(TI)..., soft=SoftData(trend, x -> x), tol=1)
@test mean(reals[1][:,1:10,:]) ≤ mean(reals[1][:,11:20,:])

# hard data is honored everywhere
TI = ones(20,20,20)
obs = rand(size(TI))
data = HardData([(i,j,k)=>obs[i,j,k] for i=1:20, j=1:20, k=1:20])
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=data)
@test reals[1] == obs

# irregular grids via hard data conditioning
TI = ones(20,20,20)
shape = HardData()
active = trues(TI)
for i=1:20, j=1:20, k=1:20
  if (i-10)^2 + (j-10)^2 + (k-10)^2 < 25
    push!(shape, (i,j,k)=>NaN)
    active[i,j,k] = false
  end
end
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape)
@test all(isnan(reals[1][!active]))
@test all(!isnan(reals[1][active]))

# irregular training images
TI = ones(20,20,20)
TI[:,5,:] = NaN
reals = iqsim(TI, 10, 10, 10, size(TI)...)
@test reals[1] == ones(TI)
TI[1,5] = 0
reals = iqsim(TI, 10, 10, 10, size(TI)..., categorical=true)
@test reals[1] == ones(TI)

# irregular training image and irregular grid
TI = ones(20,20,20)
TI[:,5,:] = NaN
trend = SoftData(ones(TI), x -> ones(TI))
shape = HardData([(i,j,k)=>NaN for i=1:20, j=5, k=1:20])
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape)
@test all(isnan(reals[1][:,5,:]))
@test all(reals[1][:,1:4,:] .== 1)
@test all(reals[1][:,6:20,:] .== 1)
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=shape, soft=trend)
@test all(isnan(reals[1][:,5,:]))
@test all(reals[1][:,1:4,:] .== 1)
@test all(reals[1][:,6:20,:] .== 1)

# no side effects with soft data
TI = ones(20,20,20)
TI[:,5,:] = NaN
aux = ones(TI)
trend = SoftData(aux, x -> aux)
iqsim(TI, 10, 10, 10, size(TI)..., soft=trend)
@test aux == ones(TI)

# raster path and dilation both take place
TI = ones(20,20,20)
data = HardData((1,1,1)=>10)
trend = SoftData(aux, x -> aux)
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=data, soft=trend)
@test reals[1][1,1,1] == 10
@test all(reals[1][2:end] .== 1)

# boundary cut in 3D
TI = ones(20,20,20)
for cut in [:dijkstra,:boykov]
  _, _, voxs = iqsim(TI, 10, 10, 10, size(TI)..., overlapx=1/3, overlapy=1/3, overlapz=1/3, cut=cut, debug=true)
  @test 0 ≤ voxs[1] ≤ 1
end

# mean voxel reuse is in range [0,1]
TI = rand(20,20,20)
@test 0 ≤ meanvoxreuse(TI, 10, 10, 10, nreal=1) ≤ 1

# multiple realizations with hard data
TI = ones(20,20,20)
data = HardData((20,20,20)=>10)
reals = iqsim(TI, 10, 10, 10, size(TI)..., hard=data, nreal=3)
for real in reals
  @test real[20,20,20] == 10
end

# different simulation paths
for kind in [:rasterup,:rasterdown,:dilation,:random]
  path = ImageQuilting.genpath((10,10,10), kind)
  @test length(path) == 1000
end

# datum is visited first if present
path = ImageQuilting.genpath((10,10,10), :datum, [(1,1,1),(10,10,10)])
@test path[1] == 1
@test path[2] == 1000


text:
module TestMulitArg

using Test
using ThreadPools

include("util.jl")

@testset "multiarg" begin

    @testset "tmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data + y
        end
        @test tmap(fn!, (TestObj(x) for x in 1:N), adders) == collect(1+1:N+1)
        @test primary
    end

    @testset "bmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data + y
        end
        @test bmap(fn!, (TestObj(x) for x in 1:N), adders) == collect(1+1:N+1)
        @inferred bmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "qmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data + y
        end
        @test qmap(fn!, (TestObj(x) for x in 1:N), adders) == collect(1+1:N+1)
        @test primary
        @inferred qmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "qbmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data + y
        end
        @test qbmap(fn!, (TestObj(x) for x in 1:N), adders) == collect(1+1:N+1)
        @inferred qbmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "logtmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data + y
        end
        p, r = logtmap(fn!, (TestObj(x) for x in 1:N), adders)
        @test r == collect(1+1:N+1)
        @test length(p.recs) == N*2
        @test primary
        @inferred logtmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "logbmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data + y
        end
        p, r = logbmap(fn!, (TestObj(x) for x in 1:N), adders)
        @test r == collect(1+1:N+1)
        @test length(p.recs) == N*2
        @inferred logbmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "logqmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data + y
        end
        p, r = logqmap(fn!, (TestObj(x) for x in 1:N), adders)
        @test r == collect(1+1:N+1)
        @test length(p.recs) == N*2
        @test primary
        @inferred logqmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "logqbmap" begin
        N = 2 * Threads.nthreads()
        adders = [1 for _ in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data + y
        end
        p, r = logqbmap(fn!, (TestObj(x) for x in 1:N), adders)
        @test r == collect(1+1:N+1)
        @test length(p.recs) == N*2
        @inferred logqbmap(fn!, (TestObj(x) for x in 1:N), adders)
    end

    @testset "tforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += y
        end
        tforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test primary
        @inferred tforeach(fn!, objs, 1:N)
    end

    @testset "bforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += y
        end
        bforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @inferred bforeach(fn!, objs, 1:N)
    end

    @testset "qforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += y
        end
        qforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test primary
        @inferred qforeach(fn!, objs, 1:N)
    end

    @testset "qbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += y
        end
        qbforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @inferred qbforeach(fn!, objs, 1:N)
    end

    @testset "logtforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += y
        end
        p = logtforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test length(p.recs) == N*2
        @test primary
        @inferred logtforeach(fn!, objs, 1:N)
    end

    @testset "logbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += y
        end
        p = logbforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test length(p.recs) == N*2
        @inferred logbforeach(fn!, objs, 1:N)
    end

    @testset "logqforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        primary = Threads.nthreads() == 1
        fn! = (x, y) -> begin
            Threads.threadid() == 1 && (primary = true)
            x.data += y
        end
        p = logqforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test length(p.recs) == N*2
        @test primary
        @inferred logqforeach(fn!, objs, 1:N)
    end

    @testset "logqbforeach" begin
        N = 2 * Threads.nthreads()
        objs = [TestObj(x) for x in 1:N]
        fn! = (x, y) -> begin
            Threads.nthreads() == 1 || Threads.threadid() == 1 && error("Task on primary")
            x.data += y
        end
        p = logqbforeach(fn!, objs, 1:N)
        @test [x.data for x in objs] == collect(1:N)*2
        @test length(p.recs) == N*2
        @inferred logqbforeach(fn!, objs, 1:N)
    end

end

end # module


text:
using Walk1DMDP, POMDPToolbox
using Base.Test

mdp = Walk1D()
random_policy = RandomGaussian(MersenneTwister(0))
hr = HistoryRecorder(; max_steps=50, show_progress=true)
h = simulate(hr, mdp, random_policy)
h.state_hist


text:
function superoperator(command, logicset::LogicalCombo; verbose=true)
    counter!()

    logicsetcopy = deepcopy(logicset)

    superset = "xor|iff|if|or|and|notthen|then"

    #println("operatoreval($command)")
    (sum(logicset[:]) == 0) && return logicset
    !occursin(Regex("(\b| )([><=|!+\\-^&]{3}|$superset)(\b| )"), command) &&
      return operatoreval(command, logicset, verbose=verbose)
    occursin(r"\{\{.*\}\}", command) &&
      return operatorspawn(command, logicset, verbose=verbose)

    m = match(Regex("(^.*)(?:\b| )([><=|!+\\-\\^&]{3}|$superset)(?:\b| )(.*?\$)"), command)
    left, operator, right = m.captures

    ℧left  = superoperator(left ,logicset, verbose=verbose)[:]
    ℧right = superoperator(right,logicset, verbose=verbose)[:]

    ℧ = logicset[:]
    ℧η = deepcopy(℧)

    if operator ∈ ["&&&", "and"]
        ℧η = ℧ .& (℧left .& ℧right)

    elseif operator ∈ ["^^^" , "xor", "!=="]
        ℧η = ℧ .& ((℧left .& .!℧right) .| (.!℧left .& ℧right))

    # this can be dangerous, false equal to false such as with previous exclusions will cause inconsistencies
    elseif operator ∈ ["<=>","===", "iff"]
        ℧η = ℧ .& (℧left .== ℧right)

    elseif operator ∈ ["|||", "or"]
        ℧η = ℧ .& (℧left .| ℧right)

    elseif operator ∈ ["==>", "then"]
        ℧η[℧left] .= ℧[℧left]  .& ℧right[℧left]

    elseif operator ∈ ["!=>"]
        ℧η[.!℧left] .= ℧[.!℧left]  .& ℧right[.!℧left]

    elseif operator ∈ ["<==", "if"]
        ℧η[℧right] = ℧[℧right] .& ℧left[℧right]

    end

    logicsetcopy[:] = ℧η
    logicsetcopy
end


text:
# All processes related to survival:
# - growth
# - density-independent mortality
# - competition

"""
    survive!(patch, mortality)

Density independent survival of individuals in a patch. The actual mortality
probability is calculated with a metabolic formula, modified by the passed `mortality`
variable and an individual's temperature adaptation.
"""
function survive!(patch::Patch, mortality::Float64)
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        if !patch.community[idx].marked
            mass = patch.community[idx].size
            deathrate = mortality * mass^(-1/4) * exp(-act/(boltz*temp))
            dieprob = (1 - exp(-deathrate))
            if rand() * patch.community[idx].tempadaptation < dieprob
                splice!(patch.community, idx)
                continue
            end
        end
        idx += 1
    end
end

"""
    survive!(world, settings)

World-wide mortality. Sounds apocalyptic, but is just a fact of life.
"""
function survive!(world::Array{Patch,1}, settings::Dict{String, Any})
    for patch in world
        (patch.isisland || !settings["static"]) && survive!(patch, settings["mortality"]) # pmap(!,patch) ???
    end
end

"""
    grow!(patch, growthrate, capgrowth)

Growth of individuals in the given patch. The actual growthrate is calculated
with a metabolic formula, modified by the passed `growthrate` variable.
"""
function grow!(patch::Patch, growthrate::Float64, capgrowth::Bool)
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        if !patch.community[idx].marked
            repsize = patch.community[idx].traits["repsize"]
            mass = patch.community[idx].size
            if mass <= repsize # stop growth if reached repsize
                growth = growthrate * mass^(3/4) * exp(-act/(boltz*temp))
                newmass = mass + growth
                if capgrowth && newmass > patch.community[idx].traits["repsize"]
                    newmass = patch.community[idx].traits["repsize"]
                end
                if newmass > 0 && mass > 0
                    patch.community[idx].size = newmass
                else
                    splice!(patch.community, idx)
                    continue
                end
            end
        end
        idx += 1
    end
end

"""
    grow!(world, settings)

Carry out growth for all patches.
"""
function grow!(world::Array{Patch,1}, settings::Dict{String, Any})
    for patch in world
        # pmap(!,patch) ???
        if patch.isisland || !settings["static"]
            grow!(patch, settings["growthrate"], settings["capgrowth"])
        end
    end
end

"""
    compete!(patch)

While there are too many organisms in a patch, pick two at random and kill the
one that is less adapted to the local precipitation levels.
"""
function compete!(patch::Patch)
    totalmass = sum(map(x -> x.size, patch.community))
    while totalmass >= patch.capacity # occupied capacity larger than available
        firstind, secondind = rand(eachindex(patch.community), 2)
        firstind == secondind && length(eachindex(patch.community)) > 1 && continue
        if patch.community[firstind].precadaptation < patch.community[secondind].precadaptation
            totalmass -= patch.community[firstind].size
            splice!(patch.community, firstind) # profiling: expensive!
        elseif patch.community[firstind].precadaptation > patch.community[secondind].precadaptation
            totalmass -= patch.community[secondind].size
            splice!(patch.community, secondind) # profiling: expensive!
        else
            victim = rand([firstind, secondind])
            totalmass -= patch.community[victim].size
            splice!(patch.community, victim) # profiling: expensive!
        end
    end
end

"""
    compete!(world, static)

Carry out competition on all patches.
"""
function compete!(world::Array{Patch,1}, static::Bool = true)
    for patch in world
        (patch.isisland || !static) && compete!(patch) # pmap(!,patch) ???
    end
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
s = open("input.txt") do file
    readlines(file);
end

function ParseRule(rule)
    m = match(r"(\w+) (\w+) bags contain", rule);

    color = m.captures[1] * " " * m.captures[2];

    bags = [];
    counts = [];
    for m in eachmatch(r"(\d+) (\w+) (\w+) bag", rule)
        push!(bags, join(m.captures[2:3], " "))
        push!(counts, parse(Int, m.captures[1]))
    end

    return (color, bags, counts);
end

function CreateColorGraph(lines)
    bagGraph = Dict();
    for line in lines
        rule = ParseRule(line);

        bagGraph[rule[1]] = (rule[2], rule[3]);
    end

    return bagGraph;
end

graph = CreateColorGraph(s);

function GetContentCount(color, seen)
    contents = graph[color];
    colors = contents[1];
    counts = contents[2];

    myCount = 0;

    for i = 1:length(colors)
        if haskey(seen, colors[i])
            myCount += counts[i] * (1 + seen[colors[i]]);
        else
            childCount = GetContentCount(colors[i], seen);
            myCount += counts[i] * (1 + childCount);
        end
    end

    seen[color] = myCount;
    return myCount;
end

seen = Dict();
GetContentCount("shiny gold", seen)


text:
using NESSie: yukawa
using NESSie.Radon
using NESSie.Radon: radoncoll!

@testset "regularyukawapot" begin
    for T in testtypes
        # x -> ξ
        x = ones(T, 3)
        yuk = T(7)
        ret = Radon.regularyukawapot(x, x, yuk)
        @test typeof(ret) == T
        @test ret == -7
        # ξ not in origin (no cancellation)
        ξ = -ones(T, 3)
        ret = Radon.regularyukawapot(x, ξ, yuk)
        @test typeof(ret) == T
        @test ret ≈ T(-0.28867513458631466433435706025773931468027671112178301167202)
        # ξ in origin (no cancellation)
        ξ = zeros(T, 3)
        ret = Radon.regularyukawapot(x, ξ, yuk)
        @test typeof(ret) == T
        @test ret ≈ T(-0.57734713663526745894988668695023080532035174827300561882957)
        # ξ in origin (potential cancellation)
        ret = Radon.regularyukawapot(T(.001) * x, ξ, yuk)
        @test typeof(ret) == T
        @test ret ≈ T(-6.95773573664079611429119468507639364064769481076591007615189)
        # ξ in origin (potential cancellation 2)
        ret = Radon.regularyukawapot(T(.0001) * x, ξ, yuk)
        @test typeof(ret) == T
        @test ret ≈ (T == Float64 ?
            -6.99575819000175052904190394524100762244862345273988759731283 :
            -yuk)
        # ξ in origin (potential cancellation 3)
        ret = Radon.regularyukawapot(T(.00001) * x, ξ, yuk)
        @test typeof(ret) == T
        @test ret ≈ (T == Float64 ?
            -6.99957566470162580591945945320355718542304458045960164651420 :
            -yuk)
    end
end

@testset "∂ₙregularyukawapot" begin
    for T in testtypes
        # x -> ξ
        x = ones(T, 3)
        n = map(T, [1, 0, 0])
        yuk = T(7)
        ret = Radon.∂ₙregularyukawapot(x, x, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ T(14.145081595145832)
        # ξ not in origin (no cancellation)
        ξ = -ones(T, 3)
        ret = Radon.∂ₙregularyukawapot(x, ξ, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ T(0.048112522396707305228639137148607498838982054839859863103104)
        # ξ in origin (no cancellation)
        ξ = zeros(T, 3)
        ret = Radon.∂ₙregularyukawapot(x, ξ, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ T(0.192436385477375021033020985138032408973537816345889754008150)
        # ξ in origin (potential cancellation)
        ret = Radon.∂ₙregularyukawapot(T(.001) * x, ξ, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ T(14.03126641709760402264474861719654017225786809145550365963283)
        # ξ in origin (potential cancellation 2)
        ret = Radon.∂ₙregularyukawapot(T(.0001) * x, ξ, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ (T == Float64 ?
            14.13365345844970855427401235115544109979343227826364222026794 :
            yuk^2 / 2 / Float32(√3))
        # ξ in origin (potential cancellation 3)
        ret = Radon.∂ₙregularyukawapot(T(.00001) * x, ξ, yuk, n)
        @test typeof(ret) == T
        @test ret ≈ (T == Float64 ?
            14.14393831379399210175378534648974379367907886588210581085651 :
            yuk^2 / 2 / Float32(√3))
    end
end

@testset "radoncoll!" begin
    for T in testtypes
        elem = [Triangle(zeros(T, 3), T[2, 0, 0], T[0, 2, 0])]
        # compute triangle area ∫dx
        res = zeros(T, 1, 1)
        radoncoll!(res, elem, [zeros(T, 3)], (_, __, ___, ____) -> one(T), zero(T))
        @test res[1] ≈ 2one(T)
        # compute prism volume
        res = zeros(T, 2, 1)
        radoncoll!(res, elem, [T[0, 0, 4], T[0, 0, 6]], (x, ξ, _, __) -> ξ[3] - x[3], zero(T))
        @test res[1] ≈ 8one(T)
        @test res[2] ≈ 12one(T)
        # compute tetrahedron volume
        res = zeros(T, 1, 1)
        radoncoll!(res, elem, [T[0, 0, 2]], (x, ξ, _, __) -> ξ[3] - x[1] - x[2], zero(T))
        @test res[1] ≈ 4one(T)/3
        # a slightly more elaborate tetrahedron volume
        elem = [Triangle(zeros(T, 3), T[3, 0, 0], T[0, 2, 0])]
        res = zeros(T, 1, 1)
        radoncoll!(res, elem, [T[0, 0, 6]], (x, ξ, _, __) -> ξ[3] - 2x[1] - 3x[2], zero(T))
        @test res[1] ≈ 6one(T)
        # ∫∫6x²-40y dA
        elem = [Triangle(T[0, 3, 0], T[1, 1, 0], T[5, 3, 0])]
        res = zeros(T, 1, 1)
        radoncoll!(res, elem, [zeros(T, 3)], (x, _, __, ___) -> 6*x[1]^2-40x[2], zero(T))
        @test res[1] ≈ T(-935/3)
    end
end