text:
struct MTL_short
	band::Int
	rad_mul::Float64
	rad_add::Float64
	rad_max::Float64
	reflect_mul::Float64
	reflect_add::Float64
	reflect_max::Float64
	sun_azim::Float64
	sun_elev::Float64
	sun_dist::Float64
	K1::Float64
	K2::Float64
end

const KW = Dict{Symbol,Any}
find_in_dict = GMT.find_in_dict

const Lsat8_desc = Dict(
	1  => "Band 1 - Coastal aerosol [0.43-0.45]",
	2  => "Band 2 - Blue [0.45-0.51]",
	3  => "Band 3 - Green [0.53-0.59]",
	4  => "Band 4 - Red [0.64-0.67]",
	5  => "Band 5 - NIR [0.85-0.88]",
	6  => "Band 6 - SWIR 1 [1.57-1.65]",
	7  => "Band 7 - SWIR 2 [2.11-2.29]",
	9  => "Band 9 - Cirrus [1.36-1.38]",
	10 => "Band 10 - Thermal IR 1 [10.6-11.19]",
	11 => "Band 11 - Thermal IR 2 [11.50-12.51]")

const Sentinel2_10m_desc = Dict(
	2  => "Band 2 - Blue [0.490]",
	3  => "Band 3 - Green [0.560]",
	4  => "Band 4 - Red [0.665]",
	8  => "Band 8 - NIR [0.842]")

const Sentinel2_20m_desc = Dict(			# Common to both 20 & 60 m
	1  => "Band 1 - Coastal aerosol [0.443]",	# Actually only available in 60 m.
	2  => "Band 2 - Blue [0.490]",
	3  => "Band 3 - Green [0.560]",
	4  => "Band 4 - Red [0.665]",
	5  => "Band 5 - Red Edge 1 [0.705]",
	6  => "Band 6 - Red Edge 2 [0.740]",
	7  => "Band 7 - Red Edge 3 [0.783]",
	8  => "Band 8A - Red Edge 4 (NIR) [0.865]",	# ~Landsat8 Band 5
	9  => "Band 9 - Water vapour [0.945]",
	10 => "Band 10 - Cirrus [1.375]",		# ~Landsat8 Band 9
	11 => "Band 11 - SWIR 1 [1.610]",		# ~Landsat8 Band 6
	12 => "Band 12 - SWIR 2 [2.190]")		# ~Landsat8 Band 7

# ----------------------------------------------------------------------------------------------------------
"""
    Irgb = truecolor(bndR, bndG, bndB)

Take three Landsat8/Sentinel2 UINT16 GMTimages or the file names of those bands and compose
an RGB true color image applying automatic histogram stretching.

Return an UInt8 RGB GMTimage

    Irgb = truecolor(cube::GMTImage, bands::Vector{Int})

Make an RGB composition of the 3 bands passed in the vector 'bands' from the layers in the multi-layered GMTimage `cube`

Return an auto-stretched UInt8 RGB GMTimage

    Irgb = truecolor(cube::String, [bands::Vector{Int}], [bandnames::Vector{String}], [raw=false])

Make an RGB composition of 3 bands from the `cube` file holding a UInt16 multi-layered array (often created with `cutcube`)
The band selection can be made with `bands` vector, case in which we will search for bands named "Band band[k]"
or where the bands description contain the contents of `bandnames`. If none of `bands` or `bandnames` is used
we search for a made up `bandnames=["red", "green", "blue"]`.

Return an auto-stretched UInt8 RGB GMTimage OR a GMTimage{UInt16,3} if the `raw` option is set to `true`.

### Example:
Make an RGB composite from data in the cube file "LC08__cube.tiff"
```julia
I = truecolor("LC08__cube.tiff");
```
"""
function truecolor(bndR, bndG, bndB)
	I = isa(bndR, GMT.GMTimage) ? bndR : gmtread(bndR)
	img = Array{UInt8}(undef, size(I,1), size(I,2), 3)
	if (eltype(I) == UInt8)
		img[:,:,1] .= I.image
	else
		_ = mat2img(I.image, stretch=true, img8=view(img,:,:,1), scale_only=1)
	end
	I = isa(bndG, GMT.GMTimage) ? bndG : gmtread(bndG)
	@assert size(I,1) == size(img,1) && size(I,2) == size(img,2)
	if (eltype(I) == UInt8)
		img[:,:,2] .= I.image
	else
		_ = mat2img(I.image, stretch=true, img8=view(img,:,:,2), scale_only=1)
	end
	I = isa(bndB, GMT.GMTimage) ? bndB : gmtread(bndB)
	@assert size(I,1) == size(img,1) && size(I,2) == size(img,2)
	if (eltype(I) == UInt8)
		img[:,:,3] .= I.image
	else
		_ = mat2img(I.image, stretch=true, img8=view(img,:,:,3), scale_only=1)
	end
	Io = mat2img(img, I);	Io.layout = "TRBa"
	Io
end

truecolor(cube::GMT.GMTimage{UInt16, 3}, layers::Vector{Int}) = truecolor(cube, layers=layers)
function truecolor(cube::GMT.GMTimage{UInt16, 3}; layers::Vector{Int}=Int[])
	(length(layers) != 3) && error("For an RGB composition 'bands' must be a 3 elements array and not $(length(layers))")
	(cube.layout[3] != 'B') && error("For an RGB composition the image object must be Band interleaved and not $(cube.layout)")
	img = Array{UInt8, 3}(undef, size(cube,1), size(cube,2), 3)
	layers = find_layers(cube, layers, 3)
	_ = mat2img(@view(cube[:,:,layers[1]]), stretch=true, img8=view(img,:,:,1), scale_only=1)
	_ = mat2img(@view(cube[:,:,layers[2]]), stretch=true, img8=view(img,:,:,2), scale_only=1)
	_ = mat2img(@view(cube[:,:,layers[3]]), stretch=true, img8=view(img,:,:,3), scale_only=1)
	Io = mat2img(img, cube);	Io.layout = "TRBa"
	Io
end

truecolor(cube::GMT.GMTgrid{Float32, 3}, layers::Vector{Int}) = truecolor(cube, layers=layers)
function truecolor(cube::GMT.GMTgrid{Float32, 3}; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[],
	               bandnames::Vector{String}=String[])
	(isempty(bands) && isempty(bandnames) && isempty(layers)) && (bandnames = ["red", "green", "blue"])
	isempty(layers) && (layers = find_layers(cube, bandnames, bands))
	(length(layers) != 3) && error("For an RGB composition 'bands' must be a 3 elements array and not $(length(layers))")
	img = Array{UInt8, 3}(undef, size(cube,1), size(cube,2), 3)
	img[:,:,1] = rescale(@view(cube[:,:,layers[1]]), stretch=true, type=UInt8)
	img[:,:,2] = rescale(@view(cube[:,:,layers[2]]), stretch=true, type=UInt8)
	img[:,:,3] = rescale(@view(cube[:,:,layers[3]]), stretch=true, type=UInt8)
	Io = mat2img(img, cube);	Io.layout = "TRBa"
	Io
end

function truecolor(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], raw::Bool=false)
	# The `raw` option returns a GMTimage{UInt16, 3} and does not convert to UInt8 with auto-stretch (default)
	(isempty(bands) && isempty(bandnames) && isempty(layers)) && (bandnames = ["red", "green", "blue"])
	rgb = subcube(cube, bands=bands, layers=layers, bandnames=bandnames)
	(raw) && return rgb
	return (eltype(rgb) <: AbstractFloat) ? truecolor(rgb, layers=[1,2,3]) : mat2img(rgb, stretch=:auto)
end

# ----------------------------------------------------------------------------------------------------------
"""
subcube(cube::String; bands=Int[], bandnames=String[], layers=Int[])

Extracts a subcube from `cube` with the layers in the `bands` vector, case in which we will search for bands
named "Band band[k]", or those whose names correspond (even partially and case insensitive) to the descriptions
in `bandnames` string vector. This means that the options `bands` and `bandnames` can only be used in 'cubes'
with bands description. The `layers` option blindly extract the `cube` planes listed in the `layer` vector.

Returns a GMTimage

### Example
Extracts the Red, Green and Blue layers from a Landsat 8 cube created with `cutcube`

```
Irgb = subcube("LC08__cube.tiff", bandnames = ["red", "green", "blue"])
```
"""
function subcube(cube::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[])
	# This is also used as a helper function in 'truecolor' and the radiometric indices functions.
	(isempty(bands) && isempty(bandnames) && isempty(layers)) && error("Must specify at least one way to select layers.")
	(isempty(layers)) && (layers = find_layers(cube, bands=bands, bandnames=bandnames)[1])
	gmtread(cube, band=layers, layout="TRBa")		# This one is still UInt16
end

# ----------------------------------------------------------------------------------------------------------
function find_layers(cube::GMT.GMTimage{UInt16, 3}, list::Vector{Int}, n_layers::Int)
	# This function is not finished
	if (maximum(list) < 200)		# The list of bands to pass to the caling fun
		(maximum(list) > size(cube,3)) && error("Not enough bands to satisfy the bands list request.")
		bands = list
	#=
	elseif (!isempty(cube.v))		# Must match the wavelength in `list` with those of cube.v
		bands, n = zeros(Int, length(list)), 0
		for w in list
			d, ind = findmin(abs.(w .- cube.v))
			if (d < tol)
				bands[n += 1] = ind
				continue
			end
		end
		(n != length(list)) && error("Some of the wavelegth in $(list) are not present in the `cube.v`")
		bands
	=#
	else
		error("The `cube` object does not have a frequencies (`v` coordinates) vector as required here.")
	end
	(length(bands) != n_layers) && error("Need $(n_layers) bands but got $(length(bands))")
	bands
end

#function find_layers(cube::GMT.GMTimage{UInt16, 3}, bandnames::Vector{String}=String[], bands::Vector{Int}=Int[])
function find_layers(cube::AbstractArray, bandnames::Vector{String}=String[], bands::Vector{Int}=Int[])
	# Find the layers corresponding to the (parts of) contents of "bandnames"
	(!isa(cube, GMT.GMTimage{UInt16, 3}) && !isa(cube, GMT.GMTgrid{Float32, 3})) && error("'cube' must be a 3D GMTimage or GMTgrid")
	(isempty(cube.names)) && error("The `cube` object does not have a `names` (band names) assigned field as required here.")
	(!isempty(bands)) && (bandnames = ["Band $(bands[k])" for k = 1:length(bands)])		# Create a bandnames vector
	helper_find_layers(lowercase.(cube.names), bandnames)
end

function find_layers(fname::String; bands::Vector{Int}=Int[], layers::Vector{Int}=Int[], bandnames::Vector{String}=String[], alllayers::Bool=false)
	# 'layers' just return itself after checking that the cube file actually contain that many layers
	# 'bands' search for "Band bands[k]"
	# 'bandnames' search the bands description for the first layer that contains bandnames[k]
	# Returns the numeric layers (1-based) corresponding to the search criteria and the bands description
	(!alllayers && isempty(bands) && isempty(layers) && isempty(bandnames)) &&
		error("Must use either the 'bands' OR the 'bandnames' option.")
	GMT.ressurectGDAL()		# Yep, shit, sometimes its needed.
	ds = GMT.Gdal.unsafe_read(fname)
	nbands = GMT.Gdal.nraster(ds)
	msg = ""
	(nbands < 2) && (msg = "This file ($fname) does not contain cube data (more than one layer).")
	(!isempty(layers) && maximum(layers) > nbands) && (msg = "Asked for more 'layers' than this cube contains.")
	(msg != "") && (GMT.Gdal.GDALClose(ds.ptr); error(msg))
	desc = fill("", nbands)
	[desc[k] = GMT.Gdal.GDALGetDescription(GMT.Gdal.GDALGetRasterBand(ds.ptr, k)) for k = 1:nbands]
	GMT.Gdal.GDALClose(ds.ptr)

	(!isempty(layers)) && return layers, desc
	(all(desc .== "")) && error("This cube file has no band descriptions so cannot use the 'band' or 'bandnames' options.")
	(alllayers) && return collect(1:nbands), desc	# OK, just return them ALL

	(!isempty(bands)) && (bandnames = ["Band $(bands[k])" for k = 1:length(bands)])		# Create a bandnames vector
	_layers = helper_find_layers(lowercase.(desc), bandnames)
	return _layers, desc
