text:
"""
A Fourier basis on the interval `[0,1]`. The basis functions are given by
`exp(2 π i k)`, with `k` ranging from `-N` to `N` for Fourier series of odd
length (`2N+1`) and from `-N+1` to `N` for even length. In the latter case, the
highest frequency basis function corresponding to frequency `N` is a cosine.

The basis functions are ordered the way they are expected by a typical FFT
implementation. The frequencies k are in the following order:
```
0 1 2 3 ... N -N -N+1 ... -2 -1,
```
for odd length and
```
0 1 2 3 ... N -N+1 ... -2 -1,
```
for even length.

The Fourier basis is orthonormal with respect to a continuous measure (for odd
length Fourier bases only) and a discrete measure.
"""
struct Fourier{T <: Real} <: Dictionary{T,Complex{T}}
	n	::	Int
end

Fourier(n::Int) = Fourier{Float64}(n)

convert(::Type{Fourier{T}}, d::Fourier) where {T} = Fourier{T}(d.n)

@deprecate Fourier(n, a::Number, b::Number) Fourier(n)→a..b
@deprecate Fourier{T}(n, a::Number, b::Number) where {T} Fourier(n)→a..b

size(b::Fourier) = (b.n,)

oddlength(b::Fourier) = isodd(length(b))
evenlength(b::Fourier) = iseven(length(b))

similar(b::Fourier, ::Type{T}, n::Int) where {T} = Fourier{T}(n)

show(io::IO, b::Fourier{Float64}) = print(io, "Fourier($(length(b)))")
show(io::IO, b::Fourier{T}) where T = print(io, "Fourier{$(T)}($(length(b)))")

# Properties

isreal(b::Fourier) = false

isbasis(b::Fourier) = true

isorthogonal(b::Fourier, μ::FourierWeight) = true
isorthogonal(b::Fourier, μ::Weight) = islebesguemeasure(μ) && support(μ) == support(b)
isorthogonal(b::Fourier, μ::DiscreteWeight) = isuniform(μ) && compatible_domains(b, points(μ))

isorthonormal(b::Fourier, μ::FourierWeight) = oddlength(b)
isorthonormal(b::Fourier, μ::Weight) = isorthogonal(b, μ) && oddlength(b)
isorthonormal(b::Fourier, μ::DiscreteWeight) =
	isorthogonal(b, μ) && isnormalized(μ) && (length(b)==length(μ) || oddlength(b))

isbiorthogonal(b::Fourier) = true

hasinterpolationgrid(b::Fourier) = true

hasextension(b::Fourier) = isodd(length(b))

# For hastransform we introduce some more functionality:
# - Check whether the given periodic equispaced grid is compatible with the FFT operators
iscompatible(dict::Fourier, grid::AbstractEquispacedGrid) =
	compatible_domains(dict, grid) && length(dict)==length(grid)
# - Fourier grids are of course okay
iscompatible(dict::Fourier, grid::FourierGrid) = length(dict)==length(grid)
# - Any non-periodic grid is not compatible
iscompatible(dict::Fourier, grid::AbstractGrid) = false
# - We have a transform if the grid is compatible
hasgrid_transform(dict::Fourier, gb, grid) = iscompatible(dict, grid)

compatible_domains(dict::Fourier, grid) = false
compatible_domains(dict::Fourier, grid::AbstractEquispacedGrid) =
	isperiodic(grid) && support(dict) ≈ coverdomain(grid)

interpolation_grid(b::Fourier{T}) where {T} = FourierGrid{T}(length(b))

gauss_rule(dict::Fourier) = NormalizedDiracComb(interpolation_grid(dict))

struct MappedFourier{T} <: MappedDict{T,Complex{T}}
	superdict	::	Fourier{T}
	map			::	ScalarAffineMap{T}
end

mapped_dict(dict::Fourier{T}, map::ScalarAffineMap{S}) where {S,T} =
	MappedFourier{promote_type(S,T)}(dict, map)

##################
# Native indices
##################

const FourierFrequency = NativeIndex{:fourier}
const FFreq = FourierFrequency

frequency(idxn::FourierFrequency) = value(idxn)

Base.show(io::IO, idx::NativeIndex{:fourier}) =
	print(io, "Fourier frequency $(value(idx))")

"""
`FFTIndexList` defines the map from native indices to linear indices
for a finite Fourier basis, when the indices are ordered in the way they
are expected in the FFT routine.
"""
struct FFTIndexList <: IndexList{FourierFrequency}
	n	::	Int
end

length(list::FFTIndexList) = list.n
size(list::FFTIndexList) = (list.n,)

# The frequency of an even Fourier basis ranges from -N+1 to N.
# The frequency of an odd Fourier basis ranges from -N to N.
# This makes for a small difference in the ordering.
function getindex(m::FFTIndexList, idx::Int)
	n = length(m)
	nhalf = n >> 1
	if idx <= nhalf+1
		FourierFrequency(idx-1)
	else
		if iseven(n)
			FourierFrequency(idx-2*nhalf-1)
		else
			FourierFrequency(idx-2*nhalf-2)
		end
	end
end

function getindex(list::FFTIndexList, idxn::FourierFrequency)
	k = value(idxn)
	k >= 0 ? k+1 : length(list)+k+1
end

ordering(b::Fourier) = FFTIndexList(length(b))

# Shorthand: compute the linear index based on the size and element type
# of an array only
linear_index(idxn::FourierFrequency, size::Tuple{Int}, T) = FFTIndexList(size[1])[idxn]

# Convenience: compute with integer frequencies, rather than FourierFrequency types
idx2frequency(b::Fourier, idx) = frequency(native_index(b, idx))
frequency2idx(b::Fourier, k) = linear_index(b, FourierFrequency(k))

nhalf(b::Fourier) = nhalf(length(b))
nhalf(n::Int) = n>>1

maxfrequency(b::Fourier) = nhalf(b)
minfrequency(b::Fourier) = oddlength(b) ? -nhalf(b) : -nhalf(b)+1

#############
# Evaluation
#############

support(b::Fourier{T}) where T = UnitInterval{T}()

hasmeasure(b::Fourier) = true
measure(b::Fourier{T}) where T = FourierWeight{T}()

period(b::Fourier{T}) where {T} = T(1)

# One has to be careful here not to match Floats and BigFloats by accident.
# Hence the conversion of pi to T here:
exponent(b::Fourier{T}, k::FourierFrequency) where {T} = 2*T(pi)*im*k
exponent(b::Fourier, i::Int) = exponent(b, native_index(b, i))

# Even-length Fourier series have a cosine at the maximal frequency.
iscosine(b::Fourier, k::FourierFrequency) = iseven(length(b)) && (k==length(b)>>1)

function unsafe_eval_element(b::Fourier, idxn::FourierFrequency, x)
	z = exp(exponent(b, idxn)*x)
	if iscosine(b, idxn)
		# The cosine is the real part of the exponential, but we make it a complex
		# number for type-stability
		z = complex(real(z))
	end
	z
end

function unsafe_eval_element_derivative(b::Fourier, idxn::FourierFrequency, x, order)
	@assert order >= 0
	f = exponent(b, idxn)
	z = f^order * exp(f * x)
	if iscosine(b, idxn)
		# We have to take the real part for the cosine, yet return a complex number
		z = complex(real(z))
	end
	z
end

function unsafe_moment(b::Fourier, idxn::FourierFrequency)
	T = codomaintype(b)
	frequency(idx) == 0 ? T(1) : T(0)
end

twiddle(T, x) = exp(T(pi)*2im*x)

