text:
# Particle selection routines

struct SelectionCache
	D::Dict{Symbol,Any}
end

SelectionCache() = SelectionCache(Dict{Symbol,Any}())

function Base.get!(f, cache::SelectionCache, groupkey::Symbol, itemkey::Symbol)
	group = get!(cache.D, groupkey) do
		Dict{Symbol,Any}()
	end
	get!(f, group, itemkey)
end

Base.delete!(cache::SelectionCache, groupkey::Symbol) =
		delete!(cache.D, groupkey)

function getchainids(model::ParticleCollection, cache::SelectionCache)
	get(model, :chainids) do
		Repeated("", length(model))
	end
end

function getelements(model::ParticleCollection, cache::SelectionCache)
	get!(cache, :system, :elements) do
		infermissingelements(model)
	end
end

function getmasses(model::ParticleCollection, cache::SelectionCache)
	get(model, :masses) do
		get!(cache, :system, :masses) do
			masses = Vector{Float64}(undef, length(model))
			elements = getelements(model, cache)
			infermasses!(masses, model.names, elements)
		end
	end
end

function getss(model::ParticleCollection, cache::SelectionCache)
	get(model, :SS) do
		get!(cache, :frame, :SS) do
			inferss(model)
		end
	end
end

function getpbccell(model::ParticleCollection, cache::SelectionCache)
	get!(cache, :frame, :cell) do
		pbccell(get(model.header, :cell, nothing))
	end
end

function getpbcpos(R::AbstractVector{Vector3D},
		cell::Union{TriclinicPBC,Nothing}, cache::SelectionCache)
	get!(cache, :frame, :pbcpos) do
		Rw, Kw = get!(cache, :buffers, :pbcpos) do
			similar(R), similar(R)
		end
		pbcpos!(Rw, Kw, R, cell)
	end
end

function getnbrtable(n::Integer, cell::Union{TriclinicPBC,Nothing},
		cache::SelectionCache)
	get!(cache, :frame, :nbrtable) do
		buffer = get!(cache, :buffers, :nbrtable) do
			NeighborTableBuffer()
		end
		NeighborTable(n, cell, buffer)
	end
end

abstract type Selector end

Base.map(s::Selector, model::ParticleCollection, cache = SelectionCache()) =
		map!(s, BitVector(undef, length(model)), model, cache)

function Base.map!(s::Selector, results::BitVector, model::ParticleCollection,
		cache = SelectionCache())
	if ! isempty(model)
		select!(results, eachindex(model), s, model, cache)
	end
	results
end

Base.view(model::ParticleCollection, s::Selector, cache = SelectionCache()) =
		MulticollectionView(model, findall(map(s, model, cache)))

Base.getindex(model::ParticleCollection, s::Selector) = view(model, s)

abstract type SelectionMode end

struct ModelMode <: SelectionMode end

Base.show(io::IO, ::ModelMode) = "model"

struct ChainMode <: SelectionMode end

Base.show(io::IO, ::ChainMode) = "chain"

function getmcrp(model::ParticleCollection, cache::SelectionCache)
	get!(cache, :system, :mcrp) do
		mcrp(model)
	end
end

geth2(model::ParticleCollection, ::ChainMode, cache::SelectionCache) =
		getmcrp(model, cache).flath1

struct ResidueMode <: SelectionMode end

Base.show(io::IO, ::ResidueMode) = "residue"

geth2(model::ParticleCollection, ::ResidueMode, cache::SelectionCache) =
		getmcrp(model, cache).flath2

struct ParticleMode <: SelectionMode end

Base.show(io::IO, ::ParticleMode) = "particle"

struct FragmentMode <: SelectionMode end

Base.show(io::IO, ::FragmentMode) = "fragment"

function gettopology(model::ParticleCollection, cache::SelectionCache)
	get(model, :topology) do
		get!(cache, :system, :topology) do
			infertopology(model)
		end
	end
end

function geth2(model::ParticleCollection, ::FragmentMode, cache::SelectionCache)
	get!(cache, :system, :mfp) do
		mfp(mfptree(gettopology(model, cache)), length(model))
	end
end

iscallable(f) = !isempty(methods(f))

function checkcallable(f)
	if ! iscallable(f)
		error("expected callable")
	end
	f
end

indexfunction(f) = checkcallable(f)

indexfunction(i::Integer) = (J -> i)

indexfunction(I::AbstractVector{<:Integer}) = (J -> I)

positionfunction(f) = checkcallable(f)

positionfunction(R::AbstractVector{<:Real}) = positionfunction(Vector3D(R))

positionfunction(R::Vector3D) = positionfunction([R])

positionfunction(R::AbstractVector{Vector3D}) = (model, cache) -> R

integerpredicate(f) = checkcallable(f)

integerpredicate(i::Integer) = (==(i))

function integerpredicate(I::AbstractVector{<:Integer})
	J = sort(I)
	i -> !isempty(searchsorted(J, i))
end

stringpredicate(f) = checkcallable(f)

stringpredicate(s::AbstractString) = namematcher(s)

stringpredicate(S::AbstractVector{<:AbstractString}) = namematcher(S)

realpredicate(f) = checkcallable(f)

realpredicate(r::Real) = (s -> isapprox(s, r))

vector3dpredicate(f) = checkcallable(f)

vector3dpredicate(V::AbstractVector{<:Real}) = vector3dpredicate(Vector3D(V))

vector3dpredicate(V::Vector3D) = (W -> isapprox(W, V))

struct NotSelector{T<:Selector} <: Selector
	s::T

	NotSelector{T}(s::T) where {T<:Selector} = new(s)
end

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::NotSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s, model, cache)
	for i in subset
		results[i] = !results[i]
	end
end

struct AndSelector{T1<:Selector,T2<:Selector} <: Selector
	s1::T1
	s2::T2

	AndSelector{T1,T2}(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			new(s1, s2)
end

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::AndSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s1, model, cache)
	subset2 = [i for i in subset if results[i]]
	results2 = similar(results)
	select!(results2, subset2, s.s2, model, cache)
	for i in subset2
		if ! results2[i]
			results[i] = false
		end
	end
end

struct OrSelector{T1<:Selector,T2<:Selector} <: Selector
	s1::T1
	s2::T2

	OrSelector{T1,T2}(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			new(s1, s2)
end

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::OrSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s1, model, cache)
	subset2 = [i for i in subset if ! results[i]]
	results2 = similar(results)
	select!(results2, subset2, s.s2, model, cache)
	for i in subset2
		if results2[i]
			results[i] = true
		end
	end
end

struct XorSelector{T1<:Selector,T2<:Selector} <: Selector
	s1::Selector
	s2::Selector

	XorSelector{T1,T2}(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			new(s1, s2)
end

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::XorSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s1, model, cache)
	results2 = similar(results)
	select!(results2, subset2, s.s2, model, cache)
	for i in subset
		results[i] = results[i] ⊻ results2[i]
	end
end

struct CachedSelector{T<:Selector} <: Selector
	s::T
	groupkey::Symbol
	itemkey::Symbol

	CachedSelector{T}(s::T, groupkey::Symbol, itemkey::Symbol) where
			{T<:Selector} = new(s, groupkey, itemkey)
end

CachedSelector(s::T, groupkey::Symbol, itemkey::Symbol) where {T<:Selector} =
		CachedSelector{T}(s, groupkey, itemkey)

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::CachedSelector, model::ParticleCollection, cache::SelectionCache)
	cached::BitVector = get!(cache, s.groupkey, s.itemkey) do
		work = BitVector(undef, length(model))
		select!(work, eachindex(model), s.s, model, cache)
		work
	end
	for i in subset
		results[i] = cached[i]
	end
end

struct MapSelector{T<:AbstractVector{Bool}} <: Selector
	map::T

	MapSelector{T}(map::T) where {T<:AbstractVector{Bool}} = new(map)
end

MapSelector(map::T) where {T<:AbstractVector{Bool}} = MapSelector{T}(map)

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::MapSelector, model::ParticleCollection, cache::SelectionCache)
	for i in subset
		results[i] = s.map[i]
	end
end

struct IndexSelector{T1,T2<:SelectionMode,T3<:SelectionMode} <:Selector
	f::T1
	by::T2
	ineach::T3

	function IndexSelector{T1,T2,T3}(f::T1, by::T2, ineach::T3) where
			{T1,T2<:SelectionMode,T3<:SelectionMode}
		if ! iscompatible(IndexSelector, by, ineach)
			error("cannot select indices by $(by) in each $(ineach)")
		end
		new(f, by, ineach)
	end
end

IndexSelector(f::T1, by::T2, ineach::T3) where
		{T1,T2<:SelectionMode,T3<:SelectionMode} =
		IndexSelector{T1,T2,T3}(f, by, ineach)

iscompatible(::Type{<:IndexSelector},
		by::SelectionMode, ineach::SelectionMode) = false

iscompatible(::Type{<:IndexSelector},
		by::ParticleMode, ineach::ModelMode) = true

iscompatible(::Type{<:IndexSelector},
		by::SelectionMode, ineach::ModelMode) = true

iscompatible(::Type{<:IndexSelector},
		by::ParticleMode, ineach::SelectionMode) = true

iscompatible(::Type{<:IndexSelector}, by::ResidueMode, ineach::ChainMode) = true

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::IndexSelector, model::ParticleCollection, cache::SelectionCache)
	work = falses(length(model))
	selectindices!(work, s.f, s.by, s.ineach, model, cache)
	for i in subset
		results[i] = work[i]
	end
end

selectindices!(results::BitVector, f, by::ParticleMode, ineach::ModelMode,
		model::ParticleCollection, cache::SelectionCache) =
		selectindices0r!(results, f(eachindex(model)))

function selectindices0r!(results::BitVector,
		I::Union{AbstractVector{<:Integer},Integer})
	for i in I
		results[i] = true
	end
end

selectindices!(results::BitVector, f, by::ParticleMode, ineach::SelectionMode,
		model::ParticleCollection, cache::SelectionCache) =
		selectindices01!(results, f, geth2(model, ineach, cache).tree)

function selectindices01!(results::BitVector, f,
		tree::AbstractVector{<:AbstractVector{<:Integer}})
	for branch in tree
		selectindices01!(results, branch, f(eachindex(branch)))
	end
end

function selectindices01!(results::BitVector, branch::AbstractVector{<:Integer},
		I::Union{AbstractVector{<:Integer},Integer})
	for i in I
		results[branch[i]] = true
	end
end

selectindices!(results::BitVector, f, by::SelectionMode, ineach::ModelMode,
		model::ParticleCollection, cache::SelectionCache) =
		selectindices1r!(results, f, geth2(model, by, cache).tree)

selectindices1r!(results::BitVector, f,
		tree::AbstractVector{<:AbstractVector{<:Integer}}) =
		selectindices1r!(results, tree, f(eachindex(tree)))

function selectindices1r!(results::BitVector,
		tree::AbstractVector{<:AbstractVector{<:Integer}},
		I::Union{AbstractVector{<:Integer},Integer})
	for i in I
		for j in tree[i]
			results[j] = true
		end
	end
end

selectindices!(results::BitVector, f, by::ResidueMode, ineach::ChainMode,
		model::ParticleCollection, cache::SelectionCache) =
		selectindices12!(results, f, getmcrp(model, cache).tree)

function selectindices12!(results::BitVector, f,
		tree::AbstractVector{<:AbstractVector{<:AbstractVector{<:Integer}}})
	for branch in tree
		selectindices12!(results, branch, f(eachindex(branch)))
	end
end

function selectindices12!(results::BitVector,
		branch::AbstractVector{<:AbstractVector{<:Integer}},
		I::Union{AbstractVector{<:Integer},Integer})
	for i in I
		for j in branch[i]
			results[j] = true
		end
	end
end

abstract type ByGroupSelector <: Selector end

struct ExpandSelector{T1<:Selector,T2<:SelectionMode} <: ByGroupSelector
	s::T1
	by::T2

	function ExpandSelector{T1,T2}(s::T1, by::T2) where
			{T1<:Selector,T2<:SelectionMode}
		if ! iscompatible(ByGroupSelector, by)
			error("cannot expand selections by $(by)")
		end
		new(s, by)
	end
end

ExpandSelector(s::T1, by::T2) where {T1<:Selector,T2<:SelectionMode} =
		ExpandSelector{T1,T2}(s, by)

iscompatible(::Type{<:ByGroupSelector}, ::SelectionMode) = false

iscompatible(::Type{<:ByGroupSelector}, ::ChainMode) = true

iscompatible(::Type{<:ByGroupSelector}, ::ResidueMode) = true

iscompatible(::Type{<:ByGroupSelector}, ::FragmentMode) = true

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::ExpandSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s, model, cache)
	work = falses(length(model))
	selectpartial!(work, subset, results, geth2(model, s.by, cache)...)
	for i in subset
		results[i] = work[i]
	end
end

function selectpartial!(results::BitVector, subset::AbstractVector{<:Integer},
		selected::BitVector, tree::AbstractVector{<:AbstractVector{<:Integer}},
		paths::AbstractVector{<:Integer})
	for i in subset
		if selected[i]
			if ! results[i]
				for j in tree[paths[i]]
					results[j] = true
				end
			end
		end
	end
end

struct RestrictSelector{T1<:Selector,T2<:SelectionMode} <: ByGroupSelector
	s::T1
	by::T2

	function RestrictSelector{T1,T2}(s::T1, by::T2) where
			{T1<:Selector,T2<:SelectionMode}
			if ! iscompatible(ByGroupSelector, by)
				error("cannot restrict selections by $(by)")
			end
		new(s, by)
	end
end

RestrictSelector(s::T1, by::T2) where {T1<:Selector,T2<:SelectionMode} =
		RestrictSelector{T1,T2}(s, by)

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::RestrictSelector, model::ParticleCollection, cache::SelectionCache)
	select!(results, subset, s.s, model, cache)
	work = trues(length(model))
	selectfull!(work, subset, results, geth2(model, s.by, cache)...)
	for i in subset
		results[i] = work[i]
	end
end