end

helper_find_layers(desc::Vector{String}, band_names::String) = helper_find_layers(desc, [band_names])
function helper_find_layers(desc::Vector{String}, band_names::Vector{String})::Vector{Int}
	# Note, 'desc' is supposed to be in lowercase already.
	bnd_names = lowercase.(band_names)
	layers, n = zeros(Int, length(bnd_names)), 0
	for bnd_name in bnd_names
		((ind = findfirst(contains.(desc, bnd_name))) !== nothing) && (layers[n += 1] = ind)
	end
	(n != length(bnd_names)) && error("All or some of the names in $(band_names) are not present in the `cube` names/description.\n\tUse the `reportbands()` function to check the bands description.")
	layers
end

# ----------------------------------------------------------------------------------------------------------
"""
    reportbands(in; [layers=Int[]])
or

    reportbands(in, layer;)

Report the Bands description of the `in` input argument. This can be a GMTimage, a GMTgrid or a file name (a String) of
a 'cube' file. Normally one made with the `cutcube` function. When the use conditions of this function are not met,
either a warning or an error message (if too deep to be caught as a warning) will be issued.

- `layers`: When this optional parameter is used, report the description of the bands in the vector `layers`
- `layer`: A scalar with a unique band number. Alternative form to `reportbands(in, layers=[layer])`

Returns a string vector.
"""
reportbands(in, layer::Int) = reportbands(in, layers=[layer])
function reportbands(in; layers::Vector{Int}=Int[])
	(!isa(in, GMT.GMTimage) && !isa(in, GMT.GMTgrid) && !isa(in, String)) &&
		error("Bad input. Must be a GMTimage, a GMTgrid or a file name. Not $(typeof(in))")
	if (isa(in, GMT.GMTimage) || isa(in, GMT.GMTgrid))
		isempty(in.names) && (println("This image object does not have a `names` assigned field"); return nothing)
		return (isempty(layers)) ? in.names : in.names[layers]
	else
		layers, desc = find_layers(in, layers=layers, alllayers=isempty(layers))
		return desc[layers]
	end
end

# ----------------------------------------------------------------------------------------------------------
"""
    cutcube(names=String[], bands=Int[], template="", region=nothing, extension=".TIF", description=String[], mtl="", sentinel2=0, save="")

Cut a 3D cube out of a Landsat/Sentinel scene within a subregion `region` and a selection of bands.

- `names`: (optional) A vector with the individual bands full file name
- `bands`: When `names` is not provided give a vector of integers corresponding to the choosen bands.
           This works well for Landsat and most of Sentinel bands. However, in later case, there are also
           bands that contain characters, for example band 8A. In this case `bands` should be a vector of
           strings including the extension. _e.g._ ["02.jp2", "8A.jp2"]
- `template`: Goes together with the `bands` option. They are both composed a template * band[n] to recreate
           the full file name of each band.
- `region` Is the region to extract and must contain the extracting region limits as [W, E, S, N] or a
           GMT style -R string (without the leading "-R").
- `extension`: In case the `bands` is numeric but file extensions are not "*.TIF" (case insensitive),
           use the extension passed by this option.
- `description`: A vector of strings (as many as bands) with a description for each band. If not provided and
           the file is recognized as a Landasat 8, band description is added automatically, otherwise
           we build one with the bands file names. This info will saved if data is written to a file.
- `mtl`:   If reading from Landsat and the MTL file is not automatically found (you get an error) use this
           option to pass the full name of the MTL file.
- `sentinel2`: ESA is just unconsistent and names change with time and band numbers can have character (e.g. 8A)
           hence we need help to recognize Sentinel files so the known description can be assigned.
           Use `sentinel=10`, or `=20` or `=60` to indicate Sentinel files at those resolutions.
- `save`:  The file name where to save the output. If not provided, a GMTimage is returned.

Return: `nothing` if the result is written in file or a GMTimage otherwise.

## Examples

```julia
# Cut a Landsat 8 scene for a small region (in UTM) and return a GMTimage with 3 bands in UInt16.
temp = "C:\\SIG_AnaliseDadosSatelite\\SIG_ADS\\DadosEx2\\LC82040332015145LGN00\\LC82040332015145LGN00_B";
cube = cutcube(bands=[2,3,4], template=temp, region=[479670,492720,4282230,4294500])

# The same example as above but save the data in a GeoTIFF disk file and use a string for `region`
cutcube(bands=[2,3,4], template=temp, region="479670/492720/4282230/4294500", save="landsat_cube.tif")
```
"""
function cutcube(; names::Vector{String}=String[], bands::AbstractVector=Int[], template::String="", region=nothing, extension::String=".TIF", description::Vector{String}=String[], save::String="", sentinel2::Int=0, mtl::String="")
	(region === nothing) && error("The `region` option cannot be empty.")
	if (isempty(names))
		(isempty(bands) || template == "") && error("When band file `names` are not provided, MUST indicate `bands` AND `template`")
		!isa(bands, Vector{<:Integer}) && !isa(bands, Vector{<:String}) && error("`bands` must a vector of Int or Strings.")
		names = Vector{String}(undef, length(bands))
		if isa(bands, Vector{Int})
			[names[k] = @sprintf("%s%d%s", template, bands[k], extension) for k = 1:length(bands)]
			!isfile(names[1]) &&	# Landsat uses "B2.TIF" and Sentinel "B02.jp2", try again.
				[names[k] = @sprintf("%s%.02d%s", template, bands[k], extension) for k = 1:length(bands)]
		else
			[names[k] = @sprintf("%s%s", template, bands[k]) for k = 1:length(bands)]
		end
	end
	# Now test if any of the file names, either passed in or generated here, do not exist
	for name in names
		!isfile(name) && error("File name $name does not exist. Must stop here.")
	end
	MTL = read_mtl(names[1], mtl, get_full=true)
	(MTL !== nothing) && (MTL = ["MTL=" * join(MTL, "\n")])

	# Little parsing of the -R string but does not test if W < E & S < N
	_region::String = isa(region, String) ? region :
	                  @sprintf("%.12g/%.12g/%.12g/%.12g", region[1], region[2], region[3], region[4])
	startswith(_region, "-R") && (_region = _region[3:end])		# Tolerate a region that starts with "-R"
	(length(findall("/", _region)) != 3) && error("Badly formed region string: $_region")

	desc = assign_description(names, description, sentinel2)
	cube = grdcut(names[1], R=_region)
	mat = cube.image
	for k = 2:length(bands)
		B = grdcut(names[k], R=_region)
		mat = cat(mat, B.image, dims=3)
	end
	cube = mat2img(mat, cube, names=desc)
	(save != "") && gdaltranslate(cube, dest=save, meta=MTL)
	return (save != "") ? nothing : cube
end

function assign_description(names::Vector{String}, description::Vector{String}, sentinel2::Int=0)
	# Create a description for each band. If 'description', the cutcube() kwarg, is provided we use it as is.
	# Next we try to find if 'names' indicate a Landsat8 origin and if yes we use the known names & frequencies
	# Otherwise we use the file names as descriptors.
	(!isempty(description) && length(names) != length(description)) &&
		error("'description' and 'names' vectors must have the same length")
	desc = (!isempty(description)) ? description : Vector{String}(undef, length(names))

	t = splitext(splitdir(names[1])[2])[1]
	if (isempty(description) && (startswith(t, "LC08_") || startswith(t, "LC8")))	# Have Landsat8 data
		for k = 1:length(names)
			t = splitext(splitdir(names[k])[2])[1]
			ind = findfirst("_B", t)
			bnd = parse(Int, t[ind[end]+1:end])
			desc[k] = Lsat8_desc[bnd]
		end
	elseif (sentinel2 == 10 || sentinel2 == 20 || sentinel2 == 60)
		for k = 1:length(names)
			t = splitdir(names[k])[2]
			ind = findfirst("_B", t)
			bnd = tryparse(Int, t[ind[end]+1:ind[end]+2])
			(bnd === nothing && t[ind[end]+1:ind[end]+2] == "8A") && (bnd = 8)	# ESA is crazzy. Fck character
			desc[k] = (sentinel2 == 10) ? Sentinel2_10m_desc[bnd] : Sentinel2_20m_desc[bnd]	# 20m applyies also for 60m
		end
	else
		if (isempty(description))
			[desc[k] = splitext(splitdir(names[k])[2])[1] for k = 1:length(names)]
		else
			desc = description
		end
	end
	return desc
end

# ----------------------------------------------------------------------------------------------------------
"""
read_mtl(band_name::String, mtl::String=""; get_full=false)

Use the `band_name` of a Landsat8 band to find the MTL file with the scene parameters at which that band
belongs and read the params needed to compute Brightness temperature, radiance at top of atmosphere, etc.
If the MTL file does not lieve next to the band file, send its name via the `mtl` argument.

The `get_full` option makes this function return a tring with contents of the MTL file or `nothing` if
the MTL file is not found.

### Returns a tuple with:

(band=band, rad_mul=rad_mul, rad_add=rad_add, rad_max=rad_max, reflect_mul=reflect_mul, reflect_add=reflect_add, reflect_max=reflect_max, sun_azim=sun_azim, sun_elev=sun_elev, sun_dis=sun_azim, K1=K1, K2=K2)

or a string with MTL contents (or nothing if MTL file is not found)
"""
function read_mtl(fname::String, mtl::String=""; get_full::Bool=false)
	pato,_fname = splitdir(fname)
	((ind = findfirst("_B", _fname)) === nothing && !get_full) &&
		error("This $(fname) is not a valid Landsat8 band file name or of a Landsat 8 cube file.")
	(ind === nothing) && return nothing		# Only happens when get_full = true and name is no Landsat
	if (mtl == "")
		mtl = join(pato, _fname[1:ind[1]] * "MTL.txt")
		if (!isfile(mtl))
			lst = filter(x -> endswith(x, "MTL.txt"), readdir((pato == "") ? "." : pato))
			if (length(lst) == 1 && startswith(lst[1], _fname[1:16]))
				mtl = joinpath(pato, lst[1])
				(!isfile(mtl) && get_full) && return nothing		# Not fatal error in this case
				!isfile(mtl) && (@warn("MTL file was not transmitted in input and I couldn't find it."); return nothing)
			end
		end
	end

	(!get_full) && (__fname = splitext(_fname)[1];	band = parse(Int, __fname[ind[1]+2:end]))
	f = open(mtl);	lines = readlines(f);	close(f)
	return (get_full) ? lines : parse_mtl(lines, band)
end

function parse_mtl(mtl::Vector{<:AbstractString}, band::Int)
	# Parse the 'mtl' string vector and extract the info relevant for 'band'
	# Make this a separate function so it can be called from read_mtl() or from the
	# MTL metadata stored in a cube file created with cutcube().

	function get_par(str)
		ind = findfirst(findfirst.(str, mtl) .!== nothing)
		parse(Float64, split(mtl[ind], "=")[2])
	end

	rad_mul = get_par("RADIANCE_MULT_BAND_$(band)")
	rad_add = get_par("RADIANCE_ADD_BAND_$(band)")
	rad_max = get_par("RADIANCE_MAXIMUM_BAND_$(band)")
	reflect_mul = (band < 10) ? get_par("REFLECTANCE_MULT_BAND_$(band)") : 1.0
	reflect_add = (band < 10) ? get_par("REFLECTANCE_ADD_BAND_$(band)") : 0.0
	reflect_max = (band < 10) ? get_par("REFLECTANCE_MAXIMUM_BAND_$(band)") : 0.0
	sun_azim = get_par("SUN_AZ")
	sun_elev = get_par("SUN_EL")
	sun_dist = get_par("EARTH_SUN")
	K1 = (band >= 10) ? get_par("K1_CONSTANT_BAND_$(band)") : 0.0
	K2 = (band >= 10) ? get_par("K2_CONSTANT_BAND_$(band)") : 0.0
	MTL_short(band, rad_mul, rad_add, rad_max, reflect_mul, reflect_add, reflect_max, sun_azim, sun_elev, sun_dist, K1, K2)
end

# ----------------------------------------------------------------------------------------------------------
function helper1_sats(fname::String, band_layer::Int)
	I::GMT.GMTimage{UInt16, 2} = (band_layer == 0) ? gmtread(fname) : gmtread(fname, layer=band_layer, layout="TRB")
	indNaN = isnodata(I)
	o = Matrix{Float32}(undef, size(I))
	return I, indNaN, o