# For large n, the routine below is slightly faster than evaluating all
# the complex exponentials. It is based on the fact that the Fourier exponentials
# are integer powers of an exponential.
function unsafe_dict_eval!(result, dict::Fourier, x)
	result[1] = 1
	n1 = nhalf(dict)
	w = twiddle(domaintype(dict), x)
	z = one(w)
	for i in 1:n1
		z *= w
		result[i+1] = z
		result[end-i+1] = conj(z)
	end
	if evenlength(dict)
		result[n1+1] = real(result[n1+1])
	end
	result
end

# By default, we preserve the odd/even property of the size when extending
extensionsize(b::Fourier) = evenlength(b) ? 2*length(b) : max(3,2*length(b)-1)

approx_length(b::Fourier, n::Int) = n

"Shift an expansion to the right by delta."
function shift(b::Fourier, coefficients, delta)
	# Only works for odd-length Fourier series for now
	@assert oddlength(b)
	S = domaintype(b)
	coef2 = copy(coefficients)
	for i in eachindex(coefficients)
		coef2[i] *= exp(2 * S(pi) * im * idx2frequency(b, i) * delta)
	end
	coef2
end

function transform_from_grid(T, src::GridBasis, dest::Fourier, grid; options...)
	@assert iscompatible(dest, grid)
	forward_fourier_operator(src, dest, T; options...)
end

function transform_to_grid(T, src::Fourier, dest::GridBasis, grid; options...)
	@assert iscompatible(src, grid)
	inverse_fourier_operator(src, dest, T; options...)
end

evaluation(::Type{T}, dict::Fourier, gb::GridBasis, grid::FourierGrid; options...) where {T} =
	resize_and_transform(T, dict, gb, grid; options...)

function evaluation(::Type{T}, dict::Fourier, gb::GridBasis, grid::PeriodicEquispacedGrid; options...) where {T}
	if coverdomain(grid)≈support(dict)
		resize_and_transform(T, dict, gb, grid; options...)
	else
		@debug "Periodic grid mismatch with Fourier basis"
		dense_evaluation(T, dict, gb; options...)
	end
end

to_periodic_grid(dict::Fourier, grid::AbstractGrid) = nothing
to_periodic_grid(dict::Fourier, grid::PeriodicEquispacedGrid{T}) where {T} =
	iscompatible(dict, grid) ? FourierGrid{T}(length(grid)) : nothing

function evaluation(::Type{T}, dict::Fourier, gb::GridBasis, grid; options...) where {T}
	grid2 = to_periodic_grid(dict, grid)
	if grid2 != nothing
		gb2 = GridBasis{T}(grid2)
		evaluation(T, dict, gb2, grid2; options...) * gridconversion(gb, gb2; options...)
	else
		@debug "Evaluation: could not convert $(string(grid)) to periodic grid"
		dense_evaluation(T, dict, gb; options...)
	end
end

# Try to efficiently evaluate a Fourier series on a regular equispaced grid
function evaluation(::Type{T}, fs::Fourier, gb::GridBasis, grid::EquispacedGrid; options...) where {T}
	# We can use the fft if the equispaced grid is a subset of the periodic grid
	if coverdomain(grid) ≈ support(fs)
		# TODO: cover the case where the EquispacedGrid is like a PeriodicEquispacedGrid
		# but with the right endpoint added
		return dense_evaluation(T, fs, gb; options...)
	elseif issubset(coverdomain(grid), support(fs))
		a, b = endpoints(coverdomain(grid))
		if a==b
			return dense_evaluation(T, fs, gb; options...)
		end
		# We are dealing with a subgrid. The main question is: if we extend it
		# to the full support, is it compatible with a periodic grid?
		h = step(grid)
		nleft = a/h
		nright = (1-b)/h
		if (nleft ≈ round(nleft)) && (nright ≈ round(nright))
			nleft_int = round(Int, nleft)
			nright_int = round(Int, nright)
			ntot = length(grid) + nleft_int + nright_int - 1
			super_grid = FourierGrid{domaintype(fs)}(ntot)
			super_gb = GridBasis{coefficienttype(gb)}(super_grid)
			E = evaluation(T, fs, super_gb; options...)
			R = IndexRestriction{T}(super_gb, gb, nleft_int+1:nleft_int+length(grid))
			return R*E
		end
		dense_evaluation(T, fs, gb; options...)
	else
		dense_evaluation(T, fs, gb; options...)
	end
end

function evaluation(::Type{T}, dict::Fourier, gb::GridBasis, grid::MidpointEquispacedGrid; options...) where {T}
	if isodd(length(grid)) && coverdomain(grid)≈support(dict)
		if length(grid) == length(dict)
			A = evaluation(T, dict, FourierGrid{domaintype(dict)}(length(dict)); options...)
			diag = zeros(T, length(dict))
			delta = step(grid)/2
			for i in 1:length(dict)
				diag[i] = exp(2 * T(pi) * im * idx2frequency(dict, i) * delta)
			end
			D = DiagonalOperator(dict, diag)
			wrap_operator(dict, gb, A*D)
		else
			dict2 = resize(dict, length(grid))
			evaluation(T, dict2, grid) * extension(T, dict, dict2)
		end
	else
		dense_evaluation(T, dict, gb; options...)
	end
end

############################
# Extension and restriction
############################

# We make special-purpose operators for the extension of a Fourier series,
# since we have to add zeros in the middle of the given Fourier coefficients.
# This can not be achieved with a single IndexExtension.
# It could be a composition of two, but this special case is widespread and hence
# we make it more efficient.
struct FourierIndexExtension{T} <: DictionaryOperator{T}
	src		::	Dictionary
	dest	::	Dictionary
	n1		::	Int
	n2		::	Int

	function FourierIndexExtension{T}(src::Dictionary, dest::Dictionary, n1::Int, n2::Int) where {T}
		@assert n1 <= n2
		new(src, dest, n1, n2)
	end
end

FourierIndexExtension(src::Dictionary, dest::Dictionary) =
	FourierIndexExtension{operatoreltype(src, dest)}(src, dest)

FourierIndexExtension{T}(src, dest) where {T} =
	FourierIndexExtension{T}(src, dest, length(src), length(dest))

string(op::FourierIndexExtension) = "Fourier series extension from length $(op.n1) to length $(op.n2) (T=$(eltype(op)))"

unsafe_wrap_operator(src, dest, op::FourierIndexExtension{T}) where {T} =
	FourierIndexExtension{T}(src, dest, op.n1, op.n2)

extension(::Type{T}, src::Fourier, dest::Fourier; options...) where {T} =
	FourierIndexExtension{T}(src, dest)

function apply!(op::FourierIndexExtension, coef_dest, coef_src)
	## @assert length(dest) > length(src)

	nh = nhalf(length(coef_src))
	# We have to distinguish between even and odd length Fourier series
	# because the memory layout is slightly different
	if iseven(length(coef_src))
		for i = 0:nh-1
			coef_dest[i+1] = coef_src[i+1]
		end
		for i = 1:nh-1
			coef_dest[end-nh+i+1] = coef_src[nh+1+i]
		end
		coef_dest[nh+1] = coef_src[nh+1]/2
		coef_dest[end-nh+1] = coef_src[nh+1]/2
		for i = nh+2:length(coef_dest)-nh
			coef_dest[i] = 0
		end
	else
		for i = 0:nh
			coef_dest[i+1] = coef_src[i+1]
		end
		for i = 1:nh
			coef_dest[end-nh+i] = coef_src[nh+1+i]
		end
		for i = nh+2:length(coef_dest)-nh
			coef_dest[i] = 0
		end
	end
	coef_dest