function selectfull!(results::BitVector, subset::AbstractVector{<:Integer},
		selected::BitVector, tree::AbstractVector{<:AbstractVector{<:Integer}},
		paths::AbstractVector{<:Integer})
	for i in subset
		if ! selected[i]
			if results[i]
				for j in tree[paths[i]]
					results[j] = false
				end
			end
		end
	end
end

abstract type WithinSelector <: Selector end

struct WithinPositionSelector{T1,T2<:SelectionMode} <: WithinSelector
	dmax::Float64
	of::T1
	by::T2

	function WithinPositionSelector{T1,T2}(dmax::Real, of::T1, by::T2) where
			{T1,T2<:SelectionMode}
		if ! iscompatible(WithinSelector, by)
			error("cannot select within a distance by $(by)")
		end
		dmax > 0.0 || error("expected strictly positive distance")
		new(dmax, of, by)
	end
end

iscompatible(::Type{<:WithinSelector}, by::SelectionMode) =
		(by == ParticleMode() || iscompatible(ByGroupSelector, by))

WithinPositionSelector(dmax::Real, of::T1, by::T2) where
		{T1,T2<:SelectionMode} = WithinPositionSelector{T1,T2}(dmax, of, by)

WithinSelector(dmax::Real, of, by::SelectionMode) =
		WithinPositionSelector(dmax, positionfunction(of), by)

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::WithinPositionSelector, model::ParticleCollection,
		cache::SelectionCache)
	work = trues(length(model))
	for i in subset
		work[i] = false
	end
	cell = getpbccell(model, cache)
	Rw, Kw = getpbcpos(model.R, cell, cache)
	nbrtable = getnbrtable(length(Rw), cell, cache)
	Rref = s.of(model, cache)
	selectwithinpos!(work, s.dmax, s.by, Rw, Kw, Rref, cell, nbrtable, model,
			cache)
	for i in subset
		results[i] = work[i]
	end
end

selectwithinpos!(results::BitVector, dmax::Real, by::ParticleMode,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rref::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable, model::ParticleCollection,
		cache::SelectionCache) =
		selectwithinpos0!(results, dmax, Rw, Kw, Rref, cell, nbrtable)

function selectwithinpos0!(results::BitVector, dmax::Real,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rref::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable)
	dmax2 = dmax^2
	J = Int[]
	for Rrefi in Rref
		Rrefiw, Krefiw = pbcpos(Rrefi, cell)
		for j in findnear!(J, Rw, Kw, Rrefiw, Krefiw, cell, dmax, nbrtable)
			if ! results[j]
				if sqmindist(Rw[j], Kw[j], Rrefiw, Krefiw, cell) <= dmax2
					results[j] = true
				end
			end
		end
	end
end

selectwithinpos!(results::BitVector, dmax::Real, by::SelectionMode,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rref::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable, model::ParticleCollection,
		cache::SelectionCache) = selectwithinpos1!(results, dmax, Rw, Kw, Rref,
		cell, nbrtable, geth2(model, by, cache)...)

function selectwithinpos1!(results::BitVector, dmax::Real,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Rref::AbstractVector{Vector3D}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable,
		tree::AbstractVector{<:AbstractVector{<:Integer}},
		paths::AbstractVector{<:Integer})
	dmax2 = dmax^2
	J = Int[]
	for Rrefi in Rref
		Rrefiw, Krefiw = pbcpos(Rrefi, cell)
		for j in findnear!(J, Rw, Kw, Rrefiw, Krefiw, cell, dmax, nbrtable)
			if ! results[j]
				if sqmindist(Rw[j], Kw[j], Rrefiw, Krefiw, cell) <= dmax2
					for k in tree[paths[j]]
						results[k] = true
					end
				end
			end
		end
	end
end

struct WithinSelectionSelector{T1<:Selector,T2<:SelectionMode} <: WithinSelector
	dmax::Float64
	of::T1
	by::T2

	function WithinSelectionSelector{T1,T2}(dmax::Real, of::T1, by::T2) where
			{T1<:Selector,T2<:SelectionMode}
		if ! iscompatible(WithinSelector, by)
			error("cannot select within a distance by $(by)")
		end
		dmax > 0.0 || error("expected strictly positive distance")
		new(dmax, of, by)
	end
end

WithinSelectionSelector(dmax::Real, of::T1, by::T2) where
		{T1<:Selector,T2<:SelectionMode} =
		WithinSelectionSelector{T1,T2}(dmax, of, by)

WithinSelector(dmax::Real, of::Selector, by::SelectionMode) =
		WithinSelectionSelector(dmax, of, by)

function select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::WithinSelectionSelector, model::ParticleCollection,
		cache::SelectionCache)
	work = trues(length(model))
	for i in subset
		work[i] = false
	end
	cell = getpbccell(model, cache)
	Rw, Kw = getpbcpos(model.R, cell, cache)
	nbrtable = getnbrtable(length(Rw), cell, cache)
	Iref = findall(map(s.of, model, cache))
	selectwithinsel!(work, s.dmax, s.by, Rw, Kw, Iref, cell, nbrtable, model,
			cache)
	for i in subset
		results[i] = work[i]
	end
end

selectwithinsel!(results::BitVector, dmax::Real, by::ParticleMode,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Iref::AbstractVector{<:Integer}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable, model::ParticleCollection,
		cache::SelectionCache) =
		selectwithinsel0!(results, dmax, Rw, Kw, Iref, cell, nbrtable)

function selectwithinsel0!(results::BitVector, dmax::Real,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Iref::AbstractVector{<:Integer}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable)
	dmax2 = dmax^2
	J = Int[]
	for i in Iref
		for j in findnear!(J, Rw, Kw, i, cell, dmax, nbrtable)
			if ! results[j]
				if sqmindist(Rw[j], Kw[j], Rw[i], Kw[i], cell) <= dmax2
					results[j] = true
				end
			end
		end
	end
end

selectwithinsel!(results::BitVector, dmax::Real, by::SelectionMode,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Iref::AbstractVector{<:Integer}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable, model::ParticleCollection,
		cache::SelectionCache) = selectwithinsel1!(results, dmax, Rw, Kw, Iref,
		cell, nbrtable, geth2(model, by, cache)...)

function selectwithinsel1!(results::BitVector, dmax::Real,
		Rw::AbstractVector{Vector3D}, Kw::AbstractVector{Vector3D},
		Iref::AbstractVector{<:Integer}, cell::Union{TriclinicPBC,Nothing},
		nbrtable::NeighborTable,
		tree::AbstractVector{<:AbstractVector{<:Integer}},
		paths::AbstractVector{<:Integer})
	dmax2 = dmax^2
	J = Int[]
	for i in Iref
		for j in findnear!(J, Rw, Kw, i, cell, dmax, nbrtable)
			if ! results[j]
				if sqmindist(Rw[j], Kw[j], Rw[i], Kw[i], cell) <= dmax2
					for k in tree[paths[j]]
						results[k] = true
					end
				end
			end
		end
	end
end

struct PropertySelector{T1,T2} <: Selector
	getter::T1
	predicate::T2

	PropertySelector{T1,T2}(getter::T1, predicate::T2) where {T1,T2} =
			new(getter, predicate)
end

PropertySelector(getter::T1, predicate::T2) where {T1,T2} =
		PropertySelector{T1,T2}(getter, predicate)

select!(results::BitVector, subset::AbstractVector{<:Integer},
		s::PropertySelector, model::ParticleCollection, cache::SelectionCache) =
		selecttrue!(results, subset, s.predicate, s.getter(model, cache))

function selecttrue!(results::BitVector, subset::AbstractVector{<:Integer},
		predicate, V::AbstractVector)
	for i in subset
		results[i] = predicate(V[i])
	end
end