end

function parse_lsat8_file(fname::String; band::Int=0, layer::Int=0, mtl::String="")
	# See if 'fname' is of a plain Landsat8 .tif file or of a cube created with cutcube().
	# Depending on the case find the MTL info from file or from the cube's Metadata. Former case
	# still accepts that the MTL file name be transmitted via the 'mtl' option.
	# If no errors so far, parse the MTL and extract the parameters concerning the wished 'band'.
	# This band number will be fetch from the band file name (full Landsat8 product name), or must be
	# transmitted via 'band' option when reading a cube file.
	# In alternative to 'band' we can use 'layer' and fetch the Band number that is stored in that layer.
	GMT.ressurectGDAL()				# Another black-hole plug attempt.
	ds = GMT.Gdal.unsafe_read(fname)
	nbands = GMT.Gdal.nraster(ds)
	meta = GMT.Gdal.GDALGetMetadata(ds.ptr, C_NULL)
	if (nbands > 1)
		desc = Vector{String}(undef,0)
		[append!(desc, [GMT.Gdal.GDALGetDescription(GMT.Gdal.GDALGetRasterBand(ds.ptr, bd))]) for bd = 1:nbands]
	else
		desc = [""]
	end
	GMT.Gdal.GDALClose(ds.ptr)

	if (nbands > 1)
		if (0 < layer < 12)			# Try to find which Band is at layer 'layer'
			((_band = tryparse(Int, split(desc[layer])[2])) === nothing) && @warn("Failed to find Band with description in layer $layout")
			(_band !== nothing) && (band = _band)
		end
		(band < 1) && error("The `band` option must contain the wished Landsat 8 band number. Use `band=N` to set it.")
		(band > 11) && throw(ArgumentError("Bad Landsat 8 band number $band"))
		MTL = ((ind = findfirst(startswith.(meta, "MTL=GROUP ="))) !== nothing) ? meta[ind][5:end] :
			error("Data in a cube (> one layer) must contain the MTL info in Metadata and this one does not.")

		((band_layer = findfirst(startswith.(desc, "Band $band"))) === nothing) && error("Band $band not found in this cube")
		pars = parse_mtl(split(MTL, "\n"), band)
	else
		pars = read_mtl(fname, mtl)		# No 'band' here because it's suposed to be findable in 'fname'
		band_layer = 0
	end
	band_layer, pars
end

# ----------------------------------------------------------------------------------------------------------
function helper_dns(fname::String, band::Int, bandname::String, bandnames::String, mtl::String)
	# Helper function for the dn2temperature, dn2radiance, dn2... functions
	(band == 0 && bandname == "" && bandnames == "") && error("No information provided on what Band to process.")
	(bandname == "" && bandnames != "") && (bandname = bandnames)	# Accept singular and plural name
	if (bandname != "")
		layers, = find_layers(fname, bandnames=[bandname])
		band_layer, pars = parse_lsat8_file(fname, layer=layers[1], mtl=mtl)
		(band_layer != layers[1]) && error("Internal error. $band_layer and $(layers[1]) should be equal.")
	else
		band_layer, pars = parse_lsat8_file(fname, band=band, mtl=mtl)
	end
	band_layer, pars
end

const dns_doc = "
- `fname`: The name of either a ``LANDSAT_PRODUCT_ID`` geotiff band, or the name of a cube file created with
  the `cutcube` function. In the first case, if the companion ``...MTL.txt`` file is not in the same directory
  as `fname` one can still pass it via the `mtl=path-to-MTL-file` option. In the second case it is mandatory
  to use one of the following two options.
- `band`: _cubes_ created with [`cutcube`](@ref) assign descriptions starting with \"Band 1 ...\" an so on
  the other bands. So when `band` is used we search for the band named \"Band N\", where N = `band`.
- `bandname`: 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 `bandname` string that will be
  matched against the cube's bands descriptions. We can use the `reportbands` function to see the bands description.
- `save`:  The file name where to save the output. If not provided, a GMTgrid is returned.

Returns a Float32 GMTgrid
"

# ----------------------------------------------------------------------------------------------------------
"""
    R = dn2temperature(fname::String; band::Int=0, mtl::String="", save::String)

Computes the brigthness temperature of Landasat8 termal band (10 or 11)

$dns_doc

# Example:
Compute the brightness temperature of Band 10 stored in a `cube`

```
T = dn2temperature(cube, band=10)
```
"""
function dn2temperature(fname::String; band::Int=0, mtl::String="", save::String="")
	band_layer, pars = parse_lsat8_file(fname, band=band, mtl=mtl)
	(pars.band < 10) && throw(ArgumentError("Brightness temperature is only for bands 10 or 11. Not this one: $(pars.band)"))
	I, indNaN, o = helper1_sats(fname, band_layer)
	@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)
		o[k] = indNaN[k] ? NaN32 : pars.K2 / (log(pars.K1 / (I.image[k] * pars.rad_mul + pars.rad_add) + 1.0)) - 273.15
	end
	G = mat2grid(o, I)
	(save != "") && gdaltranslate(G, dest=save)
	return (save != "") ? nothing : G
end

# ----------------------------------------------------------------------------------------------------------
function helper_dns_op(I::GMTimage, mul::Float64, add::Float64, indNaN::Matrix{Bool}, o::Matrix{Float32})
	# Helper function for the dn2radiance, dn2reflectance functions. Avoid IF branch if no NaNs
	if any(indNaN)
		@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)
			o[k] = indNaN[k] ? NaN32 : I.image[k] * mul + add
		end
	else
		@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)  o[k] = I.image[k] * mul + add  end
	end
	mat2grid(o, I)
end

# ----------------------------------------------------------------------------------------------------------
"""
    R = dn2radiance(fname::String, [band::Int, bandname::String, mtl::String, save::String])

Computes the radiance at TopOfAtmosphere of a Landsat 8 file

$dns_doc

# Example:
Compute the radiance TOA of Band 2 file.
```
R = dn2temperature("LC08_L1TP_204033_20210525_20210529_02_T1_B2.TIF")
```
"""
function dn2radiance(fname::String; band::Int=0, bandname::String="", bandnames::String="", mtl::String="", save::String="")
	dn2aux(fname, "radiance"; band=band, bandname=bandname, bandnames=bandnames, mtl=mtl, save=save)
end

# ----------------------------------------------------------------------------------------------------------
"""
    R = dn2reflectance(fname::String, [band::Int, bandname::String, mtl::String, save::String])

Computes the TopOfAtmosphere planetary reflectance of a Landsat8 file

$dns_doc

# Example:
Compute the reflectance TOA of Red Band stored in a `cube`
```
R = dn2reflectance(cube, bandname="red")
```
"""
function dn2reflectance(fname::String; band::Int=0, bandname::String="", bandnames::String="", mtl::String="", save::String="")
	dn2aux(fname, "reflect"; band=band, bandname=bandname, bandnames=bandnames, mtl=mtl, save=save)
end

function dn2aux(fname::String, fun::String; band::Int=0, bandname::String="", bandnames::String="", mtl::String="", save::String="")
	# Most of dn2reflectance & dn2radiance codes are similar, so gather it here under a common function.
	function inn_helper_rad(fname, band, bandname, bandnames, mtl)
		band_layer, pars = helper_dns(fname, band, bandname, bandnames, mtl)
		I, indNaN, o = helper1_sats(fname, band_layer)
		helper_dns_op(I, pars.rad_mul, pars.rad_add, indNaN, o)
	end
	function inn_helper_ref(fname, band, bandname, bandnames, mtl)
		band_layer, pars = helper_dns(fname, band, bandname, bandnames, mtl)
		(pars.band >= 10) && error("Computing Reflectance for Thermal bands is not defined.")
		I, indNaN, o = helper1_sats(fname, band_layer)
		s_elev = sin(pars.sun_elev * pi/180)
		fact_x = pars.reflect_mul / s_elev
		fact_a = pars.reflect_add / s_elev
		helper_dns_op(I, fact_x, fact_a, indNaN, o)
	end

	helper_fun = (fun == "radiance") ? inn_helper_rad : inn_helper_ref		# Which helper function to use.
	if (band == 0 && bandname == "" && bandnames == "")
		bdnames = reportbands(fname)
		if (fun == "reflect")
			((ind = findfirst(contains.(bdnames, "Band 10"))) !== nothing) && (deleteat!(bdnames, ind))
			((ind = findfirst(contains.(bdnames, "Band 11"))) !== nothing) && (deleteat!(bdnames, ind))
		end
		G = helper_fun(fname, 0, bdnames[1], "", mtl)
		mat = G.z
		for k = 2:length(bdnames)
			t = helper_fun(fname, 0, bdnames[k], "", mtl)
			mat = cat(mat, t.z, dims=3)
		end
		G = mat2grid(mat, G)
		G.names = bdnames
	else
		G = helper_fun(fname, band, bandname, bandnames, mtl)
	end
	(save != "") && gdaltranslate(G, dest=save)
	return (save != "") ? nothing : G
end

# ----------------------------------------------------------------------------------------------------------
"""
    R = reflectance_surf(fname::String, [band::Int, bandname::String, mtl::String, save::String])

Computes the radiance-at-surface of Landsat8 band using the COST model.

$dns_doc
"""
function reflectance_surf(fname::String; band::Int=0, mtl::String="", save::String="")
	band_layer, pars = parse_lsat8_file(fname, band=band, mtl=mtl)
	(pars.band >= 10) && error("Computing Surface Reflectance for Thermal bands is not defined.")
	I, indNaN, o = helper1_sats(fname, band_layer)

	s_elev = sin(pars.sun_elev * pi/180)
	Esun = (pi * pars.sun_dist ^2) * pars.rad_max / pars.reflect_max
	TAUv = 1.0;		TAUz = s_elev;		Esky = 0.0;		sun_prct=1;
	(pars.band == 6 || pars.band == 7 || pars.band == 9) && (TAUz = 1.0)
	Sun_Radiance = TAUv * (Esun * s_elev * TAUz + Esky) / (pi * pars.sun_dist ^2)

	tmp = I.image
	tmp = filter(tmp -> tmp != 0, tmp)		# This is 3 times faster then: tmp = I.image[I.image .> 0]
	darkDN = quantile!(tmp, 0.01)
	@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)
		o[k] = I.image[k] * pars.rad_mul + pars.rad_add
	end
	radiance_dark = pars.rad_mul * darkDN + pars.rad_add	# 0.01%
	LHaze = radiance_dark - sun_prct * Sun_Radiance / 100
	@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)
		o[k] = (o[k] - LHaze) / Sun_Radiance
	end

	@inbounds Threads.@threads for k = 1:size(I,1)*size(I,2)
		if     (o[k] < 0) o[k] = 0
		elseif (o[k] > 1) o[k] = 1
		end
		indNaN[k] && (o[k] = NaN32)
	end
	G = mat2grid(o, I)
	(save != "") && gdaltranslate(G, dest=save)
	return (save != "") ? nothing : G
end

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
module Generators

export generate

struct GenTask
    task::Task
end

Base.iterate(gen::GenTask) = Base.iterate(gen, nothing)
function Base.iterate(gen::GenTask, _)
    result = yieldto(gen.task)
    istaskdone(gen.task) ? nothing : (result, nothing)
end

Base.IteratorSize(_::GenTask) = Base.SizeUnknown()

"""
    generate(func::Function)

Create a generator from func. Func receives a "yield" function as its argument. You can use
this function to produce values.

Here's an example of a recursive generator:

```jldoctest; output = false
flattened(list) = generate() do yield
    fl(el::Vector) = foreach(fl, el)
    fl(el) = yield(el)
    fl(list)
end

[flattened([1, [2, [3], 4], 5])...] == [1,2,3,4,5]

# output

true
```

"""
function generate(func::Function)
    cur = current_task()
    GenTask(Task(()-> begin
        func(result-> yieldto(cur, result))
        schedule(cur)
    end))
end

end

ext: jl


