text:
using LinearAlgebra, Random
using TensorRefinement.Auxiliary

import Base: length, deepcopy, reverse!, ndims, size, ones, zeros
import Base: fill!
import Random: rand!
import LinearAlgebra: rank
import Base: permutedims!
import Base: getindex
import Base: append!, prepend!, push!, pushfirst!, pop!, popfirst!, insert!, deleteat!
import Base: vcat, hcat
import LinearAlgebra: lmul!, rmul!
import Base: *, kron, +
import LinearAlgebra: qr!, svd!

export TT
export length, ndims, size, rank, ranks
export deepcopy, reverse!, permutedims!
export fill!, rand!
export getfirstfactor, getlastfactor, getfactor, setfirstfactor!, setlastfactor!, setfactor!
export rankselect!, rankselect
export getindex
export append!, prepend!, push!, pushfirst!, pop!, popfirst!, insert!, deleteat!
export compose!, compose, composecore!, composecore, composeblock!, composeblock
export vcat, hcat, dcat
export lmul!, rmul!, mul, had, *, kron, ⊗, add, +
export qr!, svd!

struct TT{T<:Number,N}
	factors::Dec{T,N}

	function TT(factors::Dec{T,N}) where {T<:Number,N}
		if N < 2
			throw(ArgumentError("each factor is required to have two rank indices"))
		end
		U = new{T,N}(factors)
		return U
	end

end

const Path = Union{Vector{Int},Int,UnitRange{Int},StepRange{Int,K} where K,NTuple{M,Int} where M,Vector{Any},Nothing,Colon}
const Permutation = Union{NTuple{K,Int},Vector{Int}} where K

function TT(::Type{T}, d::Int, L::Int) where {T<:Number}
	checkndims(d)
	checklength(L)
	factors = Vector{Array{T,d+2}}(undef, L)
	U = TT(factors)
	return U
end

TT(d::Int, L::Int) = TT(Float64, d, L)

function TT(::Type{T}, d::Int) where {T<:Number}
	checkndims(d)
	return TT(Vector{Array{T,d+2}}(undef, 0))
end

TT(d::Int) = TT(Float64, d)

function TT(U::TT{T,N}) where {T<:Number,N}
	return U
end

function TT(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number}
	checksize(n[:,:])
	nmat = isa(n, DecSize)
	nlen = size(n, 2)
	rvec = isa(r, DecRank)
	rlen = length(r)
	if any(r .≤ 0)
		throw(ArgumentError("the rank should be a positive integer or a vector of such"))
	end
	if rvec && (rlen < 2)
		throw(ArgumentError("when the rank parameter is specified as a vector, it should contain at least two entries"))
	end
	if rvec && nmat && rlen ≠ nlen+1
		throw(ArgumentError("when the size parameter is specified as a matrix with L columns and the rank parameter is specified as a vector, the latter should contain L+1 entries"))
	end
	if first < 0
		throw(ArgumentError("the first rank, when specified separately, should be positive"))
	end
	if last < 0
		throw(ArgumentError("the last rank, when specified separately, should be positive"))
	end
	if rvec && (first > 0 || last > 0)
		throw(ArgumentError("the first and last ranks are allowed to be specified separately only if the rank parameter is specified as an integer"))
	end
	if len < 0
		throw(ArgumentError("the number of factors, when specified, should be positive"))
	end
	if len > 0
		if nmat && len ≠ nlen
			throw(ArgumentError("when the number of rows in the size matrix is larger than one and the number of factors is specified, the two should be equal"))
		end
		if rvec && len+1 ≠ rlen
			throw(ArgumentError("when the rank parameter is specified as a vector with L+1 entries and the number of factors len is specified, it is required that len=L"))
		end
	end
	L = max(nlen, rlen-1, len)
	if !nmat
		n = repeat(n[:,:], outer=(1,L))
	end
	if !rvec
		r = r*ones(Int, L+1)
		if first > 0
			r[1] = first
		end
		if last > 0
			r[L+1] = last
		end
	end
	checkrank(r; len=L)
	factors = [ Array{T}(undef, r[ℓ], n[:,ℓ]..., r[ℓ+1]) for ℓ ∈ 1:L ]
	U = TT(factors)
	return U
end