module Selectors

	using ..Dorothy
	using ..Dorothy.Utils
	using ..Dorothy.Properties

	export
			Chain, Residue, Fragment,

			Map, Index, Expand, Restrict, Within,

			Id, Name, ResId, ResName, ChainId, Element, R, V, F, Mass, Charge,
			BFactor, Occupancy, Charge, SS,

			VSite, Hydrogen, Heavy, Water, Protein, AcidResidue, BasicResidue,
			ChargedResidue, PolarResidue, HydrophobicResidue, Backbone,
			MainChain, SideChain, Calpha, Cα, Nter, Cter, NuclAcid, Lipid, Ion,
			Helix, AlphaHelix, Helix310, PiHelix, Turn, Sheet, Strand, Bridge,
			Loop, Coil, Bend

	const Model = Dorothy.ModelMode()
	const Chain = Dorothy.ChainMode()
	const Residue = Dorothy.ResidueMode()
	const Particle = Dorothy.ParticleMode()
	const Fragment = Dorothy.FragmentMode()

	Base.:(!)(s::T) where {T<:Selector} = Dorothy.NotSelector{T}(s)

	Base.:(&)(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			Dorothy.AndSelector{T1,T2}(s1, s2)

	Base.:(|)(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			Dorothy.OrSelector{T1,T2}(s1, s2)

	Base.:(⊻)(s1::T1, s2::T2) where {T1<:Selector,T2<:Selector} =
			Dorothy.XorSelector{T1,T2}(s1, s2)

	Map(map::AbstractVector{Bool}) = Dorothy.MapSelector(map)

	function Map(view::MolecularModelView)
		model = parent(view)
		map = falses(length(model))
		for i in parentindices(view)
			map[i] = true
		end
		Dorothy.MapSelector(map)
	end

	Index(X; by::SelectionMode = Particle, ineach::SelectionMode = Model) =
			Dorothy.IndexSelector(Dorothy.indexfunction(X), by, ineach)

	Expand(s::Selector; by::SelectionMode) = Dorothy.ExpandSelector(s, by)

	Restrict(s::Selector; by::SelectionMode) = Dorothy.RestrictSelector(s, by)

	Within(dmax::Real; of) = Dorothy.WithinSelector(dmax, of, Particle)

	Expand(s::Dorothy.WithinPositionSelector; by::SelectionMode) =
			Dorothy.WithinPositionSelector(s.dmax, s.of, by)

	Expand(s::Dorothy.WithinSelectionSelector; by::SelectionMode) =
			Dorothy.WithinSelectionSelector(s.dmax, s.of, by)

	Id(X...) = Dorothy.PropertySelector((model, cache) -> model.ids,
			Dorothy.integerpredicate(X...))

	Name(X...) = Dorothy.PropertySelector((model, cache) -> model.names,
			Dorothy.stringpredicate(X...))

	ResId(X...) = Dorothy.PropertySelector((model, cache) -> model.resids,
			Dorothy.integerpredicate(X...))

	ResName(X...) = Dorothy.PropertySelector((model, cache) -> model.resnames,
			Dorothy.stringpredicate(X...))

	ChainId(X...) = Dorothy.PropertySelector(Dorothy.getchainids,
			Dorothy.stringpredicate(X...))

	Element(X...) = Dorothy.PropertySelector(Dorothy.getelements,
			Dorothy.stringpredicate(X...))

	R(x) = Dorothy.PropertySelector((model, cache) -> model.R,
			Dorothy.vector3dpredicate(x))

	V(x) = Dorothy.PropertySelector((model, cache) -> model.V,
			Dorothy.vector3dpredicate(x))

	F(x) = Dorothy.PropertySelector((model, cache) -> model.F,
			Dorothy.vector3dpredicate(x))

	Mass(x) = Dorothy.PropertySelector(Dorothy.getmasses,
			Dorothy.realpredicate(x))

	Charge(x) = Dorothy.PropertySelector((model, cache) -> model.charges,
			Dorothy.realpredicate(x))

	BFactor(x) = Dorothy.PropertySelector((model, cache) -> model.bfactors,
			Dorothy.realpredicate(x))

	Occupancy(x) = Dorothy.PropertySelector((model, cache) ->
			model.occupancies, Dorothy.realpredicate(x))

	SS(X...) = Dorothy.PropertySelector(Dorothy.getss,
			Dorothy.stringpredicate(X...))

	const VSite = Dorothy.CachedSelector(Name(isvsite), :system, :vsites)

	const Hydrogen = Dorothy.CachedSelector(Element("H"), :system, :hydrogens)

	const Heavy = ! (Hydrogen | VSite)

	const Water = Dorothy.CachedSelector(ResName(iswater), :system, :water)

	const Protein =
			Dorothy.CachedSelector(ResName(isprotein), :system, :protein)

	const AcidResidue = Dorothy.CachedSelector(ResName(isacidresidue),
			:system, :acidresidues)

	const BasicResidue = Dorothy.CachedSelector(ResName(isbasicresidue),
			:system, :basicresidues)

	const ChargedResidue = Dorothy.CachedSelector(ResName(ischargedresidue),
			:system, :chargedresidues)

	const PolarResidue = Dorothy.CachedSelector(ResName(ispolarresidue),
			:system, :polarresidues)

	const HydrophobicResidue =
			Dorothy.CachedSelector(ResName(ishydrophobicresidue),
			:system, :hydrophobicresidues)

	const MainChainName =
			Dorothy.CachedSelector(Name(Properties.ismainchainname),
			:system, :mainchainnames)

	const MainChain = MainChainName & Protein

	const SideChain = !MainChainName & Protein

	const BackboneName =
			Dorothy.CachedSelector(Name(Properties.isbackbonename),
			:system, :backbonenames)

	const Backbone = BackboneName & Protein

	const CAlpha = Name("CA") & Protein

	const Cα = CAlpha

	const Nter = Index(1, by=Residue, ineach=Chain) & Protein

	const Cter = Index(last, by=Residue, ineach=Chain) & Protein

	const NuclAcid = Dorothy.CachedSelector(ResName(isnuclacid),
			:system, :nuclacid)

	const Lipid = Dorothy.CachedSelector(ResName(islipid), :system, :lipid)

	const Ion = Dorothy.CachedSelector(ResName(ision), :system, :ion)

	const Helix = Dorothy.CachedSelector(SS(ishelix), :frame, :helix)

	const AlphaHelix = SS(isalphahelix)

	const Helix310 = SS(ishelix310)

	const PiHelix = SS(ispihelix)

	const Turn = SS(isturn)

	const Sheet = Dorothy.CachedSelector(SS(issheet), :frame, :sheet)

	const Strand = SS(isstrand)

	const Bridge = SS(isbridge)

	const Loop = Dorothy.CachedSelector(SS(isloop), :frame, :loop)

	const Coil = SS(iscoil)

	const Bend = SS(isbend)

end


text:
#=
create_linked_list:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-30
=#

module LinkedListCreation

    # node class
    mutable struct Node
        data
        p_next
    end

    # head node
    head = nothing

    # method to print contents of linked list
    function print_list()
        tmp = head
        while tmp != nothing
            println(tmp.data)
            tmp = tmp.p_next
        end
    end

    function main()
        # create simple linked list with 3 nodes
        # next pointer of each node is initialized as Null
        global head = Node(1, nothing)
        second = Node(2, nothing)
        third = Node(3, nothing)

        # link first node with second
        head.p_next = second

        # link second node with third
        second.p_next = third

        print_list()
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    using .LinkedListCreation
    LinkedListCreation.main()
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
export mouseover, mouse_selection, mouseposition, hovered_scene
export select_rectangle

mouseover() = error("not implemented")

# What does this function do?
function mouse_in_scene(scene)
    p = rootparent(scene)
    lift(pixelarea(p), pixelarea(scene), events(scene).mouseposition) do pa, sa, mp
        Vec(mp) .- minimum(sa)
    end
end

"""
    mouse_selection(scene)

Returns the plot under the current mouse position in `scene`.
"""
function mouse_selection
    # TODO this needs to be implemented here via select_mouse_native
end

# What does this function do?
to_screen(scene, mpos) = Point2f0(mpos) .- Point2f0(minimum(pixelarea(scene)[]))

"""
    mouseposition(scene = hovered_scene()) -> pos
Return the current position of the mouse `pos` in _data points_ of the
given `scene`.

By default uses the `scene` that the mouse is currently hovering over.
"""
mouseposition(scene = hovered_scene()) =
to_world(scene, to_screen(scene, events(scene).mouseposition[]))

"""
    hovered_scene()
Return the `scene` that the mouse is currently hovering over.

Properly identifies the scene for a plot with multiple sub-plots.
"""
hovered_scene() = error("hoevered_scene is not implemented yet.")

"""
    select_rectangle(scene; kwargs...) -> rect
Interactively select a rectangle on a `scene` by clicking the left mouse button,
dragging and then un-clicking. The function returns an **observable** `rect` whose
value corresponds to the selected rectangle on the scene. In addition the function
_plots_ the selected rectangle on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted rectangle disappears.

The value of the returned observable is updated **only** when the user un-clicks
(i.e. when the final value of the rectangle has been decided) and only if the
rectangle has area > 0.

The `kwargs...` are propagated into `lines!` which plots the selected rectangle.
"""
function select_rectangle(scene; kwargs...)
    key = Mouse.left
    waspressed = Node(false)
    rect = Node(FRect(0, 0, 1, 1)) # plotted rectangle
    rect_ret = Node(FRect(0, 0, 1, 1)) # returned rectangle

    # Create an initially hidden rectangle
    rect_vis = lines!(
        scene, rect, raw = true, visible = false, kwargs...,
    )[end] # Why do I have to do [end] ?

    on(events(scene).mousedrag) do drag
        if ispressed(scene, key) && is_mouseinside(scene)
            mp = mouseposition(scene)
            if drag == Mouse.down
                waspressed[] = true
                rect_vis[:visible] = true # start displaying
                rect[] = FRect(mp, 0.0, 0.0)
            elseif drag == Mouse.pressed
                mini = minimum(rect[])
                rect[] = FRect(mini, mp - mini)
            end
        else
            if drag == Mouse.up && waspressed[] # User has selected the rectangle
                waspressed[] = false
                r = absrect(rect[])
                w, h = widths(r)
                if w > 0.0 && h > 0.0 # Ensure that the rectangle has non0 size.
                    rect_ret[] = r
                end
            end
            # always hide if not the right key is pressed
            rect_vis[:visible] = false # make the plotted rectangle invisible
        end
        return rect_ret
    end
    return rect_ret
end


text:
module TestJoin
    using Base.Test
    using DataTables

    name = DataTable(ID = NullableArray([1, 2, 3]),
                     Name = NullableArray(["John Doe", "Jane Doe", "Joe Blogs"]))
    job = DataTable(ID = NullableArray([1, 2, 2, 4]),
                    Job = NullableArray(["Lawyer", "Doctor", "Florist", "Farmer"]))

    # Join on symbols or vectors of symbols
    join(name, job, on = :ID)
    join(name, job, on = [:ID])

    # Soon we won't allow natural joins
    #@test_throws join(name, job)

    # Test output of various join types
    outer = DataTable(ID = NullableArray([1, 2, 2, 3, 4]),
                      Name = NullableArray(["John Doe", "Jane Doe", "Jane Doe", "Joe Blogs", Nullable()]),
                      Job = NullableArray(["Lawyer", "Doctor", "Florist", Nullable(), "Farmer"]))

    # (Tests use current column ordering but don't promote it)
    right = outer[Bool[!isnull(x) for x in outer[:Job]], [:ID, :Name, :Job]]
    left = outer[Bool[!isnull(x) for x in outer[:Name]], :]
    inner = left[Bool[!isnull(x) for x in left[:Job]], :]
    semi = unique(inner[:, [:ID, :Name]])
    anti = left[Bool[isnull(x) for x in left[:Job]], [:ID, :Name]]

    @test isequal(join(name, job, on = :ID), inner)
    @test isequal(join(name, job, on = :ID, kind = :inner), inner)
    @test isequal(join(name, job, on = :ID, kind = :outer), outer)
    @test isequal(join(name, job, on = :ID, kind = :left), left)
    @test isequal(join(name, job, on = :ID, kind = :right), right)
    @test isequal(join(name, job, on = :ID, kind = :semi), semi)
    @test isequal(join(name, job, on = :ID, kind = :anti), anti)
    @test_throws ArgumentError join(name, job)
    @test_throws ArgumentError join(name, job, on=:ID, kind=:other)

    # Join with no non-key columns
    on = [:ID]
    nameid = name[on]
    jobid = job[on]

    @test isequal(join(nameid, jobid, on = :ID), inner[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :inner), inner[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :outer), outer[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :left), left[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :right), right[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :semi), semi[on])
    @test isequal(join(nameid, jobid, on = :ID, kind = :anti), anti[on])

    # Join on multiple keys
    dt1 = DataTable(A = 1, B = 2, C = 3)
    dt2 = DataTable(A = 1, B = 2, D = 4)

    join(dt1, dt2, on = [:A, :B])

    # Test output of cross joins
    dt1 = DataTable(A = 1:2, B = 'a':'b')
    dt2 = DataTable(A = 1:3, C = 3:5)

    cross = DataTable(A = [1, 1, 1, 2, 2, 2],
                      B = ['a', 'a', 'a', 'b', 'b', 'b'],
                      C = [3, 4, 5, 3, 4, 5])

    @test isequal(join(dt1, dt2[[:C]], kind = :cross), cross)

    # Cross joins handle naming collisions
    @test size(join(dt1, dt1, kind = :cross)) == (4, 4)

    # Cross joins don't take keys
    @test_throws ArgumentError join(dt1, dt2, on = :A, kind = :cross)

    # test empty inputs
    simple_dt(len::Int, col=:A) = (dt = DataTable(); dt[col]=NullableArray(collect(1:len)); dt)
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :left),  simple_dt(0))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :left),  simple_dt(2))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :left),  simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :right), simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :right), simple_dt(2))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :right), simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :inner), simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :inner), simple_dt(0))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :inner), simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :outer), simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :outer), simple_dt(2))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :outer), simple_dt(2))
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :semi),  simple_dt(0))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :semi),  simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :semi),  simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(0), on = :A, kind = :anti),  simple_dt(0))
    @test isequal(join(simple_dt(2), simple_dt(0), on = :A, kind = :anti),  simple_dt(2))
    @test isequal(join(simple_dt(0), simple_dt(2), on = :A, kind = :anti),  simple_dt(0))
    @test isequal(join(simple_dt(0), simple_dt(0, :B), kind = :cross), DataTable(A=Int[], B=Int[]))
    @test isequal(join(simple_dt(0), simple_dt(2, :B), kind = :cross), DataTable(A=Int[], B=Int[]))
    @test isequal(join(simple_dt(2), simple_dt(0, :B), kind = :cross), DataTable(A=Int[], B=Int[]))

    # issue #960
    dt1 = DataTable(A = 1:50,
                    B = 1:50,
                    C = 1)
    categorical!(dt1, :A)
    categorical!(dt1, :B)
    join(dt1, dt1, on = [:A, :B], kind = :inner)

    # Test that Array{Nullable} works when combined with NullableArray (#1088)
    dt = DataTable(Name = Nullable{String}["A", "B", "C"],
                   Mass = [1.5, 2.2, 1.1])
    dt2 = DataTable(Name = ["A", "B", "C", "A"],
                    Quantity = [3, 3, 2, 4])
    @test join(dt2, dt, on=:Name, kind=:left) == DataTable(Name = NullableArray(["A", "B", "C", "A"]),
                                                           Quantity = NullableArray([3, 3, 2, 4]),
                                                           Mass = NullableArray([1.5, 2.2, 1.1, 1.5]))

    # Test that join works when mixing Array and NullableArray (#1151)
    dt = DataTable([collect(1:10), collect(2:11)], [:x, :y])
    dtnull = DataTable(x = NullableArray(1:10), z = NullableArray(3:12))
    @test join(dt, dtnull, on = :x) ==
        DataTable([collect(1:10), collect(2:11), NullableArray(3:12)], [:x, :y, :z])
    @test join(dtnull, dt, on = :x) ==
        DataTable([NullableArray(1:10), NullableArray(3:12), collect(2:11)], [:x, :z, :y])

    @testset "all joins" begin
        dt1 = DataTable(Any[[1, 3, 5], [1.0, 3.0, 5.0]], [:id, :fid])
        dt2 = DataTable(Any[[0, 1, 2, 3, 4], [0.0, 1.0, 2.0, 3.0, 4.0]], [:id, :fid])
        N = Nullable()

        @test join(dt1, dt2, kind=:cross) ==
            DataTable(Any[repeat([1, 3, 5], inner = 5),
                          repeat([1, 3, 5], inner = 5),
                          repeat([0, 1, 2, 3, 4], outer = 3),
                          repeat([0, 1, 2, 3, 4], outer = 3)],
                      [:id, :fid, :id_1, :fid_1])
        @test typeof.(join(dt1, dt2, kind=:cross).columns) ==
            [Vector{Int}, Vector{Float64}, Vector{Int}, Vector{Float64}]

        i(on) = join(dt1, dt2, on = on, kind = :inner)
        l(on) = join(dt1, dt2, on = on, kind = :left)
        r(on) = join(dt1, dt2, on = on, kind = :right)
        o(on) = join(dt1, dt2, on = on, kind = :outer)
        s(on) = join(dt1, dt2, on = on, kind = :semi)
        a(on) = join(dt1, dt2, on = on, kind = :anti)

        @test s(:id) ==
              s(:fid) ==
              s([:id, :fid]) == DataTable(Any[[1, 3], [1, 3]], [:id, :fid])
        @test typeof.(s(:id).columns) ==
              typeof.(s(:fid).columns) ==
              typeof.(s([:id, :fid]).columns) == [Vector{Int}, Vector{Float64}]
        @test a(:id) ==
              a(:fid) ==
              a([:id, :fid]) == DataTable(Any[[5], [5]], [:id, :fid])
        @test typeof.(a(:id).columns) ==
              typeof.(a(:fid).columns) ==
              typeof.(a([:id, :fid]).columns) == [Vector{Int}, Vector{Float64}]

        on = :id
        @test i(on) == DataTable(Any[[1, 3], [1, 3], [1, 3]], [:id, :fid, :fid_1])
        @test typeof.(i(on).columns) == [Vector{Int}, Vector{Float64}, Vector{Float64}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]),
                                 fid_1 = NullableArray([1, 3, N]))
        @test typeof.(l(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Float64}]
        @test r(on) == DataTable(id = NullableArray([1, 3, 0, 2, 4]),
                                 fid = NullableArray([1, 3, N, N, N]),
                                 fid_1 = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Float64}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 5, N, N, N]),
                                 fid_1 = NullableArray([1, 3, N, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Float64}]

        on = :fid
        @test i(on) == DataTable(Any[[1, 3], [1.0, 3.0], [1, 3]], [:id, :fid, :id_1])
        @test typeof.(i(on).columns) == [Vector{Int}, Vector{Float64}, Vector{Int}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]),
                                 id_1 = NullableArray([1, 3, N]))
        @test typeof.(l(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Int}]
        @test r(on) == DataTable(id = NullableArray([1, 3, N, N, N]),
                                 fid = NullableArray([1, 3, 0, 2, 4]),
                                 id_1 = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Int}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, N, N, N]),
                                 fid = NullableArray([1, 3, 5, 0, 2, 4]),
                                 id_1 = NullableArray([1, 3, N, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64},
                                         NullableVector{Int}]

        on = [:id, :fid]
        @test i(on) == DataTable(Any[[1, 3], [1, 3]], [:id, :fid])
        @test typeof.(i(on).columns) == [Vector{Int}, Vector{Float64}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]))
        @test typeof.(l(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64}]
        @test r(on) == DataTable(id = NullableArray([1, 3, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 5, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableVector{Int},
                                         NullableVector{Float64}]
    end

    @testset "all joins with CategoricalArrays" begin
        dt1 = DataTable(Any[CategoricalArray([1, 3, 5]),
                              CategoricalArray([1.0, 3.0, 5.0])], [:id, :fid])
        dt2 = DataTable(Any[CategoricalArray([0, 1, 2, 3, 4]),
                              CategoricalArray([0.0, 1.0, 2.0, 3.0, 4.0])], [:id, :fid])
        N = Nullable()
        DRT = CategoricalArrays.DefaultRefType

        @test join(dt1, dt2, kind=:cross) ==
            DataTable(Any[repeat([1, 3, 5], inner = 5),
                          repeat([1, 3, 5], inner = 5),
                          repeat([0, 1, 2, 3, 4], outer = 3),
                          repeat([0, 1, 2, 3, 4], outer = 3)],
                      [:id, :fid, :id_1, :fid_1])
        @test typeof.(join(dt1, dt2, kind=:cross).columns) ==
            [CategoricalVector{i, DRT} for i in [Int, Float64, Int, Float64]]

        i(on) = join(dt1, dt2, on = on, kind = :inner)
        l(on) = join(dt1, dt2, on = on, kind = :left)
        r(on) = join(dt1, dt2, on = on, kind = :right)
        o(on) = join(dt1, dt2, on = on, kind = :outer)
        s(on) = join(dt1, dt2, on = on, kind = :semi)
        a(on) = join(dt1, dt2, on = on, kind = :anti)

        @test s(:id) ==
              s(:fid) ==
              s([:id, :fid]) == DataTable(Any[[1, 3], [1, 3]], [:id, :fid])
        @test typeof.(s(:id).columns) ==
              typeof.(s(:fid).columns) ==
              typeof.(s([:id, :fid]).columns) == [CategoricalVector{Int, DRT},
                                                  CategoricalVector{Float64, DRT}]
        @test a(:id) ==
              a(:fid) ==
              a([:id, :fid]) == DataTable(Any[[5], [5]], [:id, :fid])
        @test typeof.(a(:id).columns) ==
              typeof.(a(:fid).columns) ==
              typeof.(a([:id, :fid]).columns) == [CategoricalVector{Int, DRT},
                                                  CategoricalVector{Float64, DRT}]

        on = :id
        @test i(on) == DataTable(Any[[1, 3], [1, 3], [1, 3]], [:id, :fid, :fid_1])
        @test typeof.(i(on).columns) == [CategoricalVector{Int, DRT},
                                         CategoricalVector{Float64, DRT},
                                         CategoricalVector{Float64, DRT}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]),
                                 fid_1 = NullableArray([1, 3, N]))
        @test typeof.(l(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Float64, DRT}]
        @test r(on) == DataTable(id = NullableArray([1, 3, 0, 2, 4]),
                                 fid = NullableArray([1, 3, N, N, N]),
                                 fid_1 = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Float64, DRT}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 5, N, N, N]),
                                 fid_1 = NullableArray([1, 3, N, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Float64, DRT}]

        on = :fid
        @test i(on) == DataTable(Any[[1, 3], [1.0, 3.0], [1, 3]], [:id, :fid, :id_1])
        @test typeof.(i(on).columns) == [CategoricalVector{Int, DRT},
                                         CategoricalVector{Float64, DRT},
                                         CategoricalVector{Int, DRT}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]),
                                 id_1 = NullableArray([1, 3, N]))
        @test typeof.(l(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Int, DRT}]
        @test r(on) == DataTable(id = NullableArray([1, 3, N, N, N]),
                                 fid = NullableArray([1, 3, 0, 2, 4]),
                                 id_1 = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Int, DRT}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, N, N, N]),
                                 fid = NullableArray([1, 3, 5, 0, 2, 4]),
                                 id_1 = NullableArray([1, 3, N, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT},
                                         NullableCategoricalVector{Int, DRT}]

        on = [:id, :fid]
        @test i(on) == DataTable(Any[[1, 3], [1, 3]], [:id, :fid])
        @test typeof.(i(on).columns) == [CategoricalVector{Int, DRT},
                                         CategoricalVector{Float64, DRT}]
        @test l(on) == DataTable(id = NullableArray([1, 3, 5]),
                                 fid = NullableArray([1, 3, 5]))
        @test typeof.(l(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT}]
        @test r(on) == DataTable(id = NullableArray([1, 3, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 0, 2, 4]))
        @test typeof.(r(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT}]
        @test o(on) == DataTable(id = NullableArray([1, 3, 5, 0, 2, 4]),
                                 fid = NullableArray([1, 3, 5, 0, 2, 4]))
        @test typeof.(o(on).columns) == [NullableCategoricalVector{Int, DRT},
                                         NullableCategoricalVector{Float64, DRT}]
    end