text:
export iauTtut1
"""
Time scale transformation:  Terrestrial Time, TT, to Universal Time,
UT1.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tt1,tt2    double    TT as a 2-part Julian Date
   dt         double    TT-UT1 in seconds

Returned:
   ut11,ut12  double    UT1 as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Notes:

   1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned ut11,ut12 follow
   suit.

   2. The argument dt is classical Delta T.

Reference:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTtut1(double tt1, double tt2, double dt,
#              double *ut11, double *ut12)

function iauTtut1(tt1::Real, tt2::Real, dt::Real)
   ref_tdb1 = Ref{Float64}(0.0)
   ref_tdb2 = Ref{Float64}(0.0)

   status = ccall((:iauTtut1, libsofa_c), Cint,
                   (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tt1), convert(Float64, tt2), convert(Float64, dt),
                   ref_tdb1, ref_tdb2)

   return status, ref_tdb1[], ref_tdb2[]
end


text:
module ReinforcementLearningEnvironmentAtari

using ArcadeLearningEnvironment, GR
import ArcadeLearningEnvironment: getScreen
import ReinforcementLearningBase: AbstractEnv, DiscreteSpace,
       interact!, getstate, reset!, plotenv, actionspace, sample

struct AtariEnv <: AbstractEnv
    ale::Ptr{Nothing}
    screen::Array{UInt8, 1}
    getscreen::Function
    actions::Array{Int32, 1}
    actionspace::DiscreteSpace
    noopmax::Int64
end

actionspace(env::AtariEnv) = env.actionspace

"""
    AtariEnv(name; colorspace = "Grayscale", frame_skip = 4, noopmax = 20,
                   color_averaging = true, repeat_action_probability = 0.)

