text:
const generic_docs = "
- The first form accepts inputs as matrices, or file names of the data bands.
- The second form is more versatile but also more complex to describe.
  - `cube`: Is the file name of a 'cube', a multi-layered file normally created with the [`cutcube`](@ref) function.
     If this file was created with band descriptions one can use the `bands` or the `bandnames` options.
  - `bands`: _cubes_ created with [`cutcube`](@ref) assign descriptions starting with \"Band 1 ...\" an so on
    the other bands. So when `bands` is used we search for bands named \"Band band[k]\", where band[k] loops
    over all elements of the `bands` vector. WARNING: the elements order in the vector must be sorted in increasing
    wavelength numbers, _i.e._ like the example for the first form.
  - `layers`: Use this option when you are certain of the bands order in the cube or the it doesn't have a bands
    description. The selection will be made with cube[:,:,layer[1]], etc... WARNING: same warn as above.
  - `bandnames`: When we know the common designation of a band, for example \"Green\", or any part of a band
    description, for example \"NIR\", we can use that info to create a `bandnames` string vector that will be
    matched against the cube's bands descriptions.
- `kwargs`:
  - `threshold`: When a threshold is provided we return a GMTgrid where `vals[ij] < threshold = NaN`
  - `classes`: is a vector with up to 3 elements (class separators) and we return a  UInt8 GMTimage with the
    indices categorized into vals[ij] > classes[1] = 1; vals[ij] > classes[2] = 2; vals[ij] > classes[3] = 3 and 0 otherwise.
  - `mask`: Used together with `threshold` outputs a UInt8 GMTimage mask with `vals[ij] >= threshold = 255` and 0 otherwise
     If `mask=-1` (or any other negative number) we compute instead a mask where `vals[ij] < threshold = 255` and 0 otherwise
  - `save`: Use `save=\"file_name.ext\"` to save the result in a disk file. File format is picked from file extension.

If none of `bands`, `layers` or `bandnames` is provided, we use the default band names shown in the first form.

See also https://www.indexdatabase.de/ for a list of indices and the appropriate band names per sensor.

Returns either a Float32 GMTgrid or a UInt8 GMTimage if the `mask` or `classes` options are used.
"

# ----------------------------------------------------------------------------------------------------------
function helper_si_method(cube::String, index::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[],
	                      bandnames::Vector{String}=String[], defbandnames::Vector{String}=String[], kw...)
	# Helper function to compute Spectral Indices from a 'cube' file and somehow band selection.
	(isempty(bands) && isempty(bandnames) && isempty(layers)) && (bandnames = defbandnames)
	sc = subcube(cube, bands=bands, layers=layers, bandnames=bandnames)
	sp_indices(sc, collect(1:size(sc,3)); index=index, kw...)	# Here we know that the layers are all of those in scube
end

# ----------------------------------------------------------------------------------------------------------
"""
    CLG = clg(green, redEdge3; kw...)

Green cholorphyl index. Wu et al 2012.

CLG = (redEdge3)/(green)-1
"""
clg(green, redEdge3; kw...) = sp_indices(green, redEdge3; index="CLG", kw...)
clg(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="CLG", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    CLRE = clre(redEdge1, redEdge3; kw...)

RedEdge cholorphyl index. Clevers and Gitelson 2013.

CLRE = (redEdge3)/(redEdge1)-1
"""
clre(redEdge1, redEdge3; kw...) = sp_indices(redEdge1, redEdge3; index="CLRE", kw...)
clre(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="CLRE", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    EVI = evi(blue, red, nir; kw...)
or

    EVI = evi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Enhanced vegetation index. Huete et al 1990

EVI = G * ((nir - red) / (nir + C1 * red - C2 * blue + Levi));
C1, C2, G, Levi = 6.0, 7.5, 2.5, 1.

$(generic_docs)

"""
evi(blue, red, nir; kw...) = sp_indices(blue, red, nir; index="EVI", kw...)
evi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 3); index="EVI", kw...)
evi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "EVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["blue", "red", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    EVI2 = evi2(red, nir; kw...)
or

    EVI2 = evi2(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Two-band Enhanced vegetation index. Jiang et al 2008

EVI2 = G * ((nir - red) / (nir + 2.4 * red ))

$(generic_docs)
"""
evi2(red, nir; kw...) = sp_indices(red, nir; index="EVI2", kw...)
evi2(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="EVI2", kw...)
evi2(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "EVI2"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    GNDVI = gndvi(green, nir; kw...)
or

    GNDVI = gndvi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

green Normalized diff vegetation index: more sensitive to cholorphyll than ndvi. Gitelson, A., and M. Merzlyak

GNDVI = (nir - green) / (nir + green)

$(generic_docs)
"""
gndvi(green, nir; kw...) = sp_indices(green, nir; index="GNDVI", kw...)
gndvi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="GNDVI", kw...)
gndvi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "GNDVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["green", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    MNDWI = mndwi(green, swir2; kw...)
or

    MNDWI = mndwi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Modified Normalised Difference Water Index. Xu2006

MNDWI = (green-swir2) / (green+swir2)

$(generic_docs)
"""
mndwi(green, swir2; kw...) = sp_indices(swir2, green; index="MNDWI", kw...)
mndwi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="MNDWI", kw...)
mndwi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "MNDWI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["green", "swir 2"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    MTCI = mtci(red, redEdge1, redEdge2; kw...)

Meris Terrestrial Chlorophyll Index. Clevers and Gitelson 2013, Dash and Curran 2004

MTCI = (redEdge2-redEdge1) / (redEdge1-red)
"""
mtci(red, redEdge1, redEdge2; kw...) = sp_indices(red, redEdge1, redEdge2; index="MTCI", kw...)
mtci(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 3); index="MTCI", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    MCARI = mcari(green, red, redEdge1; kw...)

Modified Chlorophyll Absorption ratio index. Daughtery et al. 2000

MCARI = (redEdge1 - red - 0.2 * (redEdge1 - green)) * (redEdge1 / red)
"""
# Sentinel-2 Band 5 (VNIR), Band 4 (Red) and Band 3 (Green).
mcari(green, red, redEdge1; kw...) = sp_indices(green, red, redEdge1; index="MCARI", kw...)
mcari(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 3); index="MCARI", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    MSAVI = msavi(red, nir; kw...)
or

    MSAVI = msavi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Modified soil adjusted vegetation index. Qi 1994

MSAVI = nir + 0.5 - (0.5 * sqrt(pow(2.0 * nir + 1.0, 2) - 8.0 * (nir - (2.0 * red))))

$(generic_docs)
"""
msavi(red, nir; kw...) = sp_indices(red, nir; index="MSAVI", kw...)
msavi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="MSAVI", kw...)
msavi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "MSAVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NBRI = nbri(nir, swir3; kw...)

Normalised Burn Ratio Index. Garcia 1991

NBRI = (nir - swir3) / (nir + swir3)
"""
nbri(nir, swir3; kw...) = sp_indices(swir3, nir; index="NBRI", kw...)
nbri(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NBRI", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NDVI = ndvi(red, nir; kw...)
or

    NDVI = ndvi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Compute the NDVI vegetation index. Input can be either the bands file names, or GMTimage objects
with the band's data.

NDVI = (nir - red) / (nir + red)

$(generic_docs)
"""
ndvi(red, nir; kw...) = sp_indices(red, nir; index="NDVI", kw...)
ndvi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NDVI", kw...)
ndvi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "NDVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NDWI = ndwi(green, nir; kw...)
or

    NDWI = ndwi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Normalized difference water index. McFeeters 1996. NDWI => (green - nir)/(green + nir)

NDWI = (green - nir)/(green + nir)