end


text:
# test CommonConvWrapper

using Test
using NLsolve
using IncrementalInference
using Statistics

import IncrementalInference: getSample

mutable struct FunctorWorks
  a::Array{Float64,2}
end

@testset "FunctorWorks" begin

function (fw::FunctorWorks)(x)
  fw.a[1,1] = -1.0
  nothing
end

global A = rand(2,3)
global At = deepcopy(A)
global At[1,1] = -1.0

global fvar = FunctorWorks(A)
fvar(0.0)
@test At == A

end

mutable struct FunctorArray{T}
  # This was a tuple, but array will like work better in the long term
  fnc!::Function
  a::Array{T, 1}
end

@testset "FunctorArray" begin

function testarray!(a1::Array{Float64, 2}, a2::Array{Float64,2})
  a1[1,1] = -1.0
  @show a1
  nothing
end

function (fi::FunctorArray)(x)
  fi.fnc!(fi.a...)
end

global A = rand(2,3)
global B = rand(2,3)
global t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
global At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorArray(testarray!, t)
fvar([0.0])
@test At == A

end

struct Test2 <: FunctorInferenceType
end

@testset "CommonConvWrapper test" begin

function (tt::Test2)(res::Vector{Float64}, userdata::FactorMetadata, idx::Int, meas::Tuple{Array{Float64,2}}, tp1::Array{Float64,2}, tp2::Array{Float64,2})
  tp1[1,1]=-2.0;
  res[:] .= 1.0
  nothing;
end

global N = 3
global tst2 = Test2()
global A = rand(2,N)
global B = rand(2,N)
global At = deepcopy(A)
global t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
t[1][1,1] = -10.0
@test A[1,1] == -10
# @show typeof(t)
# generalwrapper = GenericWrapParam{Test2}(tst2, t, 1, 1)
global ccw = CommonConvWrapper(tst2, t[1], 2, t) # generalwrapper.measurement = rand(1,1)

ccw.cpt[Threads.threadid()].factormetadata
ccw.cpt[Threads.threadid()].particleidx
ccw.measurement = (randn(1,N),)
ccw.params[ccw.cpt[Threads.threadid()].activehypo]

global x, res = zeros(2), zeros(2)
@time ccw(res, x)

At[1,1] = -2.0
At[2,1] = 0.0
@test A == At

# resource requirements

# fm = FactorMetadata()
# meas = (randn(2,3), )
#
# @time tst2(res, fm, 1, meas, A, B)

end

# abstract Nonparametric <: Function
# This is what the internmediate user would be contributing
mutable struct Pose1Pose1Test{T} <: FunctorPairwise
  Dx::T
  Pose1Pose1Test{T}() where T = new()
  Pose1Pose1Test{T}(a::T) where T = new(a)
  # Pose1Pose1Test(a::T) where T = new(a)
end
Pose1Pose1Test(a::T) where T = Pose1Pose1Test{T}(a)

getSample(pp1t::Pose1Pose1Test{T}, N::Int=1) where T = (reshape(rand(pp1t.Dx,N),1,N),)

#proposed standardized parameter list, does not have to be functor
function (Dp::Pose1Pose1Test)(res::Array{Float64},
                              userdata::FactorMetadata,
                              idx::Int,
                              meas::Tuple{Array{Float64,2}},
                              p1::Array{Float64,2},
                              p2::Array{Float64,2} )
  #
  res[1] = meas[1][1,idx] - (p2[1,idx] - p1[1,idx])
  nothing
end

@testset "Test in factor graph setting..." begin

global N = 100
global p1 = rand(1,N)
global p2 = rand(1,N)
global t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

global odo = Pose1Pose1Test{Normal}(Normal(100.0,1.0))
global ccw = CommonConvWrapper(odo, t[1], 1, t, measurement=getSample(odo, N))

freshSamples!(ccw, N)
global x, res = zeros(1), zeros(1)

@time for n in 1:N
  ccw.cpt[Threads.threadid()].particleidx = n

  ccw(res, x)
  # each point should be near 100.0
  @test 50.0 < res[1] < 150.0
end

# which variable position to solve for
ccw.varidx = 2

# do for a single point (i=1)

r = nlsolve( ccw, [0.0], inplace=true )

# for all points
@time for i in 1:N
  ccw.cpt[Threads.threadid()].particleidx = i
    # ccw(x, res)
  r = nlsolve( ccw, ccw.params[ccw.varidx][:,ccw.cpt[Threads.threadid()].particleidx], inplace=true )
  ccw.params[ccw.varidx][1,ccw.cpt[Threads.threadid()].particleidx] = r.zero[1]
end

@test abs(Statistics.mean(p1)-0.0) < 4.0
@test abs(Statistics.mean(p2)-100.0) < 4.0

end

@testset "Test with CommonConvWrapper for un-permuted root finding..." begin

global N = 110
global p1 = rand(1,N)
global p2 = rand(1,N)
global t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

global odo = Pose1Pose1Test(Normal(100.0,1.0))
# varidx=2 means we are solving for p2 relative to p1

global measurement = getSample(odo, N)
@show zDim = size(measurement[1],1)

global solvefor = 2

global ccw = CommonConvWrapper(odo, t[solvefor], zDim, t, measurement=measurement)
@show ccw.varidx = solvefor
# gwp = GenericWrapParam{Pose1Pose1Test}(odo, t, 2, 1, (zeros(0,1),) , getSample) #getSample(odo, N)

freshSamples!(ccw, N)
# and return complete fr/gwp
@time for n in 1:N
  # gwp(x, res)
  ccw.cpt[Threads.threadid()].particleidx = n
  numericRootGenericRandomizedFnc!( ccw )
end

# @show gwp.params

@test 90.0 < Statistics.mean(ccw.params[ccw.varidx]) < 110.0
@test -10.0 < Statistics.mean(ccw.params[1]) < 10.0

println("and in the reverse direction, achieved by simply changing CommonConvWrapper.varidx to 1...")

global solvefor = 1
@show ccw.varidx = solvefor
ccw.params[solvefor][:,:] = -100.0*ones(size(ccw.params[solvefor]))
ccw.cpt[Threads.threadid()].X = ccw.params[solvefor]

# @show gwp.params

# fr = FastRootGenericWrapParam{Pose1Pose1Test}(gwp.params[gwp.varidx], zDim, gwp)

freshSamples!(ccw, N)
@time for n in 1:N
  # gwp(x, res)
  ccw.cpt[Threads.threadid()].particleidx = n
  numericRootGenericRandomizedFnc!( ccw )
end

# @show gwp.params

@test -10.0 < Statistics.mean(ccw.params[1]) < 10.0
@test 90.0 < Statistics.mean(ccw.params[2]) < 110.0

end

@testset "Test with FastRootGenericWrapParam for permuted root finding..." begin

@warn "test not implemented yet"

end

# use the range only example, should give a circle with nothing in the middle

@testset "Generic convolution testing in factor graph context..." begin

global N=100
global p1 = randn(1,N)
global d1 = kde!(p1)
global p2 = randn(1,N)
global t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

global fg = initfg()
# fg.registeredModuleFunctions[:Main] = getSample

global v1=addVariable!(fg, :x1, ContinuousScalar, N=N)
global v2=addVariable!(fg, :x2, ContinuousScalar, N=N)
global bws = getBW(d1)[:,1]
global f1 = addFactor!(fg, [v1], Prior(kde!(p1, bws)) )

global odo = Pose1Pose1Test(Normal(100.0,1.0))
global f2 = addFactor!(fg, [v1;v2], odo)

global tree = wipeBuildNewTree!(fg)

global pts = getVal(getVariable(fg,:x1))
@test abs(Statistics.mean(pts)-0.0) < 10.0
global pts = getVal(getVariable(fg,:x2))
@test abs(Statistics.mean(pts)-0.0) < 10.0

# inferOverTreeR!(fg, tree, N=N)
# inferOverTreeR!(fg, tree)
# @time [inferOverTreeR!(fg, tree, N=N) for i in 1:3];
tree, smt, hist = solveTree!(fg)

# using Gadfly
# plot(y=rand(10))
# plotKDE(getVertKDE(fg,:x2))

global pts = getVal(getVariable(fg,:x1))
@test abs(Statistics.mean(pts)-0.0) < 10.0

global pts = getVal(getVariable(fg,:x2))
@test abs(Statistics.mean(pts)-100.0) < 10.0

end


text:
module A5P2

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\A5P2.jl"),
  informat="script",
  out_path = "$(pwd())\\A5P2.html",
  doctype = "md2html")
=#
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt,
  HCubature, StaticArrays, CTNumerical, Optim

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const DATE_FORMAT_STR = "m/d/y"
const SP500_NAME = "GSPC"
const TBILL_NAME = "t-bill"
const TERM_NAME = "termStructure"
const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const MAX_YEAR = 2016
const IN_SAMPLE_CUTOFF_YEAR = 1975

const DEFAULT_NUM_SIMS = 10000
const PlotContainer = Union{Plot,Gadfly.Compose.Context}
const CSRegSpec = Tuple{Symbol, Symbol, Int}
const CPRegSpec = Tuple{Symbol, Int}

const MIN_VALS_FOR_EST = 10
const Ι = Complex(0.0,1.0)

abstract type AbstractΘState end
abstract type AbstractΦFunction end

#convenience function for integration in the pure julia implementation
function I1D(f::Function, a::Float64, b::Float64)

  wrapFunc(x::SVector{1,Float64})::SVector{1,Float64} = SVector{1,Float64}(f(x[1]))
  ans::SVector{1,Float64}, err::Float64 =
    hcubature(wrapFunc, SVector{1,Float64}(a), SVector{1,Float64}(b))

  return ans[1]
end

