text:
########################################################################################
# Macro to call appropriate 3D method given a 2D one
########################################################################################

macro two_points_on_the_surface(fn)
	callermodule = __module__
	quote
		function $(esc(fn))(comm, nobs1::Point2D, nobs2::Point2D,
			los::los_direction = los_radial(); kwargs...)

			r_obs = get(kwargs, :r_obs, r_obs_default)
			xobs1 = Point3D(r_obs, nobs1)
			xobs2 = Point3D(r_obs, nobs2)
			$callermodule.$fn(comm, xobs1, xobs2, los; kwargs...)
		end
	end
end

########################################################################################
# Get the modes to be used
########################################################################################

function ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)
	ℓ_range = get(kwargs, :ℓ_range, ℓ_arr)
	ν_ind_range = get(kwargs, :ν_ind_range, axes(ν_arr, 1))
	ℓ_range, ν_ind_range
end

#######################################################################################################
# Full Frequency axis
#######################################################################################################

function pad_zeros_ν(arr::Array, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros, dim::Integer = 1)
 	ax_leading  = size(arr)[1:dim-1]
	ax_trailing = size(arr)[dim+1:end]
	N_leading_zeros_ν = ν_start_zeros + first(ν_ind_range) - 1
	N_trailing_zeros_ν = Nν_Gfn - last(ν_ind_range) + ν_end_zeros

	inds_leading = (ax_leading..., N_leading_zeros_ν, ax_trailing...) :: NTuple{ndims(arr),Int}
	inds_trailing = (ax_leading..., N_trailing_zeros_ν , ax_trailing...) :: NTuple{ndims(arr),Int}

	lead_arr = zeros(eltype(arr), inds_leading)
	trail_arr = zeros(eltype(arr), inds_trailing)

	cat(lead_arr, arr, trail_arr, dims = dim) :: typeof(arr)
end

# If an OffsetArray is passed we may infer the frequency range from its axes
function pad_zeros_ν(arr::OffsetArray, Nν_Gfn,
	ν_start_zeros, ν_end_zeros, dim::Integer = 1)

	pad_zeros_ν(parent(arr), axes(arr, dim), Nν_Gfn, ν_start_zeros, ν_end_zeros, dim)
end

########################################################################################
# Read parameters for source and observation points
########################################################################################

read_parameters_for_points(; kwargs...) = read_all_parameters(; kwargs...)

function read_parameters_for_points(xobs1::SphericalPoint, xobs2::SphericalPoint; kwargs...)
	p_Gsrc = read_all_parameters(; kwargs...)
	p_Gobs1 = read_all_parameters(; kwargs..., r_src = radius(xobs1))
	p_Gobs2 = read_all_parameters(; kwargs..., r_src = radius(xobs2))
	return p_Gsrc, p_Gobs1, p_Gobs2
end
function read_parameters_for_points((xobs1, xobs1′)::NTuple{2,SphericalPoint},
	(xobs2, xobs2′)::NTuple{2,SphericalPoint}; kwargs...)
	p_Gsrc = read_all_parameters(; kwargs...)
	p_Gobs1 = read_all_parameters(; kwargs..., r_src = radius(xobs1))
	p_Gobs2 = read_all_parameters(; kwargs..., r_src = radius(xobs2))
	return p_Gsrc, p_Gobs1, p_Gobs2
end

function read_parameters_for_points(xobs1::Point3D, xobs2_arr::Vector{<:Point3D}; kwargs...)
	p_Gsrc = read_all_parameters(; kwargs...)
	p_Gobs1 = read_all_parameters(; kwargs..., r_src = radius(xobs1))
	p_Gobs2 = [read_all_parameters(; kwargs..., r_src = radius(xobs2)) for xobs2 in xobs2_arr]
	return p_Gsrc, p_Gobs1, p_Gobs2
end

function read_parameters_for_points(::Point2D, ::Vector{<:Point2D}; kwargs...)
	p_Gsrc = read_all_parameters(; kwargs...)
	r_obs = get(kwargs, :r_obs, r_obs_default)
	p_Gobs = read_all_parameters(; kwargs..., r_src = r_obs)
	return p_Gsrc, p_Gobs, p_Gobs
end

_extremaelementrange(ℓ_ωind_iter_on_proc; dims) = UnitRange(extremaelement(ℓ_ωind_iter_on_proc, dims = dims)...)

########################################################################################
# Bipolar harmonics
########################################################################################

function computeY₀₀(::los_radial, xobs1::SphericalPoint, xobs2::SphericalPoint, ℓ_range::AbstractRange{<:Integer})
	v = collectPl(cosχ(xobs1, xobs2), lmax = maximum(ℓ_range))[ℓ_range]
	norm(ℓ) = (2ℓ+1)/4π
	@. v *= norm(ℓ_range)
	return OffsetArray(v, ℓ_range)
end

function computeY₀₀(::los_earth, xobs1::SphericalPoint, xobs2::SphericalPoint, ℓ_range::AbstractRange{<:Integer})
	B = BipolarSphericalHarmonics.monopolarharmonics(GSH(), Point2D(xobs1)..., Point2D(xobs2)..., maximum(ℓ_range))
	v = BipolarSphericalHarmonics.eltypeY(B)[biposh(B, Point2D(xobs1)..., Point2D(xobs2)..., 0, 0, ℓ, ℓ) for ℓ in ℓ_range]
	norm(ℓ) = (-1)^ℓ * √(2ℓ+1)
	@. v *= norm(ℓ_range)
	OffsetArray(v, ℓ_range)
end

function computeY₁₀(::los_radial, xobs1::SphericalPoint, xobs2::SphericalPoint, ℓ_range::AbstractRange{<:Integer})
	v = epsilon.(collectPl(cosχ(xobs1, Point2D(xobs2.θ, Dual(xobs2.ϕ, 1))), lmax = maximum(ℓ_range)))[ℓ_range]
	norm(ℓ) = (2ℓ+1)/4π
	@. v *= norm(ℓ_range)
	return OffsetArray(v, ℓ_range)
end

function computeY₁₀(::los_earth, xobs1::SphericalPoint, xobs2::SphericalPoint, ℓ_range::AbstractRange{<:Integer})
	B = BipolarSphericalHarmonics.monopolarharmonics(GSH(), Point2D(xobs1)..., Point2D(xobs2)..., maximum(ℓ_range))
	v = BipolarSphericalHarmonics.eltypeY(B)[biposh(B, Point2D(xobs1)..., Point2D(xobs2)..., 1, 0, ℓ, ℓ) for ℓ in ℓ_range]
	norm(ℓ) = (-1)^ℓ * √(2ℓ+1)
	@. v *= norm(ℓ_range)
	OffsetArray(v, ℓ_range)
end

function _multiplylos!(Y12, j1j2ind, Y12_j1j2::AbstractVector, l1l2::AbstractVector)
	for (lmind, Y12_j1j2_lm) in pairs(Y12_j1j2)
		Y12_j1j2[lmind] = Y12_j1j2_lm .* l1l2
	end
	nothing
end
function _multiplylos!(Y12, j1j2ind, Y12_j1j2::StaticVector, l1l2::AbstractVector)
	Y12[j1j2ind] = Y12_j1j2 .* l1l2
	nothing
end
los_projected_biposh!(::los_radial, Y1, l1l2) = Y1
los_projected_biposh!(::los_earth, Y12::AbstractVector{<:Number}, l1l2) = Y12 .* l1l2
function los_projected_biposh!(::los_earth, Y12, l1l2)
	for (j1j2ind, Y12_j1j2) in pairs(Y12)
		Y12_j1j2 === nothing && continue
		_multiplylos!(Y12, j1j2ind, Y12_j1j2, l1l2)
	end
	return Y12
end
function los_projected_biposh!(los, Y12, l1, l2)
	l1l2 = LinearAlgebra.kron(parent(l1), parent(l2))
	los_projected_biposh!(los, Y12, l1l2)
end

_biposh_returntype(v::AbstractArray) =  _biposh_returntype(typeof(v))
_biposh_returntype(::Type{Vector{T}}) where {T} = Vector{_biposh_returntype(T)}
_biposh_returntype(::Type{SHVector{T, Vector{T}, M}}) where {T,M} = (T2 = _biposh_returntype(T); SHVector{T2, Vector{T2}, M})
_biposh_returntype(::Type{Union{Nothing, T}}) where {T} = Union{Nothing, _biposh_returntype(T)}
_biposh_returntype(::Type{T}) where {T<:Number} = T
_biposh_returntype(::Type{SVector{9,T}}) where {T<:Number} = OffsetMatrix{T,SMatrix{2,2,T,4}}

biposh_spheroidal(Y::OffsetVector) = OffsetArray(biposh_spheroidal(parent(Y)), axes(Y))
function biposh_spheroidal(Y::SHVector)
	SHArray(biposh_spheroidal(parent(Y)), SphericalHarmonicArrays.modes(Y))
end
biposh_spheroidal(Y::Vector) = eltype(_biposh_returntype(Y))[biposh_spheroidal(Y[i]) for i in eachindex(Y)]
biposh_spheroidal(Y::Number) = Y
biposh_spheroidal(::Nothing) = nothing
function biposh_spheroidal(Y::SVector{9,<:Number})
	T00 = Y[kronindex(GSH(), 0, 0)]
	T10 = Y[kronindex(GSH(), 1, 0)] + Y[kronindex(GSH(), -1, 0)]
	T01 = Y[kronindex(GSH(), 0, 1)] + Y[kronindex(GSH(), 0, -1)]
	T11 = Y[kronindex(GSH(), -1, -1)] + Y[kronindex(GSH(), 1, -1)] +
			Y[kronindex(GSH(), -1, 1)] + Y[kronindex(GSH(), 1, 1)]
	OffsetArray(SMatrix{2,2}(T00, T10, T01, T11), 0:1, 0:1)
end

function los_projected_biposh_spheroidal(computeY, xobs1, xobs2::SphericalPoint, los, ℓ_range::AbstractUnitRange)
	Y12 = computeY(los, xobs1, xobs2, ℓ_range)
	los_projected_biposh_spheroidal(Y12, xobs1, xobs2, los)
end
function los_projected_biposh_spheroidal(Y12, xobs1, xobs2::SphericalPoint, los)
	los_projected_biposh_spheroidal!(deepcopy(Y12), xobs1, xobs2, los)
end
function los_projected_biposh_spheroidal!(Y12, xobs1, xobs2::SphericalPoint, los)
	l1, l2 = line_of_sight_covariant.((xobs1, xobs2), los)
	l1l2Y12 = los_projected_biposh!(los, Y12, l1, l2)
	biposh_spheroidal(l1l2Y12)
end
function los_projected_biposh_spheroidal_Y₀₀(xobs1::SphericalPoint, xobs2::SphericalPoint, xobs1′::SphericalPoint, xobs2′::SphericalPoint, los::los_direction, ℓ_range)
	Y12 = computeY₀₀(los, xobs1, xobs2, ℓ_range)
	los_projected_biposh_spheroidal_Y₀₀(Y12, xobs1, xobs2, xobs1′, xobs2′, los)
end
function los_projected_biposh_spheroidal_Y₀₀(Y12, xobs1::SphericalPoint, xobs2::SphericalPoint, xobs1′::SphericalPoint, xobs2′::SphericalPoint, los::los_direction)
	l1, l2 = line_of_sight_covariant.((xobs1′, xobs2′), los)
	M′ = rotation_points(HelicityCovariant(), xobs1, xobs2, xobs1′, xobs2′)
	l1l2 = M′ * LinearAlgebra.kron(l1, l2)
	l1l2Y12 = los_projected_biposh!(los, deepcopy(Y12), l1l2)
	biposh_spheroidal(l1l2Y12)
end
function los_projected_biposh_spheroidal(computeY, nobs1, nobs2_arr::Vector{<:SphericalPoint}, los, ℓ_range::AbstractUnitRange)
	l1 = line_of_sight_covariant(nobs1, los)
	l2arr = line_of_sight_covariant.(nobs2_arr, los)
	_Y12 = computeY.(los, nobs1, nobs2_arr, (ℓ_range,))
	OffsetArray(permutedims(reduce(hcat,
	biposh_spheroidal.(los_projected_biposh!.(los, _Y12, (l1,), l2arr)))), :, ℓ_range)
end

function los_projected_biposh_spheroidal(computeY, xobs1, xobs2, los, ℓ_ωind_iter_on_proc::ProductSplit)
	ℓ_range = _extremaelementrange(ℓ_ωind_iter_on_proc, dims = 1)
	los_projected_biposh_spheroidal(computeY, xobs1, xobs2, los, ℓ_range)
end

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

function allocateGfn(los::los_direction, obs_at_same_height::Bool)
	α_r₁ = zeros(ComplexF64, obsindG(los))
	α_r₂ = obs_at_same_height ? α_r₁ : zeros(ComplexF64, obsindG(los))
	(; α_r₁, α_r₂)
end

function C_FITS_header(xobs1::SphericalPoint, xobs2::SphericalPoint)
	FITSHeader(["X1R", "X1THETA", "X1PHI", "X2R", "X2THETA", "X2PHI"],
		Any[float(radius(xobs1)), float(xobs1.θ), float(xobs1.ϕ),
		float(radius(xobs2)), float(xobs2.θ), float(xobs2.ϕ)],
		["Radius of the first observation point",
		"Co-latitude of the first observation point",
		"Azimuth of the first observation point",
		"Radius of the second observation point",
		"Co-latitude of the second observation point",
		"Azimuth of the second observation point"])
end