$(generic_docs)
"""
ndwi(green, nir; kw...) = sp_indices(nir, green; index="NDWI", kw...)
ndwi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NDWI", kw...)
ndwi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "NDWI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["green", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NDWI2 = ndwi2(nir, swir2; kw...)
or

    NDWI2 = ndwi2(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Normalized difference water index. Gao 1996, Chen 2005 (also known as Normalized Difference Moisture Index
NDBI and LSWI)

NDWI2 = (nir - swir2)/(nir + swir2)

$(generic_docs)
"""
ndwi2(nir, swir2; kw...) = sp_indices(swir2, nir; index="NDWI2", kw...)
ndwi2(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NDWI2", kw...)
ndwi2(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "NDWI2"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["nir", "swir 2"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NDREI1 = ndrei1(redEdge1, redEdge2; kw...)

Normalized difference red edge index. Gitelson and Merzlyak 1994

NDREI1 = (redEdge2 - redEdge1) / (redEdge2 + redEdge1)
"""
ndrei1(redEdge1, redEdge2; kw...) = sp_indices(redEdge1, redEdge2; index="NDREI1", kw...)
ndrei1(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NDREI1", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    NDREI2 = ndrei2(redEdge1, redEdge3; kw...)

Normalized difference red edge index 2. Barnes et al 2000

NDREI2 = (redEdge3 - redEdge1) / (redEdge3 + redEdge1)
"""
ndrei2(redEdge1, redEdge3; kw...) = sp_indices(redEdge1, redEdge3; index="NDREI2", kw...)
ndrei2(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="NDREI2", kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    SATVI = satvi(red, swir2, swir3; kw...)

Soil adjusted total vegetation index. Marsett 2006

SATVI = ((swir1 - red) / (swir1 + red + L)) * (1.0 + L) - (swir2 / 2.0)
"""
satvi(red, swir2, swir3; kw...) = sp_indices(red, swir2, swir3; index="SATVI", kw...)
satvi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 3); index="SATVI", kw...)
satvi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "SATVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "swir 1", "swir 2"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    SAVI = savi(red, nir; kw...)
or

    SAVI = savi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Soil adjusted vegetation index. Huete 1988

SAVI = (nir - red) * (1.0 + L) / (nir + red + L)

$(generic_docs)
"""
savi(red, nir; kw...) = sp_indices(red, nir; index="SAVI", kw...)
savi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 2); index="SAVI", kw...)
savi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "SAVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "nir"], kw...)

# ----------------------------------------------------------------------------------------------------------
"""
    SLAVI = slavi(red, nir, swir2; kw...)
or

    SLAVI = slavi(cube::String; [bands=Int[], bandnames=String[], layers=Int[]], kwargs...)

Specific Leaf Area Vegetation Index. Lymburger 2000

SLAVI = nir / (red + swir2)

$(generic_docs)
"""
slavi(red, nir, swir2; kw...) = sp_indices(red, nir, swir2; index="SLAVI", kw...)
slavi(cube::GMT.GMTimage{UInt16, 3}, bnds; kw...) = sp_indices(cube, find_layers(cube, bnds, 3); index="SLAVI", kw...)
slavi(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], kw...) =
	helper_si_method(cube, "SLAVI"; bands=bands, layers=layers, bandnames=bandnames, defbandnames=["red", "nir", "swir 2"], kw...)

# ----------------------------------------------------------------------------------------------------------
function helper_sp_indices(kwargs...)
	# Helper function that is called by two different sp_indices methods
	d = KW(kwargs)
	mask = ((val = find_in_dict(d, [:mask], false)[1]) !== nothing)
	rev_mask = (mask && val < 0) ? true : false		# See if we want to reverse the mask
	threshold = find_in_dict(d, [:threshold], false)[1]
	classes = (threshold === nothing) ? find_in_dict(d, [:classes])[1] : nothing
	(mask && threshold === nothing) && error("The `mask` option requires the `threshold=x` option")
	(classes !== nothing && length(classes) > 3) && (classes = classes[1:3]; @warn("`classes` maximum elements is 3. Clipping the others"))
	save_name::String = ((val = find_in_dict(d, [:save], false)[1]) !== nothing) ? val : ""
	dbg = (find_in_dict(d, [:Vd :dbg])[1] !== nothing)
	dd = d		# Copy to report unused keys (errors) but can't do it in 'd' because this function may be called twice
	(haskey(dd, :mask))      && delete!(dd, :mask)
	(haskey(dd, :classes))   && delete!(dd, :classes)
	(haskey(dd, :threshold)) && delete!(dd, :threshold)
	(length(d) > 0) && println("Warning: the following options were not consumed in sp_indices => ", keys(dd))
	return mask, rev_mask, classes, threshold, save_name, dbg
end

# ----------------------------------------------------------------------------------------------------------
function sp_indices(bnd1::String, bnd2::String, bnd3::String=""; index::String="", kwargs...)
	# Compute spectral indices
	do_radTOA = any(keys(kwargs) .== :dn2radiance)
	do_refTOA = any(keys(kwargs) .== :dn2reflectance)
	do_refSrf = any(keys(kwargs) .== :reflectance_surf)
	Bnd1 = (do_radTOA) ? dn2radiance(bnd1) : (do_refTOA) ? dn2reflectance(bnd1) : (do_refSrf) ? reflectance_surf(bnd1) : gmtread(bnd1)
	Bnd2 = (do_radTOA) ? dn2radiance(bnd2) : (do_refTOA) ? dn2reflectance(bnd2) : (do_refSrf) ? reflectance_surf(bnd2) : gmtread(bnd2)
	if (bnd3 != "")
		Bnd3 = (do_radTOA) ? dn2radiance(bnd3) : (do_refTOA) ? dn2reflectance(bnd3) : (do_refSrf) ? reflectance_surf(bnd3) : gmtread(bnd3)
	else
		Bnd3 = nothing
	end
	sp_indices(Bnd1, Bnd2, Bnd3; index=index, kwargs...)
end

function sp_indices(cube::GMT.GMTimage{UInt16, 3}, bands::Vector{Int}; index::String="", kw...)
	# This method recieves the cube and a vector with the bands list and calls the worker with @view
	mask, _, classes, _, save_name, dbg = helper_sp_indices(kw...)	# Do this first because if it errors no point in continuing
	(dbg) && println(cube.names)
	if (length(bands) == 2)
		o = sp_indices(@view(cube[:,:,bands[1]]), @view(cube[:,:,bands[2]]); index=index, kw...)
	else
		o = sp_indices(@view(cube[:,:,bands[1]]), @view(cube[:,:,bands[2]]), @view(cube[:,:,bands[3]]); index=index, kw...)
	end
	if (mask || classes !== nothing)
		O = mat2img(o, cube)
		O.layout, O.range[5], O.range[6] = "BRPa", 0, (mask) ? 255 : length(classes)
	else
		O = mat2grid(o, cube)
	end
	O.names = [index * " index"]
	(save_name != "") && gmtwrite(save_name, O)
	return (save_name == "") ? O : nothing
end

function sp_indices(bnd1, bnd2, bnd3=nothing; index::String="", kwargs...)
	# This is the method who does the real work.
	(index == "") && error("Must select which index to compute")
	@assert size(bnd1) == size(bnd2)
	mask, rev_mask, classes, threshold, = helper_sp_indices(kwargs...)
	img = (mask || classes !== nothing) ? fill(UInt8(0), size(bnd1)) : fill(NaN32, size(bnd1))

	if (rev_mask) fcomp = <	else  fcomp = >  end 	# Which one to use when computing masks
	mn = size(bnd1,1) * size(bnd1,2)
	i_tmax = (eltype(bnd1) <: Integer) ? 1. / typemax(eltype(bnd1)) : 1.0	# Floats go unchanged
	C1, C2, G, L, Levi = 6.0, 7.5, 2.5, 0.5, 1.0
	if (index == "CLG" || index == "CLRE")
		# Green cholorphyl index. Wu et al 2012		CLG               (redEdge3)/(green)-1
		# RedEdge cholorphyl index. Clevers and Gitelson 2013	CLRE  (redEdge3)/(redEdge1)-1
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t = bnd2[k] / bnd1[k] - 1
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				img[k] = bnd2[k] / bnd1[k] - 1
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "EVI")			# Enhanced vegetation index. Huete et al 1990
		# G * ((nir - red) / (nir + C1 * red - C2 * blue + Levi))
		(bnd3 === nothing) && error("NIR component cannot be empty")
		@assert size(bnd3) == size(bnd2)
		blue = bnd1;	red = bnd2;		nir = bnd3
		_C2 = C2 * i_tmax;
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				t = G * ((t_nir - t_red) / (t_nir + C1 * t_red - _C2 * blue[k] + Levi))
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				img[k] = G * ((t_nir - t_red) / (t_nir + C1 * t_red - _C2 * blue[k] + Levi))
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "EVI2")		# Two-band Enhanced vegetation index. Jiang et al 2008
		# G * ((nir - red) / (nir + 2.4 * red ))
		red = bnd1;		nir = bnd2
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				t = G * (t_nir - t_red) / (t_red + 2.4 * t_nir)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				t = G * (t_nir - t_red) / (t_red + 2.4 * t_nir)
				(t >= -1 && t <= 1) && (img[k] = t)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "MTCI")		# Meris Terrestrial Chlorophyll Index. Clevers and Gitelson 2013, Dash and Curran 2004
		# (redEdge2-redEdge1) / (redEdge1-red)
		red = bnd1;		redEdge1 = bnd2;	redEdge2 = bnd3
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_redEdge1 = redEdge1[k]*i_tmax
				t = (redEdge2[k]*i_tmax - t_redEdge1) / (t_redEdge1 - red[k]*i_tmax)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_redEdge1 = redEdge1[k]*i_tmax
				img[k] = (redEdge2[k]*i_tmax - t_redEdge1) / (t_redEdge1 - red[k]*i_tmax)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "MCARI")		# Modified Chlorophyll Absorption ratio index. Daughtery et al. 2000
		# (redEdge1 - red - 0.2 * (redEdge1 - green)) * (redEdge1 / red)
		# Sentinel-2 Band 5 (VNIR), Band 4 (Red) and Band 3 (Green).
		green = bnd1;	red = bnd2;		redEdge1 = bnd3
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_redEdge1 = redEdge1[k]*i_tmax;	t_red = red[k]*i_tmax
				t = (t_redEdge1 - t_red - 0.2 * (t_redEdge1 - green[k]*i_tmax)) * (t_redEdge1 / t_red)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_redEdge1 = redEdge1[k]*i_tmax;	t_red = red[k]*i_tmax
				img[k] = (t_redEdge1 - t_red - 0.2 * (t_redEdge1 - green[k]*i_tmax)) * (t_redEdge1 / t_red)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "MSAVI")		# Modified soil adjusted vegetation index.
		# nir + 0.5 - (0.5 * sqrt(pow(2.0 * nir + 1.0, 2) - 8.0 * (nir - (2.0 * red))))
		red = bnd1;		nir = bnd2;
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_nir = nir[k]*i_tmax
				t = t_nir + 0.5 - (0.5 * sqrt((2 * t_nir + 1) ^2) - 8 * (t_nir - (2 * red[k]*i_tmax)))
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_nir = nir[k]*i_tmax
				img[k] = t_nir + 0.5 - (0.5 * sqrt((2 * t_nir + 1) ^2) - 8 * (t_nir - (2 * red[k]*i_tmax)))
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "GNDVI" || index == "MNDWI" || index == "NBRI" || index == "NDVI" || index == "NDWI" || index == "NDWI2" || index == "NDREI1" || index == "NDREI2")
		# green Normalized diff vegetation index: more sensitive to cholorphyll than ndvi. Gitelson, A., and M. Merzlyak.
		# GNDVI		=> (nir - green) / (nir + green)
		# Modified Normalised Difference Water Index. Xu2006;  MNDWI	=> (green-swir2) / (green+swir2)
		# Normalised Burn Ratio Index. NBRI => (nir - swir3) / (nir + swir3)
		# Normalized difference water index. McFeeters 1996. NDWI => (green - nir)/(green + nir)
		# NDVI Normalised Difference Vegetation Index. Rouse 1974. (nir - red) / (nir + red)
		# NDBI, LSWI. Normalized difference water index. Gao 1996, Chen 2005; NDWI2 => (nir - swir2)/(nir + swir2)
		# Normalized difference red edge index. Gitelson and Merzlyak 1994; (redEdge2 - redEdge1)/(redEdge2 + redEdge1)
		# Normalized difference red edge index 2. Barnes et al 2000; (redEdge3 - redEdge1)/(redEdge3 + redEdge1)
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t1 = bnd1[k]*i_tmax;	t2 = bnd2[k]*i_tmax
				t = (t2 - t1) / (t1 + t2)
				(fcomp(t, threshold) && t <= 1) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t1 = bnd1[k]*i_tmax;	t2 = bnd2[k]*i_tmax
				t = (t2 - t1) / (t1 + t2)
				(t >= -1 && t <= 1) && (img[k] = t)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "SATVI")		# Soil adjusted total vegetation index.
		# ((swir1 - red) / (swir1 + red + L)) * (1.0 + L) - (swir2 / 2.0)
		red = bnd1;		swir1 = bnd2;	swir2 = bnd3;
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_swir1 = swir1[k]*i_tmax
				t = ((t_swir1 - t_red) / (t_swir1 + t_red + L)) * (1.0 + L) - (swir2[k]*i_tmax / 2.0)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_swir1 = swir1[k]*i_tmax
				img[k] = ((t_swir1 - t_red) / (t_swir1 + t_red + L)) * (1.0 + L) - (swir2[k]*i_tmax / 2.0)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "SAVI")		# Soil adjusted vegetation index. Huete1988
		# (nir - red) * (1.0 + L) / (nir + red + L);
		red = bnd1;		nir = bnd2;
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				t = (t_nir - t_red) * (1.0 + L) / (t_nir - t_red + L)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				t_red = red[k]*i_tmax;	t_nir = nir[k]*i_tmax
				img[k] = (t_nir - t_red) * (1.0 + L) / (t_nir - t_red + L)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	elseif (index == "SLAVI")		# Soil Adjusted Vegetation Index Huete 1988.
		# nir / (red + swir2)
		red = bnd1;		nir = bnd2;		swir2 = bnd3
		if (mask)
			@inbounds Threads.@threads for k = 1:mn
				t = nir[k]*i_tmax / (red[k]*i_tmax + swir2[k]*i_tmax)
				(fcomp(t, threshold)) && (img[k] = 255)
			end
		else
			@inbounds Threads.@threads for k = 1:mn
				img[k] = nir[k]*i_tmax / (red[k]*i_tmax + swir2[k]*i_tmax)
			end
			helper_si!(img, threshold, classes)		# Threshold or Classes if one of them is != nothing
		end
	end
	if (isa(bnd1, GMT.GMTimage) || isa(bnd1, GMT.GMTgrid))
		if (mask || classes !== nothing)
			I = mat2img(img, proj4=bnd1.proj4, wkt=bnd1.wkt, x=bnd1.x, y=bnd1.y)
			I.epsg = bnd1.epsg;		I.layout = "BRPa"
			I.range[5], I.range[6] = 0, (mask) ? 255 : length(classes)
			return I
		else
			return mat2grid(img, bnd1)
		end
	else
		return img
	end
end

function helper_si!(img, threshold, classes)
	# Helper function to apply a threshold to a grid OR create a classes image
	mn = size(img, 1) * size(img, 2)
	if (threshold !== nothing)		# In this case we return one GMTgrid where < thres = NaN
		@inbounds Threads.@threads for k = 1:mn
			(img[k] < threshold) && (img[k] = NaN32)
		end
	elseif (classes !== nothing)	# Here we return a GMTimage with up to 4 classes (maximum allowed)
		@inbounds Threads.@threads for k = 1:mn
			(img[k] >= classes[1]) && (img[k] = 1)
			(img[k] >= classes[2]) && (img[k] = 2)
		end
		if (length(classes) == 3)
			@inbounds Threads.@threads for k = 1:mn
				(img[k] >= classes[3]) && (img[k] = 3)
			end
		end
	end
end

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
#####
##### `dot`
#####

function frule((_, Δx, Δy), ::typeof(dot), x, y)
    return dot(x, y), dot(Δx, y) + dot(x, Δy)
end