function preProcessSP500()::Void
  dateFormat::String = DATE_FORMAT_STR
  sp500DF::DataFrame = readtable("$DATA_PATH\\$SP500_NAME.csv")

  #fix the dates
  rename!(sp500DF, :Date, :DATE)
  sp500DF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(sp500DF[:,:DATE])::DataVector{Date}
  sort!(sp500DF, cols = [:DATE])
  sp500DF[:,:returns] = similar(sp500DF[:,:Adj_Close])
  sp500DF[2:end,:returns] = log.(sp500DF[2:end,:Adj_Close] ./  sp500DF[1:(end-1),:Adj_Close])

  sp500DF[:,:dayOfQuarter] = zeros(Int,size(sp500DF,1))
  for i::Int ∈ 2:size(sp500DF,1)
    if Dates.firstdayofquarter(sp500DF[i,:DATE]) ≠ Dates.firstdayofquarter(sp500DF[i-1,:DATE])
      sp500DF[i,:dayOfQuarter] = 1
    else
      sp500DF[i,:dayOfQuarter] = sp500DF[i-1,:dayOfQuarter] + 1
    end
  end
  #showcols(sp500DF)

  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls", "w")
  serialize(stream, sp500DF)
  close(stream)

  return nothing
end

struct ΘBSM <: AbstractΘState
  S₀::Float64
  τ::Float64
  X::Float64
  r::Float64
  q::Float64
  σ::Float64
  σ²::Float64

  Φ::Function
end

#constructors for the BSM parameter set
ΘBSM(S₀::Float64, τ::Float64, X::Float64, r::Float64, q::Float64, σ::Float64)::ΘBSM =
  ΘBSM(S₀, τ, X, r, q, σ, σ*σ, ΦBSM)

ΘBSM(;S₀::Float64=0, τ::Float64=0, X::Float64=0, r::Float64=0, q::Float64=0, σ::Float64=0)::ΘBSM =
  ΘBSM(S₀, τ, X, r, q, σ)

#characteristic function of the black scholes
ΦBSM(Θ::ΘBSM, ω::T where T<:Union{Complex{Float64}, Float64})::Complex{Float64}  =
  exp(im*(ω * log(Θ.S₀) + ω * Θ.τ * (Θ.r - Θ.q - 0.5*Θ.σ²)) - 0.5 * Θ.σ²  *ω*ω * Θ.τ)

#parameter set for the BATES model
struct ΘBates <: AbstractΘState
  S₀::Float64
  τ::Float64
  X::Float64
  r::Float64
  q::Float64

  σ::Float64
  ρ::Float64
  v₀::Float64

  λ::Float64
  κ::Float64
  θ::Float64
  μⱼ::Float64
  σⱼ::Float64

  σ²::Float64
  σⱼ²::Float64
  Φ::Function
end

#constructors for the Bates parameter set
ΘBates(S₀::Float64, τ::Float64, X::Float64, r::Float64, q::Float64,
  σ::Float64, ρ::Float64, v₀::Float64, λ::Float64,
 κ::Float64, θ::Float64, μⱼ::Float64, σⱼ::Float64)::ΘBates =
  ΘBates(S₀, τ, X, r, q, σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ, σ*σ, σⱼ*σⱼ, ΦBates)

ΘBates(;S₀::Float64=0, τ::Float64=0, X::Float64=0, r::Float64=0, q::Float64=0,
  σ::Float64=0, ρ::Float64=0, v₀::Float64=0, λ::Float64=0,
  κ::Float64=0, θ::Float64=0, μⱼ::Float64=0, σⱼ::Float64=0)::ΘBates =
  ΘBates(S₀, τ, X, r, q, σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ)

#this is the characteristic function for Bates
function ΦBates(Θ::ΘBates, ω::T)::Complex{Float64} where T<:Union{Complex{Float64}, Float64}
  ι::Complex{Float64} = Ι

  #caclulate helper variables
  ρσiω::Complex{Float64} = Θ.ρ*Θ.σ*ι*ω
  d::Complex{Float64} = √((ρσiω-Θ.κ)^2+Θ.σ²*(ι*ω+ω*ω))
  g::Complex{Float64} = (Θ.κ-ρσiω-d)/(Θ.κ-ρσiω+d)

  #calculate next level of helper variables
  A::Complex{Float64} = ι * ω * log(Θ.S₀) + ι*ω*(Θ.r-Θ.q) * Θ.τ
  B::Complex{Float64} = Θ.θ*Θ.κ/Θ.σ² * ((Θ.κ-ρσiω-d)*Θ.τ-2.0*log((1.0-g*e^(-d*Θ.τ))/(1.0-g)))
  C::Complex{Float64} = Θ.v₀/Θ.σ² * (Θ.κ-ρσiω-d)*(1.0-e^(-d*Θ.τ)) / (1.0-g*e^(-d*Θ.τ))
  D::Complex{Float64} = -Θ.λ*Θ.μⱼ*ι*ω*Θ.τ + Θ.λ*Θ.τ*((1.0+Θ.μⱼ)^(ι*ω)*e^(0.5*Θ.σⱼ²*ι*ω*(ι*ω-1.0))-1.0)

  return e^(A+B+C+D)
end

function ΦtestBates(Θ::ΘBates, ω::T)::Complex{Float64} where T<:Union{Complex{Float64}, Float64}
  ι::Complex{Float64} = Ι

  #caclulate helper variables
  ρσiω::Complex{Float64} = Θ.ρ*Θ.σ*ι*ω
  d::Complex{Float64} = √((ρσiω-Θ.κ)^2+Θ.σ²*(ι*ω+ω*ω))
  g::Complex{Float64} = (Θ.κ-ρσiω-d)/(Θ.κ-ρσiω+d)

  #calculate next level of helper variables
  A::Complex{Float64} = ι * ω * log(Θ.S₀) + ι*ω*(Θ.r-Θ.q) * Θ.τ
  B::Complex{Float64} = Θ.θ*Θ.κ/Θ.σ² * ((Θ.κ-ρσiω-d)*Θ.τ-2.0*log((1.0-g*e^(-d*Θ.τ))/(1.0-g)))
  C::Complex{Float64} = Θ.v₀/Θ.σ² * (Θ.κ-ρσiω-d)*(1.0-e^(-d*Θ.τ)) / (1.0-g*e^(-d*Θ.τ))
  D::Complex{Float64} = -Θ.λ*Θ.μⱼ*ι*ω*Θ.τ + Θ.λ*Θ.τ*((1.0+Θ.μⱼ)^(ι*ω)*e^(0.5*Θ.σⱼ²*ι*ω*(ι*ω-1.0))-1.0)

  println("""Bates Intermediate Results:
     A: $A
     B: $B
     C: $C
     D: $D
     d: $d
     g: $g
     Answer: $(e^(A+B+C+D))""")

  return e^(A+B+C+D)
end

struct ΘHeston <: AbstractΘState
  S₀::Float64
  τ::Float64
  X::Float64
  r::Float64
  q::Float64

  σ::Float64
  ρ::Float64
  v₀::Float64

  κ::Float64
  θ::Float64

  σ²::Float64
  Φ::Function
end

#constructors for the Bates parameter set
ΘHeston(S₀::Float64, τ::Float64, X::Float64, r::Float64, q::Float64,
  σ::Float64, ρ::Float64, v₀::Float64, κ::Float64, θ::Float64)::ΘHeston =
  ΘHeston(S₀, τ, X, r, q, σ, ρ, v₀, κ, θ, σ*σ, ΦHeston)

ΘHeston(;S₀::Float64=0, τ::Float64=0, X::Float64=0, r::Float64=0, q::Float64=0,
  σ::Float64=0, ρ::Float64=0, v₀::Float64=0,
  κ::Float64=0, θ::Float64=0)::ΘHeston =
  ΘHeston(S₀, τ, X, r, q, σ, ρ, v₀, κ, θ)

#this is the characteristic function for Bates
function ΦHeston(Θ::ΘHeston, ω::T)::Complex{Float64} where T<:Union{Complex{Float64}, Float64}
  ι::Complex{Float64} = Ι

  #caclulate helper variables
  ρσiω::Complex{Float64} = Θ.ρ*Θ.σ*ι*ω
  d::Complex{Float64} = √((ρσiω-Θ.κ)^2+Θ.σ²*(ι*ω+ω*ω))
  g::Complex{Float64} = (Θ.κ-ρσiω-d)/(Θ.κ-ρσiω+d)

  #calculate next level of helper variables
  A::Complex{Float64} = ι * ω * log(Θ.S₀) + ι*ω*(Θ.r-Θ.q) * Θ.τ
  B::Complex{Float64} = Θ.θ*Θ.κ/Θ.σ² * ((Θ.κ-ρσiω-d)*Θ.τ-2.0*log((1.0-g*e^(-d*Θ.τ))/(1.0-g)))
  C::Complex{Float64} = Θ.v₀/Θ.σ² * (Θ.κ-ρσiω-d)*(1.0-e^(-d*Θ.τ)) / (1.0-g*e^(-d*Θ.τ))

  return e^(A+B+C)
end

function priceFromΦ(Θ::T)::Float64 where {T<:AbstractΘState}
  ι::Complex{Float64} = Ι

  #Integrand for Π1 (transformed to integrate from 0 to 1)
  function Φ1(ω::Float64)::Float64
    Ω::Float64 = ω/(1.0-ω)
    return real(e^(-ι * Ω * log(Θ.X)) * Θ.Φ(Θ,Ω - ι) / (ι * Ω * Θ.Φ(Θ,-ι))) / (1.0 - ω)^2
  end

  #Integrand for Π2 (transformed to integrate from 0 to 1)
  function Φ2(ω::Float64)::Float64
    Ω::Float64 = ω/(1.0-ω)
    return real(e^(-ι * Ω * log(Θ.X)) * Θ.Φ(Θ,Ω)/(ι * Ω)) / (1.0 - ω)^2
  end

  #do the integrations

  #implement with pure julia code
  Π1::Float64 = 0.5 + 1.0/π * I1D(Φ1, 0.0, 1.0)
  Π2::Float64 = 0.5 + 1.0/π * I1D(Φ2, 0.0, 1.0)

  #use the c implementation for a small performance boost
  #Π1::Float64 = 0.5 + 1.0/π * (hquadrature(Φ1, 0.0, 1.0))[1]
  #Π2::Float64 = 0.5 + 1.0/π * (hquadrature(Φ2, 0.0, 1.0))[1]
  #hquadrature(f, a, b)

  return Θ.S₀ * e^(-Θ.τ * Θ.q) * Π1 - Θ.X * e^(-Θ.τ * Θ.r) * Π2
end

#gets the price of a normal black scholes model
function priceBSMCall(S₀::Float64, τ::Float64, X::Float64, r::Float64, q::Float64, σ::Float64)::Float64
  d1::Float64 = log(S₀/X) + (r-2.0*q+σ^2/2.0)
  d1= d1 / (τ^0.5 * σ)

  return S₀*exp(-q*τ)*CTΦ(d1)-X*exp(-r*τ)*CTΦ(d1-σ*τ^0.5)
end

#implied volatility
function getImpliedσ(price::Float64, S₀::Float64, τ::Float64, X::Float64,
  r::Float64, q::Float64)::Float64

  #need to define special functions here for the solver

  #Algorithms: :LD_LBFGS, :LN_COBYLA, :LN_BOBYQA
  #  :LD_SLSQP, :LD_TNEWTON, :LD_VAR2/:LD_VAR1,
  #  :LD_MMA, :GN_ESCH, :GN_DIRECT, GN_DIRECT_L, :GD_STOGO
  # GN_ISRES, :LD_LBFGS, :LN_SBPLX,
  #GN_CRS2_LM
  #:GN_CRS2_LM***, :GN_ISRES*** (slow), :LD_MMA**, LN_SBPLX**, LD_TNEWTON** best:62,848

  function obj(σ::Real)::Real
    d1 = (log(S₀/X) + τ*(r-q+σ^2/2.0))/ (σ*τ^0.5)

    calcPrice = S₀*exp(-q*τ)*CTΦ(d1) - X*exp(-r*τ)*CTΦ(d1-σ*τ^0.5)
    return (price-calcPrice)^2
  end

  #works well:LN_NELDERMEAD , LN_SBPLX, LD_SLSQP, LN_BOBYQA
  mod = Model(solver=NLoptSolver(algorithm=:LN_BOBYQA, maxtime=1))
  #JuMP.register(mod, :priceFromVol, 1, priceFromVol, autodiff=true)
  JuMP.register(mod, :CTΦ, 1, CTΦ, autodiff=true)
  JuMP.register(mod, :obj, 1, obj, autodiff=true)

  #constrain vol to 0 and 1000%
  @variable(mod, 1.0 >= σ >= .0001, start = 0.2)

  #set the objective function
  @NLobjective(mod, Min, obj(σ))

  σBest::Float64  = 0.0
  objBest::Float64 = -9999.0

  #try
  output = solve(mod)
  σBest  = getvalue(σ)

  return σBest

end

#helper for getting the implied volatiltiy from a parameter obhect
getImpliedσ(price::Float64, Θ::T where T<:AbstractΘState) =
  getImpliedσ(price, Θ.S₀, Θ.τ, Θ.X, Θ.r, Θ.q)

#gets the implied volatility directly
function getImpliedσ(Θ::T where T<:AbstractΘState)::Float64
  price::Float64 = priceFromΦ(Θ)
  #println("price: $price")
  return getImpliedσ(price, Θ)
end