function C_FITS_header((n1, n1′)::NTuple{2,SphericalPoint}, (n2, n2′)::NTuple{2,SphericalPoint})
	xobs1 = Point3D(r_obs_default, n1)
	xobs2 = Point3D(r_obs_default, n2)
	f = C_FITS_header(xobs1, xobs2)
	for (k, v, c) in zip(["X1PTHETA", "X1PPHI", "X2PTHETA", "X2PPHI"], [n1′..., n2′...],
		["Co-latitude of the first observation point rotated",
		"Azimuth of the first observation point rotated",
		"Co-latitude of the second observation point rotated",
		"Azimuth of the second observation point rotated"
		])
		f[k] = v
		set_comment!(f, k, c)
	end
end

function C_FITS_header(n1::Point2D, n2::Vector{<:Point2D})
	header = C_FITS_header(n1, first(n2))

	θmin = minimum(x->float(x.θ), n2)
	header["THMIN"] = θmin
	θmax = maximum(x->float(x.θ), n2)
	header["THMAX"] = θmax

	ϕmin = minimum(x->float(x.ϕ), n2)
	header["PHIMIN"] = ϕmin
	ϕmax = maximum(x->float(x.ϕ), n2)
	header["PHIMAX"] = ϕmax

	header["NTHPHI"] = length(n2)

	header
end

########################################################################################
# cross-covariances
########################################################################################

function Cωℓ(::los_radial, ω, ℓ, α_r₁::AbstractArray{<:Any, 0}, α_r₂::AbstractArray{<:Any, 0}, Pl)
	ω^2 * Powspec(ω) * conj(α_r₁[]) * α_r₂[] * Pl
end

function Cωℓ(::los_earth, ω, ℓ, α_r₁::AbstractVector{<:Complex}, α_r₂::AbstractVector{<:Complex}, Y12)
	pre = ω^2 * Powspec(ω)
	s = zero(ComplexF64)
	for α in 0:1, β in 0:1
		s += pre * conj(α_r₁[α]) * α_r₂[β] * Y12[α, β]
	end
	s
end

function Cω_partial(localtimer, ℓ_ωind_iter_on_proc,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src, (ℓ_arr, 1:Nν_Gfn),
		ℓ_ωind_iter_on_proc, NGfn_files_src)

	ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)[end]

	Cω_proc = zeros(ComplexF64, ν_ind_range)

	@timeit localtimer "biposh" begin
		Y12 = los_projected_biposh_spheroidal(computeY₀₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)
	end

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	@unpack α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]

		@timeit localtimer "FITS" begin
			read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

			if r₁_ind != r₂_ind
				read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₂_ind, 1)
			end
		end

		@timeit localtimer "Cω calculation" begin
			Cω_proc[ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y12[ℓ])
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	parent(Cω_proc)
end
function Cω_partial(localtimer, ℓ_ωind_iter_on_proc,
	(xobs1, xobs1′)::NTuple{2,SphericalPoint},
	(xobs2, xobs2′)::NTuple{2,SphericalPoint}, los::los_direction,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src, (ℓ_arr, 1:Nν_Gfn),
		ℓ_ωind_iter_on_proc, NGfn_files_src)

	ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)[end]

	Cω_proc = zeros(ComplexF64, 2, ν_ind_range)

	ℓ_range = _extremaelementrange(ℓ_ωind_iter_on_proc, dims = 1)
	Y12_Helicity = computeY₀₀(los, xobs1, xobs2, ℓ_range)
	Y12 = los_projected_biposh_spheroidal(Y12_Helicity, xobs1, xobs2, los)
	Y1′2′ = los_projected_biposh_spheroidal_Y₀₀(Y12_Helicity, xobs1, xobs2, xobs1′, xobs2′, los)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	@unpack α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]

		@timeit localtimer "FITS" begin
			read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

			if r₁_ind != r₂_ind
				read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₂_ind, 1)
			end
		end

		@timeit localtimer "Cω calculation" begin
			Cω_proc[1,ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y12[ℓ])
			Cω_proc[2,ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y1′2′[ℓ])
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	parent(permutedims(Cω_proc))
end

function Cω_partial(localtimer, ℓ_ωind_iter_on_proc,
	nobs1::Point2D, nobs2_arr::Vector{<:Point2D}, los::los_direction,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)[end]

	r₁_ind = radial_grid_index(nobs1)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
		(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	@unpack α_r₁ = allocateGfn(los, true)

	Cω_proc = zeros(ComplexF64, length(nobs2_arr), ν_ind_range)

	Y12 = los_projected_biposh_spheroidal(computeY₀₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc
		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

		for n2ind in axes(Cω_proc, 1)
			Cω_proc[n2ind, ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₁, Y12[n2ind, ℓ])
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	permutedims(parent(Cω_proc))
end

########################################################################################################
# Derivatives of cross-covariance (useful in the radial case)
########################################################################################################

function ∂ϕ₂Cω_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	Cω_proc = zeros(ComplexF64, ν_ind_range)

	Y12 = los_projected_biposh_spheroidal(computeY₁₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	@unpack α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

		if r₁_ind != r₂_ind
			read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₂_ind, 1)
		end

		Cω_proc[ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y12[ℓ])
	end
	closeGfnfits(Gfn_fits_files_src)
	parent(Cω_proc)
end

function ∂ϕ₂Cω_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	nobs1::Point2D, nobs2_arr::Vector{<:Point2D}, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, Nν_Gfn = p_Gsrc

	ℓ_arr, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	Cω_proc = zeros(ComplexF64, ν_ind_range)

	Y12 = los_projected_biposh_spheroidal(computeY₁₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)

	r_obs_ind = radial_grid_index(nobs1)

	@unpack α_r₁, α_r₂ = allocateGfn(los, true)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r_obs_ind, 1)

		for n2ind in eachindex(nobs2_arr)
			Cω_proc[n2ind, ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y12[n2ind, ℓ])
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	permutedims(parent(Cω_proc))
end

function Cω_∂ϕ₂Cω_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	Cω_proc = zeros(ComplexF64, 0:1, ν_ind_range)

	Pl_cosχ = los_projected_biposh_spheroidal(computeY₀₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)
	∂ϕ₂Pl_cosχ = los_projected_biposh_spheroidal(computeY₁₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	@unpack α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)
		if r₁_ind != r₂_ind
			read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₂_ind, 1)
		end

		f = Cωℓ(los, ω, ℓ, α_r₁, α_r₂, 1)

		Cω_proc[0, ω_ind] += f * Pl_cosχ[ℓ]
		Cω_proc[1, ω_ind] += f * ∂ϕ₂Pl_cosχ[ℓ]
	end
	closeGfnfits(Gfn_fits_files_src)
	permutedims(parent(Cω_proc))
end

function Cω_∂ϕ₂Cω_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	nobs1::Point2D, nobs2_arr::Vector{<:Point2D}, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), ℓ_ωind_iter_on_proc, NGfn_files_src)

	Cω_proc = zeros(ComplexF64, 0:1, length(nobs2_arr), ν_ind_range)

	r_obs_ind = radial_grid_index(nobs1)

	lmax = maximum(ℓ_arr)

	∂ϕ₂Pl_cosχ = zeros(0:lmax, length(nobs2_arr))
	Pl_cosχ = zeros(0:lmax, length(nobs2_arr))

	Pl_cosχ = los_projected_biposh_spheroidal(computeY₀₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)
	∂ϕ₂Pl_cosχ = los_projected_biposh_spheroidal(computeY₁₀, nobs1, nobs2_arr, los, ℓ_ωind_iter_on_proc)

	@unpack α_r₁ = allocateGfn(los, true)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r_obs_ind, 1)

		f = Cωℓ(los_radial(), ω, ℓ, α_r₁, α_r₁, 1)

		for n2ind in axes(Cω_proc, 2)
			Cω_proc[0, n2ind, ω_ind] += f * Pl_cosχ[n2ind, ℓ]
			Cω_proc[1, n2ind, ω_ind] += f * ∂ϕ₂Pl_cosχ[n2ind, ℓ]
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	permutedims(parent(Cω_proc),[3, 2, 1])
end

########################################################################################
# Functions that iterate over the modes in parallel
########################################################################################

for f in (:Cω, :∂ϕ₂Cω, :Cω_∂ϕ₂Cω)
	f_str = String(f)
	f_feeder = Symbol(f, "_feeder")
	f_partial = Symbol(f, "_partial")
	@eval function $f_feeder(comm, args...; kwargs...)
		r_src, _, c_scale = read_rsrc_robs_c_scale(kwargs)
		p_Gsrc = read_all_parameters(r_src = r_src, c_scale = c_scale)
		@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros = p_Gsrc
		iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

		C = pmapsum(comm, $f_partial, iters, args..., p_Gsrc, r_src, c_scale)
		C === nothing && return nothing
		ν_ind_range = last(iters)
		return pad_zeros_ν(C, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros)
	end
	# With or without los, 3D points
	@eval function $f(comm, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
		C = $f_feeder(comm, xobs1, xobs2, los; kwargs...)
		C === nothing && return nothing

		lostag = los isa los_radial ? "" : "_los"
		filename = $f_str * lostag * ".fits"

		if get(kwargs, :save, false)
			filepath = joinpath(SCRATCH_kerneldir[], filename)
			FITSIO.fitswrite(filepath, reinterpret_as_float(C), header = C_FITS_header(xobs1, xobs2))
		end
		return C
	end
	# Multiple 2D points
	@eval function $f(comm, nobs1::Point2D, nobs2_arr::Vector{<:Point2D},
		los::los_direction = los_radial(); kwargs...)

		C = $f_feeder(comm, nobs1, nobs2_arr, los; kwargs...)
		C === nothing && return nothing
		save_to_fits_and_return($f_str * "_n2arr.fits", C, header = C_FITS_header(nobs1, nobs2_arr))
	end
end

########################################################################################################
# Spectrum of C(ℓ, ω)
########################################################################################################

function Cωℓ_spectrum_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	r_obs, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src)
	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, Nν_Gfn, ω_arr = p_Gsrc

	ℓ_range, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src, (ℓ_arr, 1:Nν_Gfn),
						ℓ_ωind_iter_on_proc, NGfn_files_src)

	@unpack α_r₁, = allocateGfn(los, true)
	r₁_ind = radial_grid_index(r_obs)

	Cℓω = zeros(ℓ_range, ν_ind_range)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		# get the only element of a 0-dim array
		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

		# m-averaged, so divided by 2ℓ+1
		Cℓω[ℓ, ω_ind] = ω^2 * Powspec(ω) * 1/4π * abs2(α_r₁[])
	end
	closeGfnfits(Gfn_fits_files_src)
	permutedims(parent(Cℓω))
end

function Cωℓ_spectrum(comm; kwargs...)
	r_src, r_obs, c_scale = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc = read_all_parameters(r_src = r_src, c_scale = c_scale)
	@unpack ℓ_arr, ν_arr = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	Cωℓ_in_range = pmapsum(comm, Cωℓ_spectrum_partial, iters,
		r_obs, los_radial(), p_Gsrc, r_src)

	Cωℓ_in_range === nothing && return nothing
	save_to_fits_and_return("Cωℓ_in_range.fits", Cωℓ_in_range)
end

########################################################################################################
# Time-domain cross-covariance
########################################################################################################

for (_fnω, _fnt) in ((:Cω, :Ct), (:∂ϕ₂Cω, :∂ϕ₂Ct))
	_fnts = String(_fnt)
	@eval function $_fnt(comm, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
		C = $_fnω(comm, xobs1, xobs2, los; kwargs...)
		C === nothing && return nothing

		C_t = $_fnt(C; kwargs...)
		lostag = los isa los_radial ? "" : "_los"
		vtag = xobs2 isa AbstractVector ? "_n2arr" : ""
		filename = $_fnts * vtag * lostag * ".fits"
		save_to_fits_and_return(filename, C_t, header = C_FITS_header(xobs1, xobs2))
	end
	@eval function $_fnt(Cω_arr::AbstractArray{<:Complex}; kwargs...)
		dν = get(kwargs, :dν) do
			read_parameters("dν"; kwargs...)[1]
		end

		τ_ind_range = get(kwargs, :τ_ind_range, Colon())

		C_t = fft_ω_to_t(Cω_arr, dν)

		if !(τ_ind_range isa Colon)
			C_t = C_t[τ_ind_range,..]
		end

		return C_t
	end
end

########################################################################################################
# Cross-covariance at all distances on the equator, essentially the time-distance diagram
########################################################################################################

function CΔϕω_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	r₁, r₂, los::los_radial, Δϕ_arr,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((r₁, r₂))

	nϕ = length(Δϕ_arr)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src, (ℓ_arr, 1:Nν_Gfn),
						ℓ_ωind_iter_on_proc, NGfn_files_src)

	Cϕω_arr = zeros(ComplexF64, nϕ, ν_ind_range)

	n1 = Point2D(pi/2, 0)
	n2_arr = Point2D.(pi/2, Δϕ_arr)
	Y12 = los_projected_biposh_spheroidal(computeY₀₀, n1, n2_arr, los, ℓ_ωind_iter_on_proc)

	@unpack α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₁_ind, 1)

		if r₁_ind != r₂_ind
			read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, obsindFITS(los), 1, r₂_ind, 1)
		end

		for ϕ_ind in 1:nϕ
			Cϕω_arr[ϕ_ind, ω_ind] += Cωℓ(los, ω, ℓ, α_r₁, α_r₂, Y12[ϕ_ind, ℓ])
		end
	end
	closeGfnfits(Gfn_fits_files_src)
	parent(Cϕω_arr)
end

function CΔϕω(comm, r₁ = r_obs_default, r₂ = r_obs_default,
	los::los_radial = los_radial(); kwargs...)

	r_src, _, c_scale = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros = p_Gsrc
	ℓ_range, ν_ind_range = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	Δϕ_arr = get(kwargs, :Δϕ_arr, LinRange(0, π, 100))

	Cω_in_range = pmapsum(comm, CΔϕω_partial, (ℓ_range, ν_ind_range),
			r₁, r₂, los, Δϕ_arr, p_Gsrc, r_src, c_scale)

	Cω_in_range === nothing && return nothing
	CΔϕω = pad_zeros_ν(Cω_in_range, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros, 2)
	if get(kwargs, :save, true)
		filename = joinpath(SCRATCH_kerneldir[], "CΔϕω.fits")
		FITSIO.fitswrite(filename, reinterpret_as_float(CΔϕω))
	end
	return CΔϕω