end

struct FourierIndexRestriction{T} <: DictionaryOperator{T}
	src		::	Dictionary
	dest	::	Dictionary
	n1		::	Int
	n2		::	Int

	function FourierIndexRestriction{T}(src::Dictionary, dest::Dictionary, n1::Int, n2::Int) where {T}
		@assert n1 >= n2
		new(src, dest, n1, n2)
	end
end

FourierIndexRestriction(src::Dictionary, dest::Dictionary) =
	FourierIndexRestriction{operatoreltype(src, dest)}(src, dest)

FourierIndexRestriction{T}(src, dest) where {T} =
	FourierIndexRestriction{T}(src, dest, length(src), length(dest))

string(op::FourierIndexRestriction) = "Fourier series restriction from length $(op.n1) to length $(op.n2) (T=$(eltype(op)))"

unsafe_wrap_operator(src, dest, op::FourierIndexRestriction{T}) where {T} =
	FourierIndexRestriction{T}(src, dest, op.n1, op.n2)

restriction(::Type{T}, src::Fourier, dest::Fourier; options...) where {T} =
	FourierIndexRestriction{T}(src, dest)

function apply!(op::FourierIndexRestriction, coef_dest, coef_src)
	## @assert length(dest) < length(src)

	nh = nhalf(length(coef_dest))
	if isodd(length(coef_dest))
		for i = 0:nh
			coef_dest[i+1] = coef_src[i+1]
		end
		for i = 1:nh
			coef_dest[nh+1+i] = coef_src[end-nh+i]
		end
	else
		for i = 0:nh-1
			coef_dest[i+1] = coef_src[i+1]
		end
		for i = 1:nh-1
			coef_dest[nh+1+i] = coef_src[end-nh+i+1]
		end
		# if isodd(length(coef_src))
			# coef_dest[nh+1] = (coef_src[nh+1] + coef_src[end-nh+1])
		# else
			coef_dest[nh+1] = (coef_src[nh+1] + coef_src[end-nh+1]) / 2 # Added  / 2 to get adjoint correct
		# end
	end
	coef_dest
end

isdiag(::FourierIndexExtension) = true
isdiag(::FourierIndexRestriction) = true

adjoint(op::FourierIndexExtension{T}) where {T} =
	FourierIndexRestriction{T}(dest(op), src(op), op.n2, op.n1)

adjoint(op::FourierIndexRestriction{T}) where {T} =
	FourierIndexExtension{T}(dest(op), src(op), op.n2, op.n1)::DictionaryOperator

###################
# Differentiation
###################

# The Fourier series of odd length are closed under differentiation of any order.
# Fourier series of even length contain a cosine and they must be extended to
# the next odd length in order to exactly represent its derivative.

# We have a derivative of any integer order
hasderivative(b::Fourier) = true
hasderivative(b::Fourier, order::Int) = order >= 0

function derivative_dict(Φ::Fourier, order::Int; options...)
	@assert order >= 0
	if order == 0
		Φ
	else
		oddlength(Φ) ? Φ : resize(Φ, length(Φ)+1)
	end
end

# Both differentiation and antidifferentiation are diagonal operations
diff_scaling_function(b::Fourier, idx, symbol) = symbol(exponent(b, idx))
diff_scaling_function(b::Fourier, idx, order::Int) = diff_scaling_function(b,idx,x->x^order)

function antidiff_scaling_function(b::Fourier, idx, order::Int)
	T = domaintype(b)
	k = idx2frequency(b, idx)
	k==0 ? Complex{T}(0) : 1 / (k * 2 * T(pi) * im)^order
end

function differentiation(::Type{T}, src::Fourier, dest::Fourier, order::Int; options...) where {T}
	if orderiszero(order)
		@assert src==dest
		IdentityOperator{T}(src)
	else
		if oddlength(src)
			@assert src == dest
			pseudodifferential_operator(T, src, dest, x->x^order; options...)
		else
			@assert length(dest) == length(src)+1
			E = extension(T, src, dest)
			D = differentiation(T, dest, dest, order; options...)
			D*E
		end
	end
end

pseudodifferential_operator(src::Dictionary, args...; options...) =
	pseudodifferential_operator(operatoreltype(src), src, args...; options...)

pseudodifferential_operator(::Type{T}, src::Fourier, symbol::Function; options...) where {T} =
	pseudodifferential_operator(T, src, derivative_dict(src; options...), symbol; options...)

function pseudodifferential_operator(::Type{T}, src::Fourier, dest::Fourier, symbol::Function; options...) where {T}
	if oddlength(src)
		@assert src == dest
		_pseudodifferential_operator(T, src, dest, symbol; options...)
	else
		@assert length(dest) == length(src)+1
		_pseudodifferential_operator(T, src, dest, symbol; options...) * extension(T, src, dest; options...)
	end
end

_pseudodifferential_operator(::Type{T}, src, dest, symbol::Function; options...) where {T} =
	DiagonalOperator{T}([diff_scaling_function(src, idx, symbol) for idx in eachindex(src)], src, dest)

##########################
# Arithmetical operations
##########################

# Multiplication of Fourier Series
function (*)(src1::Fourier, src2::Fourier, coef_src1, coef_src2)
	T = promote_type(eltype(coef_src1), eltype(coef_src2))
	if oddlength(src1) && evenlength(src2)
	    dsrc2 = resize(src2, length(src2)+1)
	    (*)(src1, dsrc2, coef_src1, extension(T, src2, dsrc2)*coef_src2)
	elseif evenlength(src1) && oddlength(src2)
		dsrc1 = resize(src1, length(src1)+1)
	    (*)(dsrc1, src2, extension(T, src1,dsrc1)*coef_src1,coef_src2)
	elseif evenlength(src1) && evenlength(src2)
		dsrc1 = resize(src1, length(src1)+1)
	    dsrc2 = resize(src2, length(src2)+1)
		T1 = eltype(coef_src1)
		T2 = eltype(coef_src2)
	    (*)(dsrc1,dsrc2,extension(T, src1, dsrc1)*coef_src1, extension(T, src2, dsrc2)*coef_src2)
	else # they are both odd
		@assert domaintype(src1) == domaintype(src2)
	    dest = Fourier{domaintype(src1)}(length(src1)+length(src2)-1)
	    coef_src1 = [coef_src1[(nhalf(src1))+2:end]; coef_src1[1:nhalf(src1)+1]]
	    coef_src2 = [coef_src2[(nhalf(src2))+2:end]; coef_src2[1:nhalf(src2)+1]]
	    coef_dest = conv(coef_src1,coef_src2)
	    coef_dest = [coef_dest[(nhalf(dest)+1):end]; coef_dest[1:(nhalf(dest))]]
	    (dest,coef_dest)
	end
end

## Inner products

# Evaluate the inner product of two Fourier basis functions on the full domain
function innerproduct_fourier_full(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq)
	T = codomaintype(b1)
	# The outcome can be zero, one or 1/2.
	# The latter happens when integrating a cosine with another cosine, or with
	# a complex exponential that contains a cosine with the same frequency.
	# Apart from that special case, the integral is zero when the indices differ
	# and one when they are equal.
	if (iscosine(b1, i) && abs(j) == i) || (iscosine(b2, j) && abs(i) == j)
		one(T)/2
	elseif i != j
		zero(T)
	else
		one(T)
	end
end