#script function for running the GARCH problem
function A5P2Script(;refreshData::Bool=true, doGraphs::Bool = true)
  if !isfile("$DATA_PATH\\$SP500_NAME.jls") || refreshData
    preProcessSP500()
  end

  #read the binary file
  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls")
  sp500DF::DataFrame = deserialize(stream)
  close(stream)

  returns::Vector{Float64} = sp500DF[2:end,:returns]

  S₀::Float64 = 100.0
  τ::Float64 = 0.33
  X::Float64 = 120.0
  r::Float64 = 0.04
  q::Float64 = 0.0
  σ::Float64 = 0.15

  testBSMOption::ΘBSM = ΘBSM(S₀=S₀, τ=τ, X=X, r=r, q=q, σ=σ)
  testBSMPrice::Float64 = priceFromΦ(testBSMOption)
  println("Price of BSM option: $testBSMPrice")

  testImpliedσ::Float64 = getImpliedσ(testBSMPrice, testBSMOption)
  println("Implied vol: $testImpliedσ")

  #test the parameters for the Bates model
  ρ::Float64 = -0.43
  v₀::Float64 = 0.15^2

  λ::Float64 = 1.764
  κ::Float64 = 0.011
  θ::Float64 = 0.09
  μⱼ::Float64 = -0.0301
  σⱼ::Float64 = 0.11

  testBatesOption::ΘBates = ΘBates(S₀=S₀, τ=τ, X=X, r=r, q=q,
    σ=σ, ρ=ρ, v₀=v₀, λ=λ, κ=κ, θ=θ, μⱼ=μⱼ, σⱼ=σⱼ)

  testBatesPrice::Float64 = priceFromΦ(testBatesOption)
  println("Price of Bates option: $testBatesPrice")

  testImpliedσ = getImpliedσ(testBatesOption)
  println("Implied vol: $testImpliedσ")

  ΦtestBates(testBatesOption, 1.0)

  testHestonOption::ΘHeston = ΘHeston(S₀=S₀, τ=τ, X=X, r=r, q=q,
    σ=σ, ρ=ρ, v₀=v₀, κ=κ, θ=θ)

  testHestonPrice::Float64 = priceFromΦ(testHestonOption)
  println("Price of Heston option: $testHestonPrice")

  testImpliedσ = getImpliedσ(testHestonOption)
  println("Implied vol: $testImpliedσ")

  #this sets the ranges
  if doGraphs
    τs::Vector{Float64} = collect(0.01:0.01:5.0)
    Xs::Vector{Float64} = collect(50:0.1:150)

    Nτ::Int = length(τs)
    NX::Int = length(Xs)

    #setup the option parameters
    BatesVsτ::Vector{ΘBates} = ((t::Float64)->ΘBates(S₀, t, X, r, q,
      σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ)).(τs)
    BatesVsX::Vector{ΘBates} = ((x::Float64)->ΘBates(S₀, τ, x, r, q,
      σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ)).(Xs)

    HestonVsτ::Vector{ΘHeston} = ((t::Float64)->ΘHeston(S₀, t, X, r, q,
      σ, ρ, v₀, κ, θ)).(τs)
    HestonVsX::Vector{ΘHeston} = ((x::Float64)->ΘHeston(S₀, τ, x, r, q,
      σ, ρ, v₀, κ, θ)).(Xs)

    #pre-allocate the implied volatilities
    BatesVsτIV::Vector{Float64} = Vector{Float64}(length(τs))
    BatesVsXIV::Vector{Float64} = Vector{Float64}(length(Xs))
    HestonVsτIV::Vector{Float64} = Vector{Float64}(length(τs))
    HestonVsXIV::Vector{Float64} = Vector{Float64}(length(Xs))

    println("Getting $Nτ τ estimates")
    Threads.@threads for i::Int ∈ 1:Nτ
      BatesVsτIV[i] = getImpliedσ(BatesVsτ[i])
      HestonVsτIV[i] = getImpliedσ(HestonVsτ[i])
    end

    println("Getting $NX X estimates")
    Threads.@threads for i::Int ∈ 1:NX
      BatesVsXIV[i] = getImpliedσ(BatesVsX[i])
      HestonVsXIV[i] = getImpliedσ(HestonVsX[i])
    end

    #form the data so that it can be graphed
    τDF::DataFrame = DataFrame(τ=τs, σ=BatesVsτIV)
    τDF[:,:model] = :Bates
    τDF = [τDF; DataFrame(DataFrame(τ=τs, σ=HestonVsτIV, model=[:Heston for i∈1:Nτ]))]

    XDF::DataFrame = DataFrame(X=Xs, σ=BatesVsXIV)
    XDF[:,:model] = :Bates
    XDF = [XDF; DataFrame(DataFrame(X=Xs, σ=HestonVsXIV, model=[:Heston for i∈1:NX]))]

    #get the graphs
    Gadfly.push_theme(:default)
    plots::Vector{PlotContainer} = Vector{PlotContainer}()
    plotNames::Vector{String} = Vector{String}()

    push!(plotNames, "TimeSmile")
    push!(plots, plot(τDF, x=:τ, y=:σ, color=:model,
      Geom.line, style(key_position=:right),
      Guide.ylabel("Implied σ"), Guide.xlabel("Time to Maturity"), Coord.Cartesian(ymin=0.0,ymax=1.0, xmin=0.0, xmax=5.0),
      Guide.title("Volatility Smile by Time (μⱼ=-0.0301, θ=.09, σⱼ=0.11)")))

    push!(plotNames, "MoneySmile")
    push!(plots, plot(XDF, x=:X, y=:σ, color=:model,
       Geom.line, style(key_position=:right),
      Guide.ylabel("Implied σ"), Guide.xlabel("X Given S₀=100"), Coord.Cartesian(ymin=0.0,ymax=1.0, xmin=50, xmax=150),
      Guide.title("Volatility Smile by Moneyness (μⱼ=-0.0301, θ=.09, σⱼ=0.11))")))

    ####Now repeat for a wider smile
    μⱼ = -0.20
    σ = 0.4
    σⱼ = 0.25
    v₀ = 0.2^2

    #get the new option models
    BatesVsτ = ((t::Float64)->ΘBates(S₀, t, X, r, q,
      σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ)).(τs)
    BatesVsX = ((x::Float64)->ΘBates(S₀, τ, x, r, q,
      σ, ρ, v₀, λ, κ, θ, μⱼ, σⱼ)).(Xs)

    HestonVsτ = ((t::Float64)->ΘHeston(S₀, t, X, r, q,
      σ, ρ, v₀, κ, θ)).(τs)
    HestonVsX = ((x::Float64)->ΘHeston(S₀, τ, x, r, q,
      σ, ρ, v₀, κ, θ)).(Xs)

    println("Getting $Nτ τ estimates")
    Threads.@threads for i::Int ∈ 1:Nτ
      BatesVsτIV[i] = getImpliedσ(BatesVsτ[i])
      HestonVsτIV[i] = getImpliedσ(HestonVsτ[i])
    end

    println("Getting $NX X estimates")
    Threads.@threads for i::Int ∈ 1:NX
      BatesVsXIV[i] = getImpliedσ(BatesVsX[i])
      HestonVsXIV[i] = getImpliedσ(HestonVsX[i])
    end

    #form the data so that it can be graphed
    τDF = DataFrame(τ=τs, σ=BatesVsτIV)
    τDF[:,:model] = :Bates
    τDF = [τDF; DataFrame(DataFrame(τ=τs, σ=HestonVsτIV, model=[:Heston for i∈1:Nτ]))]

    XDF = DataFrame(X=Xs, σ=BatesVsXIV)
    XDF[:,:model] = :Bates
    XDF = [XDF; DataFrame(DataFrame(X=Xs, σ=HestonVsXIV, model=[:Heston for i∈1:NX]))]

    push!(plotNames, "TimeSmile Enhanced")
    push!(plots, plot(τDF, x=:τ, y=:σ, color=:model,
      Geom.line, #Theme(line_width = 0.25pt),
      Guide.ylabel("Implied σ"), Guide.xlabel("Time to Maturity"), Coord.Cartesian(ymin=0.0,ymax=1.0, xmin=0.0, xmax=5.0),
      Guide.title("Volatility Smile by Time (μⱼ=-0.15, σ=.4, σⱼ=0.25, v₀ = 0.2²)")))

    push!(plotNames, "MoneySmile Enhanced")
    push!(plots, plot(XDF, x=:X, y=:σ, color=:model,
       Geom.line, #Theme(line_width = 0.25pt),
      Guide.ylabel("Implied σ"), Guide.xlabel("X Given S₀=100"), Coord.Cartesian(ymin=0.0,ymax=1.0, xmin=50, xmax=150),
      Guide.title("Volatility Smile by Moneyness (μⱼ=-0.15, σ=.4, σⱼ=0.25, v₀ = 0.2²))")))

      for i ∈ 1:length(plots)
        draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 7.5inch, 9inch),plots[i])
      end
  end

end

@time begin
  A5P2Script(refreshData=false, doGraphs=true)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# %% Chaotic scattering
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using PyPlot, DynamicalBilliards
fig, axs = subplots(1,3)
axs[3].axis("off")
ax3 = fig.add_subplot(1, 6, 5)
ax4 = fig.add_subplot(1, 6, 6)
DynamicalBilliards.obcolor(::Obstacle) = matplotlib.colors.to_rgb(COLORS[3])

# First plot: illustration of scattering function
# I GIVE UP, I'll Make this plot in PowerPoint
off = 0.3
r = 0.2
offset = [0.0, off]
center = 2.5
R(φ) = [cos(φ) -sin(φ);
        sin(φ)  cos(φ)]
R3 = R(2π/3)
enclosing = billiard_sinai(r, 2center, 2center)
enclosing = enclosing[2:5]

# disk0 = Disk([center,0], 1.5r, "0")
# sca(axs[1])
# plot(disk0)
# axs[1].set_xlim(-1.5 + center, 0.5 + center)
# axs[1].set_ylim(-1+0.3, 1+0.3)
axs[1].axis("off")

# p0 = Particle(center-1, 0.22, 0)
# bd0 = Billiard(enclosing..., disk0)
# x,y,vx,vy = timeseries(p0, bd0, 2; dt = 0.01)
# axs[1].plot(x[1:130],y[1:130]; color = "C1")
# x0, y0 = p0.pos
# vx0, vy0 = p0.vel
# axs[1].quiver(x0, y0, 0.08vx0, 0.08vy0; angles = "xy", scale = 0.5, width = 0.018, color="C1", zorder = 99)
# axs[1].text(x0, y0-0.2, "\$\\phi\$")
# axs[1].text(x0, y0-0.2, "\$a\$")

# Okay now clarify the scattering function
disk1 = Disk([center, center] .+ offset, r, "green")
disk2 = Disk([center, center] .+ R3*offset, r, "red")
disk3 = Disk([center, center] .+ R3*R3*offset, r, "purple")

disks = Billiard(disk1, disk2, disk3)
plot(disks; ax = axs[2])
axs[2].set_xlim(1.9,3.1)
axs[2].set_ylim(2.0,3.2)
scattering = Billiard(enclosing..., disk1, disk2, disk3)

axs[2].axis("off")

# make some particles

terminate(n, τ, i, p) = i ∈ 1:4
φs = 5π/6 .+ range(0; step = 0.001, length = 3) .+ 0.17

ps = [Particle((R(φ)*[0.5, 0] .+ [center, center])..., φ+π) for φ in φs]

partcol = ( "C3", "C1", "C0",)
for (i, p) in enumerate(ps);
    x0, y0 = p.pos
    vx0, vy0 = p.vel
    axs[2].quiver(x0, y0, 0.04vx0, 0.04vy0; angles = "xy", scale = 0.5, width = 0.01,
    color = partcol[i], zorder = 99)
    x,y = timeseries(p, scattering, terminate)
    axs[2].plot(x, y; color = partcol[i], lw = 2, alpha = 0.75)
end

# Detailed computation of input-output function
φs = range(0, 2π/3; length = 100_000)
θs = zero(φs)

for (i, φ) in enumerate(φs);
    p = Particle((R(φ)*[0.5, 0] .+ [center, center])..., φ+π)
    timeseries!(p, scattering, terminate)
    θs[i] = atan(p.vel[2], p.vel[1])
end

ax3.plot(φs, θs .+ π; ls = "None", marker = "o", ms = 0.5, alpha = 0.5)
ax3.set_xticks([0, 2π/3])
ax3.set_xticklabels(["0", "2π/3"])
ax3.set_yticks([0, 2π])
ax3.set_yticklabels(["0", "2π"])
ax4.set_yticks([])
ax3.set_xlim(0, 2π/3)
ax4.set_ylim(0, 2π)
ax3.set_ylim(0, 2π)

φs = range(0.1555, 0.1567; length = 100_000)
θs = zero(φs)
for (i, φ) in enumerate(φs);
    p = Particle((R(φ)*[0.5, 0] .+ [center, center])..., φ+π)
    timeseries!(p, scattering, terminate)
    θs[i] = atan(p.vel[2], p.vel[1])
end
ax4.plot(φs, θs .+ π; ls = "None", marker = "o", ms = 0.5, alpha = 0.5)
ax4.set_xlim(φs[1], φs[end])
ax4.set_xticks([])
ax3.set_xlabel("\$\\phi_i\$", labelpad = -25)
ax3.set_ylabel("\$\\phi_o = \\mathbf{S}(\\phi_i)\$", labelpad = -25)
# ax3.axvspan(φs[1], φs[end]; color = "C1") # not even visible lol

fig.tight_layout(;pad = 0.25)
fig.subplots_adjust(wspace = 0.1)

# wsave(plotsdir("8", "chaoticscattering"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
##################################################################################################
## Unmatched test
##################################################################################################
"""
    `comp_unmat_stat_cmd(mods_g1,mods_g2)`
    Function that computes a test statistic Tcmd for each analysis region in unmatched case
    between groups.
    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_stat_cmd(mods_g1,mods_g2)
    ```
"""
function comp_unmat_stat_cmd(mods_g1::Vector{RegStruct}, mods_g2::Vector{RegStruct})::Vector{Float64}

    # Compute statistic
    cmds = []
    @inbounds for s1 = 1:length(mods_g1)
        @inbounds for s2 = 1:length(mods_g2)
            push!(cmds, comp_cmd(mods_g1[s1], mods_g2[s2]))
        end
    end

    # Return vector of average CMD between groups
    return sum(cmds) / (length(mods_g1) * length(mods_g2))