end

function CtΔϕ(comm, r₁ = r_obs_default, r₂ = r_obs_default,
	los::los_direction = los_radial(); kwargs...)

	dν = get(kwargs, :dν) do
		read_parameters("dν"; kwargs...)[1]
	end

	C = CΔϕω(comm, r₁, r₂, los; kwargs..., save = false)
	C === nothing && return nothing
	CωΔϕ = permutedims(C)
	CtΔϕ = fft_ω_to_t(CωΔϕ, dν)

	τ_ind_range = get(kwargs, :τ_ind_range, Colon())
	if !(τ_ind_range isa Colon)
		CtΔϕ = CtΔϕ[τ_ind_range, ..]
	end

	if get(kwargs, :save, true)
		tag = los isa los_radial ? "" : "_los"
		filename = joinpath(SCRATCH_kerneldir[], "C_t_phi$tag.fits")
		FITSIO.fitswrite(filename, CtΔϕ)
	end
	return CtΔϕ
end

function CΔϕt(comm, r₁ = r_obs_default, r₂ = r_obs_default,
	los::los_direction = los_radial(); kwargs...)

	_Ctϕ = CtΔϕ(comm, r₁, r₂, los; kwargs..., save = false)
	_Ctϕ === nothing && return nothing
	C = permutedims(_Ctϕ)

	if get(kwargs, :save, true)
		tag = los isa los_radial ? "" : "_los"
		filename = joinpath(SCRATCH_kerneldir[], "C_phi_t$tag.fits")
		FITSIO.fitswrite(filename, C)
	end
	return C
end

########################################################################################################
# Cross-covariance in a rotating frame
########################################################################################################

function Cτ_rotating_partial(localtimer, ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_radial,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, c_scale = 1)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = c_scale)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	ℓ_range, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src, (ℓ_arr, 1:Nν_Gfn),
						ℓ_ωind_iter_on_proc, NGfn_files_src)

	α₁conjα₂ωℓ = zeros(ComplexF64, ℓ_range, ν_ind_range)

	α_r₁, α_r₂ = allocateGfn(los, r₁_ind == r₂_ind)

	# Read all radial parts
	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		read_Gfn_file_at_index!(α_r₁, Gfn_fits_files_src,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1, 1, r₁_ind, 1)

		conjα₁_α₂ = complex(abs2(α_r₁[]))

		if r₁_ind != r₂_ind
			read_Gfn_file_at_index!(α_r₂, Gfn_fits_files_src,
				(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1, 1, r₂_ind, 1)

			conjα₁_α₂ = α_r₁[]*conj(α_r₂[])
		end

		α₁conjα₂ωℓ[ℓ, ω_ind] = conjα₁_α₂ * ω^2*Powspec(ω)
	end

	closeGfnfits(Gfn_fits_files_src)

	permutedims(parent(α₁conjα₂ωℓ))
end

# Without los, radial components, 3D points
function Cτ_rotating(comm, xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_radial = los_radial(); kwargs...)
	# Return C(Δϕ, ω) = RFFT(C(Δϕ,τ)) = RFFT(IRFFT(C0(Δϕ - Ωτ, ω))(τ))
	r_src, _, c_scale = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros, Nt, dt, dν = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)
	ℓ_range, ν_ind_range = iters

	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)
	τ_ind_range = get(kwargs, :τ_ind_range, 1:Nt)

	lmax = maximum(ℓ_range)

	Nτ = length(τ_ind_range)

	# The first step is loading in the αωℓ
	# Fourier transform the αωℓ and transpose to obtain αℓt
	Cω_in_range = pmapsum(comm, Cτ_rotating_partial, iters,
			xobs1, xobs2, los, p_Gsrc, r_src, c_scale)

	αωℓ = pad_zeros_ν(Cω_in_range, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros)

	αtℓ = fft_ω_to_t(αωℓ, dν)
	αℓt = permutedims(αtℓ)[:,τ_ind_range]

	np = min(nworkers(), Nτ)

	αℓt = distribute(αℓt, procs = workers()[1:np], dist=[1, np])

	τ_ind_tracker = RemoteChannel(()->Channel{Bool}(length(τ_ind_range)))
	progress_bar_τ = Progress(length(τ_ind_range), 1,"Cτ : ")

	Cτ_arr = @sync begin
		t = @async begin
			Cτ_arr =
			try
				DArray((Nτ,), workers()[1:np],[np]) do inds
					τ_ind_range_proc = τ_ind_range[first(inds)]
					αℓt_local = OffsetArray(αℓt[:lp], ℓ_range,τ_ind_range_proc)
					Cτ_arr = zeros(τ_ind_range_proc)
					for τ_ind in τ_ind_range_proc
						# τ goes from -T/2 to T/2-dt but is fftshifted
						τ = (τ_ind<=div(Nt, 2) ? (τ_ind-1) : (τ_ind-1- Nt)) * dt
						xobs2′ = Point3D(radius(xobs2), xobs2.θ, xobs2.ϕ-Ω_rot*τ)
						Pl_cosχ = collectPl(cosχ(xobs1, xobs2′), lmax = lmax)
						for ℓ in ℓ_range
							Cτ_arr[τ_ind] += (2ℓ+1)/4π * αℓt_local[ℓ,τ_ind] * Pl_cosχ[ℓ]
						end
						put!(τ_ind_tracker, true)
					end
					finalize(τ_ind_tracker)
					parent(Cτ_arr)
				end
			finally
				put!(τ_ind_tracker, false)
			end
		end
		while take!(τ_ind_tracker)
			next!(progress_bar_τ)
		end
		fetch(t)
	end

	OffsetArray(Array(Cτ_arr), τ_ind_range)
end

# # 2D points
# @two_points_on_the_surface Cτ_rotating

# Without los, radial components, multiple 3D points
function Cτ_rotating(comm, xobs1, xobs2_arr::Vector{<:SphericalPoint}, los::los_direction = los_radial();
	Ω_rot = 20e2/Rsun, τ_ind_range = nothing, kwargs...)

	# Return C(Δϕ, ω) = RFFT(C(Δϕ,τ)) = RFFT(IRFFT(C0(Δϕ - Ωτ, ω))(τ))
	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack Nt, dt = p_Gsrc

	if isnothing(τ_ind_range)
		τ_ind_range = [1:div(Nt, 2) for _ in xobs2_arr]
	end

	τ_ind_max_span = minimum(minimum.(τ_ind_range)):maximum(maximum.(τ_ind_range))

	Cτ_arr = [zeros(τ_inds) for τ_inds in τ_ind_range]

	@showprogress 1 "Cτ τ ind : " for τ_ind in τ_ind_max_span
		τ = (τ_ind-1) * dt
		xobs2′_arr = [shiftϕ(xobs2, -Ω_rot*τ)  for (rank,xobs2) in enumerate(xobs2_arr) if τ_ind in τ_ind_range[rank]]
		xobs2′inds_arr = [rank for rank in eachindex(xobs2_arr) if τ_ind in τ_ind_range[rank]]
		Ct_xobs2_arr = Ct(comm, xobs1, xobs2′_arr, los; kwargs...)[τ_ind, :]
		for (Ct_x2, x2ind) in zip(Ct_xobs2_arr, xobs2′inds_arr)
			Cτ_arr[x2ind][τ_ind] = Ct_x2
		end
	end

	return Cτ_arr
end

#######################################################################################################################################
# Allocate arrays
#######################################################################################################################################

function allocatearrays(::SoundSpeed, los::los_direction, obs_at_same_height)
	Gsrc = zeros(ComplexF64, 0:1, srcindG(los)..., nr)
	drGsrc = zeros(ComplexF64, 0:0, srcindG(los)..., nr)
	Gobs1 = zeros(ComplexF64, 0:1, srcindG(los)..., nr)
	drGobs1 = zeros(ComplexF64, 0:0, srcindG(los)..., nr)
	Gobs2, drGobs2 = Gobs1, drGobs1

	Gobs1_cache = Dict{Int, typeof(Gobs1)}()
	drGobs1_cache = Dict{Int, typeof(drGobs1)}()
	Gobs2_cache = obs_at_same_height ? Gobs1_cache : Dict{Int, typeof(Gobs2)}()
	drGobs2_cache = obs_at_same_height ? drGobs1_cache : Dict{Int, typeof(drGobs2)}()

	divGsrc = zeros(ComplexF64, srcindG(los)..., nr)
	divGobs = zeros(ComplexF64, srcindG(los)..., nr)

	# f_αjₒjₛ(r, rᵢ, rₛ) = -2ρc ∇⋅Gjₒ(r, rᵢ)_α ∇⋅Gjₛ(r, rₛ)_0
	_f = zeros(srcindG(los)..., nr)
	fjₒjₛ_r₁_rsrc = StructArray{ComplexF64}((_f, zero(_f)))
	fjₒjₛ_r₂_rsrc = obs_at_same_height ? fjₒjₛ_r₁_rsrc : zero(fjₒjₛ_r₁_rsrc)

	# H_βαjₒjₛ(r; r₁, r₂, rₛ) = conj(f_αjₒjₛ(r, r₁, rₛ)) Gβ0jₛ(r₂, rₛ)
	_H = zeros(obsindG(los)..., srcindG(los)..., nr)
	Hjₒjₛω_r₁r₂ = StructArray{ComplexF64}((_H, zero(_H)))
	Hjₒjₛω_r₂r₁ = obs_at_same_height ? Hjₒjₛω_r₁r₂ : zero(Hjₒjₛω_r₁r₂)

	tworealconjhωHjₒjₛω_r₁r₂ = zeros(obsindG(los)..., srcindG(los)..., nr)
	tworealconjhωconjHjₒjₛω_r₂r₁ = zero(tworealconjhωHjₒjₛω_r₁r₂)

	(; Gsrc, drGsrc, Gobs1, drGobs1, Gobs2, drGobs2,
		Gobs1_cache, drGobs1_cache, Gobs2_cache, drGobs2_cache,
		divGsrc, divGobs,
		fjₒjₛ_r₁_rsrc, fjₒjₛ_r₂_rsrc,
		Hjₒjₛω_r₁r₂, Hjₒjₛω_r₂r₁,
		tworealconjhωHjₒjₛω_r₁r₂, tworealconjhωconjHjₒjₛω_r₂r₁)
end

function allocatearrays(::Flow, los::los_direction, obs_at_same_height, _::Val{N} = Val(1)) where {N}
	Gsrc = zeros(ComplexF64, 0:1, srcindG(los)..., nr)
	drGsrc = zeros(ComplexF64, 0:1, srcindG(los)..., nr)
	Gobs1 = zeros(ComplexF64, 0:1, srcindG(los)..., nr)
	Gobs2 = Gobs1

	Gobs1_cache = Dict{Int, typeof(Gobs1)}()
	Gobs2_cache = Dict{Int, typeof(Gobs2)}()

	_G1 = zeros(0:1, srcindG(los)..., nr);
	_G2 = StructArray{ComplexF64}((_G1, zero(_G1)));
	Gparts_r₁ = OffsetVector([_G2, zero(_G2)], 0:1)
	Gparts_r₂ = obs_at_same_height ? Gparts_r₁ : zero.(Gparts_r₁)

	_G3 = zeros(0:1, srcindG(los)..., nr, 0:1);
	Gparts_r₁_2 = StructArray{ComplexF64}((_G3, zero(_G3)))
	Gparts_r₂_2 = obs_at_same_height ? Gparts_r₁_2 : zero(Gparts_r₁_2)

	# This is Gγℓjₒjₛω_α₁0(r, r₁, rₛ) as computed in the paper
	# Stored as Gγℓjₒjₛ_r₁[r,γ,αᵢ] = Gparts_r₁[γ][r, 0, αᵢ] + ζ(jₛ, jₒ, ℓ) Gparts_r₁[γ][:, 1, αᵢ]
	_G = zeros(obsindG(los)..., nr, 0:1);
	Gγℓjₒjₛ_r₁ = StructArray{ComplexF64}((_G, zero(_G)));
	Gγℓjₒjₛ_r₂ = obs_at_same_height ? Gγℓjₒjₛ_r₁ : zero(Gγℓjₒjₛ_r₁);

	# This is given by Hγℓjₒjₛω_α₁α₂(r, r₁, r₂) = conj(Gγℓjₒjₛω_α₁0(r, r₁, rₛ)) * Gjₛω_α₂0(r₂, rₛ)
	# Stored as Hγℓjₒjₛ_r₁r₂[r,γ,α₁,α₂] = conj(Gγℓjₒjₛ_r₁[r,γ,α₁]) * Gα₂r_r₂[α₂]
	# Hγℓjₒjₛ_r₁r₂ = zeros(ComplexF64, nr, 0:1, srcindG(los)..., obsindG(los)...)
	_H = zeros(obsindG(los)..., obsindG(los)..., nr, 0:1);
	Hγℓjₒjₛ_r₁r₂ = StructArray{ComplexF64}((_H, zero(_H)));
	Hγℓjₒjₛ_r₂r₁ = obs_at_same_height ? Hγℓjₒjₛ_r₁r₂ : zero(Hγℓjₒjₛ_r₁r₂);

	if N === 1
		twoimagconjhωHγℓjₒjₛ_r₁r₂ = zeros(obsindG(los)..., obsindG(los)..., nr, 0:1)
		twoimagconjhωconjHγℓjₒjₛ_r₂r₁ = zeros(obsindG(los)..., obsindG(los)..., nr, 0:1)
	else
		twoimagconjhωHγℓjₒjₛ_r₁r₂ = zeros(obsindG(los)..., obsindG(los)..., N, nr, 0:1)
		twoimagconjhωconjHγℓjₒjₛ_r₂r₁ = zeros(obsindG(los)..., obsindG(los)..., N, nr, 0:1)
	end

	# This is Gγℓjₒjₛ_r₁ for γ=1, used for the validation test
	G¹₁jj_r₁ = zeros(ComplexF64, srcindG(los)..., nr)
	G¹₁jj_r₂ = obs_at_same_height ? G¹₁jj_r₁ : zero(G¹₁jj_r₁)

	# Hjₒjₛαβ(r; r₁, r₂, rₛ) = conj(f_αjₒjₛ(r, r₁, rₛ)) Gβ0jₛ(r₂, rₛ)
	# We only use this for the validation case of ℑu⁺, so jₒ = jₛ = j
	H¹₁jj_r₁r₂ = zeros(ComplexF64, srcindG(los)..., obsindG(los)..., nr)
	H¹₁jj_r₂r₁ = obs_at_same_height ? H¹₁jj_r₁r₂ : zero(H¹₁jj_r₁r₂)

	(; Gsrc, Gobs1, Gobs2, drGsrc,
		Gobs1_cache, Gobs2_cache,
		Gparts_r₁, Gparts_r₂,
		Gparts_r₁_2, Gparts_r₂_2,
		Gγℓjₒjₛ_r₁, Gγℓjₒjₛ_r₂,
		Hγℓjₒjₛ_r₁r₂, Hγℓjₒjₛ_r₂r₁,
		twoimagconjhωconjHγℓjₒjₛ_r₂r₁,
		twoimagconjhωHγℓjₒjₛ_r₁r₂,
		G¹₁jj_r₁, G¹₁jj_r₂,
		H¹₁jj_r₁r₂, H¹₁jj_r₂r₁)
