text:
function populatekernel!(::Flow, K, Kₗₘ, Y, modes_lm, K_all_modes, Ω_ind, temp⁰, temp⁺¹)

	# The kernel is defined as
	# Kr = Σ_l K₀,ₗ₀*(r) Yₗ₀ + 2 Σ_{m=1}^{l} ℜ[K₀,ₗₘ* Yₗₘ]
	# Kθ = -√2 ℜ[ Σ_l K₁,ₗ₀*(r) Y¹ₗ₀ + Σ_{m=1}^{l} (K₁,ₗₘ* Y¹ₗₘ - K₋₁,ₗₘ* Y⁻¹ₗₘ) ]
	# Kϕ = √2 ℑ[ Σ_l K₁,ₗ₀*(r) Y¹ₗ₀ + Σ_{m=1}^{l} (K₁,ₗₘ* Y¹ₗₘ - K₋₁,ₗₘ* Y⁻¹ₗₘ) ]

	_modes_lm = SphericalHarmonicModes.ofordering(K_all_modes, modes_lm)
	for (lm_ind, (l,m)) in enumerate(_modes_lm)

		# Positive and negative m's are summed over, we don't need to loop over negative m's
		# separately
		m < 0 && continue

		Ylm = Y[(l,m)]

		@views begin

		if iszero(m)
			@. temp⁰ = real(conj(Kₗₘ[.., 0, lm_ind])*Ylm[0])
			@. temp⁺¹ = conj(Kₗₘ[.., 1, lm_ind])*Ylm[1]
		else
			# Sum of the +m and -m terms
			@. temp⁰ = 2real(conj(Kₗₘ[.., 0, lm_ind])*Ylm[0])
			@. temp⁺¹ = conj(Kₗₘ[.., 1, lm_ind])*Ylm[1] - Kₗₘ[..,-1, lm_ind]*conj(Ylm[-1])
		end

		end# views

		# r-component
		@. K[.., 1, Ω_ind] += temp⁰
		# θ-component
		@. K[.., 2, Ω_ind] += -√2*real(temp⁺¹)
		# ϕ-component
		@. K[.., 3, Ω_ind] += √2*imag(temp⁺¹)
	end
	return K
end

function Ku_longitudinal_slice_loop(localtimer, ϕθ::ProductSplit,
	Kₗₘ_fits_filename, modes_lm::LM, K_all_modes::LM)

	l_min, l_max = extrema(l_range(modes_lm))
	m_min, m_max = extrema(m_range(modes_lm))

	lm_ind_min = modeindex(K_all_modes, l_min, m_min)
	lm_ind_max = modeindex(K_all_modes, l_max, m_max)

	Kₗₘ = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename), "r") do f
		hdu = f[1]::ImageHDU{Float64,3}
		arr = reinterpret(ComplexF64, read(hdu))[.., lm_ind_min:lm_ind_max]
		OffsetArray(arr, :, -1:1, lm_ind_min:lm_ind_max)
	end

	K = zeros(size(Kₗₘ, 1), 3, length(ϕθ)) # 2nd axis is for vector components (r̂, ̂θ, ̂ϕ)

	temp⁺¹ = zeros(ComplexF64, nr)
	temp⁰ = zeros(nr)

	S = VectorSphericalHarmonics.cache(l_max)

	for (Ω_ind,(ϕ,θ)) in enumerate(ϕθ)
		VectorSphericalHarmonics.cache!(S, θ, ϕ)
		Y = genspharm(modes_lm, θ, ϕ, S)

		populatekernel!(Flow(), K, Kₗₘ, Y, modes_lm,
			K_all_modes, Ω_ind, temp⁰, temp⁺¹)
	end

	return K
end

function Ku_longitudinal_slice_from_Kₗₘ(comm, Kₗₘ_fits_filename; kwargs...)

	header = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename), "r") do f
		read_header(f[1])
	end

	K_all_modes = LM(header["L_MIN"]:header["L_MAX"],
					header["M_MIN"]:header["M_MAX"])

	l_min = get(kwargs,:l_min, header["L_MIN"]); header["L_MIN"] = l_min;
	l_max = get(kwargs,:l_max, header["L_MAX"]); header["L_MAX"] = l_max;
	m_min = get(kwargs,:m_min, min(header["M_MIN"], l_max)); header["M_MIN"] = m_min;
	m_max = get(kwargs,:m_max, min(header["M_MAX"], l_max)); header["M_MAX"] = m_max;
	kernel_modes = LM(l_min:l_max, m_min:m_max)
	nl = length(l_min:l_max)

	nθ_default = max(360, 4header["L_MAX"])
	nθ = get(kwargs,:nθ, nθ_default)
	θ_default = LinRange(0,π, nθ)
	θ = get(kwargs,:θ,θ_default)
	nθ = length(θ)
	ϕ = get(kwargs,:ϕ,(header["PHI1"] + header["PHI2"])/2)
	nl = length(l_range(kernel_modes))

	Krθ = pmapreduce(comm, Ku_longitudinal_slice_loop, (x...) -> cat(x..., dims = 3),
		(ϕ:ϕ,θ), Kₗₘ_fits_filename, kernel_modes, K_all_modes)

	Krθ === nothing && return nothing

	Krθ = permutedims(Krθ, [1, 3, 2]) # last axis will store the vector components

	# Write to fits file

	header["THSLMN"] = minimum(θ)
	set_comment!(header,"THSLMN","Minimum colatitude in slice")

	header["THSLMX"] = maximum(θ)
	set_comment!(header,"THSLMX","Maximum colatitude in slice")

	header["NTH"] = nθ
	set_comment!(header,"NTH","Number of grid points in colatitude")

	header["PHSL"] = ϕ
	set_comment!(header,"PHSL","Azimuth at which the slice is computed")

	FITS(joinpath(SCRATCH_kerneldir[],"Ku_rθslice_from_Klm.fits"),"w") do f
		FITSIO.write(f, Krθ, header=header)
	end
	Krθ
end

function Ku_latitudinal_slice_loop(localtimer, ϕθ::ProductSplit,
	Kₗₘ_fits_filename, modes_lm::LM, K_all_modes::LM)

	l_min, l_max = extrema(l_range(modes_lm))
	m_min, m_max = extrema(m_range(modes_lm))

	lm_ind_min = modeindex(K_all_modes, l_min, m_min)
	lm_ind_max = modeindex(K_all_modes, l_max, m_max)

	Kₗₘ = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename),"r") do f
		hdu = f[1]::ImageHDU{Float64,3}
		arr = reinterpret(ComplexF64, read(hdu))[.., lm_ind_min:lm_ind_max]
		OffsetArray(arr, axes(arr, 1), -1:1, lm_ind_min:lm_ind_max)
	end

	K = zeros(size(Kₗₘ, 1), 3, length(ϕθ)) # 2nd axis is for vector components (r,θ,ϕ)

	temp⁺¹ = zeros(ComplexF64, nr)
	temp⁰ = zeros(nr)

	S = VectorSphericalHarmonics.cache(l_max)

	for (Ω_ind,(ϕ,θ)) in enumerate(ϕθ)
		VectorSphericalHarmonics.cache!(S, θ, ϕ)
		Y = genspharm(modes_lm, θ, ϕ, S)

		populatekernel!(Flow(), K, Kₗₘ, Y, modes_lm,
			K_all_modes, Ω_ind, temp⁰, temp⁺¹)
	end

	return K
end

function Ku_latitudinal_slice_from_Kₗₘ(comm, Kₗₘ_fits_filename; kwargs...)
	header = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename), "r") do f
		read_header(f[1])
	end

	K_all_modes = LM(header["L_MIN"]:header["L_MAX"],
					header["M_MIN"]:header["M_MAX"])

	l_min = get(kwargs,:l_min, header["L_MIN"]); header["L_MIN"] = l_min;
	l_max = get(kwargs,:l_max, header["L_MAX"]); header["L_MAX"] = l_max;
	m_min = get(kwargs,:m_min, min(header["M_MIN"], l_max)); header["M_MIN"] = m_min;
	m_max = get(kwargs,:m_max, min(header["M_MAX"], l_max)); header["M_MAX"] = m_max;
	kernel_modes = LM(l_min:l_max, m_min:m_max)
	nl = length(l_min:l_max)

	θ = get(kwargs,:θ,(header["TH1"] + header["TH2"])/2)
	nϕ_default = max(720, 4header["L_MAX"])
	nϕ = get(kwargs,:nϕ, nϕ_default)
	ϕ_default = LinRange(0, 2π, nϕ)
	ϕ = get(kwargs,:ϕ,ϕ_default)
	nϕ = length(ϕ)
	nl = length(l_range(kernel_modes))

	Krϕ = pmapreduce(comm, Ku_latitudinal_slice_loop, (x...) -> cat(x..., dims = 3), (ϕ, θ:θ),
		Kₗₘ_fits_filename, kernel_modes, K_all_modes)

	Krϕ === nothing && return nothing

	Krϕ = permutedims(Krϕ,[1, 3, 2]) # last axis will store the vector components

	header["PHSLMN"] = minimum(ϕ)
	set_comment!(header,"PHSLMN","Minimum azimuth in slice")

	header["PHSLMX"] = maximum(ϕ)
	set_comment!(header,"PHSLMX","Maximum azimuth in slice")

	header["NPHI"] = nϕ
	set_comment!(header,"NPHI","Number of grid points in azimuth")

	header["THSL"] = θ
	set_comment!(header,"THSL","Colatitude at which the slice is computed")

	FITS(joinpath(SCRATCH_kerneldir[], "Ku_rϕslice_from_Klm.fits"),"w") do f
		FITSIO.write(f, Krϕ, header=header)
	end
	Krϕ
end

function Ku_surface_slice_loop(localtimer, ϕθ::ProductSplit,
	Kₗₘ_fits_filename, modes_lm::LM, K_all_modes::LM, r_obs_ind)

	l_min, l_max = extrema(l_range(modes_lm))
	m_min, m_max = extrema(m_range(modes_lm))

	lm_ind_min = modeindex(K_all_modes, l_min, m_min)
	lm_ind_max = modeindex(K_all_modes, l_max, m_max)

	Kₗₘ = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename),"r") do f
		hdu = f[1]::ImageHDU{Float64,3}
		arr = reinterpret_as_complex(read(hdu))[r_obs_ind, :, lm_ind_min:lm_ind_max]
		OffsetArray(arr, -1:1, lm_ind_min:lm_ind_max)
	end

	K = zeros(3, length(ϕθ))

	temp⁺¹ = zeros(ComplexF64)
	temp⁰ = zeros()

	V = VectorSphericalHarmonics.VSHCache(Float64, modes_lm)

	for (Ω_ind, (ϕ, θ)) in enumerate(ϕθ)
		Y = genspharm!(V, θ, ϕ)

		populatekernel!(Flow(), K, Kₗₘ, Y, modes_lm,
			K_all_modes, Ω_ind, temp⁰, temp⁺¹)
	end

	return K
end

function Ku_surface_slice_from_Kₗₘ(comm, Kₗₘ_fits_filename; kwargs...)

	r_obs = get(kwargs,:r_obs, r_obs_default)

	r_obs_ind = radial_grid_index(r_obs)

	header = FITS(joinpath(SCRATCH_kerneldir[], Kₗₘ_fits_filename), "r") do f
		read_header(f[1])
	end

	K_all_modes = LM(header["L_MIN"]:header["L_MAX"],
					header["M_MIN"]:header["M_MAX"])

	l_min = get(kwargs,:l_min, header["L_MIN"]); header["L_MIN"] = l_min;
	l_max = get(kwargs,:l_max, header["L_MAX"]); header["L_MAX"] = l_max;
	m_min = get(kwargs,:m_min, min(header["M_MIN"], l_max)); header["M_MIN"] = m_min;
	m_max = get(kwargs,:m_max, min(header["M_MAX"], l_max)); header["M_MAX"] = m_max;
	kernel_modes = LM(l_min:l_max, m_min:m_max)
	nl = length(l_min:l_max)

	nϕ_default = max(720, 4header["L_MAX"])
	nϕ = get(kwargs,:nϕ, nϕ_default)
	ϕ_default = LinRange(0, 2π, nϕ)
	ϕ = get(kwargs,:ϕ,ϕ_default)
	nϕ = length(ϕ)

	nθ_default = max(360, 4header["L_MAX"])
	nθ = get(kwargs,:nθ, nθ_default)
	θ_default = LinRange(0,π, nθ)
	θ = get(kwargs,:θ,θ_default)
	nθ = length(θ)
	nl = length(l_range(kernel_modes))

	K = pmapreduce(comm, Ku_surface_slice_loop, hcat, (ϕ,θ),
			Kₗₘ_fits_filename, kernel_modes, K_all_modes, r_obs_ind)

	K === nothing && return nothing

	Kϕθ = permutedims(reshape(K, 3, nϕ, nθ),[2, 3, 1]) # last axis will store the vector components

	header["PHSLMN"] = minimum(ϕ)
	set_comment!(header,"PHSLMN","Minimum azimuth in slice")

	header["PHSLMX"] = maximum(ϕ)
	set_comment!(header,"PHSLMX","Maximum azimuth in slice")

	header["NPHI"] = nϕ
	set_comment!(header,"NPHI","Number of grid points in azimuth")

	header["THSLMN"] = minimum(θ)
	set_comment!(header,"THSLMN","Minimum colatitude in slice")

	header["THSLMX"] = maximum(θ)
	set_comment!(header,"THSLMX","Maximum colatitude in slice")

	header["NTH"] = nθ
	set_comment!(header,"NTH","Number of grid points in colatitude")

	header["ROBS"] = r_obs
	set_comment!(header,"ROBS","Observation radius")

	FITS(joinpath(SCRATCH_kerneldir[], "Ku_ϕθ_from_Klm.fits"),"w") do f
		FITSIO.write(f, Kϕθ, header=header)
	end
	Kϕθ
end

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

function populatekernel!(::SoundSpeed, K, Kₗₘ, Y, modes_lm, Ω_ind)

	# The kernel is defined as K = Σ_l Kₗ₀(r) Yₗ₀ + 2 Σ_{m=1}^{l} ℜ[Kₗₘ Yₗₘ]

	for (lm_ind,(l, m)) in enumerate(modes_lm)

		# don't need to iterate over negative m's as they are already summed over
		(m < 0) && continue

		# lm_ind = modeindex(modes_lm,(l, m))

		Ylm = Y[(l,m)]

		@views begin
			if iszero(m)
				# Term is purely real
				@. K[.., Ω_ind] += real(Kₗₘ[.., lm_ind])*real(Ylm)
			else
				# Sum of the +t and -t terms
				@. K[.., Ω_ind] += 2real(Kₗₘ[.., lm_ind]*Ylm)
			end
		end
	end
end

function Kc_angular_slice_loop(localtimer, ϕθ::ProductSplit, Klm_fits_filename, st_iterator::LM)

	l_max = maximum(l_range(st_iterator))

	S = SphericalHarmonics.cache(l_max)

	Klm = FITS(joinpath(SCRATCH_kerneldir[], Klm_fits_filename),"r") do f
		hdu = f[1]::ImageHDU{Float64,2}
		reinterpret(ComplexF64, read(hdu))
	end

	K = zeros(size(Klm, 1), length(ϕθ))

	for (Ω_ind,(ϕ,θ)) in enumerate(ϕθ)
		computePlmcostheta!(S, θ)
		Ylm = computeYlm!(S, θ, ϕ)
		populatekernel!(SoundSpeed(), K, Klm, Ylm, st_iterator, Ω_ind)
	end
	return K
end

function Kc_longitudinal_slice_from_Kₗₘ(comm, Kst_fits_filename; kwargs...)

	r_src = get(kwargs,:r_src, r_src_default)

	Gfn_path_src = Gfn_path_from_source_radius(r_src)

	@load(joinpath(Gfn_path_src,"parameters.jld2"),ℓ_arr)

	header = FITS(joinpath(SCRATCH_kerneldir[], Kst_fits_filename), "r") do f
		read_header(f[1])
	end

	l_max = header["L_MAX"]
	m_max = header["M_MAX"]
	l_min = header["L_MIN"]
	m_min = header["M_MIN"]
	st_iterator = LM(l_min:l_max, m_min:m_max)

	θ = get(kwargs,:θ, LinRange(0,π, get(kwargs,:nθ, 4l_max)))
	nθ = length(θ)

	header["THSLMN"] = minimum(θ)
	set_comment!(header,"THSLMN","Minimum colatitude in slice")

	header["THSLMX"] = maximum(θ)
	set_comment!(header,"THSLMX","Maximum colatitude in slice")

	header["NTH"] = nθ
	set_comment!(header,"NTH","Number of grid points in colatitude")

	ϕ = get(kwargs,:ϕ,(header["PHI1"] + header["PHI2"])/2)

	header["PHSL"] = ϕ
	set_comment!(header,"PHSL","Azimuth at which the slice is computed")

	Krθ = pmapreduce(comm, Kc_angular_slice_loop, hcat, (ϕ:ϕ, θ),
			Kst_fits_filename, st_iterator;
			progress_str = "Angles processed : ")

	Krθ === nothing && return nothing

	FITS(joinpath(SCRATCH_kerneldir[], "Krθ_from_Kst.fits"),"w") do f
		FITSIO.write(f, Krθ, header=header)
	end
	Krθ
end

function Kc_latitudinal_slice_from_Kₗₘ(comm, Kst_fits_filename; kwargs...)

	r_src = get(kwargs,:r_src, r_src_default)

	Gfn_path_src = Gfn_path_from_source_radius(r_src)

	@load(joinpath(Gfn_path_src,"parameters.jld2"),ℓ_arr)

	header = FITS(joinpath(SCRATCH_kerneldir[], Kst_fits_filename),"r") do f
		read_header(f[1])
	end

	l_max = header["L_MAX"]
	m_max = header["M_MAX"]
	l_min = header["L_MIN"]
	m_min = header["M_MIN"]
	st_iterator = LM(l_min:l_max, m_min:m_max)

	ϕ = get(kwargs,:ϕ, LinRange(0, 2π, get(kwargs,:nϕ, 4l_max)))
	nϕ = length(ϕ)

	header["PHSLMN"] = minimum(ϕ)
	set_comment!(header,"PHSLMN","Minimum azimuth in slice")

	header["PHSLMX"] = maximum(ϕ)
	set_comment!(header,"PHSLMX","Maximum azimuth in slice")

	header["NPHI"] = nϕ
	set_comment!(header,"NPHI","Number of grid points in azimuth")

	θ = get(kwargs,:θ,(header["TH1"] + header["TH2"])/2)
	header["THSL"] = θ
	set_comment!(header,"THSL","Colatitude at which the slice is computed")

	Krϕ = pmapreduce(comm, Kc_angular_slice_loop, hcat, (ϕ, θ:θ),
			Kst_fits_filename, st_iterator)

	Krϕ === nothing && return nothing
	FITS(joinpath(SCRATCH_kerneldir[],"Krϕ_from_Kst.fits"),"w") do f
		FITSIO.write(f, Krϕ, header=header)
	end
	Krϕ
end

function Kc_surface_slice_loop(localtimer, ϕθ::ProductSplit,
	Klm_fits_filename, modes_lm, r_obs_ind)

	l_max = maximum(l_range(modes_lm))

	S = SphericalHarmonics.cache(l_max)

	Klm = FITS(joinpath(SCRATCH_kerneldir[], Klm_fits_filename),"r") do f
		hdu = f[1]::ImageHDU{Float64,2}
		reinterpret(ComplexF64, read(hdu))[r_obs_ind,:]
	end

	K = zeros(length(ϕθ))

	for (Ω_ind,(ϕ,θ)) in enumerate(ϕθ)
		computePlmcostheta!(S, θ)
		Ylm = computeYlm!(S, θ, ϕ)

		populatekernel!(SoundSpeed(), K, Klm, Ylm, modes_lm,Ω_ind)
	end
	return K
end

function Kc_surface_slice_from_Kₗₘ(comm, Kst_fits_filename; kwargs...)

	r_src = get(kwargs,:r_src, r_src_default)
	r_obs = get(kwargs,:r_obs, r_obs_default)

	r_obs_ind = radial_grid_index(r_obs)

	Gfn_path_src = Gfn_path_from_source_radius(r_src)

	@load(joinpath(Gfn_path_src,"parameters.jld2"),ℓ_arr)

	header = FITS(joinpath(SCRATCH_kerneldir[], Kst_fits_filename),"r") do f
		read_header(f[1])
	end

	l_max = header["L_MAX"]
	m_max = header["M_MAX"]
	l_min = header["L_MIN"]
	m_min = header["M_MIN"]
	st_iterator = LM(l_min:l_max, m_min:m_max)

	ϕ = get(kwargs,:ϕ, LinRange(0, 2π, get(kwargs,:nϕ, 4l_max)))
	nϕ = length(ϕ)
	θ = get(kwargs,:θ, LinRange(0, π, get(kwargs,:nθ, 4l_max)))
	nθ = length(θ)

	K = pmapreduce(comm, Kc_surface_slice_loop, vcat, (ϕ,θ),
			Kst_fits_filename, st_iterator, r_obs_ind)

	K === nothing && return nothing

	Kϕθ = reshape(K, nϕ, nθ)

	header["PHSLMN"] = minimum(ϕ)
	set_comment!(header,"PHSLMN","Minimum azimuth in slice")

	header["PHSLMX"] = maximum(ϕ)
	set_comment!(header,"PHSLMX","Maximum azimuth in slice")

	header["NPHI"] = nϕ
	set_comment!(header,"NPHI","Number of grid points in azimuth")

	header["THSLMN"] = minimum(θ)
	set_comment!(header,"THSLMN","Minimum colatitude in slice")

	header["THSLMX"] = maximum(θ)
	set_comment!(header,"THSLMX","Maximum colatitude in slice")

	header["NTH"] = nθ
	set_comment!(header,"NTH","Number of grid points in colatitude")

	header["ROBS"] = r_obs
	set_comment!(header,"ROBS","Observation radius")

	FITS(joinpath(SCRATCH_kerneldir[],"Kϕθ_from_Kst.fits"),"w") do f
		FITSIO.write(f, Kϕθ, header=header)
	end
	Kϕθ
end

#######################################################################
# The following function compute the 3D kernel using the fleshed-out formulation
# and not from the spherical harmonic coefficients
# They are only defined for radial displacements, and are used to compare
# our results with those obtained using a formalism similar to Mandal et al. (2017)
#######################################################################

function populatekernel3D!(::SoundSpeed, K,Ω_ind,
	realconjhωHjₒjₛω_r₁r₂,     Pjₒ1_Pjₛ2,
	realconjhωconjHjₒjₛω_r₂r₁, Pjₒ2_Pjₛ1,
	pre)

	P12 = Pjₒ1_Pjₛ2[Ω_ind]
	P21 = Pjₒ2_Pjₛ1[Ω_ind]

	iszero(P12) && iszero(P21) && return

	for r_ind in UnitRange(axes(K, 1))
		K[r_ind, Ω_ind] += pre * (
						realconjhωHjₒjₛω_r₁r₂[r_ind] * P12  +
						realconjhωconjHjₒjₛω_r₂r₁[r_ind] * P21 )
	end
end

function Kc_partial(localtimer, modes_iter_proc::ProductSplit,
	xobs1, xobs2, los::los_radial, θ_arr, ϕ_arr, hω_arr,
	p_Gsrc::Union{Nothing, ParamsGfn}=nothing,
	p_Gobs1::Union{Nothing, ParamsGfn}=nothing,
	p_Gobs2::Union{Nothing, ParamsGfn}=p_Gobs1,
	r_src=r_src_default, r_obs=nothing)

	p_Gsrc = read_all_parameters(p_Gsrc, r_src=r_src)
	p_Gobs1 = read_all_parameters(p_Gobs1, r_src=xobs1.r)
	p_Gobs2 = read_all_parameters(p_Gobs2, r_src=xobs2.r)

	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, dω = p_Gsrc

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

	nθ, nϕ = length(θ_arr), length(ϕ_arr)
	nΩ = nθ * nϕ
	θϕ_iter = Iterators.product(θ_arr,ϕ_arr)

	# jₛ is the angular degree of the wave from the source, and
	# jₒ is the angular degree of the wave from the observation point
	jₛ_min_proc, jₛ_max_proc = extrema(modes_iter_proc, 1)
	νind_min_proc,νind_max_proc = extrema(modes_iter_proc, 2)
	j_min, j_max = extrema(ℓ_arr)
	jₛ_range_proc = jₛ_min_proc:jₛ_max_proc
	jₒjₛ_allmodes = L₂L₁Δ(jₛ_range_proc, 2j_max, j_min:j_max)
	jₒ_range_proc = l₂_range(jₒjₛ_allmodes)

	# Need to load all the modes on this processor
	modes_obsGfn_proc = ProductSplit((jₒ_range_proc,νind_min_proc:νind_max_proc), 1, 1)

	Gfn_fits_files_src = Gfn_fits_files(Gfn_path_src,
						(ℓ_arr, 1:Nν_Gfn), modes_iter_proc, NGfn_files_src)

	# Gℓ′ω(r, robs) files
	Gfn_fits_files_obs1 = Gfn_fits_files(Gfn_path_obs1,
						(ℓ_arr, 1:Nν_Gfn), modes_obsGfn_proc, NGfn_files_obs1)

	Gfn_fits_files_obs2 = Gfn_fits_files(Gfn_path_obs2,
						(ℓ_arr, 1:Nν_Gfn), modes_obsGfn_proc, NGfn_files_obs2)

	K = zeros(nr, nΩ)

	arrs = allocatearrays(SoundSpeed(), los, obs_at_same_height)
	@unpack Gsrc, drGsrc, Gobs1, drGobs1, Gobs2, drGobs2, divGsrc, divGobs = arrs
	@unpack fjₒjₛ_r₁_rsrc, fjₒjₛ_r₂_rsrc, Hjₒjₛω_r₁r₂, Hjₒjₛω_r₂r₁ = arrs
	@unpack Gobs1_cache, drGobs1_cache, Gobs2_cache, drGobs2_cache = arrs
	@unpack tworealconjhωHjₒjₛω_r₁r₂, tworealconjhωconjHjₒjₛω_r₂r₁ = arrs

	@unpack Gobs1_cache, Gobs2_cache = arrs
	for jₒ in jₒ_range_proc
		Gobs1_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
		drGobs1_cache[jₒ] = zeros(ComplexF64, axes(drGsrc))
		if !obs_at_same_height
			Gobs2_cache[jₒ] = zeros(ComplexF64, axes(Gsrc))
			drGobs2_cache[jₒ] = zeros(ComplexF64, axes(drGsrc))
		end
	end

	@timeit localtimer "legendre" begin
	P1, P2 = zeros(0:j_max, nΩ), zeros(0:j_max, nΩ)
	for (Ω_ind,(θ,ϕ)) in enumerate(θϕ_iter)
		collectPl!(view(P1, :, Ω_ind), cosχ((θ,ϕ), xobs1), lmax = j_max)
		collectPl!(view(P2, :, Ω_ind), cosχ((θ,ϕ), xobs2), lmax = j_max)
	end

	P1 = permutedims(P1)
	P2 = permutedims(P2)
	end # timer

	Pjₒ1_Pjₛ2 = zeros(nΩ)
	Pjₒ2_Pjₛ1 = zeros(nΩ)

	ω_ind_prev = -1
	for (jₛ,ω_ind) in modes_iter_proc

		ω = ω_arr[ω_ind]
		conjhω = conj(hω_arr[ω_ind])
		dωω²Pω = dω/2π * ω^2 * Powspec(ω)

		@timeit localtimer "FITS" begin

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

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

		end # timer

		Grr_r₁_rsrc = Gsrc[r₁_ind, 0]
		Grr_r₂_rsrc = Gsrc[r₂_ind, 0]

		jₒ_range_jₛ = l₂_range(jₒjₛ_allmodes, jₛ)

		for jₒ in jₒ_range_jₛ

			pre = dωω²Pω * (2jₒ+1)/4π * (2jₛ+1)/4π

			Gobs1 = Gobs1_cache[jₒ]
			drGobs1 = drGobs1_cache[jₒ]
			Gobs2 = Gobs2_cache[jₒ]
			drGobs2 = drGobs2_cache[jₒ]

			@timeit localtimer "FITS" begin

			# Green function about receiver location
			if ω_ind != ω_ind_prev || jₒ ∉ l₂_range(jₒjₛ_allmodes, jₛ-1)
				read_Gfn_file_at_index!(Gobs1, Gfn_fits_files_obs1,
					(ℓ_arr, 1:Nν_Gfn),(jₒ,ω_ind), NGfn_files_obs1,:, 1:2, 1, 1)

				# Derivative of Green function about receiver location
				read_Gfn_file_at_index!(drGobs1, Gfn_fits_files_obs1,
					(ℓ_arr, 1:Nν_Gfn),(jₒ,ω_ind), NGfn_files_obs1,:, 1:1, 1, 2)
			end
			end # timer

			@timeit localtimer "radial term" begin
			# precompute the radial term in f
			radial_fn_δc_firstborn!(fjₒjₛ_r₁_rsrc, Gsrc, drGsrc, jₛ, divGsrc,
				Gobs1, drGobs1, jₒ, divGobs)
			end # timer

			@timeit localtimer "radial term 2" begin
			Hjₒjₛω!(Hjₒjₛω_r₁r₂, fjₒjₛ_r₁_rsrc, Grr_r₂_rsrc)
			end
			@timeit localtimer "radial term 3" begin
			@. tworealconjhωHjₒjₛω_r₁r₂ = 2real(conjhω * Hjₒjₛω_r₁r₂)
			end

			if r₁_ind != r₂_ind
				@timeit localtimer "FITS" begin
				if ω_ind != ω_ind_prev || jₒ ∉ l₂_range(jₒjₛ_allmodes, jₛ-1)
					# Green function about receiver location
					read_Gfn_file_at_index!(Gobs2, Gfn_fits_files_obs2,
						(ℓ_arr, 1:Nν_Gfn),(jₒ,ω_ind), NGfn_files_obs2,:, 1:2, 1, 1)

					# Derivative of Green function about receiver location
					read_Gfn_file_at_index!(drGobs2, Gfn_fits_files_obs2,
						(ℓ_arr, 1:Nν_Gfn),(jₒ,ω_ind), NGfn_files_obs2,:, 1:1, 1, 2)
				end
				end # timer

				@timeit localtimer "radial term" begin
				radial_fn_δc_firstborn!(fjₒjₛ_r₂_rsrc, Gsrc, drGsrc, jₛ, divGsrc,
					Gobs2, drGobs2, jₒ, divGobs)

				end # timer
			end
			@timeit localtimer "radial term 2" begin
			Hjₒjₛω!(Hjₒjₛω_r₂r₁, fjₒjₛ_r₂_rsrc, Grr_r₁_rsrc)
			end
			@timeit localtimer "radial term 3" begin
			@. tworealconjhωconjHjₒjₛω_r₂r₁ = 2real(conjhω * conj(Hjₒjₛω_r₂r₁))
			end

			@timeit localtimer "legendre" begin
			@inbounds for Ωind in eachindex(Pjₒ1_Pjₛ2)
				Pjₒ1_Pjₛ2[Ωind] = P1[Ωind, jₒ]*P2[Ωind, jₛ]
				Pjₒ2_Pjₛ1[Ωind] = P1[Ωind, jₛ]*P2[Ωind, jₒ]
			end
			end # timer

			@timeit localtimer "kernel" begin
			for Ω_ind in axes(K, 2)
				populatekernel3D!(SoundSpeed(), K, Ω_ind,
				tworealconjhωHjₒjₛω_r₁r₂, Pjₒ1_Pjₛ2,
				tworealconjhωconjHjₒjₛω_r₂r₁, Pjₒ2_Pjₛ1,
				pre)
			end
			end # timer
		end

		ω_ind_prev = ω_ind

	end

	map(closeGfnfits,(Gfn_fits_files_src, Gfn_fits_files_obs1, Gfn_fits_files_obs2))

	signaltomaster!(timers_channel, localtimer)
	map(finalize_except_wherewhence,(progress_channel, timers_channel))

	return K
end

function generatefitsheader()
    FITSHeader(["r1","th1","phi1","r2","th2","phi2",
		"l_max","m_max","jmin","jmax","nuimin","nuimax",
		"PHSLMN","PHSLMX","NPHI",
		"THSLMN","THSLMX","NTH"],
		Any[float(xobs1.r), float(xobs1.θ), float(xobs1.ϕ),
		float(xobs2.r), float(xobs2.θ), float(xobs2.ϕ), Int(l_max), Int(m_max),
		minimum(ℓ_range), maximum(ℓ_range), minimum(ν_ind_range), maximum(ν_ind_range),
		minimum(ϕ), maximum(ϕ), length(ϕ),
		minimum(θ), maximum(θ), length(θ)],
		["Radius of the first observation point",
		"Colatitude of the first observation point",
		"Azimuth of the first observation point",
		"Radius of the second observation point",
		"Colatitude of the second observation point",
		"Azimuth of the second observation point",
		"Maximum angular degree of the perturbation",
		"Maximum azimuthal order of the perturbation",
		"Minimum wave mode","Maximum wave mode",
		"Minimum wave frequency index","Maximum wave frequency index",
		"Minimum azimuth in slice","Maximum azimuth in slice",
		"Number of grid points in azimuth",
		"Minimum colatitude in slice","Maximum colatitude in slice",
		"Number of grid points in colatitude"])
end

function _Kc(comm, xobs1::Point3D, xobs2::Point3D, los::los_direction;θ,ϕ, kwargs...)
	r_src, r_obs = read_rsrc_robs_c_scale(kwargs)
	p_Gsrc, p_Gobs1, p_Gobs2 = read_parameters_for_points(xobs1, xobs2; kwargs...)

	@unpack ν_arr,ℓ_arr = p_Gsrc
	ℓ_range, ν_ind_range, modes_iter = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)

	hω_arr = get(kwargs,:hω) do
		h = hω(comm, TravelTimes(), xobs1, xobs2, los; kwargs...,
			ℓ_range=ℓ_range, ν_ind_range=ν_ind_range,
			print_timings=false)
		_broadcast(h, 0, comm)
	end

	pmapsum_timed(comm, Kc_partial, modes_iter,
		xobs1, xobs2, los,θ,ϕ, hω_arr,
		p_Gsrc, p_Gobs1, p_Gobs2, r_src, r_obs)
end

function Kc_latitudinal_slice(comm, xobs1::Point3D, xobs2::Point3D, los::los_radial=los_radial(); kwargs...)

	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack ℓ_arr,ν_arr = p_Gsrc
	ℓ_range, = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	_,ℓmax = extrema(ℓ_range)

	ϕ = get(kwargs,:ϕ, LinRange(0, 2π, get(kwargs,:nϕ, 4ℓmax)))
	θ = get(kwargs,:θ,(xobs1.θ+xobs2.θ)/2)

	Kc_3D_lat_slice = _Kc(comm, xobs1, xobs2, los; kwargs..., θ=θ, ϕ=ϕ)
	Kc_3D_lat_slice === nothing && return nothing

	header = generatefitsheader()

	filename = joinpath(SCRATCH_kerneldir[],
		"Kc_3D_lat_slice_jcutoff$(ℓmax).fits")

	FITS(filename,"w") do f
		FITSIO.write(f, Kc_3D_lat_slice, header=header)
	end
	Kc_3D_lat_slice
end

function Kc_longitudinal_slice(comm, xobs1::Point3D, xobs2::Point3D, los::los_radial=los_radial(); kwargs...)
	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack ℓ_arr,ν_arr = p_Gsrc
	ℓ_range, = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	_,ℓmax = extrema(ℓ_range)

	θ = get(kwargs,:θ, LinRange(0,π, get(kwargs,:nθ, 4ℓmax)))
	ϕ = get(kwargs,:ϕ,(xobs1.ϕ+xobs2.ϕ)/2)

	Kc_3D_long_slice = _Kc(comm, xobs1, xobs2, los; kwargs..., θ=θ, ϕ=ϕ)
	Kc_3D_long_slice === nothing && return nothing

	header = generatefitsheader()

	filename = joinpath(SCRATCH_kerneldir[],
		"Kc_3D_long_slice_jcutoff$(ℓmax).fits")

	FITS(filename,"w") do f
		FITSIO.write(f, Kc_3D_long_slice, header=header)
	end
	Kc_3D_long_slice
end

# 3D profile
function Kc(comm, xobs1::Point3D, xobs2::Point3D, los::los_radial=los_radial(); kwargs...)

	p_Gsrc = read_all_parameters(; kwargs...)
	@unpack ℓ_arr,ν_arr = p_Gsrc
	ℓ_range,_ = ℓ_and_ν_range(ℓ_arr,ν_arr; kwargs...)
	_,ℓmax = extrema(ℓ_range)
	# These may be derived form ℓmax

	l_max = get(kwargs,:l_max, 2ℓmax);

	Nnodes = get(kwargs,:nθ, ceil(Int,(l_max+1)/2))
	nodes,_ = gausslegendre(Nnodes)
	θdefault = acos.(nodes)[end:-1:1] # flipped to get to increasing θ from increasing cosθ

	θ = get(kwargs,:θ,θdefault)
	ϕ = get(kwargs,:ϕ, LinRange(0, 2π, get(kwargs,:nϕ, 2l_max+1)))

	K = _Kc(comm, xobs1, xobs2, los; kwargs...,θ=θ,ϕ=ϕ)
	K === nothing && return nothing
	Kc_3D = copy(reshape(K, nr, length(θ), length(ϕ)))

	header = generatefitsheader()

	filename = joinpath(SCRATCH_kerneldir[],
		"Kc_3D_jcutoff$(ℓmax).fits")

	FITS(filename,"w") do f
		FITSIO.write(f, Kc_3D, header=header)
	end
	Kc_3D
end


text:
"""
    Vertex

Struct containing bare values of vertex component and the respective s, t and u channel.
* `bare :: Vector{Float64}` : bare couplings on irreducible lattice sites
* `ch_s :: Channel`         : particle-particle (s) channel
* `ch_t :: Channel`         : direct particle-hole (t) channel
* `ch_u :: Channel`         : crossed particle-hole (u) channel
"""
struct Vertex
    bare :: Vector{Float64}
    ch_s :: Channel
    ch_t :: Channel
    ch_u :: Channel
end

# generate vertex dummy
function get_vertex_empty(
    r :: Reduced_lattice,
    m :: Mesh
    ) :: Vertex

    # init bare
    bare = zeros(Float64, length(r.sites))

    # init channels
    ch_s = get_channel_empty(r, m)
    ch_t = get_channel_empty(r, m)
    ch_u = get_channel_empty(r, m)

    # build vertex
    Γ = Vertex(bare, ch_s, ch_t, ch_u)

    return Γ
end

# get interpolated value of vertex in certain channel for a given frequency buffer
function get_vertex(
    site :: Int64,
    b    :: Buffer,
    Γ    :: Vertex,
    ch   :: Int64
    )    :: Float64

    val = 0.0

    if ch == 1
        val = get_channel(site, b, Γ.ch_s)
    elseif ch == 2
        val = get_channel(site, b, Γ.ch_t)
    elseif ch == 3
        val = get_channel(site, b, Γ.ch_u)
    end

    return val
end

# get interpolated value of vertex in certain channel for a given frequency buffer on all lattice sites
function get_vertex_avx!(
    r        :: Reduced_lattice,
    b        :: Buffer,
    Γ        :: Vertex,
    ch       :: Int64,
    temp     :: SubArray{Float64, 1, Array{Float64, 3}},
    exchange :: Bool,
    sgn      :: Float64
    )        :: Nothing

    if ch == 1
        get_channel_avx!(r, b, Γ.ch_s, temp, exchange, sgn)
    elseif ch == 2
        get_channel_avx!(r, b, Γ.ch_t, temp, exchange, sgn)
    elseif ch == 3
        get_channel_avx!(r, b, Γ.ch_u, temp, exchange, sgn)
    end

    return nothing
end

# replace vertex with another vertex (except for bare)
function replace_with!(
    Γ1 :: Vertex,
    Γ2 :: Vertex
    )  :: Nothing

    replace_with!(Γ1.ch_s, Γ2.ch_s)
    replace_with!(Γ1.ch_t, Γ2.ch_t)
    replace_with!(Γ1.ch_u, Γ2.ch_u)

    return nothing
end

# multiply vertex with factor (except for bare)
function mult_with!(
    Γ   :: Vertex,
    fac :: Float64
    )   :: Nothing

    mult_with!(Γ.ch_s, fac)
    mult_with!(Γ.ch_t, fac)
    mult_with!(Γ.ch_u, fac)

    return nothing
end

# multiply vertex with some factor and add to other vertex (except for bare)
function mult_with_add_to!(
    Γ2  :: Vertex,
    fac :: Float64,
    Γ1  :: Vertex
    )   :: Nothing

    mult_with_add_to!(Γ2.ch_s, fac, Γ1.ch_s)
    mult_with_add_to!(Γ2.ch_t, fac, Γ1.ch_t)
    mult_with_add_to!(Γ2.ch_u, fac, Γ1.ch_u)

    return nothing
end

"""
    get_abs_max(
        Γ :: Vertex
        ) :: Float64

Returns maximum absolute value of a vertex component across channels.
"""
function get_abs_max(
    Γ :: Vertex
    ) :: Float64

    max_s = get_abs_max(Γ.ch_s)
    max_t = get_abs_max(Γ.ch_t)
    max_u = get_abs_max(Γ.ch_u)
    Γ_max = max(max_s, max_t, max_u)

    return Γ_max
end

# resample a vertex component to new meshes via trilinear interpolation
function resample_from_to!(
    m_old :: Mesh,
    Γ_old :: Vertex,
    m_new :: Mesh,
    Γ_new :: Vertex
    )     :: Nothing

    resample_from_to!(m_old.Ωs, m_old.νs, Γ_old.ch_s, m_new.Ωs, m_new.νs, Γ_new.ch_s)
    resample_from_to!(m_old.Ωt, m_old.νt, Γ_old.ch_t, m_new.Ωt, m_new.νt, Γ_new.ch_t)
    resample_from_to!(m_old.Ωs, m_old.νs, Γ_old.ch_u, m_new.Ωs, m_new.νs, Γ_new.ch_u)

    return nothing
end

ext: jl


text:
using StaticArrays

push_str!(ex::Vector{Expr}, s::String) = push!(ex, Meta.parse(s))
append_str!(ex::Vector{Expr}, s::Vector{String}) = append!(ex, Meta.parse.(s))

# univariate monomial
@inline _m1(α::Number, x::T) where {T <: Number} =
      (@fastmath x^α)

# univariate monomial derivative
@inline _m1d(α::Number, x::T) where {T <: Number} =
      α == 0 ? zero(T) : (@fastmath α * x^(α-1))

function monomial(α, x::SVector{K, T}) where {K, T}
   m = one(T)
   for i = 1:K
      @inbounds m *= _m1(α[i], x[i])
   end
   return m
end

@generated function monomial_d(α, x::SVector{K, T}) where {K, T}
   exprs = Expr[]
   # @assert length(α) >= K
   # evaluate the scalar monomials
   #  f = SVector{...}( x[1]^α[1], ...)
   # df = SVector{...}( α[1] * x[1]^(α[1]-1), ...)
   ex_f = "f = @SVector $T["
   ex_df = "df = @SVector $T["
   for i = 1:K
      ex_f  *= " _m1(α[$i], x[$i]), "
      ex_df *= "_m1d(α[$i], x[$i]), "
   end
   ex_f  =  ex_f[1:end-2] * "]"
   ex_df = ex_df[1:end-2] * "]"

   push_str!(exprs, ex_f)
   push_str!(exprs, ex_df)
   push_str!(exprs, "m = 1.0")
   for j = 1:K
      push_str!(exprs, "m *= f[$j]")
   end

   # evaluate the derivative
   for i = 1:K
      push_str!(exprs, "dm_$i = 1.0")
      for j = 1:K
         if i == j
            push_str!(exprs, "dm_$i *= df[$j]")
         else
            push_str!(exprs, "dm_$i *= f[$j]")
         end
      end
   end
   ex_dm = "m_d = @SVector $T["
   for i = 1:K
      ex_dm *= "dm_$i, "
   end
   ex_dm = ex_dm[1:end-2] * "]"
   push_str!(exprs, ex_dm)

   quote
      $(Expr(:meta, :inline))
      @inbounds $(Expr(:block, exprs...))
      return m, m_d
   end
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
@inline function new_graph(::Type{I}, ::Type{F}) where I where F
    g = MetaGraphs.MetaDiGraph{I, F}()
    MetaGraphs.set_indexing_prop!(g, :class)
    return PharmGraph{I, F}(g)
end

@inline function build_graph(; kwargs...)
    config = Config(; kwargs...)
    return build_graph(config)
end

@inline function build_graph(config::Config)
    graph = new_graph(Int64, Float64)
    return build_graph!(graph::PharmGraph, config::Config)
end

@inline function build_graph!(graph::PharmGraph, config::Config)
    all_filenames = [
        config.rxnrel,
        config.rxnsat,
    ]
    all_numlines = countlines_filename.(all_filenames)
    total_numlines = sum(all_numlines)
    showprogress::Bool = config.showprogress
    wait_time = showprogress ? Float64(1.0) : Float64(Inf)
    p = ProgressMeter.Progress(total_numlines + 10, wait_time)

    open(config.rxnsat, "r") do io
        for line in eachline(io)
            ProgressMeter.next!(p;
                showvalues = [
                    (:Stage, "1 of 3"),
                    (:File, "RXNSAT.RRF"),
                ],
            )
            elements = split(line, "|")
            left_system = "RXCUI"
            left_value = elements[1]
            RIGHT_SYSTEM_TYPE = elements[9]
            if RIGHT_SYSTEM_TYPE == "NDC"
                add_right_node = true
                edges = :bidirectional
                right_system = "NDC"
                right_value = elements[11]
            elseif RIGHT_SYSTEM_TYPE == "ATC_LEVEL"
                add_right_node = true
                edges = :right_to_left
                right_system = string(elements[10], elements[11])
                right_value = elements[6]
            else
                add_right_node = false
            end
            left_node = PharmClass(left_system, left_value)
            PharmaceuticalClassification._add_edge!(graph, left_node => left_node)
            if add_right_node
                right_node = PharmClass(right_system, right_value)
                PharmaceuticalClassification._add_edge!(graph, right_node => right_node)
                if edges == :left_to_right || edges == :bidirectional
                    PharmaceuticalClassification._add_edge!(graph, left_node => right_node)
                end
                if edges == :right_to_left || edges == :bidirectional
                    PharmaceuticalClassification._add_edge!(graph, right_node => left_node)
                end
            end
        end
    end

    ProgressMeter.next!(p;
        showvalues = [
            (:Stage, "2 of 3"),
            (:Description, "Add edges for ATC hierarchy"),
        ],
    )
    add_atc_hierarchy_edges!(graph)

    open(config.rxnrel, "r") do io
        for line in eachline(io)
            ProgressMeter.next!(p;
                showvalues = [
                    (:Stage, "3 of 3"),
                    (:File, "RXNREL.RRF"),
                ],
            )
            elements = split(line, "|")
            if elements[11] == "RXNORM" && elements[3] == "CUI" && elements[7] == "CUI"
                relationship = elements[8]
                if relationship == "consists_of" || relationship == "has_ingredient"
                    more_specific_rxcui = elements[5]
                    element_1 = elements[1]
                    if occursin(":", element_1)
                        less_specific_rxcui = split(element_1, ':')[2]
                    else
                        less_specific_rxcui = element_1
                    end
                    more_specific_node = PharmClass("RXCUI", more_specific_rxcui)
                    less_specific_node = PharmClass("RXCUI", less_specific_rxcui)
                    PharmaceuticalClassification._add_edge!(graph, less_specific_node => more_specific_node)
                end
            end
        end
    end

    ProgressMeter.finish!(p)

    return graph
end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.
##################################################################################

import ..SimData
import ..rawdata

"""
    seriesoverlay(ts1, ts2)

Return a new [`TSeries`](@ref) over the full range of both arguments. The overlapping part
contains values from the last argument.

!!! note "Deprecation Note"

    This function will be removed in the future. Use
    `TimeSeriesEcon.overlay` instead.  Note the important difference that
    in `TimeSeriesEcon.overlay` the values in the overlay come from the *first*
    series where the value exists, as opposed to `seriesoverlay` where it was
    from the last one.

See also: [`dictoverlay`](@ref)

"""
function seriesoverlay end
export seriesoverlay
@deprecate seriesoverlay(ts1::TSeries, ts2::TSeries) overlay(ts2, ts1)

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

"""
    dictoverlay(d1, d2)

Merge two dictionaries. Common key where the values are [`TSeries`](@ref) of the
same frequency are overlayed. Otherwise, a common key takes the value of the
last Dict containing it.

!!! note "Deprecation Note"

    This function will be removed. Use `TimeSeriesEcon.overlay` instead. An
    important difference is that `TimeSeriesEcon.overlay` keeps the values from
    the first argument where the key appears, while `dictoverlay` keeps it from
    the last.

"""
function dictoverlay end
export dictoverlay
@deprecate dictoverlay(D1::Dict{String,<:Any}, D2::Dict{String,<:Any}) overlay(Workspace(D2), Workspace(D1))

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

"""
    array2data(matrix, vars, range)
    array2data(matrix, model, plan)
    array2workspace(matrix, vars, range)
    array2workspace(matrix, model, plan)
    data2array(simdata)
    data2workspace(simdata)
    workspace2array(w, vars, range)
    workspace2array(w, model, plan)
    workspace2data(w, vars, range)
    workspace2data(w, model, plan)
"""
(
    function array2data end,
    function array2workspace end,
    function data2array end,
    function data2workspace end,
    function workspace2array end,
    function workspace2data end
)
export array2data, array2workspace, data2array, data2workspace, workspace2array, workspace2data

@inline array2data(matrix::AbstractMatrix, model::Model, plan::Plan; copy = false) = array2data(matrix, model.varshks, plan.range; copy = copy)
@inline array2data(matrix::AbstractMatrix, vars, range; copy = false) = SimData(range, vars, copy ? Base.copy(matrix) : matrix)

@inline array2workspace(matrix::AbstractMatrix, model::Model, plan::Plan; copy = false) = array2workspace(matrix, model.varshks, plan.range; copy = copy)
@inline array2workspace(matrix::AbstractMatrix, vars, range; copy = false) = Workspace(Symbol(v) => TSeries(range, copy ? Base.copy(matrix[:, i]) : matrix[:, i]) for (i, v) = enumerate(vars))

@inline data2array(simdata::SimData, m::Model, p::Plan; copy = false) = data2array(simdata, m.varshks, p.range; copy = copy)
@inline data2array(simdata::SimData; copy = false) = copy ? Base.copy(rawdata(simdata)) : rawdata(simdata)
@inline data2array(simdata::SimData, vars, range; copy = false) = copy ? Base.copy(rawdata(simdata[range, vars])) : rawdata(simdata[range, vars])
@inline data2workspace(simdata::SimData, m::Model, p::Plan; copy = false) = data2workspace(simdata, m.varshks, p.range; copy = copy)
@inline data2workspace(simdata::SimData; copy = false) = data2workspace(simdata, axes(simdata,2), rangeof(simdata); copy = copy)
@inline data2workspace(simdata::SimData, vars, range; copy = false) = Workspace(Symbol(v) => copy ? Base.copy(simdata[range, v]) : simdata[range, v] for v in vars)

@inline workspace2array(w::Workspace, vars, range::AbstractUnitRange; copy = false) = hcat((w[Symbol(v)][range] for v in vars)...)
@inline workspace2array(w::Workspace, model::Model, plan::Plan; copy = false) = workspace2array(w, model.varshks, plan.range; copy = copy)

function workspace2data(w::Workspace, vars, range::AbstractUnitRange; copy = false)
    ret = SimData(range, vars)
    for v in vars
        copyto!(ret[v], range, w[Symbol(v)])
    end
    return ret
end
@inline workspace2data(w::Workspace, model::Model, plan::Plan; copy = false) = workspace2data(w, model.varshks, plan.range; copy = copy)

@inline workspace2data(w::Workspace, model::Model; copy = false) = workspace2data(w, model.varshks; copy = copy)
@inline workspace2array(w::Workspace, model::Model; copy = false) = workspace2array(w, model.varshks; copy = copy)
function workspace2array(w::Workspace, vars; copy = false)
    range = mapreduce(v -> rangeof(w[Symbol(v)]), intersect, vars)
    return hcat((w[Symbol(v)][range] for v in vars)...)
end

function workspace2data(w::Workspace, vars; copy = false)
    range = mapreduce(v -> rangeof(w[Symbol(v)]), intersect, vars)
    return hcat(MVTSeries(range); (v => w[Symbol(v)] for v in vars)...)
end

export dict2array, array2dict, dict2data, data2dict
@deprecate dict2array(d::AbstractDict, args...; kwargs...) workspace2array(Workspace(d), args...; kwargs...)
@deprecate dict2data(d::AbstractDict, args...; kwargs...) workspace2data(Workspace(d), args...; kwargs...)
@deprecate dict2array(d::Workspace, args...; kwargs...) workspace2array(d, args...; kwargs...)
@deprecate dict2data(d::Workspace, args...; kwargs...) workspace2data(d, args...; kwargs...)
@deprecate array2dict(args...; kwargs...) array2workspace(args...; kwargs...)
@deprecate data2dict(args...; kwargs...) data2workspace(args...; kwargs...)


text:
import Base: getindex

@entitytype Entity
@entitytype Project
@entitytype Folder
@entitytype File
@entitytype Link
@entitytype Schema

@dicttype DictObject  pydict_object.DictObject
@dicttype Evaluation
@dicttype Submission
@dicttype Activity
@dicttype Annotations pyannotations.Annotations
@dicttype Wiki
@dicttype Team

create(::Type{Entity},args...;kwargs...) = synapsecall(synapseclient.Entity,:create,args...;kwargs...)

@synapsefunction AbstractEntity local_state

@synapsefunction Activity used
@synapsefunction Activity executed


text:
# # [Basic SIR example](@id sir_basic)
# Sean L. Wu (@slwu89), 2021-1-8

using Individual.Sampling
using Individual.SchemaBase

using Catlab.Present, Catlab.CSetDataStructures, Catlab.Theories, Catlab.CategoricalAlgebra, Catlab.Graphics, Catlab.Graphs
using Plots, GraphViz

# ## Introduction
# The SIR (Susceptible-Infected-Recovered) model is the "hello, world!" model of for infectious disease simulations,
# and here we describe how to use the basic schema for Markov models to build it in Individual.jl. We only use the
# base `TheoryIBM` schema, because the model is a Markov chain. This tutorial largely mirrors the SIR model tutorial from
# the [R package "individual"](https://mrc-ide.github.io/individual/articles/Tutorial.html), which inspired individual.jl.

# For more information about the SIR model and stochastic epidemic models in general, please see ["A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis"](https://doi.org/10.1016/j.idm.2017.03.001).

# The schema looks like this:

to_graphviz(TheoryIBM)

# There are morphisms from the object `Person` into `State`. The first, `state` gives the current state of each individual in the
# simulation. The second `state_update` is where state updates can be queued; at the end of a time step, `state` is swapped
# with `state_update` and `state_update` is reset. This ensures that state updates obey a FIFO order, and that individuals
# update synchronously.

# There are also labels for states to make it easier to write self-documenting models.

# The basic schema could be extended with further `Attr`s if needing to model individual level heterogeneity, like immmune respose, etc.

# ## Parameters

# To start, we should define some parameters. The epidemic will be simulated in a population of 1000, where 8 persons are initially infectious, whose indices are randomly sampled.
# The effective contact rate β will be a function of the deterministic R0 and recovery rate γ. We also specify Δt, which is the size of the time step.
# Because individual’s time steps are all of unit length, we scale transition probabilities by Δt to create models with different sized steps, interpreting the discrete time model as a discretization of a continuous time model.
# If the maximum time is tmax then the overall number of time steps is tmax/Δt.

N = 1000
I0 = 8
S0 = N - I0
Δt = 0.1
tmax = 100
steps = Int(tmax/Δt)
γ = 1/10 # recovery rate
R0 = 2.5
β = R0 * γ # R0 for corresponding ODEs

initial_states = fill("S", N)
initial_states[rand(1:N, I0)] .= "I"
state_labels = ["S", "I", "R"];

# ## Model object

# The "IBM" (Individual Based Model) schema needs the type parameter `String` because it defines a
# single attribute, that giving names to the categorical set of states.

const SIR = IBM{String}()
initialize_states(SIR, initial_states, state_labels);

# ## Processes

# In order to model infection, we need a process. This is a function that takes only a single argument, `t`, for the current time step (unused here, but can model time-dependent processes, such as seasonality or school holiday).
# Within the function, we get the current number of infectious individuals, then calculate the per-capita force of infection on each susceptible person, λ=βI/N.
# Next we get the set of susceptible individuals and use the sample method to randomly select those who will be infected on this time step.
# The probability is given by $$ 1- e^{-\lambda \Delta t} $$. The method `bernoulli_sample` will automatically calculate that probability when given 3 arguments. Finally, we queue a state update for those individuals who were sampled.

function infection_process(t::Int)
    I = npeople(SIR, "I")
    N = npeople(SIR)
    λ = β * I/N
    S = get_index_state(SIR, "S")
    S = bernoulli_sample(S, λ, Δt)
    queue_state_update(SIR, S, "I")
end

# The recovery process is simpler, as the per-capita transition probability from I to R does not depend on the state of the system.

function recovery_process(t::Int)
    I = get_index_state(SIR, "I")
    I = bernoulli_sample(I, γ, Δt)
    queue_state_update(SIR, I, "R")
end

# ## Simulation

# We use `render_process` to create a rendering (output) process and
# a matrix giving state counts by time step.

# Then we draw a trajectory using `simulation_loop` and plot the results.
# Please note that `simulation_loop` is a simple wrapper around the correct updating order
# and that a user can always write a simple for loop to update the model.

state_out, render_process = render_states(SIR, steps)

simulation_loop(SIR, [infection_process, recovery_process, render_process], steps)

plot(
    (1:steps) * Δt,
    state_out,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)


text:
using VoronoiCells
using GeometryBasics
using Test

@testset "Plotting" begin
    @testset "Edges for plotting" begin
        points = [Point2(0.25, 0.25), Point2(0.75, 0.25), Point2(0.5, 0.75)]
        rect = Rectangle(Point2(0, 0), Point2(1, 1))
        tess = voronoicells(points, rect)

        p = VoronoiCells.corner_coordinates(tess)

        @test length(filter(isnan, p)) == length(points)

        expected_corners = [
            [0.0, 0.6875],
            [0.0, 0.0],
            [0.5, 0.0],
            [0.5, 0.4375],
            [0.0, 0.6875],
            [0.5, 0.4375],
            [0.5, 0.0],
            [1.0, 0.0],
            [1.0, 0.6875],
            [0.5, 0.4375],
            [0.0, 1.0],
            [0.0, 0.6875],
            [0.5, 0.4375],
            [1.0, 0.6875],
            [1.0, 1.0],
            [0.0, 1.0]
        ]

        @test filter(!isnan, p) ≈ expected_corners atol = 1e-4
    end
end


text:
if VERSION > v"0.6.3"
    using Printf
end

include("alias.jl")
include("utils.jl")
include("PGBF.jl")
include("CGBF.jl")
include("overlap.jl")

#=
using PyPlot
const plt = PyPlot

function do_plot_x( bf::PGBF, filesave::String )
    const N = 100
    x = linspace(-5.0,5.0,N)
    y = zeros(N)
    for i = 1:N
      y[i] = evaluate( bf, x[i], 0.0, 0.0 )
    end
    plt.clf()
    plt.plot( x, y, marker="o" )
    plt.grid()
    plt.savefig(filesave, dpi=300)
end
=#

function test_main()
    bf1 = PGBF(0.3, (0.0,0.0,0.0), (1,0,0))
    println(bf1)
    #do_plot_x(bf1,"bf1.png")

    bf2 = PGBF(0.4,(0.1,0.0,0.0), (1,0,0))
    println(bf2)
    #do_plot_x(bf2,"bf2.png")

    S1 = overlap(bf1,bf1)
    @printf("overlap(bf1,bf1) = %18.10f\n", S1)

    S1 = overlap(bf2,bf2)
    @printf("overlap(bf2,bf2) = %18.10f\n", S1)

    S1 = overlap(bf1,bf2)
    @printf("overlap(bf1,bf2) = %18.10f\n", S1)
end

@time test_main()


text:
#### Tendency types

# Terminology:
#
# `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
# `__Tᵤ__   ____T₁____   ______T₂______   ___S___`
#
#  - `Yᵢ` - the i-th prognostic variable
#  - `Y` - the prognostic state (column vector)
#  - `G = ∇Y` - the gradient of the prognostic state (rank 2 tensor)
#  - `F₁` - the first order tendency flux (rank 2 tensor)
#  - `F₂` - the second order tendency flux (rank 2 tensor)
#  - `Tᵤ` - the explicit time derivative (column vector)
#  - `T₁` - the first order flux divergence (column vector)
#  - `T₂` - the second order flux divergence (column vector)
#  - `S` - the non-conservative source (column vector)

export PrognosticVariable
export FirstOrder, SecondOrder
export Flux, Source
export TendencyDef

"""
    PrognosticVariable

Subtypes are used for specifying
each prognostic variable.
"""
abstract type PrognosticVariable end

"""
    AbstractOrder

Subtypes are used for dispatching
on the flux order.
"""
abstract type AbstractOrder end
struct FirstOrder <: AbstractOrder end
struct SecondOrder <: AbstractOrder end

"""
    AbstractTendencyType

Subtypes are used for specifying a
tuple of tendencies to be accumulated.
"""
abstract type AbstractTendencyType end
struct Flux{O <: AbstractOrder} <: AbstractTendencyType end
struct Source <: AbstractTendencyType end

"""
    TendencyDef

Subtypes are used for specifying
each tendency definition.
"""
abstract type TendencyDef{TT <: AbstractTendencyType, PV <: PrognosticVariable} end

"""
    eq_tends(::PrognosticVariable, ::BalanceLaw, ::AbstractTendencyType)

A tuple of `TendencyDef`s given
 - `PrognosticVariable` prognostic variable
 - `AbstractTendencyType` tendency type
 - `BalanceLaw` balance law

i.e., a tuple of `TendencyDef`s corresponding
to `F₁`, `F₂`, **or** `S` for a single
prognostic variable in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
function eq_tends end

"""
    prognostic_vars(::BalanceLaw)

A tuple of `PrognosticVariable`s given
the `BalanceLaw`.

i.e., a tuple of `PrognosticVariable`s
corresponding to the column-vector `Yᵢ` in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
prognostic_vars(::BalanceLaw) = ()

export sources
"""
    sources(bl::BalanceLaw)

A tuple of `TendencyDef{Source}`s
given the `BalanceLaw`.

i.e., a tuple of `TendencyDef{Source}`s
corresponding to the column-vector `S` in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
function sources(bl::BalanceLaw)
    tend = eq_tends.(prognostic_vars(bl), Ref(bl), Ref(Source()))
    tend = filter(x -> x ≠ nothing, tend)
    return Tuple(Iterators.flatten(tend))
end

export fluxes
"""
    fluxes(bl::BalanceLaw, order::O) where {O <: AbstractOrder}

A tuple of `TendencyDef{Flux{O}}`s
given the `BalanceLaw` and the `order::O`.

i.e., a tuple of `TendencyDef{Flux{O}}`s
corresponding to the column-vector `F₁`
or `F₂` given the flux order `order::O` in:

    `∂_t Yᵢ + (∇•F₁(Y))ᵢ + (∇•F₂(Y,G)))ᵢ = (S(Y,G))ᵢ`
"""
function fluxes(bl::BalanceLaw, order::O) where {O <: AbstractOrder}
    tend = eq_tends.(prognostic_vars(bl), Ref(bl), Ref(Flux{O}()))
    tend = filter(x -> x ≠ nothing, tend)
    return Tuple(Iterators.flatten(tend))
end


text:
#=

  Strimko puzzle in Julia.

  From
  360: A New Twist on Latin Squares
  http://threesixty360.wordpress.com/2009/08/04/a-new-twist-on-latin-squares/
  """
  The idea is simple: each row and column of an nxn grid must contain
  the number 1, 2, ... n exactly once (that is, the grid must form a
  Latin square), and each "stream" (connected path in the grid) must
  also contain the numbers 1, 2, ..., n exactly once.
  """

  For more information, see:
  * http://www.strimko.com/
  * http://www.strimko.com/rules.htm
  * http://www.strimko.com/about.htm
  * http://www.puzzlersparadise.com/Strimko.htm

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function strimko(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    streams = problem[:streams]
    placed = problem[:placed]
    n, _ = size(streams)

    @variable(model, 1 <= x[1:n,1:n] <= n, Int)

    latin_square(model, x)

    # streams
    for s in 1:n
        positions = [x[i,j] for i in 1:n, j in 1:n if s == streams[i,j]]
        @constraint(model,positions in CS.AllDifferentSet())
    end

    # places
    for p in eachrow(placed)
        @constraint(model,x[p[1], p[2]] == p[3])
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                display(x_val)

            end
        end
    else
        println("status:$status")
    end

    return status
end

strimko_problems = Dict(
    #
    # Strimko Monthly #02
    # Via http://www.hakank.org/minizinc/strimko2_002.dzn
    #
    :2  => Dict(
        :streams => resize_matrix([[1,1,2,2,2,2,2],
                   [1,1,2,3,3,3,2],
                   [1,4,1,3,3,5,5],
                   [4,4,3,1,3,5,5],
                   [4,6,6,6,7,7,5],
                   [6,4,6,4,5,5,7],
                   [6,6,4,7,7,7,7]]),
        :placed => resize_matrix([[2,1,1],
                   [2,3,7],
                   [2,5,6],
                   [2,7,4],
                   [3,2,7],
                   [3,6,1],
                   [4,1,4],
                   [4,7,5],
                   [5,2,2],
                   [5,6,6]])
    ),

    #
    # Strimko Weekly Set 067
    # Via http://www.hakank.org/minizinc/strimko2_067.dzn
    :67 => Dict(
       :streams =>  resize_matrix([[1,1,1,2,3],
                 [1,2,2,2,3],
                 [1,2,4,5,3],
                 [5,4,5,4,3],
                 [4,5,5,4,3
                 ]]),
        :placed =>  resize_matrix([[1,3,4],
                [1,4,1],
                [3,3,2],
                [3,5,3],
                [5,4,5]])
    ),

    #
    # Strimko Weekly Set 068
    # Via http://www.hakank.org/minizinc/strimko2_068.dzn
    :68 => Dict(

        :streams => resize_matrix([[1,2,2,4],
               [2,1,4,2],
               [3,4,1,3],
               [4,3,3,1]]),
        :placed =>  resize_matrix([[2,2,3],
               [2,3,2],
               [3,3,1]])
    ),

    # Strimko Weekly Set 069
    # Via http://www.hakank.org/minizinc/strimko2_069.dzn
    :69  => Dict(
        :streams => resize_matrix([[1,2,3,3,3,4],
                  [2,1,3,5,4,3],
                  [2,1,3,5,5,4],
                  [2,6,1,6,5,4],
                  [2,6,1,6,4,5],
                  [6,2,6,1,5,4]]),
        :placed =>  resize_matrix([[2,2,4],
                  [2,3,1],
                  [2,4,3],
                  [2,5,2],
                  [3,2,1],
                  [3,5,6],
                  [4,3,5],
                  [4,4,2]])
    ),

    # Strimko Weekly Set 070
    # Via http://www.hakank.org/minizinc/strimko2_070.dzn
    :70 => Dict(
        :streams =>  resize_matrix([[1,2,3,3,3],
                   [2,1,1,3,1],
                   [2,2,3,1,4],
                   [5,2,5,4,4],
                   [5,5,5,4,4]]),
        :placed =>  resize_matrix([[1,1,1],
                   [2,5,4],
                   [4,1,2],
                   [5,4,5]])
    )

)

for p in sort(collect(keys(strimko_problems)))
    println("\nProblem $p")
    @time strimko(strimko_problems[p])
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
stoichmat(rn::ReactionSystem) = prodstoichmat(rn) - substoichmat(rn)

reformat_jumps(S::Matrix, species_to_index::Dict, x::AbstractVector) =  [x .+ S[i,:] for i in 1:size(S,1)]

split_poly(p::Polynomial, x::AbstractVector, t) = prod(x .^ p.x.Z[1][1:length(x)])*t^p.x.Z[1][end]

Base.show(io::IO, p::MarkovProcess) = println(io, "MarkovProcess")

∂(p,x) = differentiate.(p,x)
∂²(p,x,y) = differentiate(∂(p,x),y)
∂(X) = [intersect(@set(X.p[i] == 0), [@set(X.p[k] >= 0) for k in 1:length(X.p) if k != i]...) for i in 1:length(X.p)]

function reformat_reactions(rxns::Vector{Reaction}, species_to_index::Dict, x::AbstractVector)
    props = []
    for r in rxns
        @unpack rate, substrates, substoich, only_use_rate = r
        a = rate*polynomial(MonomialVector(x,0))
        if !only_use_rate
            for (s, ν) in enumerate(substoich)
                idx = species_to_index[substrates[s]]
                a *= prod(x[idx] - i for i in 0:ν-1)/factorial(ν) # consistent with rxns
            end
        end
        push!(props, a)
    end
    return props
end

function stoich_bounds(rn::ReactionSystem, x0::Dict, solver)
    S = stoichmat(rn)
    B = nullspace(S)
    specs = species(rn)
    x_init = [x0[spec] for spec in specs]
    m = Model(solver)
    @variable(m, x[1:length(specs)] >= 0)
    @constraint(m, B'*x .== B'*x_init)
    scales = Dict()
    for i in 1:length(specs)
        @objective(m, Max, x[i])
        optimize!(m)
        if termination_status(m) == MOI.DUAL_INFEASIBLE
            @warn  "state space unbounded, autoscaling failed!"
            scales[specs[i]] = 1.0
        else
            scales[specs[i]] = min(1.0, objective_value(m))
        end
    end
    return scales
end

extract_primal_solution(gmp::GMPModel) = objective_value(gmp), termination_status(gmp), MOI.get(gmp.approximation_model, MOI.SolveTime())
extract_solution(gmp::GMPModel) = dual_objective_value(gmp), dual_status(gmp), MOI.get(gmp.approximation_model, MOI.SolveTime())
extract_primal_solution(model::Model) = objective_value(model), termination_status(model), MOI.get(model, MOI.SolveTime())
extract_solution(model::Model) = dual_objective_value(model), dual_status(model), MOI.get(model, MOI.SolveTime())

function remap_states!(P::ReactionProcess, bnds::Bounds)
    bnds.bounds = Dict(P.state_to_species[key] => bnds.bounds[key] for key in keys(bnds.bounds))
    bnds.info.solution_time = Dict(P.state_to_species[key] => bnds.info.solution_time[key] for key in keys(bnds.info.solution_time))
    bnds.info.termination_status = Dict(P.state_to_species[key] => bnds.info.termination_status[key] for key in keys(bnds.info.termination_status))
end

function build_approximate_model(model::GMPModel, ::MO.AbstractPrimalMode)
    degree = MO.approximation_info(model).degree
    amodel = Model()
    # initiate moments
    pvar = Dict(i => MO.moments_variable(amodel, v, degree) for (i,v) in MO.gmp_variables(model))
    # add substitutions

    # add measure condition on moments
    just = Dict{Int, Any}()
    scheme_parts = Dict{Int, Vector{MO.SchemePart}}()
    for (i, v) in MO.gmp_variables(model)
        just[i] = []
        scheme_parts[i] = MO.approximation_scheme(model, v)
        for sp in scheme_parts[i]
            cref = MO.primal_scheme_constraint(amodel, sp, pvar[i])
            push!(just[i], cref)
        end
    end

    # add constraints
    pcon = Dict{Int, Vector{ConstraintRef}}()
    for (i, con) in MO.gmp_constraints(model)
        if shape(con) isa MO.MomentConstraintShape
            cref = @constraint amodel integrate(pvar, jump_function(con)) in moi_set(con)
            pcon[i] = [cref]
        elseif shape(con) isa MeasureConstraintShape
            # add measure constraints
            refmeas = moi_set(con)
            mons = monomials(maxdegree_basis(approx_basis(refmeas), variables(refmeas), approximation_degree(model)))
            pcon[i] = @constraint amodel sum(c.*(MM.expectation.(mons, pvar[index(v)])) for (c, v) in jump_function(con)) .== integrate.(mons, refmeas)
        end
    end
    return amodel, pvar
end

function get_moment(p::Polynomial, meas, model)
    expr = AffExpr(0.0)
    for i in 1:length(p.a)
        index = findfirst(isequal(p.x[i]), meas.x)
        add_to_expression!(expr, p.a[i], VariableRef(model, MOI.VariableIndex(index)))
    end
    return expr
end

function UpperDiag(vars, n)
    M = AffExpr.(zeros(n,n))
    k = 1
    for i in 1:n, j in i:n
        M[i,j] += vars[k]
        k += 1
    end
    return M
end


text:
function smolyak_knn_points(nodes::Array{T,2},x::Array{T,1},num_points::S) where {T<:AbstractFloat,S<:Integer}

  n = size(nodes,1)

  if num_points > n
    println("smolyak_knn_points: Too many points")
    throw(ErrorException())
  end

  dist_vector    = zeros(T,n)
  index_vector_k = Array{S}(num_points)
  dist_vector_k  = Array{T}(num_points)

  for i = 1:n
    dist_vector[i] = norm(x-nodes[i,:])^2
  end

  for i = 1:num_points
    (dist,index)       = findmin(dist_vector)
    dist_vector_k[i]   = dist
    index_vector_k[i]  = index
    dist_vector[index] = Inf
  end

  return index_vector_k, dist_vector_k

end

function smolyak_knn_evaluate(y::AbstractArray{T,1},nodes::Array{T,2},x::Array{T,1},k = 2*size(nodes,2)+1) where T<:AbstractFloat

  (index_vector_k,dist_vector_k) = smolyak_knn_points(nodes,x,k)

  weights = Array{T}(k)
  for i = 1:k
    weights[i] = 1.0/(dist_vector_k[i]+eps(T))
  end

  weights = weights./sum(weights)

  y_hat = zero(T)

  for i = 1:k
    y_hat += y[index_vector_k[i]]*weights[i]
  end

  return y_hat

end


text:
using MPI, Test

include(joinpath("..", "testhelpers.jl"))

@testset "Diagnostics" begin
    runmpi(joinpath(@__DIR__, "sin_test.jl"), ntasks = 2)
    runmpi(joinpath(@__DIR__, "dm_tests.jl"), ntasks = 2)
    runmpi(joinpath(@__DIR__, "Debug/test_statecheck.jl"), ntasks = 2)
end


text:
"""
    initialize_σ_Ising(L::Int, start_type::Symbol, rng::AbstractRNG)

Initialize the Ising model lattice `σ` with values ∈ {-1, 1}.

Arguments
* `L`         : Side length of the square lattice
* `start_type`: Either `:hot` or `:cold` to initialize the lattice with a hot or cold start
* `rng`       : Random number generator
Returns
* `σ`         : 2D array representing a square lattice
"""
function initialize_σ_Ising(L::Int, start_type::Symbol, rng::AbstractRNG)
    if start_type == :cold
        σ = ones(Int8, (L, L))
    elseif start_type == :hot
        σ = rand(rng, Int8[-1, 1], (L, L))
    else
        error("invalid start_type provided, valid choices are [:cold, :hot]")
    end

    return σ
end

"""
    initialize_σ_Potts(q::Int, L::Int, start_type::Symbol, rng::AbstractRNG)

Initialize the Potts model lattice `σ` with values ∈ {1, ..., q}.

Arguments
* `q`: Number of possible spin values
* `L`         : Side length of the square lattice
* `start_type`: Either `:hot` or `:cold` to initialize the lattice with a hot or cold start
* `rng`       : Random number generator
Returns
* `σ`         : 2D array representing a square lattice
"""
function initialize_σ_Potts(q::Int, L::Int, start_type::Symbol, rng::AbstractRNG)
    if start_type == :cold
        σ = ones(Int8, (L, L))
    elseif start_type == :hot
        σ = rand(rng, UnitRange{Int8}(1, q), (L, L))
    else
        error("invalid start_type provided, valid choices are [:cold, :hot]")
    end

    return σ
end

"""
    initialize_σ_XY(L::Int, start_type::Symbol, rng::AbstractRNG)

Initialize the XY model lattice `σ` with unit vectors represented by `NTuple{2, Float64}`.

Arguments
* `L`         : Side length of the square lattice
* `start_type`: Either `:hot` or `:cold` to initialize the lattice with a hot or cold start
* `rng`       : Random number generator
Returns
* `σ`         : 2D array representing a square lattice
    represented by a two-tuple
"""
function initialize_σ_XY(L::Int, start_type::Symbol, rng::AbstractRNG)
    if start_type == :cold
        σ = [(0.0, 1.0) for i = 1:L, j = 1:L]
    elseif start_type == :hot
        σ = [random_XYVector(rng) for i = 1:L, j = 1:L]
    else
        error("invalid start_type provided, valid choices are [:cold, :hot]")
    end

    return σ
end

"""
    random_XYVector(rng::AbstractRNG)

Generates a random unit vector represented by `NTuple{2, Float64}`.

Arguments
* `rng` : Random number generator
Returns
* A two-tuple representing a unit-vector with random angular orientation
"""
function random_XYVector(rng::AbstractRNG)
    θ = rand(rng) * 2π
    return (cos(θ), sin(θ))
end


text:
struct Error <: Exception
    str::String
    T
    code
end

Error(buf::AbstractString, T, code, pos, tlen) = Error(buf, T, code)
Error(buf::AbstractVector{UInt8}, T, code, pos, tlen) = Error(String(buf[pos:(pos + tlen - 1)]), T, code)

function Error(buf::IO, T, code, pos, tlen)
    fastseek!(buf, pos - 1)
    bytes = read(buf, tlen)
    return Error(String(bytes), T, code)
end

function Base.showerror(io::IO, e::Error)
    c = e.code
    println(io, "Parsers.Error ($(codes(c))):")
    println(io, text(c))
    println(io, "attempted to parse $(e.T) from: \"$(escape_string(e.str))\"")
end

"""
A bitmask value, with various bits corresponding to different parsing signals and scenarios.

`Parsers.xparse` returns a `code` value with various bits set according to the various scenarios
encountered while parsing a value.

    * `INVALID`: there are a number of invalid parsing states, all include the INVALID bit set (check via `Parsers.invalid(code)`)
    * `OK`: signals specifically that a valid value of type `T` was parsed (check via `Parsers.ok(code)`)
    * `SENTINEL`: signals that a valid sentinel value was detected while parsing, passed via the `sentinel` keyword argument to `Parsers.Options` (check via `Parsers.sentinel(code)`)
    * `QUOTED`: a `openquotechar` from `Parsers.Options` was detected at the beginning of parsing (check via `Parsers.quoted(code)`)
    * `DELIMITED`: a `delim` character or string from `Parsers.Options` was detected while parsing (check via `Parsers.delimited(code)`)
    * `NEWLINE`: a non-quoted newline character (`'\n'`), return character (`'\r'`), or CRLF (`"\r\n"`) was detected while parsing (check via `Parsers.newline(code)`)
    * `EOF`: the end of file was reached while parsing
    * `ESCAPED_STRING`: an `escapechar` from `Parsers.Options` was encountered while parsing (check via `Parsers.escapedstring(code)`)
    * `INVALID_QUOTED_FIELD`: a `openquotechar` were detected when parsing began, but no corresponding `closequotechar` were found to correctly close a quoted field, this is usually a fatal parsing error because parsing will continue until EOF to look for the close quote character (check via `Parsers.invalidquotedfield(code)`)
    * `INVALID_DELIMITER`: a `delim` character or string were eventually detected, but not at the expected position (directly after parsing a valid value), indicating there are extra, invalid characters between a valid value and the expected delimiter (check via `Parsers.invaliddelimiter(code)`)
    * `OVERFLOW`: overflow occurred while parsing a type, like `Integer`, that have limits on valid values (check via `Parsers.overflow(code)`)

One additional convenience function is provided, `Parsers.quotednotescaped(code)`, which checks if a value was quoted,
but didn't contain any escape characters, useful to indicate if a string may be used "as-is", instead of needing to be unescaped.
"""
const ReturnCode = Int16

const SUCCESS = 0b0000000000000000 % ReturnCode
const INVALID = 0b1000000000000000 % ReturnCode

# success flags
const OK                   = 0b0000000000000001 % ReturnCode
const SENTINEL             = 0b0000000000000010 % ReturnCode

# property flags
const QUOTED               = 0b0000000000000100 % ReturnCode
const DELIMITED            = 0b0000000000001000 % ReturnCode
const NEWLINE              = 0b0000000000010000 % ReturnCode
const EOF                  = 0b0000000000100000 % ReturnCode
const ESCAPED_STRING       = 0b0000001000000000 % ReturnCode

# invalid flags
const INVALID_QUOTED_FIELD = 0b1000000001000000 % ReturnCode
const INVALID_DELIMITER    = 0b1000000010000000 % ReturnCode
const OVERFLOW             = 0b1000000100000000 % ReturnCode
const INVALID_TOKEN        = 0b1000010000000000 % ReturnCode

ok(x::ReturnCode) = (x & (OK | INVALID)) == OK
invalid(x::ReturnCode) = x < SUCCESS
sentinel(x::ReturnCode) = (x & SENTINEL) == SENTINEL
quoted(x::ReturnCode) = (x & QUOTED) == QUOTED
delimited(x::ReturnCode) = (x & DELIMITED) == DELIMITED
newline(x::ReturnCode) = (x & NEWLINE) == NEWLINE
escapedstring(x::ReturnCode) = (x & ESCAPED_STRING) == ESCAPED_STRING
invalidquotedfield(x::ReturnCode) = (x & INVALID_QUOTED_FIELD) == INVALID_QUOTED_FIELD
invaliddelimiter(x::ReturnCode) = (x & INVALID_DELIMITER) == INVALID_DELIMITER
overflow(x::ReturnCode) = (x & OVERFLOW) == OVERFLOW
quotednotescaped(x::ReturnCode) = (x & (QUOTED | ESCAPED_STRING)) == QUOTED
invalidtoken(x::ReturnCode) = (x & INVALID_TOKEN) == INVALID_TOKEN
eof(x::ReturnCode) = (x & EOF) == EOF

memcmp(a::Ptr{UInt8}, b::Ptr{UInt8}, len::Int) = ccall(:memcmp, Cint, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), a, b, len) == 0

function checksentinel(source::AbstractVector{UInt8}, pos, len, sentinel)
    sentinelpos = 0
    startptr = pointer(source, pos)
    # sentinel is an iterable of Tuple{Ptr{UInt8}, Int}, sorted from longest sentinel string to shortest
    for (ptr, ptrlen) in sentinel
        if pos + ptrlen - 1 <= len
            match = memcmp(startptr, ptr, ptrlen)
            if match
                sentinelpos = pos + ptrlen
                break
            end
        end
    end
    return sentinelpos
end

function checksentinel(source::IO, pos, len, sentinel)
    sentinelpos = 0
    origpos = position(source)
    for (ptr, ptrlen) in sentinel
        matched = match!(source, ptr, ptrlen)
        if matched
            sentinelpos = pos + ptrlen
            break
        end
    end
    fastseek!(source, origpos)
    return sentinelpos
end

function checkdelim(source::AbstractVector{UInt8}, pos, len, (ptr, ptrlen))
    startptr = pointer(source, pos)
    delimpos = pos
    if pos + ptrlen - 1 <= len
        match = memcmp(startptr, ptr, ptrlen)
        if match
            delimpos = pos + ptrlen
        end
    end
    return delimpos
end

function checkdelim(source::IO, pos, len, (ptr, ptrlen))
    delimpos = pos
    matched = match!(source, ptr, ptrlen)
    if matched
        delimpos = pos + ptrlen
    end
    return delimpos
end

@inline function match!(source::IO, ptr, ptrlen)
    b = peekbyte(source)
    c = unsafe_load(ptr)
    if b != c
        return false
    elseif ptrlen == 1
        incr!(source)
        return true
    end
    pos = position(source)
    i = 1
    while true
        incr!(source)
        if Base.eof(source)
            fastseek!(source, pos)
            return false
        end
        b = peekbyte(source)
        c = unsafe_load(ptr + i)
        if b != c
            fastseek!(source, pos)
            return false
        elseif i == ptrlen - 1
            break
        end
        i += 1
    end
    incr!(source)
    return true
end

"""
    Parsers.fastseek!(io::IO, n::Integer)

    Without valididty checks, seek an `IO` to desired byte position `n`. Used in Parsers.jl to
    seek back to a previous location already parsed.
"""
function fastseek! end

fastseek!(io::IO, n::Integer) = seek(io, n)
function fastseek!(io::IOBuffer, n::Integer)
    io.ptr = n + 1
    return
end
fastseek!(io::AbstractVector{UInt8}, n::Integer) = nothing

"""
    Parsers.readbyte(io::IO)::UInt8

    Consume a single byte from an `IO` without checking `eof(io)`.
"""
function readbyte end

"""
    Parsers.peekbyte(io::IO)::UInt8

    Return, but do not consume, the next byte from an `IO` without checking `eof(io)`.
"""
function peekbyte end

incr!(io::IO) = readbyte(io)
readbyte(from::IO) = Base.read(from, UInt8)
peekbyte(from::IO) = UInt8(Base.peek(from))

function readbyte(from::IOBuffer)
    i = from.ptr
    @inbounds byte = from.data[i]
    from.ptr = i + 1
    return byte
end

function peekbyte(from::IOBuffer)
    @inbounds byte = from.data[from.ptr]
    return byte
end

function incr!(from::IOBuffer)
    from.ptr += 1
    return
end

incr!(::AbstractVector{UInt8}) = nothing
peekbyte(from::IO, pos) = peekbyte(from)
function peekbyte(from::AbstractVector{UInt8}, pos)
    @inbounds b = from[pos]
    return b
end

eof(::AbstractVector{UInt8}, pos::Integer, len::Integer) = pos > len
eof(source::IO, pos::Integer, len::Integer) = Base.eof(source)

function text(r)
    str = ""
    if r & QUOTED > 0
        str = "encountered an opening quote character, initial value parsing "
    else
        str = "initial value parsing "
    end
    if r & OK > 0
        str *= "succeeded"
    else
        str *= "failed"
    end
    if r & (~INVALID & OVERFLOW) > 0
        str *= ", value overflowed"
    end
    if r & SENTINEL > 0
        str *= ", a sentinel value was parsed"
    end
    if r & ESCAPED_STRING > 0
        str *= ", encountered escape character"
    end
    if r & (~INVALID & INVALID_QUOTED_FIELD) > 0
        str *= ", invalid quoted value"
    end
    if r & DELIMITED > 0
        str *= ", a valid delimiter was parsed"
    end
    if r & NEWLINE > 0
        str *= ", a newline was encountered"
    end
    if r & EOF > 0
        str *= ", reached EOF"
    end
    if r & (~INVALID & INVALID_DELIMITER) > 0
        str *= ", invalid delimiter"
    end
    return str
end

codes(r) = chop(chop(string(
    ifelse(r > 0, "SUCCESS: ", "INVALID: "),
    ifelse(r & OK > 0, "OK | ", ""),
    ifelse(r & SENTINEL > 0, "SENTINEL | ", ""),
    ifelse(r & QUOTED > 0, "QUOTED | ", ""),
    ifelse(r & ESCAPED_STRING > 0, "ESCAPED_STRING | ", ""),
    ifelse(r & DELIMITED > 0, "DELIMITED | ", ""),
    ifelse(r & NEWLINE > 0, "NEWLINE | ", ""),
    ifelse(r & EOF > 0, "EOF | ", ""),
    ifelse(r & (~INVALID & INVALID_QUOTED_FIELD) > 0, "INVALID_QUOTED_FIELD | ", ""),
    ifelse(r & (~INVALID & INVALID_DELIMITER) > 0, "INVALID_DELIMITER | ", ""),
    ifelse(r & (~INVALID & OVERFLOW) > 0, "OVERFLOW | ", "")
)))

ptrlen(s::String) = (pointer(s), sizeof(s))

"""
    PosLen(pos, len, ismissing, escaped)

A custom 64-bit primitive that allows efficiently storing the byte position
and length of a value within a byte array, along with whether a sentinel
value was parsed, and whether the parsed value includes escaped characters.
Specifically, the use of 64-bits is:
  * 1 bit to indicate whether a sentinel value was encountered while parsing
  * 1 bit to indicate whether the escape character was encountered while parsing
  * 42 bits to note the byte position as an integer where a value is located in a byte array (max array size ~4.4TB)
  * 20 bits to note the length of a parsed value (max length of ~1MB)

These individual "fields" can be retrieved via dot access, like `poslen.missingvalue`, `poslen.escapedvalue`,
`poslen.pos`, `poslen.len`.

`Parsers.xparse(String, buf, pos, len, opts)` returns `Parsers.Result{PosLen}`, where the `x.val` is a `PosLen`.
"""
primitive type PosLen 64 end

const MAX_POS = 4398046511104
const MAX_LEN = 1048575
@noinline postoolarge(pos) = throw(ArgumentError("position argument to Parsers.PosLen ($pos) is too large; max position allowed is $MAX_POS"))
@noinline lentoolarge(len) = throw(ArgumentError("length argument to Parsers.PosLen ($len) is too large; max length allowed is $MAX_LEN"))

@inline function PosLen(pos::Integer, len::Integer, ismissing=false, escaped=false)
    pos > MAX_POS && postoolarge(pos)
    len > MAX_LEN && lentoolarge(len)
    pos = Int64(pos) << 20
    pos |= ifelse(ismissing, MISSING_BIT, 0)
    pos |= ifelse(escaped, ESCAPE_BIT, 0)
    return Base.bitcast(PosLen, pos | Int64(len))
end

const MISSING_BIT = Base.bitcast(Int64, 0x8000000000000000)
const ESCAPE_BIT = Base.bitcast(Int64, 0x4000000000000000)
const POS_BITS = Base.bitcast(Int64, 0x3ffffffffff00000)
const LEN_BITS = Base.bitcast(Int64, 0x00000000000fffff)

@noinline invalidproperty(nm) = throw(ArgumentError("invalid property $nm for PosLen"))

function Base.getproperty(x::PosLen, nm::Symbol)
    y = Base.bitcast(Int64, x)
    nm === :pos && return (y & POS_BITS) >> 20
    nm === :len && return y & LEN_BITS
    nm === :missingvalue && return (y & MISSING_BIT) == MISSING_BIT
    nm === :escapedvalue && return (y & ESCAPE_BIT) == ESCAPE_BIT
    invalidproperty(nm)
end

"""
    Parsers.getstring(buf_or_io, poslen::PosLen, e::UInt8) => String

When calling `Parsers.xparse` with a `String` type argument, a `Parsers.Result{PosLen}` is returned, which has 3 fields:
  * `val`: a [`PosLen`](@ref) value which stores the starting byte position and length of the parsed string value
  * `code`: a parsing return code indicating success/failure
  * `tlen`: the total number of bytes parsed, which may differ from `val.len` if delimiters or open/close quotes were parsed

If the actual parsed `String` _is_ needed, however, you can pass your source and the `res.val::PosLen` to `Parsers.getstring`
to get the actual parsed `String` value.
"""
function getstring end

_unsafe_string(p, len) = ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8}, Int), p, len)

@inline function getstring(source::Union{IO, AbstractVector{UInt8}}, x::PosLen, e::UInt8)
    x.escapedvalue && return unescape(source, x, e)
    if source isa AbstractVector{UInt8}
        return _unsafe_string(pointer(source, x.pos), x.len)
    else
        vpos, vlen = x.pos, x.len
        fastseek!(source, vpos - 1)
        str = Base.StringVector(vlen)
        readbytes!(source, str, vlen)
        return String(str)
    end
end

getstring(str::AbstractString, poslen::PosLen, e::UInt8) = getstring(codeunits(str), poslen, e)

# if a cell value of a csv file has escape characters, we need to unescape it
@noinline function unescape(origbuf, x::PosLen, e)
    n = x.len
    if origbuf isa AbstractVector{UInt8}
        source = view(origbuf, x.pos:(x.pos + x.len - 1))
    else
        origpos = position(origbuf)
        fastseek!(origbuf, x.pos - 1)
        source = origbuf
    end
    out = Base.StringVector(n)
    len = 1
    i = 1
    @inbounds begin
        while i <= n
            b = peekbyte(source, i)
            if b == e
                incr!(source)
                i += 1
                b = peekbyte(source, i)
            end
            out[len] = b
            len += 1
            incr!(source)
            i += 1
        end
    end
    if origbuf isa IO
        fastseek!(origbuf, origpos)
    end
    resize!(out, len - 1)
    return String(out)
end


text:
using Compat.Test
using Compat.LinearAlgebra
using MathProgBase

function quadprogtest(solver)
    @testset "Testing quadprog with $solver" begin
        sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf,solver)
        @test sol.status == :Optimal
        @test isapprox(sol.objval, 130/70, atol=1e-6)
        @test isapprox(norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), 0.0, atol=1e-6)

        @testset "QP1" begin
            m = MathProgBase.LinearQuadraticModel(solver)
            MathProgBase.loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)

            MathProgBase.setquadobj!(m,diagm(0 => [10.0,10.0,10.0]))
            rows = [1, 2, 2, 2, 3, 3, 3]
            cols = [1, 1, 1, 2, 2, 3, 3]
            vals = Float64[2, 0.5, 0.5, 2, 1, 1, 1]
            MathProgBase.setquadobj!(m,rows,cols,vals)
            MathProgBase.optimize!(m)
            stat = MathProgBase.status(m)
            @test stat == :Optimal
            @test isapprox(MathProgBase.getobjval(m), 130/70, atol=1e-6)
            @test isapprox(norm(MathProgBase.getsolution(m) - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), 0.0, atol=1e-6)
        end

        @testset "QP2" begin
            m = MathProgBase.LinearQuadraticModel(solver)
            MathProgBase.loadproblem!(m, [-1. 1.; 1. 1.], [0.,0.], [Inf,Inf], [1.,1.], [0.,0.], [Inf,Inf], :Max)
            MathProgBase.addquadconstr!(m, [2], [1.], [1], [1], [1.], '<', 2)
            MathProgBase.optimize!(m)
            stat = MathProgBase.status(m)
            @test stat == :Optimal
            @test isapprox(MathProgBase.getobjval(m), 2.25, atol=1e-6)
            @test isapprox(norm(MathProgBase.getsolution(m) - [0.5,1.75]), 0.0, atol=1e-3)
        end
    end
end

function qpdualtest(solver)
    @testset "Testing QP duals with $solver" begin
        # max x
        # s.t. x^2 <= 2
        @testset "QP1" begin
            m = MathProgBase.LinearQuadraticModel(solver)
            MathProgBase.loadproblem!(m, zeros(0,1), [-Inf], [Inf], [1.0], Float64[], Float64[], :Max)
            MathProgBase.addquadconstr!(m, [], [], [1], [1], [1.0], '<', 2.0)
            MathProgBase.optimize!(m)
            stat = MathProgBase.status(m)

            @test MathProgBase.numlinconstr(m) == 0
            @test MathProgBase.numquadconstr(m) == 1
            @test MathProgBase.numconstr(m) == 1
            @test stat == :Optimal
            @test isapprox(MathProgBase.getobjval(m), sqrt(2), atol=1e-6)
            @test isapprox(MathProgBase.getsolution(m)[1], sqrt(2), atol=1e-6)
            @test isapprox(MathProgBase.getquadconstrduals(m)[1], 0.5/sqrt(2), atol=1e-6)
        end

        # min -x
        # s.t. x^2 <= 2
        @testset "QP2" begin
            m = MathProgBase.LinearQuadraticModel(solver)
            MathProgBase.loadproblem!(m, zeros(0,1), [-Inf], [Inf], [-1.0], Float64[], Float64[], :Min)
            MathProgBase.addquadconstr!(m, [], [], [1], [1], [1.0], '<', 2.0)
            MathProgBase.optimize!(m)
            stat = MathProgBase.status(m)

            @test MathProgBase.numlinconstr(m) == 0
            @test MathProgBase.numquadconstr(m) == 1
            @test MathProgBase.numconstr(m) == 1
            @test stat == :Optimal
            @test isapprox(MathProgBase.getobjval(m), -sqrt(2), atol=1e-6)
            @test isapprox(MathProgBase.getsolution(m)[1], sqrt(2), atol=1e-6)
            @test isapprox(MathProgBase.getquadconstrduals(m)[1], -0.5/sqrt(2), atol=1e-6)
        end
    end
end

function socptest(solver)
    @testset "Testing SOCP interface with $solver" begin
        # min t
        # s.t. x + y >= 1
        #      x^2 + y^2 <= t^2
        #      t >= 0
        m = MathProgBase.LinearQuadraticModel(solver)
        MathProgBase.loadproblem!(m, [ 1.0 1.0 0.0 ], [-Inf,-Inf,0.0], [Inf,Inf,Inf], [0.0,0.0,1.0], [1.0],[Inf], :Min)
        MathProgBase.addquadconstr!(m, [], [], [1,2,3], [1,2,3], [1.0,1.0,-1.0],'<',0.0)
        MathProgBase.optimize!(m)
        stat = MathProgBase.status(m)

        @test stat == :Optimal
        @test isapprox(MathProgBase.getobjval(m), sqrt(1/2), atol=1e-6)
        @test isapprox(norm(MathProgBase.getsolution(m) - [0.5,0.5,sqrt(1/2)]), 0.0, atol=1e-3)
    end
end


text:
"""
    apply!(σ, changes, [changes_r])

Applies the changes `changes` to the `σ`.

If `state isa DoubleState` then single-value changes
are applied to the columns of the state (in order to
compute matrix-operator products). Otherwise it should
be a tuple with changes of row and columns.

If the state is double but there is only 1 element of changes,
it's applied to the rows.

If changes is nothing, does nothing.
"""
apply!(σ::AbstractArray, cngs::Nothing) = σ

@inline function apply!(σ::ADoubleState, cngs_l::Union{StateChanges,Nothing})
    apply!(row(σ), cngs_l)
    return σ
end

@inline function apply!(σ::ADoubleState, (cngs_l, cngs_r))
    apply!(row(σ), cngs_l)
    apply!(col(σ), cngs_r)
    return σ
end

function apply!(σ::ADoubleState, cngs_l, cngs_r)
    apply!(row(σ), cngs_l)
    apply!(col(σ), cngs_r)
    return σ
end

@inline apply!(σ::AbstractVector, cngs::Nothing) = σ
function apply!(σ::AbstractVector, cngs)
    for (site, val)=cngs
        σ[site] = val
    end
    return σ
end

"""
    apply(state, cngs)

Applies the changes `cngs` to the state `σ`, by allocating a
copy.

See also @ref(apply!)
"""
apply(σ::AbstractState, cngs) = apply!(deepcopy(σ), cngs)


text:
#=
$title Sensitivity Analysis in Economic Load Dispatch

$onText
For more details please refer to Chapter 3 (Gcode3.2), of the following book:
Soroudi, Alireza. Power System Optimization Modeling in GAMS. Springer, 2017.
--------------------------------------------------------------------------------
Model type: QCP
--------------------------------------------------------------------------------
Contributed by
Dr. Alireza Soroudi
IEEE Senior Member
Email: alireza.soroudi@gmail.com
We do request that publications derived from the use of the developed GAMS code
explicitly acknowledge that fact by citing
Soroudi, Alireza. Power System Optimization Modeling in GAMS. Springer, 2017.
DOI: doi.org/10.1007/978-3-319-62350-4
$offText

Converstion from GAMS to JuMP by Jon Duan
=#

N_Gen = 5
Gen_set = 1:N_Gen
Gen = [Symbol("G$i") for i in Gen_set]

N_counter = 11
counter = 1:N_counter

using DataFrames

Load = 400

Other_set = ["a"     ,"b",      "c",       "Pmin",  "Pmax"]
Other = [Symbol("$i") for i in Other_set]

# https://github.com/davidavdav/NamedArrays.jl
# using DataStructures
using NamedArrays
data = Float64[  [3     20     100     28    206];
                 [4.05  18.07  98.87   90    284];
                 [4.05  15.55  104.26  68    189];
                 [3.99  19.21  107.21  76    266];
                 [3.88  26.18  95.31   19    53]]

data =  NamedArray(data, ( Gen, Other ), ("Rows", "Cols"))
# df[:y2] = parse.(Int,df[:yearsAsString])

using JuMP
using Gurobi

EDsensitivity = Model(solver=GurobiSolver())

@variable(EDsensitivity,  data[i, :Pmin]<= P[i in Gen]<=data[i, :Pmax])

@variable(EDsensitivity, OF>=0)

@objective(EDsensitivity, Min, sum( data[i, :a]*P[i]^2 + data[i,:b]*P[i] for i in Gen))

@constraint(EDsensitivity, Balance, sum( P[i]for i in Gen ) >= Load)

# post Optimization
# initialized
report = DataFrame(OF=Float64[], Load=Float64[])
repGen = convert(DataFrame,Dict([i => Float64[]  for i in Gen]))

# Modifying constraints
# JuMP does not currently support changing constraint coefficients. For less-than and greater-than constraints, the right-hand-side can be changed, e.g.:
for i in counter

    # update the model with increasing Load
    Load =sum(data[i,:Pmin] for i in Gen) +
        ((i -1)/(N_counter-1))*sum(data[i,:Pmax] - data[i,:Pmin] for i in Gen)
    JuMP.setRHS(Balance, Load)  # Now

    # Solve the problem
    solve(EDsensitivity)

    # Print our final solution
    println("Objective value: ", getobjectivevalue(EDsensitivity))

    println("P = ", getvalue(P))
    push!(repGen, [getvalue(P)[i] for i in Gen])
    println("Marginal Value of Balance = ", getdual(Balance))
    # check balance
    sum(getvalue(P))
    # check margian value for P
    getdual(P)
    push!(report, [getobjectivevalue(EDsensitivity), Load])

end

# https://jump.readthedocs.io/en/latest/probmod.html
# JuMP will use the ability to modify problems exposed by the solver if possible, and will still work even if the solver does not support this functionality by passing the complete problem to the solver every time.

# Pkg.add("StatPlots")
using DataFrames, Plots, StatPlots, CSV

# http://docs.juliaplots.org/latest/tutorial/

@df report plot(:Load, :OF)

plot(report[:Load], report[:OF])

@df repGen plot(counter, cols(1:5), label = names(repGen))