end
"""
    `get_unmat_cmd_tbl(mods_g1,mods_g2)`

    Function that computes CMD between all possible pairs.

    # Examples
    ```julia-repl
    julia> CpelNano.get_unmat_cmd_tbl(mods_g1,mods_g2)
    ```
"""
function get_unmat_cmd_tbl(mods_g1::Vector{RegStruct}, mods_g2::Vector{RegStruct})::Dict{NTuple{2,Int64},Vector{Float64}}

    # Number of analysis regions
    rs = mods_g1[1]
    num_nls_reg = rs.nls_rgs.num

    # Number of replicates
    n_g1 = length(mods_g1)
    n_g2 = length(mods_g2)
    n_tot = n_g1 + n_g2

    # Init matrix
    cmd_tbl = Dict{NTuple{2,Int64},Vector{Float64}}()

    # Compute statistic within group 1
    @inbounds for s1 = 1:n_g1
        @inbounds for s2 = (s1+1):n_g1
            cmd_tbl[(s1,s2)] = comp_cmd(mods_g1[s1], mods_g1[s2])
        end
    end

    # Compute statistic within group 2
    @inbounds for s1 = 1:n_g2
        @inbounds for s2 = (s1+1):n_g2
            cmd_tbl[(n_g1+s1,n_g1+s2)] = comp_cmd(mods_g2[s1], mods_g2[s2])
        end
    end

    # Compute statistic between groups
    @inbounds for s1 = 1:n_g1
        @inbounds for s2 = 1:n_g2
            cmd_tbl[(s1,n_g1+s2)] = comp_cmd(mods_g1[s1], mods_g2[s2])
        end
    end

    # Return matrix of CMD
    return cmd_tbl

end
"""
    `comp_unmat_perm_stat_cmd(N_G1,N_G2,NUM_NLS_RGS,CMD_TBL,PERM_IDS)`

    Function that produces a permutation statistic for Tcmd in unmatched case.

    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_perm_stat_cmd(n_g1,n_g2,num_nls_rgs,cmd_tbl,ind_g1)
    ```
"""
function comp_unmat_perm_stat_cmd(n_g1::Int64,n_g2::Int64,num_nls_rgs::Int64,cmd_tbl::Dict{Tuple{Int64,Int64},Array{Float64,1}},ind_g1::Vector{Int64})::Vector{Float64}

    # Replicates
    ind_g2 = collect(1:(n_g1 + n_g2))
    ind_g2 = deleteat!(ind_g2,ind_g1)

    # Get statistic
    cmd = fill(0.0,num_nls_rgs)
    @inbounds for i in ind_g1
        @inbounds for j in ind_g2
            if haskey(cmd_tbl,(i,j))
                cmd .+= cmd_tbl[(i,j)]
            else
                cmd .+= cmd_tbl[(j,i)]
            end
        end
    end

    # Return stat for permutation
    return cmd / (n_g1 * n_g2)

end
"""
    `comp_unmat_stat_mml(μs_g1,μs_g2)`

    Function that computes MML difference between groups (unmatched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_stat_mml(μs_g1,μs_g2)
```
"""
function comp_unmat_stat_mml(μs_g1::Vector{Vector{Float64}}, μs_g2::Vector{Vector{Float64}})::Vector{Float64}

    # Compute means within groups
    mml1 = sum(μs_g1) / length(μs_g1)
    mml2 = sum(μs_g2) / length(μs_g2)

    # Return Tmml
    return mml1 - mml2

end
"""
    `comp_unmat_perm_stat_mml(μs_g1,μs_g2,PERM_IDS)`

    Function that produces a permutation statistic for Tmml (unmatched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_perm_stat_mml(μs_g1, μs_g2, perm_ids)
```
"""
function comp_unmat_perm_stat_mml(μs_g1::Vector{Vector{Float64}}, μs_g2::Vector{Vector{Float64}}, perm_ids::Vector{Int64})::Vector{Float64}

    # Get vectors for each group
    μs_g2p = vcat(μs_g1, μs_g2)
    μs_g1p = μs_g2p[perm_ids]
    deleteat!(μs_g2p, perm_ids)

    # Return
    return comp_unmat_stat_mml(μs_g1p, μs_g2p)

end
"""
    `comp_unmat_stat_nme(hs_g1,hs_g2)`

    Function that computes NME difference between groups (unmatched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_stat_nme(hs_g1,hs_g2)
```
"""
function comp_unmat_stat_nme(hs_g1::Vector{Vector{Float64}}, hs_g2::Vector{Vector{Float64}})::Vector{Float64}

    # Compute h within group
    nme1 = sum(hs_g1) / length(hs_g1)
    nme2 = sum(hs_g2) / length(hs_g2)

    # Return difference in nme
    return nme1 - nme2

end
"""
    `comp_unmat_perm_stat_nme(hs_g1,hs_g2,PERM_IDS)`

    Function that produces a permutation statistic for Tnme (unmatched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_unmat_perm_stat_nme(hs_g1,hs_g2,perm_ids)
```
"""
function comp_unmat_perm_stat_nme(hs_g1::Vector{Vector{Float64}}, hs_g2::Vector{Vector{Float64}}, perm_ids::Vector{Int64})::Vector{Float64}

    # Get vectors for each group
    hs_g2p = vcat(hs_g1, hs_g2)
    hs_g1p = hs_g2p[perm_ids]
    deleteat!(hs_g2p, perm_ids)

    # Return
    return comp_unmat_stat_nme(hs_g1p, hs_g2p)

end
"""
    `unmat_nls_reg_test(MODELS_G1,MODELS_G2)`

    Function that performs hypothesis testing in unmatched samples group comparison in an estimation region.

    # Examples
    ```julia-repl
    julia> CpelNano.unmat_nls_reg_test(ms_g1,ms_g2)
```
"""
function unmat_est_reg_test(ms_g1::Vector{RegStruct}, ms_g2::Vector{RegStruct})::RegStatTestStruct

    # Initialize output object
    test_struct = RegStatTestStruct(ms_g1[1])

    # Get stats for k-th analysis region
    μs_g1 = [x.mml for x in ms_g1]
    μs_g2 = [x.mml for x in ms_g2]
    hs_g1 = [x.nme for x in ms_g1]
    hs_g2 = [x.nme for x in ms_g2]

    # Compute observed stats
    tmml_obs = comp_unmat_stat_mml(μs_g1, μs_g2)
    tnme_obs = comp_unmat_stat_nme(hs_g1, hs_g2)
    tcmd_obs = comp_unmat_stat_cmd(ms_g1, ms_g2)

    # Init p-values
    tmml_pvals = fill(NaN, length(tmml_obs))
    tnme_pvals = fill(NaN, length(tnme_obs))
    tcmd_pvals = fill(NaN, length(tcmd_obs))

    # Compute number of possible randomizations
    n_g1 = length(μs_g1)
    n_g2 = length(μs_g2)
    L = binomial(n_g1 + n_g2, n_g1)

    # If enough data compute p-values (p<0.05)
    if L>20

        # Check if exact test
        exact = L < LMAX

        # Create iteratable object with all combinations
        comb_iter = combinations(1:(n_g1+n_g2), n_g1)

        # Get group label combinations to use
        comb_iter_used = []
        if exact
            # Use all group assignments
            comb_iter_used = comb_iter
        else
            # Use Lmax group assignments
            ind_subset = rand(1:L, LMAX)
            @inbounds for (ind, comb) in enumerate(comb_iter)
                (ind in ind_subset) && push!(comb_iter_used, comb)
            end
        end

        # Compute all CMD pairs
        cmd_tbl = get_unmat_cmd_tbl(ms_g1,ms_g2)

        # Use method for random permutation
        tmml_perms = map(x -> comp_unmat_perm_stat_mml(μs_g1, μs_g2, x), comb_iter_used)
        tnme_perms = map(x -> comp_unmat_perm_stat_nme(hs_g1, hs_g2, x), comb_iter_used)
        tcmd_perms = map(x -> comp_unmat_perm_stat_cmd(n_g1, n_g2, length(tmml_obs), cmd_tbl, x), comb_iter_used)

        # Compute p-values two-sided test
        @inbounds for k in 1:length(tmml_obs)

            # Check if data
            isnan(tmml_obs[k]) && continue

            # Get permutation stats from k-th analysis region
            tmml_perms_k = [perm[k] for perm in tmml_perms]
            tnme_perms_k = [perm[k] for perm in tnme_perms]
            tcmd_perms_k = [perm[k] for perm in tcmd_perms]

            # Get number of permutation stats equal or above observed
            tmml_pval_k = sum(abs.(tmml_perms_k) .>= abs(tmml_obs[k]))
            tnme_pval_k = sum(abs.(tnme_perms_k) .>= abs(tnme_obs[k]))
            tcmd_pval_k = sum(tcmd_perms_k .>= tcmd_obs[k])

            # Compute p-values
            tmml_pvals[k] = exact ? tmml_pval_k / length(tmml_perms_k) : (1.0 + tmml_pval_k) / (1.0 + length(tmml_perms_k))
            tnme_pvals[k] = exact ? tnme_pval_k / length(tnme_perms_k) : (1.0 + tnme_pval_k) / (1.0 + length(tnme_perms_k))
            tcmd_pvals[k] = exact ? tcmd_pval_k / length(tcmd_perms_k) : (1.0 + tcmd_pval_k) / (1.0 + length(tcmd_perms_k))

        end

    end

    # Fill return object
    @inbounds for k in 1:length(tmml_obs)
        test_struct.tests.tmml_test[k] = (tmml_obs[k],tmml_pvals[k])
        test_struct.tests.tnme_test[k] = (tnme_obs[k],tnme_pvals[k])
        test_struct.tests.tcmd_test[k] = (tcmd_obs[k],tcmd_pvals[k])
    end

    # Return output structure
    return test_struct

end
##################################################################################################
## Matched test
##################################################################################################
"""
    `comp_mat_stat_tcmd(MODS_G1,MODS_G2)`

    Function that computes CMD between matched pairs in two groups (matched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_mat_stat_tcmd(mods_g1,mods_g2)
    ```
"""
function comp_mat_stat_tcmd(mods_g1::Vector{RegStruct}, mods_g2::Vector{RegStruct})::Vector{Float64}

    # Compute stat
    cmds = []
    @inbounds for s = 1:length(mods_g1)
        push!(cmds, comp_cmd(mods_g1[s], mods_g2[s]))
    end

    # Return mean CMD
    return sum(cmds) / length(mods_g1)

end
"""
    `mat_reg_test_tcmd(mods_g1,mods_g2)`

    Function that performs hypothesis testing in matched samples group comparison for Tcmd.

    # Examples
    ```julia-repl
    julia> CpelNano.mat_reg_test_tcmd(mods_g1,mods_g2)
    ```
"""
function mat_reg_test_tcmd(mods_g1::Vector{RegStruct}, mods_g2::Vector{RegStruct})::NTuple{2,Vector{Float64}}

    # Return stat-pval pair
    return (comp_mat_stat_tcmd(mods_g1, mods_g2), fill(NaN, length(mods_g1)))

end
"""
    `comp_mat_diff_mml(MMLs_G1,MMLs_G2)`

    Function that computes MML differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_mat_diff_mml(μs_g1,μs_g2)
    ```
"""
function comp_mat_diff_mml(μs_g1::Vector{Vector{Float64}}, μs_g2::Vector{Vector{Float64}})::Vector{Vector{Float64}}

    # Compute mean differences
    diffs = []
    @inbounds for s = 1:length(μs_g1)
        push!(diffs, μs_g1[s] - μs_g2[s])
    end

    # Return vector of differences
    return diffs

end
"""
    `comp_mat_diff_nme(NMEs_G1,NMEs_G2)`

    Function that computes NME differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_mat_diff_mml(hs_g1,hs_g2)
    ```
"""
function comp_mat_diff_nme(hs_g1::Vector{Vector{Float64}}, hs_g2::Vector{Vector{Float64}})::Vector{Vector{Float64}}

    # Compute entropy differences
    diffs = []
    @inbounds for s = 1:length(hs_g1)
        push!(diffs, hs_g1[s] - hs_g2[s])
    end

    # Return vector of differences
    return diffs

end
"""
    `comp_mat_j_stat(diffs,j)`

    Function that computes permutation statistic for j-th sign assigment given vector of differences
    between matched pairs (matched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_mat_j_stat(diffs,j)
    ```
"""
function comp_mat_j_stat(diffs::Vector{Vector{Float64}}, j::Int64)::Vector{Float64}

    # Init diff in permutation
    diffs_perm = fill(0.0, length(diffs[1]))

    # Change sign if pertinent
    @inbounds for k = 1:length(diffs[1])
        j_k = j
        @inbounds for s = 1:length(diffs)
            diffs_perm[k] += Bool(j_k & 1) ? diffs[s][k] : -diffs[s][k]
            j_k >>= 1
        end
    end

    # Return permutation statistic
    return diffs_perm / length(diffs)

end
"""
    `comp_mat_perm_stats(DIFFs,Js)`

    Function that computes permutation statistics given a vector of differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> CpelNano.comp_mat_perm_stats(diffs,js)
    ```
"""
function comp_mat_perm_stats(diffs::Vector{Vector{Float64}}, js::Vector{Int64})::Vector{Vector{Float64}}

    # Return all possible signed sums
    return [comp_mat_j_stat(diffs, j) for j in js]