Returns an AtariEnv that can be used in an RLSetup of the
[ReinforcementLearning](https://github.com/jbrea/ReinforcementLearning.jl)
package. Check the deps/roms folder of the ArcadeLearningEnvironment package to
see all available `name`s.
"""
function AtariEnv(name;
                  colorspace = "Grayscale",
                  frame_skip = 4, noopmax = 20,
                  color_averaging = true,
                  actionset = :minimal,
                  repeat_action_probability = 0.)
    ale = ALE_new()
    setBool(ale, "color_averaging", color_averaging)
    setInt(ale, "frame_skip", Int32(frame_skip))
    setFloat(ale, "repeat_action_probability",
             Float32(repeat_action_probability))
    loadROM(ale, name)
    if colorspace == "Grayscale"
        screen = Array{Cuchar}(undef, 210*160)
        getscreen = getScreenGrayscale
    elseif colorspace == "RGB"
        screen = Array{Cuchar}(undef, 3*210*160)
        getscreen = getScreenRGB
    elseif colorspace == "Raw"
        screen = Array{Cuchar}(undef, 210*160)
        getscreen = getScreen
    end
    actions = actionset == :minimal ? getMinimalActionSet(ale) : getLegalActionSet(ale)
    actionspace = DiscreteSpace(length(actions), 1)
    AtariEnv(ale, screen, getscreen, actions, actionspace, noopmax)
end

function getScreen(p::Ptr, s::Array{Cuchar, 1})
    sraw = getScreen(p)
    for i in 1:length(s)
        s[i] =  sraw[i] .>> 1
    end
end

function interact!(env::AtariEnv, a)
    r = act(env.ale, env.actions[a])
    env.getscreen(env.ale, env.screen)
    (observation=env.screen, reward=r, isdone=game_over(env.ale))
end

function getstate(env::AtariEnv)
    env.getscreen(env.ale, env.screen)
    (observation=env.screen, isdone=game_over(env.ale))
end

function reset!(env::AtariEnv)
    reset_game(env.ale)
    for _ in 1:rand(0:env.noopmax) act(env.ale, Int32(0)) end
    env.getscreen(env.ale, env.screen)
    (observation=env.screen,)
end

imshowgrey(x::Array{UInt8, 2}) = imshowgrey(x[:], size(x))
imshowgrey(x::Array{UInt8, 1}, dims) = imshow(reshape(x, dims), colormap = 2)
imshowcolor(x::Array{UInt8, 3}) = imshowcolor(x[:], size(x))

function imshowcolor(x::Array{UInt8, 1}, dims)
    clearws()
    setviewport(0, dims[1]/dims[2], 0, 1)
    setwindow(0, 1, 0, 1)
    y = (zeros(UInt32, dims...) .+ 0xff) .<< 24
    img = UInt32.(x)
    @simd for i in 1:length(y)
        @inbounds y[i] += img[3*(i-1) + 1] + img[3*(i-1) + 2] << 8 + img[3*i] << 16
    end
    drawimage(0, 1, 0, 1, dims..., y)
    updatews()
end

function plotenv(env::AtariEnv)
    x = zeros(UInt8, 3 * 160 * 210)
    getScreenRGB(env.ale, x)
    imshowcolor(x, (160, 210))
end

export AtariEnv
end # module


text:
as_mpl_type(x::Colorant)                   = toRGBAtuple(x)
as_mpl_type(x::AbstractArray{<:Quantity})  = ustrip(x)
as_mpl_type(x::Quantity)                   = ustrip(x)
as_mpl_type(x::Units)                      = string(x)
as_mpl_type(x)                             = x


text:
module TIWT
export tidwt,
       tiwpt,
       tiwpd,
       itidwt,
       itiwpt,
       itiwpd

using Wavelets

using ..Utils

# ========== Single Step Translation-Invariant Wavelet Transform ==========
"""
    tidwt_step(v, d, h, g)

Performs one level of the translation-invariant discrete wavelet transform (TWDWT) on the
vector `v`, which is the `d`-th level scaling coefficients (Note the 0-th level scaling
coefficients is the raw signal). The vectors `h` and `g` are the detail and scaling filters.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of TIDWT
tidwt_step(v, 0, h, g)
```

**See also:** [`tidwt_step!`](@ref)
"""
function tidwt_step(v::AbstractVector{T}, d::Integer, h::Vector{S}, g::Vector{S}) where
                   {T<:Number, S<:Number}
    n = length(v)
    w₁ = zeros(T, n)
    w₂ = zeros(T, n)

    tidwt_step!(w₁, w₂, v, d, h, g)
    return w₁, w₂
end

"""
    tidwt_step!(w₁, w₂, v, d, h, g)

Same as `tidwt_step` but without array allocation.

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)
w₁ = zeros(8)
w₂ = zeros(8)

# One step of TIDWT
tidwt_step!(w₁, w₂, v, 0, h, g)
```

**See also:** [`tidwt_step`](@ref)
"""
function tidwt_step!(w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     v::AbstractVector{T},
                     d::Integer,
                     h::Vector{S},
                     g::Vector{S}) where {T<:Number, S<:Number}
    # Variable allocations to be used in Wavelets.Transforms.filtdown!
    n = length(v)                           # Signal length
    nc = Utils.nodelength(n, d+1)           # Node length at level d+1
    filtlen = length(h)                     # Filter length
    si = Vector{S}(undef, filtlen-1)        # Temp. filter vector

    # One step of translation-invariant transform
    for b in 0:(1<<d-1)
        # ----- Transform original block -----
        block = @view v[Utils.packet(d,b,n)]
        st = 2*b*nc + 1                     # Starting index in w₁, w₂
        # Decompose using low and high pass filter
        Transforms.filtdown!(g, si, w₁, st, nc, block, 1, 0, false)
        Transforms.filtdown!(h, si, w₂, st, nc, block, 1, -filtlen+1, true)
        # ----- Transform circularly shifted block -----
        block = circshift(block, 1)
        st = (2*b+1)*nc + 1                 # Starting index in w₁, w₂
        # Decompose using low and high pass filter
        Transforms.filtdown!(g, si, w₁, st, nc, block, 1, 0, false)
        Transforms.filtdown!(h, si, w₂, st, nc, block, 1, -filtlen+1, true)
    end
    return w₁, w₂
end

"""
    itidwt_step(w₁, w₂, d, g, h)

Perform one level of the inverse translation-invariant discrete wavelet transform (ITIDWT)
on the vector `w₁` and `w₂`, which are the `d+1`-th level scaling and detail coefficients
respectivaly (Note the 0-th level scaling coefficients correspond to the raw signal). The
vectors `h` and `g` are the high and low pass filters.

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: `d+1`-level scaling coefficients.
- `w₂::AbstractVector{T} where T<:Number`: `d+1`-level detail coefficients.
- `d::Integer`: Depth level of output coefficients.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
`v::Vector{T}`: Output of inverse translation-invariant wavelet transform.

# Examples
```julia
using WaveletsExt

# Setup
w₁, w₂ = randn(8), randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, false)

# One step of iTIDWT
tidwt_step(w₁, w₂, 0, h, g)
```

**See also:** [`itidwt_step!`](@ref)
"""
function itidwt_step(w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     d::Integer,
                     h::Vector{S},
                     g::Vector{S}) where {T<:Number, S<:Number}
    v = similar(w₁)
    itidwt_step!(v, w₁, w₂, d, h, g)
    return v
end

"""
    itidwt_step!(v, w₁, w₂, d, h, g)

Same as `itidwt_step` but without array allocation.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Output vector allocation.
- `w₁::AbstractVector{T} where T<:Number`: `d+1`-level scaling coefficients.
- `w₂::AbstractVector{T} where T<:Number`: `d+1`-level detail coefficients.
- `d::Integer`: Depth level of output coefficients.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
`v::AbstractVector{T}`: Output from inverse translation-invariant wavelet transform of `w₁`
and `w₂`.

# Examples
```julia
using WaveletsExt

# Setup
w₁, w₂ = randn(8), randn(8)
v = similar(w₁)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, false)

# One step of iTIDWT
tidwt_step(v, w₁, w₂, 0, h, g)
```

**See also:** [`itidwt_step`](@ref)
"""
function itidwt_step!(v::AbstractVector{T},
                      w₁::AbstractVector{T},
                      w₂::AbstractVector{T},
                      d::Integer,
                      h::Vector{S},
                      g::Vector{S}) where {T<:Number, S<:Number}
    # Variable allocations to be used in Wavelets.Transforms.filtup!
    n = length(v)                           # Signal length
    np = Utils.nodelength(n, d)             # Node length at level d
    nc = Utils.nodelength(n, d+1)           # Node length at level d+1
    filtlen = length(h)                     # Filter length
    si = Vector{S}(undef, filtlen-1)        # Temp. filter vector
    v₁, v₂ = similar(v), similar(v)         # Temp. output storages

    # One step of translation-invariant transform
    for b in 0:(1<<d-1)
        sp = b*np + 1                       # Starting index in v
        # ----- Inverse transform original block -----
        sc = 2*b*nc + 1                     # Starting index in w₁, w₂
        # Reconstruct using low and high pass filter
        Transforms.filtup!(false, g, si, v₁, sp, np, w₁, sc, -filtlen+1, false)
        Transforms.filtup!(true, h, si, v₁, sp, np, w₂, sc, 0, true)
        # ----- Inverse transform circularly shifted block -----
        sc = (2*b+1)*nc + 1                 # Starting index in w₁, w₂
        # Reconstruct using low and high pass filter
        Transforms.filtup!(false, g, si, v₂, sp, np, w₁, sc, -filtlen+1, false)
        Transforms.filtup!(true, h, si, v₂, sp, np, w₂, sc, 0, true)
        # Adjust shift
        v₂[Utils.packet(d,b,n)] = circshift(v₂[Utils.packet(d,b,n)],-1)
    end
    # Compute result from v₁ and v₂
    for i in eachindex(v)
        @inbounds v[i] = (v₁[i] + v₂[i])/2
    end
    return v
end

# ========== Translation-Invariant DWT ==========
@doc raw"""
    tidwt(x, wt[, L])

Computes the translation-invariant discrete wavelet transform (TIDWT) for `L` levels.

# Arguments
- `x::AbstractVector{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
  \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Matrix{T}`: Output from TIDWT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIDWT
xw = tidwt(x, wt)
```

**See also:** [`itidwt`](@ref)
"""
function tidwt(x::AbstractVector{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)
    n = length(x)
    w = Array{T,2}(undef, (n,L))
    v = copy(x)

    # TIDWT for L levels
    @inbounds begin
        for d in 0:(L-1)
            v, w[:,L-d] = tidwt_step(v, d, h, g)
        end
    end
    return [v w]
end

"""
    itidwt(xw, wt)

Computes the inverse translation-invariant discrete wavelet transform (iTIDWT) on `xw`.

# Arguments
- `xw::AbstractArray{T,2} where T<:Number`: TIDWT-transformed array.
- `wt::OrthoFilter`: Orthogonal wavelet filter.

# Returns
`::Vector{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIDWT
xw = tidwt(x, wt)

# iTIDWT
x̂ = itidwt(xw, wt)
```

**See also:** [`tidwt`](@ref)
"""
function itidwt(xw::AbstractArray{T,2}, wt::OrthoFilter) where T<:Number
    # Setup
    _, k = size(xw)
    g, h = WT.makereverseqmfpair(wt, false)
    x = xw[:,1]             # Output
    L = k-1                 # Number of decomposition levels

    # iTIDWT
    for i in 1:(k-1)
        @inbounds x = itidwt_step(x, xw[:,i+1], L-i, h, g)
    end
    return x
end

# ========== Translation-Invariant WPT ==========
@doc raw"""
    tiwpt(x, wt[, L])

Computes `L` levels of translation-invariant wavelet packet transform (TIWPT) on `x`.

# Arguments
- `x::AbstractVector{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
  \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Matrix{T}`: Output from TIWPT on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIWPT
xw = tiwpt(x, wt)
```

**See also:** [`itiwpt`](@ref), [`tidwt`](@ref)
"""
function tiwpt(x::AbstractVector{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)
    n = length(x)
    xw = Matrix{T}(undef, (n, 1<<L))
    xw[:,1] = x

    # TIWPT for L levels
    for d in 0:(L-1)
        nblock = 1<<d
        for b in 0:(nblock-1)
            # Parent index
            np = (1<<L)÷nblock
            i = b*np + 1
            # Children indices
            nc = np÷2
            j₁ = 2*b*nc + 1
            j₂ = (2*b+1)*nc + 1
            # Overwrite output of TIDWT directly onto xw
            xw[:,j₁], xw[:,j₂] = tidwt_step(xw[:,i], d, h, g)
        end
    end
    return xw
end

"""
    itiwpt(xw, wt)

Computes the inverse translation-invariant wavelet packet transform (iTIWPT) on `xw`.

# Arguments
- `xw::AbstractArray{T,2} where T<:Number`: TIDWT-transformed array.
- `wt::OrthoFilter`: Orthogonal wavelet filter.

# Returns
`::Vector{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIWPT
xw = tiwpt(x, wt)

# iTIWPT
x̂ = itiwpt(xw, wt)
```

**See also:** [`tiwpt`](@ref), [`itidwt`](@ref)
"""
function itiwpt(xw::AbstractArray{T,2}, wt::OrthoFilter) where T<:Number
    # Sanity check
    n, m = size(xw)
    @assert isdyadic(m) || throw(ArgumentError("Number of columns of xw is not dyadic."))
    @assert m ≤ maxtransformlevels(n) ||
            throw(ArgumentError("Number of nodes in `xw` is more than possible number of nodes at any depth for signal of length `n`"))

    # Setup
    g, h = WT.makereverseqmfpair(wt, false)
    temp = copy(xw)                         # Temp. array to store intermediate outputs
    L = ndyadicscales(m)                    # Number of decompositions from xw

    # iTIWPT
    for d in (L-1):-1:0
        nblock = 1<<d
        for b in 0:(nblock-1)
            # Parent index
            np = (1<<L)÷nblock
            i = b*np + 1
            # Children indices
            nc = np÷2
            j₁ = 2*b*nc + 1
            j₂ = (2*b+1)*nc + 1
            # Overwrite output of iTIWPT directly onto temp
            temp[:,i] = itidwt_step(temp[:,j₁], temp[:,j₂], d, h, g)
        end
    end
    return temp[:,1]
end

# ========== Translation-Invariant WPD ==========
@doc raw"""
    tiwpd(x, wt[, L])

Computes `L` levels of translation-invariant wavelet packet decomposition (TIWPD) on `x`.

# Arguments
- `x::AbstractVector{T} where T<:Number`: Original signal, preferably of size 2ᴷ where ``K
  \in \mathbb{N}``.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of decomposition.

# Returns
`::Matrix{T}`: Output from TIWPD on `x`.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIWPD
xw = tiwpd(x, wt)
```

**See also:** [`itiwpd`](@ref), [`tiwpt`](@ref)
"""
function tiwpd(x::AbstractVector{T},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)
    n = length(x)                       # Size of signal
    n₀ = 1<<(L+1)-1                     # Total number of nodes
    n₁ = n₀ - (1<<L)                    # Total number of nodes excluding leaf nodes
    xw = Matrix{T}(undef, (n, n₀))      # Output allocation
    xw[:,1] = x

    # TIWPD
    for i in 1:n₁
        d = floor(Integer, log2(i))     # Current depth
        xw[:,left(i)], xw[:,right(i)] = tidwt_step(xw[:,i], d, h, g)
    end
    return xw
end

"""
    itiwpd(xw, wt[, L])
    itiwpd(xw, wt, tree)

Computes the inverse translation-invariant wavelet packet decomposition (iTIWPD) on `xw`.

# Arguments
- `xw::AbstractArray{T,2} where T<:Number`: TIDWT-transformed array.
- `wt::OrthoFilter`: Orthogonal wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(size(xw,1))`) Depth of inverse decomposition.
- `tree::BitVector`: Representation of binary tree used as basis for the signal `x`.

# Returns
`::Vector{T}`: Inverse transformed signal.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
x = generatesignals(:heavysine)
wt = wavelet(WT.haar)

# TIWPD
xw = tiwpd(x, wt)

# iTIWPD
x̂ = itiwpd(xw, wt)
```

**See also:** [`tiwpd`](@ref), [`itiwpt`](@ref)
"""
function itiwpd(xw::AbstractArray{T,2},
                wt::OrthoFilter,
                L::Integer = maxtransformlevels(size(xw,1))) where T<:Number
    # Sanity check
    @assert L ≤ maxtransformlevels(x) ||
            throw(ArgumentError("Too many transform levels (length(x) < 2ᴸ)"))
    @assert L ≥ 1 || throw(ArgumentError("L must be ≥ 1"))
    return itiwpd(xw, wt, maketree(size(xw,1), L, :full))
end

function itiwpd(xw::AbstractArray{T,2}, wt::OrthoFilter, tree::BitVector) where
                T<:Number
    # Sanity check
    @assert isvalidtree(xw[:,1], tree)

    # Setup
    g, h = WT.makereverseqmfpair(wt, false)
    tmp = copy(xw)

    # iTIWPD
    for (i, haschild) in Iterators.reverse(enumerate(tree))
        # Current node has child => Compute one step of inverse transform
        if haschild
            # Current depth
            d = floor(Integer, log2(i))
            # Parent node
            v = @view tmp[:,i]
            # Children nodes
            w₁ = @view tmp[:,left(i)]
            w₂ = @view tmp[:,right(i)]
            # Inverse transform
            itidwt_step!(v, w₁, w₂, d, h, g)
        end
    end
    return tmp[:,1]
end

end


text:
import Gtk.GLib: g_type_from_name
type DataFrameStoreImpl
  parent :: GObjectImpl
  private :: Ptr{Void}
end
type DataFrameStoreClassImpl
  parent :: GObjectClass
end
type DataFrameStore <: Gtk.GLib.GObject
  handle :: Ptr{DataFrameStoreImpl}
  #handle :: Ptr{GObject}
  stamp :: Int32
  df :: AbstractDataFrame
  function DataFrameStore()
    obj = g_object_new(DataFrameStoreGetType())
    n = new(obj.handle)
    println(n.handle)
    println(obj.handle)
    self = Gtk.gobject_move_ref(n, obj)
    dat = unsafe_load(self.handle)
    dat.private = pointer_from_objref(self)
    unsafe_store!(self.handle, dat)
    self.stamp = rand(Int32)
    self.df = DataFrame()
    self
  end
end
function DataFrameStore(df :: AbstractDataFrame)
  retval = DataFrameStore()
  retval.df = df
  retval
end

function DataFrameStoreTreeModelInit(iface :: Ptr{GtkTreeModelIface})
  i = unsafe_load(iface)
  i.get_flags = cfunction(DataFrameStoreGetFlags, Cint, (Ptr{DataFrameStoreImpl},))
  i.get_n_columns = cfunction(DataFrameStoreGetNColumns, Cint, (Ptr{DataFrameStoreImpl},))
  i.get_column_type = cfunction(DataFrameStoreGetColumnType, Gtk.GLib.GType, (Ptr{DataFrameStoreImpl}, Cint))
  i.get_iter = cfunction(DataFrameStoreGetIter, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}, GtkTreePath))
  i.get_path = cfunction(DataFrameStoreGetPath, Ptr{GtkTreeView}, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}))
  i.get_value = cfunction(DataFrameStoreGetValue, Void, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}, Cint, Ptr{Gtk.GLib.GValue}))
  i.iter_next = cfunction(DataFrameStoreIterNext, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}))
  i.iter_previous = cfunction(DataFrameStoreIterPrev, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}))
  i.iter_children = cfunction(DataFrameStoreIterChildren, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}, Ptr{GtkTreeIterImpl}))
  i.iter_has_child = cfunction(DataFrameStoreIterHasChild, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}))
  i.iter_n_children = cfunction(DataFrameStoreIterNChildren, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}))
  i.iter_nth_child = cfunction(DataFrameStoreIterNthChild, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}, Ptr{GtkTreeIterImpl}, Cint))
  i.iter_parent = cfunction(DataFrameStoreIterParent, Cint, (Ptr{DataFrameStoreImpl}, Ptr{GtkTreeIterImpl}, Ptr{GtkTreeIterImpl}))
  unsafe_store!(iface, i)
  return
end
DataFrameStoreType = g_type_from_name(:DataFrameStore) :: Gtk.GLib.GType
function DataFrameStoreGetType()
  #global DataFrameStoreType:: Gtk.GLib.GType
  global DataFrameStoreType :: Gtk.GLib.GType
  if DataFrameStoreType == 0
    info = GTypeInfo(
      UInt16(sizeof(DataFrameStoreClassImpl)),
      C_NULL,
      C_NULL,
      C_NULL,
      C_NULL,
      C_NULL,
      UInt16(sizeof(DataFrameStoreImpl)),
      UInt16(0),
      C_NULL,
      C_NULL
    )
    ifaceInfo = GInterfaceInfo(
      cfunction(DataFrameStoreTreeModelInit, Void, (Ptr{GtkTreeModelIface},)),
      C_NULL,
      C_NULL
    )
    DataFrameStoreType = g_type_register_static(g_type_from_name(:GObject), "DataFrameStore", info, Int32(0))
    modelType = ccall((:gtk_tree_model_get_type, Gtk.libgtk), Gtk.GLib.GType, ())
    g_type_add_interface_static(DataFrameStoreType, modelType, ifaceInfo)
  end
  DataFrameStoreType
end

function DataFrameStoreGetFlags(tree_model :: Ptr{DataFrameStoreImpl})
  Cint(3)
end
function DataFrameStoreGetNColumns(tree_model :: Ptr{DataFrameStoreImpl})
  println(tree_model)
  println(unsafe_load(tree_model))
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  Cint(length(self.df))
end
function DataFrameStoreGetColumnType(tree_model :: Ptr{DataFrameStoreImpl}, index :: Cint)
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  typ = eltypes(self.df)[index+1]
  JuliaTypeToGType(typ)

end

function DataFrameStoreGetIter(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl}, path :: GtkTreePath)
  index = unsafe_load(GAccessor.indices(path))
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  it = unsafe_load(iter)
  (rows, cols) = size(self.df)
  if index > length(self.df)
    it.stamp = 0
    unsafe_store!(iter, it)
    return Cint(false)
  else
    it.stamp = self.stamp
    it.user1 = index + 1
    unsafe_store!(iter, it)
    return Cint(true)
  end
end

function DataFrameStoreGetPath(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl})
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  it = unsafe_load(iter)
  ret :: Ptr{GtkTreeView} = C_NULL
  (rows, cols) = size(self.df)
  if it.stamp != self.stamp
    return ret
  elseif it.user1 > rows
    return ret
  end
  ret = ccall((:gtk_tree_path_new, Gtk.libgtk), Ptr{GtkTreePath}, ())
  ccall((:gtk_tree_path_append_index, Gtk.libgtk), Void, (GtkTreePath, Cint), ret, Cint(it.user1))
  return ret
end

function DataFrameStoreGetValue(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl}, col :: Cint, value :: Ptr{Gtk.GLib.GValue})
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  it = unsafe_load(iter)
  row = it.user1
  (rows, cols) = size(self.df)
  if col+1 > cols
    return
  end
  if it.stamp != self.stamp
    return
  end
  item = self.df[row, col+1]
  gv = gvalue(item)
  unsafe_store!(value, gv)
  return
end

function DataFrameStoreIterNext(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl})
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  it = unsafe_load(iter)
  (rows, cols) = size(self.df)
  if it.stamp != self.stamp
    return Cint(false)
  end
  it.user1 = it.user1 + 1
  if it.user1 > rows
    it.stamp = 0
    unsafe_store!(iter, it)
    return Cint(false)
  end
  unsafe_store!(iter, it)
  return Cint(true)
end
function DataFrameStoreIterPrev(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl})
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  it = unsafe_load(iter)
  if it.stamp != self.stamp
    return Cint(false)
  end
  it.user1 = it.user1 - 1
  if it.user1 <= 0
    it.stamp = 0
    unsafe_store!(iter, it)
    return Cint(false)
  end
  unsafe_store!(iter, it)
  return Cint(true)
end

function DataFrameStoreIterChildren(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl}, parent :: Ptr{GtkTreeIterImpl})
  it = unsafe_load(iter)
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  if parent != C_NULL
    it.stamp = 0
    unsafe_store!(iter, it)
    return Cint(false)
  end
  (rows, cols) = size(self.df)
  if rows > 0
    it.stamp = self.stamp
    it.user1 = 1
    unsafe_store(iter, it)
    return Cint(true)
  end
  it.stamp = 0
  unsafe_store!(iter, it)
  return Cint(false)
end

function DataFrameStoreIterHasChild(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl})
  return Cint(false)
end

function DataFrameStoreIterNChildren(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl})
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  (rows, cols) = size(self.df)
  if iter == C_NULL
    return Cint(rows)
  end
  it = unsafe_load(iter)
  if self.stamp != it.stamp
    return Cint(-1)
  end
  return Cint(0)
end

function DataFrameStoreIterNthChild(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl}, parent :: Ptr{GtkTreeIterImpl}, n :: Cint)
  if parent != C_NULL
    return Cint(false)
  end
  self = unsafe_pointer_to_objref(unsafe_load(tree_model).private)
  (rows, cols) = size(self.df)
  if n+1 < 1 || n+1 > rows
    return Cint(false)
  end
  it = unsafe_load(iter)
  it.stamp = self.stamp
  it.user1 = n+1
  unsafe_store!(iter, it)
  return Cint(true)
end
function DataFrameStoreIterParent(tree_model :: Ptr{DataFrameStoreImpl}, iter :: Ptr{GtkTreeIterImpl}, child :: Ptr{GtkTreeIterImpl})
  it = unsafe_load(iter)
  it.stamp = 0
  unsafe_store!(iter, it)
  return Cint(false)
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SingleMass <: Shaft
        H::Float64
        D::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of single mass shaft model. Typically represents the rotor mass.

# Arguments
- `H::Float64`: Rotor inertia constant in MWs/MVA, validation range: (0, nothing)
- `D::Float64`: Rotor natural damping in pu, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SingleMass <: Shaft
    "Rotor inertia constant in MWs/MVA"
    H::Float64
    "Rotor natural damping in pu"
    D::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SingleMass(H, D, ext=Dict{String, Any}(), )
    SingleMass(H, D, ext, [:δ, :ω], 2, InfrastructureSystemsInternal(), )
end

function SingleMass(; H, D, ext=Dict{String, Any}(), )
    SingleMass(H, D, ext, )
end

# Constructor for demo purposes; non-functional.
function SingleMass(::Nothing)
    SingleMass(;
        H=0,
        D=0,
        ext=Dict{String, Any}(),
    )
end

"""Get SingleMass H."""
get_H(value::SingleMass) = value.H
"""Get SingleMass D."""
get_D(value::SingleMass) = value.D
"""Get SingleMass ext."""
get_ext(value::SingleMass) = value.ext
"""Get SingleMass states."""
get_states(value::SingleMass) = value.states
"""Get SingleMass n_states."""
get_n_states(value::SingleMass) = value.n_states
"""Get SingleMass internal."""
get_internal(value::SingleMass) = value.internal


text:
#
#   BoostFractor
#
#   V: 2019-11-29
#
#   Stefan Knirck
#
#
#   This file defines common interfaces for all algorithms, such as (x,y)
#   coordinate system, the SetupBoundaries structure, Propagators, ...
#

export SetupBoundaries, SeedSetupBoundaries, propagator, propagator1D, CoordinateSystem, SeedCoordinateSystem, wavelength

function wavelength(frequency::Float64)
    speed_of_light = 299792458. # [m]
    return speed_of_light / frequency
end

struct CoordinateSystem
    # Real Space Coordinates
    X::Array{Float64,1}
    Y::Array{Float64,1}
    # K space Coordinates
    kX::Array{Float64,1}
    kY::Array{Float64,1}
end

function SeedCoordinateSystem(;X = -0.5:0.01:0.5, Y = -0.5:0.01:0.5) # units [m]
    kX = get_kspace_coords(X)
    kY = get_kspace_coords(Y)
    return CoordinateSystem(X, Y, kX, kY)
end

"""
    get_kspace_coords(RealSpaceCoords)

Calculate k space coordinate system from real space one. Helper function for SeedCoordinateSystem.
"""
function get_kspace_coords(RealSpaceCoords)
    minimum_Kspace = 2*pi/(maximum(RealSpaceCoords)*2)
    maximum_Kspace = minimum_Kspace * (length(RealSpaceCoords)-1)/2
    coordsKspace = -maximum_Kspace:minimum_Kspace:maximum_Kspace
    return coordsKspace
end

@doc raw"""
# Summary
    mutable struct SetupBoundaries <: Any

Define properties of dielectric boundaries. Coordinate system?

# Fields:
- `distance::Array{Float64,1}` ```> 0```: Distance in z direction to boundary
- `r::Array{Complex{Float64},1}` ```[0, 1]```: Boundary reflection coefficient for right-propagating wave
- `eps::Array{Complex{Float64},1}`: Dielectric permittivity to the right of each boundary"
- `relative_tilt_x` ```> 0```: Tilt in x direction [rad?]
- `relative_tilt_y` ```> 0```: Tilt in y direction [rad?]
- `relative_surfaces::Array{Complex{Float64},3}` ```?```: Surface roughness. z offset (1st dim) at x,y (2nd, 3rd dims)
"""
mutable struct SetupBoundaries
    distance::Array{Float64,1} # = [15e-3, 5e-3,0]
    r::Array{Complex{Float64},1}   # = [1,-0.5,0.5,0]
    eps::Array{Complex{Float64},1}   # = [1,9,1]
    relative_tilt_x # = [0,0]
    relative_tilt_y # = [0,0]
    relative_surfaces::Array{Complex{Float64},3} # = [z, x,y ]
    # etc.
end

## Convenient tools ################################################################################
@doc raw"""
    SeedSetupBoundaries(diskno=3)

Initialize `mutable struct SetupBoundaries` with sensible values.

# Arguments
- `diskno::Int` ```> 0```: Number of dielectric discs
"""
function SeedSetupBoundaries(coords::CoordinateSystem; diskno=3, distance=nothing, epsilon=nothing, relative_tilt_x=zeros(2*diskno+2), relative_tilt_y=zeros(2*diskno+2), relative_surfaces=zeros(2*diskno+2 , length(coords.X), length(coords.Y)))

    # Initialize SetupBoundaries entries with default values given diskno. Rest was already defined in function definition.
    if distance === nothing # [m]
        distance = [0.0]
        append!(distance, [ x % 2 == 1 ? 8e-3 : 1e-3 for x in 1:2*(diskno) ])
        append!(distance, 0e-3) #8e-3)
    end

    if epsilon === nothing
        epsilon = [NaN]
        append!(epsilon, [ x % 2 == 1.0 ? 1.0 : 9.0 for x in 1:2*(diskno) ])
        append!(epsilon, 1.0)
    end
    epsilon = Array{Complex{Float64}}(epsilon)
    epsilon[imag.(epsilon) .== 0.0] .= real.(epsilon[imag.(epsilon) .== 0.0]) .- 0.0im#Make sure that a real epsilon is casted to Re(eps)-0.0im.

    reflectivities = complex([1.0])
    R = [(sqrt(epsilon[i-1]) - sqrt(epsilon[i])) / (sqrt(epsilon[i-1]) + sqrt(epsilon[i])) for i in 3:length(epsilon)]
    append!(reflectivities, R)

    # Check if initialization was self-consistent
    length(distance) == length(reflectivities)+1 == length(epsilon) == length(relative_tilt_x) == length(relative_tilt_y) == size(relative_surfaces, 1) || throw(DimensionMismatch("the arrays in your SetupBoundaries objects don't fit together!"))

    return SetupBoundaries(distance, Array{Complex{Float64}}(reflectivities), epsilon, relative_tilt_x, relative_tilt_y, relative_surfaces)
end

## The heart of it #################################################################################

## Propagators ########################################################################################
#TODO: propagator and propagator NoTilts are in-place, julia convention is name! instead of name
"""
    propagator(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda)

Do the FFT of E0 on a disk, propagate the beam a given distance and do the iFFT.
Note that this method is in-place. If it should be called more than one time on the
same fields, use propagator(copy(E0), ...).

Assume: Tilt is small, additional phase is obtained by propagating all fields just
with k0 to the tilted surface (only valid if diffraction effects are small).

# Arguments
- `E0::Array{Float64,2}`: Electric field before propagation
- `dz`: Distance propagated in z direction
- `diskR`: Radius of discs
- `eps`: Dielectric permittivity
- `tilt_x`: Disc tilt in x direction
- `tilt_y`: Disc tilt in y direction
- `surface`: Surface roughness of disc (only at end of propagation)
- `lambda`: Wavelength of electric field

See also: [`propagatorMomentumSpace`](@ref)
"""
function propagator(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda, coords::CoordinateSystem)
    #For real eps we need to make sure it's represented as real(eps)-0.0im so that sqrt(k0^2- Kx^2 - Ky^2)[propagatorNoTilts] has negative imaginary part for
    #Kx^2 + Ky^2>k0 and these modes propagate with a loss:
    imag(eps) == 0.0 ? eps = real(eps) - 0.0im : nothing
    k0 = 2*pi/lambda*sqrt(eps)
    # Call the propagator and add a phase imposed by the tilt
    E0 = propagatorNoTilts(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda, coords)
    # Tilts:
    E0 .*= [exp(-1im*k0*tilt_x*x) * exp(-1im*k0*tilt_y*y) for x in coords.X, y in coords.Y]
    # More general: Any surface misalignments:
    E0 .*= exp.(-1im*k0*surface) #(the element wise (exp.) is important, otherwise "surface" is treated as a matrix!)
    return E0
end

"""
    propagatorNoTilts(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda)

Wrapped by [`propagator`](@ref). Go there for documentation. Tilt arguments to be
compatible with other propagators.
"""
function propagatorNoTilts(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda, coords::CoordinateSystem)
    # Diffract at the Disk. Only the disk is diffracting.
    E0 .*= [abs(x^2 + y^2) < diskR^2 for x in coords.X, y in coords.Y]
    # FFT the E-Field to spatial frequencies
    # fft! and ifft! in the current release (1.2.2) only work with type ComplexF32 and ComplexF64
    # fft and ifft seem more stable
    FFTW.fft!(E0)
    E0 = FFTW.fftshift(E0)

    # TODO: If maximum k is higher than k0, then it is not defined
    #       what happens with this mode
    #       We should give a warning and handle this here
    #       At the moment the script will just also propagate with a loss for those components
    # Propagate through space
    k0 = 2*pi/lambda*sqrt(eps)
    k_prop = [sqrt(k0^2 - Kx^2 - Ky^2) for Kx in coords.kX, Ky in coords.kY]

    E0 = E0 .* exp.(-1im*k_prop*dz)
    # Backtransform
    E0 = FFTW.ifftshift(E0)
    FFTW.ifft!(E0)
    return E0
end

"""
    propagatorMomentumSpace(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda)

Propagator that assumes E0 is already in momentum space. Mix between [`propagator`](@ref)
and [`propagatorNoTilts`](@ref). Go to [`propagator`](@ref) for documentation.
"""
function propagatorMomentumSpace(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda, coords::CoordinateSystem)
    # Propagate through space
    k0 = 2*pi/lambda*sqrt(eps)
    k_prop = [sqrt(k0^2 - Kx^2 - Ky^2) for Kx in coords.kX, Ky in coords.kY]
    E0 = E0 .* exp.(-1im*k_prop*dz)

    # Transform to position space
    E0 = FFTW.ifftshift(E0)
    FFTW.ifft!(E0)

    # Diffract at the Disk. Only the disk is diffracting.
    E0 .*= [abs(x^2 + y^2) < diskR^2 for x in coords.X, y in coords.Y]

    # Kick (tilt)
    if tilt_x != 0 || tilt_y != 0
        E0 .*= [exp(-1im*k0*tilt_x*x) * exp(-1im*k0*tilt_y*y) for x in coords.X, y in coords.Y]
    end

    # FFT the E-Field to spatial frequencies / momentum space
    FFTW.fft!(E0)
    E0 = FFTW.fftshift(E0)

    return E0
end

"""
    propagator1D(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda)

This propagator just does the phase propagation. Go to [`propagator`](@ref)
for documentation. 3D arguments to be compatible with other propagators.
"""
function propagator1D(E0, dz, diskR, eps, tilt_x, tilt_y, surface, lambda, coords::CoordinateSystem)
    # Version of the propagator without the fft
    # should be faster and easy to check consistency with 1D calc

    # Propagate through space
    k0 = 2*pi/lambda*sqrt(eps)
    e1 = E0.*exp(-1im*k0*dz)
    return e1

end

#TODO: This function is no longer supported. reflectivity_transmissivity_1d does not exist; DiskDefinition has merged into SetuoBoundaries.
@doc raw"""
    initialize_reflection_transmission(freq::Float64, bdry::SetupBoundaries, disk::DiskDefiniton)

OUTDATED! Does not work, do not use!
Calculate reflection and transmission coefficients.

# Arguments
- `freq::Float64` ```> 0```: Frequency of EM radiation
- `bdry::SetupBoundaries`: Properties of dielectric boundaries
- `disk::DiskDefiniton`: Properties of dielectric discs
"""
# function initialize_reflection_transmission(freq::Float64, bdry::SetupBoundaries, coords::CoordinateSystem)#, disk::DiskDefiniton)
#     if disk === nothing
#         # Iniatilize reflection coefficients according to epsilon
#         r_left = ones(length(bdry.eps))
#         r_left[1] = -1
#         for i in 2:length(bdry.eps)
#             # The reflectivity at this point
#             r_left = (sqrt(bdry.eps[i-1])-sqrt(bdry.eps[i]))/(sqrt(bdry.eps[i])+sqrt(bdry.eps[i-1]))
#         end
#         r_right = -r_left
#         t_left = 1. + r_left
#         t_right = 1 .+ r_right
#     else
#         # Initailize reflection coefficients according to disk model
#         ref, trans = reflectivity_transmissivity_1d(freq, disk.thickness)
#         r_left = ones(length(bdry.eps),length(coords.kX),length(coords.kY)).*ref
#         r_right = r_left
#         t_left = ones(length(bdry.eps),length(coords.kX),length(coords.kY)).*trans
#         t_right = t_left
#     end
#     return r_left, r_right, t_left, t_right
# end


text:
# This file contains the tests for the `Both` module

module TestBoth

using PDFHighlights
using Test

# Print the header
println("\e[1;32mRUNNING\e[0m: TestBoth.jl")

const pdf = joinpath(@__DIR__, "..", "pdf", "TestPDF.pdf")

@testset "get_authors" begin

    dir = joinpath(@__DIR__, "..")

    @test get_authors(dir) == String["Pavel Sobolev",]

    csv = "oof.csv"
    import_highlights(csv, pdf; quiet=true)

    @test get_authors(csv) == repeat(["Pavel Sobolev",], 7)

    isfile(csv) && rm(csv)

    @test_throws PDFHighlights.Internal.Exceptions.NotCSVorDir("oof") get_authors("oof")

end

@testset "get_highlights" begin

    # With concatenation
    @test get_highlights(pdf) == String[
        "Highlight 1",
        "Highlight 2 Highlight 3",
        "Highlight 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    # Without concatenation
    @test get_highlights(pdf; concatenate=false) == String[
        "Highlight 1",
        "Highlight 2",
        "Highlight 3",
        "High-",
        "light 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    dir = joinpath(@__DIR__, "..")

    # With concatenation
    @test get_highlights(dir) == String[
        "Highlight 1",
        "Highlight 2 Highlight 3",
        "Highlight 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    # Without concatenation
    @test get_highlights(dir; concatenate=false) == String[
        "Highlight 1",
        "Highlight 2",
        "Highlight 3",
        "High-",
        "light 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    csv = "oof.csv"
    import_highlights(csv, pdf; quiet=true)

    @test get_highlights(csv) == String[
        "Highlight 1",
        "Highlight 2 Highlight 3",
        "Highlight 4",
        "Highhighlight 5",
        "6th Highhigh light-",
        "High light 7",
        "8th Highlight-",
    ]

    isfile(csv) && rm(csv)

    @test_throws(
        PDFHighlights.Internal.Exceptions.NotCSVorPDForDir("oof"),
        get_highlights("oof")
    )

end

@testset "get_titles" begin

    dir = joinpath(@__DIR__, "..")

    @test get_titles(dir) == String["A dummy PDF for tests",]

    csv = "oof.csv"
    import_highlights(csv, pdf; quiet=true)

    @test get_titles(csv) == repeat(["A dummy PDF for tests",], 7)

    isfile(csv) && rm(csv)

    @test_throws PDFHighlights.Internal.Exceptions.NotCSVorDir("oof") get_titles("oof")

end

end


text:
# Parametric type examples

abstract type Asset end
abstract type Investment <: Asset end
abstract type Equity <: Investment end

struct Stock <: Equity
    symbol::String
    name::String
end

struct StockHolding{T <: Real}
    stock::Stock
    quantity::T
end

#= REPL

julia> stock = Stock("AAPL", "Apple, Inc.");

julia> holding = StockHolding(stock, 100)
StockHolding{Int64}(Stock("AAPL", "Apple, Inc."), 100)

julia> holding = StockHolding(stock, 100.00)
StockHolding{Float64}(Stock("AAPL", "Apple, Inc."), 100.0)

julia> holding = StockHolding(stock, 100 // 3)
StockHolding{Rational}(Stock("AAPL", "Apple, Inc."), 100//3)
=#

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

struct StockHolding2{T <: Real, P <: AbstractFloat}
    stock::Stock
    quantity::T
    price::P
    marketvalue::P
end

#= REPL
julia> holding = StockHolding2(stock, 100, 180.00, 18000.00)
StockHolding2{Int64,Float64}(Stock("AAPL", "Apple, Inc."), 100, 180.0, 18000.0)

julia> holding = StockHolding2(stock, 100, 180.00, 18000)
ERROR: MethodError: no method matching StockHolding2(::Stock, ::Int64, ::Float64, ::Int64)
Closest candidates are:
  StockHolding2(::Stock, ::T<:Real, ::P<:AbstractFloat, ::P<:AbstractFloat) where {T<:Real, P<:AbstractFloat} at REPL[77]:2
=#

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

abstract type Holding{P} end

mutable struct StockHolding3{T, P} <: Holding{P}
    stock::Stock
    quantity::T
    price::P
    marketvalue::P
end

mutable struct CashHolding{P} <: Holding{P}
    currency::String
    amount::P
    marketvalue::P
end

#= REPL
julia> stock = Stock("AAPL", "Apple, Inc.");

julia> certificate_in_the_safe = StockHolding3(stock, 100, 180.00, 18000.00)
StockHolding3{Int64,Float64}(Stock("AAPL", "Apple, Inc."), 100, 180.0, 18000.0)

julia> StockHolding3{Int64,Float64} <: Holding{Float64}
true

julia> certificate_in_the_safe isa Holding{Float64}
true

julia> Holding{Float64} <: Holding
true

julia> Holding{Int} <: Holding
true
=#


text:
# This file provides reproducible code for the extraction of SigmaRidgeRegression.CLLData
# from the R/Bioconductor MOFA package.

# A R installation (it will be called through `RCall`) is required
# with an installation of the `MOFAdata` package.
# This package may be installed from within `R` as follows:

# ```r
# if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("MOFAdata")
#```

using JLD2
using RCall

R"""
	data("CLL_data", package="MOFAdata")

	# use methylation data, gene expression data and drug responses as predictors
	CLL_data <- CLL_data[1:3]
	CLL_data <- lapply(CLL_data,t)
	ngr <- sapply(CLL_data,ncol)
	CLL_data <- Reduce(cbind, CLL_data)

	#only include patient samples profiles in all three omics
	CLL_data2 <- CLL_data[apply(CLL_data,1, function(p) !any(is.na(p))),]
	dim(CLL_data2)

	# prepare design matrix and response
	X <- CLL_data2[,!grepl("D_002", colnames(CLL_data))]
	y <- rowMeans(CLL_data2[,grepl("D_002", colnames(CLL_data))])
	annot <- rep(1:3, times = ngr-c(5,0,0)) # group annotations to drugs, meth and RNA
	ngr_prime <- ngr-c(5,0,0)
"""

# run with seed from Velten & Huber
R"""
set.seed(9876)
foldid <- sample(rep(seq(10), length=nrow(X)))
"""

@rget foldid
@rget X
@rget y
@rget ngr_prime

cll_data = (X = X, y = y, ngr = Int.(ngr_prime), foldid = foldid)
JLD2.@save "cll_data.jld2" {compress=true} cll_data


text:
# this is an experimental version that allocates less. It is
# about 30% faster, for a very simple likelihood, but the
# time is so short that it's not really important. For more complex
# likelihoods, it seems likely that making the likelihood computation fast is more beneficial than making the sampler fast (?)
using Distributions

function Prior!(p, θ)
    p .= pdf.(LogNormal(1.0,1.0), θ[1])
end

function lnL!(ll, θ, y)
    ll .= sum(logpdf.(Exponential(θ[1]), y))
end

function Proposal!(θᵗ, θᶜ, tuning)
    θᵗ .= rand(LogNormal(log(θᶜ[1]),tuning))
end

struct model
    data::Array{Float64}
    prior!::Function
    lnL!::Function
    proposal!::Function
end

function mcmc(θᶜ, model, reps, burnin)
    tuning = 0.5
    θᵗ = copy(θᶜ) # θᵗ params
    chain = copy(θᶜ) # need defined at this scope
    ℒᵗ = [-Inf]
    ℒᶜ = [-Inf]
    pt = [-Inf]
    pc = [1.]
    Prior!(pc, θᶜ)
    model.lnL!(ℒᶜ, θᶜ, model.data) # get initial likelihood
    for i = 1:reps + burnin
        model.proposal!(θᵗ, θᶜ, tuning)
        model.lnL!(ℒᵗ, θᵗ, model.data)
        Prior!(pt, θᵗ)
        # if accepted, update the current quantities
        if rand() < exp(ℒᵗ[1] - ℒᶜ[1])*pt[1]/pc[1]
            θᶜ .= θᵗ
            pc .= pt
            ℒᶜ .= ℒᵗ
        end
        # if out of burning, push current to chain
        i > burnin + 1 ? push!(chain, θᶜ[1]) :
            (i == burnin + 1 ? chain = copy(θᶜ) : nothing)
    end
    acceptance = (size(unique(chain)) ./ size(chain))[1]
return chain, acceptance
end

function main()
y = rand(Exponential(3.0),30)
θ = [rand(LogNormal(1.0, 1.0))]   # initial value draw from prior
m = model(y, Prior!, lnL!, Proposal!)
chain, acceptance = mcmc(θ, m, 100000, 10000)
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Type{NamedTuple{(:X, :Z, :W), T} where T<:Tuple},Core.Tuple{Base.AbstractVector{T} where T, Base.AbstractVector{T} where T, Base.AbstractVector{T} where T}})   # time: 0.0056522
end


text:
using Documenter, NelderMead

makedocs(sitename="NelderMead.jl")

deploydocs(
    repo = "github.com/jwscook/NelderMead.jl.git",
)


text:
const global True  = ConstantTerm("⊤")
const global False = ConstantTerm("⊥")

"""
   FOOL2FOL!(a::AbstractLogicTerm, e::Env)

Convert a FOOL logic expression with associated Environment `e` to first order form. This
may modify the environment by introducing new typed `FunctionTerm`s. See [A First Class
Boolean Sort in First-Order Theorem Proving and TPTP](https://arxiv.org/pdf/1505.01682.pdf)
for details.
"""
function FOOL2FOL!(a, e=Env())
  # f2f! replaces all FOOL formula with FOL formula and adds new terms
  function f2f!(fp::Union{FunctionTerm,PredicateTerm}, defs, e::Env, boolContext)
    newArgs = similar(fp.args, Union{Variable,FunctionTerm})
    for (i,a) in enumerate(fp.args)
      a = f2f!(a, defs, e, false)  # recursively apply f2f!
      if typeof(a) <: Union{PredicateTerm,JunctionTerm,QuantifierTerm,NegationTerm}
        fv = freeVar(a)
        F = FunctionTerm(newFunctor(e), fv);  newArgs[i] = F
        eq = equalTerm(F, True)
        ts = getindex.(Ref(e), getproperty.(fv, :name))  # get argument types
        e[F.name.name] = tuple(getindex.(ts,1)..., BoolType())
        newFormula = AndTerm(OrTerm(NegationTerm(a), eq), OrTerm(NegationTerm(eq), a))
        push!(defs, isempty(fv) ? newFormula : AQuantifierTerm(fv, newFormula))
      else
        newArgs[i] = a
      end
    end
    rootType(fp)(fp.name, newArgs)
  end
  f2f!(n::NegationTerm, defs, e::Env, boolContext) = NegationTerm(f2f!(n.scope, defs, e, true))
  function f2f!(j::JunctionTerm, defs, e::Env, boolContext)
    newJuncts = map(x -> isa(x, Variable) ? equalTerm(x, True) : x, f2f!.(j.juncts, Ref(defs), Ref(e), true))
    rootType(j)(convert(Vector{typejoin(typeof.(newJuncts)...)}, newJuncts))
  end
  f2f!(q::QuantifierTerm, defs, e::Env, boolContext) = rootType(q)(q.variables, f2f!(q.scope, defs, e, true))
  f2f!(v::Variable, defs, e::Env, boolContext) = boolContext ? equalTerm(v, True) : v

  defs = Vector{SententialTerm}(undef, 0)
  ret = f2f!(a, defs, e, true)
  if isempty(defs)
    ret
  else
    push!(defs, ret)
    e[True.name.name] = BoolType();  e[False.name.name] = BoolType()
    AndTerm(convert(Vector{typejoin(typeof.(defs)...)}, defs))
  end
end

"""
    pullCommonQuantifiers(s::SententialTerm)

Pull out common quantifiers if possible; this is only useful before variable
renaming.
"""
pullCommonQuantifiers(p::PredicateTerm) = p
pullCommonQuantifiers(n::NegationTerm) = NegationTerm(pullCommonQuantifiers(n.scope))
function pullCommonQuantifiers(j::JunctionTerm)
  newJuncts = pullCommonQuantifiers.(j.juncts)
  QJ = pairedQuantifierType(j)
  if all(rootType.(newJuncts) .== QJ)
    commonVars = reduce((cs,vs) -> intersect(vs,cs), getproperty.(newJuncts[2:end], :variables);
    init = newJuncts[1].variables)
    if ~isempty(commonVars)
      newScope = [begin
                    diffVars = setdiff(j.variables, commonVars)
                    isempty(diffVars) ? j.scope :  QJ(diffVars, j.scope)
                  end for j in newJuncts]
      QJ(commonVars, rootType(j)(newScope))
    else
      rootType(j)(newJuncts...)
    end
  else
    rootType(j)(newJuncts...)
  end
end
pullCommonQuantifiers(q::QuantifierTerm) = rootType(q)(q.variables, pullCommonQuantifiers(q.scope))

"""
    simplify(s::SententialTerm)

Simplify the structure of s by removing double negations, aggregating nested
ands/ors and collecting quantified variabes in nested quantifiers

"""
simplify(p::Union{Variable,PredicateTerm}) = p
simplify(a::J) where {J<:JunctionTerm{<:LiteralTerm}} = a
function simplify(j::JunctionTerm)
  if length(j.juncts) == 1
    j.juncts[1]
  else
    J = rootType(j);  newJuncts = SententialTerm[]
    for c in j.juncts
      cs = simplify(c)
      rootType(cs) == J ? append!(newJuncts, cs.juncts) : push!(newJuncts, cs)
    end
  end
  J(newJuncts...)
end
function simplify(n::NegationTerm)
  ns = simplify(n.scope)
  isa(ns, NegationTerm) ? ns.scope : NegationTerm(ns)
end
function simplify(q::QuantifierTerm)
  Q = rootType(q);  s = simplify(q.scope)
  rootType(s) == Q ? Q(union(q.variables, s.variables), s.scope) : Q(q.variables, s)
end

"""
    negationNormalForm(s::SententialTerm)

Obtain the negation normal form by pushing negations inward.
"""
negationNormalForm(p::PredicateTerm)= p
negationNormalForm(a::JunctionTerm) = rootType(a)(negationNormalForm.(a.juncts))
negationNormalForm(q::QuantifierTerm) = rootType(q)(q.variables, negationNormalForm(q.scope))
function negationNormalForm(n::NegationTerm)
  nnf(p::PredicateTerm) = NegationTerm(p)
  nnf(n::NegationTerm) = negationNormalForm(n.scope)
  nnf(a::OrTerm) = AndTerm(negationNormalForm.(NegationTerm.(a.juncts)))
  nnf(a::AndTerm) = OrTerm(negationNormalForm.(NegationTerm.(a.juncts)))
  nnf(q::EQuantifierTerm) = AQuantifierTerm(q.variables, negationNormalForm(NegationTerm(q.scope)))
  nnf(q::AQuantifierTerm) = EQuantifierTerm(q.variables, negationNormalForm(NegationTerm(q.scope)))

  nnf(n.scope)
end

"""
    renameVars(s::SententialTerm, e=Env())

Standardize the names of all quantified variables so they are unique. Use new variable names
drawn from the environment `e`.
"""
function renameVars(s::SententialTerm, e=Env())
  rename_(p::PredicateTerm, seen, replacements) =
    PredicateTerm(p.name, [isa(a,Variable) ? get(replacements, a, a) : a for a in p.args])
  rename_(n::NegationTerm, seen, replacements) = NegationTerm(rename_(n.scope, seen, replacements))
  rename_(j::JunctionTerm, seen, replacements) = rootType(j)(rename_.(j.juncts, Ref(seen), Ref(replacements)))
  function rename_(q::QuantifierTerm, seen, replacements)
    newVars = [v ∈ seen ? replacements[v] = newVariable(v, e) : v for v in q.variables]
    rootType(q)(newVars, rename_(q.scope, union!(seen, Set(newVars)), replacements))
  end

  rename_(s, Set{Variable}(), Dict{Variable,Variable}())
end

"""
   pullQuantifiers(s::SententialTerm)

Pull quantifiers outwards assuming that s has had it's variables standardized.
"""
pullQuantifiers(p::PredicateTerm) = p
pullQuantifiers(n::NegationTerm) = NegationTerm(pullQuantifiers(n.scope))
function pullQuantifiers(j::JunctionTerm)
  function pullThroughJunction(js, J) # pull quantifiers through junction J
    es, as = isa.(js, EQuantifierTerm), isa.(js, AQuantifierTerm)
    if any(es)  # pull existential quantifiers first to simplify Skolemization
      EVars = collect(Iterators.flatten(getproperty.(js[es], :variables)))
      scope = map((e,j) -> e ? j.scope : j, es, js)
      EQuantifierTerm(EVars, pullThroughJunction(scope, J))
    elseif any(as)
      AVars = collect(Iterators.flatten(getproperty.(js[as], :variables)))
      scope = map((a,j) -> a ? j.scope : j, as, js)
      AQuantifierTerm(AVars, pullThroughJunction(scope, J))
    else
      J(js)
    end
  end

  pullThroughJunction(pullQuantifiers.(j.juncts), rootType(j))
end
pullQuantifiers(q::QuantifierTerm) = rootType(q)(q.variables, pullQuantifiers(q.scope))

"""
    skolemize(s::SententialTerm, e=Env())

Remove all existential quantifiers from `s` (which is assumed to be in prenex normal form)
by Skolemization. We use the environment `e` to define new variables/functors.
"""
function skolemize(s::SententialTerm, e=Env())
  skolemize_(v::Variable, replacements, _) = get(replacements, v, v)
  skolemize_(fp::Union{FunctionTerm,PredicateTerm}, replacements, previousAVars) =
    rootType(fp)(fp.name, skolemize_.(fp.args, Ref(replacements), Ref(previousAVars)))
  skolemize_(n::NegationTerm, replacements, previousAVars) =
    NegationTerm(skolemize_(n.scope, replacements, previousAVars))
  skolemize_(j::JunctionTerm, replacements, previousAVars) =
    rootType(j)(skolemize_.(j.juncts, Ref(replacements), Ref(previousAVars)))
  skolemize_(a::AQuantifierTerm, replacements, previousAVars) =
    AQuantifierTerm(a.variables, skolemize_(a.scope, replacements, append!(previousAVars, a.variables)))
  function skolemize_(eq::EQuantifierTerm, replacements, previousAVars)
    setindex!.(Ref(replacements), FunctionTerm.(newFunctor.(eq.variables, e), Ref(previousAVars)), eq.variables)
    skolemize_(e.scope, replacements, previousAVars)
  end

  skolemize_(s, Dict{Variable,FunctionTerm}(), Variable[])
end

"""
    flatten(ls)

Utility to flatten a list of lists into a list.
"""
flatList(ls) = collect(Iterators.flatten(ls))

using IterTools

"""
    distributeOrOverAnd(s::SententialTerm)

Distribute ∨ over ∧ assuming we are in prenex normal form. This is used when
converting a sentence `s` to clausal form.
"""
distributeOrOverAnd(l::LiteralTerm) = AndTerm(OrTerm([l]))
function distributeOrOverAnd(a::AndTerm)
  orTerms = mapreduce(j->distributeOrOverAnd(j).juncts, vcat, a.juncts; init=OrTerm[])
  AndTerm(orTerms)
end
function distributeOrOverAnd(o::OrTerm)
  clauseLiterals = reduce(o.juncts; init=Vector{Vector{LiteralTerm}}[]) do l,j
    jJuncts = getproperty.(distributeOrOverAnd(j).juncts, :juncts)
    push!(l, jJuncts)
  end
  clauses = [OrTerm(flatList(c)) for c in IterTools.ivec(Iterators.product(clauseLiterals...))]
  AndTerm(clauses)
end
distributeOrOverAnd(q::QuantifierTerm) = rootType(q)(q.variables, distributeOrOverAnd(q.scope))

"""
    prenexNormalForm(s::SententialTerm, e=Env())

Convert `s` to prenex normal form. Use the environment `e` to define new variables/functors.
"""
prenexNormalForm(s, e=Env()) = s |> negationNormalForm |> pullCommonQuantifiers |>
  Base.Fix2(renameVars, e) |> pullQuantifiers |> distributeOrOverAnd |> simplify

"""
    skolemNormalForm(s::SententialTerm, e=Env())

Return the Skolem normal form of `s` using the environment `e` to define new
variables/functors. The result is in CNF form.
"""
skolemNormalForm(s::SententialTerm, e=Env()) = s |> Base.Fix2(prenexNormalForm, e) |> skolemize |> simplify

"""
    stripQuantifiers(s::SententialTerm)

Remove all quantifiers from sentence s. This is useful when constructing the
universally quantified conjunction normal form.
"""
stripQuantifiers(p::PredicateTerm) = p
stripQUantifiers(n::NegationTerm) = NegationTerm(stripQuantifiers(n.scope))
stripQuantifiers(j::JunctionTerm) = rootType(j)(stripQuantifiers.(j.juncts))
stripQuantifiers(q::QuantifierTerm) = stripQuantifiers.(q.scope)