function rrule(::typeof(dot), x::AbstractArray, y::AbstractArray)
    project_x = ProjectTo(x)
    project_y = ProjectTo(y)
    function dot_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        x̄ = @thunk(project_x(reshape(y .* ΔΩ', axes(x))))
        ȳ = @thunk(project_y(reshape(x .* ΔΩ, axes(y))))
        return (NoTangent(), x̄, ȳ)
    end
    return dot(x, y), dot_pullback
end

#####
##### 3-arg `dot`
#####

function frule(
    (_, Δx, ΔA, Δy),
    ::typeof(dot),
    x::AbstractVector{<:Number},
    A::AbstractMatrix{<:Number},
    y::AbstractVector{<:Number},
)
    return dot(x, A, y), dot(Δx, A, y) + dot(x, ΔA, y) + dot(x, A, Δy)
end

function rrule(
    ::typeof(dot),
    x::AbstractVector{<:Number},
    A::AbstractMatrix{<:Number},
    y::AbstractVector{<:Number},
)
    project_x = ProjectTo(x)
    project_A = ProjectTo(A)
    project_y = ProjectTo(y)
    Ay = A * y
    z = adjoint(x) * Ay
    function dot_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        dx = @thunk project_x(conj(ΔΩ) .* Ay)
        dA = @thunk project_A(ΔΩ .* x .* adjoint(y))
        dy = @thunk project_y(ΔΩ .* (adjoint(A) * x))
        return (NoTangent(), dx, dA, dy)
    end
    dot_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent(), ZeroTangent())
    return z, dot_pullback
end

function rrule(
    ::typeof(dot),
    x::AbstractVector{<:Number},
    A::Diagonal{<:Number},
    y::AbstractVector{<:Number},
)
    project_x = ProjectTo(x)
    project_A = ProjectTo(A)
    project_y = ProjectTo(y)
    z = dot(x, A, y)
    function dot_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        dx = @thunk project_x(conj(ΔΩ) .* A.diag .* y)  # A*y is this broadcast, can be fused
        dA = @thunk project_A(Diagonal(ΔΩ .* x .* conj(y)))  # calculate N not N^2 elements
        dy = @thunk project_y(ΔΩ .* conj.(A.diag) .* x)
        return (NoTangent(), dx, dA, dy)
    end
    dot_pullback(::ZeroTangent) = (NoTangent(), ZeroTangent(), ZeroTangent(), ZeroTangent())
    return z, dot_pullback
end

#####
##### `cross`
#####

function frule((_, Δa, Δb), ::typeof(cross), a::AbstractVector, b::AbstractVector)
    return cross(a, b), cross(Δa, b) .+ cross(a, Δb)
end

# TODO: support complex vectors
function rrule(::typeof(cross), a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
    project_a = ProjectTo(a)
    project_b = ProjectTo(b)
    Ω = cross(a, b)
    function cross_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        da = @thunk(project_a(cross(b, ΔΩ)))
        db = @thunk(project_b(cross(ΔΩ, a)))
        return (NoTangent(), da, db)
    end
    return Ω, cross_pullback
end

#####
##### `det`
#####

function frule((_, Δx), ::typeof(det), x::AbstractMatrix)
    Ω = det(x)
    # TODO Performance optimization: probably there is an efficent
    # way to compute this trace without during the full compution within
    return Ω, Ω * tr(x \ Δx)
end
frule((_, Δx), ::typeof(det), x::Number) = (det(x), Δx)

function rrule(::typeof(det), x::Union{Number,AbstractMatrix})
    Ω = det(x)
    function det_pullback(ΔΩ)
        ∂x = x isa Number ? ΔΩ : Ω * ΔΩ * inv(x)'
        return (NoTangent(), ∂x)
    end
    return Ω, det_pullback
end

#####
##### `logdet`
#####

function frule((_, Δx), ::typeof(logdet), x::Union{Number,StridedMatrix{<:Number}})
    Ω = logdet(x)
    return Ω, tr(x \ Δx)
end

function rrule(::typeof(logdet), x::Union{Number,StridedMatrix{<:Number}})
    Ω = logdet(x)
    function logdet_pullback(ΔΩ)
        ∂x = x isa Number ? ΔΩ / x' : ΔΩ * inv(x)'
        return (NoTangent(), ∂x)
    end
    return Ω, logdet_pullback
end

#####
##### `logabsdet`
#####

function frule((_, Δx), ::typeof(logabsdet), x::AbstractMatrix)
    Ω = logabsdet(x)
    (y, signy) = Ω
    b = tr(x \ Δx)
    ∂y = real(b)
    ∂signy = eltype(x) <: Real ? ZeroTangent() : im * imag(b) * signy
    ∂Ω = Tangent{typeof(Ω)}(∂y, ∂signy)
    return Ω, ∂Ω
end

function rrule(::typeof(logabsdet), x::AbstractMatrix)
    Ω = logabsdet(x)
    function logabsdet_pullback(ΔΩ)
        (Δy, Δsigny) = ΔΩ
        (_, signy) = Ω
        f = signy' * Δsigny
        imagf = f - real(f)
        g = real(Δy) + imagf
        ∂x = g * inv(x)'
        return (NoTangent(), ∂x)
    end
    return Ω, logabsdet_pullback
end

#####
##### `trace`
#####

function frule((_, Δx), ::typeof(tr), x)
    return tr(x), tr(Δx)
end

function rrule(::typeof(tr), x)
    # This should really be a FillArray
    # see https://github.com/JuliaDiff/ChainRules.jl/issues/46
    d = size(x, 1)
    function tr_pullback(ΔΩ)
        return (NoTangent(), Diagonal(fill(ΔΩ, d)))
    end
    return tr(x), tr_pullback
end

#####
##### `pinv`
#####

@scalar_rule pinv(x) -(Ω^2)

function frule(
    (_, ẋ), ::typeof(pinv), x::AbstractVector{T}, tol::Real=0
) where {T<:Union{Real,Complex}}
    y = pinv(x, tol)
    Δx = unthunk(ẋ)
    ∂y′ = sum(abs2, parent(y)) .* Δx .- 2real(y * Δx) .* parent(y)
    ∂y = y isa Transpose ? transpose(∂y′) : adjoint(∂y′)
    return y, ∂y
end

function frule(
    (_, ẋ), ::typeof(pinv), x::LinearAlgebra.AdjOrTransAbsVec{T}, tol::Real=0
) where {T<:Union{Real,Complex}}
    y = pinv(x, tol)
    Δx = unthunk(ẋ)
    ∂y = sum(abs2, y) .* vec(Δx') .- 2real(Δx * y) .* y
    return y, ∂y
end

# Formula for derivative adapted from Eq 4.12 of
# Golub, Gene H., and Victor Pereyra. "The Differentiation of Pseudo-Inverses and Nonlinear
# Least Squares Problems Whose Variables Separate."
# SIAM Journal on Numerical Analysis 10(2). (1973). 413-432. doi: 10.1137/0710036
function frule((_, Ȧ), ::typeof(pinv), A::AbstractMatrix{T}; kwargs...) where {T}
    Y = pinv(A; kwargs...)
    m, n = size(A)
    # contract over the largest dimension
    ΔA = unthunk(Ȧ)
    if m ≤ n
        ∂Y = -Y * (ΔA * Y)
        ∂Y = add!!(∂Y, (ΔA' - Y * (A * ΔA')) * (Y' * Y))  # (I - Y A) ΔA' Y' Y
        ∂Y = add!!(∂Y, Y * (Y' * ΔA') * (I - A * Y))  # Y Y' ΔA' (I - A Y)
    else
        ∂Y = -(Y * ΔA) * Y
        ∂Y = add!!(∂Y, (I - Y * A) * (ΔA' * Y') * Y)  # (I - Y A) ΔA' Y' Y
        ∂Y = add!!(∂Y, (Y * Y') * (ΔA' - (ΔA' * A) * Y))  # Y Y' ΔA' (I - A Y)
    end
    return Y, ∂Y
end

function rrule(
    ::typeof(pinv), x::Union{AbstractVector{T},LinearAlgebra.AdjOrTransAbsVec{T}}
) where {T<:Union{Real,Complex}}
    y, full_pb = rrule(pinv, x, 0)
    pinv_pullback(Δy) = return full_pb(Δy)[1:2]
    return y, pinv_pullback
end

function rrule(
    ::typeof(pinv), x::AbstractVector{T}, tol::Real
) where {T<:Union{Real,Complex}}
    y = pinv(x, tol)
    function pinv_pullback(ȳ)
        Δy = unthunk(ȳ)
        ∂x = sum(abs2, parent(y)) .* vec(Δy') .- 2real(y * Δy') .* parent(y)
        return (NoTangent(), ∂x, ZeroTangent())
    end
    return y, pinv_pullback
end

function rrule(
    ::typeof(pinv), x::LinearAlgebra.AdjOrTransAbsVec{T}, tol::Real
) where {T<:Union{Real,Complex}}
    y = pinv(x, tol)
    function pinv_pullback(ȳ)
        Δy = ȳ
        ∂x′ = sum(abs2, y) .* Δy .- 2real(y' * Δy) .* y
        ∂x = x isa Transpose ? transpose(conj(∂x′)) : adjoint(∂x′)
        return (NoTangent(), ∂x, ZeroTangent())
    end
    return y, pinv_pullback
end

function rrule(::typeof(pinv), A::AbstractMatrix{T}; kwargs...) where {T}
    Y = pinv(A; kwargs...)
    function pinv_pullback(Ȳ)
        ΔY = unthunk(Ȳ)
        m, n = size(A)
        # contract over the largest dimension
        if m ≤ n
            ∂A = (Y' * -ΔY) * Y'
            ∂A = add!!(∂A, (Y' * Y) * (ΔY' - (ΔY' * Y) * A)) # Y' Y ΔY' (I - Y A)
            ∂A = add!!(∂A, (I - A * Y) * (ΔY' * Y) * Y') # (I - A Y) ΔY' Y Y'
        elseif m > n
            ∂A = Y' * (-ΔY * Y')
            ∂A = add!!(∂A, Y' * (Y * ΔY') * (I - Y * A)) # Y' Y ΔY' (I - Y A)
            ∂A = add!!(∂A, (ΔY' - A * (Y * ΔY')) * (Y * Y')) # (I - A Y) ΔY' Y Y'
        end
        return (NoTangent(), ∂A)
    end
    return Y, pinv_pullback
end

#####
##### `sylvester`
#####

# included because the primal uses `schur`, for which we don't have a rule

function frule(
    (_, ΔA, ΔB, ΔC),
    ::typeof(sylvester),
    A::StridedMatrix{T},
    B::StridedMatrix{T},
    C::StridedMatrix{T},
) where {T<:BlasFloat}
    RA, QA = schur(A)
    RB, QB = schur(B)
    D = QA' * (C * QB)
    Y, scale = LAPACK.trsyl!('N', 'N', RA, RB, D)
    Ω = rmul!(QA * (Y * QB'), -inv(scale))
    ∂D = QA' * (mul!(muladd(ΔA, Ω, ΔC), Ω, ΔB, true, true) * QB)
    ∂Y, scale2 = LAPACK.trsyl!('N', 'N', RA, RB, ∂D)
    ∂Ω = rmul!(QA * (∂Y * QB'), -inv(scale2))
    return Ω, ∂Ω
end

# included because the primal mutates and uses `schur` and LAPACK

function rrule(
    ::typeof(sylvester), A::StridedMatrix{T}, B::StridedMatrix{T}, C::StridedMatrix{T}
) where {T<:BlasFloat}
    RA, QA = schur(A)
    RB, QB = schur(B)
    D = QA' * (C * QB)
    Y, scale = LAPACK.trsyl!('N', 'N', RA, RB, D)
    Ω = rmul!(QA * (Y * QB'), -inv(scale))
    function sylvester_pullback(ΔΩ)
        ∂Ω = T <: Real ? real(ΔΩ) : ΔΩ
        ∂Y = QA' * (∂Ω * QB)
        trans = T <: Complex ? 'C' : 'T'
        ∂D, scale2 = LAPACK.trsyl!(trans, trans, RA, RB, ∂Y)
        ∂Z = rmul!(QA * (∂D * QB'), -inv(scale2))
        return NoTangent(), @thunk(∂Z * Ω'), @thunk(Ω' * ∂Z), @thunk(∂Z * inv(scale))
    end
    return Ω, sylvester_pullback
end

#####
##### `lyap`
#####

# included because the primal uses `schur`, for which we don't have a rule

function frule(
    (_, Ȧ, Ċ), ::typeof(lyap), A::StridedMatrix{T}, C::StridedMatrix{T}
) where {T<:BlasFloat}
    ΔA = unthunk(Ȧ)
    ΔC = unthunk(Ċ)
    R, Q = schur(A)
    D = Q' * (C * Q)
    Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
    Ω = rmul!(Q * (Y * Q'), -inv(scale))
    ∂D = Q' * (mul!(muladd(ΔA, Ω, ΔC), Ω, ΔA', true, true) * Q)
    ∂Y, scale2 = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, ∂D)
    ∂Ω = rmul!(Q * (∂Y * Q'), -inv(scale2))
    return Ω, ∂Ω
end

# included because the primal mutates and uses `schur` and LAPACK

function rrule(
    ::typeof(lyap), A::StridedMatrix{T}, C::StridedMatrix{T}
) where {T<:BlasFloat}
    R, Q = schur(A)
    D = Q' * (C * Q)
    Y, scale = LAPACK.trsyl!('N', T <: Complex ? 'C' : 'T', R, R, D)
    Ω = rmul!(Q * (Y * Q'), -inv(scale))
    function lyap_pullback(ΔΩ)
        ∂Ω = T <: Real ? real(ΔΩ) : ΔΩ
        ∂Y = Q' * (∂Ω * Q)
        ∂D, scale2 = LAPACK.trsyl!(T <: Complex ? 'C' : 'T', 'N', R, R, ∂Y)
        ∂Z = rmul!(Q * (∂D * Q'), -inv(scale2))
        return NoTangent(),
        @thunk(mul!(∂Z * Ω', ∂Z', Ω, true, true)),
        @thunk(∂Z * inv(scale))
    end
    return Ω, lyap_pullback
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RegularGrid(dims, origin, spacing)

A regular grid with dimensions `dims`, lower left corner at `origin`
and cell spacing `spacing`. The three arguments must have the same length.

    RegularGrid(start, finish, dims=dims)

Alternatively, construct a regular grid from a `start` point (lower left)
to a `finish` point (upper right).

    RegularGrid{T}(dims)
    RegularGrid{T}(dim1, dim2, ...)

Finally, a regular grid can be constructed by only passing the dimensions
`dims` as a tuple, or by passing each dimension `dim1`, `dim2`, ... separately.
In this case, the origin and spacing default to (0,0,...) and (1,1,...).

## Examples

Create a 3D grid with 100x100x50 locations:

```julia
julia> RegularGrid(100,100,50)
```

Create a 2D grid with 100x100 locations and origin at (10.,20.) units:

```julia
julia> RegularGrid((100,100),(10.,20.),(1.,1.))
```

Create a 1D grid from -1 to 1 with 100 locations:

```julia
julia> RegularGrid((-1.,),(1.,), dims=(100,))
```
"""
struct RegularGrid{T,N} <: SpatialDomain{T,N}
  dims::Dims{N}
  origin::SVector{N,T}
  spacing::SVector{N,T}

  function RegularGrid{T,N}(dims, origin, spacing) where {N,T}
    @assert all(dims .> 0) "dimensions must be positive"
    @assert all(spacing .> 0) "spacing must be positive"
    new(dims, origin, spacing)
  end
end

RegularGrid(dims::Dims{N}, origin::Vec{N,T}, spacing::Vec{N,T}) where {N,T} =
  RegularGrid{T,N}(dims, origin, spacing)

RegularGrid(dims::Dims{N}, origin::NTuple{N,T}, spacing::NTuple{N,T}) where {N,T} =
  RegularGrid{T,N}(dims, origin, spacing)

RegularGrid(start::Vec{N,T}, finish::Vec{N,T};
            dims::Dims{N}=ntuple(i->100, N)) where {N,T} =
  RegularGrid{T,N}(dims, start, @. (finish - start) / (dims - 1))

RegularGrid(start::NTuple{N,T}, finish::NTuple{N,T};
            dims::Dims{N}=ntuple(i->100, N)) where {N,T} =
  RegularGrid(SVector(start), SVector(finish); dims=dims)

RegularGrid{T}(dims::Dims{N}) where {N,T} =
  RegularGrid{T,N}(dims, ntuple(i->zero(T), N), ntuple(i->one(T), N))

RegularGrid{T}(dims::Vararg{Int,N}) where {N,T} = RegularGrid{T}(dims)

RegularGrid(dims::Dims{N}) where {N} = RegularGrid{Float64}(dims)

RegularGrid(dims::Vararg{Int,N}) where {N} = RegularGrid{Float64}(dims)

Base.size(grid::RegularGrid) = grid.dims
origin(grid::RegularGrid)    = grid.origin
spacing(grid::RegularGrid)   = grid.spacing

nelms(grid::RegularGrid) = prod(grid.dims)

function coordinates!(buff::AbstractVector{T}, grid::RegularGrid{T,N},
                      location::Int) where {N,T}
  intcoords = CartesianIndices(grid.dims)[location]
  @inbounds @simd for i in 1:N
    buff[i] = grid.origin[i] + (intcoords[i] - 1)*grid.spacing[i]
  end
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, grid::RegularGrid{T,N}) where {N,T}
  dims = join(grid.dims, "×")
  print(io, "$dims RegularGrid{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", grid::RegularGrid{T,N}) where {N,T}
  println(io, grid)
  println(io, "  origin:  ", Tuple(grid.origin))
  print(  io, "  spacing: ", Tuple(grid.spacing))
end


text:
# Model structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

const IntTuple = Union{Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer},Array{Int64,1},Array{Int32,1}}
const RealTuple = Union{Tuple{Real,Real}, Tuple{Real,Real,Real},Array{Float64,1},Array{Float32,1}}

export Model, PhysicalParameter, get_dt

###################################################################################################
# PhysicalParameter abstract vector

"""
    PhysicalParameter
        n::IntTuple
        d::RealTuple
        o::RealTuple
        data::Union{Array, Number}

PhysicalParameter structure for physical space parameter.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`data`: the array of the parameter values of size n

Constructor
===========

A `PhysicalParameter` can be constructed in various ways but always require the origin `o` and grid spacing `d` that
cannot be infered from the array.

    PhysicalParameter(v::Array{vDT}, d, o) where `v` is an n-dimensional array and n=size(v)

    PhysicalParameter(n, d, o; vDT=Float32) Creates a zero PhysicalParameter

    PhysicalParameter(v::Array{vDT}, A::PhysicalParameter) Creates a PhysicalParameter from the Array `v` with n, d, o from `A`

    PhysicalParameter(v::Array{vDT, N}, n::Tuple, d::Tuple, o::Tuple) where `v` is a vector or nd-array that is reshaped into shape `n`

    PhysicalParameter(v::vDT, n::Tuple, d::Tuple, o::Tuple) Creates a constant (single number) PhyicalParameter

"""
mutable struct PhysicalParameter{vDT} <: AbstractVector{vDT}
    n::Tuple
    d::Tuple
    o::Tuple
    data::Union{Array, vDT}
end

mutable struct PhysicalParameterException <: Exception
    msg :: String
end

PhysicalParameter(v::BitArray{N}, args...) where N = v

function PhysicalParameter(v::Array{vDT}, d::Tuple, o::Tuple) where {vDT}
    length(size(v)) != length(o) && throw(PhysicalParameterException("Input array should be $(length(o))-dimensional"))
    return PhysicalParameter{vDT}(size(v), d, o, v)
end

function PhysicalParameter(n::Tuple, d::Tuple, o::Tuple; vDT=Float32)
    return PhysicalParameter{vDT}(n, d, o, zeros(vDT, n))
end

function PhysicalParameter(v::Array{vDT}, A::PhysicalParameter) where {vDT}
    return PhysicalParameter{vDT}(A.n, A.d, A.o, v)
end

function PhysicalParameter(v::Array{vDT, N}, n::Tuple, d::Tuple, o::Tuple) where {vDT, N}
    length(v) != prod(n) && throw(PhysicalParameterException("Incompatible number of element in input $(length(v)) with n=$(n)"))
    N == 1 && (v = reshape(v, n))
    return PhysicalParameter{vDT}(n, d, o, v)
end

function PhysicalParameter(v::vDT, n::Tuple, d::Tuple, o::Tuple) where vDT<:Number
    return PhysicalParameter{vDT}(n, d, o, v)
end

PhysicalParameter(p::PhysicalParameter, n::Tuple, d::Tuple, o::Tuple) = p
PhysicalParameter(p::PhysicalParameter) = p

# transpose and such.....
conj(x::PhysicalParameter{vDT}) where vDT = x
transpose(x::PhysicalParameter{vDT}) where vDT = PhysicalParameter{vDT}(x.n[end:-1:1], x.d[end:-1:1], x.o[end:-1:1], permutedims(x.data, length(x.n):-1:1))
adjoint(x::PhysicalParameter{vDT}) where vDT = transpose(x)

# Basic overloads
size(A::PhysicalParameter) = (prod(A.n), 1)
length(A::PhysicalParameter) = prod(A.n)

function norm(A::PhysicalParameter, order::Real=2)
    return norm(vec(A.data), order)
end

dot(A::PhysicalParameter, B::PhysicalParameter) = dot(vec(A.data), vec(B.data))
dot(A::PhysicalParameter, B::Array) = dot(vec(A.data), vec(B))
dot(A::Array, B::PhysicalParameter) = dot(vec(A), vec(B.data))

display(A::PhysicalParameter) = println("$(typeof(A)) of size $(A.n) with origin $(A.o) and spacing $(A.d)")
show(A::PhysicalParameter) = show(A.data)
showarg(io::IO, A::PhysicalParameter, toplevel) = print(io, typeof(A), " with size $(A.n), spacing $(A.d) and origin $(A.o)")

# Indexing
firstindex(A::PhysicalParameter) = 1
lastindex(A::PhysicalParameter) = length(A)
lastindex(A::PhysicalParameter, dim::Int) = A.n[dim]

function promote_shape(p::PhysicalParameter, A::Array{vDT, N}) where {vDT, N}
    (size(A) != p.n && N>1) && return promote_shape(p.data, A)
    (length(A) == prod(p.n) && N==1) && return size(A)
    return promote_shape(A, A)
end

promote_shape(A::Array{vDT, N}, p::PhysicalParameter) where {vDT, N} = promote_shape(p, A)
reshape(p::PhysicalParameter, n::Tuple{Vararg{Int64,N}}) where N = (prod(n)==prod(p.n) && return p)

dotview(A::PhysicalParameter{vDT}, I::Vararg{Union{Function, Int, UnitRange{Int}}, N}) where {vDT, N} = dotview(A.data, I...)
Base.dotview(m::PhysicalParameter, i) = Base.dotview(m.data, i)

getindex(A::PhysicalParameter, i::Int) = A.data[i]
getindex(A::PhysicalParameter, i::Colon) = A.data[:]

function getindex(A::PhysicalParameter{T}, I::Vararg{Union{Int, BitArray, Function, StepRange{Int}, UnitRange{Int}}, N}) where {N, T}
    new_v = getindex(A.data, I...)
    length(size(new_v)) != length(A.n) && (return new_v)
    s = [i == (:) ? 0 : i[1]-1 for i=I]
    st = [getattr(i, :step, 1) for i=I]
    new_o = [ao+i*d for (ao, i, d)=zip(A.o, s, A.d)]
    new_d = [d*s for (d, s)=zip(A.d, st)]
    PhysicalParameter{T}(size(new_v), tuple(new_d...), tuple(new_o...), new_v)
end

setindex!(A::PhysicalParameter, v, I::Vararg{Union{Int, Function, UnitRange{Int}}, N}) where {N} = setindex!(A.data, v, I...)
setindex!(A::PhysicalParameter, v, i::Int) = (A.data[i] = v)

# Constructiors by copy
similar(x::PhysicalParameter{vDT}) where {vDT} = vDT(0) .* x
copy(x::PhysicalParameter{vDT}) where {vDT} = PhysicalParameter{vDT}(x.n, x.d, x.o, x.data)

# Equality
isequal(A::PhysicalParameter, B::PhysicalParameter) = (A.data == B.data && A.o == B.o && A.d == B.d)
isapprox(A::PhysicalParameter, B::PhysicalParameter; kwargs...) = (isapprox(A.data, B.data) && A.o == B.o && A.d == B.d)
isapprox(A::PhysicalParameter, B::AbstractArray; kwargs...) = isapprox(A.data, B)
isapprox(A::AbstractArray, B::PhysicalParameter; kwargs...) = isapprox(A, B.data)

# # Arithmetic operations
function compare(A::PhysicalParameter, B::PhysicalParameter)
    A.o != B.o && throw(PhysicalParameterException("Incompatible origins $(A.o) and $(B.o)"))
    A.d != B.d && throw(PhysicalParameterException("Incompatible spacing $(A.d) and $(B.d)"))
    A.n != B.n && throw(PhysicalParameterException("Incompatible sizes $(A.n) and $(B.n)"))
end

function +(A::PhysicalParameter{vDT}, B::PhysicalParameter{vDT}) where {vDT}
    compare(A, B)
    return PhysicalParameter(A.data + B.data, A)
end

function -(A::PhysicalParameter{vDT}, B::PhysicalParameter{vDT}) where {vDT}
    compare(A, B)
    return PhysicalParameter(A.data - B.data, A)
end

function *(A::PhysicalParameter{vDT}, B::PhysicalParameter{vDT}) where {vDT}
    compare(A, B)
    return PhysicalParameter(A.data .* B.data, A)
end

function /(A::PhysicalParameter{vDT}, B::PhysicalParameter{vDT}) where {vDT}
    compare(A, B)
    return PhysicalParameter(A.data ./ B.data, A)
end

# Brodacsting
BroadcastStyle(::Type{<:PhysicalParameter}) = Broadcast.ArrayStyle{PhysicalParameter}()

function similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{PhysicalParameter}}, ::Type{ElType}) where ElType
    # Scan the inputs for the ArrayAndChar:
    A = find_aac(bc)
    # Use the char field of A to create the output
    ElType == Bool ? Ad = similar(Array{ElType}, axes(A.data)) : Ad = similar(A.data)
    PhysicalParameter(Ad, A.d, A.o)
end

"`A = find_aac(As)` returns the first PhysicalParameter among the arguments."
find_aac(bc::Base.Broadcast.Broadcasted) = find_aac(bc.args)
find_aac(args::Tuple) = find_aac(find_aac(args[1]), Base.tail(args))
find_aac(x) = x
find_aac(::Tuple{}) = nothing
find_aac(a::PhysicalParameter, rest) = a
find_aac(::Any, rest) = find_aac(rest)

function broadcasted(f::Function, A::AbstractArray{vDT, N}, p::PhysicalParameter{RDT}) where {vDT, RDT, N}
    if size(A) != p.n
        if length(A) == length(p)
            return PhysicalParameter(materialize(broadcasted(f, A, vec(p.data))), p.n, p.d, p.o)
        end
        throw(PhysicalParameterException("Incompatible sizes $(size(A)) and $(p.n)"))
    end
    return PhysicalParameter(materialize(broadcasted(f, A, p.data)), p.d, p.o)
end

function broadcasted(f::Function, p::PhysicalParameter{RDT}, A::AbstractArray{vDT, N}) where {vDT, RDT, N}
    if size(A) != p.n
        if length(A) == length(p)
            return PhysicalParameter(materialize(broadcasted(f, vec(p.data), A)), p.n, p.d, p.o)
        end
        throw(PhysicalParameterException("Incompatible sizes $(size(A)) and $(p.n)"))
    end
    return PhysicalParameter(materialize(broadcasted(f, p.data, A)), p.d, p.o)
end

function broadcasted(f::Function, p1::PhysicalParameter{RDT}, p2::PhysicalParameter{RDT}) where {RDT}
    p1.n != p2.n && throw(PhysicalParameterException("Incompatible sizes $(p1.n) and $(p2.n)"))
    return PhysicalParameter(materialize(broadcasted(f, p1.data, p2.data)), p1.d, p1.o)
end

broadcasted(f::Function, p::PhysicalParameter, bc::Base.Broadcast.Broadcasted) = broadcasted(f, p, materialize(bc))
broadcasted(f::Function, bc::Base.Broadcast.Broadcasted, p::PhysicalParameter) = broadcasted(f, materialize(bc), p)

function *(A::Union{joMatrix, joLinearFunction, joLinearOperator, joCoreBlock}, p::PhysicalParameter{RDT}) where {RDT}
    @warn "JOLI linear operator, returning julia Array"
    return A*vec(p.data)
end

materialize!(p::PhysicalParameter{RDT}, b::Array{<:Number, N}) where{RDT, N} = (p.data[:] .= b)
materialize!(p::PhysicalParameter{RDT}, b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{N}}) where{RDT, N} = materialize!(p, collect(b))
materialize!(p::Array, b::Broadcast.Broadcasted{Broadcast.ArrayStyle{PhysicalParameter}}) = materialize!(p, reshape(collect(b), size(p)))

# Linalg Extras
mul!(x::PhysicalParameter, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::Array) = mul!(view(x.data, :), F, y)
mul!(x::PhysicalParameter, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::PhysicalParameter) = mul!(view(x.data, :), F, view(y.data, :))
mul!(x::Array, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::PhysicalParameter) = mul!(x, F, y[1:end])

# For ploting
array2py(p::PhysicalParameter{vDT}, i::Int64, I::CartesianIndex{N}) where {vDT, N} = array2py(p.data, i, I)

###################################################################################################
# Isotropic acoustic

"""
    Model
        n::IntTuple
        d::RealTuple
        o::RealTuple
        nb::Integer
        params::Dict
        rho::Array

Model structure for seismic velocity models.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`nb`: number of absorbing boundary points in each direction

`params`: Physical parameters such has squared slowness, denisty or THomesne parameters

Constructor
===========

The parameters `n`, `d`, `o` and `m` are mandatory, whith `nb` and other physical parameters being optional input arguments.

    Model(n, d, o, m; nb=40, rho=1, epsilon=0, delta=0, theta=0, phi=0)

where

`m`: velocity model in slowness squared (s^2/km^2)

`epsilon`: Epsilon thomsen parameter ( between -1 and 1)

`delta`: Delta thomsen parameter ( between -1 and 1 and delta < epsilon)

`theta`: Anisotopy dip in radian

`phi`: Anisotropy asymuth in radian

`rho`: density (g / m^3)

"""
mutable struct Model
    n::IntTuple
    d::RealTuple
    o::RealTuple
    nb::Integer # number of absorbing boundaries points on each side
    params::Dict
end

###################################################################################################
# Constructors

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m;
               epsilon=nothing, delta=nothing, theta=nothing,
               phi=nothing, rho=nothing, nb=40)

    params = Dict(:m => PhysicalParameter(Float32.(m), d, o))
    for (name, val)=zip([:rho, :epsilon, :delta, :theta, :phi], [rho, epsilon, delta, theta, phi])
        ~isnothing(val) && (params[name] = PhysicalParameter(Float32.(val), n, d, o))
    end

    return Model(n, d, o, nb, params)
end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho; nb=40)
    return Model(n, d, o, m; rho=rho, nb=nb)
end

get_dt(m::Model; dt=nothing) = calculate_dt(m; dt=dt)
getindex(m::Model, sym::Symbol) = m.params[sym]

Base.setproperty!(m::Model, s::Symbol, p::PhysicalParameter{Float32}) = (m.params[s] = p)

function Base.getproperty(obj::Model, sym::Symbol)
    if sym == :params
        return getfield(obj, sym)
    elseif sym in keys(obj.params)
        return obj.params[sym]
    else
         # fallback to getfield
        return getfield(obj, sym)
    end
end

similar(x::PhysicalParameter{vDT}, m::Model) where {vDT} = PhysicalParameter(m.n, m.d, m.o; vDT=vDT)
similar(x::Array, m::Model) where {vDT} = similar(x, m.n)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
"""
	function [H,LVpen] = histoptimal(Y, Dmax)
 PERFORMS DENSITY ESTIMATION using PENALIZED REGULAR HISTOGRAM, author Yves ROZENHOLC
 Compute automaticaly the number of bins of the classical density estimator,
 the regular histogram, and return description of the obtained histogram(s). The
 number(s) of bins chose by the algorithm derived from the works of
 MODEL and COMMANDE LINE
 MODEL :
 Y n-sample of an UNKNOWN density f
 COMMANDE LINE :
 H = histoptimal(Y, Dmax, withfig)
 [H,LVpen] = histoptimal(Y, Dmax, withfig)
 INPUT
      Y, the observations - an n sample of an unknown density f
      Dmax, the maximum number of bins of the histogram
           Default value [n/log(n)].
           For very large n, use smaller value to reduce the computation time.
      withfig, if 1 plot the selected histogram
           Default 1.
 OUTPUT
      H, the selected histogram
      LVpen, the penalized Log-Likelihood of the selected histogram
 CONNECTED PAPERS
 HOW MANY BINS SHOULD BE PUT IN A REGULAR HISTOGRAM.
 by L. BIRG\'E and Y. ROZENHOLC in ESAIM P&S
 Other references :
 	A. Barron, L. Birg'e, P.Massart 1995
 	G. Castellan 1998
 for details you can contact
 	Yves ROZENHOLC
 	yves.rozenholc@univ-paris5.fr
 version du 17/03/00
"""
function histoptimal(Y, Dmax::Int = floor(Int, length(Y)/log(length(Y))))
	Y = Y[:];
	n = length(Y)
	D = 2:Dmax;
	pen = (D .- 1) .+ log.(D) .^ 2.5;

	mini, maxi =minimum(Y), maximum(Y);
	mini ≈ maxi && @error "almost zero range of data"
	if (mini<0)||(maxi>1)
		Y = (Y .- mini)./(maxi - mini);
	else
   mini, maxi = 0, 1;
	end;

	LV = zeros(length(D)); #Le vecteur des max de la log-vraisemblance
	for (i, d) = enumerate(D)
		Nd = bindata(Y, d); # valeur de la fct de repartition aux bornes
		ind = findall(Nd .!=0 );
    LV[i] = sum(Nd[ind] .* log.(Nd[ind] .* d ./ n));
	end
	ii = argmax(LV .- pen)
	LV[ii] - pen[ii], D[ii]
end

function bindata(x, k::Int)
	δ = (maximum(x) - minimum(x)) / k
	(δ == 0) && @error "cannot bin interval of zero length"
	o = zeros(k);
	bindata!(o, x, δ)
end

function bindata(x, δ::AbstractFloat)
	o = zeros(ceil(Int,(maximum(x) - minimum(x)) / δ))
	bindata!(o, x, δ)
end

function bindata!(o, x, δ::AbstractFloat)
	foreach(z -> o[min(length(o),floor(Int, z/δ) + 1)] += 1, x)
	o
end

struct EqualHistogram{T}
	p :: Vector{T}
	mi::T
	mx::T
	δ::T
end

function EqualHistogram(x)
	mi, mx = minimum(x), maximum(x)
	mx - mi == 0 && @error("The span of histogram should be greater than zero")

	x = (x .- mi) ./ (mx - mi)

	k = histoptimal(x)[2]
	δ = 1 / k
	p = bindata(x, δ) ./ (length(x)*(mx - mi))
	EqualHistogram(p, mi, mx, δ)
end

function (m::EqualHistogram)(x)
	ii = floor.(Int, (x .- m.mi) ./ ((m.mx - m.mi) * m.δ))
	o = similar(x)
	fill!(o, log(10*eps(eltype(x))))
	mask = (ii .>0 ) .& (ii .< length(m.p))
	o[mask] .= -log.( m.p[ii[mask]] .+ 10*eps(eltype(x)))
	o
end

struct Loda{V,E}
	ws::Vector{V}
	histograms::Vector{E}
end

criterion(fs) = sum(abs.(fs[end] .- fs[end - 1])) / sum(abs.(fs[2] .- fs[1]))

function Loda(x, τ::AbstractFloat = 0.01, sparse::Bool = false)
	sparse == true && @error("Not implemented")
	Loda(optimize_loda_params(x, τ)[1:2]...)
end

function optimize_loda_params(x,τ)
	d = size(x,1)
	ws = []
	fs = []
	histograms = []
	k = 1
	cumf = zeros(size(x,2))
	while true
		w = randn(1, d)
		p = w * x
		h = EqualHistogram(p)
		cumf .+= h(p)[:]
		push!(fs, cumf ./ k)
		push!(ws, w)
		push!(histograms, h)

		length(fs) > 2 && criterion(fs) < τ && break
		k += 1
	end
	(ws, histograms, fs)
end

function (m::Loda)(x)
	o = similar(x, size(x,2))
	o .= 0
	for i in 1:length(m.ws)
		h, w = m.histograms[i], m.ws[i]
		o .-= h(w * x)[:]
	end
	o ./= length(m.ws)
end


text:
#
# vectors.jl -
#
# Implement basic operations for *vectors*.  Here arrays of any rank are
# considered as *vectors*, the only requirements are that, when combining
# *vectors*, they have the same list of axes (i.e. the same dimensions for
# most arrays).  These methods are intended to be used for numerical
# optimization.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

"""
    vnorm2([T,] v)

yields the Euclidean (L2) norm of `v`.  The floating point type of the result
can be imposed by optional argument `T`.  Also see [`vnorm1`](@ref) and
[`vnorminf`](@ref).

"""
function vnorm2(v::AbstractArray{<:Floats})
    s = zero(real(eltype(v)))
    @inbounds @simd for i in eachindex(v)
        s += abs2(v[i])
    end
    return sqrt(s)
end

"""
    vnorm1([T,] v)

yields the L1 norm of `v`, that is the sum of the absolute values of its
elements.  The floating point type of the result can be imposed by optional
argument `T`.  For a complex valued argument, the result is the sum of the
absolute values of the real part and of the imaginary part of the elements
(like BLAS `asum`).

See also [`vnorm2`](@ref) and [`vnorminf`](@ref).

"""
function vnorm1(v::AbstractArray{<:Reals})
    s = zero(real(eltype(v)))
    @inbounds @simd for i in eachindex(v)
        s += abs(v[i])
    end
    return s
end

function vnorm1(v::AbstractArray{<:Complexes})
    si = sr = zero(real(eltype(v)))
    @inbounds @simd for i in eachindex(v)
        z = v[i]
        sr += abs(real(z))
        si += abs(imag(z))
    end
    return sr + si
end

"""
    vnorminf([T,] v)

yields the infinite norm of `v`, that is the maximum absolute value of its
elements.  The floating point type of the result can be imposed by optional
argument `T`.  Also see [`vnorm1`](@ref) and [`vnorm2`](@ref).

"""
function vnorminf(v::AbstractArray{<:Reals})
    absmax = zero(real(eltype(v)))
    @inbounds @simd for i in eachindex(v)
        absmax = max(absmax, abs(v[i]))
    end
    return absmax
end

function vnorminf(v::AbstractArray{<:Complexes})
    abs2max = zero(real(eltype(v)))
    @inbounds @simd for i in eachindex(v)
        abs2max = max(abs2max, abs2(v[i]))
    end
    return sqrt(abs2max)
end

# Versions with forced type of output result.
for func in (:vnorm2, :vnorm1, :vnorminf)
    @eval $func(::Type{T}, v) where {T<:AbstractFloat} =
        convert(T, $func(v))::T
end

#------------------------------------------------------------------------------

"""
    vcreate(x)

yields a new variable instance similar to `x`.  If `x` is an array, the
element type of the result is a floating-point type.

Also see [`similar`](@ref).

"""
vcreate(x::AbstractArray{T}) where {R<:Real,T<:Union{R,Complex{R}}} =
    similar(x, float(T))

#------------------------------------------------------------------------------

"""
    vcopy!(dst, src) -> dst

copies the contents of `src` into `dst` and returns `dst`.  This function
checks that the copy makes sense (for instance, for array arguments, the
`copyto!` operation does not check that the source and destination have the
same dimensions).

Also see [`copyto!`](@ref), [`vcopy`](@ref), [`vswap!`](@ref).

"""
function vcopy!(dst::AbstractArray{<:Real,N},
                src::AbstractArray{<:Real,N}) where {N}
    if dst !== src
        axes(dst) == axes(src) || arguments_have_incompatible_axes()
        copyto!(dst, src)
    end
    return dst
end

function vcopy!(dst::AbstractArray{<:Complex{<:Real},N},
                src::AbstractArray{<:Complex{<:Real},N}) where {N}
    if dst !== src
        axes(dst) == axes(src) || arguments_have_incompatible_axes()
        copyto!(dst, src)
    end
    return dst
end

"""
    vcopy(x)

yields a fresh copy of the *vector* `x`.  If `x` is is an array, the element
type of the result is a floating-point type.

Also see [`copy`](@ref), [`vcopy!`](@ref), [`vcreate!`](@ref).

"""
vcopy(x) = vcopy!(vcreate(x), x)

"""
    vswap!(x, y)

exchanges the contents of `x` and `y` (which must have the same element type
and axes if they are arrays).

Also see [`vcopy!`](@ref).

"""
vswap!(x::T, y::T) where {T<:AbstractArray} =
    x === y || _swap!(x, y)

vswap!(x::AbstractArray{T,N}, y::AbstractArray{T,N}) where {T,N} =
    _swap!(x, y)

# Forced swapping.
_swap!(x::AbstractArray{T,N}, y::AbstractArray{T,N}) where {T,N} =
    @inbounds @simd for i in all_indices(x, y)
        temp = x[i]
        x[i] = y[i]
        y[i] = temp
    end

#------------------------------------------------------------------------------

"""
    vfill!(x, α) -> x

sets all elements of `x` with the scalar value `α` and return `x`.  The default
implementation just calls `fill!(x,α)` but this method may be specialized for
specific types of variables `x`.

Also see [`vzero!`](@ref), [`fill!`](@ref).

"""
vfill!(x, α) = fill!(x, α)

"""
    vzero!(x) -> x

fills `x` with zeros and returns it.  The default implementation just calls
`fill!(x,zero(eltype(x)))` but this method may be specialized for specific
types of variables `x`.

Also see [`vfill!`](@ref).

"""
vzero!(x) = vfill!(x, zero(eltype(x)))

"""
    vzeros(x)

yields a *vector* like `x` filled with zeros.

Also see [`vones`](@ref), [`vcreate`](@ref), [`vfill!`](@ref).

"""
vzeros(x) = vzero!(vcreate(x))

"""
    vones(x)

yields a *vector* like `x` filled with ones.

Also see [`vzeros`](@ref), [`vcreate`](@ref), [`vfill!`](@ref).

"""
vones(x) = vfill!(vcreate(x), 1)
vones(x::AbstractArray{T}) where {T} = vfill!(vcreate(x), one(T))

#------------------------------------------------------------------------------

"""
    vscale!(dst, α, src) -> dst

overwrites `dst` with `α*src` and returns `dst`.  Computations are done at the
numerical precision of `promote_eltype(src,dst)`.  The source argument may be
omitted to perform *in-place* scaling:

    vscale!(x, α) -> x

overwrites `x` with `α*x` and returns `x`.  The convention is that the result
is zero-filled if `α=0` (whatever the values in the source).

Methods are provided by default so that the order of the factor `α` and the
source vector may be reversed:

    vscale!(dst, src, α) -> dst
    vscale!(α, x) -> x

Also see [`vscale`](@ref), [`LinearAlgebra.rmul!](@ref).

"""
function vscale!(dst::AbstractArray{<:Floats,N},
                 α::Number,
                 src::AbstractArray{<:Floats,N}) where {N}
    if α == 1
        vcopy!(dst, src)
    elseif α == 0
        axes(dst) == axes(src) || arguments_have_incompatible_axes()
        vzero!(dst)
    elseif α == -1
        @inbounds @simd for i in all_indices(dst, src)
            dst[i] = -src[i]
        end
    else
        alpha = promote_multiplier(α, src)
        @inbounds @simd for i in all_indices(dst, src)
            dst[i] = alpha*src[i]
        end
    end
    return dst
end

# In-place scaling.
function vscale!(x::AbstractArray{<:Floats}, α::Number)
    if α == 0
        vzero!(x)
    elseif α == -1
        @inbounds @simd for i in eachindex(x)
            x[i] = -x[i]
        end
    elseif α != 1
        alpha = promote_multiplier(α, x)
        @inbounds @simd for i in eachindex(x)
            x[i] *= alpha
        end
    end
    return x
end

# In-place scaling with reverse order of arguments.
vscale!(α::Number, x) = vscale!(x, α)

# Scaling for other *vector* types.
vscale!(x, α::Number) = vscale!(x, α, x)
vscale!(dst, src, α::Number) = vscale!(dst, α, src)

# The following methods are needed to avoid looping forever.
vscale!(::Number, ::Number) = error("bad argument types")
vscale!(::Any, ::Number, ::Number) = error("bad argument types")
vscale!(::Number, ::Any, ::Number) = error("bad argument types")
vscale!(::Number, ::Number, ::Any) = error("bad argument types")
vscale!(::Number, ::Number, ::Number) = error("bad argument types")

"""
    vscale(α, x)

or

    vscale(x, α)

yield a new *vector* whose elements are those of `x` multiplied by the scalar
`α`.

Also see [`vscale!`](@ref), [`vcreate`](@ref).

"""
vscale(α::Number, x) = vscale!(vcreate(x), α, x)
vscale(x, α::Number) = vscale(α, x)

#------------------------------------------------------------------------------
# ELEMENT-WISE MULTIPLICATION

"""
    vproduct(x, y) -> z

yields the element-wise multiplication of `x` by `y`.  To avoid allocating the
result, the destination array `dst` can be specified with the in-place version
of the method:

    vproduct!(dst, [sel,] x, y) -> dst

which overwrites `dst` with the elementwise multiplication of `x` by `y`.
Optional argument `sel` is a selection of indices to which apply the operation.

"""
vproduct(x::V, y::V) where {V} = vproduct!(vcreate(x), x, y)

vproduct(x::AbstractArray{<:Any,N}, y::AbstractArray{<:Any,N}) where {N} =
    vproduct!(similar(x, promote_eltype(x,y)), x, y)

for Td in (AbstractFloat, Complex{<:AbstractFloat}),
    Tx in (AbstractFloat, Complex{<:AbstractFloat}),
    Ty in (AbstractFloat, Complex{<:AbstractFloat})

    if Td <: Complex || (Tx <: Real && Ty <: Real)

        @eval function vproduct!(dst::AbstractArray{<:$Td,N},
                                 x::AbstractArray{<:$Tx,N},
                                 y::AbstractArray{<:$Ty,N}) where {N}
            @inbounds @simd for i in all_indices(dst, x, y)
                dst[i] = x[i]*y[i]
            end
            return dst
        end

        @eval function vproduct!(dst::AbstractArray{<:$Td,N},
                                 sel::AbstractVector{Int},
                                 x::AbstractArray{<:$Tx,N},
                                 y::AbstractArray{<:$Ty,N}) where {N}
            if checkselection(sel, dst, x, y)
                @inbounds @simd for j in eachindex(sel)
                    i = sel[j]
                    dst[i] = x[i]*y[i]
                end
            end
            return dst
        end

    end

end

@doc @doc(vproduct) vproduct!

#------------------------------------------------------------------------------
# VECTOR UPDATE

"""
    vupdate!(y, [sel,] α, x) -> y

overwrites `y` with `α*x + y` and returns `y`.  The code is optimized for some
specific values of the multiplier `α`.  For instance, if `α` is zero, then `y`
is left unchanged without using `x`.  Computations are performed at the
numerical precision of `promote_eltype(x,y)`.

Optional argument `sel` is a selection of indices to which apply the operation.
Note that if an index is repeated, the operation will be performed several
times at this location.

See also: [`vscale!`](@ref), [`vcombine!](@ref).

"""
function vupdate!(y::AbstractArray{<:Number,N},
                  α::Number,
                  x::AbstractArray{<:Number,N}) where {N}
    axes(x) == axes(y) || arguments_have_incompatible_axes()
    if α == 1
        @inbounds @simd for i in eachindex(x, y)
            y[i] += x[i]
        end
    elseif α == -1
        @inbounds @simd for i in eachindex(x, y)
            y[i] -= x[i]
        end
    elseif α != 0
        alpha = promote_multiplier(α, x)
        @inbounds @simd for i in eachindex(x, y)
            y[i] += alpha*x[i]
        end
    end
    return y
end

function vupdate!(y::AbstractArray{<:Floats,N},
                  sel::AbstractVector{Int},
                  α::Number,
                  x::AbstractArray{<:Floats,N}) where {N}
    if checkselection(sel, x, y)
        if α == 1
            @inbounds @simd for j in eachindex(sel)
                i = sel[j]
                y[i] += x[i]
            end
        elseif α == -1
            @inbounds @simd for j in eachindex(sel)
                i = sel[j]
                y[i] -= x[i]
            end
        elseif α != 0
            alpha = promote_multiplier(α, x)
            @inbounds @simd for j in eachindex(sel)
                i = sel[j]
                y[i] += alpha*x[i]
            end
        end
    end
    return y
end

#------------------------------------------------------------------------------
# LINEAR COMBINATION

"""
    vcombine(α, x, β, y) -> dst

yields the linear combination `dst = α*x + β*y`.

----

To avoid allocating the result, the destination array `dst` can be specified
with the in-place version of the method:

    vcombine!(dst, α, x, β, y) -> dst

The code is optimized for some specific values of the multipliers `α` and `β`.
For instance, if `α` (resp. `β`) is zero, then the prior contents of `x`
(resp. `y`) is not used.

The source(s) and the destination can be the same.  For instance, the two
following lines of code produce the same result:

    vcombine!(dst, 1, dst, α, x)
    vupdate!(dst, α, x)

See also: [`vscale!`](@ref), [`vupdate!](@ref).

"""
vcombine(α::Number, x::V, β::Number, y::V) where {V} =
    vcombine!(vcreate(x), α, x, β, y)

function vcombine!(dst::AbstractArray{<:Number,N},
                   α::Number, x::AbstractArray{<:Number,N},
                   β::Number, y::AbstractArray{<:Number,N}) where {N}
    axes(dst) == axes(x) == axes(y) || arguments_have_incompatible_axes()
    if α == 0
        if β == 0
            vzero!(dst)
        elseif β == 1
            _vcombine!(dst, axpby_yields_y,     0,x, 1,y)
        elseif β == -1
            _vcombine!(dst, axpby_yields_my,    0,x,-1,y)
        else
            b = promote_multiplier(β, y)
            _vcombine!(dst, axpby_yields_by,    0,x, b,y)
        end
    elseif α == 1
        if β == 0
            _vcombine!(dst, axpby_yields_x,     1,x, 0,y)
        elseif β == 1
            _vcombine!(dst, axpby_yields_xpy,   1,x, 1,y)
        elseif β == -1
            _vcombine!(dst, axpby_yields_xmy,   1,x,-1,y)
        else
            b = promote_multiplier(β, y)
            _vcombine!(dst, axpby_yields_xpby,  1,x, b,y)
        end
    elseif α == -1
        if β == 0
            _vcombine!(dst, axpby_yields_mx,   -1,x, 0,y)
        elseif β == 1
            _vcombine!(dst, axpby_yields_ymx,  -1,x, 1,y)
        elseif β == -1
            _vcombine!(dst, axpby_yields_mxmy, -1,x,-1,y)
        else
            b = promote_multiplier(β, y)
            _vcombine!(dst, axpby_yields_bymx, -1,x, b,y)
        end
    else
        a = promote_multiplier(α, x)
        if β == 0
            _vcombine!(dst, axpby_yields_ax,    a,x, 0,y)
        elseif β == 1
            _vcombine!(dst, axpby_yields_axpy,  a,x, 1,y)
        elseif β == -1
            _vcombine!(dst, axpby_yields_axmy,  a,x,-1,y)
        else
            b = promote_multiplier(β, y)
            _vcombine!(dst, axpby_yields_axpby, a,x, b,y)
        end
    end
    return dst
end

function _vcombine!(dst::AbstractArray{<:Number,N},
                    f::Function,
                    α::Number, x::AbstractArray{<:Number,N},
                    β::Number, y::AbstractArray{<:Number,N}) where {N}
    @inbounds @simd for i in eachindex(dst, x, y)
        dst[i] = f(α, x[i], β, y[i])
    end
end

@doc @doc(vcombine) vcombine!

#------------------------------------------------------------------------------
# INNER PRODUCT

"""
    vdot([T,] [w,] x, y)

yields the inner product of `x` and `y`; that is, the sum of `conj(x[i])*y[i]`
or, if `w` is specified, the sum of `w[i]*conj(x[i])*y[i]` (`w` must have
real-valued elements), for all indices `i`.  Optional argument `T` is the
floating point type of the result.

Another possibility is:

    vdot([T,] sel, x, y)

with `sel` a selection of indices to restrict the computation of the inner
product to some selected elements.  This yields the sum of `x[i]*y[i]` for all
`i ∈ sel`.

If the arguments have complex-valued elements and `T` is specified as a
floating-point type, complexes are considered as vectors of pairs of reals and
the result is:

    vdot(T::Type{AbstractFloat}, x, y)
    -> ((x[1].re*y[1].re + x[1].im*y[1].im) +
        (x[2].re*y[2].re + x[2].im*y[2].im) + ...)

"""
vdot(::Type{T}, x, y) where {T<:AbstractFloat} = convert(T,vdot(x,y))::T
vdot(::Type{T}, w, x, y) where {T<:AbstractFloat} = convert(T,vdot(w,x,y))::T

function vdot(x::AbstractArray{<:AbstractFloat,N},
              y::AbstractArray{<:AbstractFloat,N}) where {N}
    s = zero(promote_eltype(x, y))
    @inbounds @simd for i in all_indices(x, y)
        s += x[i]*y[i]
    end
    return s
end

function vdot(x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(promote_eltype(x, y))
    @inbounds @simd for i in all_indices(x, y)
        s += conj(x[i])*y[i]
    end
    return s
end

# This one yields the real part of the dot product, just as if complexes were
# pairs of reals.
function vdot(T::Type{<:AbstractFloat},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(real(promote_eltype(x, y)))
    @inbounds @simd for i in all_indices(x, y)
        xi = x[i]
        yi = y[i]
        s += real(xi)*real(yi) + imag(xi)*imag(yi)
    end
    return convert(T, s)::T
end

function vdot(T::Type{Complex{<:AbstractFloat}},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    return convert(T, vdot(x, y))::T
end

function vdot(w::AbstractArray{<:AbstractFloat,N},
              x::AbstractArray{<:AbstractFloat,N},
              y::AbstractArray{<:AbstractFloat,N}) where {N}
    s = zero(promote_eltype(w, x, y))
    @inbounds @simd for i in all_indices(w, x, y)
        s += w[i]*x[i]*y[i]
    end
    return s
end

function vdot(T::Type{<:AbstractFloat},
              w::AbstractArray{<:AbstractFloat,N},
              x::AbstractArray{<:AbstractFloat,N},
              y::AbstractArray{<:AbstractFloat,N}) where {N}
    return convert(T, vdot(w, x, y))::T
end

function vdot(w::AbstractArray{<:AbstractFloat,N},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(promote_eltype(w, x, y))
    @inbounds @simd for i in all_indices(w, x, y)
        s += w[i]*conj(x[i])*y[i]
    end
    return s
end

function vdot(T::Type{<:AbstractFloat},
              w::AbstractArray{<:AbstractFloat,N},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(real(promote_eltype(w, x, y)))
    @inbounds @simd for i in all_indices(w, x, y)
        xi = x[i]
        yi = y[i]
        s += (real(xi)*real(yi) + imag(xi)*imag(yi))*w[i]
    end
    return convert(T, s)::T
end

function vdot(T::Type{Complex{<:AbstractFloat}},
              w::AbstractArray{<:AbstractFloat,N},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    return convert(T, vdot(w, x, y))::T
end

function vdot(sel::AbstractVector{Int},
              x::AbstractArray{<:AbstractFloat,N},
              y::AbstractArray{<:AbstractFloat,N}) where {N}
    s = zero(promote_eltype(x, y))
    if checkselection(sel, x, y)
        @inbounds @simd for j in eachindex(sel)
            i = sel[j]
            s += x[i]*y[i]
        end
    end
    return s
end

function vdot(T::Type{<:AbstractFloat},
              sel::AbstractVector{Int},
              x::AbstractArray{<:AbstractFloat,N},
              y::AbstractArray{<:AbstractFloat,N}) where {N}
    return convert(T, vdot(sel, x, y))::T
end

function vdot(sel::AbstractVector{Int},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(promote_eltype(x, y))
    if checkselection(sel, x, y)
        @inbounds @simd for j in eachindex(sel)
            i = sel[j]
            s += conj(x[i])*y[i]
        end
    end
    return s
end

function vdot(T::Type{<:AbstractFloat},
              sel::AbstractVector{Int},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    s = zero(real(promote_eltype(x, y)))
    if checkselection(sel, x, y)
        @inbounds @simd for j in eachindex(sel)
            i = sel[j]
            xi = x[i]
            yi = y[i]
            s += real(xi)*real(yi) + imag(xi)*imag(yi)
        end
    end
    return convert(T, s)::T
end

function vdot(T::Type{Complex{<:AbstractFloat}},
              sel::AbstractVector{Int},
              x::AbstractArray{<:Complex{<:AbstractFloat},N},
              y::AbstractArray{<:Complex{<:AbstractFloat},N}) where {N}
    return convert(T, vdot(sel, x, y))::T
end

# Check compatibility os selected indices with other specifed array(s) and
# return whether the selection is non-empty.
@inline function checkselection(sel::AbstractVector{Int},
                                A::AbstractArray{<:Any,N}) where {N}
    @certify IndexStyle(sel) === IndexLinear()
    @certify IndexStyle(A) === IndexLinear()
    flag = !isempty(sel)
    if flag
        imin, imax = extrema(sel)
        I = eachindex(IndexLinear(), A)
        ((first(I) ≤ imin) & (imax ≤ last(I))) || out_of_range_selection()
    end
    return flag
end

@inline function checkselection(sel::AbstractVector{Int},
                                A::AbstractArray{<:Any,N},
                                B::AbstractArray{<:Any,N}) where {N}
    @certify IndexStyle(B) === IndexLinear()
    axes(A) == axes(B) || arguments_have_incompatible_axes()
    checkselection(sel, A)
end

@inline function checkselection(sel::AbstractVector{Int},
                                A::AbstractArray{<:Any,N},
                                B::AbstractArray{<:Any,N},
                                C::AbstractArray{<:Any,N}) where {N}
    @certify IndexStyle(B) === IndexLinear()
    @certify IndexStyle(C) === IndexLinear()
    axes(A) == axes(B) == axes(C) || arguments_have_incompatible_axes()
    checkselection(sel, A)
end

@noinline out_of_range_selection() =
    bad_argument("some selected indices are out of range")


text:
struct BackwardLocal <: MainPathAlgorithm end

struct BackwardLocalResult{T<:Integer, U<:Real} <: MainPathResult
  edges::Vector{Edge{T}}
  weights::Vector{U}
end

function backward_local(
  g::AbstractGraph{T},
  s::Vector{R},
  weights::AbstractMatrix{U}
  ) where {T,R <: Integer} where U <: Real

  source = ifelse.(indegree(g) .== 0, true, false)
  visited = falses(nv(g))
  S = copy(s)
  mp_e = Vector{Edge{T}}()

  while !all(source[S])
    next = Vector{T}()
    for u in S
      (source[u] || visited[u]) && continue
      maxnb = max_inneighbors(g, u, weights)
      append!(next, maxnb)
      [push!(mp_e, Edge(nb, u)) for nb in maxnb]
      visited[u] = true
    end
    S = next
  end

  from = src.(mp_e)
  to = dst.(mp_e)
  w_e = [weights[s,t] for (s, t) in zip(from, to)]
  return BackwardLocalResult(mp_e, w_e)
end

function main_path(
  g::AbstractGraph{T},
  s::AbstractVector{S},
  weights::AbstractMatrix{U},
  ::BackwardLocal
  ) where {S, T <: Integer} where {U <: Real}

  return backward_local(g, s, weights)
end

function main_path(
  g::AbstractGraph{T},
  s::AbstractVector{S},
  ::BackwardLocal
  ) where {S, T <: Integer}

  return backward_local(g, s, weights(g))
end

function main_path(
  g::AbstractGraph{T},
  s::Integer,
  ::BackwardLocal
  ) where {T <: Integer}

   return backward_local(g, [s], weights(g))
 end

function main_path(
  g::AbstractGraph{T},
  ::BackwardLocal
  ) where {T <: Integer}

  w = weights(g)
  sinks = get_sinks(g)
  sweights = [max_inweight(g, v, w) for v in sinks]
  s = sinks[sweights .== maximum(sweights)]

  return backward_local(g, s, w)
end


text:
using Base.Meta

export emit

function emit(arg::ParsedArgument)
  local out::Expr = Expr(:(::), arg.name, arg.typ)
  if arg.varargs
    out = Expr(:(...), out)
  end
  if isdefined(arg, :default)
    out = Expr(:kw, out, arg.default)
  end
  return out
end

emit(args::Vector{ParsedArgument}) = map(emit, args)

function emit_args(func::ParsedFunction)
  if isdefined(func, :kwargs)
    return [Expr(:parameters, emit(func.kwargs)...), emit(func.args)...]
  else
    return emit(func.args)
  end
end

make_quotenode(s::Symbol) = eval(quot(quot(s)))

function emit_name(namespace::Vector{Symbol}, name::Symbol)
  if length(namespace)>0
    return Expr(:(.), emit_name(namespace[1:end-1], namespace[end]), make_quotenode(name))
  else
    return name
  end
end

function emit_name(func::ParsedFunction)
  if !isdefined(func, :name)
    return []
  elseif isdefined(func, :types)
    if isdefined(func, :namespace)
      return [Expr(:curly, emit_name(func.namespace, func.name), func.types...)]
    else
      return [Expr(:curly, func.name, func.types...)]
    end
  elseif isdefined(func, :namespace)
    return [emit_name(func.namespace, func.name)]
  else
    return [func.name]
  end
end

function emit(func::ParsedFunction)
  const head = isdefined(func, :name) ? :call : :tuple # :call for a method else :tuple
  return Expr(:function, Expr(head, emit_name(func)..., emit_args(func)...), func.body)
end

macro commutative(ex::Expr)
  # Parse the function
  pfunc = ParsedFunction(ex)

  if (!isdefined(pfunc, :name) ||
      length(pfunc.args) != 2 ||
      pfunc.args[1].varargs || pfunc.args[2].varargs ||
      isdefined(pfunc.args[1], :default) || isdefined(pfunc.args[2], :default) ||
      (isdefined(pfunc, :kwargs) && length(pfunc.kwargs) != 0))
    error("Expected a method of 2 simple arguments only")
  end

  if pfunc.args[1].typ == pfunc.args[2].typ
    error("Arguments must be of different types")
  end

  # commutative so we can just swap the arguments
  reverse!(pfunc.args)

  if pfunc.args[1].typ == :Any || pfunc.args[2].typ == :Any
    # If one of the arguments is nothing we'll get warning unless we
    # define a method with symetric types first
    pfunc2 = deepcopy(pfunc)
    if pfunc2.args[1].typ == :Any
      pfunc2.args[1].typ = pfunc2.args[2].typ
    else
      pfunc2.args[2].typ = pfunc2.args[1].typ
    end
    # produce three forms of the method
    return esc(quote
      $(emit(pfunc2))
      $ex
      $(emit(pfunc))
    end)
  else
    # produce two forms of the method
    return esc(quote
      $ex
      $(emit(pfunc))
    end)
  end
end


text:
module TestProgress

export LoudTestSet

using Test
using Test: AbstractTestSet, DefaultTestSet, Pass, finish, record

const DOT = '•'

struct LoudTestSet{T<:AbstractTestSet} <: AbstractTestSet
    wrapped::T

    LoudTestSet{T}(desc) where T = new(T(desc))
end

LoudTestSet(desc::AbstractString) = LoudTestSet{DefaultTestSet}(desc)

Test.record(ts::LoudTestSet, res) = record(ts.wrapped, res)
function Test.record(ts::LoudTestSet, res::Pass)
    printstyled(DOT, color=:green)
    record(ts.wrapped, res)
    return res
end

function Test.finish(ts::LoudTestSet)
    # TODO: only print nested testsets up to some maximum depth?
    printstyled("\nFinished " * ts.wrapped.description, color=:cyan)
    print("\n")
    finish(ts.wrapped)
end

end  # module


text:
"""
	readcampbell(filein)
Load data stored in Campbell logger file format to DataFrame

**Input**
* filein: full name of the input file

**Output**
* dataframe containing columns + timevector in DateTime format
* data units

**Example**
```
data,units = readcampbell("../test/input/campbell_data.tsf");
```
"""
function readcampbell(filein::String)
	channels = Array{SubString{String},1}();
	units = Array{SubString{String},1}();
	# Read header
	open(filein,"r") do fid
		readline(fid);# firts line is not important
		# Read channel names
		temp = split(replace(readline(fid),r" |\r|\""=>""),','); # remove funny symbols
		channels = temp[2:end];
		# Read units
		temp = split(replace(readline(fid),r" |\r|\""=>""),','); # remove funny symbols
		units = temp[2:end];
	end
	# read data
	dataall = readdlm(filein,',',skipstart=4,header=false);
	data = DataFrame(datetime=DateTime.(dataall[:,1],"yyyy-mm-dd HH:MM:SS"))
	for (i,v) in enumerate(channels)
		data[Symbol(v)] = dataall[:,i+1];
	end
	return data,units;
end


text:
using LeastSquaresOptim

function rosenbrock(x)
    [1 - x[1], 100 * (x[2]-x[1]^2)]
end
x0 = zeros(2)
optimize(rosenbrock, x0, Dogleg(), lower = fill(0.0, length(x0)))
optimize(rosenbrock, x0, LevenbergMarquardt(), lower = fill(0.0, length(x0)))

max_stars_repo_name: matthieugomez/LeastSquares

max_issues_repo_name: matthieugomez/LeastSquares

max_forks_repo_name: matthieugomez/LeastSquares


text:
using Distributions, Random
using Statistics
using SpecialFunctions
const __TRUNC = 0.64;
const __TRUNC_RECIP = 1.0 / __TRUNC;
"""
    PolyaGamma(b::Int, c::Real)

## Arguments
- `b::Int`
- `c::Real` exponential tilting

## Keyword Arguments
- `nmax::Int` Sampling parameter
- `trunc::Int` Sampling parameter

Create a PolyaGamma sampler with parameters `b` and `c`
"""
struct PolyaGamma{Tc,A} <: Distributions.ContinuousUnivariateDistribution
    # For sum of Gammas.
    b::Int
    c::Tc
    trunc::Int
    nmax::Int
    bvec::A
    #Constructor
    function PolyaGamma{T}(b::Int, c::T, trunc::Int, nmax::Int) where {T<:Real}
        if trunc < 1
            @warn "trunc < 1. Setting trunc=1."
            trunc = 1
        end
        bvec = [convert(T, (twoπ * (k - 0.5))^2) for k in 1:trunc]
        return new{typeof(c),typeof(bvec)}(b, c, trunc, nmax, bvec)
    end
end

Statistics.mean(d::PolyaGamma) = d.b / (2 * d.c) * tanh(0.5 * d.c)

function PolyaGamma(b::Int, c::T; nmax::Int=10, trunc::Int=200) where {T<:Real}
    return PolyaGamma{T}(b, c, trunc, nmax)
end

function Distributions.pdf(d::PolyaGamma, x)
    return cosh(d.c / 2)^d.b * 2.0^(d.b - 1) / gamma(d.b) * sum(
        ((-1)^n) * gamma(n + d.b) / gamma(n + 1) * (2 * n + b) / (sqrt(2 * π * x^3)) *
        exp(-(2 * n + b)^2 / (8 * x) - c^2 / 2 * x) for n in 0:(d.nmax)
    )
end

## Sampling
function Distributions.rand(rng::AbstractRNG, d::PolyaGamma{T}) where {T<:Real}
    if iszero(d.b)
        return zero(T)
    end
    return sum(Base.Fix1(draw_like_devroye, rng), d.c * ones(d.b))
end

## Utility functions
function a(n::Int, x::Real)
    k = (n + 0.5) * π
    if x > __TRUNC
        return k * exp(-0.5 * k^2 * x)
    elseif x > 0
        expnt = -1.5 * (log(0.5 * π) + log(x)) + log(k) - 2.0 * (n + 0.5)^2 / x
        return exp(expnt)
    end
end

function mass_texpon(z::Real)
    t = __TRUNC

    fz = 0.125 * π^2 + 0.5 * z^2
    b = sqrt(1.0 / t) * (t * z - 1)
    a = sqrt(1.0 / t) * (t * z + 1) * -1.0

    x0 = log(fz) + fz * t
    xb = x0 - z + logcdf(Distributions.Normal(), b)
    xa = x0 + z + logcdf(Distributions.Normal(), a)

    qdivp = 4 / π * (exp(xb) + exp(xa))

    return 1.0 / (1.0 + qdivp)
end

function rtigauss(rng::AbstractRNG, z::Real)
    z = abs(z)
    t = __TRUNC
    x = t + 1.0
    if __TRUNC_RECIP > z
        alpha = 0.0
        rate = 1.0
        d_exp = Exponential(1.0 / rate)
        while (rand(rng) > alpha)
            e1 = rand(rng, d_exp)
            e2 = rand(rng, d_exp)
            while e1^2 > 2 * e2 / t
                e1 = rand(rng, d_exp)
                e2 = rand(rng, d_exp)
            end
            x = 1 + e1 * t
            x = t / x^2
            alpha = exp(-0.5 * z^2 * x)
        end
    else
        mu = 1.0 / z
        while (x > t)
            y = randn(rng)^2
            half_mu = 0.5 * mu
            mu_Y = mu * y
            x = mu + half_mu * mu_Y - half_mu * sqrt(4 * mu_Y + mu_Y^2)
            if rand(rng) > mu / (mu + x)
                x = mu^2 / x
            end
        end
    end
    return x
end

# ////////////////////////////////////////////////////////////////////////////////
# 				  // Sample //
# ////////////////////////////////////////////////////////////////////////////////

function draw_like_devroye(rng::AbstractRNG, c::Real)
    # Change the parameter.
    c = 0.5 * abs(c)

    # Now sample 0.25 * J^*(1, Z := Z/2).
    fz = 0.125 * π^2 + 0.5 * c^2
    # ... Problems with large Z?  Try using q_over_p.
    # double p  = 0.5 * __PI * exp(-1.0 * fz * __TRUNC) / fz;
    # double q  = 2 * exp(-1.0 * Z) * pigauss(__TRUNC, Z);

    x = 0.0
    s = 1.0
    y = 0.0
    # int iter = 0; If you want to keep track of iterations.
    d_exp = Exponential()
    while true
        if rand(rng) < mass_texpon(c)
            x = __TRUNC + rand(rng, d_exp) / fz
        else
            x = rtigauss(rng, c)
        end
        s = a(0, x)
        y = rand(rng) * s
        n = 0
        go = true

        # Cap the number of iterations?
        while (go)
            n = n + 1
            if isodd(n)
                s = s - a(n, x)
                if y <= s
                    return 0.25 * x
                end
            else
                s = s + a(n, x)
                if y > s
                    go = false
                end
            end
        end
        # Need Y <= S in event that Y = S, e.g. when x = 0.
    end
end # draw_like_devroye


text: include("./arithmetic.jl")