# Evaluate the inner product on the interval [a,b]
function innerproduct_fourier_part(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq, a, b)
	S = domaintype(b1)
	T = codomaintype(b1)
	# convert 2π to type S and 2πi to type T
	twopi = 2*S(pi)
	tpi = 2im*T(pi)
	if iscosine(b1, i)
		if iscosine(b2, j)
			if abs(i) == abs(j)
				T(-cos(twopi*i*a)*sin(twopi*i*a)+cos(twopi*i*b)*sin(twopi*i*b)-twopi*i*(a-b))/(2*twopi*i)
			else
				T((-i-j)*sin(twopi*(i-j)*a)+(i+j)*sin(twopi*(i-j)*b)-(sin(twopi*(i+j)*a)-sin(twopi*(i+j)*b))*(i-j))/(2*i^2*twopi-2*j^2*twopi)
			end
		else
			if i == j
				((-cos(twopi*i*a)-sin(twopi*i*a))*exp(tpi*a)+exp(tpi*b)*(sin(twopi*i*b)+cos(twopi*i*b)))/(2*twopi*i)
			else
				(-T(im)*j*exp(tpi*j*a)*cos(twopi*i*a)+T(im)*j*exp(tpi*j*b)*cos(twopi*i*b)-i*exp(tpi*j*a)*sin(twopi*i*a)+i*exp(tpi*j*b)*sin(twopi*i*b))/((2*i^2-2*j^2)*S(pi))
			end
		end
	else
		if iscosine(b2, j)
			if abs(i) == abs(j)
				(-T(im)*cos(twopi*i*a)^2-cos(twopi*i*a)*sin(twopi*i*a)+T(im)*cos(twopi*i*b)^2+cos(twopi*i*b)*sin(twopi*i*b)-twopi*i*(a-b))/(2*twopi*i)
			else
				(-T(im)*i*exp(-tpi*i*a)*cos(twopi*j*a)+T(im)*i*exp(-tpi*i*b)*cos(twopi*j*b)+j*exp(-tpi*i*a)*sin(twopi*j*a)-j*exp(-tpi*i*b)*sin(twopi*j*b))/((2*i^2-2*j^2)*pi)
			end
		else
			if i == j
				T(b-a)
			else
				-T(1im)/(twopi*(j-i)) * (exp(tpi*b*(j-i)) - exp(tpi*a*(j-i)))
			end
		end
	end
end

# For the uniform measure on [0,1], invoke innerproduct_fourier_full
innerproduct_native(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq, m::FourierWeight; options...) =
	innerproduct_fourier_full(b1, i, b2, j)

function innerproduct_native(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq, m::LebesgueMeasure; options...)
	d = support(m)
	if d isa AbstractInterval
		if infimum(d) == 0 && supremum(d) == 1
			innerproduct_fourier_full(b1, i, b2, j)
		else
			innerproduct_fourier_part(b1, i, b2, j, infimum(d), supremum(d))
		end
	else
		default_dict_innerproduct(b1, i, b2, j, m)
	end
end

function gram(::Type{T}, dict::Fourier, measure::FourierWeight; options...) where {T}
	@assert isorthogonal(dict, measure) # some robustness.
	if iseven(length(dict))
		CoefficientScalingOperator{T}(dict, (length(dict)>>1)+1, one(T)/2)
	else
		@assert isorthonormal(dict, measure)
		IdentityOperator{T}(dict, dict)
	end
end

function gram(::Type{T}, dict::Fourier, measure::DiscreteWeight, grid::AbstractEquispacedGrid, weights::FillArrays.AbstractFill; options...) where {T}
	if coverdomain(grid) ≈ support(dict) && isperiodic(grid)
		if isorthonormal(dict, measure)
			IdentityOperator{T}(dict)
		elseif isorthogonal(dict, measure)
			if isodd(length(dict)) || (length(dict)==length(grid))
				ScalingOperator{T}(dict, unsafe_weight(measure, 1)*length(grid))
			else
				CoefficientScalingOperator{T}(dict, (length(dict)>>1)+1, one(T)/2)*ScalingOperator{T}(dict, weights[1]*length(grid))
			end
		end
	else
		default_mixedgram_discretemeasure(T, dict, dict, measure, grid, weights; options...)
	end
end


text:
"""
This formulation of the duration constraints, adds over the start times looking backwards.

"""
function device_duration_retrospective(ps_m::CanonicalModel,
                                        duration_data::Tuple{Vector{String},Vector{UpDown}},
                                        initial_duration_on::Vector{InitialCondition},
                                        initial_duration_off::Vector{InitialCondition},
                                        time_range::UnitRange{Int64},
                                        cons_name::Symbol,
                                        var_names::Tuple{Symbol,Symbol,Symbol})

    set_name = duration_data[1]

    name_up = Symbol(cons_name,:_up)
    name_down = Symbol(cons_name,:_down)

    ps_m.constraints[name_up] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[name_down] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

        for t in time_range, (ix,name) in enumerate(set_name)

                if t - duration_data[2][ix].up >= 1
                    tst = duration_data[2][ix].up
                else
                    tst = max(1.0, duration_data[2][ix].up - initial_duration_on[ix].value)
                end

                if t - duration_data[2][ix].down >= 1
                    tsd = duration_data[2][ix].down
                else
                    tsd = max(1.0, duration_data[2][ix].down - initial_duration_off[ix].value)
                end

                ps_m.constraints[name_up][name, t] = JuMP.@constraint(ps_m.JuMPmodel, sum([ps_m.variables[var_names[2]][name,i] for i in ((t - round(tst) + 1) :t) if i > 0 ]) <= ps_m.variables[var_names[1]][name,t])
                ps_m.constraints[name_down][name, t] = JuMP.@constraint(ps_m.JuMPmodel, sum([ps_m.variables[var_names[3]][name,i] for i in ((t - round(tsd) + 1) :t) if i > 0]) <= (1 - ps_m.variables[var_names[1]][name,t]))

        end

    return

end

ext: jl


text:
using SparseArrays, FillArrays, Mill

"""
	ExtractOneHot(ks, k, v)

	Converts a Vector of `Dict` items to one-hot encoding by using key
	`k` to identify a name of item out of `ks` and values
	of key `v` as values.

```juliadoctest
julia> samples = ["{\"name\": \"a\", \"count\" : 1}",
		"{\"name\": \"b\", \"count\" : 2}",];
julia> samples = JSON.parse.(samples);
julia> e = ExtractOneHot(["a","b"], "name", "count");
julia> e(vs).data
3×1 SparseArrays.SparseMatrixCSC{Int64,Int64} with 2 stored entries:
  [1, 1]  =  1
  [2, 1]  =  2
```

	If `v` is equal to `nothing`, then it boils down to one-hot encoding
```juliadoctest
julia> e = ExtractOneHot(["a","b"], "name", nothing);
julia> e(vs).data
3×1 SparseArrays.SparseMatrixCSC{Int64,Int64} with 2 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
```

	If there is key in the data which is not known (it was not part of `vs`),
	than it is assigned to an special designed key serving as "unknown`
```juliadoctest
julia> vs = JSON.parse.(["{\"name\": \"c\", \"count\" : 1}"]);
julia> e = ExtractOneHot(["a","b"], "name", nothing);
julia> e(vs).data
3×1 SparseArrays.SparseMatrixCSC{Bool,Int64} with 1 stored entry:
  [3, 1]  =  1

```
"""
struct ExtractOneHot{K,I,V} <: AbstractExtractor
	k::K
	v::V
	key2id::Dict{I,Int}
	n::Int
end

extractsmatrix(s::ExtractOneHot) = false