end

#######################################################################################################################################
# First Born approximation for flows
#######################################################################################################################################

function δCωℓ_FB!(δC_r, ω_ind, ω, ℓ, Y12ℓ, Gsrc, G¹₁jj_r₂, G¹₁jj_r₁, r₁_ind, r₂_ind, ::los_radial)
	G_r₁_rsrc = Gsrc[0, r₁_ind]
	G_r₂_rsrc = Gsrc[0, r₂_ind]

	pre = ω^3 * Powspec(ω) * Y12ℓ

	for r_ind in UnitRange(axes(δC_r, 1))
		δC_r[r_ind, ω_ind] += pre *
			(conj(G_r₁_rsrc) * G¹₁jj_r₂[r_ind] + G_r₂_rsrc * conj(G¹₁jj_r₁[r_ind]))
	end
	return δC_r
end

function δCωℓ_FB!(δC_r, ω_ind, ω, ℓ, Y12ℓ, Gsrc, G¹₁jj_r₂, G¹₁jj_r₁, r₁_ind, r₂_ind, ::los_earth)
	pre = ω^3 * Powspec(ω)* Ω(ℓ, 0) * √(1/6)
	for r_ind in UnitRange(axes(δC_r, 1))
		s = zero(eltype(δC_r))
		for β in 0:1, γ in 0:1
			Gγ_r₁_rsrc = Gsrc[γ, 0, r₁_ind]
			Gβ_r₂_rsrc = Gsrc[β, 0, r₂_ind]

			gγ_r₁_rsrc = G¹₁jj_r₁[γ, r_ind]
			gβ_r₂_rsrc = G¹₁jj_r₂[β, r_ind]

			s += (conj(Gγ_r₁_rsrc) * gβ_r₂_rsrc + Gβ_r₂_rsrc * conj(gγ_r₁_rsrc)) * Y12ℓ[γ, β]
		end

		δC_r[r_ind, ω_ind] += pre * s
	end
	return δC_r
end

_integrate_δCrω(δC_r, ::los_radial) = -2im*integrate(r, (@. r^2 * ρ * δC_r))
_integrate_δCrω(δC_r, ::los_earth) = -4integrate(r, (@. r^2 * ρ * δC_r))

function δCω_uniform_rotation_firstborn_integrated_over_angle_partial(
	localtimer,
	ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, r_obs = nothing, Ω_rot = 20e2/Rsun)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src)
	p_Gobs1 = read_all_parameters(p_Gobs1, r_src = radius(xobs1))
	p_Gobs2 = read_all_parameters(p_Gobs2, r_src = radius(xobs2))

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	Gfn_path_obs1, NGfn_files_obs1 =  p_Gobs1.path, p_Gobs1.num_procs
	Gfn_path_obs2, NGfn_files_obs2 = p_Gobs2.path, p_Gobs2.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2 =
			Gfn_fits_files.((Gfn_path_src, Gfn_path_obs1, Gfn_path_obs2),
			((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
			(NGfn_files_src, NGfn_files_obs1, NGfn_files_obs2))

	δC_r = zeros(ComplexF64, nr, ν_ind_range)

	arrs = allocatearrays(Flow(), los, r₁_ind == r₂_ind)
	@unpack Gsrc, Gobs1, Gobs2, G¹₁jj_r₁, G¹₁jj_r₂ = arrs;

	Y12 = los_projected_biposh_spheroidal(computeY₁₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		# Green function about source location
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		# Green function about receiver location
		read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)

		radial_fn_uniform_rotation_firstborn!(G¹₁jj_r₁, Gsrc, Gobs1, ℓ, los)

		if r₁_ind != r₂_ind
			# Green function about receiver location
	    	read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)

			radial_fn_uniform_rotation_firstborn!(G¹₁jj_r₂, Gsrc, Gobs2, ℓ, los)
		end

		δCωℓ_FB!(δC_r, ω_ind, ω, ℓ, Y12[ℓ], Gsrc, G¹₁jj_r₂, G¹₁jj_r₁, r₁_ind, r₂_ind, los)
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))

	δCω = Ω_rot * _integrate_δCrω(δC_r, los)

	parent(δCω)
end

# 3D points
function δCω_uniform_rotation_firstborn_integrated_over_angle(comm,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction = los_radial(); kwargs...)

	r_src, r_obs = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc = read_all_parameters(r_src = r_src)
	p_Gobs1 = read_all_parameters(r_src = radius(xobs1))
	p_Gobs2 = read_all_parameters(r_src = radius(xobs2))
	@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)

	Cω_in_range = pmapsum(comm,
		δCω_uniform_rotation_firstborn_integrated_over_angle_partial,
		iters, xobs1, xobs2, los, p_Gsrc, p_Gobs1, p_Gobs2, r_src, r_obs, Ω_rot)

	Cω_in_range === nothing && return nothing
	ν_ind_range = last(iters)
	δCω_flows_FB = pad_zeros_ν(Cω_in_range, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros)

	save_to_fits_and_return("δCω_flows_FB.fits", δCω_flows_FB)
end

# # 2D points
# @two_points_on_the_surface δCω_uniform_rotation_firstborn_integrated_over_angle

#######################################################################################################################################
# δC(ω) = C(ω) - C0(ω) for flows
########################################################################################################################################

# Linearized, without los, radial components, 3D points
function δCω_uniform_rotation_rotatedwaves_linearapprox(comm, xobs1, xobs2,
	los::los_radial = los_radial(); kwargs...)

	# We compute δC(xobs1, xobs2) = -iΩ ∑ℓ (2ℓ+1)/4π ∂ω (ω^2 P(ω) αℓω*(r₂)αℓω(r₁))∂Δϕ Pl(cos(Δϕ))

	r_src = get(kwargs, :r_src, r_src_default)
	Gfn_path_src = Gfn_path_from_source_radius(r_src, c_scale = get(kwargs, :c_scale, 1))
	@load joinpath(Gfn_path_src,"parameters.jld2") dω

	∂ϕC = ∂ϕ₂Cω(comm, xobs1, xobs2, los; kwargs...)
	∂ϕC === nothing && return nothing

	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)
	∂ω∂ϕC = D(size(∂ϕC, 1))*∂ϕC ./ dω
	δCω_flows_rotated_linear = @. -im*Ω_rot*∂ω∂ϕC
	save_to_fits_and_return("δCω_flows_rotated_linear.fits", δCω_flows_rotated_linear)
end

# With or without los, 3D points
function δCω_uniform_rotation_rotatedwaves(comm, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	dt = get(kwargs, :dt) do
		read_parameters("dt"; kwargs...)[1]
	end
	δCt = δCt_uniform_rotation_rotatedwaves(comm, xobs1, xobs2, los; kwargs...)
	δCt === nothing && return nothing
	δCω_flows_rotated = fft_t_to_ω(δCt, dt)
	save_to_fits_and_return("δCω_flows_rotated.fits", δCω_flows_rotated)
end

# With or without los, using precomputed time-domain cross-covariance
function δCω_uniform_rotation_rotatedwaves(δCt::AbstractArray; kwargs...)
	dt = get(kwargs, :dt) do
		read_parameters("dt"; kwargs...)[1]
	end
	δCω_flows_rotated = fft_t_to_ω(δCt, dt)
	return save_to_fits_and_return("δCω_flows_rotated.fits", δCω_flows_rotated)
end

#######################################################################################################################################
# δC(t) = C(t) - C0(t) for flows
########################################################################################################################################

# With or without los
function δCt_uniform_rotation_rotatedwaves(comm, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	dν = get(kwargs, :dν) do
		read_parameters("dν"; kwargs...)[1]
	end

	C′_t = Cτ_rotating(comm, xobs1, xobs2, los; kwargs...)
	C0_t = get(kwargs, :Ct) do
		C0_ω = get(kwargs, :Cω) do
			Cω(comm, xobs1, xobs2, los; kwargs...)
		end
		C0_ω === nothing ? nothing : Ct(C0_ω; dν=dν, kwargs...)
	end
	C′_t === nothing && return nothing
	δCt_flows_rotated = parent(C′_t) .- parent(C0_t)
	save_to_fits_and_return("δCt_flows_rotated.fits", δCt_flows_rotated)
end

# With or without los, multiple points
function δCt_uniform_rotation_rotatedwaves(comm, xobs1, xobs2_arr::Vector,
	los::los_direction = los_radial(); kwargs...)

	C′_t = Cτ_rotating(comm, xobs1, xobs2_arr, los; kwargs...)
	C0_t = get(kwargs, :Ct) do
			[(ind,τ_ind) -> begin
					C = Ct(comm, xobs1, xobs2_arr, los; kwargs...)
					C === nothing && return nothing
					τ_ind = axes(C′_t[ind],1)
					C[τ_ind, ind]
				end
				for ind in eachindex(xobs2_arr)]
	end
	C′_t === nothing && return nothing
	return parent(C′_t) .- C0_t
end

# Linearized, with or without los
function δCt_uniform_rotation_rotatedwaves_linearapprox(comm, xobs1, xobs2,
	los::los_direction = los_radial(); kwargs...)

	C = ∂ϕ₂Ct(comm, xobs1, xobs2, los; kwargs...)
	C === nothing && return nothing
	Nt = read_parameters("Nt"; kwargs...)
	τ_ind_range = get(kwargs, :τ_ind_range, 1:Nt)
	t = fftfreq(Nt, dt*Nt)
	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)
	@views -Ω_rot .* t[τ_ind_range] .* C[τ_ind_range]
end

# Linearized, with or without los, multiple 3D points
function δCt_uniform_rotation_rotatedwaves_linearapprox(comm, xobs1::Point3D,
	xobs2_arr::Vector{<:Point3D}, los::los_direction = los_radial(); kwargs...)

	_, Nt = read_parameters("Nt"; kwargs...)

	t = fftfreq(Nt, dt*Nt)
	τ_ind_range = get(kwargs, :τ_ind_range, [1:Nt for _ in xobs2_arr])

	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)

	C = ∂ϕ₂Ct(comm, xobs1, xobs2_arr, los; kwargs...)
	C === nothing && return nothing
	δCt = -Ω_rot .* t .* C
	map(zip(τ_ind_range, eachindex(xobs2_arr))) do (τ_ind, x2ind)
		δCt[τ_ind, x2ind]
	end
end

# Linearized, with or without los, multiple 2D points
function δCt_uniform_rotation_rotatedwaves_linearapprox(comm, nobs1::Point2D,
	nobs2_arr::Vector{<:Point2D}, los::los_direction = los_radial(); kwargs...)

	C = ∂ϕ₂Ct(nobs1, nobs2_arr, los; kwargs...)
	C === nothing && return nothing
	δCt_uniform_rotation_rotatedwaves_linearapprox(C; kwargs...)
end

# Linearized, with or without los, using precomputed time-domain cross-covariance
function δCt_uniform_rotation_rotatedwaves_linearapprox(∂ϕ₂Ct_arr; kwargs...)
	Nt, dt = read_parameters("Nt","dt"; kwargs...)

	τ_ind_range = get(kwargs, :τ_ind_range, 1:Nt)
	t = fftfreq(Nt, dt*Nt)

	Ω_rot = get(kwargs, :Ω_rot, 20e2/Rsun)

	@views @. -Ω_rot * t[τ_ind_range] * ∂ϕ₂Ct_arr[τ_ind_range]
end

#######################################################################################################################################
# First Born approximation for sound speed perturbation
#######################################################################################################################################

function δC_FB_soundspeed!(δC_rω, ω_ind, ω, Y12ℓ, H¹₁jj_r₂r₁, H¹₁jj_r₁r₂, ::los_radial)
	pre = ω^2*Powspec(ω)* Y12ℓ
	@. δC_rω[:, ω_ind] += pre * (conj(H¹₁jj_r₂r₁) + H¹₁jj_r₁r₂)
	return δC_rω
end

function δC_FB_soundspeed!(δC_rω, ω_ind, ω, Y12ℓ, Hjₒjₛω_r₂r₁, Hjₒjₛω_r₁r₂, ::los_earth)
	pre = ω^2 * Powspec(ω)
	for r_ind in eachindex(r)
		YH = zero(eltype(Y12ℓ))
		for β in 0:1, α in 0:1
			YH += Y12ℓ[α, β] * ( Hjₒjₛω_r₁r₂[α, β, r_ind] + conj(Hjₒjₛω_r₂r₁[β, α, r_ind]) )
		end
		δC_rω[r_ind, ω_ind] += pre * YH
	end
	return δC_rω
end

function δCω_isotropicδc_firstborn_integrated_over_angle_partial(localtimer,
	ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_direction,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs1::Union{Nothing, ParamsGfn} = nothing,
	p_Gobs2::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, r_obs = nothing, c_scale = 1+1e-5)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = 1)
	p_Gobs1 = read_all_parameters(p_Gobs1, r_src = radius(xobs1), c_scale = 1)
	p_Gobs2 = read_all_parameters(p_Gobs2, r_src = radius(xobs2), c_scale = 1)

	Gfn_path_src, NGfn_files_src = p_Gsrc.path, p_Gsrc.num_procs
	Gfn_path_obs1, NGfn_files_obs1 =  p_Gobs1.path, p_Gobs1.num_procs
	Gfn_path_obs2, NGfn_files_obs2 = p_Gobs2.path, p_Gobs2.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	ϵ = c_scale-1

	Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2 =
		Gfn_fits_files.((Gfn_path_src, Gfn_path_obs1, Gfn_path_obs2),
		((ℓ_arr, 1:Nν_Gfn),), (ℓ_ωind_iter_on_proc,),
		(NGfn_files_src, NGfn_files_obs1, NGfn_files_obs2))

	δC_rω = zeros(ComplexF64, nr, ν_ind_range)

	arrs = allocatearrays(SoundSpeed(), los, r₁_ind == r₂_ind)
	@unpack Gsrc, drGsrc, divGsrc, Gobs1, drGobs1, Gobs2, divGobs = arrs
	@unpack Hjₒjₛω_r₁r₂, Hjₒjₛω_r₂r₁ = arrs
	@unpack fjₒjₛ_r₁_rsrc, fjₒjₛ_r₂_rsrc = arrs

	Y12 = los_projected_biposh_spheroidal(computeY₀₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		# Green function about the source location
		read_Gfn_file_at_index!(Gsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:2, srcindFITS(los), :, 1)

		Gγr_r₁_rsrc = αrcomp(Gsrc, r₁_ind, los)
		Gγr_r₂_rsrc = αrcomp(Gsrc, r₂_ind, los)

		# Derivative of Green function about the source location
		read_Gfn_file_at_index!(drGsrc, Gfn_fits_files_src,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_src, 1:1, srcindFITS(los), :, 2)

		# Green function about the receiver location
		read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs1, 1:2, srcindFITS(los), :, 1)

		# Derivative of Green function about the receiver location
		read_Gfn_file_at_index!(drGobs1, Gfn_fits_files_obs1,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs1, 1:1, srcindFITS(los), :, 2)

		radial_fn_isotropic_δc_firstborn!(fjₒjₛ_r₁_rsrc,
			Gsrc, drGsrc, divGsrc, Gobs1, drGobs1, divGobs, ℓ)

		Hjₒjₛω!(Hjₒjₛω_r₁r₂, fjₒjₛ_r₁_rsrc, Gγr_r₂_rsrc)

		if r₁_ind != r₂_ind
			# Green function about the receiver location
    		read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
    			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs2, 1:2, srcindFITS(los), :, 1)

    		# Derivative of Green function the about receiver location
    		read_Gfn_file_at_index!(drGobs2, Gfn_fits_files_obs2,
    			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_obs2, 1:1, srcindFITS(los), :, 2)

			radial_fn_isotropic_δc_firstborn!(fjₒjₛ_r₂_rsrc,
				Gsrc, drGsrc, divGsrc, Gobs2, drGobs2, divGobs, ℓ)

			Hjₒjₛω!(Hjₒjₛω_r₂r₁, fjₒjₛ_r₂_rsrc, Gγr_r₁_rsrc)
		end

		δC_FB_soundspeed!(δC_rω, ω_ind, ω, Y12[ℓ], Hjₒjₛω_r₂r₁, Hjₒjₛω_r₁r₂, los)
	end

	foreach(closeGfnfits, (Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))

	C = parent(δC_rω)
	integrate(r, (@. r^2 * (ϵ*c) * C))
end

function _δCω_isotropicδc_firstborn(comm, xobs1, xobs2, args...; kwargs...)
	r_src, r_obs, c_scale = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs..., c_scale = 1)
	@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	δCω_in_range = pmapsum(comm,
		δCω_isotropicδc_firstborn_integrated_over_angle_partial,
		iters, xobs1, xobs2, args...,
		p_Gsrc, p_Gobs1, p_Gobs2, r_src, r_obs, c_scale)

	δCω_in_range === nothing && return nothing

	ν_ind_range = last(iters)
	pad_zeros_ν(δCω_in_range, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros)
end

# With our without los, 3D points
function δCω_isotropicδc_firstborn_integrated_over_angle(comm, xobs1::SphericalPoint, xobs2::SphericalPoint,
	los::los_direction = los_radial(); kwargs...)

	δCω_isotropicδc_FB = _δCω_isotropicδc_firstborn(comm, xobs1, xobs2, los; kwargs...)
	δCω_isotropicδc_FB === nothing && return nothing

	tag = los isa los_radial ? "" : "_los"
	filename = "δCω_isotropicδc_FB$tag.fits"
	save_to_fits_and_return(filename, δCω_isotropicδc_FB)
end

# # 2D points
# @two_points_on_the_surface δCω_isotropicδc_firstborn_integrated_over_angle

function δCt_isotropicδc_firstborn_integrated_over_angle(comm, xobs1, xobs2,
	los::los_direction = los_radial(); kwargs...)

	dν = read_parameters("dν"; kwargs...)[1]

	δCω = δCω_isotropicδc_firstborn_integrated_over_angle(comm, xobs1, xobs2, los; kwargs...)
	δCω === nothing && return nothing
	δCt_isotropicδc_FB = fft_ω_to_t(δCω, dν)

	tag = los isa los_radial ? "" : "_los"
	filename = "δCt_isotropicδc_FB$tag.fits"
	save_to_fits_and_return(filename, δCt_isotropicδc_FB)
end

#######################################################################################################################################
# δC(ω) = C(ω) - C0(ω) for sound speed perturbations
#######################################################################################################################################