end
"""
    `mat_est_reg_test(MODELS_G1,MODELS_G2)`

    Function that performs hypothesis testing in unmatched samples group comparison in an estimation region.

    # Examples
    ```julia-repl
    julia> CpelNano.mat_est_reg_test(ms_g1,ms_g2)
    ```
"""
function mat_est_reg_test(ms_g1::Vector{RegStruct}, ms_g2::Vector{RegStruct})::RegStatTestStruct

    # Initialize output object
    test_struct = CpelNano.RegStatTestStruct(ms_g1[1])

    # Get stats for k-th analysis region
    μs_g1 = [m.mml for m in ms_g1]
    μs_g2 = [m.mml for m in ms_g2]
    hs_g1 = [m.nme for m in ms_g1]
    hs_g2 = [m.nme for m in ms_g2]

    # Compute observed differences
    mml_diffs = CpelNano.comp_mat_diff_mml(μs_g1, μs_g2)
    nme_diffs = CpelNano.comp_mat_diff_nme(hs_g1, hs_g2)

    # Compute observed stats
    tmml_obs = sum(mml_diffs) / length(ms_g1)
    tnme_obs = sum(nme_diffs) / length(ms_g1)
    tcmd_obs = CpelNano.comp_mat_stat_tcmd(ms_g1, ms_g2)

    # Init p-values
    tmml_pvals = fill(NaN, length(tmml_obs))
    tnme_pvals = fill(NaN, length(tnme_obs))

    # If enough data compute p-values
    if (0.5^length(μs_g1))<0.05

        # Get group label combinations to use
        exact = 2^length(ms_g1) < LMAX
        js = exact ? collect(0:2^length(ms_g1) - 1) : rand(0:2^length(ms_g1) - 1, LMAX)

        # Compute permutation stats
        tmml_perms = CpelNano.comp_mat_perm_stats(mml_diffs, js)
        tnme_perms = CpelNano.comp_mat_perm_stats(nme_diffs, js)

        # Compute p-values two-sided test
        @inbounds for k in 1:length(tmml_obs)

            # Check if data
            isnan(tmml_obs[k]) && continue

            # Get permutation stats from k-th analysis region
            tmml_perms_k = [perm[k] for perm in tmml_perms]
            tnme_perms_k = [perm[k] for perm in tnme_perms]

            # Get number of permutation stats equal or above observed
            tmml_pval_k = sum(abs.(tmml_perms_k) .>= abs(tmml_obs[k]))
            tnme_pval_k = sum(abs.(tnme_perms_k) .>= abs(tnme_obs[k]))

            # Compute p-values
            tmml_pvals[k] = exact ? tmml_pval_k / length(tmml_perms_k) : (1.0 + tmml_pval_k) / (1.0 + length(tmml_perms_k))
            tnme_pvals[k] = exact ? tnme_pval_k / length(tnme_perms_k) : (1.0 + tnme_pval_k) / (1.0 + length(tnme_perms_k))

        end

    end

    # Fill return object
    @inbounds for k in 1:length(tmml_pvals)
        test_struct.tests.tmml_test[k] = (tmml_obs[k], tmml_pvals[k])
        test_struct.tests.tnme_test[k] = (tnme_obs[k], tnme_pvals[k])
        test_struct.tests.tcmd_test[k] = (tcmd_obs[k], NaN)
    end

    # Return output structure
    return test_struct

end
##################################################################################################
## User functions
##################################################################################################
"""
    `read_in_grp_mods(MOD_FILES,CHR,FA_REC,CONFIG)`

    Function that reads in model files of group in chromosome `CHR`.

    # Examples
    ```julia-repl
    julia> CpelNano.read_in_grp_mods(mod_fls,chr,fa_rec,config)
    ```
"""
function read_in_grp_mods(mod_fls::Vector{String}, chr::String, fa_rec::FASTA.Record, config::CpelNanoConfig)::Vector{Dict{String,RegStruct}}

    # Read in models
    ms = Vector{Dict{String,RegStruct}}()
    for f in mod_fls
        push!(ms, read_mod_file_chr(f, chr, fa_rec, config))
    end

    # Return models
    return ms

end
"""
    `find_comm_regs(MODS_G1,MODS_G2)`

    Function that finds common regions between the two groups.

    # Examples
    ```julia-repl
    julia> CpelNano.find_comm_regs(ms_g1,ms_g2)
    ```
"""
function find_comm_regs(ms_g1::Vector{Dict{String,RegStruct}}, ms_g2::Vector{Dict{String,RegStruct}})::Vector{String}

    # Deal with group 1
    g1_dic = Dict{String,Int64}()
    for rep in ms_g1
        for k in keys(rep)
            g1_dic[k] = k in keys(g1_dic) ? g1_dic[k] + 1 : 1
        end
    end

    # Deal with group 2
    g2_dic = Dict{String,Int64}()
    for rep in ms_g2
        for k in keys(rep)
            g2_dic[k] = k in keys(g2_dic) ? g2_dic[k] + 1 : 1
        end
    end

    # Common regions
    comm_regs = Vector{String}()
    for k in keys(g1_dic)
        k in keys(g2_dic) && push!(comm_regs, k)
    end

    # Return common regions
    return comm_regs

end
"""
    `pmap_diff_grp_comp(REG_ID,MODS_G1,MODS_G2,CONFIG)`

    Function that performs differential methylation analysis between groups given two groups of
    CpelNano model files at a given estimation region `REG_ID`.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_diff_grp_comp(reg_id,mods_g1,mods_g2,config)
    ```
"""
function pmap_diff_grp_comp(reg_id::String, ms_g1::Vector{Dict{String,RegStruct}}, ms_g2::Vector{Dict{String,RegStruct}}, config::CpelNanoConfig)::RegStatTestStruct

    ## Perform testing

    print_log("Working on: $(reg_id)")

    # Matched vs Unmatched
    if config.matched

        ## Get available pairs for estimation region
        reg_ms_g1 = Vector{RegStruct}()
        reg_ms_g2 = Vector{RegStruct}()
        for rep_ind in 1:length(ms_g1)
            if haskey(ms_g1[rep_ind],reg_id) && haskey(ms_g2[rep_ind],reg_id)
                # Get sample dictionaries
                rep_dic_g1 = ms_g1[rep_ind]
                rep_dic_g2 = ms_g2[rep_ind]
                # Push pair
                push!(reg_ms_g1, rep_dic_g1[reg_id])
                push!(reg_ms_g2, rep_dic_g2[reg_id])
            end
        end

        # Matched test
        out_test = mat_est_reg_test(reg_ms_g1, reg_ms_g2)

    else

        ## Get models for specific region

        # Group 1
        reg_ms_g1 = Vector{RegStruct}()
        for dict in ms_g1
            haskey(dict, reg_id) && push!(reg_ms_g1, dict[reg_id])
        end

        # Group 2
        reg_ms_g2 = Vector{RegStruct}()
        for dict in ms_g2
            haskey(dict, reg_id) && push!(reg_ms_g2, dict[reg_id])
        end

        ## Unmatched test
        out_test = unmat_est_reg_test(reg_ms_g1, reg_ms_g2)

    end

    # Return test struct
    return out_test

end
"""
    `diff_grp_comp(MODEL_FILES_G1,MODEL_FILES_G2,FASTA,CONFIG)`

    Function that performs differential methylation analysis between groups given two groups of
    CpelNano model files.

    # Examples
    ```julia-repl
    julia> CpelNano.diff_grp_comp(mod_fls_g1,mod_fls_g2,fasta,config)
    ```
"""
function diff_grp_comp(mod_fls_g1::Vector{String}, mod_fls_g2::Vector{String}, fasta::String, config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names, chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Get output file names
    config.out_diff_files = OutputDiffFiles(config.out_dir, config.out_prefix)

    # Check if output file exists
    check_diff_output_exists(config) && return nothing

    # Loop over chromosomes
    for (i, chr) in enumerate(chr_names)

        print_log("Processing chromosome $(chr)")

        ## Get fasta record

        fa_rec = get_chr_fa_rec(chr, fasta)

        ## Read in models

        ms_g1 = read_in_grp_mods(mod_fls_g1, chr, fa_rec, config)
        ms_g2 = read_in_grp_mods(mod_fls_g2, chr, fa_rec, config)

        ## Find regions with at least one sample per group

        comm_regs = find_comm_regs(ms_g1, ms_g2)

        # If no data in chromosome continue
        length(comm_regs) > 0 || continue
        config.verbose && print_log("Comm regions: $(comm_regs)")

        ## Process each region
        pmap_out = pmap(x -> pmap_diff_grp_comp(x, ms_g1, ms_g2, config), comm_regs)

        ## Write chromosome output
        write_diff_output(pmap_out, config)

    end

    # Perform multiple hypothesis testing
    print_log("Multiple hypothesis testing")
    mult_hyp_corr(config)

    # Return nothing
    return nothing

end


text:
module Kokkos
  using CxxWrap
  using MPI
  import ..libjltrilinos

  using CustomUnitRanges: filename_for_zerorange
  include(filename_for_zerorange)

  @wrapmodule(libjltrilinos, :register_kokkos)

  """
  Expose a view using its raw data pointer
  """
  struct PtrView{ScalarT,N,LayoutT} <: AbstractArray{ScalarT,N}
    ptr::Ptr{ScalarT}
    size::NTuple{N,Int}
    stored_view # prevent GC of the actual view
  end

  make_dimensions(v,::Type{Val{N}}) where {N} = ([Int(Kokkos.dimension(v,i-1)) for i in 1:N]...,)

  function view(v::View3{Ptr{Ptr{ST}},LayoutT,SpaceT}, ndims::Type{Val{N}}) where {ST,LayoutT,SpaceT,N}
    return PtrView{ST,N,LayoutT}(Kokkos.ptr_on_device(v),make_dimensions(v,ndims), v)
  end

  function view(v::View4{Ptr{Ptr{ST}},LayoutT,SpaceT,Nothing}, ndims::Type{Val{N}}) where {ST,LayoutT,SpaceT,N}
    return PtrView{ST,N,LayoutT}(Kokkos.ptr_on_device(v),make_dimensions(v,ndims), v)
  end

  # Array interface, 0-based for consistency with global-to-local index mapping
  Base.IndexStyle(::Type{PtrView{ScalarT,N,LayoutT}}) where {ScalarT,N,LayoutT} = IndexLinear()
  Base.axes(v::PtrView{ScalarT,N,LayoutLeft}) where {ScalarT,N} = map(ZeroRange, v.size)
  Base.axes(v::PtrView{ScalarT,N,LayoutLeft},d) where {ScalarT,N} = ZeroRange(v.size[d])
  Base.getindex(v::PtrView{ScalarT,1,LayoutLeft}, i::Integer) where {ScalarT} = unsafe_load(v.ptr,i+1)
  Base.setindex!(v::PtrView{ScalarT,1,LayoutLeft}, value, i::Integer) where {ScalarT} = unsafe_store!(v.ptr,value,i+1)
  Base.getindex(v::PtrView{ScalarT,N,LayoutLeft}, i::Integer) where {ScalarT,N} = unsafe_load(v.ptr,i)
  Base.setindex!(v::PtrView{ScalarT,N,LayoutLeft}, value, i::Integer) where {ScalarT,N} = unsafe_store!(v.ptr,value,i)
  Base.similar(v::PtrView{ScalarT,N,LayoutT}) where {ScalarT,N,LayoutT} = Array{ScalarT,N}(undef, length.(axes(v)))
  Base.size(v::PtrView) = v.size

  function Base.copy!(dest::AbstractArray{T,N}, src::PtrView{T,N,LayoutT}) where {T,N,LayoutT}
    println("copying")
    @show maximum(src)
    @boundscheck size.(axes(src)) == size.(axes(dest))
    for (I,J) in zip(eachindex(IndexStyle(dest), dest), eachindex(IndexStyle(src), src))
      @inbounds dest[I] = src[J]
    end
    return dest
  end

  struct StaticCrsGraph
    rowmap::Vector{Int}
    entries::Vector{Int}

    function StaticCrsGraph(nb_rows, nb_entries)
      rowmap = Vector{Int}(undef, nb_rows+1)
      entries = Vector{Int}(undef, nb_entries)
      rowmap[end] = nb_entries
      return new(rowmap,entries)
    end
  end

  get_datatype(::Type{View3{Ptr{DT},LayoutT,SpaceT}}) where {DT,LayoutT,SpaceT} = DT
  get_datatype(::Type{View3{CxxWrap.ConstPtr{DT},LayoutT,SpaceT}}) where {DT,LayoutT,SpaceT} = DT

  function to_cpp(g::StaticCrsGraph, ::Type{GraphT}) where {GraphT}
    entries_t = entries_type(GraphT)
    rowmap_t = row_map_type(GraphT)
    return GraphT(makeview(entries_t, "entries", convert(Vector{get_datatype(entries_t)}, g.entries)), makeview(rowmap_t, "rowmap", convert(Vector{get_datatype(rowmap_t)}, g.rowmap)))
  end
end


text:
using QNaNs
using Test

@test isnan(qnan(Int64(3))) == true
@test isnan(qnan(Int32(3))) == true
@test isnan(qnan(Int16(3))) == true

@test typeof(qnan(Int64(5))) == Float64
@test typeof(qnan(Int32(5))) == Float32
@test typeof(qnan(Int16(5))) == Float16

@test qnan(qnan(Int64(17))) == 17
@test qnan(qnan(Int32(17))) == 17
@test qnan(qnan(Int16(17))) == 17

@test qnan(qnan(Int64(-22))) == -22
@test qnan(qnan(Int32(-22))) == -22
@test qnan(qnan(Int16(-22))) == -22

@test qnan(zero(Int64)) === NaN
@test qnan(zero(Int32)) === NaN32


text:
include("../src/FluxFederated.jl")
using .FluxFederated
using Sockets
using Flux # Note: this takes awhile to load

# Connect to localhost port 3000
con = connect(3000)
println("Connected to localhost port 3000")

# Establish shared AES key via RSA
shared_key = handshake_client(con)

# Communicate securely
received_msg = secure_receive(con, shared_key) # --> "Wanna train some models?"
secure_send(con, "Oh yes", shared_key)

println("Received: $received_msg")

# Receive a Flux AI model via socket connection
model_bytes = secure_receive_bytes(con, shared_key)
model = model_from_bytes(model_bytes) # --> Chain(Dense(10,5,relu),Dense(5,2),softmax)

println("Received model: $model\n")

println("Model parameters:\n", params(model))

close(con)


text:
function interpp2(X,Y,V,Xq,Yq)
    itp=interpolate((X, Y), V', Gridded(Linear()))
    sz = size(Xq)
    y=zeros(sz[1],sz[2])
    for i=1:sz[1]
        for j=1:sz[2]
            y[i,j]=itp(Xq[i,j],Yq[i,j])
        end
    end
    return y
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel