text:
################################################################################
# Utility functions
################################################################################

# Indexing

srcindFITS(::los_radial) = 1
srcindFITS(::los_earth) = 1:2

srcindG(::los_radial) = ()
srcindG(::los_earth) = (0:1,)

obsindFITS(::los_radial) = 1
obsindFITS(::los_earth) = 1:2

obsindG(::los_radial) = ()
obsindG(::los_earth) = (0:1,)

sizeindG(::los_radial) = ()
sizeindG(::los_earth) = (2,)

αrcomp(G::AbstractArray{ComplexF64, 2}, r_ind, α) = G[α, r_ind]
αrcomp(G::AbstractArray{ComplexF64, 3}, r_ind, α) = G[α, 0, r_ind]
αrcomp(G::AbstractArray{ComplexF64, 2}, r_ind, ::los_radial) = αrcomp(G, r_ind, 0)
αrcomp(G::AbstractArray{ComplexF64, 3}, r_ind, ::los_radial) = αrcomp(G, r_ind, 0)
function αrcomp(G::AbstractArray{ComplexF64, 2}, r_ind, los::los_earth)
	αrcomp(G, r_ind, Base.IdentityUnitRange(0:1))
end
function αrcomp(G::AbstractArray{ComplexF64, 3}, r_ind, los::los_earth)
	αrcomp(G, r_ind, Base.IdentityUnitRange(0:1))
end

##################################################################################
# Compute the processor range for a particular range of modes
##################################################################################

function ℓ′ω_range_from_ℓω_range(ℓ_ωind_iter_on_proc::ProductSplit, s_max, ℓ_arr)
	ω_ind_min, ω_ind_max = extremaelement(ℓ_ωind_iter_on_proc, dims = 2)
	ℓ′min_ωmin = minimum(minimum(intersect(ℓ_arr, max(ℓ-s_max, 0):ℓ+s_max))
		for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc if ω_ind==ω_ind_min)
	ℓ′max_ωmax = maximum(maximum(intersect(ℓ_arr, max(ℓ-s_max, 0):ℓ+s_max))
		for (ℓ, ω_ind) in ℓ_ωind_iter_on_proc if ω_ind==ω_ind_max)
	return (ℓ′min_ωmin, ω_ind_min),(ℓ′max_ωmax, ω_ind_max)
end

function Gfn_processor_range((ℓ_arr, ω_ind_arr), ℓ_ωind_iter_on_proc::ProductSplit, s_max, num_procs)
	modes = ℓ′ω_range_from_ℓω_range(ℓ_ωind_iter_on_proc, s_max, ℓ_arr)
	proc_id_min = ParallelUtilities.whichproc((ℓ_arr, ω_ind_arr), first(modes), num_procs)
	proc_id_max = ParallelUtilities.whichproc((ℓ_arr, ω_ind_arr), last(modes), num_procs)
	return proc_id_min:proc_id_max
end

###################################################################################
## Read Green function slices (real and imaginary parts) and assign to arrays
###################################################################################

function read_Gfn_file_at_index!(G::AbstractArray{<:Complex},
	G_fits_dict::Dict, (ℓ_arr, ω_ind_arr)::Tuple{Vararg{AbstractUnitRange, 2}},
	ℓω::NTuple{2, Int}, num_procs::Integer, I...)

	v = ParallelUtilities.whichproc_localindex((ℓ_arr, ω_ind_arr), ℓω, num_procs)
	v === nothing && error("could not locate ℓω = $ℓω")
	proc_id_mode, ℓω_index_in_file = v
	G_file = G_fits_dict[proc_id_mode].hdu
	read_Gfn_file_at_index!(G, G_file, I..., ℓω_index_in_file)
	return G
end

function read_Gfn_file_at_index!(G::AbstractArray{<:Complex,0}, G_hdu::ImageHDU, I...)
	G1D = reshape(G, 1)
	read!(G_hdu, reinterpret_as_float(G1D), :, I...)
	return G
end
function read_Gfn_file_at_index!(G, G_hdu::ImageHDU, I...)
	read!(G_hdu, reinterpret_as_float(G), :, I...)
	return G
end

function Gfn_fits_files(path::String, proc_id_range::AbstractUnitRange)
	function f(procid, path)
		filepath = joinpath(path, @sprintf "Gfn_proc_%03d.fits" procid)
		f_FITS = FITS(filepath,"r") :: FITS
		f_HDU = f_FITS[1] :: ImageHDU{Float64,6}
		(fits = f_FITS, hdu = f_HDU)
	end
	Dict(procid => f(procid, path) for procid in proc_id_range)
end

function Gfn_fits_files(path::String, (ℓ_arr, ω_ind_arr), ℓ_ωind_iter_on_proc::ProductSplit, num_procs::Integer)
	proc_range = ParallelUtilities.procrange_recast((ℓ_arr, ω_ind_arr), ℓ_ωind_iter_on_proc, num_procs)
	Gfn_fits_files(path, proc_range)
end

function Gfn_fits_files(path::String, (ℓ_arr, ω_ind_arr), ℓ_ωind_iter_on_proc::ProductSplit, s_max::Integer, num_procs::Integer)
	proc_range = Gfn_processor_range((ℓ_arr, ω_ind_arr), ℓ_ωind_iter_on_proc, s_max, num_procs)
	Gfn_fits_files(path, proc_range)
end

closeGfnfits(d::Dict) = map(x -> close(x.fits), values(d))

# Directory to save output to
function Gfn_path_from_source_radius(r_src::Real; c_scale = 1)
	dir = "Greenfn_src$((r_src/Rsun > 0.99 ?
			(@sprintf "%dkm" (r_src-Rsun)/1e5) :
			(@sprintf "%.2fRsun" r_src/Rsun)))"
	if c_scale != 1
		dir *= "_c_scale_$(@sprintf "%g" c_scale)"
	end
	dir *= "_flipped"
	joinpath(SCRATCH[], dir)
end
Gfn_path_from_source_radius(x::Point3D; kwargs...) = Gfn_path_from_source_radius(x.r; kwargs...)

#####################################################################################
# Solve for components
#####################################################################################

function load_solar_model()

	modelS_meta = readdlm("$(@__DIR__)/ModelS.meta", comments = true, comment_char='#')::Matrix{Float64};
	Msun, Rsun = modelS_meta;

	modelS = readdlm("$(@__DIR__)/ModelS", comments = true, comment_char='#')::Matrix{Float64};
	modelS_detailed = readdlm("$(@__DIR__)/ModelS.detailed", comments = true, comment_char='#')::Matrix{Float64};

	HMIfreq = readdlm("$(@__DIR__)/m181q.1216")::Matrix{Float64};
	ℓ_HMI::Vector{Int}, ν_HMI::Vector{Float64}, γ_HMI::Vector{Float64} = HMIfreq[:, 1], HMIfreq[:, 3], HMIfreq[:, 5];

	# Fit modes above ℓ=11 and 2mHz<ν<4mHz
	mode_filter = @. (ℓ_HMI > 11) & (ν_HMI > 1.5e3) & (ν_HMI < 4.5e3);

	ν_HMI = ν_HMI[mode_filter];
	ℓ_HMI = ℓ_HMI[mode_filter];
	γ_HMI = γ_HMI[mode_filter];

	# Fit γ(ω) in Hz, the HMI frequencies are in μHz
	γ_damping = polyfit(ν_HMI.*(2π*1e-6), γ_HMI.*(2π*1e-6), 3);

	# invert the grid to go from inside to outside,
	# and leave out the center to avoid singularities
	# Start from r = 0.2Rsun to compare with Mandal. et al 2017

	#r_start_ind_mandal = searchsortedfirst(modelS[:, 1], 0.2, rev = true)
	r_start_ind_skipzero::Int = searchsortedfirst(modelS[:, 1], 2e-2, rev = true)
	#r_start_ind_somewhere_inside = searchsortedfirst(modelS[:, 1], 0.5, rev = true)

	#set r_start_ind to r_start_ind_mandal to compare with Mandal. et al 2017
	r_start_ind = r_start_ind_skipzero

	r::Vector{Float64} = modelS[r_start_ind:-1:1, 1]*Rsun;

	nr = length(r);

	dr = D(nr; stencil_gridpts = Dict(7=>3, 43=>5))*r;
	ddr = sparse(dbydr(dr));
	c = modelS[r_start_ind:-1:1, 2]::Vector{Float64};
	ρ = modelS[r_start_ind:-1:1, 3]::Vector{Float64};

	G = 6.67428e-8 # cgs units
	m = @. Msun*exp(modelS_detailed[r_start_ind:-1:1, 2]::Vector{Float64})

	g = @. G*m/r^2

	N2 = @. g * modelS_detailed[r_start_ind:-1:1, 15]::Vector{Float64} / r

	return Rsun, nr, r, dr, ddr, c,ρ, g, N2, γ_damping
end

const Rsun, nr, r, dr, ddr, c, ρ, g, N2, γ_damping = load_solar_model()

radial_grid_index(x::SphericalPoint) = radial_grid_index(radius(x))
radial_grid_index(r_pt::Real) = radial_grid_index(r, r_pt)
radial_grid_index(r, r_pt::Real) = searchsortedfirst(r, r_pt)
radial_grid_closest(x::SphericalPoint) = radial_grid_closest(radius(x))
radial_grid_closest(r_pt::Real) = radial_grid_closest(r, r_pt)
radial_grid_closest(r, r_pt::Real) = r[radial_grid_index(r, r_pt)]

const r_src_default = radial_grid_closest(r, Rsun - 75e5)
const r_obs_default = radial_grid_closest(r, Rsun + 150e5)

radius(x::Point3D) = x.r
radius(::Point2D) = r_obs_default

function read_rsrc_robs_c_scale(kwargs)
	r_src = get(kwargs, :r_src, r_src_default)
	r_obs = get(kwargs, :r_obs, r_obs_default)
	c_scale = get(kwargs, :c_scale, 1)
	r_src, r_obs, c_scale
end

# Source components in Hansen VSH basis
function σsrc_grid(r_src = r_src_default)
	r_src_ind = radial_grid_index(r, r_src)
	r_src_on_grid = r[r_src_ind];
	max(r_src_on_grid - r[max(1, r_src_ind-2)], r[min(nr, r_src_ind+2)] - r_src_on_grid)
end

delta_fn_gaussian_approx(r_src, σsrc) = @. exp(-(r - r_src)^2 / 2σsrc^2) / √(2π*σsrc^2) /r^2;

Ω(ℓ, N) = √((ℓ+N)*(ℓ-N+1)/2)
ζjₒjₛs(j₁, j₂, ℓ) = (Ω(j₁, 0)^2 + Ω(j₂, 0)^2 - Ω(ℓ, 0)^2)/(Ω(j₁, 0)*Ω(j₂, 0))
Njₒjₛs(jₒ, jₛ, s) = √((2jₒ+1)*(2jₛ+1)/(4π*(2s+1)))

# Functions to write the frequency grid and other parameters to the data directory

function frequency_grid(Gfn_save_directory; kwargs...)
	ν_low = get(kwargs, :ν_low, 2.0e-3)
	ν_high = get(kwargs, :ν_high, 4.5e-3)
	num_ν=get(kwargs, :num_ν, 1250)
	@assert(num_ν>1,"Need at least two points in frequency to construct the grid")
	ν_Nyquist = get(kwargs, :ν_Nyquist, 16e-3)

	ℓ_arr = get(kwargs, :ℓ_arr, 1:1)

	dν = (ν_high - ν_low)/(num_ν-1); dω = 2π*dν

	# choose values on a grid
	ν_low_index = Int64(floor(ν_low/dν)); ν_low = ν_low_index*dν
	ν_high_index = num_ν + ν_low_index - 1; ν_high = ν_high_index*dν;
	Nν_Gfn = num_ν
	ν_Nyquist_index = Int64(ceil(ν_Nyquist/dν)); ν_Nyquist = ν_Nyquist_index*dν

	Nν = ν_Nyquist_index + 1; Nt = 2*(Nν-1)
	ν_full = (0:ν_Nyquist_index).*dν;
	ν_arr = (ν_low_index:ν_high_index).*dν ;
	T = 1/dν; dt = T/Nt;
	ν_start_zeros = ν_low_index # index starts from zero
	ν_end_zeros = ν_Nyquist_index - ν_high_index

	ω_arr = 2π .* ν_arr;

	if !isdir(Gfn_save_directory)
		mkdir(Gfn_save_directory)
	end

	@save(joinpath(Gfn_save_directory,"parameters.jld2"),
		ν_arr, ω_arr, ν_full, dν, dω, ℓ_arr,
		ν_start_zeros, ν_end_zeros, Nν, Nt, dt, T, Nν_Gfn, ν_Nyquist)

	ℓ_arr, ω_arr
end

function append_parameters(Gfn_save_directory; kwargs...)
	paramfile = joinpath(Gfn_save_directory,"parameters.jld2")
	params = jldopen(paramfile,"a+")
	for (k, v) in Dict(kwargs)
		params[string(k)] = v
	end
	close(params)
end

function update_parameters(Gfn_save_directory; kwargs...)
	paramfile = joinpath(Gfn_save_directory,"parameters.jld2")
	params = load(paramfile)
	for (k, v) in Dict(kwargs)
		params[string(k)] = v
	end
	save(paramfile, params)
end

struct ParamsGfn
	path::String
	Nν_Gfn::Int
	ν_Nyquist::Float64
	ℓ_arr::UnitRange{Int}
	dν::Float64
	dω::Float64
	Nν::Int
	ν_full::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
	ν_arr::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
	T::Float64
	dt::Float64
	Nt::Int
	ν_start_zeros::Int
	ν_end_zeros::Int
	ω_arr::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
	num_procs::Int
end

function read_all_parameters(; kwargs...)
	r_src = get(kwargs, :r_src, r_src_default)
	c_scale = get(kwargs, :c_scale, 1)
	Gfn_path_src = Gfn_path_from_source_radius(r_src, c_scale = c_scale)
	# load all parameters from the file at the path
	params_all = load(joinpath(Gfn_path_src,"parameters.jld2"))
	# pack the parameters including the path into a struct
	ParamsGfn(Gfn_path_src,
		[params_all[String(k)] for k in fieldnames(ParamsGfn) if k != :path]...)
end

read_all_parameters(p::ParamsGfn; kwargs...) = p
function read_all_parameters(::Nothing; kwargs...)
	read_all_parameters(; kwargs...) :: ParamsGfn
end

function read_parameters(params...; kwargs...)
	isempty(params) && return ()
	params_all = read_all_parameters(; kwargs...)
	Tuple(getproperty(params_all, Symbol(p)) for p in params)
end

get_numprocs(path) = load(joinpath(path, "parameters.jld2"), "num_procs")

# Functions to compute Green function components in the helicity basis
# We solve in the Hansen VSH basis first and change over to the PB basis

function solve_for_components!(M, S, α, β)
	H = M\S # solve the equation

	@views @. α = H[1:nr]
	@views @. β = H[nr+1:end]
    return α, β
end

function source!(Sh, ω, ℓ, delta_fn)
    @views @. Sh[nr+1:end-1] = -√(ℓ*(ℓ+1))/ω^2 * delta_fn[1:nr-1] /(r[1:nr-1]*ρ[1:nr-1])
	Sh[end] = 0
    return Sh;
end

function ℒr!(L12, L21, c′; stencil_gridpts = Dict(6=>3, 42=>5), kw...)
    # Dirichlet condition β(r_out)=0 on the right edge
    dbydr!(L12, dr, stencil_gridpts = stencil_gridpts, left_edge_npts = 2, right_edge_npts = 2)
	L12[1:2] .= 0 # dirichlet boundary condition on G at the inner boundary

    dinds = diagind(L12)
    @. L12[dinds] += g/c′^2

    # Right edge ghost to use point α(r_out) in computing derivatives
    # Dirichlet condition α(r_in)=0 on the left edge
    dbydr!(L21, dr, stencil_gridpts = stencil_gridpts, left_edge_npts = 2, right_edge_npts = 2)
	L21[end-1:end] .= 0

    dinds = diagind(L21)
    @. L21[dinds] += 2/r - g/c′^2

	return L12, L21
end

function wave_operator_diagonal!(Mdiag, Mdiagfixed, ω, ℓ, c′, invρr²; kw...)
	ωC = ω - im * γ_damping(ω)

	Mdiag .= Mdiagfixed

	for (dind, rind) in zip(1:nr, 1:nr)
		Mdiag[dind] += -ωC^2 * ρ[rind]
	end
	for (dind, rind) in zip(nr+1:2nr, 1:nr)
		Mdiag[dind] += -ℓ*(ℓ+1) / ωC^2 * invρr²[rind]
	end

	return Mdiag
end

function greenfn_components_onemode!(ω, ℓ, αrℓω, βrℓω, αhℓω, βhℓω, M_ωℓ, Mdiagfixed, Mdiag,
	Sr, Sh, delta_fn, c′, invρr²; kwargs...)

	tangential_source = get(kwargs, :tangential_source, true)

	source!(Sh, ω, ℓ, delta_fn)

	wave_operator_diagonal!(Mdiag, Mdiagfixed, ω, ℓ, c′, invρr²; kwargs...);
	M_ωℓ[diagind(M_ωℓ)] .= Mdiag

    M_sparse = sparse(M_ωℓ)
	M_lu = lu(M_sparse)

	# radial source
	solve_for_components!(M_lu, Sr, αrℓω, βrℓω)

	if tangential_source
		solve_for_components!(M_lu, Sh, αhℓω, βhℓω)
	end
    return nothing
end

function greenfn_components_somemodes_serial_oneproc(ℓ_ωind_proc::ProductSplit,
	r_src, c_scale, ω_arr,
	tangential_source::Bool = true,
	tracker::Union{Nothing, RemoteChannel} = nothing,
	)

	# save real and imaginary parts separately
	β = 0:1
	γ = tangential_source ? (0:1) : (0:0)
	G = zeros(ComplexF64, nr, β, γ, 0:1, length(ℓ_ωind_proc))

	r_src_on_grid = radial_grid_closest(r, r_src)
	σsrc = σsrc_grid(r_src_on_grid)
	delta_fn = delta_fn_gaussian_approx(r_src_on_grid, σsrc)
	Sr, Sh = zeros(2nr), zeros(2nr);
	@views @. Sr[2:nr] = delta_fn[2:nr];

	αr, βr = zeros(ComplexF64, nr), zeros(ComplexF64, nr);
	αh, βh = zeros(ComplexF64, nr), zeros(ComplexF64, nr);

	# temporary arrys used to compute the derivative operator
	M_ωℓ = zeros(ComplexF64, 2nr, 2nr);
	L12 = zeros(nr, nr);
	L21 = zero(L12);
	Mdiagfixed = zeros(2nr); # the terms that do not depend on ω and ℓ
	Mdiag = zeros(ComplexF64, 2nr); # the terms that do not depend on ω and ℓ

	stencil_gridpts = Dict(6=>3, 42=>5);

	if c_scale != 1
		c′= c .* c_scale
		N2′= @. N2 + g^2/c^2 * (1-1/c_scale^2)
	else
		c′ = c;
		N2′ = N2;
	end
	invρr² = @. 1 / (ρ * r^2);

	@. Mdiagfixed[1:nr] = ρ * N2′;
	@. Mdiagfixed[nr+1:end] = 1 / (ρ * c′^2);

	Gωℓ, drGωℓ = zeros(ComplexF64, nr), zeros(ComplexF64, nr)

	ℒr!(L12, L21, c′, stencil_gridpts = stencil_gridpts)

	# top right, ∂ᵣβ, (nr x nr)
	M_ωℓ[1:nr, nr+1:end] = L12

	# bottom left ∂ᵣα, (nr x nr)
	M_ωℓ[nr+1:end, 1:nr] = L21

	for (ℓωind, (ℓ, ω_ind)) in enumerate(ℓ_ωind_proc)

		ω = ω_arr[ω_ind]

		greenfn_components_onemode!(ω, ℓ, αr,
			βr, αh, βh, M_ωℓ, Mdiagfixed, Mdiag, Sr, Sh, delta_fn,
			c′, invρr²,
			r_src = r_src, tangential_source = tangential_source,
			c_scale = c_scale);

		# radial component for radial source
		@. G[:, 0, 0, 0, ℓωind] = αr

		# tangential component for radial source
		@. G[:, 1, 0, 0, ℓωind] = Ω(ℓ, 0)/(ρ*r*ω^2) * βr

		if tangential_source
			# radial component for tangential source
			@. G[:, 0, 1, 0, ℓωind] = αh/√2

			# tangential component for tangential source
			@. G[:, 1, 1, 0, ℓωind] = Ω(ℓ, 0)/(ρ*r*ω^2) * (βh - delta_fn)
		end

		for srcind in UnitRange(axes(G,3)), obsind in UnitRange(axes(G,2))
			Gωℓ .= @view G[:, obsind, srcind, 0, ℓωind]
			mul!(drGωℓ, ddr, Gωℓ)
			G[:, obsind, srcind, 1, ℓωind] = drGωℓ
		end

		(tracker isa RemoteChannel) && put!(tracker, true)
	end

	permutedims(reshape(reinterpret_as_float(G), 2, nr, 2, 2, 2, :), [1, 3, 4, 2, 5, 6])
end

function greenfn_components_somemodes_serial_oneproc_fits(ℓ_ωind_proc::ProductSplit,
	r_src, c_scale,
	ω_arr, Gfn_save_directory,
	tangential_source::Bool = true,
	tracker::Union{Nothing, RemoteChannel} = nothing
    )

	rank = ParallelUtilities.workerrank(ℓ_ωind_proc)
	save_path = joinpath(Gfn_save_directory, @sprintf "Gfn_proc_%03d.fits" rank)

	FITS(save_path, "w") do file
		GR = greenfn_components_somemodes_serial_oneproc(
			ℓ_ωind_proc, r_src, c_scale, ω_arr, tangential_source, tracker)

		write(file, GR)
	end

	nothing
end

function greenfn_components(r_src = r_src_default; kwargs...)
	tangential_source = get(kwargs, :tangential_source, true)
	c_scale = get(kwargs, :c_scale, 1)
	Gfn_save_directory = Gfn_path_from_source_radius(r_src, c_scale = c_scale)

	if !isdir(Gfn_save_directory)
		@info "Creating $Gfn_save_directory"
		mkdir(Gfn_save_directory)
	end

	@info "Saving output to $Gfn_save_directory"

	ℓ_arr, ω_arr = frequency_grid(Gfn_save_directory; kwargs...);
	Nν_Gfn = length(ω_arr); ω_ind_arr = 1:Nν_Gfn
	@info "Using ℓ_arr = $(ℓ_arr) and $(Nν_Gfn) frequency bins"

	modes_iter = Iterators.product(ℓ_arr, ω_ind_arr)
	num_tasks = length(modes_iter)
	num_procs = min(num_tasks, nworkers())

	append_parameters(Gfn_save_directory, num_procs = num_procs)

	tracker = RemoteChannel(() -> Channel{Bool}(num_tasks + 1))
	prog_bar = Progress(num_tasks, 2, "Green functions computed : ")

	wrapper(x) = greenfn_components_somemodes_serial_oneproc_fits(x, r_src, c_scale, ω_arr,
	Gfn_save_directory, tangential_source, tracker)

	@sync begin
	 	@async begin
	 		try
	 			pmapbatch_productsplit(wrapper, ℓ_arr, ω_ind_arr);
	 			println("Finished computing Green functions")
	 		finally
	 			put!(tracker, false)
	 			finish!(prog_bar)
	 		end
	 	end
	 	while take!(tracker)
	 		 next!(prog_bar)
	 	end
	end;
	nothing # to suppress the task done message
end

function Gfn_reciprocity_partial(localtimer, ℓ_ωind_iter_on_proc, r_src, r_obs)
	ℓ_arr, ν_arr = ParallelUtilities.getiterators(ℓ_ωind_iter_on_proc)
	Nν_Gfn = length(ν_arr)

	r_src_ind = radial_grid_index(r, r_src)
	r_obs_ind = radial_grid_index(r, r_obs)
	Gfn_path_src = Gfn_path_from_source_radius(r_src)
	Gfn_path_obs = Gfn_path_from_source_radius(r_obs)
	num_procs_src = get_numprocs(Gfn_path_src)
	num_procs_obs = get_numprocs(Gfn_path_obs)

	G_reciprocity = zeros(2, ℓ_arr, ν_arr)

	Gfn_fits_files_src, Gfn_fits_files_obs =
		Gfn_fits_files.((Gfn_path_src, Gfn_path_obs), ((ℓ_arr, 1:Nν_Gfn),),
			(ℓ_ωind_iter_on_proc,), (num_procs_src, num_procs_obs))

	G01 = zeros(ComplexF64)
	G10 = zeros(ComplexF64)

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

		read_Gfn_file_at_index!(G10,
			Gfn_fits_files_src, (ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), num_procs_src,
			2, 1, r_obs_ind, 1)

		read_Gfn_file_at_index!(G01,
			Gfn_fits_files_obs, (ℓ_arr, 1:Nν_Gfn), (ℓ, ω_ind), num_procs_obs,
			1, 2, r_src_ind, 1)

		G_reciprocity[1, ℓ, ω_ind] = abs(G10[])
		G_reciprocity[2, ℓ, ω_ind] = abs(G01[])
	end

	return parent(G_reciprocity)
end

function Gfn_reciprocity(comm; kwargs...)
	r_src = get(kwargs, :r_src, r_src_default)
	r_obs = get(kwargs, :r_obs, r_obs_default)
	Gfn_path_src = Gfn_path_from_source_radius(r_src)
	@load(joinpath(Gfn_path_src,"parameters.jld2"),
		ν_arr, Nν_Gfn, ℓ_arr, num_procs)

	ℓ_range = get(kwargs, :ℓ_range, ℓ_arr)
	ν_ind_range = get(kwargs, :ν_ind_range, 1:Nν_Gfn)

	G_reciprocity = pmapsum(comm, Gfn_reciprocity_partial, (ℓ_range, ν_ind_range), r_src, r_obs)
	permutedims(G_reciprocity, [3,2,1])
end

#################################################################################
# Functions that are used in computing the kernels
#################################################################################

_mulR(a::Complex, b::Complex) = real(a)*real(b) - imag(a)*imag(b)
_mulI(a::Complex, b::Complex) = real(a)*imag(b) + imag(a)*real(b)

_structarrayparent(S::StructArray{<:Complex}) = StructArray{eltype(S)}(map(parent, (S.re, S.im)))
_structarrayparent(S::StructArray{<:Complex}, ::Type{T}) where {T} = StructArray{eltype(S)}(map(T∘parent, (S.re, S.im)))
_structarrayparent(S::AbstractArray, ::Type{T}) where {T} = T(no_offset_view(parent(S)))
_structarray(S::StructArray{<:Complex}, ::Type{T}) where {T} = StructArray{eltype(S)}(map(T, (S.re, S.im)))
_structarray(S::AbstractArray, ::Type{T}) where {T} = T(no_offset_view(S))

function divG_radial!(divG::AbstractVector, ℓ::Integer, G::AbstractMatrix, drG::AbstractMatrix)
	# components in PB VSH basis
	pre = -2Ω(ℓ, 0)

	T = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	GP = T(parent(G))
	drGP = T(parent(drG))

	@turbo for r_ind in eachindex(divG)
		divG[r_ind] = pre * GP[2, r_ind]/r[r_ind] + drGP[1, r_ind] + 2/r[r_ind]*GP[1, r_ind]
	end

	divG
end

function divG_radial!(divG::AbstractMatrix, ℓ::Integer, G::AbstractArray{<:Any, 3}, drG::AbstractArray{<:Any, 3})
	# components in PB VSH basis
	pre = -2Ω(ℓ, 0)

	Tdiv = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	divGP = Tdiv(parent(divG))

	TG = HybridArray{Tuple{2,2,StaticArrays.Dynamic()}}
	GP = TG(parent(G))
	drGP = TG(parent(drG))

	@turbo for r_ind in UnitRange(axes(divG, 2))
		invr = 1/r[r_ind]
		for β in UnitRange(axes(GP, 2))
			divGP[β, r_ind] = pre * GP[2, β, r_ind]*invr + drGP[1, β, r_ind] + 2invr*GP[1, β, r_ind]
		end
	end

	divG
end

# Radial components of dG for sound-speed perturbations
function radial_fn_δc_firstborn!(f::StructArray{<:Complex,1},
	divGsrc::AbstractVector{<:Complex}, divGobs::AbstractVector{<:Complex})

	@turbo for I in eachindex(f)
		pre = -ρ[I] * 2c[I]
		f.re[I] = pre * _mulR(divGobs[I], divGsrc[I])
		f.im[I] = pre * _mulI(divGobs[I], divGsrc[I])
	end
	return f
end

function radial_fn_δc_firstborn!(f::AbstractMatrix{<:Complex},
	divGsrc::AbstractVector{<:Complex}, divGobs::AbstractMatrix{<:Complex})

	HT = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	divGobsH = HT(parent(divGobs))

	fR = HT(parent(f.re))
	fI = HT(parent(f.im))

	@turbo for r_ind in eachindex(r)
		twoρcdivGsrc = -ρ[r_ind] * 2c[r_ind] * divGsrc[r_ind]
		for obsind in axes(fR,1)
			fR[obsind, r_ind] = _mulR(divGobsH[obsind, r_ind], twoρcdivGsrc)
			fI[obsind, r_ind] = _mulI(divGobsH[obsind, r_ind], twoρcdivGsrc)
		end
	end
	return f
end

function radial_fn_isotropic_δc_firstborn!(f, Gsrc::AA, drGsrc::AA, divGsrc,
	Gobs::AA, drGobs::AA, divGobs, ℓ) where {AA<:AbstractArray{ComplexF64}}
	radial_fn_δc_firstborn!(f, Gsrc, drGsrc, ℓ, divGsrc, Gobs, drGobs, ℓ, divGobs)
end

# Only radial component
function radial_fn_δc_firstborn!(f::AbstractVector{<:Complex},
	Gsrc::AA, drGsrc::AA, ℓ::Integer, divGsrc::BB,
	Gobs::AA, drGobs::AA, ℓ′::Integer, divGobs::BB) where
		{AA<:AbstractArray{<:Complex, 2}, BB<:AbstractArray{<:Complex, 1}}

	# G is a tensor with one trailing vector index, the first axis is r
	# In this case the divG arrays are 1D (r)
	divG_radial!(divGobs, ℓ′, Gobs, drGobs)
	divG_radial!(divGsrc, ℓ, Gsrc, drGsrc)

	radial_fn_δc_firstborn!(f, divGsrc, divGobs)
end

# All components
function radial_fn_δc_firstborn!(f::AbstractMatrix{<:Complex},
	Gsrc::AA, drGsrc::AA, ℓ::Integer, divGsrc::BB,
	Gobs::AA, drGobs::AA, ℓ′::Integer, divGobs::BB) where
		{AA<:AbstractArray{<:Complex, 3}, BB<:AbstractArray{<:Complex, 2}}

	# G is a tensor with two vector indices
	# In this case the divG arrays are 2D (r, vec_ind)
	divG_radial!(divGobs, ℓ′, Gobs, drGobs)
	divG_radial!(divGsrc, ℓ, Gsrc, drGsrc)

	divGsrc_0 = view(divGsrc, 0, :) # Source is radial

	radial_fn_δc_firstborn!(f, divGsrc_0, divGobs)
end

# All components
# H_βαjₒjₛ(r; r₁, r₂, rₛ) = conj(f_αjₒjₛ(r, r₁, rₛ)) Gβ0jₛ(r₂, rₛ)
function Hjₒjₛω!(Hjₒjₛω_r₁r₂::AbstractArray{<:Any, 3},
	fjₒjₛ_r₁_rsrc::AbstractMatrix, Grjₛ_r₂_rsrc::AbstractVector)

	HT = HybridArray{Tuple{2,2,StaticArrays.Dynamic()}}
	HP = _structarrayparent(Hjₒjₛω_r₁r₂, HT)

	FT = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	FP = _structarrayparent(fjₒjₛ_r₁_rsrc, FT)

	Grjₛ_r₂_rsrcP = no_offset_view(Grjₛ_r₂_rsrc)

	@turbo for r_ind in axes(HP, 3), γ in axes(HP, 2)
		Gγrjₛ_r₂_rsrc = Grjₛ_r₂_rsrcP[γ]
		for α in axes(HP, 1)
			HP[α, γ, r_ind] = conj(FP[α, r_ind]) * Gγrjₛ_r₂_rsrc
		end
	end
	return Hjₒjₛω_r₁r₂
end

# Only radial component
# H_00jₒjₛ(r; r₁, r₂, rₛ) = conj(f_0jₒjₛ(r, r₁, rₛ)) G00jₛ(r₂, rₛ)
function Hjₒjₛω!(Hjₒjₛω_r₁r₂::AbstractVector, fjₒjₛ_r₁_rsrc::AbstractVector, Grrjₛ_r₂_rsrc)
	@turbo for I in CartesianIndices(Hjₒjₛω_r₁r₂)
		Hjₒjₛω_r₁r₂[I] = conj(fjₒjₛ_r₁_rsrc[I]) * Grrjₛ_r₂_rsrc
	end
	return Hjₒjₛω_r₁r₂
end

# Radial components of dG for flows
# Only tangential (+) component
function radial_fn_uniform_rotation_firstborn!(G::AbstractVector,
	Gsrc::AA, Gobs::AA, j, ::los_radial) where {AA<:AbstractMatrix}

	for r_ind in UnitRange(axes(Gsrc,2))
		G[r_ind] = Gsrc[0, r_ind] * Gobs[0, r_ind] -
					Gsrc[0, r_ind] * Gobs[1, r_ind]/Ω(j, 0) -
					Gsrc[1, r_ind]/Ω(j, 0) * Gobs[0, r_ind] +
					ζjₒjₛs(j, j, 1) * Gsrc[1, r_ind] * Gobs[1, r_ind]
	end
	return G
end

function radial_fn_uniform_rotation_firstborn!(G::AbstractMatrix,
	Gsrc::AA, Gobs::AA, j, ::los_earth) where {AA<:AbstractArray{<:Any, 3}}

	for r_ind in UnitRange(axes(Gsrc,3)), α₂ in UnitRange(axes(Gobs, 2))
		G[α₂, r_ind] = Gsrc[0, 0, r_ind] * Gobs[0, α₂, r_ind] -
					Gsrc[0, 0, r_ind] * Gobs[1, α₂, r_ind]/Ω(j, 0) -
					Gsrc[1, 0, r_ind]/Ω(j, 0) * Gobs[0, α₂, r_ind] +
					ζjₒjₛs(j, j, 1) * Gsrc[1, 0, r_ind] * Gobs[1, α₂, r_ind]
	end
	return G
end

function Gⱼₒⱼₛω_u⁺_firstborn!(Gparts⁺, Gsrc, jₛ, Gobs, jₒ, ::los_radial)

	# The actual G¹ₗⱼ₁ⱼ₂ω_00(r, rᵢ, rₛ) =  G[:, 0] + ζ(jₛ, jₒ, ℓ)G[:, 1]
	# We store the terms separately and add them up for each ℓ

	GT = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	GsrcH = GT(parent(Gsrc))
	GobsH = GT(parent(Gobs))

	GH = _structarrayparent(Gparts⁺, GT)

	Ωjₒ0 = Ω(jₒ, 0)
	Ωjₛ0 = Ω(jₛ, 0)

	@turbo for r_ind in axes(GH, 2)
		GH[1, r_ind] = GsrcH[1, r_ind] * GobsH[1, r_ind] -
						GsrcH[1, r_ind] * GobsH[2, r_ind]/Ωjₒ0 -
						GsrcH[2, r_ind] * GobsH[1, r_ind]/Ωjₛ0

		GH[2, r_ind] = GsrcH[2, r_ind] * GobsH[2, r_ind]
	end
	return Gparts⁺
end

function Gⱼₒⱼₛω_u⁺_firstborn!(Gparts⁺, Gsrc, jₛ, Gobs, jₒ, ::los_earth)

	# The actual G¹ₗⱼ₁ⱼ₂ω_α₁0(r, rᵢ, rₛ) = G[:, 0, α₁] + ζ(jₛ, jₒ, ℓ)G[:, 1, α₁]
	# We store the terms separately and add them up for each ℓ

	GT = HybridArray{Tuple{2,2,StaticArrays.Dynamic()}}
	GsrcH = GT(parent(Gsrc))
	GobsH = GT(parent(Gobs))

	GH = _structarrayparent(Gparts⁺, GT)

	Ωjₒ0 = Ω(jₒ, 0)
	Ωjₛ0 = Ω(jₛ, 0)

	@turbo for r_ind in axes(GH, 3), α₁ in axes(GH, 2)
		GH[1, α₁, r_ind] = GsrcH[1, 1, r_ind] * GobsH[1, α₁, r_ind] -
							GsrcH[1, 1, r_ind] * GobsH[2, α₁, r_ind]/Ωjₒ0 -
							GsrcH[2, 1, r_ind] * GobsH[1, α₁, r_ind]/Ωjₛ0

		GH[2, α₁, r_ind] = GsrcH[2, 1, r_ind] * GobsH[2, α₁, r_ind]
	end
	return Gparts⁺
end

function Gⱼₒⱼₛω_u⁰_firstborn!(Gparts⁰, drGsrc, jₛ, Gobs, jₒ, ::los_radial)

	# The term corresponding to u⁰ is G⁰₀ ∂ᵣG⁰₀ + ζ(jₛ, jₒ, ℓ) G¹₀ ∂ᵣG¹₀
	# We express this as G⁰ₗⱼ₁ⱼ₂ω_00(r, rᵢ, rₛ) =  G[:, 0] + ζ(jₛ, jₒ, ℓ)G[:, 1]
	# We store the terms separately and add them up for each ℓ

	GT = HybridArray{Tuple{2,StaticArrays.Dynamic()}}
	drGsrcH = GT(parent(drGsrc))
	GobsH = GT(parent(Gobs))

	GH = _structarrayparent(Gparts⁰, GT)

	@turbo for I in CartesianIndices(GH)
		GH[I] = GobsH[I] * drGsrcH[I]
	end
	return Gparts⁰
end

function Gⱼₒⱼₛω_u⁰_firstborn!(Gparts⁰, drGsrc, jₛ, Gobs, jₒ, ::los_earth)

	# The term corresponding to u⁰ is G⁰α₁ ∂ᵣG⁰₀ + ζ(jₛ, jₒ, ℓ) G¹α₁ ∂ᵣG¹₀
	# We express this as G⁰ₗⱼ₁ⱼ₂ω_α₁0(r, rᵢ, rₛ) = G[0, α₁, :] + ζ(jₛ, jₒ, ℓ)G[1, α₁, :]
	# We store the terms separately and add them up for each ℓ

	GT = HybridArray{Tuple{2,2,StaticArrays.Dynamic()}}
	GobsH = GT(parent(Gobs))
	drGsrcH = GT(parent(drGsrc))

	GH = _structarrayparent(Gparts⁰, GT)

	@turbo for r_ind in axes(GH, 3), α₁ in axes(GH, 2), ind in axes(GH, 1)
		GH[ind, α₁, r_ind] = GobsH[ind, α₁, r_ind] * drGsrcH[ind, 1, r_ind]
	end
	return Gparts⁰
end

# This function computes Gparts
# Components (0) and (+)
function Gⱼₒⱼₛω_u⁰⁺_firstborn!(Gparts, Gparts2, Gsrc, drGsrc, jₛ, Gobs, jₒ, los::los_direction)
	Gⱼₒⱼₛω_u⁰_firstborn!(Gparts[0], drGsrc, jₛ, Gobs, jₒ, los)
	Gⱼₒⱼₛω_u⁺_firstborn!(Gparts[1],  Gsrc, jₛ, Gobs, jₒ, los)
	copyto!(Gparts2, 1, Gparts[0], 1, length(Gparts[0]))
	copyto!(Gparts2, length(Gparts[0]) + 1, Gparts[1], 1, length(Gparts[1]))
	return Gparts2
end

# We evaluate Gγₗⱼ₁ⱼ₂ω_00(r, rᵢ, rₛ) = Gsum[γ, r] = Gparts[0, :] + ζ(jₛ, jₒ, ℓ)G[1, :]
# This function is run once for each ℓ using cached values of G
function Gγₗⱼₒⱼₛω_α₁_firstborn!(Gγₗⱼₒⱼₛω_α₁::AbstractMatrix, Gparts2, jₒ, jₛ, ℓ)

	coeff = ζjₒjₛs( jₒ, jₛ, ℓ)

	GT = HybridArray{Tuple{StaticArrays.Dynamic(),2}}
	G = _structarrayparent(Gγₗⱼₒⱼₛω_α₁, GT)

	GPT = HybridArray{Tuple{2,StaticArrays.Dynamic(),2}}
	GP = _structarrayparent(Gparts2, GPT)

	@turbo for I in CartesianIndices(G)
		G[I] = GP[1, I] + coeff * GP[2, I]
	end
	return Gγₗⱼₒⱼₛω_α₁
end

# We evaluate Gγₗⱼ₁ⱼ₂ω_α₁0(r, rᵢ, rₛ) = Gsum[r, γ, α₁] = G[:, 0, α₁] + ζ(jₛ, jₒ, ℓ)G[:, 1, α₁]
# This function is run once for each ℓ using cached values of G
function Gγₗⱼₒⱼₛω_α₁_firstborn!(Gγₗⱼₒⱼₛω_α₁::StructArray{<:Complex, 3}, Gparts2, jₒ, jₛ, ℓ)

	coeff = ζjₒjₛs( jₒ, jₛ, ℓ)

	GT = HybridArray{Tuple{2,StaticArrays.Dynamic(),2}}
	G = _structarrayparent(Gγₗⱼₒⱼₛω_α₁, GT)

	GPT = HybridArray{Tuple{2,2,StaticArrays.Dynamic(),2}}
	GP = _structarrayparent(Gparts2, GPT)

	@turbo for I in CartesianIndices(G)
		G[I] = GP[1, I] + coeff * GP[2, I]
	end
	return Gγₗⱼₒⱼₛω_α₁
end

# We compute Hγₗⱼ₁ⱼ₂ω_00(r, r₁, r₂, rₛ) = conj(Gγₗⱼ₁ⱼ₂ω_00(r, r₁, rₛ)) * Gⱼ₂ω_00(r₂, rₛ)
function Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(H::StructArray{<:Complex, 2}, Gparts,
	Gγₗⱼₒⱼₛω_α₁::AbstractMatrix{<:Complex}, Grr::Complex, jₒ, jₛ, ℓ)

	HT = HybridArray{Tuple{StaticArrays.Dynamic(),2}}
	HR = HT(parent(H.re))
	HI = HT(parent(H.im))

	GT = HybridArray{Tuple{StaticArrays.Dynamic(),2}}
	GR = GT(parent(Gγₗⱼₒⱼₛω_α₁.re))
	GI = GT(parent(Gγₗⱼₒⱼₛω_α₁.im))

	Gγₗⱼₒⱼₛω_α₁_firstborn!(Gγₗⱼₒⱼₛω_α₁, Gparts, jₒ, jₛ, ℓ)
	# @. H = conj(Gγₗⱼₒⱼₛω_α₁) * Grr
	@turbo for I in CartesianIndices(HR)
		HR[I] = GR[I] * real(Grr) + GI[I] * imag(Grr)
		HI[I] = GR[I] * imag(Grr) - GI[I] * real(Grr)
	end
	return H
end

# We compute Hγₗⱼ₁ⱼ₂ω_α₁α₂(r, r₁, r₂, rₛ) = conj(Gγₗⱼ₁ⱼ₂ω_α₁0(r, r₁, rₛ)) * Gⱼ₂ω_α₂0(r₂, rₛ)
function Hγₗⱼₒⱼₛω_α₁α₂_firstborn!(H::StructArray{<:Complex, 4}, Gparts,
	Gγₗⱼₒⱼₛω_α₁::StructArray{<:Complex, 3}, G::AbstractVector{<:Complex}, jₒ, jₛ, ℓ)

	Gγₗⱼₒⱼₛω_α₁_firstborn!(Gγₗⱼₒⱼₛω_α₁, Gparts, jₒ, jₛ, ℓ)

	HT = HybridArray{Tuple{2,2,StaticArrays.Dynamic(),2}}
	HS = _structarrayparent(H, HT)

	GT = HybridArray{Tuple{2,StaticArrays.Dynamic(),2}}
	Gγₗⱼₒⱼₛω_α₁sH = _structarrayparent(Gγₗⱼₒⱼₛω_α₁, GT)

	Gp = parent(G)

	@turbo for γ in axes(HS, 4), rind in axes(HS, 3), obsind2 in axes(HS, 2)
		Gsrc_obs2 = Gp[obsind2]
		for obsind1 in axes(HS, 1)
			HS[obsind1, obsind2, rind, γ] = conj(Gγₗⱼₒⱼₛω_α₁sH[obsind1, rind, γ]) * Gsrc_obs2
		end
	end
	return H
end


text:
@theme VimTheme Dict(
    :name => "Vim",
    :description => "A theme similar to the vim 7.0 default theme.",
    :comments => "Based on the theme from Pygments.",
    :style => S"bg: 000000; fg: cccccc",
    :tokens => Dict(
        COMMENT             => S"fg: 000080",
        COMMENT_PREPROC     => S"",
        COMMENT_SPECIAL     => S"bold; fg: cd0000",

        ERROR               => S"bg: ff0000",

        GENERIC_DELETED     => S"fg: cd0000",
        GENERIC_EMPH        => S"italic",
        GENERIC_ERROR       => S"fg: ff0000",
        GENERIC_HEADING     => S"bold; fg: 000080",
        GENERIC_INSERTED    => S"fg: 00cd00",
        GENERIC_OUTPUT      => S"fg: 888",
        GENERIC_PROMPT      => S"bold; fg: 000080",
        GENERIC_STRONG      => S"bold",
        GENERIC_SUBHEADING  => S"bold; fg: 800080",
        GENERIC_TRACEBACK   => S"fg: 04d",

        KEYWORD             => S"fg: cdcd00",
        KEYWORD_DECLARATION => S"fg: 00cd00",
        KEYWORD_NAMESPACE   => S"fg: cd00cd",
        KEYWORD_PSEUDO      => S"",
        KEYWORD_TYPE        => S"fg: 00cd00",

        NAME                => S"",
        NAME_BUILTIN        => S"fg: cd00cd",
        NAME_CLASS          => S"fg: 00cdcd",
        NAME_EXCEPTION      => S"bold; fg: 666699",
        NAME_VARIABLE       => S"fg: 00cdcd",

        NUMBER              => S"fg: cd00cd",

        OPERATOR            => S"fg: 3399cc",
        OPERATOR_WORD       => S"fg: cdcd00",

        STRING              => S"fg: cd0000",

        TEXT                => S"",

        WHITESPACE          => S"",
    )
)

ext: jl


text:
using Plots.Dates
using StatsPlots

json_files_names = [
    ["../../../nim/godot_test/_godotapi/api.json","godot_api.json"],
    ["../../../nim/godot-nim-stub/_godotapi/api.json","godot_api2.json"],
    ["../../../impala/anydsl/llvm_build/compile_commands.json","llvm_compile_commands.json"],
    ["../../../web/vuegl-pool/package-lock.json","vuegl_pkg_lock.json"],
    ["../../../c/simdjson/jsonexamples/twitter.json","simdjson_twitter.json"],
]
json_files = map(first, json_files_names)
json_names = map(last, json_files_names)

parsers_names = [
    ["bin/zig_json", "zig std lib"],
    ["../../JSONTestSuite/simdjson/simdjson", "c++ simdjson"],
    ["../zig-out/bin/simdjzon", "zig simdjzon"],
    ["../../JSONTestSuite/parsers/test_nim/test_json", "nim json"],
    ["../../JSONTestSuite/parsers/test_go/test_json", "go json"],
]
parsers = map(first, parsers_names)
parser_names = map(last, parsers_names)
simdjzon_parser_idx = findfirst(x -> x == "zig simdjzon", parser_names)

# begin
fileslen = length(json_files)
parserslen = length(parsers)
bench_results = Dict(parser => [[] for _ in 1:fileslen] for parser in parsers)

function bench(parser, json_file, i, j; dry_run=false, debug=false)
    out = Pipe()
    proc = run(pipeline(`./time.sh $parser $json_file`, stdout=out), wait=true)
    close(out.in)
    output = String(read(out))

    debug && println("$(parser) $json_file exitcode $(proc.exitcode) stdout $output")
    if !dry_run
        @assert(proc.exitcode == 0, "parse validation failure")
        push!(bench_results[parser][i], parse(Float16, output))
    end
end

# dry run to warm up cache
for parser in parsers, json_file in json_files
    bench(parser, json_file, 1, 1; dry_run=true)
end
num_runs = 10
# gather data
for j in 1:num_runs, (i, json_file) in enumerate(json_files)
    for parser in parsers
        bench(parser, json_file, i, j; debug=false)
    end
end

# average the results and
# append simdjzon results to csv file
_mn = [[] for _ in parsers]
parser_idx = 1
open("benchmark_results.csv"; append=true) do results_file
    write(results_file, join(["Date", json_names...], ", "), "\n")
    now = Dates.now()
    write(results_file, Dates.format(now, "Y-m-d"), ", ")
    for parser in parsers
        parser_results = bench_results[parser]
        for result_set in parser_results
            avg = sum(result_set) / length(result_set)
            if parser_idx == simdjzon_parser_idx
                write(results_file, "$avg, ")
            end
            push!(_mn[parser_idx], avg)
        end
        global parser_idx += 1
    end
    write(results_file, "\n")
end

# prepare for plotting
lenp = length(parsers)
lenf = length(json_files)
json_file_sizes_kb = map(x -> filesize(x) ÷ 1000, json_files)
json_names_lens = map(x -> json_names[x] * " ($(json_file_sizes_kb[x])Kb)", 1:lenf)
nam = repeat(json_names_lens, inner=lenp)
sx = repeat(parser_names, outer=lenf)

# reshape, flipping rows/cols
dim1 = length(_mn)
dim2 = length(_mn[1])
matrix = Array{eltype(_mn[1]),2}(undef, dim1, dim2)
for i in 1:dim1, j in 1:dim2
    matrix[i,j] = _mn[i][j]
end
mn = collect(Iterators.flatten(reshape(matrix, dim2, dim1)))

# plot
l = @layout [a{0.01h}; grid(1, 1)]
p = fill(plot(), 2, 1)
p[1] = plot(title="JSON Validation Times", framestyle=nothing, showaxis=false, xticks=false, yticks=false, margin=0Plots.mm)
p[2] = groupedbar(nam, mn, group=sx, ylabel="milliseconds",
            title="average of $num_runs runs", xrotation=-7, titlefontsize=8,
            margin=0Plots.mm, framestyle=nothing, margin_right=10Plots.mm
    )
plot(p..., layout=l)
savefig("validation_grouped.png")

# old line plot
if false
    using Plots
    x = 1:fileslen
    p = plot(x, _mn[1], title="validation time",
                label=parser_names[1], lw=3,
                ylabel="milliseconds", yrotation=0,
                xlabel="file", xticks=(x, json_names_lens), xrotation=-7
        )
    for i in 2:length(parsers)
        line = _mn[i]
        plot!(p, x, line, label=parser_names[i], lw=3)
    end
    savefig("validation_line.png")
end

max_stars_repo_name: travisstaloch/simdjzon

max_issues_repo_name: travisstaloch/simdjzon

max_forks_repo_name: travisstaloch/simdjzon


text:
function add_def(s::Solver, x::Z3Expr)
    name = bool_const(ctx(s), "$(id(x))")
    add(s, name == x)
    name
end

function relax_core(s::Solver, core::ExprVector, cs::Vector{<:Z3Expr})
    prefix = bool_val(ctx(s), true)
    for c in core
        deleteat!(cs, findfirst(isequal(c), cs))
    end
    for i in 1:length(core)-1
        prefix = add_def(s, and(core[i], prefix))
        push!(cs, add_def(s, or(prefix, core[i+1])))
    end
    cs
end

function maxsat(s::Solver, cs::Vector{<:Z3Expr})
    cost = 0
    cs0 = copy(cs)
    cs1 = copy(cs)
    push(s)
    while (res = check(s, cs1)) == Z3.unsat && cost < length(cs0)
        cost += 1
        cs1 = relax_core(s, unsat_core(s), cs1)
    end
    @debug "Maxsat" cost res==Z3.sat reason_unknown(s)
    res != Z3.sat && return res, -1, nothing
    satisfied = [c for c in cs0 if is_true(Z3.eval(get_model(s), c, false))]
    m = get_model(s)
    pop(s, 1)
    @debug "Maxsat" satisfied m
    res, cost, m
end


text:
using ArgParse, Distributed

function parse_cli(args)
    settings = ArgParseSettings("ABC-MK CLI")

    @add_arg_table! settings begin
        # CMD 1
        "rates", "R"
            help = "Function to solve fixation and polymorphic rates analitically. The function will create N random models from prior values. Use the arguments to defined the input range for each parameter.\nIf rho and/or theta are set to nothing, the function will input random values given the range 0.0005:0.0005:0.01. Otherwise you can fix the values.\nIf gL is set to nothing, the function will not account the role of the weakly selected alleles in the estimation.\nThe function returns a HDF5 file containing models solved and rates. The rates will be used to compute summary statistics required at ABC.\nPlease check the documentation to get more info about models parameters or detailed arguments description https://jmurga.github.io/Analytical.jl/dev/cli/ to check model
            "
            action = :command
        # CMD 2
        "parse_data", "P"
            help = "Function to parse polymorphic and divergence data from Uricchio et. al (2019) and Murga-Moreno et al (2019). Please input a path to create a new analysis folder. You can filter the dataset using a file containing a list of Ensembl IDs. The function returns files containing raw polymorphic and divergence data, parsed SFS and parsed divegence required to estimate summary statistics. Please check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli/"
            action = :command
        # CMD 3
        "summaries", "S"
            help = "Function to parse polymorphic and divergence data from Uricchio et. al (2019) and Murga-Moreno et al (2019). Please input a path to create a new analysis folder. You can filter the dataset using a file containing a list of Ensembl IDs. The function returns files containing raw polymorphic and divergence data, parsed SFS and parsed divegence required to estimate summary statistics. Please check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli/"
            action = :command
        # CMD 4
        "inference", "I"
            help = "ABCreg inference. The function returns posterior distributions from ABC inference. Each posterior file contains information about alpha_w, alpha_s, alpha, gamNeg and shape parameter. The number of posterior distributions will depend on the number of bootstrap replicas. Check the documentation to get more info https://jmurga.github.io/Analytical.jl/dev/cli"
            action = :command

    end

    add_arg_table!(settings["rates"],
        ["--pop_size"],
        Dict(
            :help => "Population size",
            :arg_type => Int64,
            :default => 1000
        ),
        ["--sample_size"],
        Dict(
            :help => "Sample size",
            :arg_type => Int64,
            :required => true
        ),
        ["--dac"],
        Dict(
            :help => "Derived Allele Count",
            :arg_type => String,
            :required => true
        ),
        ["--gam_neg"],
        Dict(
            :help => "Selection coefficient for deleterious alleles",
            :arg_type => String,
            :required => true
        ),
        ["--positive_strong"],
        Dict(
            :help => "Selection coefficient for strongly beneficial alleles",
            :arg_type => String,
            :required => true
        ),
        ["--positive_weak"],
        Dict(
            :help => "Selection coefficient for weakly beneficial alleles",
            :arg_type => String,
            :required => true
        ),
        ["--shape"],
        Dict(
            :help => "Shape value modeling Gamma distribution for deleterious alleles",
            :arg_type => Float64,
            :default => 0.184
        ),
        ["--rho"],
        Dict(
            :help => "Recombination rate",
            :arg_type => Float64,
            :default => 0.001
        ),
        ["--theta"],
        Dict(
            :help => "Mutation rate on coding locus",
            :arg_type => Float64,
            :default => 0.001
        ),
        ["--solutions"],
        Dict(
            :help => "Mutation rate on coding locus",
            :arg_type => Int64,
            :default => 100_000
        ),
        ["--output"],
        Dict(
            :help => "Output file",
            :arg_type => String,
            :default => "rates.jld2"
        ),
        ["--scheduler"],
        Dict(
            :help => "Select scheduler manager",
            :arg_type => String,
            :default => "local"
        ),
        ["--nthreads"],
        Dict(
            :help => "Select number of threads to parallelize",
            :arg_type => Int64,
            :default => 1
        ),
    )

    add_arg_table!(settings["parse_data"],
        "folder",
        Dict(
            :help => "a positional argument",
            :required => true,
            :arg_type => String
        ),
        ["--dataset","-d"],
        Dict(
            :help => "a positional argument",
            :default => "tgp",
            :arg_type => String
        ),
        ["--gene_list","-g"],
        Dict(
            :help => "a positional argument",
            :arg_type => Union{Bool,String},
            :default => false
        ),
        ["--bins","-b"],
        Dict(
            :help => "a positional argument",
            :arg_type => Union{Bool,Int64},
            :default => false
        )
    )

    add_arg_table!(settings["summaries"],
        "folder",
        Dict(
            :help => "Folder path containing SFS and divergence files to run the analysis",
            :required => true,
            :arg_type => String
        ),
        ["--rates"],
        Dict(
            :help => "H5 file containing precomputed rates",
            :required => true,
            :arg_type => String
        ),
        ["--sample_size"],
        Dict(
            :help => "Sample size",
            :required => true,
            :arg_type => Int64
        ),
        ["--dac"],
        Dict(
            :help => "Derived allele count",
            :required => true,
            :arg_type => String
        ),
        ["--summstat_size"],
        Dict(
            :help => "Define number of summary estatistic to perform ABC",
            :required => true,
            :arg_type => Int64
        ),
        ["--bootstrap"],
        Dict(
            :help => "Allow bootstrap following polyDFE manual",
            :arg_type => Bool,
            :default => false
        ),
        ["--replicas"],
        Dict(
            :help => "Number of bootstrap replicas",
            :arg_type => Int64,
            :default => 1
        ),
        ["--nthreads"],
        Dict(
            :help => "Select scheduler manager",
            :arg_type => Int64,
            :default => 1
        ),
        ["--scheduler"],
        Dict(
            :help => "Select scheduler manager",
            :arg_type => String,
            :default => "local"
        )
    )
    add_arg_table!(settings["inference"],
        "folder",
        Dict(
            :help => "Folder path containing SFS and divergence files to run the analysis",
            :required => true,
            :arg_type => String
        ),
        ["--S"],
        Dict(
            :help => "Define number of summary estatistic to perform ABC",
            :required => true,
            :arg_type => Int64
        ),
        ["--tol"],
        Dict(
            :help => "Tolerance",
            :required => true,
            :arg_type => Float64
        ),
        ["--abcreg"],
        Dict(
            :help => "ABCreg static binary",
            :required => true,
            :arg_type => String
        ),
        ["--nthreads"],
        Dict(
            :help => "Select scheduler manager",
            :arg_type => Int64,
            :default => 1
        ),
        ["--scheduler"],
        Dict(
            :help => "Select scheduler manager",
            :arg_type => String,
            :default => "local"
        )
    )

    return parse_args(settings)
end

cli = parse_cli(ARGS)

for cmd in keys(cli)
    if(cmd == "parse_data")

        @eval using Analytical, DataFrames, CSV

        folder = cli[cmd]["folder"]

        run(`mkdir -p $folder`)

        dataset = lowercase(cli[cmd]["dataset"])
        data    = folder * "/" * dataset * ".txt"

        @eval download("https://raw.githubusercontent.com/jmurga/Analytical.jl/master/data/"* $dataset * ".txt",$data)

        # Check if bins or gene_list are defined
        gene_list = cli[cmd]["gene_list"]
        @eval if $gene_list != false
            @eval gList = CSV.read($gene_list,DataFrame,header=false) |> Array
        else
            @eval gList = nothing
        end

        bins = cli[cmd]["bins"]

        @eval if $bins != 0
            @eval bins_size = $bins
        else
            bins_size = nothing
        end

        # Parsing TGP data
        if dataset == "tgp"
            @eval α,sfs, divergence = Analytical.parse_sfs(sample_size=661,data=$data,gene_list=$gList,bins=$bins_size)
        elseif occursin("zi",dataset)
            @eval α,sfs, divergence = Analytical.parse_sfs(sample_size=154,data=$data,gene_list=$gList,bins=$bins_size,isolines=true)
        elseif occursin("ral",dataset)
            @eval α,sfs, divergence = Analytical.parse_sfs(sample_size=160,data=$data,gene_list=$gList,bins=$bins_size,isolines=true)
        end
        # Writting data to folder
        @eval sName = $folder * "/sfs.tsv"
        @eval dName = $folder * "/div.tsv"

        @eval CSV.write($sName,DataFrame($sfs,:auto),delim='\t',header=false)
        @eval CSV.write($dName,DataFrame($divergence',:auto),delim='\t',header=false)
    elseif (cmd == "rates")

        neg    = parse.(Int,split(cli[cmd]["gam_neg"],":"))
        strong = parse.(Int,split(cli[cmd]["positive_strong"],":"))

        if (cli[cmd]["positive_weak"] == false)
            weak = nothing
        else
            weak = parse.(Int,split(cli[cmd]["positive_weak"],":"))
        end

        dac = parse.(Int,split(cli[cmd]["dac"],","))

        if (cli[cmd]["rho"] == false)
            rho = nothing
        else
            rho = cli[cmd]["rho"]
        end

        if (cli[cmd]["theta"] == "nothing")
            theta = nothing
        else
            theta = cli[cmd]["theta"]
        end

        scheduler = cli[cmd]["scheduler"];nthreads = cli[cmd]["nthreads"]
        if scheduler == "slurm"
            @eval using ClusterManagers
            @eval addprocs_slurm($nthreads)
        elseif scheduler == "htcondor"
            @eval using ClusterManagers
            @eval addprocs_htc($nthreads)
        else
            @eval addprocs($nthreads)
        end

        ne = cli[cmd]["pop_size"];samples= cli[cmd]["sample_size"];shape = cli[cmd]["shape"]
        @eval @everywhere using Analytical, ParallelUtilities
        @eval adap = Analytical.parameters(N=$ne,n=$samples,dac=$dac,al=$shape)

        @eval cnv = Analytical.binomial_dict()
        @eval Analytical.binomOp!($adap,$cnv.bn);

        solutions = cli[cmd]["solutions"]; output = cli[cmd]["output"]
        @eval Analytical.rates(param = $adap,convoluted_samples=$cnv,gH=$strong[1]:$strong[2],gL=$weak[1]:$weak[2],gamNeg=$neg[1]:$neg[2],iterations = $solutions,rho=$rho,theta=$theta,shape=$adap.al,output=$output,scheduler=$scheduler);

        for i in workers()
            rmprocs(i)
        end
    elseif (cmd == "summaries")

        scheduler = cli[cmd]["scheduler"];nthreads = cli[cmd]["nthreads"];

        if scheduler == "slurm"
            @eval using ClusterManagers
            @eval addprocs_slurm($nthreads)
        elseif scheduler == "htcondor"
            @eval using ClusterManagers
            @eval addprocs_htc($nthreads)
        else
            @eval addprocs($nthreads)
        end

        samples       = cli[cmd]["sample_size"]
        dac           = parse.(Int,split(cli[cmd]["dac"],","))
        bootstrap     = cli[cmd]["bootstrap"]
        replicas      = cli[cmd]["replicas"]
        rates         = cli[cmd]["rates"]
        summstat_size = cli[cmd]["summstat_size"]
        folder        = cli[cmd]["folder"]

        @eval  using JLD2, DataFrames, CSV, ProgressMeter
        @eval @everywhere using  Analytical, ParallelUtilities
        @eval adap    = Analytical.parameters(n=$samples,dac = $dac)

        @eval if ($bootstrap == true)
            @eval summstat  = Analytical.summary_statistics(param=$adap,h5_file=$rates,analysis_folder=$folder,summstat_size=$summstat_size,replicas=$replicas,bootstrap=true)
        else
            @eval summstat  = Analytical.summary_statistics(param=$adap,h5_file=$rates,analysis_folder=$folder,summstat_size=$summstat_size,replicas=$replicas,bootstrap=false)
        end

        for i in workers()
            rmprocs(i)
        end
    elseif (cmd == "inference")

        scheduler = cli[cmd]["scheduler"]
        nthreads = cli[cmd]["nthreads"]
        folder   = cli[cmd]["folder"]
        S        = cli[cmd]["S"]
        tol      = cli[cmd]["tol"]
        abcreg   = cli[cmd]["abcreg"]

        if scheduler == "slurm"
            @eval using ClusterManagers
            @eval addprocs_slurm($nthreads)
        elseif scheduler == "htcondor"
            @eval using ClusterManagers
            @eval addprocs_htc($nthreads)
        else
            @eval addprocs($nthreads)
        end

        @eval @everywhere using Analytical,ParallelUtilities

        @eval Analytical.ABCreg(analysis_folder=$folder,S=$S,tol=$tol,abcreg=$abcreg)

        for i in workers()
            rmprocs(i)
        end
    end
end


text:
""
function build_QCModel_result(qcm::QuantumCircuitModel, solve_time::Number)
    # try-catch is needed until solvers reliably support ResultCount()
    result_count = 1
    try
        result_count = JuMP.result_count(qcm.model)
    catch
        Memento.warn(_LOGGER, "the given optimizer does not provide the ResultCount() attribute, assuming the solver returned a solution which may be incorrect.");
    end

    solution = Dict{String,Any}()

    if result_count > 0
        solution = build_QCModel_solution(qcm)
    else
        Memento.warn(_LOGGER, "Quantum circuit model has no results, solution cannot be built")
    end

    result = Dict{String,Any}(
        "optimizer" => JuMP.solver_name(qcm.model),
        "termination_status" => JuMP.termination_status(qcm.model),
        "primal_status" => JuMP.primal_status(qcm.model),
        "objective" => get_objective_value(qcm.model),
        "objective_lb" => get_objective_bound(qcm.model),
        "solve_time" => solve_time,
        "solution" => solution,
    )

    # Objective slack Penalty
    if qcm.data["decomposition_type"] == "approximate"
        result["objective_slack_penalty"] = qcm.options.objective_slack_penalty
    end

    return result
end

""
function get_objective_value(model::JuMP.Model)
    obj_val = NaN

    try
        obj_val = JuMP.objective_value(model)
    catch
        Memento.warn(_LOGGER, "Objective value is unbounded. Problem may be infeasible or not constrained properly");
    end

    return obj_val
end

""
function get_objective_bound(model::JuMP.Model)
    obj_lb = -Inf

    try
        obj_lb = JuMP.objective_bound(model)
    catch
    end

    return obj_lb
end

function build_QCModel_solution(qcm::QuantumCircuitModel)
    solution = Dict{String,Any}()
    for i in keys(qcm.variables)
        solution[String(i)] = JuMP.value.(qcm.variables[i])
    end
    return solution
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct VariableReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        timeframe::Float64
        requirement::Float64
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure for the procurement products for system simulations.

# Arguments
- `name::String`
- `timeframe::Float64`: the relative saturation timeframe, validation range: (0, nothing), action if invalid: error
- `requirement::Float64`: the required quantity of the product should be scaled by a Forecast
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct VariableReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    "the relative saturation timeframe"
    timeframe::Float64
    "the required quantity of the product should be scaled by a Forecast"
    requirement::Float64
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function VariableReserve{T}(name, timeframe, requirement, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), ) where T <: ReserveDirection
    VariableReserve{T}(name, timeframe, requirement, ext, forecasts, InfrastructureSystemsInternal(), )
end

function VariableReserve{T}(; name, timeframe, requirement, ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), ) where T <: ReserveDirection
    VariableReserve{T}(name, timeframe, requirement, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function VariableReserve{T}(::Nothing) where T <: ReserveDirection
    VariableReserve{T}(;
        name="init",
        timeframe=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::VariableReserve) = value.name
"""Get VariableReserve timeframe."""
get_timeframe(value::VariableReserve) = value.timeframe
"""Get VariableReserve requirement."""
get_requirement(value::VariableReserve) = value.requirement
"""Get VariableReserve ext."""
get_ext(value::VariableReserve) = value.ext

InfrastructureSystems.get_forecasts(value::VariableReserve) = value.forecasts
"""Get VariableReserve internal."""
get_internal(value::VariableReserve) = value.internal


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using LineSearches
using Flux
using StatsBase
using SparseArrays
using Distributed
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs
include("../../../main/compute.jl")

#Bruno-like denoising using wst

##5) Iterative: using empirical noisy dbn (for non-Omega) and true: with pixelwise regularization
ARGS_buffer = ["reg", "nonapd", "noiso", "scratch_NM/StandardizedExp/Nx64/", "full_3losstest", "Full+Eps"]
ENV_buffer= "1000"
numfile = Base.parse(Int, ENV_buffer)
println(numfile, ARGS_buffer[1], ARGS_buffer[2])

if ARGS_buffer[1]=="log"
    logbool=true
else
    if ARGS_buffer[1]!="reg" error("Invalid log arg") end
    logbool=false
end

if ARGS_buffer[2]=="apd"
    apdbool = true
else
    if ARGS_buffer[2]!="nonapd" error("Invalid apd arg") end
    apdbool=false
end

if ARGS_buffer[3]=="iso"
    isobool = true
else
    if ARGS_buffer[3]!="noiso" error("Invalid iso arg") end
    isobool=false
end

direc = ARGS_buffer[4] #"../StandardizedExp/Nx64/noisy_stdtrue/" #Change
datfile = direc * "Data_" * string(numfile) * ".jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]

fname_save = direc * "scratch_NM/NewWrapper/5-30/" * string(numfile) * "_try1"  #Change
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense

optim_settings = Dict([("iterations", 20), ("norm", false), ("minmethod", ConjugateGradient())])
recon_settings = Dict([("log", logbool), ("Invcov_matrix", ARGS_buffer[6]), ("optim_settings", optim_settings), ("eps_value_sfd", 1e-5), ("eps_value_init", 1e-10)]) #Add constraints

recon_settings["datafile"] = datfile

if dhc_args[:iso]
    error("Not constructed for iso")
else #Not iso
    coeffmask = falses(2+Nf+Nf^2)
    coeffmask[Nf+3:end] .= Diagonal(trues(Nf))[:] #Diagonal(triu(trues(Nf, Nf)))[:]
end

#Weight given to terms
lval2 = 0.0
lval3 = 0.0

println("Regularizer Lambda=", round(lval3, sigdigits=3))
#input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov, dhc_args, LossFunc, dLossFunc;
#FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), lambda=0.001, func_specific_params=nothing

#Generate fixed noise maps
Nr=10000
sigma = loaddf["std"]
noisemaps = []
for n=1:Nr
    noisyim = randn((Nx, Nx)).*sigma
    push!(noisemaps, noisyim)
end

s_noisy = DHC_compute_wrapper(loaddf["init"], filter_hash, norm=false; dhc_args...)[coeffmask]
starget = s_noisy #DHC_compute_wrapper(loaddf["true_img"], filter_hash, norm=false; dhc_args...)[coeffmask] #

@time pmap(image -> DHC_compute_wrapper(image, filter_hash; norm=false, dhc_args...)[coeffmask], noisemaps)

function ens_noise()
    empsn = []
    for n=1:Nr
        push!(empsn,  DHC_compute_wrapper(noisemaps[n], filter_hash, norm=false; dhc_args...)[coeffmask])
    end
    return empsn
end

@time ens_noise()

function Loss_noisyfast()
    empsn = []

function Loss_noisy(img_curr, filter_hash, dhc_args; noisemaps=nothing, coeffmask=nothing, starget=nothing)
    #func_specific_params should contain: coeff_mask2, target2, invcov2
    empsn = []
    M = 500
    totdiff = 0
    for n=1:M
        noisy = img_curr + noisemaps[n]
        diffsq = sum((DHC_compute_wrapper(noisy, filter_hash; norm=false, dhc_args...)[coeffmask] .- starget).^2)
        totdiff += diffsq
    end
    totdiff = totdiff/M
    return totdiff
end

function dLoss_noisy!(storage_grad, img_curr, filter_hash, dhc_args; noisemaps=nothing, coeffmask=nothing, starget=nothing, dA=nothing)
    #func_specific_params should contain: coeff_mask2, target2, invcov2
    empsn = []
    M = 500 #length(noisemaps)
    totdiff = zeros(Nx^2, 1)
    for n=1:M
        noisy = img_curr + noisemaps[n]
        diffsq = 2*(DHC_compute_wrapper(noisy, filter_hash; norm=false, dhc_args...)[coeffmask] .- starget)
        wts = zeros(2+Nf+Nf^2)
        wts[coeffmask] .= diffsq
        wts = reshape(wts[Nf+3:end], (Nf, Nf))
        #wts = ReconFuncs.augment_weights_S20(convert(Array{Float64,1}, diffsq), filter_hash, dhc_args, coeffmask)
        derve = Deriv_Utils_New.wst_S20_deriv_sum(noisy, filter_hash, wts)
        totdiff .+= derve
    end
    totdiff ./= M
    storage_grad .= reshape(totdiff, (Nx, Nx))
end

func_specific_params = Dict([(:coeffmask=> coeffmask), (:starget=>starget), (:noisemaps=> noisemaps)])

recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings

if logbool
    error("Not implemented here")
end

res, recon_img = image_recon_derivsum_custom(init, filter_hash, dhc_args, Loss_noisy, dLoss_noisy!; optim_settings=optim_settings, func_specific_params)
heatmap(init)
heatmap(recon_img, title="Using true shift")
heatmap(true_img)

function calc_1dps_local(image, kbins::Array{Float64,1})
    #Assumes uniformly spaced kbins
    Nx = size(image)[1]
    fft_zerocenter = fftshift(fft(image))
    impf = abs2.(fft_zerocenter)
    x = (collect(1:Nx) * ones(Nx)') .- (Nx/2.0)
    y = (ones(Nx) * collect(1:Nx)') .- (Nx/2.0)
    krad  = (x.^2 + y.^2).^0.5
    meanpk = zeros(size(kbins))
    kdel = kbins[2] - kbins[1]
    #println(size(meanpk), " ", kdel)
    for k=1:size(meanpk)[1]
        filt = findall((krad .>= (kbins[k] - kdel./2.0)) .& (krad .<= (kbins[k] + kdel./2.0)))
        meanpk[k] = mean(impf[filt])
    end
    return meanpk
end

function J_hashindices(J_values, fhash)
    jindlist = []
    for jval in J_values
        push!(jindlist, findall(fhash["J_L"][:, 1].==jval))
    end
    return vcat(jindlist'...)
end

function J_S1indices(J_values, fhash)
    #Assumes this is applied to an object of length 2+Nf+Nf^2 or 2+Nf
    return J_hashindices(J_values, fhash) .+ 2
end

kbins=convert(Array{Float64, 1}, collect(1:32))
apdsmoothed = imfilter(init, Kernel.gaussian(1.0))
smoothps = calc_1dps_local(apdsmoothed, kbins)
true_ps = calc_1dps_local(true_img, kbins)
initps = calc_1dps_local(init, kbins)
recps = calc_1dps_local(recon_img, kbins)

JS1ind = J_S1indices([0, 1, 2, 3], filter_hash)
clim= (minimum(true_img), maximum(true_img))
p1 = heatmap(recon_img, title="Recon", clim=clim)
p3 = plot(log.(kbins), log.(true_ps), label="True")
plot!(log.(kbins), log.(recps), label="Recon")
plot!(log.(kbins), log.(initps), label="Init")
plot!(log.(kbins), log.(smoothps), label="Smoothed Init")
plot!(title="P(k): Denoising using Shift(init)")
xlabel!("lnk")
ylabel!("lnP(k)")
p2 = heatmap(true_img, title="True Img", clim=clim)
p4= heatmap(init, title="Init Img", clim=clim)
p5= heatmap(apdsmoothed, title="Smoothed init", clim=clim)
residual = recon_img- true_img
rlims = (minimum(residual), maximum(residual))
#symmax = maximum([abs(minimum(residual)), maximum(residual)])
#rg = cgrad(:bwr, [-symmax/2.0, symmax/2.0])
p6 = heatmap(residual, title="Residual: Recon - True", clims=rlims, c=:bwr)
p7 = heatmap(apdsmoothed- true_img, title="Residual: SmoothedInit - True", clims=rlims, c=:bwr)

struesel = Data_Utils.fnlog(DHC_compute_wrapper(true_img, filter_hash, norm=false; dhc_args...))
sinitsel = Data_Utils.fnlog(DHC_compute_wrapper(init, filter_hash, norm=false; dhc_args...))
ssmoothsel = Data_Utils.fnlog(DHC_compute_wrapper(apdsmoothed, filter_hash, norm=false; dhc_args...))
sreconsel = Data_Utils.fnlog(DHC_compute_wrapper(recon_img, filter_hash, norm=false; dhc_args...))
slims = (minimum(struesel[JS1ind]), maximum(struesel[JS1ind]))
cg = cgrad([:blue, :white, :red])
truephi, trueomg = round(struesel[2+filter_hash["phi_index"]], sigdigits=3), round(struesel[2+filter_hash["Omega_index"]], sigdigits=3)
reconphi, reconomg = round(sreconsel[2+filter_hash["phi_index"]], sigdigits=3), round(sreconsel[2+filter_hash["Omega_index"]], sigdigits=3)
smoothphi, smoothomg = round(ssmoothsel[2+filter_hash["phi_index"]], sigdigits=3), round(ssmoothsel[2+filter_hash["Omega_index"]], sigdigits=3)

p8 = heatmap(struesel[JS1ind], title="True Coeffs ϕ=" * string(truephi) * "Ω=" * string(trueomg) , clims=slims, c=cg)
p9 = heatmap(sreconsel[JS1ind], title="Recon Coeffs ϕ=" * string(reconphi) * "Ω=" * string(reconomg), clims=slims, c=cg)
p10 = heatmap(ssmoothsel[JS1ind], title="Smooth Init ϕ=" * string(smoothphi) * "Ω=" * string(smoothomg), clims=slims, c=cg)
p = plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, layout=(5, 2), size=(1800, 2400))
savefig(p, "scratch_NM/NewWrapper/5-30/wphlike_s2r_new.png")
fracres = (init .- true_img)./true_img
fps = (initps .- true_ps)./true_ps
println("Mean Abs Frac, Init = ", round(mean(abs.(fracres)), digits=3), "Smoothed = ", round(mean(abs.((apdsmoothed .- true_img)./true_img)), digits=3), "Recon = ", round(mean(abs.((recon_img .- true_img)./true_img)), digits=3))
println("MSE, Init = ", round(mean((init .- true_img).^2), digits=6), "Smoothed = ", round(mean((apdsmoothed .- true_img).^2), digits=6), "Recon = ", round(mean((recon_img .- true_img).^2), digits=6))
println("Power Spec Frac Res, Init = ", round(mean(abs.(fps)), digits=3), "Smoothed = ", round(mean(abs.((smoothps .- true_ps)./true_ps)), digits=3), "Recon = ", round(mean(abs.((recps .- true_ps)./true_ps)), digits=3))

##5) Iterative: using empirical noisy dbn (for non-Omega) and true: with pixelwise regularization
ARGS_buffer = ["reg", "nonapd", "noiso", "scratch_NM/StandardizedExp/Nx64/", "full_3losstest", "Full+Eps"]
ENV_buffer= "1000"
numfile = Base.parse(Int, ENV_buffer)
println(numfile, ARGS_buffer[1], ARGS_buffer[2])

if ARGS_buffer[1]=="log"
    logbool=true
else
    if ARGS_buffer[1]!="reg" error("Invalid log arg") end
    logbool=false
end

if ARGS_buffer[2]=="apd"
    apdbool = true
else
    if ARGS_buffer[2]!="nonapd" error("Invalid apd arg") end
    apdbool=false
end

if ARGS_buffer[3]=="iso"
    isobool = true
else
    if ARGS_buffer[3]!="noiso" error("Invalid iso arg") end
    isobool=false
end

direc = ARGS_buffer[4] #"../StandardizedExp/Nx64/noisy_stdtrue/" #Change
datfile = direc * "Data_" * string(numfile) * ".jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]

fname_save = direc * "scratch_NM/NewWrapper/5-30/" * string(numfile) * "_try1"  #Change
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, t=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])
dhc_args = Dict(:doS2=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense

optim_settings = Dict([("iterations", 20), ("norm", false), ("minmethod", ConjugateGradient())])
recon_settings = Dict([("log", logbool), ("Invcov_matrix", ARGS_buffer[6]), ("optim_settings", optim_settings), ("eps_value_sfd", 1e-5), ("eps_value_init", 1e-10)]) #Add constraints

recon_settings["datafile"] = datfile

if dhc_args[:iso]
    error("Not constructed for iso")
else #Not iso
    coeffmask = falses(2+Nf+Nf^2)
    coeffmask[Nf+3:end] .= Diagonal(trues(Nf))[:] #Diagonal(triu(trues(Nf, Nf)))[:]
end

#Weight given to terms
lval2 = 0.0
lval3 = 0.0

println("Regularizer Lambda=", round(lval3, sigdigits=3))
#input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov, dhc_args, LossFunc, dLossFunc;
#FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), lambda=0.001, func_specific_params=nothing

#Generate fixed noise maps
Nr=10000
sigma = loaddf["std"]
noisemaps = []
for n=1:Nr
    noisyim = randn((Nx, Nx)).*sigma
    push!(noisemaps, noisyim)
end

s_noisy = DHC_compute_wrapper(loaddf["init"], filter_hash, norm=false; dhc_args...)[coeffmask]
starget = s_noisy #DHC_compute_wrapper(loaddf["true_img"], filter_hash, norm=false; dhc_args...)[coeffmask] #

@time pmap(image -> DHC_compute_wrapper(image, filter_hash; norm=false, dhc_args...)[coeffmask], noisemaps)

function ens_noise()
    empsn = []
    for n=1:Nr
        push!(empsn,  DHC_compute_wrapper(noisemaps[n], filter_hash, norm=false; dhc_args...)[coeffmask])
    end
    return empsn
end

@time ens_noise()

function Loss_noisyfast()
    empsn = []

function Loss_noisy(img_curr, filter_hash, dhc_args; noisemaps=nothing, coeffmask=nothing, starget=nothing)
    #func_specific_params should contain: coeff_mask2, target2, invcov2
    empsn = []
    M = 500
    totdiff = 0
    for n=1:M
        noisy = img_curr + noisemaps[n]
        diffsq = sum((DHC_compute_wrapper(noisy, filter_hash; norm=false, dhc_args...)[coeffmask] .- starget).^2)
        totdiff += diffsq
    end
    totdiff = totdiff/M
    return totdiff
end

function dLoss_noisy!(storage_grad, img_curr, filter_hash, dhc_args; noisemaps=nothing, coeffmask=nothing, starget=nothing, dA=nothing)
    #func_specific_params should contain: coeff_mask2, target2, invcov2
    empsn = []
    M = 500 #length(noisemaps)
    totdiff = zeros(Nx^2, 1)
    for n=1:M
        noisy = img_curr + noisemaps[n]
        diffsq = 2*(DHC_compute_wrapper(noisy, filter_hash; norm=false, dhc_args...)[coeffmask] .- starget)
        wts = zeros(2+Nf+Nf^2)
        wts[coeffmask] .= diffsq
        wts = reshape(wts[Nf+3:end], (Nf, Nf))
        #wts = ReconFuncs.augment_weights_S20(convert(Array{Float64,1}, diffsq), filter_hash, dhc_args, coeffmask)
        derve = Deriv_Utils_New.wst_S20_deriv_sum(noisy, filter_hash, wts)
        totdiff .+= derve
    end
    totdiff ./= M
    storage_grad .= reshape(totdiff, (Nx, Nx))
end

func_specific_params = Dict([(:coeffmask=> coeffmask), (:starget=>starget), (:noisemaps=> noisemaps)])

recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings

if logbool
    error("Not implemented here")
end

res, recon_img = image_recon_derivsum_custom(init, filter_hash, dhc_args, Loss_noisy, dLoss_noisy!; optim_settings=optim_settings, func_specific_params)
heatmap(init)
heatmap(recon_img, title="Using true shift")
heatmap(true_img)

function calc_1dps_local(image, kbins::Array{Float64,1})
    #Assumes uniformly spaced kbins
    Nx = size(image)[1]
    fft_zerocenter = fftshift(fft(image))
    impf = abs2.(fft_zerocenter)
    x = (collect(1:Nx) * ones(Nx)') .- (Nx/2.0)
    y = (ones(Nx) * collect(1:Nx)') .- (Nx/2.0)
    krad  = (x.^2 + y.^2).^0.5
    meanpk = zeros(size(kbins))
    kdel = kbins[2] - kbins[1]
    #println(size(meanpk), " ", kdel)
    for k=1:size(meanpk)[1]
        filt = findall((krad .>= (kbins[k] - kdel./2.0)) .& (krad .<= (kbins[k] + kdel./2.0)))
        meanpk[k] = mean(impf[filt])
    end
    return meanpk
end

function J_hashindices(J_values, fhash)
    jindlist = []
    for jval in J_values
        push!(jindlist, findall(fhash["J_L"][:, 1].==jval))
    end
    return vcat(jindlist'...)
end

function J_S1indices(J_values, fhash)
    #Assumes this is applied to an object of length 2+Nf+Nf^2 or 2+Nf
    return J_hashindices(J_values, fhash) .+ 2
end

kbins=convert(Array{Float64, 1}, collect(1:32))
apdsmoothed = imfilter(init, Kernel.gaussian(1.0))
smoothps = calc_1dps_local(apdsmoothed, kbins)
true_ps = calc_1dps_local(true_img, kbins)
initps = calc_1dps_local(init, kbins)
recps = calc_1dps_local(recon_img, kbins)

JS1ind = J_S1indices([0, 1, 2, 3], filter_hash)
clim= (minimum(true_img), maximum(true_img))
p1 = heatmap(recon_img, title="Recon", clim=clim)
p3 = plot(log.(kbins), log.(true_ps), label="True")
plot!(log.(kbins), log.(recps), label="Recon")
plot!(log.(kbins), log.(initps), label="Init")
plot!(log.(kbins), log.(smoothps), label="Smoothed Init")
plot!(title="P(k): Denoising using Shift(init)")
xlabel!("lnk")
ylabel!("lnP(k)")
p2 = heatmap(true_img, title="True Img", clim=clim)
p4= heatmap(init, title="Init Img", clim=clim)
p5= heatmap(apdsmoothed, title="Smoothed init", clim=clim)
residual = recon_img- true_img
rlims = (minimum(residual), maximum(residual))
#symmax = maximum([abs(minimum(residual)), maximum(residual)])
#rg = cgrad(:bwr, [-symmax/2.0, symmax/2.0])
p6 = heatmap(residual, title="Residual: Recon - True", clims=rlims, c=:bwr)
p7 = heatmap(apdsmoothed- true_img, title="Residual: SmoothedInit - True", clims=rlims, c=:bwr)

struesel = Data_Utils.fnlog(DHC_compute_wrapper(true_img, filter_hash, norm=false; dhc_args...))
sinitsel = Data_Utils.fnlog(DHC_compute_wrapper(init, filter_hash, norm=false; dhc_args...))
ssmoothsel = Data_Utils.fnlog(DHC_compute_wrapper(apdsmoothed, filter_hash, norm=false; dhc_args...))
sreconsel = Data_Utils.fnlog(DHC_compute_wrapper(recon_img, filter_hash, norm=false; dhc_args...))
slims = (minimum(struesel[JS1ind]), maximum(struesel[JS1ind]))
cg = cgrad([:blue, :white, :red])
truephi, trueomg = round(struesel[2+filter_hash["phi_index"]], sigdigits=3), round(struesel[2+filter_hash["Omega_index"]], sigdigits=3)
reconphi, reconomg = round(sreconsel[2+filter_hash["phi_index"]], sigdigits=3), round(sreconsel[2+filter_hash["Omega_index"]], sigdigits=3)
smoothphi, smoothomg = round(ssmoothsel[2+filter_hash["phi_index"]], sigdigits=3), round(ssmoothsel[2+filter_hash["Omega_index"]], sigdigits=3)

p8 = heatmap(struesel[JS1ind], title="True Coeffs ϕ=" * string(truephi) * "Ω=" * string(trueomg) , clims=slims, c=cg)
p9 = heatmap(sreconsel[JS1ind], title="Recon Coeffs ϕ=" * string(reconphi) * "Ω=" * string(reconomg), clims=slims, c=cg)
p10 = heatmap(ssmoothsel[JS1ind], title="Smooth Init ϕ=" * string(smoothphi) * "Ω=" * string(smoothomg), clims=slims, c=cg)
p = plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, layout=(5, 2), size=(1800, 2400))
savefig(p, "scratch_NM/NewWrapper/5-30/wphlike_s1only_new.png")
fracres = (init .- true_img)./true_img
fps = (initps .- true_ps)./true_ps
println("Mean Abs Frac, Init = ", round(mean(abs.(fracres)), digits=3), "Smoothed = ", round(mean(abs.((apdsmoothed .- true_img)./true_img)), digits=3), "Recon = ", round(mean(abs.((recon_img .- true_img)./true_img)), digits=3))
println("MSE, Init = ", round(mean((init .- true_img).^2), digits=6), "Smoothed = ", round(mean((apdsmoothed .- true_img).^2), digits=6), "Recon = ", round(mean((recon_img .- true_img).^2), digits=6))
println("Power Spec Frac Res, Init = ", round(mean(abs.(fps)), digits=3), "Smoothed = ", round(mean(abs.((smoothps .- true_ps)./true_ps)), digits=3), "Recon = ", round(mean(abs.((recps .- true_ps)./true_ps)), digits=3))

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
module STREAMBenchmark

using BenchmarkTools
using Statistics, Downloads
using Base.Threads: nthreads, @threads
using LoopVectorization

include("kernels_allthreads.jl")
include("kernels_nothreads.jl")
include("kernels_nthreads.jl")
include("benchmarks.jl")
include("original.jl")

export memory_bandwidth, benchmark, scaling_benchmark

end


text:
"""

$(TYPEDEF)

Structures to contain the volumes obtained from calculations.

$(TYPEDFIELDS)

"""
@with_kw mutable struct Volume
    total::Float64
    bulk::Float64
    domain::Float64
    shell::Vector{Float64}
end

Volume(nbins::Int) = Volume(0.0, 0.0, 0.0, zeros(Float64, nbins))

function reset!(v::Volume)
    v.total = 0.0
    v.bulk = 0.0
    v.domain = 0.0
    @. v.shell = 0.0
    return nothing
end

#function Base.show(io::IO, v::Volume)
#  n = length(v.shell)
#  println(" Mean total box volume: $(v.total) ")
#  println(" Mean bulk volume: $(v.bulk) ")
#  println(" Mean solute domain volume: $(v.domain) ")
#  println(" Volumes of first, medium, and last solvation shells: $(v.shell[1]), $(v.shell[round(Int,n/2)]), $(v.shell[n])")
#end


text:
@everywhere using Evolution
using Base.Test
using ValueHistories
using UnicodePlots

@everywhere function rosenbrock(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

history = VectorUnivalueHistory(Float64)

population = [MinimizerPhenotype(rand(2) * 30, rosenbrock) for i = 1:1000]
for i = 1:30
    map!(fitness!, population)
    mutate!(population, Elitism(RandomNormal(), 0.1))
    crossover!(population, Elitism(Random(), 0.2))

    push!(history, i, fitness(best(population)))
end

solution = best(population)
println(solution)
print(lineplot(get(history)..., height = 10, title = "Best Fitness over Generations"))


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(Blosc_jll.__init__)})
end


text:
@testset "weeks test"  begin
    @test isa(Weeks(s -> 1 / s), Weeks)
    @test isa(WeeksErr(s -> 1 / s), WeeksErr)

    fl = Weeks(s -> 1 / s^2)

    @test isapprox(fl(1), 1)
    @test isapprox(fl([1 2; 3 4]), [1 2; 3 4])

    fl = Weeks(s -> s / (1 + s^2), 64)
    @test isapprox(fl([1 2; 3 4]), broadcast(cos,[1.0 2.0; 3.0 4.0]))

    fle = WeeksErr(s -> s/(1+s^2),  64)
    @test isapprox(fle([1 2; 3 4])[1], cos.([1 2; 3 4]))

    ## Optim is broken.
    # e1 = abs(fl(10.0) - cos(10.0))
    # InverseLaplace.optimize(fl, 10.0)
    # e2 = abs(fl(10.0) - cos(10.0))
    # @test e2 < e1

    # e1 = abs(fle(10.0)[1] - cos(10.0))
    # optimize(fle,10.0)
    # e2 = abs(fl(10.0)[1] - cos(10.0))
    # @test e2 < e1

    fle = WeeksErr(s -> s/(1+s^2),  64)
    (t1,e1) = fle(10.0)
    @test isapprox([t1,e1], [-0.8390715290763104,1.6069233509475447e-11]; atol = 1e-9)

    setNterms(fle, 8)
    e1 = abs(fle(10.0)[1] - cos(10.0))
    setNterms(fle, 80)
    e2 = abs(fle(10.0)[1] - cos(10.0))
    @test e2 < e1

    fle = WeeksErr(s -> s/(1+s^2),  64)
    c1 = copy(fle.coefficients)
    @test string(fle) == "WeeksErr{Float64}(Nterms=64,sigma=1.0,b=1.0)"
    setparameters(fle,2.0,2.0,80)
    @test string(fle) == "WeeksErr{Float64}(Nterms=80,sigma=2.0,b=2.0)"

    # Check that Laguerre coefficients have been recomputed
    @test c1 != fle.coefficients

    ### Complex

    function Fcomplex(s)
        # Laplace domain
        α = complex(-0.3, 6.0)
        return 1 / (s - α)
    end

    function fcomplex(t)
        # Time domain
        α = complex(-0.3, 6.0)
        return exp(α * t)
    end

    let Fc = Weeks(Fcomplex, 1024, datatype=Complex),  trange = range(0.0, stop=30.0, length=1000)
        @test isapprox(Fc.(trange), fcomplex.(trange), atol=0.001)
    end
end # @testset


text:
using ApproxFun, MultivariateOrthogonalPolynomials, StaticArrays, FillArrays, BlockArrays, Test
import ApproxFunBase: checkpoints, plan_transform!, fromtensor
import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendreconic!, legendre2duffyconic!, c_plan_rottriangle, plan_transform,
                        c_execute_tri_hi2lo, c_execute_tri_lo2hi, duffy2legendrecone_column_J!, duffy2legendrecone!, legendre2duffycone!,
                        duffy2legendreconic!, ConicTensorizer, pointsize, tensortransform, _coefficients!

@testset "Conic" begin
    @testset "tensorizer derivation" begin
        A = [1 2 3 5 6;
             4 7 8 0 0;
             9 0 0 0 0]
        B = PseudoBlockArray(A, Ones{Int}(3), [1; Fill(2,2)])

        a = Vector{eltype(A)}()
        for N = 1:blocksize(B,2), K=1:N
            append!(a, vec(B[Block(K,N-K+1)]))
        end

        @test a == fromtensor(ConicTensorizer(),A) == fromtensor(LegendreConic(),A) == 1:9

        N = isqrt(length(a))
        M = 2N-1
        Ã = zeros(eltype(a), N, M)
        B = PseudoBlockArray(Ã, Ones{Int}(3), [1; Fill(2,2)])
        k = 1
        for N = 1:blocksize(B,2), K=1:N
            V = view(B, Block(K,N-K+1))
            for j = 1:length(V)
                V[j] = a[k]
                k += 1
            end
        end
        @test Ã == totensor(ConicTensorizer(),a) == totensor(LegendreConic(),a) == A
    end

    @testset "DuffyConic" begin
        f = Fun((t,x,y) -> 1, DuffyConic(), 10)
        @test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1

        f = Fun((t,x,y) -> t, DuffyConic(), 10)
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)

        f = Fun((t,x,y) -> x, DuffyConic(), 10)

        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 0.1

        f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyConic(), 1000)
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)

        f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyConic())
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ exp(cos(sqrt(0.1^2+0.2^2)*0.1)*0.2)

        m,ℓ = (1,1)
        f = (txy) -> ((t,x,y) = txy;  θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
        g = Fun(f, DuffyConic())
        t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
        @test g(t,x,y) ≈ f((t,x,y))
    end

    @testset "Legendre<>DuffyConic" begin
        for k = 0:10
            a = [zeros(k); 1.0; zeros(5)]
            F = totensor(rectspace(DuffyConic()), a)
            F = pad(F, :, 2size(F,1)-1)
            T = eltype(a)
            P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
            @test legendre2duffyconic!(P, duffy2legendreconic!(P, copy(F))) ≈ F

            b = coefficients(a, LegendreConic(), DuffyConic())
            @test a ≈ coefficients(b, DuffyConic(), LegendreConic())[1:length(a)]
        end
    end

    @testset "LegendreConicPlan" begin
        p = points(LegendreConic(),10)
        @test length(p) == 6
        v = fill(1.0,length(p))
        n = length(v)
        N = (1 + isqrt(1+8n)) ÷ 4
        M = 2N-1
        D = plan_transform!(rectspace(DuffyConic()), reshape(v,N,M))
        N,M = D.plan[1][2],D.plan[2][2]
        V=reshape(v,N,M)
        @test D*V ≈ [[1 zeros(1,2)]; zeros(1,3)]
        T = Float64
        P = c_plan_rottriangle(N, zero(T), zero(T), zero(T))
        duffy2legendreconic!(P,V)
        @test V ≈ [[1 zeros(1,2)]; zeros(1,3)]
        @test_broken fromtensor(DuffyConic(), V) ≈ [1; Zeros(5)]

        v = fill(1.0,length(p))
        @test plan_transform(LegendreConic(),v)*v ≈ [1; Zeros(3)]
        @test v == fill(1.0,length(p))
    end

    @testset "LegendreConic" begin
        f = Fun((t,x,y) -> 1, LegendreConic(), 10)
        @test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1

        f = Fun((t,x,y) -> t, LegendreConic(), 10)
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ sqrt(0.1^2+0.2^2)
        f = Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 10)
        @test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1+sqrt(0.1^2+0.2^2)+0.1+0.2

        @time Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 1000)

        f = Fun((t,x,y) -> exp(x*cos(t+y)), LegendreConic())

        for (m,ℓ) in ((0,0), (0,1), (0,2), (1,1), (1,2), (2,2))
            f = (txy) -> ((t,x,y) = txy;  θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
            g = Fun(f, LegendreConic())
            t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
            @test g(t,x,y) ≈ f((t,x,y))
        end
    end
end

@testset "Cone" begin
    @testset "ConeTensorizer" begin
        ts = MultivariateOrthogonalPolynomials.ConeTensorizer()
        @test totensor(ts,1:10) == [1 2 3 5 6 7; 4 8 9 0 0 0; 10 0 0 0 0 0]
    end

    @testset "DuffyCone" begin
        p = points(DuffyCone(), 10)
        @test p isa Array{SVector{3,Float64},3}
        V = fill(1.0, size(p))
        M,_,N = size(V)
        R = Array{Float64}(undef,M,sum(1:M))
        P = plan_transform(DuffyCone(), Array{Float64}(undef, size(p)))
        for k = 1:M
            R[k,:] = P.plan[2]*V[k,:,:]
        end
        for j = 1:size(M,2)
            R[:,j] = P.plan[1]*R[:,j]
        end
        @test R ≈ [1/sqrt(2) 0 0; 0 0 0]

        @test P * V ≈ [1/sqrt(2); Zeros(9)] ≈ [Fun((x,y) -> 1, ZernikeDisk()).coefficients; Zeros(9)]

        f = Fun((t,x,y) -> 1, DuffyCone(), 10)
        @test f.coefficients ≈ [1/sqrt(2); Zeros(9)]
        @test f(0.3,0.1,0.2) ≈ 1

        f = Fun((t,x,y) -> t, DuffyCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.3

        g = Fun(ZernikeDisk(),[0,1.0])
        p = points(DuffyCone(), 10)

        gg = txy -> ((t,x,y) = txy; g(x/t,y/t))
        V = gg.(p)
        P = plan_transform(DuffyCone(), V)
        @test P*V ≈ [0;1;Zeros(8)] ≈ Fun((t,x,y) -> g(x/t,y/t), DuffyCone(), 10).coefficients

        f = Fun((t,x,y) -> x, DuffyCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.1

        f = Fun((t,x,y) -> y, DuffyCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.2

        f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyCone(), 2000)
        @test f(0.3,0.1,0.2) ≈ exp(cos(0.3*0.1)*0.2)

        f = Fun((t,x,y) -> exp(cos(t*x)*y), DuffyCone())
        @test f(0.3,0.1,0.2) ≈ exp(cos(0.3*0.1)*0.2)
    end

    @testset "TriTransform" begin
        m = 1
        p = Fun(NormalizedJacobi(0,2m+2,0..1), [1])
        q = Fun(x -> p(x)*(1-x)^m * 2^m, NormalizedJacobi(0,2,0..1))
        F = [zeros(2)  q.coefficients]
        F = pad(F, :, 2size(F,1)-1)
        T = Float64
        P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
        c_execute_tri_lo2hi(P, F)
        @test F ≈ (B = zeros(size(F)); B[1,m+1] = 1; B)
    end

    @testset "Legendre<>DuffyCone" begin
        for (m,k,ℓ) in ((0,0,0), (0,0,1), (0,0,2), (1,0,0), (1,0,1), (1,1,0), (1,1,1), (1,1,2))
            Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
            f = (txy) -> ((t,x,y) = txy;  θ = atan(y,x); Fun(NormalizedJacobi(0,2m+2,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * Y(x/t,y/t))
            g = Fun(f, DuffyCone())
            a = g.coefficients
            F = totensor(rectspace(DuffyCone()), a)
            for j = 1:size(F,2)
                F[:,j] = coefficients(F[:,j], NormalizedJacobi(0,1,Segment(1,0)), NormalizedJacobi(0,2,Segment(1,0)))
            end

            T = eltype(a)
            P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))

            Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
            for J = 1:size(F,2)
                duffy2legendrecone_column_J!(P, Fc, F, J)
            end

            @test F ≈ (B = zeros(size(F)); B[ℓ+1,sum(1:m)+k+1] = 1; B)
        end
        for k = 0:10
            a = [zeros(k); 1.0; zeros(5)]
            F = totensor(rectspace(DuffyCone()), a)
            F = pad(F, :, 2size(F,1)-1)

            T = eltype(a)
            P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))

            @test legendre2duffycone!(P, duffy2legendrecone!(P, copy(F))) ≈ F

            b = coefficients(a, LegendreCone(), DuffyCone())
            @test a ≈ coefficients(b, DuffyCone(), LegendreCone())[1:length(a)]
        end
    end

    @testset "LegendreConePlan" begin
        p = points(LegendreCone(),10)
        @test length(p) == 12 == 2*2*3
        @test size(p) == (2,2,3)
        V = fill(1.0,size(p))
        M,_,N = size(V)
        D = plan_transform(DuffyCone(), V)
        C = tensortransform(D,V)
        @test C ≈ [[1/sqrt(2) zeros(1,2)]; zeros(1,3)]
        T = Float64
        P = c_plan_rottriangle(M, zero(T), zero(T), one(T))
        _coefficients!(P,C,DuffyCone(),LegendreCone())
        @test C ≈ [[0.8164965809277259 zeros(1,2)]; zeros(1,3)]
        @test fromtensor(LegendreCone(), C) ≈ [0.8164965809277259; Zeros(9)]

        P = plan_transform(LegendreCone(),V)
        @test P*V ≈ [0.8164965809277259; Zeros(9)]
        @test V == fill(1.0,size(p))

        p = points(LegendreCone(),20)
        V = fill(1.0,size(p))
        P = plan_transform(LegendreCone(),V)
        @test P*V ≈ [0.8164965809277259; Zeros(55)]
    end
    @testset "iduffy" begin
        p = points(LegendreCone(),20)
        V = fill(1.0,size(p))
        D = plan_transform(DuffyCone(),V)
        P = plan_transform(LegendreCone(),V)

    end

    @testset "LegendreCone" begin
        f = Fun((t,x,y) -> 1, LegendreCone(), 10)
        @test f.coefficients ≈ [0.8164965809277259; zeros(ncoefficients(f)-1)]
        @test f(0.3,0.1,0.2) ≈ 1

        f = Fun((t,x,y) -> t, LegendreCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.3
        f = Fun((t,x,y) -> x, LegendreCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.1
        f = Fun((t,x,y) -> y, LegendreCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 0.2
        f = Fun((t,x,y) -> 1+t+x+y, LegendreCone(), 10)
        @test f(0.3,0.1,0.2) ≈ 1+0.3+0.1+0.2
        f = Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 1000)
        @test f(0.3,0.1,0.2) ≈ 1+0.3+0.1+0.2
        @test ncoefficients(f) == 1771

        for (m,k,ℓ) in ((0,0,0), (1,0,0), (0,0,1), (1,1,0), (2,1,1))
            Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
            f = (txy) -> ((t,x,y) = txy; Fun(NormalizedJacobi(0,2m+2,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * Y(x/t,y/t))
            g = Fun(f, LegendreCone(),100)
            @test sum(g.coefficients) ≈ 1
            t,x,y = 0.3,0.1,0.2
            @test g(t,x,y) ≈ f((t,x,y))
        end
    end
end


text:
function zoo_min(obj::Objective, par::Parameter)
    obj_clean_history(obj)
    if isnull(par.constraint)
        if par.algorithm == "pracos"
            algorithm = pracos_opt!
            zoolog("using algorithm: pracos")
        elseif par.algorithm == "psracos"
            algorithm = psracos_opt!
            zoolog("using algorithm: psracos")
        else
            algorithm = asracos_opt!
            zoolog("using algorithm: asracos")
        end
    else
        algorithm = pposs_opt!
        zoolog("using algorith: pposs")
    end
    solution = algorithm(obj, par)
    return solution
end


text:
function generate_synthetic_data(n, param_frac, sigma; signal_to_noise=10,
    fixed_rank_sparse=false, fixed_rank=1, fixed_sparse=0)
    """
    This function samples synthetic data by adding noise to the sum of a
    random low rank matrix and a random sparse matrix.

    :param n: The row and column dimension of the data to be sampled (Int64).
    :param param_frac: Parameter that controls the rank and sparsity of the
                       sampled low rank and sparse matrix when
                       fixed_rank_sparse=false (Float64).
    :param sigma: Parameter that controls the absolute magnitude of the noise
                  (Float64).
    :param signal_to_noise: Parameter that controls the signal to noise ratio of
                            the sampled data (Float64).
    :param fixed_rank_sparse: If true, use fixed_rank and fixed_sparse as the
                              rank and sparsity respectively of the low rank
                              matrix and the sparse matrix to be sampled (Bool).
    :param fixed_rank: Rank of the low rank matrix to be sampled when
                       fixed_rank_sparse=true.
    :param fixed_sparse: Sparsity of the sparse matrix to be sampled when
                         fixed_rank_sparse=true

    :return: This function returns 5 values:
             1) The sampled n x n data matrix.
             2) The sampled n x n low rank matrix.
             3) The sampled n x n sparse matrix.
             4) The sparsity level of the sparse matrix.
             5) The rank of the low rank matrix.
    """

    if fixed_rank_sparse
        k_sparse = fixed_sparse
        k_rank = fixed_rank
    else
        k_sparse = Int64(floor(param_frac * n^2))
        k_rank = Int64(floor(param_frac * n))
    end

    # Sample a rank k_rank matrix
    d = Normal(0, signal_to_noise * sigma / n^0.5)
    X = rand(d, (n, k_rank))
    Y = rand(d, (n, k_rank))
    L_0 = X * Y'

    # Sample a sparsity k_sparse matrix
    d = Uniform(-5, 5)
    E = rand(d, (n, n))
    S_0 = zeros(n, n)

    indexes = []
    for i = 1:n
        for j = 1:n
            append!(indexes, [(i, j)])
        end
    end

    select_indexes = sample(indexes, k_sparse, replace=false)
    for index in select_indexes
        S_0[index[1], index[2]] = E[index[1], index[2]]
    end

    data = L_0 + S_0

    # Add noise to the sum of the sampled low rank and sparse matrices
    d = Normal(0, sigma^2)
    noise = rand(d, (n, n))
    data = data + noise

    return data, L_0, S_0, k_sparse, k_rank

end;

function generate_synthetic_data_symmetric(n, param_frac, sigma;
    signal_to_noise=10, fixed_rank_sparse=false, fixed_rank=1, fixed_sparse=0)
    """
    This function samples symmetric synthetic data by adding symmetric noise to
    the sum of a random symmetric low rank matrix and a random symmetric sparse
    matrix.

    :param n: The row and column dimension of the data to be sampled (Int64).
    :param param_frac: Parameter that controls the rank and sparsity of the
                       sampled low rank and sparse matrix when
                       fixed_rank_sparse=false (Float64).
    :param sigma: Parameter that controls the absolute magnitude of the noise
                  (Float64).
    :param signal_to_noise: Parameter that controls the signal to noise ratio of
                            the sampled data (Float64).
    :param fixed_rank_sparse: If true, use fixed_rank and fixed_sparse as the
                              rank and sparsity respectively of the low rank
                              matrix and the sparse matrix to be sampled (Bool).
    :param fixed_rank: Rank of the low rank matrix to be sampled when
                       fixed_rank_sparse=true.
    :param fixed_sparse: Sparsity of the sparse matrix to be sampled when
                         fixed_rank_sparse=true

    :return: This function returns 5 values:
             1) The sampled n-by-n data matrix.
             2) The sampled n-by-n low rank matrix.
             3) The sampled n-by-n sparse matrix.
             4) The sparsity level of the sparse matrix.
             5) The rank of the low rank matrix.
    """

    if fixed_rank_sparse
        k_sparse = fixed_sparse
        k_rank = fixed_rank
    else
        k_sparse = Int64(floor(param_frac * n^2))
        k_rank = Int64(floor(param_frac * n))
    end

    # Sample a rank k_rank symmetric matrix
    d = Normal(0, signal_to_noise * sigma / n^0.5)
    X = rand(d, (n, k_rank))
    L_0 = X * X'

    # Sample a sparsity k_sparse symmetric matrix
    d = Uniform(-5, 5)
    E = rand(d, (n, n))
    S_0 = zeros(n, n)

    indexes = []
    for i = 1:(n-1)
        for j = (i+1):n
            append!(indexes, [(i, j)])
        end
    end

    select_indexes = sample(indexes, Int64(floor(k_sparse/2)), replace=false)
    for index in select_indexes
        S_0[index[1], index[2]] = E[index[1], index[2]]
        S_0[index[2], index[1]] = S_0[index[1], index[2]]
    end

    if k_sparse % 2 == 1
        index = sample(1:n, 1)
        S_0[index, index] = E[index, index]
    end

    data = L_0 + S_0

    # Add symmetric noise to the sum of the sampled low rank and sparse matrices
    d = Normal(0, sigma^2)
    noise = rand(d, (n, n))
    for i = 1:(n-1)
        for j = (i+1):n
            noise[j, i] = noise[i, j]
        end
    end
    data = data + noise

    return data, L_0, S_0, k_sparse, k_rank

end;