function ExtractOneHot(ks::V, k, v = nothing) where {V<:Union{Vector, Set}}
	key2id = Dict(zip(ks, 1:length(ks)))
	ExtractOneHot(k, v, key2id, length(key2id) + 1)
end

(e::ExtractOneHot{K,V})(v::Dict) where {K, V <:Nothing} = ArrayNode(sparse([e.key2id[v[e.k]]], [1f0], true, e.n))
(e::ExtractOneHot{K,V})(v::Dict) where {K, V} = ArrayNode(sparse([e.key2id[v[e.k]]],  [1f0], [get(v, e.v, 0)], e.n))

function (e::ExtractOneHot{K,I,V})(vs::Vector) where {K, I, V<:Nothing}
	isempty(vs) && return(ArrayNode(spzeros(Float32, e.n, 1)))
	ids = [get(e.key2id, v[e.k], e.n) for v in vs]
	ArrayNode(sparse(ids, Ones(length(ids)), 1f0, e.n, 1))
end

function (e::ExtractOneHot{K,I,V})(vs::Vector) where {K, I, V}
	isempty(vs) && return(ArrayNode(spzeros(Float32, e.n, 1)))
	ids = [get(e.key2id, v[e.k], e.n) for v in vs]
	x = [Float32(v[e.v]) for v in vs]
	ArrayNode(sparse(ids, Ones(length(ids)), x, e.n, 1))
end

(e::ExtractOneHot)(::Nothing) = ArrayNode(spzeros(Float32, e.n, 1))

function Base.show(io::IO, m::ExtractOneHot; pad = [], key::String="")
	c = COLORS[(length(pad)%length(COLORS))+1]
	key *= isempty(key) ? "" : ": ";
	paddedprint(io,"$(key)OneHot d = $(m.n)\n", color = c)
end

Base.hash(e::ExtractOneHot, h::UInt) = hash((e.k, e.v, e.key2id, e.n), h)
Base.:(==)(e1::ExtractOneHot, e2::ExtractOneHot) = e1.k === e2.k && e1.v === e2.v && e1.key2id == e2.key2id && e1.n === e2.n


text:
module BigMat

using LinearAlgebra, Base, Random, Optim
    import LinearAlgebra: norm

    # functions that are available for the user
    export MPT, size, order, length, rank, norm,
           getelement, getrow, getcol, outerprod,
           MPT_SVD, TTm_ALS, TT_ALS, eye, rkrp2tn,
           mps2vec, mpo2mat,
           KathriRao, transpose, roundTT,
           mpo2mps, mps2mpo,
           unfold, diag, matrixbyvector, vectorbymatrix,
           TNrSVD,
           approxpseudoinverse

    # package code
    include("MPT.jl")

    include("functions_computeTN/MPT_SVD.jl")
    include("functions_computeTN/TT_SVD.jl")
    include("functions_computeTN/TT_ALS.jl")
    include("functions_computeTN/TTm_ALS.jl")
    include("functions_computeTN/eyeMPO.jl")
    include("functions_computeTN/rkrp2tn.jl")

    include("functions_tools/BasicAlgebra.jl")
    include("functions_tools/Contractions.jl")
    include("functions_tools/Unfoldings.jl")
    include("functions_tools/MatrixAlgebra.jl")
    include("functions_tools/roundTT.jl")
    include("functions_tools/TNrSVD.jl")

end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
"""
    maximum_xor_subarray(array::Array{Int64,1})

For finding the maximum XOR value of a subarray by running two nested for-loops. The first
for-loop is intializing the current XOR value, the second for-loop is generating the [BIT-XOR](https://en.wikipedia.org/wiki/Bitwise_operation#XOR)
and compared the current value with the temporary largest value.

# Arguments
- `array::Array{Int64,1}`: Unsorted array

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [8, 1, 2, 12, 45, 100, 43]
julia> ClassicAlgorithmsCollections.maximum_xor_subarray(arr)
110
```
"""
function maximum_xor_subarray(array::Array{Int64,1})
    result = -21345944
    n = length(array)

    for i in 1:n
        curr_xor = 0
        for j in i:n
            curr_xor = curr_xor ⊻ array[j]
            if result < curr_xor
                result = curr_xor
            end
        end
    end
    return result
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# Problem from https://exercism.io/my/tracks/julia#exercise-binary-search
# Based on my published solution https://exercism.io/tracks/julia/exercises/binary-search/solutions/5dcd7f36433245a28484889c42a957c0

# Problem Instructions:

"""
    binary_search(list, query; rev=false, lt=<, by=identity)

Implement a binary search algorithm.
Searching a sorted collection is a common task. A dictionary is a sorted list of word definitions. Given a word, one can find its definition. A telephone book is a sorted list of people's names, addresses, and telephone numbers. Knowing someone's name allows one to quickly find their telephone number and address.

If the list to be searched contains more than a few items (a dozen, say) a binary search will require far fewer comparisons than a linear search, but it imposes the requirement that the list be sorted.

In computer science, a binary search or half-interval search algorithm finds the position of a specified input value (the search "key") within an array sorted by key value.

In each step, the algorithm compares the search key value with the key value of the middle element of the array.

If the keys match, then a matching element has been found and the range of indices that equal the search key value are returned.

Otherwise, if the search key is less than the middle element's key, then the algorithm repeats its action on the sub-array to the left of the middle element or, if the search key is greater, on the sub-array to the right.

If the remaining array to be searched is empty, then the key cannot be found in the array and a special "not found" indication is returned. Search methods in Julia typically return an empty range located at the insertion point in this case.

A binary search halves the number of items to check with each iteration, so locating an item (or determining its absence) takes logarithmic time. A binary search is a dichotomic divide and conquer search algorithm.

Bonus task:
Implement keyword arguments by, lt and rev so that by specifies a transformation applied to all elements of the list, lt specifies a comparison and rev specifies if the list is ordered in reverse.

Contributed By:- [Soc Virnyl Estela](https://github.com/uncomfyhalomacro)
"""
function binary_search(list, query; rev = false, lt = <, by = identity)
    if issorted(list) || issorted(list; rev = true)
        low = !rev ? 1 : length(list)
        high = !rev ? length(list) : 1
        middle(l, h) = round(Int, (l + h) // 2)
        query = by(query)

        while !rev ? low <= high : high <= low
            mid = middle(low, high)
            by(list[mid]) === query && return mid:mid
            if lt(by(list[mid]), query)
                low = !rev ? mid + 1 : mid - 1
            else
                high = !rev ? mid - 1 : mid + 1
            end
        end
        return !rev ? (low:high) : (high:low)

    else
        throw(error("List not sorted, unable to search value"))
    end
end

# See Issue https://github.com/TheAlgorithms/Julia/issues/34
"""
    binary_search(arr::AbstractArray{T,1}, l::T, r::T, x::T) where {T<:Real}

The implementation of this binary Search is recursive and requires O(Log n) space. With iterative Binary Search, we need only O(1) space. Useful for the implementation of `exponential_search`.

Contributed By:- [Ash](https://github.com/ashwani-rathee)
"""
function binary_search(
    arr::AbstractArray{T,1},
    l::T,
    r::T,
    x::T,
) where {T<:Real}
    if (r >= l)
        mid = Int(ceil(l + (r - l) / 2))
        # println(mid)
        if (arr[mid] == x)
            return "Element present at index $mid"
        elseif (arr[mid] > x)
            binary_search(arr, l, mid - 1, x)
        else
            binary_search(arr, mid + 1, r, x)
        end
    else
        return "Element not present in array"
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
__precompile__(true)

module InvariantDistribution

using Simplices, SimplexSplitting, Parameters,
    StaticArrays, Distributions, InplaceOps

installed = Pkg.installed()
if !("Pycall" in keys(installed))
    Pkg.add("PyCall")
    ENV["PYTHON"]= ""; Pkg.build("PyCall")
end

if !("Conda" in keys(installed))
    Pkg.add("Conda")
    using Conda; Conda.add("scipy")
end

if !("Simplices" in keys(installed))
    Pkg.clone("https://github.com/kahaaga/Simplices.jl")
else
    using Simplices
end

include("invariant_maps.jl")
include("statespace.jl")
include("example_embeddings.jl")
include("invdist.jl")
include("invariantize_embedding.jl")
include("mm_exact/mm_sparse.jl")
include("mm_exact/mm_sparse_parallel.jl")
include("mm_exact/mm_parallel.jl")
include("mm_approx/mm_dd.jl")
include("mm_approx/mm_dd2.jl")
include("mm_approx/mm_dd4.jl")

export markovmatrix,
    markovmatrixp,
    estimate_invdist,
    StateSpace,
    indexin_rows,
    invariant_gaussian_embedding,
    te_test,
    mm_sparse,
    mm_sparse_parallel,
    mm_p,
    invariantize_embedding,
    invariant_bullmap_on_cube,
    prepare_mm_dd,
	mm_dd,
    mm_dd2,
    mm_dd3,
    mm_dd4

t = triang_from_embedding(Embedding(invariant_gaussian_embedding(npts = 8, covariance = 0.5, tau = 1)))
mm_dd4(t)

end # module

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
# Define logistic loss
# TODO
# [ ] Memoize expit(-yₚ[i] y[i])

import Base: length

# Define data for logistic regression problem
struct LogitData{T <: Real} <: AbstractDataset{T}
    # Features x_ij
    X::Array{T, 2}
    # Labels y_i
    y::Vector{T}
    # Current prediction. Should be equal to dot(X, p) with p
    # current parameters.
    y_pred::Vector{T}
end

function LogitData(X::Array{T, 2}, y::Vector{T}) where T
    @assert size(X, 1) == size(y, 1)
    return LogitData(X, y, zeros(T, size(y, 1)))
end

function LogitData(libsmv_file::String; scale_data=false)
    dataset = parse_libsvm(libsmv_file)
    # Convert to dense matrix
    X = to_dense(dataset)
    if scale_data
        scale!(NormalScaler(), X)
    end
    y = dataset.labels
    # Special preprocessing for covtype
    format_label!(y)
    return LogitData(X, y, zeros(size(y, 1)))
end

nfeatures(dat::LogitData) = size(dat.X, 2)
ndata(dat::LogitData) = length(dat.y)

function predict!(data::LogitData{T}, x::AbstractVector{T}, intercept=NaN) where T
    mul!(data.y_pred, data.X, x)
    if isfinite(intercept)
        data.y_pred .+= intercept
    end
    return
end

"""
Compute logistic loss of given vector parameter `ω ∈  R^p`.

## Complexity
O(n)

"""
function loss(ω::AbstractVector{T}, data::LogitData{T}) where T
    # Compute dot product
    res = zero(T)
    @inbounds for i in 1:ndata(data)
        res += loss(Logit(), data.y_pred[i], data.y[i])
    end
    return one(T) / ndata(data) * res
end

"""
Compute gradient of logistic loss for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p)

"""
function gradient!(grad::AbstractVector{T}, ω::AbstractVector{T},
                   data::LogitData{T}, fit_intercept::Bool=false) where T
    nfeats = nfeatures(data)
    n = ndata(data)
    invn = -one(T) / n
    @inbounds for j in 1:n
        tmp = invn * data.y[j] * expit(-data.y_pred[j] * data.y[j])
        for i in 1:nfeats
            grad[i] += tmp * data.X[j, i]
        end
        if fit_intercept
            grad[end] += tmp
        end
    end
    return nothing
end

"""
Compute Hessian of logistic loss for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p * (p-1) /2)

"""
function hessian!(hess::AbstractVector{T}, ω::AbstractVector{T}, data::LogitData{T},
                  fit_intercept::Bool=false) where T
    # Sanity check
    p = nfeatures(data)
    n = ndata(data)
    invn = one(T) / n
    @inbounds for i in 1:n
        σz = expit(-data.y_pred[i] * data.y[i])
        cst = invn * σz * (one(T) - σz)
        count = 1
        for j in 1:p
            for k in j:p
                hess[count] += cst * data.X[i, j] * data.X[i, k]
                count += 1
            end
            if fit_intercept
                hess[count] += cst * data.X[i, j]
                count += 1
            end
        end
        if fit_intercept
            hess[count] += cst
        end
    end
    return nothing
end

"""
Compute diagonal of Hessian for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p)

"""
function diaghess!(diagh::AbstractVector{T}, ω::AbstractVector{T},
                   data::LogitData{T}, fit_intercept::Bool=false) where T
    # Sanity check
    n = ndata(data)
    p = nfeatures(data)
    @assert length(ω) == length(diagh)
    invn = one(T) / n
    @inbounds for i in 1:n
        σz = expit(-data.y_pred[i] * data.y[i])
        cst = invn * σz * (one(T) - σz)
        for j in 1:p
            diagh[j] += cst * data.X[i, j] * data.X[i, j]
        end
        if fit_intercept
            diagh[end] += cst
        end
    end
    return nothing
end

"""
Compute Hessian-vector product of logistic loss for given vector
parameter `ω ∈  R^p`.

## Complexity
O(n x 2 x p)

"""
function hessvec!(hessvec::AbstractVector{T}, ω::AbstractVector{T},
                  vec::AbstractVector{T}, data::LogitData{T}, fit_intercept::Bool=false) where T
    p = nfeatures(data)
    n = ndata(data)
    invn = one(T) / n
    @inbounds for i in 1:n
        σz = expit(-data.y_pred[i] * data.y[i])
        cst = invn * σz * (T(1) - σz)
        acc = zero(T)
        @simd for j in 1:p
            acc += data.X[i, j] * vec[j]
        end
        pσ = cst * acc
        for j in 1:p
            hessvec[j] += pσ * data.X[i, j]
        end

        if fit_intercept
            for j in 1:p
                hessvec[j] += cst * data.X[i, j] * vec[end]
            end
            hessvec[end] += pσ + cst * vec[end]
        end
    end
    return nothing
end


text:
n = rand(1:10)
m = n + rand(0:5)
A = rand(m, n)

r = polarfact(A, alg = :svd);

# Test unitary matrix U

U = r.U
H = r.H

@test_approx_eq_eps U'*U eye(n) 1e-7

# Test Hermitian positive semifefinite matrix H

@test issym(H)

for i in eigvals(H)
    @test i >= 0.
end

@test_approx_eq A U*H

println("SVD method passed test...")


text:
using AutoHashEquals

"""
Abstract type representing a pattern used in all the various pattern matching backends.
You can use the `Pattern` constructor to recursively convert an `Expr` (or any type satisfying [`Metatheory.TermInterface`](@ref)) to a [`Pattern`](@ref).
"""
abstract type Pattern end

import Base.==
==(a::Pattern, b::Pattern) = false
TermInterface.arity(p::Pattern) = 0

"""
Pattern variables will first match on any subterm
and instantiate the substitution to that subterm.
"""
mutable struct PatVar <: Pattern
    name::Symbol
    idx::Int
end
function ==(a::PatVar, b::PatVar)
    # (a.name == b.name)
    a.idx == b.idx
end
PatVar(var) = PatVar(var, -1)

"""
A pattern literal will match only against an instance of itself.
Example:
```julia
PatLiteral(2)
```
Will match only against values that are equal (using `Base.(==)`) to 2.

```julia
PatLiteral(:a)
```
Will match only against instances of the literal symbol `:a`.
"""
@auto_hash_equals struct PatLiteral{T} <: Pattern
    val::T
end

"""
Type assertions on a [`PatVar`](@ref), will match if and only if
the type of the matched term for the pattern variable `var` is a subtype
of `type`.
"""
@auto_hash_equals struct PatTypeAssertion <: Pattern
    var::PatVar
    type::Type
end

@auto_hash_equals struct PatSplatVar <: Pattern
    var::PatVar
end

"""
This type of pattern will match if and only if
the two subpatterns exist in the same equivalence class,
in the e-graph on which the matching is performed.
**Can be used only in the e-graphs backend**
"""
@auto_hash_equals struct PatEquiv <: Pattern
    left::Pattern
    right::Pattern
end

"""
Term patterns will match
on terms of the same `arity` and with the same
function symbol (`head`).
"""
@auto_hash_equals struct PatTerm <: Pattern
    head::Any
    args::Vector{Pattern}
    metadata::NamedTuple
end
TermInterface.arity(p::PatTerm) = length(p.args)
PatTerm(head, args) = PatTerm(head, args, (;))

"""
This pattern type matches on a function application
but instead of strictly matching on a head symbol,
it has a pattern variable as head. It can be used for
example to match arbitrary function calls.
"""
@auto_hash_equals struct PatAllTerm <: Pattern
    head::PatVar
    args::Vector{Pattern}
    metadata::NamedTuple
end
TermInterface.arity(p::PatAllTerm) = length(p.args)
PatAllTerm(head, args) = PatAllTerm(head, args, (;))

"""
Collects pattern variables appearing in a pattern into a vector of symbols
"""
patvars(p::PatLiteral, s) = s
patvars(p::PatVar, s) = push!(s, p.name)
patvars(p::PatTypeAssertion, s) = patvars(p.var, s)
patvars(p::PatSplatVar, s) = patvars(p.var, s)
function patvars(p::PatEquiv, s)
    patvars(p.left, s)
    patvars(p.right, s)
end

function patvars(p::PatTerm, s)
    for x ∈ p.args
        patvars(x, s)
    end
    return s
end

function patvars(p::PatAllTerm, s)
    push!(s, p.head.name)
    for x ∈ p.args
        patvars(x, s)
    end
    return s
end

function patvars(p::Pattern)
    unique(patvars(p, Symbol[]))
end

setindex!(p::PatLiteral, pvars) = nothing
function setindex!(p::PatVar, pvars)
    p.idx = findfirst((==)(p.name), pvars)
end
setindex!(p::PatTypeAssertion, pvars) = setindex!(p.var, pvars)
setindex!(p::PatSplatVar, pvars) = setindex!(p.var, pvars)

function setindex!(p::PatEquiv, pvars)
    setindex!(p.left, pvars)
    setindex!(p.right, pvars)
end

function setindex!(p::PatTerm, pvars)
    for x ∈ p.args
        setindex!(x, pvars)
    end
end

function setindex!(p::PatAllTerm, pvars)
    setindex!(p.head, pvars)
    for x ∈ p.args
        setindex!(x, pvars)
    end
end


text:
taskid(task::Task = current_task()) = objectid(task)

const Recalls = try
    Base.require(Base.PkgId(Base.UUID(0x30af7cf3eb4344c7afa733725b72a81e), "Recalls"))
catch
    nothing
end

istracing() = false

macro trace(args...)
    Recalls === nothing && return nothing
    note = :($Recalls.@note(taskid = $(taskid)(), $(args...),))
    @assert note.head === :macrocall
    @assert note.args[2] isa LineNumberNode
    note.args[2] = __source__
    quote
        if $istracing()
            $note
        end
        nothing
    end |> esc
end
# TODO: Use Preferences.jl to controll if Recalls.jl should be loaded? Since
# `args` may affect what variables are captured in closures, the existence of
# Recalls.jl changes the program slightly.

function enable_tracing()
    if Recalls === nothing
        error(
            "Tracing requires Recalls.jl to be installed in the same environment.",
            " See: https://github.com/tkf/Recalls.jl",
        )
    end
    prev = istracing()
    @eval istracing() = true
    return prev
end

function disable_tracing()
    prev = istracing()
    @eval istracing() = false
    return prev
end


text:
module Architectures
using Base.Iterators
using MLDataUtils

export MLP, batcheval

include("mlp.jl")

"""
    batcheval(model, bsize, X...)

Evaluate `model` on inputs `X` using batches of size `bsize` allocating the output automatically.
"""
function batcheval(model, bsize, X...)
    n = nobs(X[1])
    ranges = partition(1:n, bsize)
    idxs, state = iterate(ranges)
    Y1 = model((X_[:, idxs] for X_ in X)...).data # old Flux
    # Y1 = model((X_[:, idxs] for X_ in X)...) # Zygote Flux
    Y = similar(Y1, size(Y1)[1:end-1]..., n)
    Y[:, idxs] .= Y1
    while true
        res = iterate(ranges, state)
        isnothing(res) && break
        idxs, state = res
        Y[:,idxs] .= model((X_[:, idxs] for X_ in X)...).data # old Flux
        # Y[:,idxs] .= model((X_[:, idxs] for X_ in X)...) # Zygote Flux
    end
    Y
end

end

max_stars_repo_name: dhonza/FluxGoodies

max_issues_repo_name: dhonza/FluxGoodies

max_forks_repo_name: dhonza/FluxGoodies


text:
module EdgeDetect

using ForwardDiff
using Plots
using Statistics

function approx_diff(f, x, delta)
    h=1e-6
    return (f(x+delta*h) - f(x-delta*h))/(2*h)
end

function find(f,  N=110; a = -1.0, b = 1.0, deriv=:forwarddiff)

#    println("a")

    discont_f0, stddev_f0 = find_edges(f, N, a=a, b = b, deriv=deriv)

#    println("EDGES F0 ", discont_f0)

    begin
        if deriv == :approx
            df = x -> approx_diff(f,x, b-a)
        else
#            println("ForwardDiff")
            df = x -> ForwardDiff.derivative(f, x);
        end

        discont_f1, stddev_f1 = find_edges(df, N, a=a, b = b, deriv=deriv)
#        println("EDGES F1 ", discont_f1)
    end

    edges = [discont_f0;discont_f1]

    if size(edges)[1] == 0
        return zeros(0,2), Bool[]
    end

    #
    #println("EDGES ", edges)

    inds = sortperm(edges[:,1])
    edges = edges[inds, :]

    #println("EDGES2 ", edges)

    good = [1]
    for i = 2:size(edges)[1]
        if abs(edges[i-1,1] - edges[i,1]) > 1e-3 *(b-a)
            push!(good,i)
        end
    end
    edges = edges[good,:]

    contin = Bool[]
    for i = 1:size(edges)[1]

#        println("contin ", abs(f(edges[i,1]) - f(edges[i,2]) ), " std ", stddev_f0, " r ", abs(f(edges[i,1]) - f(edges[i,2]) ) / stddev_f0 , " r2 ", f(edges[i,1]) / stddev_f0)

        if  abs(f(edges[i,1]) - f(edges[i,2]) ) / stddev_f0 < 0.5e-2 && f(edges[i,1]) / stddev_f0 < 1e3
            push!(contin, true)
        else
            push!(contin, false)
        end
    end
    return edges, contin

end

function find_edges(f, N = 110; a=-1.0, b = 1.0, deriv=:forwarddiff)

    if deriv == :approx
        df = x -> approx_diff(f,x, b-a)
    else
        df = x -> ForwardDiff.derivative(f, x);
    end

    #    f(0.1)
#    df(0.1)
    grid, grid2, fg, dfg, df_approx, df_err = get_derivs(f,df, N, a=a, b=b)

    stddev = std(fg)

    local_max = []
    for i in 2:N-1
        if abs(df_err[i])+1e-10 > abs(df_err[i-1]) + 1e-10 && abs(df_err[i]) > abs(df_err[i+1]) + 1e-10
            push!(local_max, i)
        end
    end

#    println("local_max ", local_max, grid[local_max])

    true_edge = zeros(0,2)

    for i in local_max
#        println([i, grid[i], grid[i+1], fg[i], fg[i+1], df_err[i] ])

        gridX, grid2X, fgX, dfgX, df_approxX, df_errX = get_derivs(f,df, 20, a=grid[i], b=grid[i+1])
        ind =argmax(abs.(df_errX))
        new_max_err = abs.(df_errX[ind])
        old_max_err = abs(df_err[i])

        if new_max_err >= old_max_err *0.1 || new_max_err > 1e2 #true discont doesn't shrink (much) as we get closer

            for j = 1:7
                gridX, grid2X, fgX, dfgX, df_approxX, df_errX = get_derivs(f,df, 20, a=gridX[ind], b=gridX[ind+1])
                ind =argmax(abs.(df_errX))
                new_max_err = abs.(df_errX[ind])
#                println("iter $j ", [gridX[ind], gridX[ind+1]], " " , gridX[ind+1] - gridX[ind], " " , [new_max_err, old_max_err])
                old_max_err = new_max_err
            end

            if new_max_err >= old_max_err *0.1 && new_max_err > 1e3
                true_edge = vcat(true_edge, [gridX[ind]  gridX[ind+1]])
#                println("TRUE dis ", [gridX[ind], gridX[ind+1]])
#                println()
            end

        end

    end

#    println("TRUE EDGE ", true_edge)

#    plot(grid2[2:end-1], df_err[2:end-1], label="df err ")
#    plot!(grid2, dfg, label="FD", linewidth=3)
#    plot!(grid2, df_approx, label="df appprox", linewidth=2)
#    plot!(grid, fg, label="f")

    return true_edge, stddev

end

function get_derivs(f, df, N; a=-1, b = 1)

    begin

        aa = a + (b-a)/N * 0.001
        bb = b - (b-a)/N * 0.0015

        h = (bb - aa)/(N) * (1.0 - 10^-12)
        grid = (aa):h:(bb)
        grid2 = (aa:h:bb)[1:end-1] .+ h/2
    end

    #    N = length(grid)-2

    #println("grid ")

    fg = f.(grid)
    dfg = df.(grid2)

    #fg = zeros(length(grid))
    #for i = 1:length(grid)
    #    fg[i] = f(grid[i])
    #end
    #dfg = zeros(length(grid2))
    #for i = 1:length(grid2)
    #    dfg[i] = df(grid2[i])
    #end

#    begin
#        fg = zeros(size(grid))
#        dfg = zeros(size(grid2))
#        for i = 1:length(grid)
#        fg[i] = f(grid[i])
#        end
#        for i = 1:length(grid2)
#            dfg[i] = df(grid2[i])
#        end
#    end

#    println("size fg ", size(fg))
#    println("N $N")

    begin
        df_approx = zeros(N)
        for i in 1:(N)

            df_approx[i] = (fg[i+1] - fg[i]) / h

        end
        df_err = df_approx - dfg
    end

    return grid, grid2, fg, dfg, df_approx, df_err

end

end #end module


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

__precompile__(true)

module IntervalRootFinding

using IntervalArithmetic
using ForwardDiff
using StaticArrays

using LinearAlgebra: I, Diagonal

import Base: ⊆, show, big, \

import Polynomials: roots

## Root finding
export
    derivative, jacobian,  # reexport derivative from ForwardDiff
    Root, is_unique,
    roots, find_roots,
    bisect, newton1d, quadratic_roots,
    gauss_seidel_interval, gauss_seidel_interval!,
    gauss_seidel_contractor, gauss_seidel_contractor!,
    gauss_elimination_interval, gauss_elimination_interval!,
    slope

export isunique, root_status

import IntervalArithmetic: interval, wideinterval

const derivative = ForwardDiff.derivative
const D = derivative

const where_bisect = IntervalArithmetic.where_bisect ## 127//256

include("root_object.jl")

include("newton.jl")
include("krawczyk.jl")

include("complex.jl")
include("contractors.jl")
include("roots.jl")
include("newton1d.jl")
include("quadratic.jl")
include("linear_eq.jl")
include("slopes.jl")

gradient(f) = X -> ForwardDiff.gradient(f, X[:])

ForwardDiff.jacobian(f, X::IntervalBox) = ForwardDiff.jacobian(f, X.v)

const ∇ = gradient

export ∇

function find_roots(f::Function, a::Interval{T}, method::Function = newton;
                    tolerance = eps(T), debug = false, maxlevel = 30) where {T}

    method(f, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
end

function find_roots(f::Function, f_prime::Function, a::Interval{T}, method::Function=newton;
                    tolerance=eps(T), debug=false, maxlevel=30) where {T}

    method(f, f_prime, a; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
end

function find_roots(f::Function, a::Real, b::Real, method::Function=newton;
           tolerance=eps(1.0*a), debug=false, maxlevel=30, precision::Int=-1)

    if precision >= 0
        setprecision(Interval, precision) do
            find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)
        end

    else  # use current precision

        find_roots(f, @interval(a, b), method; tolerance=tolerance, debug=debug, maxlevel=maxlevel)

    end
end

function find_roots_midpoint(f::Function, a::Real, b::Real, method::Function=newton;
           tolerance=eps(1.0*a), debug=false, maxlevel=30, precision=-1)

    roots = find_roots(f, a, b, method; tolerance=tolerance, debug=debug, maxlevel=maxlevel, precision=precision)

    T = eltype(roots[1].interval)

    midpoints = T[]
    radii = T[]
    root_symbols = Symbol[]  # :unique or :unknown

    if length(roots) == 0
        return (midpoints, radii, root_symbols)  # still empty
    end

    for root in roots
        midpoint, radius = midpoint_radius(root.interval)
        push!(midpoints, midpoint)
        push!(radii, radius)

        push!(root_symbols, root.status)

    end

    (midpoints, radii, root_symbols)

end

function clean_roots(f, roots)

    # order, remove duplicates, and include intervals X only if f(X) contains 0
    sort!(roots, lt=lexless)
    roots = unique(roots)
    roots = filter(x -> 0 ∈ f(x.interval), roots)

    # merge neighbouring roots if they touch:

    if length(roots) < 2
        return roots
    end

    new_roots = eltype(roots)[]

    base_root = roots[1]

    for i in 2:length(roots)
        current_root = roots[i]

        if isempty(base_root.interval ∩ current_root.interval) ||
                (base_root.status != current_root.status)

            push!(new_roots, base_root)
            base_root = current_root
        else
            base_root = Root(hull(base_root.interval, current_root.interval), base_root.status)
        end
    end

    push!(new_roots, base_root)

    return new_roots

end

end