function δCω_isotropicδc_C_minus_C0_partial(localtimer,
	ℓ_ωind_iter_on_proc::ProductSplit,
	xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_earth,
	p_Gsrc::Union{Nothing, ParamsGfn} = nothing,
	p_G′src::Union{Nothing, ParamsGfn} = nothing,
	r_src = r_src_default, r_obs = nothing, c_scale = 1+1e-5)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src = r_src, c_scale = 1)
	p_G′src = read_all_parameters(p_G′src, r_src = r_src, c_scale = c_scale)

	Gfn_path_src_c0, NGfn_files_c0 = p_Gsrc.path, p_Gsrc.num_procs
	Gfn_path_src_c′, NGfn_files_c′ =  p_G′src.path, p_G′src.num_procs
	@unpack ℓ_arr, ω_arr, Nν_Gfn = p_Gsrc

	_, ν_ind_range = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	Gfn_fits_files_src_c0, Gfn_fits_files_src_c′ =
	Gfn_fits_files.((Gfn_path_src_c0, Gfn_path_src_c′), ((ℓ_arr, 1:Nν_Gfn),),
		(ℓ_ωind_iter_on_proc,), (NGfn_files_c0, NGfn_files_c′))

	Cω = zeros(ComplexF64, ν_ind_range)

	r₁_ind, r₂_ind = radial_grid_index.((xobs1, xobs2))

	G1_c′ = zeros(ComplexF64, 0:1)
	G1_c0 = zeros(ComplexF64, 0:1)
	if r₁_ind != r₂_ind
		G2_c′ = zeros(ComplexF64, 0:1)
		G2_c0 = zeros(ComplexF64, 0:1)
	else
		G2_c′ = G1_c′
		G2_c0 = G1_c0
	end

	Y12 = los_projected_biposh_spheroidal(computeY₀₀, xobs1, xobs2, los, ℓ_ωind_iter_on_proc)

	for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc

		ω = ω_arr[ω_ind]

		# Green function about source location at observation point 1
		read_Gfn_file_at_index!(G1_c0, Gfn_fits_files_src_c0,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_c0, 1:2, 1, r₁_ind, 1)

		read_Gfn_file_at_index!(G1_c′, Gfn_fits_files_src_c′,
		(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_c′, 1:2, 1, r₁_ind, 1)

		# Green function about source location at observation point 2
		if r₁_ind != r₂_ind
			read_Gfn_file_at_index!(G2_c0, Gfn_fits_files_src_c0,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_c0, 1:2, 1, r₂_ind, 1)

			read_Gfn_file_at_index!(G2_c′, Gfn_fits_files_src_c′,
			(ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), NGfn_files_c′, 1:2, 1, r₂_ind, 1)
		end

		Y12ℓ = Y12[ℓ]

		s = zero(eltype(Cω))
		for β in 0:1, α in 0:1
			GG_c0 = conj(G1_c0[α]) * G2_c0[β]
			GG_c′ = conj(G1_c′[α]) * G2_c′[β]
			δGG = GG_c′ - GG_c0
			s += δGG * Y12ℓ[α, β]
		end
		Cω[ω_ind] += ω^2 * Powspec(ω) * s
	end

	foreach(closeGfnfits, (Gfn_fits_files_src_c′, Gfn_fits_files_src_c0))

	parent(Cω)
end

function _δCω_isotropicδc_C_minus_C0(comm, xobs1, xobs2,
	los::los_earth, args...; kwargs...)
	r_src, r_obs, c_scale = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc = read_all_parameters(r_src = r_src, c_scale = 1)
	p_G′src = read_all_parameters(r_src = r_src, c_scale = c_scale)
	@unpack Nν_Gfn, ν_arr, ℓ_arr, ν_start_zeros, ν_end_zeros = p_Gsrc
	iters = ℓ_and_ν_range(ℓ_arr, ν_arr; kwargs...)

	Cω_in_range = pmapsum(comm,
		δCω_isotropicδc_C_minus_C0_partial,
		iters, xobs1, xobs2, los, args..., p_Gsrc, p_G′src,
		r_src, r_obs, c_scale)

	Cω_in_range === nothing && return nothing

	ν_ind_range = last(iters)
	pad_zeros_ν(Cω_in_range, ν_ind_range, Nν_Gfn, ν_start_zeros, ν_end_zeros)
end

# Without los, radial components, 3D points
function δCω_isotropicδc_C_minus_C0(comm, xobs1::SphericalPoint, xobs2::SphericalPoint,
	los::los_radial = los_radial(); kwargs...)

	c_scale = get(kwargs, :c_scale, 1+1e-5)

	C′ = Cω(comm, xobs1, xobs2, los; kwargs..., c_scale = c_scale)
	C0 = Cω(comm, xobs1, xobs2, los; kwargs..., c_scale = 1)

	C′ === nothing && return nothing

	δCω_isotropicδc_CmC0 = @. C′- C0
	save_to_fits_and_return("δCω_isotropicδc_CmC0.fits", δCω_isotropicδc_CmC0)
end

# With los, 3D points
function δCω_isotropicδc_C_minus_C0(comm, xobs1::SphericalPoint, xobs2::SphericalPoint, los::los_earth; kwargs...)
	δCω_isotropicδc_CmC0_los = _δCω_isotropicδc_C_minus_C0(comm, xobs1, xobs2, los; kwargs...)
	δCω_isotropicδc_CmC0_los === nothing && return nothing
	save_to_fits_and_return("δCω_isotropicδc_CmC0_los.fits", δCω_isotropicδc_CmC0_los)
end

# # 2D points
# @two_points_on_the_surface δCω_isotropicδc_C_minus_C0

########################################################################################
# δC(t) = C(t) - C0(t) for sound speed perturbations
########################################################################################

# With or without los, 3D points
function δCt_isotropicδc_C_minus_C0(comm, xobs1, xobs2,
	los::los_direction = los_radial(); kwargs...)

	dν = read_parameters("dν"; kwargs...)[1]

	δCω = δCω_isotropicδc_C_minus_C0(comm, xobs1, xobs2, los; kwargs...)
	δCω === nothing && return nothing
	δCt_isotropicδc_CmC0 = fft_ω_to_t(δCω, dν)

	tag = los isa los_radial ? "" : "_los"
	filename = "δCt_isotropicδc_CmC0$tag.fits"
	save_to_fits_and_return(filename, δCt_isotropicδc_CmC0)
end

# With or without los, precomputed arrays
function δCt_isotropicδc_C_minus_C0(δCω::AbstractArray; kwargs...)
	dν = read_parameters("dν"; kwargs...)[1]
	δCt_isotropicδc_CmC0 = fft_ω_to_t(δCω, dν)
	save_to_fits_and_return("δCt_isotropicδc_CmC0.fits", δCt_isotropicδc_CmC0)
end

########################################################################################################################
# Window function
########################################################################################################################

function ridgefilterspectra(Cωℓ, ν_arr, nγ = 10)
	Nν_Gfn = length(ν_arr)
	map(axes(Cωℓ, 2)) do ℓ
		Cℓ = @view Cωℓ[:, ℓ];
		Cℓ = diff(sign.(diff(Cℓ)));
		ν_inds = findall(isequal(-2), Cℓ)
		ν_ℓ = ν_arr[ν_inds]
		γ_ℓ = γ_damping.(ν_ℓ * 2pi) ./ 2pi
		ν_low = ν_ℓ .- nγ * γ_ℓ
		ν_low_inds = clamp.(searchsortedfirst.((ν_arr,), ν_low), 1, Nν_Gfn)
		ν_high = ν_ℓ .+ nγ * γ_ℓ
		ν_high_inds = clamp.(searchsortedfirst.((ν_arr,), ν_high), 1, Nν_Gfn)
		inds = UnitRange{Int}[]
		iters = zip(ν_low_inds, ν_high_inds)
		for (ind, (ν_low_ind, ν_high_ind)) in enumerate(iters)
			if ind == length(iters)
				push!(inds, ν_low_ind:ν_high_ind)
				continue
			end
			next_low_ind = minimum(@view ν_low_inds[ind+1:end])
			if ν_high_ind < next_low_ind
				# isolated peak
				push!(inds, ν_low_ind:ν_high_ind)
			else
				push!(inds, ν_low_ind:ν_high_inds[end])
				break
			end
		end
		inds
	end
end

function bounce_filter(Δϕ, n)
	nparams = 5
	coeffs = Dict()
	for i in [1, 2, 4]
		coeffs[i] = Dict("high"=>zeros(nparams),"low"=>zeros(nparams))
	end

	coeffs[1]["high"] = [2509.132334896018, 12792.508296270391,-13946.527195127102, 8153.75242742649,-1825.7567469552703]
	coeffs[1]["low"] = [40.821191938380714, 11410.21390421857,-11116.305124138207, 5254.244817703224,-895.0009393800744]

	coeffs[2]["high"] = [4083.6946001848364, 14924.442447995087,-13880.238239469609, 7562.499279468063,-1622.5318939228978]
	coeffs[2]["low"] = [2609.4406668522433, 10536.81683213881,-7023.811081076518, 2586.7238222832298,-348.9245124332354]

	coeffs[4]["high"] = [6523.103468645263, 16081.024611219753,-7085.7174198723405, 973.4990690666436, 236.95568587146957]
	coeffs[4]["low"] = [5150.314633252216, 15040.045600508669,-8813.047362534506, 3878.5398150601663,-870.3633232120256]

	τ_low,τ_high = 0., 0.
	for (i, c) in enumerate(coeffs[n]["high"])
		τ_high += c*Δϕ^(i-1)
	end

	for (i, c) in enumerate(coeffs[n]["low"])
		τ_low += c*Δϕ^(i-1)
	end

	return τ_low,τ_high
end

function gaussian_fit(x, y)
	# use the fact that if y = Gaussian(x), log(y) = quadratic(x)
	# quadratic(x) = ax² + bx + c
	# Gaussian(x) = A*exp(-(x-x0)²/2σ²)
	# the parameters are σ=√(-1/2a), x0 = -b/2a, A = exp(c-b^2/4a)
	c, b, a = polyfit(x, log.(y), 2).a
	A = exp(c-b^2/4a)
	x0 = -b/2a
	σ = √(-1/2a)
	return A, x0, σ
end

function time_window_bounce_filter(xobs1, xobs2, dt, bounce_no = 1)
	time_window_bounce_filter(acos(cosχ(xobs1, xobs2)), dt, bounce_no)
end

function time_window_bounce_filter(Δϕ::Real, dt, bounce_no = 1)
	τ_low,τ_high = bounce_filter(Δϕ, bounce_no)
	τ_low_ind = floor(Int64,τ_low/dt);
	τ_high_ind = ceil(Int64,τ_high/dt)
	return τ_low_ind,τ_high_ind
end

function time_window_indices_by_fitting_bounce_peak(C_t::AbstractVector{<:Real},
	xobs1, xobs2; dt, bounce_no = 1, kwargs...)

	τ_low_ind,τ_high_ind = time_window_bounce_filter(xobs1, xobs2, dt, bounce_no)
	time_window_indices_by_fitting_bounce_peak(C_t,τ_low_ind,τ_high_ind; dt = dt, kwargs...)
end

hilbertenvelope(C_t, Nt = size(C_t, 1)) = abs.(hilbert(C_t[1:div(Nt, 2)]))

function time_window_indices_by_fitting_bounce_peak(C_t::AbstractVector{<:Real},
	τ_low_ind::Int64, τ_high_ind::Int64;
	dt, Nt = size(C_t, 1), kwargs...)

	env = hilbertenvelope(C_t)
	peak_center = argmax(env[τ_low_ind:τ_high_ind]) + τ_low_ind - 1
	points_around_max = env[peak_center-2:peak_center+2]
	_, t0,σt = gaussian_fit(peak_center-2:peak_center+2, points_around_max)

	max(1, floor(Int64, t0 - 2σt)):min(ceil(Int64, t0 + 2σt), Nt)
end

function time_window_indices_by_fitting_bounce_peak(C_t::AbstractMatrix{<:Real},
	xobs1, xobs2_arr::Vector, args...; kwargs...)

	t_inds_range = Vector{UnitRange}(undef, size(C_t, 2))
	for (x2ind, xobs2) in enumerate(xobs2_arr)
		t_inds_range[x2ind] = time_window_indices_by_fitting_bounce_peak(C_t[:, x2ind],
								xobs1, xobs2; kwargs...)
	end
	return t_inds_range
end

function time_window_indices_by_fitting_bounce_peak(C_t::AbstractMatrix{<:Real},
	(xobs1, xobs1′)::NTuple{2,SphericalPoint},
	(xobs2, xobs2′)::NTuple{2,SphericalPoint}, args...; kwargs...)

	@assert size(C_t, 2) == 2 "function is only defined for two sets of points related by a rotation"
	t_inds_range = Vector{UnitRange}(undef, size(C_t, 2))
	t_inds_range[1] = time_window_indices_by_fitting_bounce_peak(C_t[:, 1],
							xobs1, xobs2; kwargs...)
	t_inds_range[2] = time_window_indices_by_fitting_bounce_peak(C_t[:, 2],
							xobs1′, xobs2′; kwargs...)
	return t_inds_range
end

function time_window(a::Vector, τ_ind_range)
	b = zero(a)
	b[τ_ind_range] .= 1
	return b
end

function time_window(a::Matrix, τ_ind_range)
	b = zero(a)
	for (idx, τ_inds) in enumerate(τ_ind_range)
		b[τ_inds, idx] .= 1
	end
	return b
end

function ht(::TravelTimes, Cω_x1x2::AbstractArray, xobs1, xobs2; bounce_no = 1, kwargs...)
	Cω_x1x2 === nothing && return nothing
	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack ν_full, Nt, dt, dν = p_Gsrc

	ω_full = 2π*ν_full

	C_t = fft_ω_to_t(Cω_x1x2, dν)
	∂tCt_ω = @. Cω_x1x2 * im * ω_full
	∂tCt = fft_ω_to_t(∂tCt_ω, dν)

	τ_ind_range = get(kwargs, :τ_ind_range, nothing)
	if isnothing(τ_ind_range)
		τ_ind_range = time_window_indices_by_fitting_bounce_peak(C_t, xobs1, xobs2;
					dt = dt, Nt = Nt, bounce_no = bounce_no)
	end

	f_t = time_window(∂tCt, τ_ind_range)

	(@. -f_t * ∂tCt) ./ sum((@. f_t * ∂tCt^2 * dt), dims = 1)
end

function ht(::Amplitudes, Cω_x1x2::AbstractArray, xobs1, xobs2; bounce_no = 1, kwargs...)
	Cω_x1x2 === nothing && return nothing
	Nt, dt, dν = read_parameters("Nt","dt","dν"; kwargs...)

	C_t = fft_ω_to_t(Cω_x1x2, dν)

	τ_ind_range = get(kwargs, :τ_ind_range, nothing)
	if isnothing(τ_ind_range)
		τ_ind_range = time_window_indices_by_fitting_bounce_peak(C_t, xobs1, xobs2;
					dt = dt, Nt = Nt, bounce_no = bounce_no)
	end

	f_t = time_window(C_t, τ_ind_range)

	(@. f_t * C_t) ./ sum((@. f_t*C_t^2 * dt), dims = 1)
end

function ht(comm, m::SeismicMeasurement, xobs1, xobs2, los::los_direction = los_radial(); kwargs...)
	Cω_x1x2 = Cω(comm, xobs1, xobs2, los; kwargs...)
	Cω_x1x2 === nothing && return nothing
	ht(m, Cω_x1x2, xobs1, xobs2; kwargs...)
end

function hω(args...; kwargs...)
	h_t = ht(args...; kwargs...)
	h_t === nothing && return nothing
	dt, ν_start_zeros, ν_arr = read_parameters("dt", "ν_start_zeros", "ν_arr"; kwargs...)
	h_ω = fft_t_to_ω(h_t, dt)
	ν_ind_range = get(kwargs, :ν_ind_range, axes(ν_arr, 1)) .+ ν_start_zeros
	h_ω[ν_ind_range, ..]
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

sym=Dict("$(Char(i))"=>i for i=97:122)
A = WFST(sym)
W = typeofweight(A)
add_arc!(A,1,2,"a","a",1)
add_arc!(A,1,3,"a","a",2)
add_arc!(A,2,2,"b","b",3)
add_arc!(A,3,3,"b","b",3)
add_arc!(A,2,4,"c","c",5)
add_arc!(A,3,4,"d","d",6)
initial!(A,1)
final!(A,4)
println(A)
@test is_deterministic(A) == false
@test is_acceptor(A)
D = determinize_fsa(A)
#println(D)
@test is_deterministic(D)
@test size(D) == (3,4)
# TODO add more tests

ext: jl


text:
@testset "Iterations" begin

    @testset "NMI" begin
        # This is the example of MI(X, Y)/H(X, Y) from:
        #
        # Gao, H., Dou, Y., Yang, J., & Wang, J. (2011).
        # New methods to measure residues coevolution in proteins.
        # BMC bioinformatics, 12(1), 206.

        aln = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA)
        result = Float64[ 0     0     0     0     0     0
                          0     0     0     0     0     0
                          0     0     0     1     1     0.296
                          0     0     1     0     1     0.296
                          0     0     1     1     0     0.296
                          0     0     0.296 0.296 0.296 0 ]

        nmi = mapcolpairfreq!(normalized_mutual_information, aln,
                              Counts(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                              Val{false})
        nmi_mat = convert(Matrix{Float64}, getarray(nmi))
        @test isapprox(nmi_mat, result, rtol=1e-4)

        nmi_t = mapseqpairfreq!(normalized_mutual_information, transpose(aln),
                                Counts(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                                Val{false})
        @test nmi_mat == convert(Matrix{Float64},getarray(nmi_t))
    end

    @testset "Gaps" begin

        function _gaps(table::Union{Probabilities{Float64,1,GappedAlphabet},
                                    Counts{Float64,1,GappedAlphabet}})
            table[GAP]
        end

        table = ContingencyTable(Float64, Val{1}, GappedAlphabet())

        gaps = read(joinpath(pwd(), "data", "gaps.txt"), Raw)

        # THAYQAIHQV 0
        # THAYQAIHQ- 0.1
        # THAYQAIH-- 0.2
        # THAYQAI--- 0.3
        # THAYQA---- 0.4
        # THAYQ----- 0.5
        # THAY------ 0.6
        # THA------- 0.7
        # TH-------- 0.8
        # T--------- 0.9

        ngaps = Float64[0,1,2,3,4,5,6,7,8,9]

        colcount = mapcolfreq!(_gaps, gaps, Counts(table))
        @test all((vec(getarray(colcount)) .- ngaps) .== 0.0)

        colfract = mapcolfreq!(_gaps, gaps, Probabilities(table))
        @test all((vec(getarray(colfract)) .- ngaps./10.0) .== 0.0)

        seqcount = mapseqfreq!(_gaps, gaps, Counts(table))
        @test all((vec(getarray(seqcount)) .- ngaps) .== 0.0)

        seqfract = mapseqfreq!(_gaps, gaps, Probabilities(table))
        @test all((vec(getarray(seqfract)) .- ngaps./10.0) .== 0.0)
    end
end


text:
isless{N, ES, mode}(lhs::Sigmoid{N, ES, mode}, rhs::Sigmoid{N, ES, mode}) = @s(lhs) < @s(rhs)

Base.:<{N, ES, mode}(lhs::Sigmoid{N, ES, mode}, rhs::Sigmoid{N, ES, mode}) = (@u(rhs) == (@signbit)) | isless(lhs, rhs)
Base.:<={N, ES, mode}(lhs::Sigmoid{N, ES, mode}, rhs::Sigmoid{N, ES, mode}) = (@u(rhs) == (@signbit)) | (@s(lhs) <= @s(rhs))

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
# using SuperResolution
using Distributions
using ProgressMeter
include("fctSR.jl")
include("segmentSR.jl")

###4.. Setup KLClustering

##	Generate all possible clustering models for a _specific segment_
#	For final use, create such ModelSR for every valid SEG.idx
mutable struct ModelSR
		#Segment-restricted features
	idx::Tuple{Int64,Int64}
	modus::Int64
	dim::Int64
	loc::Array{Array{Float64,2},1}	#localisations of cluster populations
	frame::Array{Int64,1}
	model::Array{MixtureModel,1}	#MixtureModel
	member::Array{Array{Int64,1},1}	#cluster ID for all localisations
	clus_ID::Array{Array{Int64,1},1}	#cluster ID for remaining populations (needed for dim2WKLDfull; might replace with unique() and findall() for indices)
	mol_ID::Array{Int64,1}		# ground truth member ID

	cov::Array{Array{Array{Float64,2},1},1}	#covariance matrices
	comp_mix::Array{Array{Float64,1},1}	#convex coefficient of all MixtureModel components
	wkld::Array{Array{Float64,2},1}	#distance matrices KLD & comp_mix
	weight::Array{Array{Float64,2},1}	#distance matrices KLD & comp_mix
	geom::Array{Array{Geometric{Float64},1},1}
	param::Array{Array{Float64,1},1}
	KLD::Array{Array{Float64,2},1}	#distance matrices KLD & comp_mix

	K::Array{Int64,1}	#number of molecules in segment (according to ground_truth)
##Function to initialize structure
	function ModelSR(SEG::SuperResolution.SegmentSR,idx::Tuple{Int64,Int64},WKLD_modINPUT=1) #idx from list SEG.idx[j]
		obj=new()
		obj.idx=idx
		obj.frame=SEG.frame[CartesianIndex(idx)]
		#obj.modus=["KLD","WKLD"][sum(WKLD_modINPUT==1)+1]	#WKLD_modINPUT==1.."WKLD", else "KLD"
		obj.modus=WKLD_modINPUT
		obj.dim=size(SEG.all_loc,2)
		#obj.dim=size(SEG.all_loc[:,1:(end-1)])[2]	#all_loc last col is framenumbers
		#obj.loc=Array{Array{Float64,2},1}(undef,SEG.N[CartesianIndex(idx)]) #
		obj.member=Array{Array{Int64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.clus_ID=Array{Array{Int64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.mol_ID=SEG.mol_ID[CartesianIndex(idx)]
		#obj.model=Array{MixtureModel,1}(undef,SEG.N[CartesianIndex(idx)])
		obj.cov=Array{Array{Array{Float64,2},1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.comp_mix=Array{Array{Float64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.wkld=Array{Array{Float64,2},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.KLD=Array{Array{Float64,2},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.weight=Array{Array{Float64,2},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.geom=Array{Array{Geometric{Float64},1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.param=Array{Array{Float64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.K=Array{Int64,1}(undef,SEG.N[CartesianIndex(idx)])

		obj.member[1]=collect(1:SEG.N[CartesianIndex(idx)])#Array{Array{Int64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.clus_ID[1]=collect(1:SEG.N[CartesianIndex(idx)])#Array{Array{Int64,1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.cov[1]=repeat([1.0*eye(obj.dim)],SEG.N[CartesianIndex(idx)])#Array{Array{Array{Float64,2},1},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.weight[1]=ones(SEG.N[CartesianIndex(idx)],SEG.N[CartesianIndex(idx)])
		obj.comp_mix[1]=1/(SEG.N[CartesianIndex(idx)]-1)*ones(1,SEG.N[CartesianIndex(idx)])[1,:]#Array{Array{Float64,2},1}(undef,SEG.N[CartesianIndex(idx)])
		obj.geom[1]=repeat([Geometric(0.999)],SEG.N[CartesianIndex(idx)])
		sKLD=[dim2adj(SEG.loc[CartesianIndex(idx)]),]	#dim2adj: 452×3 Array{Float64,2} --> 452-element Array{Array,1}
		#N=size(sKLD[1])[1]

		#global i_inc=1
		#sKLD[1],obj.comp_mix[1],obj.cov[1],obj.clus_ID[1],obj.member[1],obj.wkld[1],obj.weight[i_inc+1],obj.geom[1]=dim2WKLDfull(sKLD[i_inc],obj.cov[i_inc],obj.comp_mix[i_inc],SEG.frame[CartesianIndex(idx)],obj.clus_ID[i_inc],obj.member[i_inc])
		# sKLD[1],ωall[1],Σall[1],Aall[1],Acoll[1],KLDeval[1]=dim2WKLDfull(sKLD[i_inc],Σall[i_inc],ωall[i_inc],Aall[i_inc],Acoll[i_inc])
##1..Extract needed features
		global obj.loc=[dim2adj(dim2adj(SEG.loc[CartesianIndex(idx)])),] #x,y coordinates of all observations [MANDATORY FOR MODEL SELECTION/EVALUATION]
		obj.model=[dim2mix(SEG.loc[CartesianIndex(idx)][:,1:2]),]
		obj.wkld[1]=deepcopy(SEG.loc[CartesianIndex(idx)])

		###4.. KLClustering
		print("number of observations = $(SEG.N[CartesianIndex(idx)])\n")
		@time(
			@showprogress(for i_inc=2:SEG.N[CartesianIndex(idx)]#while size(obj.clus_ID[i_inc],1)>1
				tmpsKLD,obj.comp_mix[i_inc],obj.cov[i_inc],obj.clus_ID[i_inc],obj.member[i_inc],obj.wkld[i_inc],obj.weight[i_inc],obj.geom[i_inc]=dim2WKLDfull(
				# sKLD[i_inc-1],obj.cov[i_inc-1],obj.comp_mix[i_inc-1],SEG.frame[CartesianIndex(idx)],obj.clus_ID[i_inc-1],obj.member[i_inc-1],WKLD_modINPUT)
				sKLD[i_inc-1],obj.cov[i_inc-1],obj.comp_mix[i_inc-1],SEG.frame[CartesianIndex(idx)],obj.clus_ID[i_inc-1],obj.member[i_inc-1],WKLD_modINPUT)
				push!(sKLD,tmpsKLD)#x,y values + pdf value of Framenumber
				push!(obj.loc,dim2adj(sKLD[i_inc]))#x,y values from sKLD variable
				push!(obj.model,dim2mix(dim2adj(sKLD[i_inc]),obj.comp_mix[i_inc]))#MixtureModel (PDF)
				#global i_inc+=1
			end)
		)
		for j=1:(length(obj.geom))
			tete=map(i->params(obj.geom[j][i])[1],1:length(obj.geom[j]))
			obj.param[j]=tete
		end
		obj.K=map(i->length(unique(obj.member[i])),1:(length(obj.member)))
		obj
	end
end
export ModelSR
##
# #x,y coordinates of cluster centers, for each iteration
# obj.loc=map(i->obj.loc[i][:,1:2],1:length(obj.loc))
#
# ###5.. Cluster Evaluation
# #sMMsim=dim2mix(𝛍sim,𝛚sim)#GMM-Distribution based on mean ("RealLoc") mol_IDs, for comparison
#
# print("##########################\n")
# print("Compute Model Selection...\n")
# print("##########################\n")
#
# xRangeSel=collect(1:(i_inc))
# Knum=map(i->size(obj.loc[i],1),collect(1:size(obj.loc,1))) #number of clusters for each interation
#
# #compute number, index and locations of unique molecules
# Acoll_unique=map(i->length(unique(Acoll[i])),1:length(Acoll))
# #print("Acoll_unique:\n"); show(Acoll_unique); print("\n")
# #compute min|max|mean KLD values
# #2:length(obj.wkld) because the 1st entry is [Y,] !!!!
# minKLDeval=map(i->minimum(obj.wkld[i][1e-10 .<obj.wkld[i]]),1:length(obj.wkld)-1)
# push!(minKLDeval,sort(obj.wkld[end-2][minKLDeval[end-1].<obj.wkld[end-2]])[1]) #pseudo minKLD value for 1 Cluster
#
# #compute absolute minKLD changes
# abs_minKLDeval=abs.(minKLDeval[1:end-1]-minKLDeval[2:end])
# push!(abs_minKLDeval,0)
# abs_minKLDeval[2:end]=abs_minKLDeval[1:end-1]
# abs_minKLDeval[1]=0
# #print("abs_minKLDeval:\n"); show(abs_minKLDeval); print("\n")
#
# #obj.loc[end-solvIDXend] is the clustering solution
# global solvIDX=findall(abs_minKLDeval.==maximum(abs_minKLDeval))[1] #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# global solvNumClust=size(obj.loc[solvIDX],1)
# print("estimated cluster for partition $(partN)_$(partN2): $(solvNumClust)\n")
# # print("length(minKLDeval)=$(length(minKLDeval))\n")
# # print("length(abs_delta_minKLDeval)=$(length(abs_minKLDeval))\n")

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>5.0, :alg=>TMJets(orderT=8, orderQ=2, maxsteps=1_000)),
    "X0"   => Hyperrectangle(low=[4.8, 1.8], high=[5.2, 2.2]),
    "info" => Dict(
        "name" => "lotka_volterra",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using Dispersal, Test, Random

Random.seed!(1234)

init = [0 0 0 0 0;
        0 0 0 0 0;
        0 0 1 0 0;
        0 0 0 0 0;
        0 0 0 0 0]

rules = Ruleset(JumpDispersal(prob_threshold=1.0, spotrange=1))
output = ArrayOutput(init; tspan=1:100)
sim!(output, rules)

output[2] =
    [0 0 0 0 0;
     0 0 0 0 0;
     0 0 1 0 0;
     0 0 0 1 0;
     0 0 0 0 0]

output[3] =
    [0 0 0 0 0;
     0 0 0 0 0;
     0 0 1 0 0;
     0 0 1 1 0;
     0 0 1 0 0]

output[end] =
    [1 1 1 1 1;
     1 1 1 1 1;
     1 1 1 1 1;
     1 1 1 1 1;
     1 1 1 1 1]


text:
include("mnist.jl")

function get_args()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--seed"
        arg_type = Int
        default = 0
        "--log"
        arg_type = String
        required = true
        "--ea"
        arg_type = String
        required = true
        "--chromosome"
        arg_type = String
        required = true
        "--cfg"
        arg_type = String
        required = true
    end

    CGP.Config.add_arg_settings!(s)
end

args = parse_args(get_args())
CGP.Config.init(args["cfg"])
srand(args["seed"])
Logging.configure(filename=args["log"], level=INFO)
ea = eval(parse(args["ea"]))
ctype = eval(parse(args["chromosome"]))

maxfit, best = ea(ctype, 4, 2, stdp_mnist)
Logging.info(@sprintf("E%0.6f", -maxfit))


text:
"""
    Certify(optimizer)

Certify uses semidefinite programming to compute over-approximated certificates for a neural network with only one hidden layer.

# Problem requirement
1. Network: one hidden layer, any activation that is differentiable almost everywhere whose derivative is bound by 0 and 1
2. Input: hypercube
3. Output: halfspace

# Return
`BasicResult`

# Method
Semindefinite programming.
Default `optimizer` is `SCSSolver()`.

# Property
Sound but not complete.

# Reference
[A. Raghunathan, J. Steinhardt, and P. Liang,
"Certified Defenses against Adversarial Examples,"
*ArXiv Preprint ArXiv:1801.09344*, 2018.](https://arxiv.org/abs/1801.09344)
"""
@with_kw struct Certify{O<:AbstractMathProgSolver}
    optimizer::O  = SCSSolver()
end

# can pass keyword args to optimizer
# Certify(optimizer::DataType = SCSSolver; kwargs...) =  Certify(optimizer(;kwargs...))

function solve(solver::Certify, problem::Problem)
    @assert length(problem.network.layers) == 2 "Certify only handles Networks that have one hidden layer. Got $(length(problem.network.layers)) total layers"
    model = Model(solver = solver.optimizer)
    c, d = tosimplehrep(problem.output)
    v = c * problem.network.layers[2].weights
    W = problem.network.layers[1].weights
    M = get_M(v[1, :], W)
    n = size(M, 1)
    P = @variable(model, [1:n, 1:n], SDP)
    # Compute value
    output = c * compute_output(problem.network, problem.input.center) .- d[1]
    epsilon = problem.input.radius[1]
    o = output + epsilon/4 * tr(M*P)
    # Specify problem
    @constraint(model, diag(P) .<= ones(n))
    @objective(model, Max, o[1])
    status = solve(model, suppress_warnings = true)
    return interpret_result(solver, status, o[1])
end

# True if o < 0
# Undertermined if otherwise
function interpret_result(solver::Certify, status, o)
    # println("Upper bound: ", getvalue(o[1]))
    if getvalue(o) <= 0
        return BasicResult(:SAT)
    else
        return BasicResult(:Unknown)
    end
end

# M is used in the semidefinite program
function get_M(v::Vector{Float64}, W::Matrix{Float64})
    m = W' * Diagonal(v)
    mxs, mys = size(m)
    o = ones(size(W, 2), 1)
    # TODO Mrow2 and Mrow3 look suspicious (dims don't match in the general case).
    Mrow1 = [zeros(1, 1+mxs)    o'*m]
    Mrow2 = [zeros(mxs, 1+mxs)     m]
    Mrow3 = [m'*o  m' zeros(mys, mys)]

    M = [Mrow1; Mrow2; Mrow3]
    return M
end


text:
"""
    install{S<:AbstractString}(config::Config, version::VersionNummber; labels::Array{S}=String[])

Downloads the latest binary version for your platform and symlinks the binary with the appropriate labels.
"""
function install{S<:AbstractString}(config::Config, version::VersionNumber; labels::Array{S}=String[])
    init(config)

    # download the julia version
    download_url = julia_url(version)
    base_name = "julia-$(version.major).$(version.minor)_$(Dates.today())"
    tmp_dest = join(config.tmp, base_name)

    info(logger, "Downloading julia $(version.major).$(version.minor) from $download_url ...")
    download(download_url, tmp_dest, false)

    # Perform some verification on the download
    if stat(tmp_dest).size < 1024
        # S3 responds with an error message when a URL doesn't exist
        unavailable = open(tmp_dest, "r") do f
            @mock contains(readstring(f), "key does not exist")
        end

        if unavailable
            error(logger, "Julia version $version does not exist")
        else
            error(logger, "Aborting install. Download appears corrupt")
        end
    end

    bin_path = Playground.install_bin(tmp_dest, config, base_name, false)
    Playground.link_julia(bin_path, config, labels)

    label_path = join(config.bin, "julia-$(version.major).$(version.minor)")
    debug(logger, "Creating link $label_path -> $bin_path")
    symlink(bin_path, label_path; exist_ok=true, overwrite=true)
end

"""
    install{S<:AbstractString}(config::Config, executable::AbstractPath; labels::Array{S}=String[])

Creates symlinks from an existing julia install.
"""
function install{S<:AbstractString}(config::Config, executable::AbstractPath; labels::Array{S}=String[])
    info(logger, "Adding julia labels $labels to $executable")
    if exists(executable)
        init(config)

        exe = abs(executable)
        while islink(exe)
            exe = join(parent(exe), readlink(exe))
        end

        exe = abs(exe)
        if isfile(exe)
            Playground.link_julia(exe, config, labels)
        else
            error(logger, "$exe is not a valid executable.")
        end
    else
        error(logger, "$executable is not a valid path")
    end
end

# @doc doc"""
#    This option
#         1. clones the git repo from the provided url.
#         2. checks out the supplied revision into the local branch_name.
#         3. attempts to build julia from scratch.
# """ ->
# function gitinstall(config::Config; url::ASCIIString="", dir::AbstractString="", revision="", labels=[])

#     error("Installing from a git repo isn't implemented yet.")

#     info("Cloning the julia repository into the playground")
#     run(`git clone $(url) $(name)` |> gitlog)

#     # Handle the cd into and out of src directory cause cd() in base
#     # seems to be broken.
#     cwd = pwd()
#     cd(dest)
#     build_julia(julia, , log_path)
#     cd(cwd)
# end

function link_julia{S<:AbstractString}(bin_path::AbstractPath, config::Config, labels::Array{S}=[])
    bin_path = abs(bin_path)
    for label in labels
        label_path = join(config.bin, label)
        debug(logger, "Creating link $label_path -> $bin_path")
        symlink(bin_path, label_path; exist_ok=true, overwrite=true)
    end

    @compat if is_unix()
        ret = readstring(`$(bin_path) -e 'versioninfo()'`)
        lines = split(ret, "\n")

        versionstr = split(lines[1], " ")[3]
        debug(logger, "Creating link julia-$(versionstr) -> $bin_path")
        symlink(bin_path, join(config.bin, "julia-$(versionstr)"); exist_ok=true, overwrite=true)

        commit_sha = strip(split(lines[2], " ")[2], '*')
        debug(logger, "Creating link julia-$(commit_sha) -> $bin_path")
        symlink(bin_path, join(config.bin, "julia-$(commit_sha)"); exist_ok=true, overwrite=true)
    end
end

# function build_julia(target, root_path, log_path)
#     info("Building julia ( $(target) )...")
#     gitlog = joinpath(log_path, "git.log")
#     buildlog = joinpath(log_path, "build.log")

#     run(`git checkout $(target)` >> gitlog)
#     info("checking out $(target)")

#     # Write the different prefix to the Make.user file before
#     # building and installing.
#     info("setting prefix in Make.user")
#     fstrm = open("Make.user","w")
#     write(fstrm, "prefix=$(root_path)")

#     info("Building julia")
#     # Build and install.
#     # TODO: log the build output properly in root_dir/log
#     run(`make` |> buildlog)
#     run(`make install` >> buildlog)
#     println("Julia has been built and installed.")
# end

@compat if is_apple()
    function install_bin(src::AbstractPath, config::Config, base_name, force)
        src_path = abs(join(config.src, base_name))
        bin_path = abs(join(config.bin, base_name))
        exe_path = abs(join(src_path, "Contents/Resources/julia/bin/julia"))

        function install_dmg(mountdir::AbstractPath)
            app_path = nothing
            try
                debug(logger, "Mounting $src @ $mountdir")
                run(`hdiutil attach -mountpoint $mountdir $src`)
                for f in readdir(mountdir)
                   if extension(f) == "app"
                        app_path = join(mountdir, f)
                        break
                   end
                end
                if app_path != nothing
                    debug(logger, "Copying $app_path to $src_path")
                    copy(app_path, src_path)
                    chmod(exe_path, Playground.JULIA_BIN_MODE)
                    debug(logger, "Permissions set to $(Playground.JULIA_BIN_MODE)")
                end
            finally
                debug(logger, "Detaching $src from $mountdir")
                run(`hdiutil detach $mountdir`)
            end
        end

        dmg_tmp_dir = mktmpdir(config.tmp)
        try
            # Don't bother running this if src_path already exists
            if !exists(src_path) || force
                info(logger, "Installing $src ...")
                install_dmg(dmg_tmp_dir)
            end

            debug(logger, "Creating link $bin_path -> $exe_path")
            symlink(exe_path, bin_path; exist_ok=true, overwrite=true)
        finally
            remove(dmg_tmp_dir)
        end

        return bin_path
    end
elseif is_unix()
    function install_bin(src::AbstractPath, config::Config, base_name, force)
        src_path = abs(join(config.src, base_name))
        bin_path = abs(join(config.bin, base_name))

        info(logger, "Installing $src ...")
        if !exists(src_path) || force
            mkdir(src_path; recursive=true)
            try
                debug(logger, "Extracting $src to $src_path")
                run(`tar -xzf $src -C $src_path`)
            catch
                remove(src_path, recursive=true)
            end
        end

        julia_bin_path = join(src_path, readdir(src_path)[1], p"bin/julia")
        chmod(julia_bin_path, Playground.JULIA_BIN_MODE)
        debug(logger, "Permissions set to $(Playground.JULIA_BIN_MODE)")
        debug(logger, "Creating link $bin_path -> $julia_bin_path")
        symlink(julia_bin_path, bin_path; exist_ok=true, overwrite=true)

        return bin_path
    end
elseif is_windows()
    function install_bin(src::AbstractPath, config::Config)
        # not sure what to do here yet.
        error(logger, "installing Julia EXEs on Windows not implemented yet.")
    end
end


text:
module Dec4Part1
ref_key = ["byr","iyr","eyr","hgt","hcl","ecl","pid","cid"]
ref_keys = filter(x -> x ≠ "cid", ref_key)

function read_passport(file_path::String)::String
    return read(file_path, String)
end

function extract_chunks(x::String, pattern)::Array{String,1}
    match_iterator = eachmatch(pattern, x)
    return map(x -> x.match, match_iterator)
end

function extract_key_from_chunk(x::String, delimiter=":")::String
    substring = split(x, delimiter)
    string = String.(substring)
    return string[1]
end

function extract_keys(chunk::String)::Array{String,1}
    chunks = extract_chunks(chunk, r"((?:[^\s])+)")
    keys_chunks = extract_key_from_chunk.(chunks)
    return keys_chunks
end

function check_validity(keys::Array{String,1})::Bool
    global ref_keys
    for key in ref_keys
        if key ∉ keys
            return false
        end
    end
    return true
end

function get_valid_passports(passports::String)::Int64
    passports_array = extract_chunks(passports, r"(?s)((?:[^\n][\n]?)+)")
    passport_keys = extract_keys.(passports_array)
    bools = check_validity.(passport_keys)
    return sum(bools)
end

function main(file_path::String)::Int64
    passports = read_passport(file_path)
    return get_valid_passports(passports)
end
end


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

# this is just one big texture
texture_atlas = loadasset("doge.png")
w, h = size(texture_atlas)
const n = 40
xrange = linspace(0, w, n)
yrange = linspace(0, h, n)
scale  = Vec2f0(step(xrange), step(yrange))

# position in a grid
positions = map(timesignal) do t
    vec(Point2f0[(x+(sin(t*2*pi)*400),y+(sin(0+y*t*0.01)*200)+(cos(t*2*pi)*200)) for x=xrange, y=yrange])
end

# uv coordinates are normalized coordinates into the texture_atlas
# they need the start point and the width of each rectangle (sprites are rectangles)
# so you will not actually index with the circle primitive, but rather with
# with the quad of the particle (the rest of the quad is transparent)
# note, that for uv coordinates, the origin is on the top left corner
uv_offset_width = vec(Vec4f0[(x,y,x+(1/n),y+(1/n)) for x=linspace(0, 1, n), y=linspace(1, 0, n)])

# when position and scale are defined, We can leave the middle and radius of
# Circle undefined, so just passing the type.
distfield = visualize((Circle, positions),
    scale=scale,
    stroke_width=1f0,
    uv_offset_width=uv_offset_width,
    stroke_color = RGBA{Float32}(0.9,0.9,0.9,1.0),
    image=texture_atlas
)
_view(distfield, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
import SparseArrays: nnz

export AbstractFspVector, FspVectorSparse, sum, get_states, get_values, get_statedict , nnz

abstract type AbstractFspVector end

"""
    FspVectorSparse{NS, IntT<:Integer, RealT<:AbstractFloat} <: AbstractFspVector

Sparse multi-indexed vector (i.e., sparse "tensor") with entries of type `RealT` and multi-indices of type `IntT` for storing CME solution and sensitivity vectors output by the Finite State Projection.
"""
Base.@kwdef mutable struct FspVectorSparse{NS,IntT<:Integer,RealT<:AbstractFloat} <: AbstractFspVector
    states::Vector{MVector{NS,IntT}}
    values::Vector{RealT}
    state2idx::Dict{MVector{NS,IntT},IntT}
end
get_states(v::FspVectorSparse) = v.states
get_values(v::FspVectorSparse) = v.values
get_statedict(v::FspVectorSparse) = v.state2idx
nnz(p::FspVectorSparse) = length(get_states(p))

"""
    FspVectorSparse(states::Vector{MVector{NS, IntT}}, values::Vector{RealT})

Construct a sparse representation of a FSP-truncated probability distribution with support `states` and proability values `values`.
"""
function FspVectorSparse(states::Vector{MVector{NS,IntT}}, values::Vector{RealT}; checksizes::Bool = true) where {NS,IntT<:Integer,RealT<:AbstractFloat}
    checksizes && (length(states) ≠ length(values)) && throw(ArgumentError("State and value lists must have equal lengths."))
    return FspVectorSparse{NS,IntT,RealT}(
        states = deepcopy(states),
        values = deepcopy(values),
        state2idx = Dict([x => i for (i, x) in enumerate(states)])
    )
end

"""
    FspVectorSparse(statespace::AbstractStateSpaceSparse{NS, NR, IntT}, statevalpairs::Vector{Pair{VecT, RealT}})

Construct a sparse representation of a FSP-truncated probability distribution on the state space `statespace` and proability values `values`.
"""
function FspVectorSparse(statespace::AbstractStateSpaceSparse{NS,NR,IntT}, statevalpairs::Vector{Pair{VecT,RealT}}) where {NS,NR,IntT<:Integer,VecT<:AbstractVector,RealT<:AbstractFloat}
    states = deepcopy(get_states(statespace))
    state2idx = deepcopy(get_statedict(statespace))
    values = zeros(RealT, get_state_count(statespace))
    for (x, v) in statevalpairs
        idx = get(state2idx, x, 0)
        (idx ≠ 0) && (values[idx] = v)
    end
    return FspVectorSparse{NS,IntT,RealT}(
        states = states,
        values = values,
        state2idx = state2idx
    )
end

function Base.sum(p::FspVectorSparse)
    Base.sum(p.values)
end

"""
    sum(p::FspVectorSparse, dims)

Reduce the FSP vector `p` by summing over species specified in `dims` and return a new FSP vector over a reduced-dimension state space.
"""
function Base.sum(p::FspVectorSparse{NS,IntT,RealT}, dims::AbstractVector{<:Integer})::FspVectorSparse where {NS,IntT<:Integer,RealT<:AbstractFloat}
    dims = unique(dims)
    !((min(dims...) >= 1) && (max(dims...) <= NS)) && throw(ArgumentError("Input dimensions must be between 1 and $(NS)."))

    keepdims = setdiff(1:NS, dims)

    full_states = get_states(p)
    full_vals = get_values(p)

    reduced_species_count = NS - length(dims)
    reduced_states = Vector{MVector{reduced_species_count,IntT}}()
    reduced_state2idx = Dict{MVector{reduced_species_count,IntT}, IntT}()
    reduced_vals = Vector{RealT}()

    rvec_length = 0
    for (i, eachstate) in enumerate(full_states)
        reduced_state = eachstate[keepdims]
        reduced_idx = get(reduced_state2idx, reduced_state, 0)
        if reduced_idx == 0
            rvec_length += 1
            reduced_state2idx[reduced_state] = rvec_length
            push!(reduced_states, copy(reduced_state))
            push!(reduced_vals, full_vals[i])
        else
            reduced_vals[reduced_idx] += full_vals[i]
        end
    end

    return FspVectorSparse{reduced_species_count,IntT,RealT}(
        states = reduced_states,
        values = reduced_vals,
        state2idx = reduced_state2idx
    )
end

"""
    Array(p::FspVectorSparse)

Construct a new dense `N`-dimensional array from the FSP vector.
"""
function Base.Array(p::FspVectorSparse{NS,IntT,RealT}) where {NS, IntT<:Integer, RealT<:AbstractFloat}
    states = get_states(p)
    if isempty(states)
        throw(ArgumentError("Cannot construct dense array from empty FspVector instance."))
    end
    values = get_values(p)

    dimcount = length(states[1])
    bounds = zeros(Int, dimcount)
    for eachstate in states
        for i in 1:dimcount
            bounds[i] = max(bounds[i], eachstate[i]+1)
        end
    end

    out = zeros(RealT, bounds...)
    tmp = MVector{NS,IntT}([0 for i in 1:NS])
    for i in 1:length(values)
        tmp .= states[i] .+ 1
        out[tmp...] = values[i]
    end
    return out
end

using Printf
function Base.show(io::IO, vec::FspVectorSparse)
    nz = length(vec.values)
    if nz == 0
        println(io, "Sparse multi-indexed vector with no nonzero entries.")
        return nothing
    end
    println("Sparse multi-indexed vector with $(nz) nonzero entries:")
    for i in 1:min(nz, 5)
        print(io, "$(vec.states[i]) ↦ ")
        @printf(io, "%.2e \n", vec.values[i])
    end
    (nz > 10) && println(repeat(" ", length(vec.states[1]) + 5), "...")
    for i in max(nz - 5, 6):nz
        print(io, "$(vec.states[i]) ↦ ")
        @printf(io, "%.2e \n", vec.values[i])
    end
    nothing
end