TT(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = TT(Float64, n, r; first=first, last=last, len=len)

function TT!(V::Factor{T,N}) where {T<:Number,N}
	factors = [V]
	U = TT(factors)
	return U
end

function TT(V::Factor{T,N}; len::Int=1) where {T<:Number,N}
	if len < 0
		throw(ArgumentError("the number of factors, when specified, should be positive"))
	end
	if len == 0
		return TT(T, N-2)
	end
	if len > 1
		p,q = factorranks(V)
		if p ≠ q
			throw(ArgumentError("the two ranks of V should be equal when the number of factors is specified as larger than one"))
		end
	end
	factors = [ copy(V) for ℓ ∈ 1:len ]
	U = TT(factors)
	return U
end

function length(U::TT{T,N}) where {T<:Number,N}
	L = length(U.factors)
	checklength(L)
	return L
end

function ndims(U::TT{T,N}) where {T<:Number,N}
	d = [ factorndims(V) for V ∈ U.factors ]
	checkndims(d)
	return d[1]
end

function size(U::TT{T,N}) where {T<:Number,N}
	L = length(U)
	d = ndims(U)
	n = [ size(U.factors[ℓ], 1+k) for k ∈ 1:d, ℓ ∈ 1:L ]
	n = n[:,:]
	checksize(n)
	return n
end

function ranks(U::TT{T,N}) where {T<:Number,N}
	L = length(U)
	d = ndims(U)
	p = [ size(U.factors[ℓ], 1) for ℓ ∈ 1:L ]
	q = [ size(U.factors[ℓ], d+2) for ℓ ∈ 1:L ]
	return p,q
end

function rank(U::TT{T,N}) where {T<:Number,N}
	p,q = ranks(U)
	try checkranks(p,q) catch e
		isa(e, DimensionMismatch) && throw(DimensionMismatch("the factors have inconsistent ranks"))
	end
	return [p..., q[end]]
end

deepcopy(U::TT{T,N}) where {T<:Number,N} = TT(deepcopy(U.factors))

function reverse!(W::TT{T,N}) where {T<:Number,N}
	L = length(W)
	reverse!(W.factors)
	for ℓ ∈ 1:L
		W.factors[ℓ] = factorranktranspose(W.factors[ℓ])
	end
	return W
end

function permutedims!(U::TT{T,N}, τ::Permutation) where {T<:Number,N}
	d = N-2
	if d == 0
		throw(ArgumentError("the decomposition should have at least one mode dimension"))
	end
	if length(τ) ≠ d || !isperm(τ)
		throw(ArgumentError("τ is not a valid permutation of the mode dimensions of U"))
	end
	isa(τ, Vector{Int}) || (τ = collect(τ))
	L = length(U)
	for ℓ ∈ 1:L
		U.factors[ℓ] = factormodetranspose(U.fators[ℓ], τ)
	end
	return U
end

function fill!(U::TT{T,N}, v::T) where {T<:Number,N}
	L = length(U)
	for ℓ ∈ 1:L
		fill!(U.factors[ℓ], v)
	end
	return U
end

function rand!(rng::AbstractRNG, U::TT{T,N}) where {T<:Number,N}
	L = length(U)
	for ℓ ∈ 1:L
		rand!(rng, U.factors[ℓ])
	end
	return U
end
rand!(U::TT{T,N}) where {T<:Number,N} = rand!(Random.GLOBAL_RNG, U)

zeros(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = fill!(TT(T, n, r; first=first, last=last, len=len), convert(T, 0))
zeros(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = ttzeros(Float64, n, r; first=first, last=last, len=len)

ones(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = fill!(TT(T, n, r; first=first, last=last, len=len), convert(T, 1))
ones(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = ttzeros(Float64, n, r; first=first, last=last, len=len)

rand(rng::AbstractRNG, ::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = rand!(rng, TT(T, n, r; first=first, last=last, len=len))
rand(rng::AbstractRNG, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = ttrand(rng, Float64, n, r; first=first, last=last, len=len)
rand(::Type{T}, n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) where {T<:Number} = ttrand(Random.GLOBAL_RNG, T, n, r; first=first, last=last, len=len)
rand(n::Union{DecSize,FactorSize}, r::Union{Int,DecRank}; first::Int=0, last::Int=0, len::Int=0) = ttrand(Random.GLOBAL_RNG, Float64, n, r; first=first, last=last, len=len)

getfirstfactor(U::TT{T,N}) where {T<:Number,N} = U.factors[1]

getlastfactor(U::TT{T,N}) where {T<:Number,N} = U.factors[end]

function getfactor(U::TT{T,N}, ℓ::Int) where {T<:Number,N}
	L = length(U)
	if ℓ ⊈ 1:L
		throw(ArgumentError("ℓ is out of range"))
	end
	return U.factors[ℓ]
end

function setfirstfactor!(U::TT{T,N}, F::Array{T,N}) where {T<:Number,N}
	U.factors[1] = F
	return U
end

function setlastfactor!(U::TT{T,N}, F::Array{T,N}) where {T<:Number,N}
	U.factors[end] = F
	return U
end

function setfactor!(U::TT{T,N}, F::Array{T,N}, ℓ::Int) where {T<:Number,N}
	L = length(U)
	if ℓ ⊈ 1:L
		throw(ArgumentError("ℓ is out of range"))
	end
	U.factors[ℓ] = F
	return U
end

function rankselect!(U::TT{T,N}, α::Indices, β::Indices) where {T<:Number,N}
	# if isa(α, Int) || isa(β, Int)
	# 	throw(ArgumentError("for consistency with Base.selectdim, scalar α and β are not accepted; use α:α or β:β instead of α or β to select a subtensro of the factor whose first or second rank is one"))
	# end
	L = length(U)
	r = rank(U); p,q = r[1],r[L+1]
	α = indvec(α; min=1, max=p)
	β = indvec(β; min=1, max=q)
	if length(α) == 0
		throw(ArgumentError("the range for the first rank is empty"))
	end
	if α ⊈ 1:p
		throw(ArgumentError("the range for the first rank is incorrect"))
	end
	if length(β) == 0
		throw(ArgumentError("the range for the second rank is empty"))
	end
	if β ⊈ 1:q
		throw(ArgumentError("the range for the second rank is incorrect"))
	end
	setlastfactor!(U, factorrankselect(getlastfactor(U), :, β))
	setfirstfactor!(U, factorrankselect(getfirstfactor(U), α, :))
	return U
end

function rankselect(U::TT{T,N}, α::Indices, β::Indices) where {T<:Number,N}
	V = deepcopy(U)
	rankselect!(V, α, β)
	return V
end

getindex(U::TT{T,N}, α::Indices, β::Indices) where {T<:Number,N} = rankselect(U, α, β)

function append!(U::TT{T,N}, V::TT{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if ndims(U) ≠ ndims(V)
		throw(DimensionMismatch("U and V are inconsistent in the number of dimensions"))
	end
	p,q = ranks(U)
	r,s = ranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
		try checkranks(r,s) catch
			throw(ArgumentError("the factors of V have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		append!(p, r); append!(q, s);
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return append!(U.factors, V.factors)
end

function prepend!(U::TT{T,N}, V::TT{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if ndims(U) ≠ ndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = ranks(U); r,s = ranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
		try checkranks(r,s) catch
			throw(ArgumentError("the factors of V have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		prepend!(p, r); prepend!(q, s);
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return prepend!(U.factors, V.factors)
end

function push!(U::TT{T,N}, V::Factor{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if ndims(U) ≠ factorndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = ranks(U)
	r,s = factorranks(V)
	if rankprecheck
		try checkranks(p,q) catch
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); push!(p, r); push!(q, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return push!(U.factors, V)
end

function pushfirst!(U::TT{T,N}, V::Factor{T,N}; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	if ndims(U) ≠ factorndims(V)
		throw(DimensionMismatch("U and V have different numbers of dimensions"))
	end
	p,q = ranks(U)
	r,s = factorranks(V)
	if rankprecheck
		try
			checkranks(p,q)
		catch e
			throw(ArgumentError("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); pushfirst!(p, r); pushfirst!(q, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return pushfirst!(U.factors, V)
end

function pop!(U::TT{T,N}) where {T<:Number,N}
	return pop!(U.factors)
end

function popfirst!(U::TT{T,N}) where {T<:Number,N}
	return popfirst!(U.factors)
end

function insert!(U::TT{T,N}, ℓ::Int, V::Factor{T,N}; path::String="", rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	L = length(U)
	if path ∉ ("forward","backward")
		throw(ArgumentError("path should be either \"forward\" or \"backward\""))
	end
	if ℓ ∉ 1:L
		throw(ArgumentError("ℓ is required to be from 1:L, where L is the number of factors in U"))
	end
	(path == "forward") && (ℓ = ℓ+1)
	(ℓ == 1) && return pushfirst!(U, V; rankprecheck=rankprecheck, rankpostcheck=rankpostcheck)
	(ℓ == L+1) && return push!(U, V; rankprecheck=rankprecheck, rankpostcheck=rankpostcheck)
	p,q = ranks(U)
	if rankprecheck
		try checkranks(p,q) catch
			throw(DimensionMismatch("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		r,s = factorranks(V); insert!(p, ℓ, r); insert!(q, ℓ, s)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U and V are inconsistent for this operation"))
		end
	end
	return insert!(U.factors, ℓ, V)
end

function deleteat!(U::TT{T,N}, Λ::Path; rankprecheck::Bool=true, rankpostcheck::Bool=true) where {T<:Number,N}
	L = length(U)
	isa(Λ, Colon) && (Λ = collect(2:L); path = "backward")
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if unique(Λ) ≠ Λ
		throw(ArgumentError("the entries of Λ should be unique"))
	end
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ should be an element or a subset of 1:L with unique entries, where L is the number of factors in U"))
	end
	p,q = ranks(U)
	if rankprecheck
		try checkranks(p,q) catch
			throw(DimensionMismatch("the factors of U have incorrect or inconsistent ranks"))
		end
	end
	if rankpostcheck
		deleteat!(p, Λ); deleteat!(q, Λ)
		try checkranks(p,q) catch
			throw(DimensionMismatch("the ranks of U are inconsistent for this operation"))
		end
	end
	return deleteat!(U.factors, Λ)
end

function compose!(W::TT{T,N}, Λ::Path; path::String="", major::String="last") where {T<:Number,N}
	if path ∉ ("","forward","backward")
		throw(ArgumentError("the value of the keyword argument path should be \"\" (default, accepted only for empty Λ and for Λ=:), \"forward\" or \"backward\""))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	L = length(W)
	rank(W)
	isa(Λ, Colon) && (Λ = collect(2:L); path = "backward")
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if path ∉ ("forward","backward")
		throw(ArgumentError("when Λ is neither empty nor a colon, path should be either \"forward\" or \"backward\""))
	end
	(path == "backward") || (Λ .+= 1; path = "backward")
	if unique(Λ) ≠ Λ
		throw(ArgumentError("Λ has duplicate entries"))
	end
	if Λ ⊈ 2:L
		throw(ArgumentError("Λ, when nonempty, should be a colon or a Vector/NTuple/UnitRange of Int, from 1:L-1 for path=\"forward\" and from 2:L for path=\"backward\", where L is the number of factors in W"))
	end
	sort!(Λ; rev=true)
	for ℓ ∈ Λ
		F = getfactor(W, ℓ-1)
		F = factorcontract(F, getfactor(W, ℓ); major=major)
		setfactor!(W, F, ℓ-1)
		deleteat!(W, ℓ; rankprecheck=false, rankpostcheck=false)
	end
	return W
end

compose!(W::TT{T,N}; path::String="", major::String="last") where {T<:Number,N} = compose!(W, :; path=path, major=major)

function compose(W::TT{T,N}, Λ::Path; path::String="", major::String="last") where {T<:Number,N}
	U = deepcopy(W)
	return compose!(U, Λ; path=path, major=major)
end

compose(W::TT{T,N}; path::String="", major::String="last") where {T<:Number,N} = compose(W, :; path=path, major=major)

function composecore(W::TT{T,N}; major::String="last") where {T<:Number,N}
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end
	L = length(W); rank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	F = copy(getfactor(W, 1))
	for ℓ ∈ 2:L
		F = factorcontract(F, getfactor(W, ℓ); major=major)
	end
	return F
end

function composecore!(U::TT{T,N}; major::String="last") where {T<:Number,N}
	compose!(U; major=major)
	return getfirstfactor(U)
end

function composeblock!(U::TT{T,N}, α::Int, β::Int; major::String="last") where {T<:Number,N}
	r = rank(U); p,q = r[1],r[end]
	if α ∉ 1:p
		throw(ArgumentError("the first rank index is out of range"))
	end
	if β ∉ 1:q
		throw(ArgumentError("the second rank index is out of range"))
	end
	rankselect!(U, α:α, β:β)
	return compose(U; major=major)
end

function composeblock(U::Dec{T,N}, α::Int, β::Int; major::String="last") where {T<:Number,N}
	V = deepcopy(U)
	V = composeblock!(V, α, β, major=major)
	return V
end

function vcat(U::TT{T,N}, V::TT{T,N}, W::Vararg{TT{T,N},M}) where {T<:Number,N,M}
	L = length(U); d = ndims(U); m = size(U); p = rank(U)
	W = (V,W...)
	for V ∈ W
		if length(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if size(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
		q = rank(V)
		if q[L+1] ≠ p[L+1]
			throw(ArgumentError("the decompositions are incompatible in the last rank"))
		end
	end
	Z = TT(T, d, 0)
	for ℓ ∈ 1:L-1
		push!(Z, factordcat(getfactor(U, ℓ), [ getfactor(V, ℓ) for V ∈ W ]...))
	end
	push!(Z, factorvcat(getlastfactor(U), [ getlastfactor(V) for V ∈ W ]...))
	return Z
end

function hcat(U::TT{T,N}, V::TT{T,N}, W::Vararg{TT{T,N},M}) where {T<:Number,N,M}
	L = length(U); m = size(U); p = rank(U)
	W = (V,W...)
	for V ∈ W
		if length(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if size(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
		q = rank(V)
		if q[1] ≠ p[1]
			throw(ArgumentError("the decompositions are incompatible in the first rank"))
		end
	end
	Z = TT(T, d, 0)
	push!(Z, factorhcat(getfirstfactor(U), [ getfirstfactor(V) for V ∈ W ]...))
	for ℓ ∈ 2:L
		push!(Z, factordcat(getfactor(U, ℓ), [ getfactor(V, ℓ) for V ∈ W ]...))
	end
	return Z
end

function dcat(U::TT{T,N}, V::TT{T,N}, W::Vararg{TT{T,N},M}) where {T<:Number,N,M}
	L = length(U); m = decsize(U)
	W = (V,W...)
	for V ∈ W
		if length(V) ≠ L
			throw(ArgumentError("the decompositions are incompatible in the number of factors"))
		end
		if size(V) ≠ m
			throw(ArgumentError("the decompositions are incompatible in mode size"))
		end
	end
	Z = TT(T, d, 0)
	for ℓ ∈ 1:L
		push!(Z, factordcat(getfactor(U, ℓ), [ getfactor(V, ℓ) for V ∈ W ]...))
	end
	return Z
end

function lmul!(α::T, U::TT{T,N}) where {T<:Number,N}
	F = getlastfactor(U)
	F .= α * F
	return U
end

rmul!(U::TT{T,N}, α::T) where {T<:Number,N} = lmul!(α, U)

function *(α::T, U::TT{T,N}) where {T<:Number,N}
	V = deepcopy(U)
	return lmul!(α, U)
end

*(U::TT{T,N}, α::T) where {T<:Number,N} = *(α, U)

function mul(U₁::TT{T,N₁}, σ₁::Indices, U₂::TT{T,N₂}, σ₂::Indices) where {T<:Number,N₁,N₂}
	n₁ = size(U₁); d₁ = ndims(U₁)
	n₂ = size(U₂); d₂ = ndims(U₂)
	L = length(U₁)
	if length(U₂) ≠ L
		throw(ArgumentError("U₁ and U₂ differ in the number of factors"))
	end
	if isa(σ₁, Vector{Any}) && length(σ₁) > 0
		throw(ArgumentError("if σ₁ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	if isa(σ₂, Vector{Any}) && length(σ₂) > 0
		throw(ArgumentError("if σ₂ is passed as a vector, it should be a vector of the type Vector{Int} or an empty vector of the type Vector{Any}"))
	end
	isa(σ₁, Vector{Int}) || (σ₁ = indvec(σ₁; max=d₁))
	isa(σ₂, Vector{Int}) || (σ₂ = indvec(σ₂; max=d₂))
	if length(σ₁) ≠ length(σ₂)
		throw(ArgumentError("the specified sets of modes of σ₁ and σ₂ are inconsistent"))
	end
	if σ₁ ⊈ 1:d₁ || unique(σ₁) ≠ σ₁
		throw(ArgumentError("the set of modes of U₁ is specified incorrectly"))
	end
	if σ₂ ⊈ 1:d₂ || unique(σ₂) ≠ σ₂
		throw(ArgumentError("the set of modes of U₂ is specified incorrectly"))
	end
	if n₁[σ₁,:] ≠ n₂[σ₂,:]
		throw(ArgumentError("U₁ and U₂ are inconsistent with respect to the specified modes"))
	end
	factors = [ factormp(getfactor(U₁, ℓ), σ₁, getfactor(U₂, ℓ), σ₂) for ℓ ∈ 1:L ]
	return TT(factors)
end

function had(U::TT{T,N}, V::TT{T,N}) where {T<:Number,N}
	L = length(U)
	if length(V) ≠ L
		throw(ArgumentError("U and V differ in the number of factors"))
	end
	if size(U) ≠ size(V)
		throw(ArgumentError("U and V are inconsistent in mode size"))
	end
	factors = [ factorhp(getfactor(U, ℓ), getfactor(V, ℓ)) for ℓ ∈ 1:L ]
	return TT(factors)
end

*(U::TT{T,N}, V::TT{T,N}) where {T<:Number,N} = had(U, V)

function kron(U₁::Dec{T,N}, U₂::Dec{T,N}) where {T<:Number,N}
	L = length(U₁)
	if length(U₂) ≠ L
		throw(ArgumentError("U₁ and U₂ differ in the number of factors"))
	end
	if ndims(U₁) ≠ ndims(U₂)
		throw(ArgumentError("U₁ and U₂ are inconsistent in the number of mode dimensions"))
	end
	factors = [ factorkp(U₁[ℓ], U₂[ℓ]) for ℓ ∈ 1:L ]
	return TT(factors)
end

⊗(U₁::Dec{T,N}, U₂::Dec{T,N}) where {T<:Number,N} = kron(U₁, U₂)

function add(U::TT{T,N}, V::TT{T,N}) where {T<:Number,N}
	m = size(U); L = length(U)
	if length(V) ≠ L
		throw(ArgumentError("U and V differ in the number of factors"))
	end
	(L == 0) && return TT(T, d, 0)
	p = rank(U); q = rank(V)
	if size(V) ≠ m
		throw(ArgumentError("U and V are inconsistent in mode size"))
	end
	if q[1] ≠ p[1]
		throw(ArgumentError("the decompositions are incompatible in the first rank"))
	end
	if q[L+1] ≠ p[L+1]
		throw(ArgumentError("the decompositions are incompatible in the last rank"))
	end
	W = TT(T, d, 0)
	if L == 1
		push!(W, getfirstfactor(U) + getfirstfactor(V))
	elseif L > 1
		push!(W, factorhcat(getfirstfactor(U), getfirstfactor(V)))
		for ℓ ∈ 2:L-1
			push!(W, factordcat(getfactor(U, ℓ), getfactor(V, ℓ)))
		end
		push!(W, factorvcat(getlastfactor(U), getlastfactor(V)))
	end
	return W
end

+(U::TT{T,N}, V::TT{T,N}) where {T<:Number,N} = add(U, V)

function qr!(W::TT{T,N}, Λ::Path; path::String="") where {T<:FloatRC,N}
	L = length(W); rank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if path ∉ ("","forward","backward")
		throw(ArgumentError("path should be either \"\" (default, accepted only when path can be deduced from Λ), \"forward\" or \"backward\""))
	end
	if path == "" && (isa(Λ, Colon) || length(unique(Λ)) == 1)
		throw(ArgumentError("path cannot be deduced from Λ and should be specified as either \"forward\" or \"backward\""))
	end
	isa(Λ, Colon) && (Λ = (path == "forward") ? (1:L) : (L:-1:1))
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if unique(Λ) ≠ Λ
		throw(ArgumentError("Λ has duplicate entries"))
	end
	if (minimum(Λ):maximum(Λ)) ⊈ Λ
		throw(ArgumentError("the entries of Λ should form a set of contiguous integers"))
	end
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ is out of range"))
	end
	fw,bw = issorted(Λ, rev=false),issorted(Λ, rev=true)
	if !fw && !bw
		throw(ArgumentError("Λ is not sorted in ascending or descending order, so it is inconsistent with any forward or backward path"))
	end
	@assert !(fw && bw && path == "")
	fw && path == "" && (path = "forward")
	bw && path == "" && (path = "backward")
	if path == "forward" && !fw
		throw(ArgumentError("Λ is not sorted in ascending order, so it is inconsistent with any forward path"))
	end
	if path == "backward" && !bw
		throw(ArgumentError("Λ is not sorted in descending order, so it is inconsistent with any backward path"))
	end
	M = length(Λ)
	for λ ∈ 1:M
		ℓ = Λ[λ]
		Q,R = factorqr!(getfactor(W, ℓ); rev=(path == "backward"))
		setfactor!(W, Q, ℓ)
		if λ < M
			ν = Λ[λ+1]
			F = factorcontract(R, getfactor(W ,ν), rev=(path == "backward"))
			setfactor!(W, F, ν)
		else
			insert!(W, ℓ, R; path=path, rankprecheck=false, rankpostcheck=true)
		end
	end
	return W
end

qr!(W::TT{T,N}; path::String="") where {T<:FloatRC,N} = qr!(W, :; path=path)

function svd!(W::TT{T,N}, Λ::Path, n::Union{Colon,DecSize}; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, maxrank::Int2=0, major::String="last") where {T<:FloatRC,N}
	# the decomposition is assumed to be orthogonal
	L = length(W); rank(W)
	if L == 0
		throw(ArgumentError("the decomposition is empty"))
	end
	if path ∉ ("","forward","backward")
		throw(ArgumentError("path should be either \"\" (default, accepted only when path can be deduced from Λ), \"forward\" or \"backward\""))
	end
	if path == "" && (isa(Λ, Colon) || length(unique(Λ)) == 1)
		throw(ArgumentError("path cannot be deduced from Λ and should be specified as either \"forward\" or \"backward\""))
	end
	isa(Λ, Colon) && (Λ = (path == "forward") ? (1:L) : (L:-1:1))
	isa(Λ, Vector{Int}) || (Λ = indvec(Λ))
	(length(Λ) == 0) && return W
	if (minimum(Λ):maximum(Λ)) ⊈ Λ
		throw(ArgumentError("the entries of Λ should form a set of contiguous integers"))
	end
	if Λ ⊈ 1:L
		throw(ArgumentError("Λ is out of range"))
	end
	fw,bw = issorted(Λ, rev=false),issorted(Λ, rev=true)
	if !fw && !bw
		throw(ArgumentError("Λ is not sorted in ascending or descending order, so it is inconsistent with any forward or backward path"))
	end
	fw && path == "" && (path = "forward")
	bw && path == "" && (path = "backward")
	if path == "forward" && !fw
		throw(ArgumentError("Λ is not sorted in ascending order, so it is inconsistent with any forward path"))
	end
	if path == "backward" && !bw
		throw(ArgumentError("Λ is not sorted in descending order, so it is inconsistent with any backward path"))
	end
	m = size(W)
	isa(n, Colon) && (n = m[:,Λ])
	if size(n,1) ≠ size(m,1)
		throw(ArgumentError("the number of rows in n should be equal to the number of dimensions in each factor of W"))
	end
	if size(n,2) ≠ length(Λ)
		throw(ArgumentError("the number of columns in n should be equal to the number of elements in Λ"))
	end
	for ℓ ∈ unique(Λ)
		λ = (Λ .== ℓ)
		if m[:,ℓ] ≠ prod(n[:,λ]; dims=2)[:]
			throw(ArgumentError("n and Λ are incompatible with the size of the factors of W"))
		end
	end
	K = length(Λ)

	if any(aTol .< 0)
		throw(ArgumentError("aTol should be a nonnegative Float64 or a vector of such"))
	end
	if isa(aTol, Float64)
		aTol = aTol*ones(K)/sqrt(K)
	elseif length(aTol) ≠ K
		throw(ArgumentError("aTol, passed as a vector, has incorrect length"))
	end

	if any(aTolDistr .< 0)
		throw(ArgumentError("aTolDistr should be a nonnegative Float64 or a vector of such"))
	end
	if isa(aTolDistr, Float64)
		aTolDistr = aTolDistr*ones(K)/sqrt(K)
	elseif length(aTolDistr) ≠ K
		throw(ArgumentError("aTolDistr, passed as a vector, has incorrect length"))
	end

	if any(rTol .< 0)
		throw(ArgumentError("rTol should be a nonnegative Float64 or a vector of such"))
	end
	if isa(rTol, Float64)
		rTol = rTol*ones(K)/sqrt(K)
	elseif length(rTol) ≠ K
		throw(ArgumentError("rTol, passed as a vector, has incorrect length"))
	end

	if any(rTolDistr .< 0)
		throw(ArgumentError("rTolDistr should be a nonnegative Float64 or a vector of such"))
	end
	if isa(rTolDistr, Float64)
		rTolDistr = rTolDistr*ones(K)/sqrt(K)
	elseif length(rTolDistr) ≠ K
		throw(ArgumentError("rTolDistr, passed as a vector, has incorrect length"))
	end

	if any(maxrank .< 0)
		throw(ArgumentError("maxrank should be a nonnegative Int or a vector of such"))
	end
	if isa(maxrank, Int)
		maxrank = maxrank*ones(Int, K)
	elseif length(maxrank) ≠ K
		throw(ArgumentError("maxrank, passed as a vector, has incorrect length"))
	end
	if major ∉ ("first","last")
		throw(ArgumentError("major should be either \"last\" (default) or \"first\""))
	end

	ε = zeros(Float64, K); δ = zeros(Float64, K)
	σ = Vector{Vector{Float64}}(undef, K)
	ρ = zeros(Int, K)
	μ = 0.0
	aTolAcc = 0.0; rTolAcc = 0.0
	for λ ∈ 1:K
		F = getfactor(W, Λ[λ])
		if λ == 1
			ε₁ = [aTol[λ],aTolDistr[λ]]; ε₁ = ε₁[ε₁ .> 0]; ε₁ = isempty(ε₁) ? 0.0 : minimum(ε₁)
			δ₁ = [rTol[λ],rTolDistr[λ]]; δ₁ = δ₁[δ₁ .> 0]; δ₁ = isempty(δ₁) ? 0.0 : minimum(δ₁)
			U,V,ε[1],δ[1],μ,ρ[1],σ[1] = factorsvd!(F, n[:,λ], :; aTol=ε₁, rTol=δ₁, rank=maxrank[λ], rev=(path == "backward"), major=major)
		else
			(aTolDistr[λ] > 0) && (aTolDistr[λ] = sqrt(aTolDistr[λ]^2+aTolAcc^2); aTolAcc = 0.0)
			(rTolDistr[λ] > 0) && (rTolDistr[λ] = sqrt(rTolDistr[λ]^2+rTolAcc^2); rTolAcc = 0.0)
			ε₁ = [aTol[λ],aTolDistr[λ],μ*rTol[λ],μ*rTolDistr[λ]]; ε₁ = ε₁[ε₁ .> 0]
			ε₁ = isempty(ε₁) ? 0.0 : minimum(ε₁)
			U,V,ε[λ],_,_,ρ[λ],σ[λ] = factorsvd!(F, n[:,λ], :; aTol=ε₁, rank=maxrank[λ], rev=(path == "backward"), major=major)
			δ[λ] = (μ > 0) ? ε[λ]/μ : 0.0
		end
		setfactor!(W, U, Λ[λ])
		if λ < K && Λ[λ+1] ≠ Λ[λ]
			F = getfactor(W, Λ[λ+1])
			F = factorcontract(V, F, rev=(path == "backward"), major=major)
			setfactor!(W, F, Λ[λ+1])
		else
			insert!(W, Λ[λ], V; path=path, rankprecheck=false, rankpostcheck=true)
			(path == "forward") && (Λ[λ:end] .+= 1)
		end
	end
	return W,ε,δ,μ,ρ,σ
end

svd!(W::TT{T,N}, Λ::Path; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, maxrank::Int2=0, major::String="last") where {T<:FloatRC,N} = svd!(W, Λ, :; path=path, aTol=aTol, aTolDistr=aTolDistr, rTol=rTol, rTolDistr=rTolDistr, maxrank=maxrank, major=major)

svd!(W::TT{T,N}; path::String="", aTol::Float2=0.0, aTolDistr::Float2=0.0, rTol::Float2=0.0, rTolDistr::Float2=0.0, maxrank::Int2=0, major::String="last") where {T<:FloatRC,N} = svd!(W, :, :; path=path, aTol=aTol, aTolDistr=aTolDistr, rTol=rTol, rTolDistr=rTolDistr, maxrank=maxrank, major=major)


text:
function FieldStopOnEvent()
    return FieldStopOnEvent(())
end

function event_occurred(obj::FieldStopOnEvent, arg0::FieldSpacecraftState, arg1::FieldEventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (FieldSpacecraftState, FieldEventDetector, jboolean), arg0, arg1, arg2)
end

ext: jl


text:
using Test
using StaticArrays
using Random
using Nystrom
using LinearAlgebra
Random.seed!(1)

@testset "IFGF test" begin
    k    = 2π       # wavenumber
    λ    = 2π/k     # wavelength
    ppw  = 4        # points per wavelength
    dx   = λ/ppw    # distance between points
    nmax = 20       # max points per leaf box
    p    = (3,5,5)  # interpolation points per dimension
    pde  = Helmholtz(dim=3,k=k)

    # geometry
    Geometry.clear_entities!()
    geo = ParametricSurfaces.Sphere(;radius=1)
    Ω   = ParametricSurfaces.Domain(geo)
    Γ   = boundary(Ω)
    np  = ceil(Int,2/dx)
    M   = ParametricSurfaces.meshgen(Γ,(np,np))
    msh = NystromMesh(M,Γ;order=1)
    nx = length(msh.dofs)
    Xpts = msh.dofs
    B = rand(ComplexF64, nx)

    for op in [SingleLayerOperator, DoubleLayerOperator]
        iop = op(pde,msh)
        compress = IFGFCompressor(;p,nmax)
        A = compress(iop)
        exa = iop*B
        @test norm(A*B-exa)/norm(exa) < 1e-2
        # @info "" norm(A*B-exa)/norm(exa)
    end
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
using Random, Distributions, LinearAlgebra

struct LabeledPoint
    point
    label
    LabeledPoint(point, label) = new(point, label)
end

function parse_iris_row(line)
    """
    sepal_length, sepal_width, petal_length, petal_width, class
    """
    row = split(line, ",")
    measurements = [parse(Float64, value) for value in row[1:end-1]]
    # class is e.g. "Iris-virginica"; we just want "virginica"
    label = split(row[end], "-")[end]

    return LabeledPoint(measurements, label)
end

parse_iris_row("5.1,3.5,1.4,0.2,Iris-setosa")


text:
function show(io::IO, obj::ConvergenceFailed)
    println(io, "failed to find a solution: $(obj.reason)")
end

"""
    RobustLinResp

Robust linear response structure.

Solve the following minimization problem:

```math
\\min \\sum_i \\rho\\left(\\dfrac{r_i}{\\hat{\\sigma}}\right)
```

# Fields

- `est`: estimator used for the model
- `y`: response vector
- `μ`: mean response vector
- `offset`: offset added to `Xβ` to form `μ`. Can be of length 0
- `wts`: prior case weights.  Can be of length 0.
- `σ`: current estimate of the scale or dispersion
- `devresid`: the deviance residuals
- `wrkwt`: working case weights for the Iteratively Reweighted Least Squares (IRLS) algorithm
- `wrkres`: working residuals for IRLS
- `wrkscaledres`: scaled residuals for IRLS

"""
mutable struct RobustLinResp{T<:AbstractFloat,V<:AbstractVector{T},M<:AbstractEstimator} <:
               RobustResp{T}
    "`est`: estimator used for the model"
    est::M
    "`y`: response vector"
    y::V
    "`μ`: mean response vector"
    μ::V
    "`offset`: offset added to `Xβ` to form `μ`. Can be of length 0"
    offset::V
    "`wts`: prior case weights.  Can be of length 0."
    wts::V
    "`σ`: current estimate of the scale or dispersion"
    σ::T
    "`devresid`: the deviance residuals"
    devresid::V
    "`wrkwt`: working case weights for the Iteratively Reweighted Least Squares (IRLS) algorithm"
    wrkwt::V
    "`wrkres`: working residuals for IRLS"
    wrkres::V
    "`wrkscaledres`: scaled residuals for IRLS"
    wrkscaledres::V

    function RobustLinResp{T,V,M}(
        est::M,
        y::V,
        μ::V,
        off::V,
        wts::V,
        σ::Real=1,
    ) where {V<:AbstractVector{T},M<:AbstractEstimator} where {T<:AbstractFloat}
        n = length(y)
        length(μ) == n ||
            throw(DimensionMismatch("lengths of μ and y are not equal: $(length(μ)!=n)"))
        ll = length(off)
        ll in (0, n) || throw(DimensionMismatch("length of offset is $ll, must be $n or 0"))
        ll = length(wts)
        ll in (0, n) || throw(DimensionMismatch("length of wts is $ll, must be $n or 0"))
        σ > 0 || throw(DomainError(σ, "scale/dispersion must be positive"))
        new{T,V,M}(
            est,
            y,
            μ,
            off,
            wts,
            float(σ),
            similar(y),
            similar(y),
            similar(y),
            similar(y),
        )
    end
end

"""
    RobustLinResp(est::M, y::V, off::V, wts::V, σ::Real=1)
            where {V<:FPVector, M<:AbstractEstimator}

Create a response structure by copying the value of `y` to the mean `μ`.

"""
function RobustLinResp(
    est::M,
    y::V,
    off::V,
    wts::V,
    σ::Real=1,
) where {V<:FPVector,M<:AbstractEstimator}
    r = RobustLinResp{eltype(y),V,M}(est, y, similar(y), off, wts, float(σ))
    initresp!(r)
    return r
end

"""
    RobustLinResp(est::M, y, off, wts) where {M<:AbstractEstimator}

Convert the arguments to float arrays
"""
function RobustLinResp(est::M, y, off, wts) where {M<:AbstractEstimator}
    RobustLinResp(
        est,
        float(collect(y)),
        float(collect(off)),
        float(collect(wts)),
        float(1),
    )
end

function initresp!(r::RobustResp{T}) where {T}
    # Set μ to 0
    fill!(r.μ, 0)

    # Reset the factor of the TauEstimator
    if isa(r.est, TauEstimator)
        update_weight!(r.est, 0)
    end

    # Set residual (without offset)
    broadcast!(-, r.wrkres, r.y, r.μ)

    # Set working weights to the data weights
    if !isempty(r.wts)
        copyto!(r.wrkwt, r.wts)
    else
        fill!(r.wrkwt, 1)
    end
end

function Base.getproperty(r::RobustLinResp, s::Symbol)
    if s ∈ (:mu, :η, :eta)
        r.μ
    elseif s ∈ (:sigma, :scale)
        r.σ
    else
        getfield(r, s)
    end
end

function dispersion(
    r::RobustResp,
    dof_residual::Real=(nobs(r) - 1),
    sqr::Bool=false,
    robust::Bool=true,
)
    wrkwt, wrkres, wrkscaledres = r.wrkwt, r.wrkres, r.wrkscaledres
    if robust
        s = (r.σ)^2 *
            sum(i -> wrkwt[i] * abs2(wrkscaledres[i]), eachindex(wrkwt, wrkres)) /
            dof_residual
    else
        s = sum(i -> wrkwt[i] * abs2(wrkres[i]), eachindex(wrkwt, wrkres)) / dof_residual
    end

    return sqr ? s : sqrt(s)
end

"""
    location_variance(r::RobustLinResp, sqr::Bool=false)

Compute the part of the variance of the coefficients `β` that is due to the encertainty
from the location. If `sqr` is `false`, return the standard deviation instead.

From Maronna et al., Robust Statistics: Theory and Methods, Equation 4.49
"""
function location_variance(
    r::RobustLinResp,
    dof_residual::Real=(nobs(r) - 1),
    sqr::Bool=false,
)
    lpsi(x) = psi(r.est, x)
    lpsider(x) = psider(r.est, x)

    if isa(r.est, UnionL1)
        @warn "coefficient variance is not well defined for L1Estimator."
        return Inf
    end

    v = if isempty(r.wts)
        v = mean((lpsi.(r.wrkscaledres)) .^ 2)
        v /= (mean(lpsider.(r.wrkscaledres)))^2
    else
        wts = weights(r.wts)
        v = mean((lpsi.(r.wrkscaledres)) .^ 2, wts)
        v /= (mean(lpsider.(r.wrkscaledres), wts))^2
    end
    v *= r.σ^2
    v *= (nobs(r) / dof_residual)

    return sqr ? v : sqrt(v)
end

deviance(r::RobustResp) = sum(r.devresid)

function nulldeviance(r::RobustResp)
    ## TODO: take wts into account
    y, σ = r.y, r.σ
    μi = mean(y)

    dev = 0
    @inbounds for i in eachindex(y)
        dev += 2 * rho(r.est, (y[i] - μi) / σ)
    end
    dev
end

## TODO: define correctly the loglikelihood of the full model
fullloglikelihood(r::RobustLinResp) = -log(r.scale) - log(estimator_norm(r.est))

loglikelihood(r::RobustResp) = fullloglikelihood(r) - deviance(r) / 2

nullloglikelihood(r::RobustResp) = fullloglikelihood(r) - nulldeviance(r) / 2

response(r::RobustResp) = r.y

Estimator(r::RobustResp) = r.est

function weights(r::RobustResp{T}) where {T<:AbstractFloat}
    if isempty(r.wts)
        weights(ones(T, length(r.y)))
    else
        weights(r.wts)
    end
end

fitted(r::RobustLinResp) = r.μ

residuals(r::RobustLinResp) = r.wrkres

workingweights(r::RobustLinResp) = r.wrkwt

"""
    nobs(obj::RobustResp)
For linear and generalized linear models, returns the number of elements of the response
"""
function nobs(r::RobustResp{T}) where {T}
    if !isempty(r.wts)
        ## Suppose that the weights are probability weights
        count(!iszero, r.wts)
        # convert(T, sum(r.wts))
    else
        size(r.y, 1)
    end
end

function applyoffset!(r::RobustLinResp)
    # Subtract the offset because it will be added in the updateres! step
    if size(r.offset, 1) == size(r.y, 1)
        @inbounds @simd for i in eachindex(r.μ, r.offset)
            r.μ[i] -= r.offset[i]
        end
    end
end

"""
    setμ!{T<:FPVector}(r::RobustResp{T}, μ::T)
Update the mean to the given value `μ`.
If no value is provided, μ is set to the mean response `mean(r.y)`.
"""
function setμ! end

function setμ!(r::RobustLinResp, μ::T) where {T<:Real}
    fill!(r.μ, μ)

    # Subtract the offset because it will be added in the updateres! step
    applyoffset!(r)
end

function setμ!(r::RobustLinResp, μ::V) where {V<:FPVector}
    copyto!(r.μ, μ)

    # Subtract the offset because it will be added in the updateres! step
    applyoffset!(r)
end

function setμ!(r::RobustLinResp)
    if !isempty(r.wts)
        # Weighted mean
        m = mean(r.y, weights(r.wts))
    else
        m = mean(r.y)
    end
    setμ!(r, m)
end

function updateres!(r::RobustLinResp; updatescale=true, kwargs...)
    if updatescale
        _updateres_and_scale!(r; kwargs...)
    else
        _updateres!(r)
    end
end

"""
    _updateres!{T<:FPVector}(r::RobustResp{T})
Update the mean, working weights, working residuals and deviance residuals,
of the response `r`. The loss function for location estimation is used.
"""
function _updateres!(
    r::RobustLinResp{T,V,M},
) where {T<:AbstractFloat,V<:FPVector,M<:AbstractEstimator}
    ## Add offset to the linear predictor, if offset is defined
    ## and copy the linear predictor to the mean
    if length(r.offset) != 0
        broadcast!(+, r.μ, r.μ, r.offset)
    end

    y, μ, σ, devresid = r.y, r.μ, r.σ, r.devresid
    wrkres, wrkscaledres, wrkwt = r.wrkres, r.wrkscaledres, r.wrkwt

    invσ = inv(σ) # reduce #allocations using * instead of /.
    @inbounds for i in eachindex(y, μ, wrkres, wrkscaledres, wrkwt, devresid)
        wrkscaledres[i] = wrkres[i] = y[i] - μ[i]
        wrkscaledres[i] *= invσ
        wrkwt[i] = weight(r.est, wrkscaledres[i])
        devresid[i] = 2 * rho(r.est, wrkscaledres[i])
    end

    ## Multiply by the observation weights
    if !isempty(r.wts)
        broadcast!(*, r.devresid, r.devresid, r.wts)
        broadcast!(*, r.wrkwt, r.wrkwt, r.wts)
    end
end

function _updateres_and_scale!(
    r::RobustLinResp{T,V,M};
    kwargs...,
) where {T<:AbstractFloat,V<:FPVector,M<:Union{SEstimator,MMEstimator}}
    ## Add offset to the linear predictor, if offset is defined
    ## and copy the linear predictor to the mean
    if length(r.offset) != 0
        broadcast!(+, r.μ, r.μ, r.offset)
    end

    y, μ, devresid = r.y, r.μ, r.devresid
    wrkres, wrkscaledres, wrkwt = r.wrkres, r.wrkscaledres, r.wrkwt

    # First pass, compute the residuals
    @inbounds for i in eachindex(y, μ, wrkres)
        wrkres[i] = y[i] - μ[i]
    end

    # Second pass, compute the scale
    updatescale!(r; kwargs...)

    # Third pass, compute the scaled residuals, weights and deviance residuals
    invσ = inv(r.σ) # reduce #allocations using * instead of /.
    @inbounds for i in eachindex(wrkres, wrkscaledres, wrkwt, devresid)
        wrkscaledres[i] = wrkres[i] * invσ
        wrkwt[i] = weight(r.est, wrkscaledres[i])
        devresid[i] = 2 * rho(r.est, wrkscaledres[i])
    end

    ## Multiply by the observation weights
    if !isempty(r.wts)
        broadcast!(*, r.devresid, r.devresid, r.wts)
        broadcast!(*, r.wrkwt, r.wrkwt, r.wts)
    end
end

function _updateres_and_scale!(
    r::RobustLinResp{T,V,M};
    kwargs...,
) where {T<:AbstractFloat,V<:FPVector,M<:TauEstimator}
    ## Add offset to the linear predictor, if offset is defined
    ## and copy the linear predictor to the mean
    if length(r.offset) != 0
        broadcast!(+, r.μ, r.μ, r.offset)
    end

    y, μ, devresid = r.y, r.μ, r.devresid
    wrkres, wrkscaledres, wrkwt = r.wrkres, r.wrkscaledres, r.wrkwt

    # First pass, compute the residuals
    @inbounds for i in eachindex(y, μ, wrkres)
        wrkres[i] = y[i] - μ[i]
    end

    # Second pass, compute the scale
    updatescale!(r; kwargs..., sigma0=:mad)

    # Third pass, compute the scaled residuals
    invσ = inv(r.σ) # reduce #allocations using * instead of /.
    @inbounds for i in eachindex(wrkres, wrkscaledres)
        wrkscaledres[i] = wrkres[i] * invσ
    end

    # Fourth pass, compute the weights of the τ-estimator
    update_weight!(r.est, wrkscaledres; wts=r.wts)

    # Fifth pass, compute the weights and deviance residuals
    @inbounds for i in eachindex(wrkscaledres, wrkwt, devresid)
        wrkwt[i] = weight(r.est, wrkscaledres[i])
        devresid[i] = 2 * rho(r.est, wrkscaledres[i])
    end

    ## Multiply by the observation weights
    if !isempty(r.wts)
        broadcast!(*, r.devresid, r.devresid, r.wts)
        broadcast!(*, r.wrkwt, r.wrkwt, r.wts)
    end
end

"""
    updatescale!(
        r::RobustLinResp;
        sigma0::Union{Symbol,AbstractFloat}=r.σ,
        fallback::Union{Nothing, AbstractFloat}=nothing,
        verbose::Bool=false,
        kwargs...
    )

Compute the M-scale estimate from the estimator and residuals.
`sigma0` : initial value of the scale (default to r.σ),
           or `:mad` if it is estimated from mad(res).
`fallback` : if the algorithm does not converge, return the fallback value.
             If nothing, raise a ConvergenceFailed error.
`kwargs` : keyword arguments for the scale_estimate function.
"""
function updatescale!(
    r::RobustLinResp;
    sigma0::Union{Symbol,AbstractFloat}=r.σ,
    fallback::Union{Nothing,AbstractFloat}=nothing,
    verbose::Bool=false,
    kwargs...,
)
    est, res = r.est, r.wrkres
    σ0 = isa(sigma0, AbstractFloat) ? sigma0 : madresidualscale(r)

    if !isbounded(est)
        @warn "scale/dispersion is not changed because the estimator ($(est)) does not"*
              " allow scale estimation, a bounded estimator should be used, like Tukey."
        return r
    end

    verbose && print("Update scale: $(σ0) ")
    σ = try
        scale_estimate(est, res; σ0=σ0, wts=r.wts, kwargs...)
    catch e
        if isa(e, ConvergenceFailed) && !isnothing(fallback)
            verbose && print("  x>  fallback")
            σ = fallback
        else
            verbose && println("  x>  error")
            rethrow(e)
        end
    end
    verbose && println("  ->  $(σ)")
    r.σ = σ
    r
end

scale(r::RobustLinResp, sqr::Bool=false) =
    if sqr
        (r.σ)^2
    else
        r.σ
    end

function tauscale(
    r::RobustLinResp{T,V,M},
    sqr::Bool=false;
    verbose::Bool=false,
    bound::AbstractFloat=0.5,
    sigma0::Union{Symbol,AbstractFloat}=r.σ,
    updatescale::Bool=false,
) where {T,V,M<:TauEstimator}
    # No need to recompute s, as it was called in updateres! and stored in r.σ
    if updatescale
        updatescale!(r; sigma0=sigma0, verbose=verbose)
    end
    tau_scale_estimate(r.est, r.wrkres, r.σ, sqr; wts=r.wts, bound=bound)
end

"""
    madresidualscale(res)

Compute the scale using the MAD of the residuals.
"""
function madresidualscale(
    res::AbstractArray;
    factor::AbstractFloat=1.0,
    wts::AbstractArray=[],
)
    factor > 0 || error("factor should be positive")

    σ = if length(wts) == length(res)
        factor * mad(wts .* abs.(res); normalize=true)
    else
        factor * mad(abs.(res); normalize=true)
    end
    return σ
end
function madresidualscale(r::RobustLinResp; factor::AbstractFloat=1.0)
    madresidualscale(r.wrkres; factor=factor, wts=r.wts)
end

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


text:
function relabelBondsToString(
            lattice :: L
        ) :: Lattice{S,Bond{String,N},U} where {LS,LB,D,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U,L<:AbstractLattice{S,B,U}}

    # create a list of new bonds
    bondlist = Bond{String,N}[
        newBond(Bond{String,N}, from(b), to(b), string(label(b)), wrap(b))
        for b in bonds(lattice)
    ]

    # return a new lattice object
    return newLattice(
        Lattice{S,Bond{String,N},U},
        deepcopy(latticeVectors(lattice)),
        deepcopy(sites(lattice)),
        bondlist,
        deepcopy(unitcell(lattice))
    )
end

function relabelBondsToString(
            lattice :: L
        ) :: L where {LS,D,N,S<:AbstractSite{LS,D},U,L<:AbstractLattice{S,Bond{String,N},U}}
    # act as the identy as the bonds are already string labeled
    return lattice
end

# export the relabeling function
export relabelBondsToString


text:
function showpixel(a::Image, x::Int, y::Int)
    c = color(a, x, y)
    hex(c.r, 2)*hex(c.g, 2)*hex(c.b, 2)
end

w = 5
h = 7
a = Image(w, h)

purple = Color(0xff, 0, 0xff)
green = Color(0, 0xff, 0)

fill!(a, purple)
splat!(a, 2, 3, green)

for i in 1:h
    for j in 1:w
        print(showpixel(a, j, i), " ")
    end
    println()
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
################################################################################
#
#    NfOrd/NfOrdGenFracIdl.jl : Fractional ideals of generic
#                               orders in number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export frac_ideal

################################################################################
#
#  Construction
#
################################################################################

doc"""
***
    frac_ideal(O::NfOrd, A::FakeFmpqMat) -> NfOrdFracIdl

> Creates the fractional ideal of $\mathcal O$ with basis matrix $A$.
"""
function frac_ideal(O::NfOrdGen, x::FakeFmpqMat)
  y = hnf(x)
  z = NfOrdGenFracIdl(O, y)
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, A::fmpz_mat, b::fmpz) -> NfOrdFracIdl

> Creates the fractional ideal of $\mathcal O$ with basis matrix $A/b$.
"""
function frac_ideal(O::NfOrdGen, x::fmpz_mat, y::fmpz)
  y = FakeFmpqMat(x, y)
  z = NfOrdGenFracIdl(O, y)
  return z
end

frac_ideal(O::NfOrdGen, x::fmpz_mat, y::Integer) = frac_ideal(O, x, fmpz(y))

doc"""
***
    frac_ideal(O::NfOrd, I::NfOrdIdl) -> NfOrdFracIdl

> Turns the ideal $I$ into a fractional ideal of $\mathcal O$.
"""
function frac_ideal(O::NfOrdGen, x::NfOrdGenIdl)
  z = NfOrdGenFracIdl(O, x, fmpz(1))
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, I::NfOrdIdl, b::fmpz) -> NfOrdFracIdl

> Creates the fractional ideal $I/b$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrdGen, x::NfOrdGenIdl, y::fmpz)
  z = NfOrdGenFracIdl(O, x, deepcopy(y)) # deepcopy x?
  return z
end

frac_ideal(O::NfOrdGen, x::NfOrdGenIdl, y::Integer) = frac_ideal(O, x, fmpz(y))

doc"""
***
    frac_ideal(O::NfOrd, a::nf_elem) -> NfOrdFracIdl

> Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrdGen, x::nf_elem)
  z = NfOrdGenFracIdl(O, deepcopy(x))
  return z
end

doc"""
***
    frac_ideal(O::NfOrd, a::NfOrdElem) -> NfOrdFracIdl

> Creates the principal fractional ideal $(a)$ of $\mathcal O$.
"""
function frac_ideal(O::NfOrdGen, x::NfOrdElem{NfOrdGen})
  z = NfOrdGenFracIdl(O, elem_in_nf(x))
  return z
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::NfOrdGenFracIdl)
  print(io, "Fractional ideal of $(order(x)) with basis matrix\n")
  print(io, "$(basis_mat(x))\n")
end


text:
"""
```julia
    AbstractStopCondition
```

An abstract type that all stopping conditions for the iterative method
should inherit from.

Each type should define the () operator on itself, and return true from that
method when the stopping condition is met, and false if the condition is not
met.
"""
abstract type AbstractStopCondition end

struct StoppingConditionIterable{I, SC}
    iter::I
    condVec::Vector{SC}
end

function iterate( iter::StoppingConditionIterable )
    next = iterate( iter.iter )
    return dispatch( iter, next )
end

function iterate( iter::StoppingConditionIterable, (inst, next) )
    if inst == :halt
        return nothing
    end

    next = iterate( iter.iter, next )
    return dispatch( iter, next )
end

function dispatch( iter::StoppingConditionIterable, next )
    if next === nothing
        return nothing
    end

    nextInst = :continue

    for s in iter.condVec
        if s( next[2] )
            nextInst = :halt
        end
    end

    return (next[1], (nextInst, next[2]))
end

"""
```julia
    halt( iter::I, condVec::Vector{SC} ) where {I, SC}
```

Construct an iterator that will check the state of `iter` against the stopping conditions
`condVec`.
"""
halt( iter::I, condVec::Vector{SC} ) where {I, SC} = StoppingConditionIterable{I, eltype(condVec)}( iter, condVec )

"""
```julia
    halt( iter::I, cond::SC ) where {I, SC <: AbstractStopCondition}
```

Construct an iterator that will check the state of `iter` against the stopping condition
`condVec`.
"""
halt( iter::I, cond::SC ) where {I, SC <: AbstractStopCondition} = StoppingConditionIterable{I, SC}( iter, [cond] )


text:
function eom_lorenz_lorenz(u, p, t)
    c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
    x1, x2, x3, y1, y2, y3 = (u...,)

    dx1 = -a₁*(x1 - x2) + c_yx*(y1 - x1)
    dx2 = -x1*x3 + a₂*x1 - x2
    dx3 = x1*x2 - a₃*x3
    dy1 = -b₁*(y1 - y2) + c_xy*(x1 - y1)
    dy2 = -y1*y3 + b₂*y1 - y2
    dy3 = y1*y2 - b₃*y3

    return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
end

function lorenz_lorenz(; u0 = rand(6),
        c_xy = 0.2, c_yx = 0.0,
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 9/3)
    ContinuousDynamicalSystem(eom_lorenz_lorenz, u0, [c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃])
end

# Create an instance of the system with default parameters
sys = lorenz_lorenz()

# Analysis setup. We'll check for an influence from variable x1 to y1, which are the
# 1st and 4st components of the orbit. We'll generate orbits consisting of 500 points.
setup = ContinuousSystemSetup(source = 1, target = 4, n_pts = 500)

test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)

# Run test
@test causality(sys, setup, test_ccm) isa Vector{Vector{T}} where {T <: Real}
@test causality(sys, setup, test_cm) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_vf) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_tog) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_jdd) isa Vector{T} where {T <: Real}
@test causality(sys, setup, test_jddt) isa OneSampleTTest


text:
module Banjori

function map2lowcase(s)
    return 97 + (s-97)%26
end

function next_domain(seed)
    sl = Array{UInt8,1}(seed)
    sl[1] = map2lowcase(sl[1]+sl[4])
    sl[2] = map2lowcase(sl[1]+2*sl[2])
    sl[3] = map2lowcase(sl[1]+sl[3]-1)
    sl[4] = map2lowcase(sl[2]+sl[3]+sl[4])
    return join([Char(x) for x in sl])
end

function gen_one(seed="cckuailong")
    if length(seed) < 4
        print("Banjori seed length must >= 4")
        return missing
    else
        return next_domain(seed)
    end
end

function gen_many(cnt=100)
    res = []
    seeds = ["somestring", "firetruck", "bulldozer", "airplane", "racecar",
            "apartment", "laptop", "laptopcomp", "malwareisbad", "crazytrain",
            "thepolice", "fivemonkeys", "hockey", "football", "baseball",
            "basketball", "trackandfield", "fieldhockey", "softball", "redferrari",
            "blackcheverolet", "yellowelcamino", "blueporsche", "redfordf150",
            "purplebmw330i", "subarulegacy", "hondacivic", "toyotaprius",
            "sidewalk", "pavement", "stopsign", "trafficlight", "turnlane",
            "passinglane", "trafficjam", "airport", "runway", "baggageclaim",
            "passengerjet", "delta1008", "american765", "united8765", "southwest3456",
            "albuquerque", "sanfrancisco", "sandiego", "losangeles", "newyork",
            "atlanta", "portland", "seattle", "washingtondc"]
    seed_len = length(seeds)
    each_seed_items = trunc(Int64, cnt/seed_len)
    remain_items = cnt % seed_len
    for seed in seeds
        for i in 1:each_seed_items
            seed = next_domain(seed)
            res = [res..., seed]
        end
    end
    for seed in seeds[1:remain_items]
        seed = next_domain(seed)
        res = [res..., seed]
    end

    return res
end

end  # module Banjori


text:
A = rand(UInt8, 50, 20)
img_x = grayim(A)
img_y = permutedims(img_x, [2, 1])

@testset "crop" begin
    # test types stability
    @test typeof(Augmentor.crop(img_x, 1:50, 1:20)) <: typeof(img_x)
    @test typeof(Augmentor.crop(img_y, 1:50, 1:20)) <: typeof(img_y)
    @test_throws BoundsError Augmentor.crop(img_x, 0:50, 1:20)
    @test_throws BoundsError Augmentor.crop(img_x, 1:50, 0:20)
    @test_throws BoundsError Augmentor.crop(img_x, 1:51, 1:20)
    @test_throws BoundsError Augmentor.crop(img_x, 1:50, 1:21)

    # test values for some sample matrices
    @test Augmentor.crop(img_x, 20:40, 5:15) == img_x[20:40, 5:15]
    @test Augmentor.crop(img_y, 20:40, 5:15) == img_y[5:15, 20:40]
end


text:
import Pkg

"""
    is_loaded(mod::AbstractString)

Return `true` if `mod` is loaded.
"""
function is_loaded(mod::AbstractString)
    mods = Base.loaded_modules
    return any(m -> m.name == mod, keys(mods))
end

function module_version(mod::Module)
    toml = Base.parsed_toml(joinpath(dirname(dirname(Pkg.pathof(mod))), "Project.toml"))
    return toml["version"]
end

function pycall_version()
    if isdefined(@__MODULE__, :PyCall)
        return module_version(PyCall)
    end
    return nothing
end

function pythoncall_version()
    if isdefined(@__MODULE__, :PythonCall)
        return module_version(PythonCall)
    end
    return nothing
end

function parse_project()
    return parse_project(Pkg.project().path)
end

function parse_project(path::AbstractString)
    fullpath = nothing
    if ! endswith(path, "Project.toml")
        for f in ("Project.toml", "JuliaProject.toml")
            _fullpath = joinpath(path, f)
            if isfile(_fullpath)
                fullpath = _fullpath
                break
            end
            if isnothing(fullpath)
                error("Can't find project file in $path")
            end
        end
    else
        fullpath = path
    end
    return Base.parsed_toml(fullpath)
end

function  get_pycall_libpython()
    pycall_jl = Base.find_package("PyCall")
    if isnothing(pycall_jl)
        return (nothing, "not installed")
    end
    deps_jl = joinpath(dirname(dirname(Base.find_package("PyCall"))), "deps", "deps.jl")
    if ! isfile(deps_jl)
        return (nothing, "not built")
    end
    include(deps_jl) # not a great way to do this
    return (libpython, "ok")
end


text:
# shows.jl

using MIRT: @shows
using Test: @test

var = ones(3,4)
@test (@shows var) isa Nothing


text:
function removefreevariables!(sdp::SparseSDP)

end


text:
# For dispatching on the function name.
immutable Fn{name}
end

fnname(::Type{Fn}) = nothing
fnname{name}(::Type{Fn{name}}) = name

# Type constructor shortcut.
Fn(names...) =
    isempty(names) ? Union{} :
    length(names) == 1 ? Type{Fn{names[1]}} :
    Union{[Type{Fn{name}} for name in names]...}

reqcomposable(base::Scope, gs...) =
    for g in gs
        domain(base) <: idomain(g) || error("expected an expression defined on type $(domain(base)): $g")
    end

reqcomposable(f, gs...) =
    for g in gs
        odomain(f) <: idomain(g) || error("expected composable expressions: $f and $g")
    end

reqsingular(fs...) =
    for f in fs
        issingular(f) || error("expected a singular expression: $f")
    end

reqplural(fs...) =
    for f in fs
        !issingular(f) || error("expected a plural expression: $f")
    end

reqnonempty(fs...) =
    for f in fs
        isnonempty(f) || error("expected a complete expression: $f")
    end

reqpartial(fs...) =
    for f in fs
        !isnonempty(f) || error("expected a partial expression: $f")
    end

reqmonic(fs...) =
    for f in fs
        ismonic(f) || error("expected an exclusive expression: $f")
    end

reqodomain(dom, fs...) =
    for f in fs
        odomain(f) <: dom || error("expected an expression of type $dom: $f")
    end

reqidentity(fs...) =
    for f in fs
        !isnull(lookup(f, :__id)) || error("expected an expression with identity: $f")
    end

reqfields(fs...) =
    for f in fs
        !isnull(f.scope.items) && !isempty(get(f.scope.items)) || error("expected a composite expression: $f")
    end

reqtag(fs...) =
    for f in fs
        !isnull(f.scope.tag) || error("expected a tagged expression: $f")
    end

function >>(f::Query, g::Query)
    reqcomposable(f, g)
    pipe = f.pipe >> g.pipe
    syntax =
        !isnull(f.syntax) && !isnull(g.syntax) ?
            Nullable{AbstractSyntax}(ComposeSyntax(get(f.syntax), get(g.syntax))) :
        !isnull(f.syntax) ? f.syntax :
        !isnull(g.syntax) ? g.syntax :
        Nullable{AbstractSyntax}()
    return Query(g, pipe=pipe, syntax=syntax)
end

select(q::Query) =
    let out = lookup(q, (:__out, 0)),
        tag = q.scope.tag
        q = !isnull(out) ? q >> compile(get(out), q.scope, AbstractSyntax[]) : q
        if !isnull(tag)
            q = Query(q, scope=settag(q.scope, get(tag)))
        end
        q
    end

identify(q::Query) =
    let id = lookup(q, (:__id, 0))
        !isnull(id) ? q >> compile(get(id), q.scope, AbstractSyntax[]) : q
    end

record(base::Query, ops::Query...) =
    record(scope(base), ops...)

function record(base::Scope, ops::Query...)
    reqcomposable(base, ops...)
    pipe = TuplePipe([op.pipe for op in ops])
    scope = nest(base, odomain(pipe))
    fields = Query[]
    items = AbstractBinding[]
    hasid = false
    hasout = false
    for (k, op) in enumerate(ops)
        field = Query(op, pipe=ItemPipe(pipe, k))
        push!(fields, field)
        binding = SimpleBinding(field)
        push!(items, binding)
        if !isnull(field.scope.tag)
            scope = addlocal(scope, (get(field.scope.tag), 0), binding)
        end
        hasid = hasid || !isnull(lookup(op, :__id))
        hasout = hasout || !isnull(lookup(op, :__out))
    end
    scope = setitems(scope, tuple(items...))
    if hasid
        scope = addlocal(scope, (:__id, 0), record(scope, [identify(field) for field in fields]...))
    end
    if hasout
        scope = addlocal(scope, (:__out, 0), record(scope, [select(field) for field in fields]...))
    end
    return Query(scope, pipe)
end

compile(base::Query, syntax::AbstractSyntax) =
    compile(scope(base), syntax)

function compile(base::Scope, syntax::LiteralSyntax)
    T = domain(base)
    pipe = isa(syntax.val, Void) ? NullPipe(T) : ConstPipe(T, syntax.val)
    scope = nest(base, odomain(pipe))
    return Query(scope, pipe, syntax=syntax)
end

function compile(base::Scope, syntax::ApplySyntax)
    binding = lookup(base, (syntax.fn, length(syntax.args)))
    !isnull(binding) ?
        Query(compile(get(binding), base, syntax.args), syntax=syntax) :
        error("undefined combinator $(syntax.fn)")
end

function compile(base::Scope, syntax::ComposeSyntax)
    f = compile(base, syntax.f)
    g = compile(f, syntax.g)
    return Query(f >> g, syntax=syntax)
end

compile{name}(fn::Type{Fn{name}}, base::Scope, arg1::AbstractSyntax, args::AbstractSyntax...) =
    compile(fn, base, compile(base, arg1), [compile(base, arg) for arg in args]...)

function compile(::Fn(:here), base::Scope)
    T = domain(base)
    pipe = HerePipe(T)
    return Query(base, pipe)
end

function compile(::Fn(:unlink), base::Scope, arg::AbstractSyntax)
    scope = nest(base, Unit)
    I = domain(base)
    pipe = UnitPipe(I)
    root_base = Query(scope, pipe)
    return root_base >> compile(root_base, arg)
end

compile(::Fn(:link), base::Scope, pred::AbstractSyntax, arg::AbstractSyntax) =
    let unlink = compile(Fn{:unlink}, base, arg),
        mix = compile(Fn{:mix}, base, compile(Fn{:here}, base), unlink),
        condition = compile(mix, pred),
        filter = compile(Fn{:filter}, base, mix, condition),
        right = compile(Fn{:right}, mix.scope)
        filter >> right
    end

function compile(fn::Fn(:before, :and_before, :after, :and_after, :around, :and_around), base::Scope, ops::Query...)
    reqcomposable(base, ops...); reqsingular(ops...); reqnonempty(ops...)
    T = domain(base)
    before = fn in (Fn{:before}, Fn{:and_before}, Fn{:around}, Fn{:and_around})
    self = fn in (Fn{:and_before}, Fn{:and_after}, Fn{:and_around})
    after = fn in (Fn{:after}, Fn{:and_after}, Fn{:around}, Fn{:and_around})
    pipe =
        isempty(ops) ? RelativePipe(T, before, self, after) :
        length(ops) == 1 ? RelativeByPipe(identify(ops[1]).pipe, before, self, after) :
            RelativeByPipe(identify(record(base, ops...)).pipe, before, self, after)
    return Query(base, pipe)
end

compile(::Fn(:record), base::Scope, ops::Query...) = record(base, ops...)

compile(fn::Fn(:select), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow)
        compile(fn, base, flow, [compile(flow, op) for op in ops]...)
    end

function compile(::Fn(:select), base::Scope, flow::Query, ops::Query...)
    reqcomposable(base, flow)
    reqcomposable(flow, ops...)
    left = record(flow, ops...)
    right = Query(flow, pipe=HerePipe(odomain(flow)))
    pipe = TuplePipe(left.pipe, right.pipe)
    scope = nest(base, odomain(pipe))
    left_field = ItemPipe(pipe, 1)
    right_field = ItemPipe(pipe, 2)
    items = AbstractBinding[]
    defs = Dict{Symbol,Query}()
    for (k, op) in enumerate(ops)
        field = Query(op, pipe=left_field >> ItemPipe(left.pipe, k))
        push!(items, SimpleBinding(field))
        if !isnull(op.scope.tag)
            scope = addlocal(scope, get(op.scope.tag), field)
        end
    end
    if !isnull(flow.scope.tag)
        scope = addlocal(scope, get(flow.scope.tag), Query(flow, pipe=right_field))
    end
    maybe_id = lookup(flow, :__id)
    if !isnull(maybe_id)
        identity = compile(get(maybe_id), flow.scope)
        identity = Query(identity, pipe=right_field >> identity.pipe)
    else
        identity = Query(flow, pipe=right_field)
    end
    scope = addlocal(scope, :__id, identity)
    maybe_out = lookup(left, :__out)
    if !isnull(maybe_out)
        selector = compile(get(maybe_out), left.scope)
        selector = Query(selector, pipe=left_field >> selector.pipe)
    else
        selector = Query(left, pipe=left_field)
    end
    scope = addlocal(scope, :__out, selector)
    return flow >> Query(scope, pipe)
end

compile(fn::Fn(:filter), base::Scope, flow::AbstractSyntax, op::AbstractSyntax) =
    let flow = compile(base, flow)
        compile(fn, base, flow, compile(flow, op))
    end

function compile(::Fn(:filter), base::Scope, flow::Query, op::Query)
    reqcomposable(base, flow); reqplural(flow)
    reqcomposable(flow, op); reqsingular(op); reqodomain(Bool, op)
    O = odomain(flow)
    pipe = FilterPipe(flow.pipe, op.pipe)
    return Query(flow, pipe=pipe)
end

function compile(::Fn(:reverse), base::Scope, flow::Query)
    reqcomposable(base, flow); reqplural(flow)
    pipe = ReversePipe(flow.pipe)
    return Query(flow, pipe=pipe)
end

function compile(::Fn(:array), base::Scope, op1::Query, ops::Query...)
    reqcomposable(base, op1, ops...); reqsingular(op1, ops...); reqnonempty(op1, ops...)
    reqodomain(odomain(op1), ops...)
    scope = op1.scope
    pipe = VectorPipe([op.pipe for op in [op1, ops...]]...)
    return Query(scope, pipe)
end

function compile(::Fn(:range), base::Scope, start::Query, step::Query, stop::Query)
    reqcomposable(base, start, step, stop); reqsingular(start, step, stop); reqodomain(Int, start, step, stop)
    scope = nest(base, Int)
    if isnonempty(start) && isnonempty(step) && isnonempty(stop)
        pipe = RangePipe(start.pipe, step.pipe, stop.pipe)
    else
        pipe =
            (start.pipe * (step.pipe * stop.pipe)) >>
            RangePipe(
                ItemPipe(Tuple{Int,Tuple{Int,Int}}, 1),
                ItemPipe(Tuple{Int,Tuple{Int,Int}}, 2) >> ItemPipe(Tuple{Int,Int}, 1),
                ItemPipe(Tuple{Int,Tuple{Int,Int}}, 2) >> ItemPipe(Tuple{Int,Int}, 2))
    end
    return Query(scope, pipe)
end

function compile(fn::Fn(:first, :last), base::Scope, flow::Query)
    reqcomposable(base, flow); reqplural(flow)
    rev = flow.scope.rev
    rev = fn == Fn{:first} ? rev : !rev
    pipe = FirstPipe(flow.pipe, rev)
    return Query(flow, pipe=pipe)
end

compile(fn::Fn(:first, :last), base::Scope, flow::AbstractSyntax, op::AbstractSyntax) =
    let flow = compile(base, flow)
        compile(fn, base, flow, compile(flow, op))
    end

function compile(fn::Fn(:first, :last), base::Scope, flow::Query, op::Query)
    reqcomposable(base, flow); reqplural(flow)
    reqcomposable(flow, op); reqsingular(op); reqnonempty(op)
    rev = op.scope.rev
    rev = fn == Fn{:first} ? rev : !rev
    op = identify(op)
    pipe = FirstByPipe(flow.pipe, op.pipe, rev)
    return Query(flow, pipe=pipe)
end

function compile(fn::Fn(:take, :skip), base::Scope, flow::Query, size::Query)
    reqcomposable(base, flow); reqplural(flow)
    reqcomposable(base, size); reqsingular(size); reqnonempty(size); reqodomain(Int, size)
    pipe = TakePipe(flow.pipe, size.pipe, fn == Fn{:take} ? false : true)
    return Query(flow, pipe=pipe)
end

function compile(::Fn(:get), base::Scope, flow::Query, idx::Query)
    reqcomposable(base, flow); reqplural(flow); reqmonic(flow); reqidentity(flow)
    reqcomposable(base, idx); reqsingular(idx); reqnonempty(idx)
    maybe_id = lookup(flow, :__id)
    identity = compile(get(maybe_id), flow.scope)
    reqodomain(odomain(identity), idx)
    pipe = GetPipe(flow.pipe, identity.pipe, idx.pipe)
    return Query(flow, pipe=pipe)
end

function compile(::Fn(:asc), base::Scope, op::Query)
    return Query(op, scope=setrev(op.scope, false))
end

function compile(::Fn(:desc), base::Scope, op::Query)
    return Query(op, scope=setrev(op.scope, true))
end

compile(fn::Fn(:sort), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow)
        compile(fn, base, flow, [compile(flow, op) for op in ops]...)
    end

function compile(fn::Fn(:sort), base::Scope, flow::Query, ops::Query...)
    reqcomposable(base, flow); reqplural(flow)
    reqcomposable(flow, ops...); reqsingular(ops...); reqnonempty(ops...)
    if isempty(ops)
        maybe_id = lookup(flow, :__id)
        if !isnull(maybe_id)
            identity = Query(compile(get(maybe_id), flow.scope))
            identity = Query(identity, scope=setrev(identity.scope, flow.scope.rev))
            return compile(fn, base, flow, identity)
        end
        pipe = SortPipe(flow.pipe, flow.scope.rev)
    else
        pipe = SortByPipe(
            flow.pipe,
            [(identify(op).pipe, op.scope.rev) for op in ops]...)
    end
    return Query(flow, pipe=pipe)
end

function compile(fn::Fn(:connect, :and_connect), base::Scope, op::Query)
    reqcomposable(base, op); reqpartial(op); reqodomain(idomain(op), op)
    pipe = ConnectPipe(op.pipe, fn == Fn{:and_connect})
    return Query(op, pipe=pipe)
end

function compile(::Fn(:depth), base::Scope, op::Query)
    reqcomposable(base, op); reqpartial(op); reqodomain(idomain(op), op)
    scope = nest(base, Int)
    pipe = DepthPipe(op.pipe)
    return Query(scope, pipe)
end

compile(fn::Fn(:sort_connect), base::Scope, flow::AbstractSyntax, op::AbstractSyntax) =
    let flow = compile(base, flow)
        compile(fn, base, flow, compile(flow, op))
    end

function compile(::Fn(:sort_connect), base::Scope, flow::Query, op::Query)
    reqcomposable(base, flow); reqplural(flow)
    reqcomposable(flow, op); reqpartial(op); reqodomain(idomain(op), op)
    maybe_id = lookup(op, :__id)
    if isnull(maybe_id)
        pipe = SortConnectPipe(flow.pipe, op.pipe)
    else
        identity = compile(get(maybe_id), op.scope)
        pipe = SortConnectPipe(flow.pipe, identity.pipe, op.pipe)
    end
    return Query(flow, pipe=pipe)
end

function compile(::Fn(:unique), base::Scope, flow::Query)
    reqcomposable(base, flow); reqplural(flow)
    maybe_id = lookup(flow, :__id)
    if !isnull(maybe_id)
        identity = compile(get(maybe_id), flow.scope)
        PipeType = ismonic(flow) ? SortByPipe : UniquePipe
        pipe = PipeType(flow.pipe, identity.pipe, flow.scope.rev)
    else
        PipeType = ismonic(flow) ? SortPipe : UniquePipe
        pipe = PipeType(flow.pipe, flow.scope.rev)
    end
    return Query(flow, pipe=pipe)
end

compile(fn::Fn(:group, :group_cube),
        base::Scope, flow::AbstractSyntax, op1::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow)
        compile(fn, base, flow, compile(flow, op1), [compile(flow, op) for op in ops]...)
    end

function compile(
        fn::Fn(:group, :group_cube, :partition, :partition_cube),
        base::Scope, flow::Query, ops::Query...)
    iscube = (fn == Fn{:group_cube} || fn == Fn{:partition_cube})
    ispartition = (fn == Fn{:partition} || fn == Fn{:partition_cube})
    reqcomposable(base, flow); reqplural(flow)
    if !ispartition
        reqcomposable(flow, ops...); reqsingular(ops...); reqnonempty(ops...)
    else
        reqcomposable(base, ops...); reqplural(ops...); reqtag(ops...)
    end
    if !ispartition
        groups = []
        for op in ops
            maybe_id = lookup(op, :__id)
            if isnull(maybe_id)
                push!(groups, (op.pipe, op.scope.rev))
            else
                identity = compile(get(maybe_id), op.scope)
                push!(groups, (op.pipe, identity.pipe, op.scope.rev))
            end
        end
    else
        groups = []
        for op in ops
            img = lookup(flow, get(op.scope.tag))
            !isnull(img) || error("undefined attribute $(get(op.tag)): $flow")
            img = compile(get(img), base)
            reqsingular(img); reqnonempty(img); reqodomain(odomain(op), img)
            maybe_id = lookup(op, :__id)
            if isnull(maybe_id)
                push!(groups, (op.pipe, img.pipe, op.scope.rev))
            else
                identity = compile(get(maybe_id), op.scope)
                push!(groups, (op.pipe, img.pipe, identity.pipe, op.scope.rev))
            end
        end
    end
    pipe = (!ispartition ? GroupPipe : PartitionPipe)(flow.pipe, iscube, groups...)
    scope = nest(base, odomain(pipe))
    fields = Query[]
    items = AbstractBinding[]
    kernel_pipe = ItemPipe(pipe, 1)
    for (k, op) in enumerate(ops)
        item_pipe = ItemPipe(pipe, 1, k)
        field = Query(op, pipe=item_pipe)
        push!(fields, field)
        item = SimpleBinding(field)
        push!(items, item)
        if !isnull(field.scope.tag)
            scope = addlocal(scope, get(field.scope.tag), item)
        end
    end
    kernel_field = record(scope, fields...)
    scope = addlocal(scope, :__id, identify(kernel_field))
    flow_field = Query(flow, pipe=ItemPipe(pipe, 2))
    push!(fields, flow_field)
    flow_item = SimpleBinding(flow_field)
    push!(items, flow_item)
    if !isnull(flow_field.scope.tag)
        scope = addlocal(scope, get(flow_field.scope.tag), flow_item)
    end
    scope = setitems(scope, tuple(items...))
    scope = addlocal(scope, :__out, select(record(scope, fields...)))
    return Query(scope, pipe)
end

function compile(::Fn(:mix), base::Scope, ops::Query...)
    reqcomposable(base, ops...); reqtag(ops...)
    if isempty(ops)
        scope = next(base, Tuple{})
        I = odomain(base)
        pipe = ConstPipe(I, ())
        return Query(scope, pipe)
    else
        pipe = MixPipe([op.pipe for op in ops])
        T = odomain(pipe)
        scope = nest(base, T)
        fields = Query[]
        items = AbstractBinding[]
        for (k, op) in enumerate(ops)
            field = Query(op, pipe=ItemPipe(T, k))
            push!(fields, field)
            item = SimpleBinding(field)
            push!(items, item)
            scope = addlocal(scope, get(op.scope.tag), item)
        end
        scope = setitems(scope, tuple(items...))
        identity = record(scope, [identify(field) for field in fields]...)
        selector = record(scope, [select(field) for field in fields]...)
        scope = addlocal(scope, :__id, identity)
        scope = addlocal(scope, :__out, selector)
        return Query(scope, pipe)
    end
end

function compile(fn::Fn(:left, :right), base::Scope)
    reqfields(Query(base, HerePipe(base.domain)))
    return compile(get(base.items)[fn == Fn{:left} ? 1 : end], base)
end

function compile(::Fn(:pack), base::Scope, ops::Query...)
    reqcomposable(base, ops...); reqtag(ops...)
    I = domain(base)
    Fs = Vector{Pair{Symbol,AbstractPipe}}()
    isisnonempty = false
    identitymap = Dict{Symbol,AbstractPipe}()
    selectormap = Dict{Symbol,AbstractPipe}()
    for op in ops
        T = odomain(op)
        tag = get(op.scope.tag)
        push!(Fs, Pair{Symbol,AbstractPipe}(tag, op.pipe))
        if isnonempty(op)
            isisnonempty = true
        end
        maybe_id = lookup(op, :__id)
        identitymap[tag] =
            !isnull(maybe_id) ? compile(get(maybe_id), op.scope).pipe : HerePipe(T)
        maybe_out = lookup(op, :__out)
        selectormap[tag] =
            !isnull(maybe_out) ? compile(get(maybe_out), op.scope).pipe : HerePipe(T)
    end
    pipe = PackPipe(Fs...)
    scope = nest(base, odomain(pipe))
    U = odomain(pipe).parameters[2]
    identity = Query(scope, CasePipe(U, true, identitymap...))
    selector = Query(scope, CasePipe(U, true, selectormap...))
    scope = addlocal(scope, :__id, identity)
    scope = addlocal(scope, :__out, selector)
    for op in ops
        tag = get(op.scope.tag)
        T = odomain(op)
        scope = addlocal(scope, tag, Query(op, pipe=CasePipe(U, false, tag => HerePipe(T))))
    end
    return Query(scope, pipe)
end

function compile(::Fn(:as), base::Scope, op::AbstractSyntax, ident::AbstractSyntax)
    (isa(ident, ApplySyntax) && isempty(ident.args)) || error("expected an identifier: $ident")
    q = compile(base, op)
    scope = settag(q.scope, ident.fn)
    return Query(q, scope=scope)
end

compile(::Fn(:(=>)), base::Scope, ident::AbstractSyntax, op::AbstractSyntax) =
    compile(Fn{:as}, base, op, ident)

compile(fn::Fn(:define), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow)
        for op in ops
            flow = compile(fn, base, flow, compile(flow, op))
        end
        flow
    end

function compile(::Fn(:define), base::Scope, flow::Query, op::Query)
    reqcomposable(base, flow)
    reqcomposable(flow, op); reqtag(op)
    return Query(flow, scope=addlocal(flow.scope, get(op.scope.tag), op))
end

function compile(::Fn(:!), base::Scope, op::Query)
    reqcomposable(base, op); reqsingular(op); reqodomain(Bool, op)
    scope = nest(base, Bool)
    pipe = ~op.pipe
    return Query(scope, pipe)
end

function compile(fn::Fn(:&, :|), base::Scope, op1::Query, op2::Query)
    reqcomposable(base, op1, op2); reqsingular(op1, op2); reqodomain(Bool, op1, op2)
    scope = nest(base, Bool)
    if fn == Fn{:&}
        pipe = op1.pipe & op2.pipe
    else
        pipe = op1.pipe | op2.pipe
    end
    return Query(scope, pipe)
end

function compile(fn::Fn(:(==), :(!=)), base::Scope, op1::Query, op2::Query)
    reqcomposable(base, op1, op2);
    typejoin(odomain(op1), odomain(op2)) != Any || reqodomain(odomain(op1), op2)
    ismonic(op1) || reqsingular(op2)
    T = odomain(op1)
    TT = Tuple{T,T}
    scope = nest(base, Bool)
    output = Output(Bool, max(omode(op1), omode(op2)))
    pipe = fn == Fn{:(==)} ? op1.pipe .== op2.pipe : op1.pipe .!= op2.pipe
    return Query(scope, pipe)
end

function compile(fn::Fn(:in), base::Scope, op1::Query, op2::Query)
    reqcomposable(base, op1, op2);
    typejoin(odomain(op1), odomain(op2)) != Any || reqodomain(odomain(op1), op2)
    reqsingular(op1); reqnonempty(op1); reqplural(op2)
    scope = nest(base, Bool)
    pipe = InPipe(op1.pipe, op2.pipe)
    return Query(scope, pipe)
end

const SCALAR_OPS = (
    :(+), :(-), :(*), :(/),
    :(<), :(<=), :(>), :(>=),
    :startswith, :endswith, :contains, :substr,
    :date,
    :day, :_day, :days, :_days,
    :month, :_month, :months, :_months,
    :year, :_year, :years, :_years,
    :_usd)

const AGGREGATE_OPS = (
    :sum, :max, :min, :mean,
    :count, :exists,
    :all, :any)

function compile(fn::Fn(SCALAR_OPS...), base::Scope, ops::Query...)
    reqcomposable(base, ops...)
    for j = eachindex(ops)
        issingular(ops[j]) || reqsingular(ops[1:j-1]...)
    end
    Ts = Type[odomain(op) for op in ops]
    method_exists(polydomain, Tuple{Type{fn}, [Type{T} for T in Ts]...}) ||
        error("expected compatible operands for operator $(fnname(fn)): $(join(ops, ", "))")
    method = polymethod(fn, Ts...)
    domain = polydomain(fn, Ts...)
    scope = nest(base, domain)
    pipe = OpPipe(method, Ts, domain, [op.pipe for op in ops])
    return Query(scope, pipe)
end

function compile(fn::Fn(AGGREGATE_OPS...), base::Scope, flow::Query)
    reqcomposable(base, flow); reqplural(flow)
    T = odomain(flow)
    method = polymethod(fn, T)
    domain = polydomain(fn, T)
    haszero = polyhaszero(fn, T)
    scope = nest(base, domain)
    pipe = AggregateOpPipe(method, T, domain, haszero, flow.pipe)
    scope = nest(base, odomain(pipe))
    return Query(scope, pipe)
end

polymethod{name}(::Type{Fn{name}}, ::Type...) = eval(name)
polyhaszero{name}(::Type{Fn{name}}, ::Type) = true

polydomain(::Fn(:(+), :(-)), ::Type{Int}) = Int
polydomain(::Fn(:(+), :(-)), ::Type{Int}, ::Type{Int}) = Int
polydomain(::Fn(:(+), :(-)), ::Type{Int}, ::Type{Float64}) = Float64
polydomain(::Fn(:(+), :(-)), ::Type{Float64}, ::Type{Int}) = Float64
polydomain(::Fn(:(+), :(-)), ::Type{Float64}, ::Type{Float64}) = Float64
polydomain{M<:Monetary}(::Fn(:(+), :(-)), ::Type{M}, ::Type{M}) = M
polydomain(::Fn(:(*)), ::Type{Int}, ::Type{Int}) = Int
polymethod(::Fn(:(/)), ::Type{Int}, ::Type{Int}) = div
polydomain(::Fn(:(/)), ::Type{Int}, ::Type{Int}) = Int
polydomain(::Fn(:(*), :(/)), ::Type{Int}, ::Type{Float64}) = Float64
polydomain(::Fn(:(*), :(/)), ::Type{Float64}, ::Type{Int}) = Float64
polydomain(::Fn(:(*), :(/)), ::Type{Float64}, ::Type{Float64}) = Float64
polydomain{T<:Number,M<:Monetary}(::Fn(:(*)), ::Type{T}, ::Type{M}) = M
polydomain{T<:Number,M<:Monetary}(::Fn(:(*)), ::Type{M}, ::Type{T}) = M
polydomain{T<:Number,M<:Monetary}(::Fn(:(*)), ::Type{M}, ::Type{T}, ::Type{T}) = M
polydomain{M<:Monetary}(::Fn(:(/)), ::Type{M}, ::Type{M}) = Float64
polydomain{T<:Number,M<:Monetary}(::Fn(:(/)), ::Type{M}, ::Type{T}) = M
polydomain{T1<:Number,T2<:Number}(::Fn(:(<), :(<=), :(>), :(>=)), ::Type{T1}, ::Type{T2}) = Bool
polydomain{M<:Monetary}(::Fn(:(<), :(<=), :(>), :(>=)), ::Type{M}, ::Type{M}) = Bool

_regex_contains(s, r) = ismatch(r, s)
_substr(s, i, j) = s[i:j]
polydomain{S1<:AbstractString,S2<:AbstractString}(
    ::Fn(:startswith, :endswith, :contains), ::Type{S1}, ::Type{S2}) = Bool
polymethod{S<:AbstractString}(::Fn(:contains), ::Type{S}, ::Type{Regex}) = _regex_contains
polydomain{S<:AbstractString}(::Fn(:contains), ::Type{S}, ::Type{Regex}) = Bool
polymethod{S<:AbstractString}(::Fn(:substr), ::Type{S}, ::Type{Int}, ::Type{Int}) = _substr
polydomain{S<:AbstractString}(::Fn(:substr), ::Type{S}, ::Type{Int}, ::Type{Int}) = S

_dates_date(s) = Date(s)
_dates_day() = Dates.Day(1)
_dates_month() = Dates.Month(1)
_dates_year() = Dates.Year(1)
polymethod(::Fn(:year), ::Type{Date}) = Dates.year
polymethod(::Fn(:month), ::Type{Date}) = Dates.month
polymethod(::Fn(:day), ::Type{Date}) = Dates.day
polydomain(::Fn(:year, :month, :day), ::Type{Date}) = Int
polymethod{S<:AbstractString}(::Fn(:date), ::Type{S}) = _dates_date
polydomain{S<:AbstractString}(::Fn(:date), ::Type{S}) = Date
polymethod(::Fn(:day, :_day, :days, :_days)) = _dates_day
polymethod(::Fn(:month, :_month, :months, :_months)) = _dates_month
polymethod(::Fn(:year, :_year, :years, :_years)) = _dates_year
polydomain(::Fn(:day, :_day, :days, :_days)) = Dates.Day
polydomain(::Fn(:month, :_month, :months, :_months)) = Dates.Month
polydomain(::Fn(:year, :_year, :years, :_years)) = Dates.Year
polydomain{P<:Dates.Period}(::Fn(:(+), :(-)), ::Type{P}) = Dates.Period
polydomain{P1<:Dates.Period,P2<:Dates.Period}(::Fn(:(+), :(-)), ::Type{P1}, ::Type{P2}) = Dates.Period
polydomain{P<:Dates.Period}(::Fn(:(+), :(-)), ::Type{Date}, ::Type{P}) = Date
polydomain{P<:Dates.Period}(::Fn(:(*)), ::Type{Int}, ::Type{P}) = Dates.Period
polydomain(::Fn(:(<), :(<=), :(>), :(>=)), ::Type{Date}, ::Type{Date}) = Bool

_usd() = Monetary{:USD}(1, 0)
polymethod(::Fn(:_usd)) = _usd
polydomain(::Fn(:_usd)) = Monetary{:USD}
polydomain{T<:Number}(::Fn(:sum), ::Type{T}) = T
polydomain{M<:Monetary}(::Fn(:sum), ::Type{M}) = M
polymethod(::Fn(:max), ::Type) = maximum
polymethod(::Fn(:min), ::Type) = minimum
polydomain{T}(::Fn(:max, :min), ::Type{T}) = T
polyhaszero(::Fn(:max, :min, :mean), ::Type) = false
polydomain(::Fn(:mean), ::Type{Int}) = Float64
polydomain(::Fn(:mean), ::Type{Float64}) = Float64
polydomain{M<:Monetary}(::Fn(:mean), ::Type{M}) = M

_exists(a) = !isempty(a)
polydomain(::Fn(:all, :any), ::Type{Bool}) = Bool
polymethod(::Fn(:exists), ::Type) = _exists
polydomain(::Fn(:exists), ::Type) = Bool
polymethod(::Fn(:count), ::Type) = length
polydomain(::Fn(:count), ::Type) = Int

compile(fn::Fn(:json), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...) =
    let flow = compile(base, flow),
        ops = [compile(flow, op) for op in ops]
        isempty(ops) ? compile(fn, base, flow) : compile(fn, base, compile(Fn{:select}, base, flow, ops...))
    end

function compile(fn::Fn(:json), base::Scope, flow::Query)
    flow = select(flow)
    if !isnull(flow.scope.items)
        parts = ()
        fields = ()
        for item in get(flow.scope.items)
            field = compile(fn, flow.scope, compile(item, flow.scope))
            if !isnull(field.scope.tag)
                fields = (fields..., Pair{Symbol,AbstractPipe}(get(field.scope.tag), field.pipe))
            end
        end
        scope = nest(base, Dict{Any,Any})
        if !isnull(flow.scope.tag)
            scope = settag(scope, get(flow.scope.tag))
        end
        pipe = flow.pipe >> DictPipe(fields...)
        if ispartial(flow)
            pipe = NullToVoidPipe(pipe)
        end
        flow = Query(scope, pipe)
    else
        if ispartial(flow)
            pipe = NullToVoidPipe(flow.pipe)
            flow = Query(flow, pipe=pipe)
        end
    end
    return flow
end

function compile(::Fn(:frame), base::Scope, flow::Query)
    return Query(flow, pipe=BindRelPipe(flow.pipe))
end

function compile(::Fn(:given), base::Scope, flow::AbstractSyntax, ops::AbstractSyntax...)
    scope = base
    qs = Query[]
    tags = Symbol[]
    for op in ops
        q = compile(scope, op)
        reqtag(q)
        push!(qs, q)
        tag = gensym(get(q.scope.tag))
        push!(tags, tag)
        binding = ParamBinding(tag, output(q))
        scope = addglobal(scope, get(q.scope.tag), binding)
    end
    query = compile(scope, flow)
    pipe = query.pipe
    reverse!(qs)
    reverse!(tags)
    for (q, tag) in zip(qs, tags)
        pipe = BindEnvPipe(pipe, tag, q.pipe)
    end
    # FIXME: remove parameters from the query scope.
    return Query(query, pipe=pipe)
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module TestIteration

import Compat
using Test, DataFrames

@testset "eachrow and eachcol" begin
    df = DataFrame(A = Vector{Union{Int, Missing}}(1:2), B = Vector{Union{Int, Missing}}(2:3))

    @test size(eachrow(df)) == (size(df, 1),)
    @test parent(eachrow(df)) === df
    @test names(eachrow(df)) == names(df)
    @test IndexStyle(eachrow(df)) == IndexLinear()
    @test sprint(summary, eachrow(df)) == "2-element DataFrameRows"
    @test Base.IndexStyle(eachrow(df)) == IndexLinear()
    @test eachrow(df)[1] == DataFrameRow(df, 1, :)
    @test collect(eachrow(df)) isa Vector{<:DataFrameRow}
    @test eltype(eachrow(df)) <: DataFrameRow
    for row in eachrow(df)
        @test isa(row, DataFrameRow)
        @test (row[:B] - row[:A]) == 1
        # issue #683 (https://github.com/JuliaData/DataFrames.jl/pull/683)
        @test collect(pairs(row)) isa Vector{Pair{Symbol, Int}}
    end

    @test size(eachcol(df)) == (size(df, 2),)
    @test parent(eachcol(df)) === df
    @test names(eachcol(df)) == names(df)
    @test IndexStyle(eachcol(df)) == IndexLinear()
    @test Base.IndexStyle(eachcol(df)) == IndexLinear()
    @test length(eachcol(df)) == size(df, 2)
    @test eachcol(df)[1] == df[:, 1]
    @test collect(eachcol(df)) isa Vector{AbstractVector}
    @test collect(eachcol(df)) == [[1, 2], [2, 3]]
    @test eltype(eachcol(df)) == AbstractVector
    for col in eachcol(df)
        @test isa(col, AbstractVector)
    end

    @test map(x -> minimum(convert(Vector, x)), eachrow(df)) == [1,2]
    @test map(Vector, eachrow(df)) == [[1, 2], [2, 3]]
    @test mapcols(minimum, df) == DataFrame(A = [1], B = [2])
    @test map(minimum, eachcol(df)) == [1, 2]
    @test eltype.(eachcol(mapcols(Vector{Float64}, df))) == [Float64, Float64]
    @test eltype(map(Vector{Float64}, eachcol(df))) == Vector{Float64}

    @test_throws MethodError for x in df; end
end

@testset "mapcols" begin
    df_mapcols = DataFrame(a=1:10, b=11:20)
    @test mapcols(sum, df_mapcols) == DataFrame(a=55, b=155)
    @test mapcols(x -> 1, df_mapcols) == DataFrame(a=1, b=1)
    @test_throws ArgumentError mapcols(x -> x[1] == 1 ? 0 : [0], df_mapcols)
    @test_throws ArgumentError mapcols(x -> x[1] == 1 ? x : 0, df_mapcols)
    @test_throws ArgumentError mapcols(x -> x[1] != 1 ? x : 0, df_mapcols)
    df_mapcols2 = mapcols(x -> x, df_mapcols)
    @test df_mapcols2 == df_mapcols
    @test df_mapcols2.a !== df_mapcols.a
    @test df_mapcols2.b !== df_mapcols.b
end

@testset "SubDataFrame" begin
    df = DataFrame([11:16 21:26 31:36 41:46])
    sdf = view(df, [3,1,4], [3,1,4])
    @test sdf == df[[3,1,4], [3,1,4]]
    @test eachrow(sdf) == eachrow(df[[3,1,4], [3,1,4]])
    @test size(eachrow(sdf)) == (3,)
    @test eachcol(sdf) == eachcol(df[[3,1,4], [3,1,4]])
    @test size(eachcol(sdf)) == (3,)
end

@testset "parent mutation" begin
    df = DataFrame([11:16 21:26 31:36 41:46])
    sdf = view(df, [3,1,4], [3,1,4])
    erd = eachrow(df)
    erv = eachrow(sdf)
    rename!(df, Symbol.(string.("y", 1:4)))
    df[!, 1] = 51:56
    @test df[1, :] == erd[1]
    @test copy(erv[1]) == (y3=33, y1=53, y4=43)
    df[!, :z] .= 1
    @test length(erd[1]) == 5 # the added column is reflected
    select!(df, Not([4,5]))
    @test copy(erd[1]) == (y1 = 51, y2 = 21, y3 = 31) # the removed columns are reflected
end

@testset "getproperty and propertynames" begin
    df_base = DataFrame([11:16 21:26 31:36 41:46])
    for df in (df_base, view(df_base, 1:3, 1:3))
        for x in (eachcol(df), eachrow(df))
            @test propertynames(x) == propertynames(df)
            for n in names(df)
                @test getproperty(x, n) === getproperty(df, n)
            end
            @test_throws ArgumentError x.a
        end
    end
end

@testset "overload Compat functions" begin
    @testset "DataFrames.$f === Compat.$f" for f in intersect(names(DataFrames), names(Compat))
        @test getproperty(DataFrames, f) === getproperty(Compat, f)
    end
end

@testset "keys and pairs for eachcol" begin
    df = DataFrame([11:16 21:26 31:36 41:46])

    cols = eachcol(df)

    @test keys(cols) == names(df)
    for (a, b, c) in zip(keys(cols), cols, pairs(cols))
        @test (a => b) == c
    end

    for (i, n) in enumerate(keys(cols))
        @test cols[i] === cols[n]
    end
    @test_throws ArgumentError cols[:non_existent]
end

end # module


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm
import Base.show

function show(filename)
	V, EV = Plasm.svg2lar(filename)
	#Plasm.view(Plasm.numbering(.075)((V,[[[k] for k=1:size(V,2)], EV])))
	V, EV = Plasm.svg2lar(filename, flag=false)
	#Plasm.view(Plasm.numbering(2)((V,[[[k] for k=1:size(V,2)], EV])))
	return V, EV
end

#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/new.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/curved.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/twopaths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/paths.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/boundarytest2.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/tile.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/interior.svg")
#show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/holes.svg")

V,EV = show("/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/holes.svg")

V = Plasm.normalize(V,flag=true)
Plasm.view(Plasm.numbering(.125)((V,[[[k] for k=1:size(V,2)], EV])))

# subdivision of input edges
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV = Lar.planar_arrangement_1(W::Lar.Points, cop_EW::Lar.ChainOp)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# visualization of component graphs
EW = Lar.cop2lar(copEV)
W = convert(Lar.Points, V')
hpcs = [ Plasm.lar2hpc(W,[EW[e] for e in comp]) for comp in bicon_comps ]
Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k]) for k=1:(length(hpcs)) ])
Plasm.view(Plasm.numbering(.125)((W,[[[k] for k=1:size(W,2)], EW])))

# computation of 2-cells and 2-boundary
V, copEV, copFE = Lar.planar_arrangement_2(V, copEV, bicon_comps)

# visualization of numbered arrangement
Plasm.view( Plasm.numbering1(0.125)((V, copEV, copFE)) )

# final solid visualization
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})

EVs = Lar.FV2EVs(copEV, copFE) # polygonal face boundaries
EVs = convert(Array{Array{Array{Int64,1},1},1}, EVs)
Plasm.viewcolor(V::Lar.Points, EVs::Array{Lar.Cells})

model = V,EVs
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))


text:
#=
Stack
=#

using SimpleDataStructures

ss = SimpleStack{Int}()

push!(ss, 1)
push!(ss, 2)
push!(ss, 3)
pop!(ss)
pop!(ss)
pop!(ss)
pop!(ss)


text:
"""
    Kasai's algorithm for linear-time construction of LCP array from Suffix Array
"""
function lcparray{T<:Integer}(sa::Vector{Int}, text::Vector{T})
    n = length(sa)
    lcps = similar(sa)
    rank = similar(sa)
    for i = 1:n
        rank[sa[i]] = i
    end

    lcp = 0
    for i = 1:n
        rank[i] > 1 || continue
        j = sa[rank[i]-1]
        while i+lcp <= n && j+lcp <= n && text[i+lcp] == text[j+lcp]
            lcp += 1
        end
        lcps[rank[i]] = lcp
        lcp > 0 && (lcp -= 1)
    end
    lcps[1] = 0
    lcps
end


text:
"""
`LoaderError` should be thrown when loader library code fails, and other libraries should
be given the chance to recover from the error.  Reports the library name and an error message:
LoaderError("ImageMagick", "Foo not available")
"""
struct LoaderError <: Exception
    library::String
    msg::String
end
Base.showerror(io::IO, e::LoaderError) = println(io, e.library, " load error: ",
                                                 e.msg, "\n  Will try next loader.")

"""
`WriterError` should be thrown when writer library code fails, and other libraries should
be given the chance to recover from the error.  Reports the library name and an error message:
WriterError("ImageMagick", "Foo not available")
"""
struct WriterError <: Exception
    library::String
    msg::String
end
Base.showerror(io::IO, e::WriterError) = println(
    io, e.library, " writer error: ",
    e.msg, "\n  Will try next writer."
)

"""
`NotInstalledError` should be thrown when a library is currently not installed.
"""
struct NotInstalledError <: Exception
    library::Symbol
    message::String
end
Base.showerror(io::IO, e::NotInstalledError) = println(io, e.library, " is not installed.")

"""
`UnknownFormat` gets thrown when FileIO can't recognize the format of a file.
"""
struct UnknownFormat{T <: Formatted} <: Exception
    format::T
end
Base.showerror(io::IO, e::UnknownFormat) = println(io, e.format, " couldn't be recognized by FileIO.")

"""
Handles error as soon as they get thrown while doing IO
"""
function handle_current_error(e, library, islast::Bool)
    bt = catch_backtrace()
    bts = sprint(io->Base.show_backtrace(io, bt))
    message = islast ? "" : "\nTrying next loading library! Please report this issue on the Github page for $library"
    @warn string(e, bts, message)
end
handle_current_error(e::NotInstalledError) = @warn string("lib ", e.library, " not installed, trying next library")

"""
Handles a list of thrown errors after no IO library was found working
"""
function handle_exceptions(exceptions::Vector, action)
    # first show all errors when there are more then one
    multiple = length(exceptions) > 1
    println(stderr, "Error$(multiple ? "s" : "") encountered while $action.")
    if multiple
        println("All errors:")
        println("===========================================")
        for (err, file) in exceptions
            showerror(stdout, err)
            println("\n===========================================")
        end
    end
    # then handle all errors.
    # this way first fatal exception throws and user can still see all errors
    # TODO, don't throw, if it contains a NotInstalledError?!
    println(stderr, "\nFatal error:")
    for exception in exceptions
        continue_ = handle_error(exception...)
        continue_ || break
    end
end

handle_error(e, q) = throw(e)

function handle_error(e::NotInstalledError, q)
    println("Library \"", e.library, "\" is not installed but is recommended as a library to load format: \"", file_extension(q), "\"")
    !isinteractive() && rethrow(e) # if we're not in interactive mode just throw
    while true
        println("Should we install \"", e.library, "\" for you? (y/n):")
        input = lowercase(chomp(strip(readline(stdin))))
        if input == "y"
            @info(string("Start installing ", e.library, "..."))
            Pkg.add(string(e.library))
            return false # don't continue
        elseif input == "n"
            @info(string("Not installing ", e.library))
            return true # User does not install, continue going through errors.
        else
            println("$input is not a valid choice. Try typing y or n")
        end
    end
    true # User does not install, continue going through errors.
end


text:
# Reference:
# [Kobayashi 2014] Kobayashi, Kazuhiro, et al. "Statistical Singing Voice
# Conversion with Direct Waveform Modification based on the Spectrum
# Differential." Fifteenth Annual Conference of the International Speech
# Communication Association. 2014.

# DIFFGMM construction based on [Kobayashi 2014]
# TODO: tests
function diffgmm(params::GMMMapParam)
    μˣ = params.μˣ
    μʸ = params.μʸ
    Σˣˣ = params.Σˣˣ
    Σˣʸ = params.Σˣʸ
    Σʸˣ = params.Σʸˣ
    Σʸʸ = params.Σʸʸ

    GMMMapParam(params.weights,
                μˣ,
                μʸ-μˣ, # Eq. (6) in [Kobayashi 2014]
                Σˣˣ,
                Σˣʸ - Σˣˣ, # Eq. (7)
                (Σˣʸ - Σˣˣ)',
                Σˣˣ + Σʸʸ - Σˣʸ - Σʸˣ # Eq. (8)
                )
end