text:
"""
    gdaltranslate(indata, opts=String[]; dest="/vsimem/tmp", kwargs...)

Convert raster data between different formats and other operations also provided by the GDAL
'gdal_translate' tool. Namely sub-region extraction and resampling.
The `kwargs` options accept the GMT region (-R), increment (-I), target SRS (-J). Any of the keywords
`outgrid`, `outfile` or `save` = outputname options to make this function save the result in disk
in the file 'outputname'. The file format is picked from the 'outputname' file extension.
When no output file name is provided it returns a GMT object (either a grid or an image, depending
on the input type). To force the return of a GDAL dataset use the option `gdataset=true`.

- `indata`: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset
- `opts`:   List of options. The accepted options are the ones of the gdal_translate utility.
            This list can be in the form of a vector of strings, or joined in a simgle string.
- `kwargs`: Besides what was mentioned above one can also use `meta=metadata`, where `metadata`
            is a string vector with the form "NAME=...." foe each of its elements. This data
            will be recognized by GDAL as Metadata.

### Returns
A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).
"""
function gdaltranslate(indata, opts=String[]; dest="/vsimem/tmp", kwargs...)
	helper_run_GDAL_fun(gdaltranslate, indata, dest, opts, "", kwargs...)
end
# ---------------------------------------------------------------------------------------------------

"""
    gdalwarp(datasets::Vector{Dataset}, options=String[]; dest="/vsimem/tmp", kw...)

Image reprojection and warping function.

### Parameters
* `datasets` The list of input datasets.
* `options` List of options (potentially including filename and open
	options). The accepted options are the ones of the gdalwarp utility.
* `kw` are kwargs that may contain the GMT region (-R), proj (-J), inc (-I) and `save=fname` options

### Returns
A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).
"""
function gdalwarp(indata, opts=String[]; dest="/vsimem/tmp", kwargs...)
	helper_run_GDAL_fun(gdalwarp, indata, dest, opts, "", kwargs...)
end

# ---------------------------------------------------------------------------------------------------
"""
    gdaldem(dataset, method, options=String[]; dest="/vsimem/tmp", colorfile=name|GMTcpt, kw...)

Tools to analyze and visualize DEMs.

### Parameters
* `dataset` The source dataset.
* `method` the processing to apply (one of "hillshade", "slope",
    "aspect", "color-relief", "TRI", "TPI", "Roughness").
* `options` List of options (potentially including filename and open options).
    The accepted options are the ones of the gdaldem utility.

# Keyword Arguments
* `colorfile` color file (mandatory for "color-relief" processing, should be empty otherwise).
* `kw...` keyword=value arguments when `method` is hillshade.

### Returns
A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).
"""
function gdaldem(indata, method::String, opts::Vector{String}=String[]; dest="/vsimem/tmp", kwargs...)
	opts = gdal_opts2vec(opts)		# Guarantied to return a Vector{String}
	if (method == "hillshade")		# So far the only method that accept kwarg options
		d = GMT.KW(kwargs)
		band = ((val = GMT.find_in_dict(d, [:band])[1]) !== nothing) ? string(val) : "1"
		append!(opts, ["-compute_edges", "-b", band])
		if ((val = GMT.find_in_dict(d, [:scale])[1]) === nothing)
			if (isa(indata, GMT.GMTgrid) && (occursin("longlat", indata.proj4) || occursin("latlong", indata.proj4)) ||
											grdinfo(indata, C="n")[1].data[end] == 1)
				append!(opts, ["-s", "111120"])
			end
		else
			append!(opts, ["-s", string(val)])
		end
		((val = GMT.find_in_dict(d, [:zfactor :zFactor])[1]) !== nothing) && append!(opts, ["-z", string(val)])
		((val = GMT.find_in_dict(d, [:azim :azimuth])[1]) !== nothing) && append!(opts, ["-az", string(val)])
		((val = GMT.find_in_dict(d, [:elev :altitude])[1]) !== nothing) && append!(opts, ["-alt", string(val)])
		((val = GMT.find_in_dict(d, [:combined :combine])[1]) !== nothing) && append!(opts, ["-combined"])
		((val = GMT.find_in_dict(d, [:multi :multidir :multiDirectional])[1]) !== nothing) && append!(opts, ["-multidirectional"])
		((val = GMT.find_in_dict(d, [:igor])[1]) !== nothing) && append!(opts, ["-igor"])
		((val = GMT.find_in_dict(d, [:alg])[1]) !== nothing) && append!(opts, ["-alg", string(val)])
		((val = GMT.find_in_dict(d, [:Horn])[1]) !== nothing) && append!(opts, ["-alg", "Horn"])
		((val = GMT.find_in_dict(d, [:Zeven :Zevenbergen])[1]) !== nothing) && append!(opts, ["-alg", "ZevenbergenThorne"])
		helper_run_GDAL_fun(gdaldem, indata, dest, opts, method, d...)
	else
		helper_run_GDAL_fun(gdaldem, indata, dest, opts, method, kwargs...)
	end
end

"""
    function ogr2ogr(indata, opts=String[]; dest="/vsimem/tmp", kwargs...)

### Parameters
* `dataset` The source dataset.
* `options` List of options (potentially including filename and open
	options). The accepted options are the ones of the gdalwarp utility.
* `kw` are kwargs that may contain the GMT region (-R), proj (-J), inc (-I) and `save=fname` options

### Returns
A GMT dataset, or a GDAL dataset (or nothing if file was writen on disk).
"""
function gdalvectortranslate(indata, opts=String[]; dest="/vsimem/tmp", kwargs...)
	helper_run_GDAL_fun(gdalvectortranslate, indata, dest, opts, "", kwargs...)
end

# ---------------------------------------------------------------------------------------------------
function helper_run_GDAL_fun(f::Function, indata, dest::String, opts, method::String="", kwargs...)
	# Helper function to run the GDAL function under 'some protection' and returning obj or saving in file

	GMT.ressurectGDAL()				# Another black-hole plug attempt.
	opts = gdal_opts2vec(opts)		# Guarantied to return a Vector{String}
	d, opts, got_GMT_opts = GMT_opts_to_GDAL(opts, kwargs...)
	((val = GMT.find_in_dict(d, [:Vd])[1]) !== nothing) && println(opts)

	# For gdaldem color-relief we need a further arg that is the name of a cpt. So save one on disk
	_cmap = C_NULL
	if (f == gdaldem && ((cmap = GMT.find_in_dict(d, GMT.CPTaliases)[1])) !== nothing)
		_cmap = tempdir() * "/GMTtmp_cpt.cpt"
		if ( (isa(cmap, String) && (lowercase(splitext(cmap)[2][2:end]) == "cpt")) || isa(cmap, GMT.GMTcpt) )
			save_cpt4gdal(cmap, _cmap)	# GDAL pretend to recognise CPTs but it almost doesn't
		else
			_cmap = cmap				# Risky, assume it's something GDAL can read
		end
	end

	dataset, needclose = get_gdaldataset(indata, opts)
	((outname = GMT.add_opt(d, "", "", [:outgrid :outfile :save])) != "") && (dest = outname)
	default_gdopts!(dataset, opts, dest)	# Assign some default options in function of the driver and data type
	((val = GMT.find_in_dict(d, [:meta])[1]) !== nothing && isa(val,Vector{String})) &&
		Gdal.GDALSetMetadata(dataset.ptr, val, C_NULL)		# Metadata must be in the form NAME=.....

	CPLPushErrorHandler(@cfunction(CPLQuietErrorHandler, Cvoid, (UInt32, Cint, Cstring)))
	o = (method == "") ? f(dataset, opts; dest=dest, gdataset=true) : f(dataset, method, opts; dest=dest, gdataset=true, colorfile=_cmap)
	(o !== nothing && o.ptr == C_NULL) && @warn("$(f) returned a NULL pointer.")
	if (o !== nothing)
		# If not explicitly stated to return a GDAL datase, return a GMT type
		if (f == ogr2ogr)
			deletedatasource(o, "/vsimem/tmp")		# WTF do I need to do this?
			!(haskey(d, :gdataset) && !d[:gdataset]) && (o = gd2gmt(o))
		else
			n_bands = (got_GMT_opts && !haskey(d, :gdataset) && isa(o, AbstractRasterBand)) ? 1 : nraster(o)
			(!haskey(d, :gdataset)) && (o = gd2gmt(o, bands=collect(1:n_bands)))
		end
		haskey(d, :gdataset) && delete!(d, :gdataset)
	end
	(needclose) && GDALClose(dataset.ptr)
	CPLPopErrorHandler();
	(length(d) > 0) && println("Warning: the following options were not consumed in $f => ", keys(d))
	o
end

# Because the GDAL reconnaissance of GMT cpts is very very weak, we must re-write CPTs in a way that it can swallow
save_cpt4gdal(cpt::String, outname::String) = save_cpt4gdal(GMT.gmtread(cpt), outname)
function save_cpt4gdal(cpt::GMT.GMTcpt, outname::String)
	fid = open(outname, "w")
	try
		println(fid, "#COLOR_MODEL = RGB\n#")
		for k = 1:size(cpt.range,1)
			@printf(fid, "%.12g\t%.0f %.0f %.0f\t%.12g\t%.0f %.0f %.0f\n", cpt.range[k,1], cpt.cpt[k,1]*255, cpt.cpt[k,2]*255, cpt.cpt[k,3]*255, cpt.range[k,2], cpt.cpt[k,4]*255, cpt.cpt[k,5]*255, cpt.cpt[k,6]*255)
		end
	catch err
		@warn("Error $(err) saving a GDAL pleasing CPT file $(outname)")
	end
	close(fid)
end

# ---------------------------------------------------------------------------------------------------
function GMT_opts_to_GDAL(opts::Vector{String}, kwargs...)
	# Helper function to process some GMT options and turn them into GDAL syntax
	d = GMT.init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	((opt_R = GMT.parse_R(d, "")[1]) != "") && append!(opts, ["-projwin", split(opt_R[4:end], '/')[[1,4,2,3]]...])	# Ugly
	((opt_J = GMT.parse_J(d, "", " ")[1]) != " ") && append!(opts, ["-a_srs", opt_J[4:end]])
	if ((opt_I = GMT.parse_I(d, "", [:I :inc :increment :spacing], 'I')) != "")	# Need the 'I' to not fall into parse_I() exceptions
		t = split(opt_I[4:end], '/')
		(length(t) == 1) ? append!(opts, ["-tr", t[1], t[1]]) : append!(opts, ["-tr", t[1], t[2]])
	end
	return d, opts, (opt_R != "" || length(opt_J) > 1 || opt_I != "")
end

# ---------------------------------------------------------------------------------------------------
function gdal_opts2vec(opts)::Vector{String}
	# Break up a string of options into a vector string as it's needed by GDAL lower level functions
	(opts == "") && return String[]
	(isempty(opts) || (isa(opts, Vector{<:AbstractString}) && length(opts) > 1)) && return opts	# if already a vec
	(eltype(opts) != Char && !(eltype(opts) <: AbstractString)) &&
		error("Options for GDAL must be a string or a vector of one string")
	_opts = (isa(opts, Vector{<:AbstractString})) ? opts[1] : opts

	ind = helper_opts2vec(_opts)
	isempty(ind) && return split(_opts)			# Perfect, just a 'simple' options list, split and go
	_opts = _opts[1:ind[1][1]-1] * replace(_opts[ind[1][1]+1 : ind[2][1]-1], ' ' => '\x7f') * _opts[ind[2][1]+1 : end]

	while (!isempty(helper_opts2vec(_opts)))	# See if we have more
		ind = helper_opts2vec(_opts)
		_opts = _opts[1:ind[1][1]-1] * replace(_opts[ind[1][1]+1 : ind[2][1]-1], ' ' => '\x7f') * _opts[ind[2][1]+1 : end]
	end

	v = split(_opts)
	[v[i] = replace(v[i], '\x7f' => ' ') for i = 1:length(v)]		# Undo the Char(127) char
end

function helper_opts2vec(opts::String)
	ind = findall("\"", opts)
	isempty(ind) && (ind = findall("'", opts))
	(!isempty(ind) && rem(length(ind), 2) != 0) && error("Delimiting characters like $(opts[ind[1][1]]) must come in pairs")
	return ind
end

# ---------------------------------------------------------------------------------------------------
function default_gdopts!(ds, opts::Vector{String}, dest::String)
	# Assign some default options in function of the driver and data type
	startswith(dest, "/vsimem") && return nothing	# In this case we won't set any defaults

	driver = shortname(getdriver(ds))
	dt = GDALGetRasterDataType(ds.ptr)
	# For some reason when MEM driver (only it?) dt comes == 1, even when data is float. So check again.
	(startswith(lowercase(driver), "mem") && dt == 1 && isa(ds, Gdal.AbstractRasterBand)) && (dt = GDALGetRasterDataType(getband(ds,1).ptr))

	ext = lowercase(splitext(dest)[2])
	isTiff = (ext == ".tif" || ext == ".tiff")
	isNC   = (driver == "netCDF" || ext == ".nc"  || ext == ".grd") && (width(ds) > 128 && height(ds > 128))
	(ext == ".grd") && append!(opts, ["-of", "netCDF"])		# Accept .grd as meaning netcdf and not Surfer ascii (GDAL default)
	((dt == 1 || isTiff) && !any(startswith.(opts, "COMPRESS"))) && append!(opts, ["-co", "COMPRESS=DEFLATE", "-co", "PREDICTOR=2"])
	((dt == 1 || isTiff) && !any(startswith.(opts, "TILED"))) && append!(opts, ["-co", "TILED=YES"])
	(dt >= 6 && !any(startswith.(opts, "a_nodata"))) && append!(opts, ["-a_nodata","NaN"])
	(isNC) && append!(opts,["-co", "FORMAT=NC4", "-co", "COMPRESS=DEFLATE", "-co", "ZLEVEL=4"])
end

# ---------------------------------------------------------------------------------------------------
function get_gdaldataset(data, opts)
	# Get a GDAL dataset from either a file name, a GMT grid or image, or a dataset itself
	# In case of a file name we must be careful and deal with possible "+b" band requests from GMT.
	needclose = false
	if isa(data, AbstractString)			# Check also for remote files (those that start with a @). MAY SCREW VIOLENTLY
		name, ext = splitext(data)			# Be carefull, the name may carry a bands request. e.g. "LC08__cube.tiff+b3,2,1"
		name = ((ind = findfirst("+", ext)) === nothing) ? data : name * ext[1:ind[1]-1]
		if (ind !== nothing && ext[ind[1]+1] == 'b')	# So we must convert the "+b3,2,1" into GDAL syntax
			o = split(ext[ind[1]+2:end], ",")			# But we must add 1 because GDAL band count is 1-based and GMT 0-based
			p = parse.(Int, o) .+ 1
			o = [string(c) for c in p]
			if (isa(opts, Vector{String}))
				[append!(opts, ["-b", o[k]]) for k = 1:length(o)]
			else
				opts *= " -b" * join(o, " -b ")
			end
		end
		ds = (name[1] == '@') ? Gdal.unsafe_read(gmtwhich(name)[1].text[1]) : Gdal.unsafe_read(name)
		needclose = true					# For some reason file remains open and we must close it explicitly
	elseif (isa(data, GMTgrid) || isa(data, GMTimage) || GMT.isGMTdataset(data) || isa(data, Matrix{<:Real}))
		ds = gmt2gd(data)
	elseif (isa(data, AbstractGeometry))	# VERY UNSATISFACTORY. SHOULD BE ABLE TO GETPARENT (POSSIBLE?)
		ds = wrapgeom(data)
	else
		ds = data							# If it's not a GDAL dataset or convenient descendent, shit will follow soon
	end
	(ds === nothing) && error("Error fetching the GDAL dataset from input $(typeof(data))")
	return ds, needclose
end

# ---------------------------------------------------------------------------------------------------
"""
    dither(indata; n_colors=256, save="", gdataset=false)

Convert a 24bit RGB image to 8bit paletted.
- Use the `save=fname` option to save the result to disk in a GeoTiff file `fname`. If `fname`
  has no extension a `.tif` one will be appended. Alternatively give file names with extension
  `.png` or `.nc` to save the file in one of those formats.
- `gdataset=true`: to return a GDAL dataset. The default is to return a GMTimage object.
- `n_colors`: Select the number of colors in the generated color table. Defaults to 256.
"""
function dither(indata, opts=String[]; n_colors::Integer=256, save::String="", gdataset::Bool=false)
	# ...
	src_ds, needclose = get_gdaldataset(indata, "")
	(nraster(src_ds) < 3) && error("Input image must have at least 3 bands")
	(isa(indata, GMTimage) && !startswith(indata.layout, "TRB")) &&
		error("Image memory layout must be `TRB` and not $(indata.layout). Load image with gdaltranslate()")
	r_band, g_band, b_band = getband(src_ds, 1), getband(src_ds, 2), getband(src_ds, 3)

	drv_name = "MEM"
	if (save != "")
		fn, ext = splitext(save)
		ext = lowercase(ext)
		if     (ext == "" || ext == ".tif" || ext == ".tiff")  drv_name = "GTiff"
		elseif (ext == ".png") drv_name = "PNG"
		elseif (ext == ".nc")  drv_name = "netCDF"
		else
			@warn("Format not supported. Only TIF, PNG or netCDF are allowed. Resorting to TIF")
			drv_name = "GTiff"
			save = fn * ".tif"
		end
		(ext == "") && (save *= ".tif")
	end

	(drv_name == "GTiff") && append!(opts, ["TILED=YES", "TILED=YES", "COMPRESS=DEFLATE", "PREDICTOR=2"])
	dst_ds = create(save, driver=getdriver(drv_name), width=width(src_ds), height=height(src_ds),
	                nbands=1, dtype=UInt8, options=opts)
	try
		setproj!(dst_ds, getproj(src_ds))
		setgeotransform!(dst_ds, getgeotransform(src_ds))
	catch
	end
	dst_band = getband(dst_ds, 1)
	ct = createcolortable(UInt32(1))

	ComputeMedianCutPCT(r_band, g_band, b_band, n_colors, ct)
	setcolortable!(dst_band, ct)
	DitherRGB2PCT(r_band, g_band, b_band, dst_band, ct)
	(needclose) && GDALClose(src_ds.ptr)
	if (save != "")						# Because a file was writen
		destroy(dst_ds)
		return nothing
	end
	(gdataset) ? dst_ds : gd2gmt(dst_ds)
end


text:
context("Mesh Types") do

context("Merging Mesh") do
    baselen = 0.4f0
    dirlen = 2f0
    a = [
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(dirlen, baselen, baselen)), RGBA(1f0,0f0,0f0,1f0)),
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, dirlen, baselen)), RGBA(0f0,1f0,0f0,1f0)),
        (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, baselen, dirlen)), RGBA(0f0,0f0,1f0,1f0))
    ]

    am = map(GLNormalMesh, a)
    axis = merge(am)
    axis2 = merge(am)
    @fact axis --> axis2 # test ==
    @fact typeof(am[1]) --> GLNormalColorMesh

    @fact typeof(axis) --> GLNormalAttributeMesh
    @fact typeof(GLPlainMesh(axis)) --> GLPlainMesh

    @fact vertextype(axis) --> Point{3, Float32}
    @fact normaltype(axis) --> Normal{3, Float32}
    @fact facetype(axis) --> Face{3, Cuint, -1}
    @fact hasvertices(axis) --> true
    @fact hasfaces(axis) --> true
    @fact hasnormals(axis) --> true
    @fact hascolors(axis) --> false
    end

context("Show") do
    baselen = 0.4f0
    dirlen = 2f0
    m = GLNormalMesh(HyperRectangle{3,Float32}(Vec3f0(baselen),
                     Vec3f0(dirlen, baselen, baselen)))
    io = IOBuffer()
    show(io,m)
    seekstart(io)
    s = "HomogenousMesh(\n    faces: 12xGeometryTypes.Face{3,UInt32,-1},     normals: 24xGeometryTypes.Normal{3,Float32},     vertices: 24xFixedSizeArrays.Point{3,Float32}, )\n"
    #@fact readall(io) --> s #Win32 and Win64 have different ordering it seems.
end

context("Primitives") do
    # issue #16
    #m = HomogenousMesh{Point{3,Float64},Face{3,Int,0}}(Sphere(Point(0,0,0), 1))
    #@fact length(vertices(m)) --> 145
    #@fact length(faces(m)) --> 288
end

context("Slice") do
    test_mesh = HomogenousMesh(Point{3,Float64}[
     Point{3,Float64}(0.0,0.0,10.0),
     Point{3,Float64}(0.0,10.0,10.0),
     Point{3,Float64}(0.0,0.0,0.0),
     Point{3,Float64}(0.0,10.0,0.0),
     Point{3,Float64}(10.0,0.0,10.0),
     Point{3,Float64}(10.0,10.0,10.0),
     Point{3,Float64}(10.0,0.0,0.0),
     Point{3,Float64}(10.0,10.0,0.0),
    ],Face{3,Int,0}[
     Face{3,Int,0}(3,7,5)
     Face{3,Int,0}(1,3,5)
     Face{3,Int,0}(1,2,3)
     Face{3,Int,0}(2,4,3)
     Face{3,Int,0}(1,5,6)
     Face{3,Int,0}(2,1,6)
     Face{3,Int,0}(4,8,3)
     Face{3,Int,0}(3,8,7)
     Face{3,Int,0}(7,8,6)
     Face{3,Int,0}(7,6,5)
     Face{3,Int,0}(2,6,4)
     Face{3,Int,0}(4,6,8)]
    )
    s1 = slice(test_mesh, 1.0)
    s2 = slice(test_mesh, 2.0)
    s3 = slice(test_mesh, 0.0)

    @fact length(s1) --> 8
    @fact length(s1) --> length(s2)
    @fact length(s3) --> 0
    exp1 = [([0.0,1.0],[0.0,0.0]),([0.0,0.0],[1.0,0.0]),
                      ([1.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,1.0]),
                      ([10.0,1.0],[10.0,10.0]),([10.0,10.0],[1.0,10.0]),
                      ([1.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,1.0])]
    exp2 = [([0.0,2.0],[0.0,0.0]),([0.0,0.0],[2.0,0.0]),
                      ([2.0,0.0],[10.0,0.0]),([10.0,0.0],[10.0,2.0]),
                      ([10.0,2.0],[10.0,10.0]),([10.0,10.0],[2.0,10.0]),
                      ([2.0,10.0],[0.0,10.0]),([0.0,10.0],[0.0,2.0])]
    @fact length(s1) --> length(exp1)
    @fact length(s2) --> length(exp2)
end

context("checkbounds") do
    m1 = HomogenousMesh([Point{3,Float64}(0.0,0.0,10.0),
                         Point{3,Float64}(0.0,10.0,10.0),
                         Point{3,Float64}(0.0,0.0,0.0)],
                        [Face{3,Int,0}(1,2,3)])
    @fact checkbounds(m1) --> true
    m2 = HomogenousMesh([Point{3,Float64}(0.0,0.0,10.0),
                         Point{3,Float64}(0.0,10.0,10.0),
                         Point{3,Float64}(0.0,0.0,0.0)],
                        [Face{3,Int,-1}(1,2,3)])
    @fact checkbounds(m2) --> false
    # empty case
    m3 = HomogenousMesh([Point{3,Float64}(0.0,0.0,10.0),
                         Point{3,Float64}(0.0,10.0,10.0),
                         Point{3,Float64}(0.0,0.0,0.0)],
                        Face{3,Int,-1}[])
    @fact checkbounds(m3) --> true
end

context("vertex normals") do
    test_mesh = HomogenousMesh(Point{3,Float64}[Point{3,Float64}(0.0,0.0,10.0),
     Point{3,Float64}(0.0,10.0,10.0),
     Point{3,Float64}(0.0,0.0,0.0),
     Point{3,Float64}(0.0,10.0,0.0),
     Point{3,Float64}(10.0,0.0,10.0),
     Point{3,Float64}(10.0,10.0,10.0),
     Point{3,Float64}(10.0,0.0,0.0),
     Point{3,Float64}(10.0,10.0,0.0),
    ],Face{3,Int,0}[
     Face{3,Int,0}(3,7,5)
     Face{3,Int,0}(1,3,5)
     Face{3,Int,0}(1,2,3)
     Face{3,Int,0}(3,2,4)
     Face{3,Int,0}(1,5,6)
     Face{3,Int,0}(2,1,6)
     Face{3,Int,0}(4,8,3)
     Face{3,Int,0}(3,8,7)
     Face{3,Int,0}(7,8,6)
     Face{3,Int,0}(5,7,6)
     Face{3,Int,0}(2,6,4)
     Face{3,Int,0}(4,6,8)]
    )
    ns = normals(test_mesh.vertices, test_mesh.faces)
    @fact length(ns) --> length(test_mesh.vertices)
    expect = [Normal(-0.408248290463863,-0.408248290463863,0.816496580927726),
              Normal(-0.816496580927726,0.408248290463863,0.408248290463863),
              Normal(-0.5773502691896257,-0.5773502691896257,-0.5773502691896257),
              Normal(-0.408248290463863,0.816496580927726,-0.408248290463863),
              Normal(0.408248290463863,-0.816496580927726,0.408248290463863),
              Normal(0.5773502691896257,0.5773502691896257,0.5773502691896257),
              Normal(0.816496580927726,-0.408248290463863,-0.408248290463863),
              Normal(0.408248290463863,0.408248290463863,-0.816496580927726)]
    @fact ns --> expect
end

end

ext: jl


text:
"""
    AbstractArchSpace

Abstract callable struct representing an architecture search space.

Architecture spaces define a range of possible hyperparameters for a model architecture. Used to create new models or parts of models.

Return an `AbstractVertex` from the search space when invoked with an input `AbstractVertex`.
"""
abstract type AbstractArchSpace end

"""
    AbstractLayerSpace

Abstract callable struct representing a search space of layers.

Return a layer in the search space given an input size and (optionally) a random number generator.
"""
abstract type AbstractLayerSpace end

"""
    AbstractParSpace{N, T}

Abstract callable struct representing a seach space of `N`D parameters of type `T`.

Return a tuple of parameters from the search space when invoked. Random number generator may be supplied.
"""
abstract type AbstractParSpace{N, T} end

"""
    AbstractWeightInit

Abstract type for weight initialization strategies
"""
abstract type AbstractWeightInit end

"""
    DefaultWeightInit <: AbstractWeightInit

Use the layer default weight initialization.
"""
struct DefaultWeightInit <: AbstractWeightInit end

"""
    IdentityWeightInit <: AbstractWeightInit
    IdentityWeightInit()

Initialize weights with an identity mapping.
"""
struct IdentityWeightInit <: AbstractWeightInit end

"""
    PartialIdentityWeightInit <: AbstractWeightInit
    PartialIdentityWeightInit(inoffset, outoffset)

Initialize weights with an identity mapping.

Parameter `inoffset` and `outoffset` shifts the weigths in case `nin != nout`. This is so that a concatenation can also be an identity mapping as a whole.
"""
struct PartialIdentityWeightInit <: AbstractWeightInit
    inoffset::Int
    outoffset::Int
end

"""
    ZeroWeightInit <: AbstractWeightInit
    ZeroWeightInit()

Initialize weights as zeros.

Main use case is to let residual layers be and identity mapping
"""
struct ZeroWeightInit <: AbstractWeightInit end

"""
    BaseLayerSpace
    BaseLayerSpace(outsizes, activationfunctions)

Search space for basic attributes common to all layers.

`outsizes` is the output sizes (number of neurons).
`activationfunctions` is the activation functions.
"""
struct BaseLayerSpace{T}
    nouts::AbstractParSpace{1, <:Integer}
    acts::AbstractParSpace{1, T}
end
BaseLayerSpace(outsizes, acts) = BaseLayerSpace(parspaceof(outsizes), parspaceof(acts))
outsize(s::BaseLayerSpace, rng=rng_default) = s.nouts(rng)
activation(s::BaseLayerSpace, rng=rng_default) = s.acts(rng)

"""
    SingletonParSpace{N, T} <:AbstractParSpace{N, T}
    SingletonParSpace(p::T...)
    Singleton2DParSpace(p::T)

Singleton search space. Has exactly one value per dimension.

`Singleton2DParSpace` is a convenience constructor for a 2D `SingletonParSpace` of `[p, p]`.
"""
struct SingletonParSpace{N, T} <:AbstractParSpace{N, T}
    p::NTuple{N, T}
end
SingletonParSpace(p::T...) where T = SingletonParSpace(p)
Singleton2DParSpace(p::T) where T = SingletonParSpace(p,p)
(s::SingletonParSpace)(rng=nothing) = s.p
(s::SingletonParSpace{1, T})(rng=nothing) where T = s.p[1]

"""
    ParSpace{N, T} <:AbstractParSpace{N, T}
    ParSpace(p::AbstractVector{T}...)
    ParSpace1D(p...)
    ParSpace2D(p::AbstractVector)

Search space for parameters.

Return independent uniform random values for all `N` dimensions from the search space when invoked.

`ParSpace1D` is a convenience constructor for an 1D `ParSpace` of `p`.
`ParSpace2D` is a convenience constructor for a 2D `ParSpace` of `[p, p]``.
"""
struct ParSpace{N, T} <:AbstractParSpace{N, T}
    p::NTuple{N, AbstractVector{T}}
end
ParSpace(p::AbstractVector{T}...) where T = ParSpace(p)
ParSpace1D(p...) = ParSpace(collect(p))
ParSpace2D(p::AbstractVector) = ParSpace(p,p)
(s::ParSpace)(rng=rng_default) = rand.((rng,), s.p)
(s::ParSpace{1})(rng=rng_default) = rand(rng, s.p[1])

"""
    CoupledParSpace{N, T} <:AbstractParSpace{N, T}
    CoupledParSpace(p::AbstractParSpace{1, T}, N)
    CoupledParSpace(p::AbstractVector{T}, N)

Search space for parameters.

Return the same uniformly sampled value for all `N` dimensions from the search space when invoked.
"""
struct CoupledParSpace{N, T} <:AbstractParSpace{N, T}
    p::AbstractParSpace{1, T}
end
CoupledParSpace(p::AbstractParSpace{1, T}, N) where T = CoupledParSpace{N}(p)
CoupledParSpace(p::AbstractVector{T}, N) where T = CoupledParSpace{N, T}(ParSpace(p))
function(s::CoupledParSpace{N})(rng=rng_default) where N
    val = s.p(rng)
    return ntuple(i -> val, N)
end

parspaceof(x) = parspaceof(1, x)
parspaceof(N, x) = parspaceof(Val(N), x)
parspaceof(::Val{N}, s::AbstractParSpace{N}) where N = s
parspaceof(::Val{N}, s::AbstractParSpace{1}) where N = parspaceof(N, s.p[1])
parspaceof(::Val{N}, s::CoupledParSpace{1}) where N = CoupledParSpace(s.p, N)

parspaceof(n::Val{N}, x) where N = parspaceof(n, ntuple(i -> x, N))
parspaceof(::Val{N}, x::NTuple{N}) where N = SingletonParSpace(x...)
parspaceof(::Val{N}, x::NTuple{N, AbstractVector}) where N = ParSpace(x...)

"""
    NamedLayerSpace <:AbstractLayerSpace
    NamedLayerSpace(name::String, s::AbstractLayerSpace)

Adds a `name` to an `AbstractLayerSpace`.
"""
struct NamedLayerSpace <:AbstractLayerSpace
    name::String
    s::AbstractLayerSpace
end
NaiveNASlib.name(::AbstractLayerSpace) = ""
NaiveNASlib.name(s::NamedLayerSpace) = s.name
function (s::NamedLayerSpace)(in::Integer,rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    ismissing(outsize) && return s.s(in, rng)
    return s.s(in, rng, outsize=outsize, wi=wi)
end

"""
    LoggingLayerSpace <: AbstractLayerSpace
    LoggingLayerSpace(s::AbstractLayerSpace; level=Logging.Debug, nextlogfun=() -> PrefixLogger("   "))
    LoggingLayerSpace(msgfun, s::AbstractLayerSpace; level = Logging.Debug, nextlogfun = () -> PrefixLogger("   "))
    LoggingLayerSpace(level::LogLevel, msgfun, nextlogfun, s::AbstractLayerSpace)

Logs `msgfun(layer)` at loglevel `level` after creating a `layer` from `s`.

Calling `nextlogfun()` produces an `AbstractLogger` which will be used when creating `layer` from `s`.

By default, this is used to add a level of indentation to subsequent logging calls which makes logs of hierarchical archspaces easier to read. Set `nextlogfun = () -> current_logger()` to remove this behaviour.
"""
struct LoggingLayerSpace{F,L<:LogLevel,LF,T <: AbstractLayerSpace}  <: AbstractLayerSpace
    msgfun::F
    level::L
    nextlogfun::LF
    s::T
end
LoggingLayerSpace(s::AbstractLayerSpace; level=Logging.Debug, nextlogfun=() -> PrefixLogger("   ")) = LoggingLayerSpace(l -> "Create $l from $(name(s))", s, level=level, nextlogfun=nextlogfun)
LoggingLayerSpace(msgfun, s::AbstractLayerSpace; level = Logging.Debug, nextlogfun = () -> PrefixLogger("   ")) = LoggingLayerSpace(msgfun, level, nextlogfun, s)
NaiveNASlib.name(s::LoggingLayerSpace) = name(s.s)
function (s::LoggingLayerSpace)(in::Integer,rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    layer = with_logger(s.nextlogfun()) do
        ismissing(outsize) ? s.s(in, rng) : s.s(in, rng, outsize=outsize, wi=wi)
    end
    msg = s.msgfun(layer)
    @logmsg s.level msg
    return layer
end

"""
    DenseSpace <:AbstractLayerSpace
    DenseSpace(base::BaseLayerSpace)
    DenseSpace(outsizes, activations)

Search space of Dense layers.
"""
struct DenseSpace <:AbstractLayerSpace
    base::BaseLayerSpace
end
DenseSpace(outsizes, activations) = DenseSpace(BaseLayerSpace(outsizes, activations))
(s::DenseSpace)(in::Integer,rng=rng_default; outsize=outsize(s.base,rng), wi=DefaultWeightInit(), densefun=Dense) = densefun(in, outsize, activation(s.base,rng); denseinitW(wi)...)

denseinitW(::DefaultWeightInit) = ()
denseinitW(::IdentityWeightInit) = (init = Flux.identity_init,)
denseinitW(wi::PartialIdentityWeightInit) = (init = (args...) -> Flux.identity_init(args...;shift=(wi.outoffset, wi.inoffset)),)
denseinitW(::ZeroWeightInit) = (init = zeros,)

"""
    ConvSpace{N} <:AbstractLayerSpace

    ConvSpace{N}(;outsizes, kernelsizes, activations=identity, strides=1, dilations=1, paddings=SamePad(), convfuns=Conv)
    ConvSpace(convfun::AbstractParSpace, base::BaseLayerSpace, ks::AbstractParSpace, stride::AbstractParSpace, dilation::AbstractParSpace, pad)

Search space of `N`D convolutional layers.

Constructor with keyword arguments takes scalars, vectors or `AbstractParSpace`s as inputs.
"""
struct ConvSpace{N} <:AbstractLayerSpace
    convfun::AbstractParSpace
    base::BaseLayerSpace
    kernelsize::AbstractParSpace{N, <:Integer}
    stride::AbstractParSpace
    dilation::AbstractParSpace
    pad::AbstractParSpace
end

ConvSpace{N}(;outsizes, kernelsizes, activations=identity, strides=1, dilations=1, paddings=SamePad(), convfuns=Conv) where N =
ConvSpace(
    parspaceof(convfuns),
    BaseLayerSpace(parspaceof(outsizes), parspaceof(activations)),
    parspaceof(N, kernelsizes),
    parspaceof(N, strides),
    parspaceof(N, dilations),
    parspaceof(paddings))

function (s::ConvSpace)(insize::Integer, rng=rng_default; outsize = outsize(s.base, rng), wi=DefaultWeightInit())
    convfun = s.convfun(rng)
    ks = Tuple(s.kernelsize(rng))
    act = activation(s.base, rng)
    pad = s.pad(rng)
    stride = s.stride(rng)
    dilation = s.dilation(rng)
    return convfun(ks, insize=>outsize, act; pad=pad, stride=stride, dilation=dilation, convinitW(wi)...)
end

convinitW(::DefaultWeightInit) = ()
convinitW(::IdentityWeightInit) = (init=Flux.identity_init,)
convinitW(wi::PartialIdentityWeightInit) = (init = (args...) -> Flux.identity_init(args...;shift=(0,0,wi.inoffset, wi.outoffset)),)
convinitW(::ZeroWeightInit) = (init=(args...) -> zeros(Float32,args...),)

"""
    BatchNormSpace <:AbstractLayerSpace
    BatchNormSpace(activationfunctions...)
    BatchNormSpace(activationfunctions::AbstractVector)

Search space of BatchNorm layers.
"""
struct BatchNormSpace <:AbstractLayerSpace
    acts::AbstractParSpace
end
BatchNormSpace(act::Function) = BatchNormSpace(SingletonParSpace(act))
BatchNormSpace(act, acts...) = BatchNormSpace(ParSpace1D(act,acts...))
BatchNormSpace(acts::AbstractVector) = BatchNormSpace(acts...)
(s::BatchNormSpace)(in::Integer, rng=rng_default;outsize=nothing, wi=nothing) = BatchNorm(in, s.acts(rng))

"""
    PoolSpace{N} <:AbstractLayerSpace
    PoolSpace{N}(;windowsizes, strides=1, paddings=SamePad(), poolfun=[MaxPool, MeanPool])

Search space of `N`D pooling layers.

Constructor with keyword arguments takes scalars/tuples, vectors or `AbstractParSpace`s as inputs.
"""
struct PoolSpace{N} <:AbstractLayerSpace
    poolfun::AbstractParSpace
    ws::AbstractParSpace{N, <:Integer}
    stride::AbstractParSpace
    pad::AbstractParSpace
end
PoolSpace{N}(;windowsizes, strides=1, paddings=0, poolfuns=[MaxPool, MeanPool]) where N = PoolSpace{N}(
    parspaceof(poolfuns),
    parspaceof(N, windowsizes),
    parspaceof(N, strides),
    parspaceof(paddings)
)

function (s::PoolSpace)(in::Integer, rng=rng_default;outsize=nothing, wi=nothing)
    poolfun = s.poolfun(rng)
    ws = Tuple(s.ws(rng))
    stride = s.stride(rng)
    pad = s.pad(rng)
    poolfun(ws, stride=stride, pad=pad)
end

default_logging() = logged(level=Logging.Debug)
"""
    LayerVertexConf
    LayerVertexConf(layerfun, traitfun)

Generic configuration template for computation graph vertices with `Flux` layers as their computation.

Both `layerfun` and `traitfun` are forwarded to `NaiveNASflux.fluxvertex`.

Intention is to make it easy to add logging, validation and pruning metrics in an uniform way.
"""
struct LayerVertexConf{F, T}
    layerfun::F
    traitfun::T
end
LayerVertexConf() = LayerVertexConf(ActivationContribution ∘ LazyMutable, validated() ∘ default_logging())

"""
    Shielded(base=LayerVertexConf(); allowed = tuple())

Create a [`LayerVertexConf`](@ref) which is shielded from mutation.

Keyword `allowed` can be used to supply a tuple (or array) of `AbstractMutation` types to allow.
"""
Shielded(base=LayerVertexConf(); allowed = tuple()) = let Shield(t) = MutationShield(t, allowed...)
    LayerVertexConf(base.layerfun, Shield ∘ base.traitfun)
end

(c::LayerVertexConf)(in::AbstractVertex, l) = fluxvertex(l,in,layerfun=c.layerfun, traitfun=c.traitfun)
(c::LayerVertexConf)(name::String, in::AbstractVertex, l) = fluxvertex(name, l,in,layerfun=c.layerfun, traitfun=c.traitfun)

Base.Broadcast.broadcastable(c::LayerVertexConf) = Ref(c)

"""
    ConcConf
    ConcConf(layerfun, traitfun)

Generic configuration template for concatenation of vertex outputs.

Both `layerfun` and `traitfun` are forwarded to `NaiveNASflux.concat`.
"""
struct ConcConf{F, T}
    layerfun::F
    traitfun::T
end
ConcConf() = ConcConf(ActivationContribution, validated() ∘ default_logging())

(c::ConcConf)(in::AbstractVector{<:AbstractVertex}) = c(in...)
(c::ConcConf)(in::AbstractVertex) = in
(c::ConcConf)(ins::AbstractVertex...) = concat(ins..., traitfun = c.traitfun, layerfun=c.layerfun)
(c::ConcConf)(name::String, in::AbstractVector{<:AbstractVertex}) = c(name, in...)
(c::ConcConf)(name::String, in::AbstractVertex) = in
(c::ConcConf)(name::String, ins::AbstractVertex...) = concat(ins..., traitfun = c.traitfun ∘ named(name), layerfun=c.layerfun)

"""
    LoggingArchSpace <: AbstractArchSpace
    LoggingArchSpace(s::AbstractArchSpace; level=Logging.Debug, nextlogfun=in -> PrefixLogger("   "))
    LoggingArchSpace(msgfun, s::AbstractArchSpace; level=Logging.Debug, nextlogfun=in -> PrefixLogger("   "))
    LoggingArchSpace(msgfun::Function, level::LogLevel, nextlogfun, s::AbstractArchSpace)

Logs `msgfun(vertex)` at loglevel `level` after creating a `vertex` from `s`.

Calling `nextlogfun(in)` where `in` is the input vertex produces an `AbstractLogger` which will be used when creating `vertex` from `s`.

By default, this is used to add a level of indentation to subsequent logging calls which makes logs of hierarchical archspaces easier to read. Set `nextlogfun = e -> current_logger()` to remove this behaviour.
"""
struct LoggingArchSpace{F,L<:LogLevel,LF,T <: AbstractArchSpace} <: AbstractArchSpace
    msgfun::F
    level::L
    nextlogfun::LF
    s::T
end
LoggingArchSpace(s::AbstractArchSpace; level=Logging.Debug, nextlogfun=in -> PrefixLogger("   ")) = LoggingArchSpace(v -> "Created $(name(v))", level, nextlogfun, s)
LoggingArchSpace(msgfun, s::AbstractArchSpace; level=Logging.Debug, nextlogfun=in -> PrefixLogger("   ")) = LoggingArchSpace(msgfun, level, nextlogfun, s)
function (s::LoggingArchSpace)(in::AbstractVertex,rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    layer = with_logger(s.nextlogfun(in)) do
        ismissing(outsize) ? s.s(in, rng) : s.s(in, rng, outsize=outsize, wi=wi)
    end
    @logmsg s.level s.msgfun(layer)
    return layer
end
function (s::LoggingArchSpace)(namestr::String, in::AbstractVertex,rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    layer = with_logger(s.nextlogfun(in)) do
        ismissing(outsize) ? s.s(namestr, in, rng) : s.s(namestr, in, rng, outsize=outsize, wi=wi)
    end
    @logmsg s.level s.msgfun(layer)
    return layer
end

"""
    VertexSpace <:AbstractArchSpace
    VertexSpace([conf::LayerVertexConf], lspace::AbstractLayerSpace)

Search space of one `AbstractVertex` with compuation drawn from `lspace`.

`conf` is used to attach other metadata to the vertex. See [`LayerVertexConf`](@ref).
"""
struct VertexSpace <:AbstractArchSpace
    conf::LayerVertexConf
    lspace::AbstractLayerSpace
end
VertexSpace(lspace::AbstractLayerSpace) = VertexSpace(LayerVertexConf(), lspace)

(s::VertexSpace)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = s.conf(in,  create_layer(outsize, nout(in), s.lspace, wi, rng))
(s::VertexSpace)(namestr::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = s.conf(join(filter(!isempty, [namestr, name(s.lspace)]), "."), in, create_layer(outsize, nout(in), s.lspace, wi, rng))
create_layer(::Missing, insize::Integer, ls::AbstractLayerSpace, wi, rng) = ls(insize, rng, wi=wi)
create_layer(outsize::Integer, insize::Integer, ls::AbstractLayerSpace, wi, rng) = ls(insize, rng, outsize=outsize, wi=wi)

"""
    NoOpArchSpace <: AbstractArchSpace

Returns input vertex without any modification.
"""
struct NoOpArchSpace <: AbstractArchSpace end
(::NoOpArchSpace)(v::AbstractVertex, args...; kwargs...) = v
(::NoOpArchSpace)(name::String, v::AbstractVertex, args...; kwargs...) = v

"""
    ArchSpace <:AbstractArchSpace
    ArchSpace(ss::AbstractLayerSpace...; conf=LayerVertexConf())
    ArchSpace(ss::AbstractArchSpace...)

Search space of `AbstractArchSpace`s.

Draws one vertex from one of `ss` (uniformly selected) when invoked.
"""
struct ArchSpace <:AbstractArchSpace
    s::AbstractParSpace{1, <:AbstractArchSpace}
end
ArchSpace(l::AbstractLayerSpace; conf=LayerVertexConf()) = ArchSpace(VertexSpace(conf,l))
ArchSpace(l::AbstractLayerSpace, ls::AbstractLayerSpace...; conf=LayerVertexConf()) = ArchSpace(VertexSpace.(conf,[l,ls...])...)
ArchSpace(s::AbstractArchSpace) = ArchSpace(SingletonParSpace(s))
ArchSpace(s::AbstractArchSpace, ss::AbstractArchSpace...) = ArchSpace(ParSpace([s, ss...]))

(s::ArchSpace)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = s.s(rng)(in, rng, outsize=outsize, wi=wi)
(s::ArchSpace)(name::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = s.s(rng)(name, in, rng, outsize=outsize, wi=wi)

"""
    ConditionalArchSpace{P, S1, S2} <: AbstractArchSpace
    ConditionalArchSpace(predicate, iftrue, iffalse=NoOpArchSpace())
    ConditionalArchSpace(;predicate, iftrue, iffalse=NoOpArchSpace())

Use `iftrue` if `predicate(invertex)` returns true, `iffalse` otherwise.
"""
Base.@kwdef struct ConditionalArchSpace{P, S1, S2} <: AbstractArchSpace
    predicate::P
    iftrue::S1
    iffalse::S2 = NoOpArchSpace()
    ConditionalArchSpace(predicate::P, iftrue::S1, iffalse::S2=NoOpArchSpace()) where {P, S1, S2} =new{P, S1, S2}(predicate, iftrue, iffalse)
end

(s::ConditionalArchSpace)(in::AbstractVertex, args...; kwargs...) = s.predicate(in) ? s.iftrue(in, args...; kwargs...) : s.iffalse(in, args...; kwargs...)
(s::ConditionalArchSpace)(name::String, in::AbstractVertex, args...; kwargs...) = s.predicate(in) ? s.iftrue(name, in, args...; kwargs...) : s.iffalse(name, in, args...; kwargs...)

"""
    RepeatArchSpace <:AbstractArchSpace
    RepeatArchSpace(s::AbstractArchSpace, r::Integer)
    RepeatArchSpace(s::AbstractArchSpace, r::AbstractVector{<:Integer})

Search space of repetitions of another `AbstractArchSpace` where number of repetitions is uniformly drawn from `r`.

Output of each generated candidate is input to next and the last output is returned.
"""
struct RepeatArchSpace <:AbstractArchSpace
    s::AbstractArchSpace
    r::AbstractParSpace{1, <:Integer}
end
RepeatArchSpace(s::AbstractArchSpace, r::Integer) = RepeatArchSpace(s, SingletonParSpace(r))
RepeatArchSpace(s::AbstractArchSpace, r::AbstractVector{<:Integer}) = RepeatArchSpace(s, ParSpace(r))

(s::RepeatArchSpace)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = foldl((next, i) -> s.s(next, rng, outsize=outsize, wi=repeatinitW(wi, next, outsize)), 1:s.r(rng), init=in)
(s::RepeatArchSpace)(name::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = foldl((next, i) -> s.s(join([name,".", i]), next, rng, outsize=outsize, wi=repeatinitW(wi, next, outsize)), 1:s.r(rng), init=in)

repeatinitW(wi::AbstractWeightInit, invertex, outsize) = wi
repeatinitW(wi::PartialIdentityWeightInit, invertex, outsize) = nout(invertex) == outsize ? IdentityWeightInit() : wi
repeatinitW(wi::PartialIdentityWeightInit, invertex, ::Missing) = wi

"""
    ArchSpaceChain <:AbstractArchSpace
    ArchSpaceChain(s::AbstractArchSpace...)

Chains multiple `AbstractArchSpace`s after each other.

Input vertex will be used to generate an output vertex from the first `AbstractArchSpace` in the chain which is then used to generate a next output vertex from the next `AbstractArchSpace` in the chain and so on. The output from the last `AbstractArchSpace` is returned.
"""
struct ArchSpaceChain{S<:AbstractVector{<:AbstractArchSpace}} <:AbstractArchSpace
    s::S
end
ArchSpaceChain(s::AbstractArchSpace...) = ArchSpaceChain(collect(s))
(s::ArchSpaceChain)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = foldl((next, ss) -> ss(next, rng, outsize=outsize, wi=repeatinitW(wi, next, outsize)), s.s, init=in)
(s::ArchSpaceChain)(name::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit()) = foldl((next, i) -> s.s[i](join([name,".", i]), next, rng, outsize=outsize, wi=repeatinitW(wi, next, outsize)), eachindex(s.s), init=in)

"""
    ForkArchSpace <:AbstractArchSpace
    ForkArchSpace(s::AbstractArchSpace, r::Integer; conf=ConcConf())
    ForkArchSpace(s::AbstractArchSpace, r::AbstractVector{<:Integer}; conf=ConcConf())

Search space of parallel paths from another `AbstractArchSpace` where number of paths is uniformly drawn from `r`.

Input vertex is input to a number of paths drawn from an `AbstractParSpace`. Concatenation of paths is output.
"""
struct ForkArchSpace <:AbstractArchSpace
    s::AbstractArchSpace
    p::AbstractParSpace{1, <:Integer}
    c::ConcConf
end
ForkArchSpace(s::AbstractArchSpace, r::Integer; conf=ConcConf()) = ForkArchSpace(s, SingletonParSpace(r), conf)
ForkArchSpace(s::AbstractArchSpace, r::AbstractVector{<:Integer}; conf=ConcConf()) = ForkArchSpace(s, ParSpace(r), conf)

function (s::ForkArchSpace)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit())
     # Make sure there are no paths with size 0, which is what happens if np > outsize
     np=min_nomissing(s.p(rng), outsize)
     np == 0 && return in
     outsizes = eq_split(outsize, np)
     return s.c(map(i -> s.s(in, rng, outsize=outsizes[i], wi=forkinitW(wi, outsizes, i)), 1:np))
 end
function (s::ForkArchSpace)(name::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    # Make sure there are no paths with size 0, which is what happens if np > outsize
    np=min_nomissing(s.p(rng), outsize)
    np == 0 && return in
    outsizes = eq_split(outsize, np)
    return s.c(name * ".cat", map(i -> s.s(join([name, ".path", i]), in, rng,outsize=outsizes[i], wi=forkinitW(wi, outsizes, i)), 1:np))
end
min_nomissing(x, ::Missing) = x
min_nomissing(x, y) = min(x,y)

function eq_split(x, n)
    rem = x
    out = Vector(undef, n)
    for i = 1:n
        out[i] = rem ÷ (n-i+1)
        rem -= out[i]
    end
    return out
end

forkinitW(wi::AbstractWeightInit, outsizes, i) = wi
forkinitW(wi::IdentityWeightInit, outsizes, i) = PartialIdentityWeightInit(mapreduce(ii -> outsizes[ii], + , 1:i-1, init=0), 0)

"""
    ResidualArchSpace <:AbstractArchSpace
    ResidualArchSpace(s::AbstractArchSpace, [conf::VertexConf])

Turns the wrapped `AbstractArchSpace` into a residual.

Return `x = y + in` where `y` is drawn from the wrapped `AbstractArchSpace` when invoked with `in` as input vertex.

`conf` is used to decorate trait and wrap the computation (e.g in an `ActivationContribution`).
"""
struct ResidualArchSpace <:AbstractArchSpace
    s::AbstractArchSpace
    conf::VertexConf
end
ResidualArchSpace(s::AbstractArchSpace) = ResidualArchSpace(s, VertexConf(traitdecoration = validated() ∘ default_logging(), outwrap = ActivationContribution))
ResidualArchSpace(l::AbstractLayerSpace) = ResidualArchSpace(VertexSpace(l))

function (s::ResidualArchSpace)(in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    add = s.conf >> in + s.s(in, rng,outsize=nout(in), wi=resinitW(wi))
    return resscale(wi, s.conf, add)
end
function (s::ResidualArchSpace)(name::String, in::AbstractVertex, rng=rng_default; outsize=missing, wi=DefaultWeightInit())
    conf = s.conf
    nconf = @set conf.traitdecoration = conf.traitdecoration ∘ named(name * ".add")
    add = nconf >> in + s.s(join([name, ".res"]), in, rng,outsize=nout(in), wi=resinitW(wi))

    sconf = @set conf.traitdecoration = conf.traitdecoration ∘ named(name * ".add.scale")
    return resscale(wi, sconf, add)
 end

resinitW(wi::AbstractWeightInit) = wi

resscale(wi, conf, v) = v
resscale(::Union{IdentityWeightInit, PartialIdentityWeightInit}, conf, v) = invariantvertex(conf.outwrap(x -> convert(eltype(x), 0.5) .* x), v; traitdecoration= conf.traitdecoration)

"""
    FunctionSpace <: AbstractArchSpace
    FunctionSpace(funs...; namesuff::String, conf=LayerVertexConf(ActivationContribution, validated() ∘ default_logging()))

Return a `SizeInvariant` vertex representing `fun(x)` when invoked with `in` as input vertex where `x` is output of `in` where `fun` is uniformly selected from `funs`.
"""
struct FunctionSpace{C<:LayerVertexConf, P<:AbstractParSpace} <: AbstractArchSpace
    conf::C
    funspace::P
    namesuff::String
end
funspace_default_conf() = LayerVertexConf(ActivationContribution, validated() ∘ default_logging())

FunctionSpace(funs...; namesuff::String, conf=funspace_default_conf()) = FunctionSpace(conf, ParSpace1D(funs...), namesuff)

(s::FunctionSpace)(in::AbstractVertex, rng=rng_default; outsize=nothing, wi=nothing) = funvertex(s, in, rng)
(s::FunctionSpace)(name::String, in::AbstractVertex, rng=rng_default; outsize=nothing, wi=nothing) = funvertex(join([name,s.namesuff]), s, in, rng)

funvertex(s::FunctionSpace, in::AbstractVertex, rng) = invariantvertex(s.conf.layerfun(s.funspace(rng)), in, traitdecoration = s.conf.traitfun)

funvertex(name::String, s::FunctionSpace, in::AbstractVertex, rng) =
invariantvertex(s.conf.layerfun(s.funspace(rng)), in, traitdecoration = s.conf.traitfun ∘ named(name))

"""
    GlobalPoolSpace(Ts...)
    GlobalPoolSpace(conf::LayerVertexConf, Ts...)

Short for `FunctionSpace` with global average or global max pooling.

Also adds a `MutationShield` to prevent the vertex from being removed by default.
"""
GlobalPoolSpace(Ts...) = GlobalPoolSpace(LayerVertexConf(ActivationContribution, MutationShield ∘ validated() ∘ default_logging()), Ts...)
GlobalPoolSpace(conf::LayerVertexConf, Ts...=(MaxPool, MeanPool)...) = FunctionSpace(GlobalPool.(Ts)..., namesuff = ".globpool", conf=conf)


text:
using ExtraFrech
using Test

@testset "ExtraFrech.jl" begin
    # Write your tests here.
    @test ExtraFrech.impurity(zeros(10)) == 0

end


text:
using Literate

function preprocess(str)
    str = replace(str, "# PREAMBLE" => "")
    str = replace(str, "# PKG_SETUP" =>
    """
    using Pkg
    Pkg.activate(@__DIR__)
    Pkg.instantiate()
    """)
    return str
end

exampledir = @__DIR__
for subdir in ["simulation", "trajectory_optimization", "reinforcement_learning", "gradients", "mechanics"]
    root = joinpath(exampledir, subdir)
    isdir(root) || continue
    @show subdir
    for file in readdir(root)
        name, ext = splitext(file)
        lowercase(ext) == ".jl" || continue
        absfile = joinpath(root, file)
        @show absfile
        Literate.notebook(absfile, root, execute=false, preprocess=preprocess)
    end
end


text:
@doc raw"""
    change_representer(M::Hyperbolic{n}, ::EuclideanMetric, p, X)

Change the Eucliden representer `X` of a cotangent vector at point `p`.
We only have to correct for the metric, which means that the sign of the last entry changes, since
for the result ``Y``  we are looking for a tangent vector such that

```math
    g_p(Y,Z) = -y_{n+1}z_{n+1} + \sum_{i=1}^n y_iz_i = \sum_{i=1}^{n+1} z_ix_i
```

holds, which directly yields ``y_i=x_i`` for ``i=1,\ldots,n`` and ``y_{n+1}=-x_{n+1}``.
"""
change_representer(::Hyperbolic, ::EuclideanMetric, ::Any, ::Any)

function change_representer!(M::Hyperbolic, Y, ::EuclideanMetric, p, X)
    copyto!(M, Y, p, X)
    Y[end] *= -1
    return Y
end

function change_metric!(::Hyperbolic, ::Any, ::EuclideanMetric, ::Any, ::Any)
    return error(
        "Changing metric from Euclidean to Hyperbolic is not possible (see Sylvester's law of inertia).",
    )
end

function check_point(M::Hyperbolic, p; kwargs...)
    if !isapprox(minkowski_metric(p, p), -1.0; kwargs...)
        return DomainError(
            minkowski_metric(p, p),
            "The point $(p) does not lie on $(M) since its Minkowski inner product is not -1.",
        )
    end
    return nothing
end

function check_vector(M::Hyperbolic, p, X; kwargs...)
    if !isapprox(minkowski_metric(p, X), 0.0; kwargs...)
        return DomainError(
            abs(minkowski_metric(p, X)),
            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal (with respect to the Minkowski inner product) in the embedding.",
        )
    end
    return nothing
end

function convert(::Type{HyperboloidTVector}, X::T) where {T<:AbstractVector}
    return HyperboloidTVector(X)
end
function convert(
    ::Type{HyperboloidTVector},
    p::P,
    X::T,
) where {P<:AbstractVector,T<:AbstractVector}
    return HyperboloidTVector(X)
end
convert(::Type{AbstractVector}, X::HyperboloidTVector) = X.value
function convert(
    ::Type{T},
    p::HyperboloidPoint,
    X::HyperboloidTVector,
) where {T<:AbstractVector}
    return X.value
end

function convert(::Type{HyperboloidPoint}, p::T) where {T<:AbstractVector}
    return HyperboloidPoint(p)
end
convert(::Type{AbstractVector}, p::HyperboloidPoint) = p.value

@doc raw"""
    convert(::Type{HyperboloidPoint}, p::PoincareBallPoint)
    convert(::Type{AbstractVector}, p::PoincareBallPoint)

convert a point [`PoincareBallPoint`](@ref) `x` (from $ℝ^n$) from the
Poincaré ball model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ to a [`HyperboloidPoint`](@ref) $π(p) ∈ ℝ^{n+1}$.
The isometry is defined by

````math
π(p) = \frac{1}{1-\lVert p \rVert^2}
\begin{pmatrix}2p_1\\⋮\\2p_n\\1+\lVert p \rVert^2\end{pmatrix}
````

Note that this is also used, when the type to convert to is a vector.
"""
function convert(::Type{HyperboloidPoint}, p::PoincareBallPoint)
    return HyperboloidPoint(convert(AbstractVector, p))
end
function convert(::Type{AbstractVector}, p::PoincareBallPoint)
    return 1 / (1 - norm(p.value)^2) .* vcat(2 .* p.value, 1 + norm(p.value)^2)
end

@doc raw"""
    convert(::Type{HyperboloidPoint}, p::PoincareHalfSpacePoint)
    convert(::Type{AbstractVector}, p::PoincareHalfSpacePoint)

convert a point [`PoincareHalfSpacePoint`](@ref) `p` (from $ℝ^n$) from the
Poincaré half plane model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ to a [`HyperboloidPoint`](@ref) $π(p) ∈ ℝ^{n+1}$.

This is done in two steps, namely transforming it to a Poincare ball point and from there further on to a Hyperboloid point.
"""
function convert(t::Type{HyperboloidPoint}, p::PoincareHalfSpacePoint)
    return convert(t, convert(PoincareBallPoint, p))
end
function convert(t::Type{AbstractVector}, p::PoincareHalfSpacePoint)
    return convert(t, convert(PoincareBallPoint, p))
end

@doc raw"""
    convert(::Type{HyperboloidTVector}, p::PoincareBallPoint, X::PoincareBallTVector)
    convert(::Type{AbstractVector}, p::PoincareBallPoint, X::PoincareBallTVector)

Convert the [`PoincareBallTVector`](@ref) `X` from the tangent space at `p` to a
[`HyperboloidTVector`](@ref) by computing the push forward of the isometric map, cf.
[`convert(::Type{HyperboloidPoint}, p::PoincareBallPoint)`](@ref).

The push forward $π_*(p)$ maps from $ℝ^n$ to a subspace of $ℝ^{n+1}$, the formula reads

````math
π_*(p)[X] = \begin{pmatrix}
    \frac{2X_1}{1-\lVert p \rVert^2} + \frac{4}{(1-\lVert p \rVert^2)^2}⟨X,p⟩p_1\\
    ⋮\\
    \frac{2X_n}{1-\lVert p \rVert^2} + \frac{4}{(1-\lVert p \rVert^2)^2}⟨X,p⟩p_n\\
    \frac{4}{(1-\lVert p \rVert^2)^2}⟨X,p⟩
\end{pmatrix}.
````
"""
function convert(::Type{HyperboloidTVector}, p::PoincareBallPoint, X::PoincareBallTVector)
    return HyperboloidTVector(convert(AbstractVector, p, X))
end
function convert(
    ::Type{T},
    p::PoincareBallPoint,
    X::PoincareBallTVector,
) where {T<:AbstractVector}
    t = (1 - norm(p.value)^2)
    den = 4 * dot(p.value, X.value) / (t^2)
    c1 = (2 / t) .* X.value + den .* p.value
    return vcat(c1, den)
end

@doc raw"""
    convert(
        ::Type{Tuple{HyperboloidPoint,HyperboloidTVector}}.
        (p,X)::Tuple{PoincareBallPoint,PoincareBallTVector}
    )
    convert(
        ::Type{Tuple{P,T}},
        (p, X)::Tuple{PoincareBallPoint,PoincareBallTVector},
    ) where {P<:AbstractVector, T <: AbstractVector}

Convert a [`PoincareBallPoint`](@ref) `p` and a [`PoincareBallTVector`](@ref) `X`
to a [`HyperboloidPoint`](@ref) and a [`HyperboloidTVector`](@ref) simultaneously,
see [`convert(::Type{HyperboloidPoint}, ::PoincareBallPoint)`](@ref) and
[`convert(::Type{HyperboloidTVector}, ::PoincareBallPoint, ::PoincareBallTVector)`](@ref)
for the formulae.
"""
function convert(
    ::Type{Tuple{HyperboloidPoint,HyperboloidTVector}},
    (p, X)::Tuple{PoincareBallPoint,PoincareBallTVector},
)
    return (convert(HyperboloidPoint, p), convert(HyperboloidTVector, p, X))
end

@doc raw"""
    convert(::Type{HyperboloidTVector}, p::PoincareHalfSpacePoint, X::PoincareHalfSpaceTVector)
    convert(::Type{AbstractVector}, p::PoincareHalfSpacePoint, X::PoincareHalfSpaceTVector)

convert a point [`PoincareHalfSpaceTVector`](@ref) `X` (from $ℝ^n$) at `p` from the
Poincaré half plane model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$ to a
[`HyperboloidTVector`](@ref) $π(p) ∈ ℝ^{n+1}$.

This is done in two steps, namely transforming it to a Poincare ball point and from there further on to a Hyperboloid point.
"""
function convert(
    t::Type{HyperboloidTVector},
    p::PoincareHalfSpacePoint,
    X::PoincareHalfSpaceTVector,
)
    return convert(t, convert(Tuple{PoincareBallPoint,PoincareBallTVector}, (p, X))...)
end
function convert(
    t::Type{T},
    p::PoincareHalfSpacePoint,
    X::PoincareHalfSpaceTVector,
) where {T<:AbstractVector}
    return convert(t, convert(Tuple{PoincareBallPoint,PoincareBallTVector}, (p, X))...)
end

@doc raw"""
    convert(
        ::Type{Tuple{HyperboloidPoint,HyperboloidTVector},
        (p,X)::Tuple{PoincareHalfSpacePoint, PoincareHalfSpaceTVector}
    )
    convert(
        ::Type{Tuple{T,T},
        (p,X)::Tuple{PoincareHalfSpacePoint, PoincareHalfSpaceTVector}
    ) where {T<:AbstractVector}

convert a point [`PoincareHalfSpaceTVector`](@ref) `X` (from $ℝ^n$) at `p` from the
Poincaré half plane model of the [`Hyperbolic`](@ref) manifold $\mathcal H^n$
to a tuple of a [`HyperboloidPoint`](@ref) and a [`HyperboloidTVector`](@ref) $π(p) ∈ ℝ^{n+1}$
simultaneously.

This is done in two steps, namely transforming it to the Poincare ball model and from there
further on to a Hyperboloid.
"""
function convert(
    t::Type{Tuple{HyperboloidPoint,HyperboloidTVector}},
    (p, X)::Tuple{PoincareHalfSpacePoint,PoincareHalfSpaceTVector},
)
    return convert(t, convert(Tuple{PoincareBallPoint,PoincareBallTVector}, (p, X)))
end

@doc raw"""
    distance(M::Hyperbolic, p, q)
    distance(M::Hyperbolic, p::HyperboloidPoint, q::HyperboloidPoint)

Compute the distance on the [`Hyperbolic`](@ref) `M`, which reads

````math
d_{\mathcal H^n}(p,q) = \operatorname{acosh}( - ⟨p, q⟩_{\mathrm{M}}),
````

where $⟨\cdot,\cdot⟩_{\mathrm{M}}$ denotes the [`MinkowskiMetric`](@ref) on the embedding,
the [`Lorentz`](@ref)ian manifold.
"""
distance(::Hyperbolic, p, q) = acosh(max(-minkowski_metric(p, q), 1.0))

embed(M::Hyperbolic, p::HyperboloidPoint) = embed(M, p.value)
embed!(M::Hyperbolic, q, p::HyperboloidPoint) = embed!(M, q, p.value)
function embed(M::Hyperbolic, p::HyperboloidPoint, X::HyperboloidTVector)
    return embed(M, p.value, X.value)
end
function embed!(M::Hyperbolic, Y, p::HyperboloidPoint, X::HyperboloidTVector)
    return embed!(M, Y, p.value, X.value)
end

function exp!(M::Hyperbolic, q, p, X)
    vn = sqrt(max(inner(M, p, X, X), 0.0))
    vn < eps(eltype(p)) && return copyto!(q, p)
    return copyto!(q, cosh(vn) * p + sinh(vn) / vn * X)
end

# overwrite the default construction on level 2 (dispatching on basis)
# since this function should not call get_vector (that relies on get_basis itself on H2)
function _get_basis(
    M::Hyperbolic,
    p,
    B::DefaultOrthonormalBasis{ℝ,TangentSpaceType};
    kwargs...,
)
    return get_basis_orthonormal(M, p, ℝ)
end

function get_basis_orthonormal(M::Hyperbolic{n}, p, r::RealNumbers) where {n}
    V = [
        _hyperbolize(M, p, [i == k ? one(eltype(p)) : zero(eltype(p)) for k in 1:n]) for
        i in 1:n
    ]
    return CachedBasis(DefaultOrthonormalBasis(r), gram_schmidt(M, p, V))
end

function get_basis_diagonalizing(M::Hyperbolic, p, B::DiagonalizingOrthonormalBasis)
    n = manifold_dimension(M)
    X = B.frame_direction
    V = [
        _hyperbolize(M, p, [i == k ? one(eltype(p)) : zero(eltype(p)) for k in 1:n]) for
        i in 1:n
    ]
    κ = -ones(n)
    if norm(M, p, X) != 0
        placed = false
        for i in 1:n
            if abs(inner(M, p, X, V[i])) ≈ norm(M, p, X) # is X a multiple of V[i]?
                V[i] .= V[1]
                V[1] .= X
                placed = true
                break
            end
        end
        if !placed
            V[1] .= X
        end
        κ[1] = 0.0
    end
    V = gram_schmidt(M, p, V; atol=4 * eps(eltype(V[1])))
    return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, κ, V))
end

@doc raw"""
    get_coordinates(M::Hyperbolic, p, X, ::DefaultOrthonormalBasis)

Compute the coordinates of the vector `X` with respect to the orthogonalized version of
the unit vectors from $ℝ^n$, where $n$ is the manifold dimension of the [`Hyperbolic`](@ref)
 `M`, utting them intop the tangent space at `p` and orthonormalizing them.
"""
get_coordinates(M::Hyperbolic, p, X, ::DefaultOrthonormalBasis)

function get_coordinates_orthonormal(M::Hyperbolic, p, X, r::RealNumbers)
    return get_coordinates(M, p, X, get_basis_orthonormal(M, p, r))
end
function get_coordinates_orthonormal!(M::Hyperbolic, c, p, X, r::RealNumbers)
    c = get_coordinates!(M, c, p, X, get_basis_orthonormal(M, p, r))
    return c
end
function get_coordinates_diagonalizing!(
    M::Hyperbolic,
    c,
    p,
    X,
    B::DiagonalizingOrthonormalBasis,
)
    c = get_coordinates!(M, c, p, X, get_basis_diagonalizing(M, p, B))
    return c
end

@doc raw"""
    get_vector(M::Hyperbolic, p, c, ::DefaultOrthonormalBasis)

Compute the vector from the coordinates with respect to the orthogonalized version of
the unit vectors from $ℝ^n$, where $n$ is the manifold dimension of the [`Hyperbolic`](@ref)
 `M`, utting them intop the tangent space at `p` and orthonormalizing them.
"""
get_vector(M::Hyperbolic, p, c, ::DefaultOrthonormalBasis)

function get_vector_orthonormal!(M::Hyperbolic, X, p, c, r::RealNumbers)
    X = get_vector!(M, X, p, c, get_basis(M, p, DefaultOrthonormalBasis(r)))
    return X
end
function get_vector!(M::Hyperbolic, X, p, c, B::DiagonalizingOrthonormalBasis)
    X = get_vector!(M, X, p, c, get_basis(M, p, B))
    return X
end

@doc raw"""
    _hyperbolize(M,q)

Given the [`Hyperbolic`](@ref)`(n)` manifold using the hyperboloid model, a point from the
$q\in ℝ^n$ can be set onto the manifold by computing its last component such that for the
resulting `p` we have that its [`minkowski_metric`](@ref) is ``⟨p,p⟩_{\mathrm{M}} = - 1``,
i.e. ``p_{n+1} = \sqrt{\lVert q \rVert^2 - 1}``
"""
_hyperbolize(::Hyperbolic, q) = vcat(q, sqrt(norm(q)^2 + 1))

@doc raw"""
    _hyperbolize(M, p, Y)

Given the [`Hyperbolic`](@ref)`(n)` manifold using the hyperboloid model and a point `p`
thereon, we can put a vector $Y\in ℝ^n$  into the tangent space by computing its last
component such that for the
resulting `p` we have that its [`minkowski_metric`](@ref) is $⟨p,X⟩_{\mathrm{M}} = 0$,
i.e. $X_{n+1} = \frac{⟨\tilde p, Y⟩}{p_{n+1}}$, where $\tilde p = (p_1,\ldots,p_n)$.
"""
_hyperbolize(::Hyperbolic, p, Y) = vcat(Y, dot(p[1:(end - 1)], Y) / p[end])

@doc raw"""
    inner(M::Hyperbolic{n}, p, X, Y)
    inner(M::Hyperbolic{n}, p::HyperboloidPoint, X::HyperboloidTVector, Y::HyperboloidTVector)

Cmpute the inner product in the Hyperboloid model, i.e. the [`minkowski_metric`](@ref) in
the embedding. The formula reads

````math
g_p(X,Y) = ⟨X,Y⟩_{\mathrm{M}} = -X_{n}Y_{n} + \displaystyle\sum_{k=1}^{n-1} X_kY_k.
````
This employs the metric of the embedding, see [`Lorentz`](@ref) space.
"""
inner(M::Hyperbolic, p, X, Y)

function log!(M::Hyperbolic, X, p, q)
    scp = minkowski_metric(p, q)
    w = q + scp * p
    wn = sqrt(max(scp .^ 2 - 1, 0.0))
    wn < eps(eltype(p)) && return zero_vector!(M, X, p)
    X .= acosh(max(1.0, -scp)) / wn .* w
    return X
end

function minkowski_metric(a::HyperboloidPoint, b::HyperboloidPoint)
    return minkowski_metric(a.value, b.value)
end
function minkowski_metric(a::HyperboloidTVector, b::HyperboloidPoint)
    return minkowski_metric(a.value, b.value)
end
function minkowski_metric(a::HyperboloidPoint, b::HyperboloidTVector)
    return minkowski_metric(a.value, b.value)
end
function minkowski_metric(a::HyperboloidTVector, b::HyperboloidTVector)
    return minkowski_metric(a.value, b.value)
end

project!(::Hyperbolic, Y, p, X) = (Y .= X .+ minkowski_metric(p, X) .* p)
function project!(
    ::Hyperbolic,
    Y::HyperboloidTVector,
    p::HyperboloidPoint,
    X::HyperboloidTVector,
)
    return (Y.value .= X.value .+ minkowski_metric(p.value, X.value) .* p.value)
end

function Random.rand!(M::Hyperbolic{N}, pX; vector_at=nothing, σ=one(eltype(pX))) where {N}
    if vector_at === nothing
        a = σ .* randn(N)
        pX[1:(end - 1)] .= a
        pX[end] = sqrt(1 + dot(a, a))
    else
        Y = σ * randn(eltype(vector_at), size(vector_at))
        project!(M, pX, vector_at, Y)
    end
    return pX
end
function Random.rand!(
    rng::AbstractRNG,
    M::Hyperbolic{N},
    pX;
    vector_at=nothing,
    σ=one(eltype(pX)),
) where {N}
    if vector_at === nothing
        a = σ .* randn(rng, N)
        pX[1:(end - 1)] .= a
        pX[end] = sqrt(1 + dot(a, a))
    else
        Y = σ * randn(rng, eltype(vector_at), size(vector_at))
        project!(M, pX, vector_at, Y)
    end
    return pX
end

function parallel_transport_to!(M::Hyperbolic, Y, p, X, q)
    w = log(M, p, q)
    wn = norm(M, p, w)
    wn < eps(eltype(p + q)) && return copyto!(Y, X)
    return copyto!(Y, X - (inner(M, p, w, X) * (w + log(M, q, p)) / wn^2))
end

max_stars_repo_name: JuliaNLSolvers/ManifoldMuseum

max_issues_repo_name: JuliaNLSolvers/ManifoldMuseum

max_forks_repo_name: JuliaNLSolvers/ManifoldMuseum


text:
# some commentary

# docstring...

mutable struct WeaveMol <: AbstractAtoms # maybe we want to add AbstractAtomGraph <: AbstractAtoms that this and AtomGraph would subtype?
    smiles::String
    elements::Vector{String}
    graph::SimpleWeightedGraph{<:Integer,<:Real}
    encoded_atom_features::Union{Matrix{<:Real},Nothing}
    encoded_pair_features::Union{Matrix{Union{Vector{<:Real},Nothing}},Nothing}
    featurization::WeaveFeaturization
    id::Any # probably makes sense to have this in addition to smiles? Maybe?
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.7700820821577 <= q <= 37.6713142442229)

                     add_NL_constraint(m, :(sigmoid(0.9266287961121966 + 0.9630341445616355*sigmoid(0.9080159995418975 + 0.10769169268185319*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.8197021439185934*sigmoid(-0.9264374759333807 + -0.7028105280156569*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.8564198456770362*sigmoid(-0.7118497063916944 + 0.36586202274128254*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7868467121279954*sigmoid(0.8429853377330745 + -0.6377098680254001*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + sigmoid(-0.9209290943615045 + -0.04560727661330999*sigmoid(0.9080159995418975 + 0.10769169268185319*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.04187997574232316*sigmoid(-0.9264374759333807 + -0.7028105280156569*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.009897273687461627*sigmoid(-0.7118497063916944 + 0.36586202274128254*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.4138607422831795*sigmoid(0.8429853377330745 + -0.6377098680254001*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + sigmoid(0.41055131347137674 + 0.5074173875637427*sigmoid(0.9080159995418975 + 0.10769169268185319*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7865169798440914*sigmoid(-0.9264374759333807 + -0.7028105280156569*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.6322955384333278*sigmoid(-0.7118497063916944 + 0.36586202274128254*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.03598168048021799*sigmoid(0.8429853377330745 + -0.6377098680254001*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + sigmoid(-0.7045693885767101 + 0.765322888364711*sigmoid(0.9080159995418975 + 0.10769169268185319*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.9087169409161477*sigmoid(-0.9264374759333807 + -0.7028105280156569*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.29202574365958966*sigmoid(-0.7118497063916944 + 0.36586202274128254*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.14636407330293322*sigmoid(0.8429853377330745 + -0.6377098680254001*sigmoid(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*sigmoid(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*sigmoid(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*sigmoid(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using MVApp

data_folder = "data/0429"
main(data_folder,97,1,10)
main(data_folder,97,-1,-1)

# main(data_folder,95,1,10)
push!(ARGS, data_folder, "1","10", "97", "56")
MVApp.julia_main()


text:
# TODO:
# * defm, set_defm!         >>> decide what to do about those
# * Base.read, Base.write   >>> rename as read_xyz, write_xyz
# * if we make Atoms immutable, then we could make the calculator
#   etc a type parameter, but it has the disadvantage that we would need
#   to create a new object every time we change it significantly.
#   maybe a good thing.
#
# * `JData` has a hard-coded floating point type. this seems necessary
#       because otherwise the kw-constructor doesn't work. probably this is
#       not important, but could be fixed at some point (low priority)
#
# * the kwconstructor only allows Float64, this is a bad shortcoming that
#   needs to be fixed I think
#
# * decide whether == should check equality of data
#
# * implement a softer `isequal` which checks for "physical" equality, i.e.,
#   it accounts for PBC, ignores any data, and allows for a controllable
#   error in the positions, momenta and masses
#
# * add fields for charges and magnetic moments?
#

const CH = JuLIP.Chemistry

import Base.Dict

export Atoms,
       bulk,
       cell_vecs,
       chemical_symbols,
       atomic_numbers, numbers,
       chemical_symbols

"""
`JData`: datatype for storing any data

some data which needs to be updated if the configuration (positions only!) has
changed too much.
"""
mutable struct JData
   max_change::Float64     # how much X may change before recomputing
   accum_change::Float64   # how much has it changed already
   data::Any
end

"""
`Atoms{T <: AbstractFloat} <: AbstractAtoms`

The main type storing information about atomic configurations. This is
not normally constructed directly by calling `Atoms`, but using one of the
utility functions such as `bulk`.

## Constructors

TODO

## Convenience functions

* `bulk`

# Getters and Setters

* `positions`, `set_positions!`
* `momenta`, `set_momenta!`
* `masses`, `set_masses!`
* `numbers`, `set_numbers!`
* TODO: `chemical_symbols`, `set_chemical_symbols!`
* `cell`, `set_cell!`
* `pbc`, `set_pbc!`
* `calculator`, `set_calculator!`
* `constraint`, `set_constraint!`
* `get_data`, `set_data!`
"""
@with_kw mutable struct Atoms{T <: AbstractFloat, TI <: Integer} <: AbstractAtoms
   X::Vector{JVec{T}} = JVec{Float64}[]       # positions
   P::Vector{JVec{T}} = JVec{Float64}[]       # momenta (or velocities?)
   M::Vector{T} = Float64[]             # masses
   Z::Vector{TI} = Int[]           # atomic numbers
   cell::JMat{T} = zero(JMat{T})                   # cell
   pbc::JVec{Bool} = JVec(false, false, false)     # boundary condition
   calc::Union{Void, AbstractCalculator} = nothing
   cons::Union{Void, AbstractConstraint} = nothing
   data::Dict{Any,JData} = Dict{Any,JData}()
end

_auto_pbc(pbc::Tuple{Bool, Bool, Bool}) = pbc
_auto_pbc(pbc::Bool) = (pbc, pbc, pbc)
_auto_pbc(pbc::AbstractVector) = tuple(pbc...)

_auto_cell(cell) = JMat(cell)
_auto_cell(C::Vector{Any}) = JMatF([ C[1] C[2] C[3] ])

Atoms(X, P, M, Z, cell, pbc; calc=NullCalculator(),
      cons = NullConstraint(), data = Dict{Any,JData}()) =
   Atoms(X, P, M, Z, _auto_cell(cell), _auto_pbc(pbc), calc, cons, data)

Atoms(Z::Vector{TI}, X::Vector{JVec{T}}; kwargs...
     ) where {TI <: Integer, T <: AbstractFloat} =
  Atoms(;Z=Z, X=X, P = zeros(JVec{T}, length(X)), M = zeros(length(X)), kwargs...)

# derived properties
length(at::Atoms) = length(at.X)

# ------------------------ access to struct fields ----------------------

# ----- getters

positions(at::Atoms) = copy(at.X)
momenta(at::Atoms) = copy(at.P)
masses(at::Atoms) = copy(at.M)
numbers(at::Atoms) = copy(at.Z)
atomic_numbers(at::Atoms) = copy(at.Z)
cell(at::Atoms) = at.cell
pbc(at::Atoms) = at.pbc
calculator(at::Atoms) = at.calc
constraint(at::Atoms) = at.cons

chemical_symbols(at::Atoms) = Chemistry.chemical_symbol.(at.Z)

# ----- setters

function set_positions!(at::Atoms{T}, X::Vector{JVec{T}})  where T
   update_data!(at, dist(at, X))
   at.X .= X
   return at
end
function set_momenta!(at::Atoms{T}, P::Vector{JVec{T}})  where T
   at.P .= P
   return at
end
function set_masses!(at::Atoms{T}, M::Vector{T})  where T
   update_data!(at, Inf)
   at.M .= M
   return at
end
function set_numbers!(at::Atoms{T, TI}, Z::Vector{TI}) where T where TI
   update_data!(at, Inf)
   at.Z .= Z
   return at
end
function set_cell!(at::Atoms{T}, C::AbstractMatrix) where T
   at.cell = JMat{T}(C)
   return at
end
function set_pbc!(at::Atoms, p::Union{AbstractVector, Tuple})
   at.pbc = JVec{Bool}(p...)
   return at
end
function set_calculator!(at::Atoms, calc::AbstractCalculator)
   at.calc = calc
   return at
end
function set_constraint!(at::Atoms, cons::AbstractConstraint)
   at.cons = cons
   return at
end

# --------------- equality tests -------------------

# we need to sort vectors of vectors to test equality of configurations
# note this is a strict notion of equality
function Base.isless(x::JVec, y::JVec)
   if x[1] < y[1]
      return true
   elseif x[1] == y[1]
      if x[2] < y[2]
         return true
      elseif x[2] == y[2]
         if x[3] < y[3]
            return true
         end
      end
   end
   return false
end

import Base.==
==(at1::Atoms{T,TI}, at2::Atoms{T,TI}) where T where TI =
   isapprox(at1, at2, tol = 0.0) && (at1.data == at2.data)

import Base.isapprox
function isapprox(at1::Atoms{T,TI}, at2::Atoms{T,TI}; tol = sqrt(eps(T))) where T where TI
   if length(at1) != length(at2)
      return false
   end
   if tol > 0
      ndigits = floor(Int, abs(log10(tol)))
      X1 = [round.(x, ndigits) for x in at1.X]
      X2 = [round.(x, ndigits) for x in at2.X]
   else
      X1, X2 = at1.X, at2.X
   end
   p1 = sortperm(X1)
   p2 = sortperm(X2)
   return (maxdist(at1.X[p1], at2.X[p2]) <= tol) &&
          (maxdist(at1.P[p1], at2.P[p2]) <= tol) # &&
          (vecnorm(at1.M[p1] - at2.M[p2], Inf) <= tol) &&
          (at1.Z[p1] == at2.Z[p2]) &&
          (vecnorm(at1.cell - at2.cell, Inf) <= tol) &&
          (at1.pbc == at2.pbc) &&
          (at1.calc == at2.calc) &&
          (at1.cons == at2.cons)
end

# --------------- some aliases for easier / direct access --------------
#                 and some access not covered above

# access positions by direct indexing
Base.getindex(at::Atoms, i::Integer) = at.X[i]
function Base.setindex!(at::Atoms{T}, x::JVec, i::Integer) where T <: AbstractFloat
   at.X[i] = JVec{T}(x)
   return at.X[i]
end

"""
`cell_vecs(at::Atoms)` : return the three cell vectors
"""
cell_vecs(at::Atoms) = at.cell[1,:], at.cell[2,:], at.cell[3,:]

function neighbourlist(at::Atoms{T}, cutoff::T; recompute=false, kwargs...) where T <: AbstractFloat
   # TODO: re-design this from scratch . . .
   PairList(positions(at), cutoff, cell(at), pbc(at); kwargs...)
   # if !has_data(at, (:nlist, cutoff)) || recompute
   #    set_transient!(at, (:nlist, cutoff),
   #          PairList(positions(at), cutoff, cell(at), pbc(at))
   #       )
   # end
   # return get_data(at, (:nlist, cutoff))
end

neighbourlist(at::Atoms) = neighbourlist(at, cutoff(at))

"""
`static_neighbourlist(at::Atoms, cutoff)`

This function first checks whether a static neighbourlist already exists
with cutoff `cutoff` and if it does then it returns the existing list.
If it does not, then it computes a new neighbour list with the current
configuration, stores it for later use and returns it.
"""
function static_neighbourlist(at::Atoms{T}, cutoff::T) where T
   if !has_data(at, (:snlist, cutoff))
      set_data!( at, (:snlist, cutoff),
            PairList(positions(at), cutoff, cell(at), pbc(at))
         )
   end
   return get_data(at, (:snlist, cutoff))
end

static_neighbourlist(at::Atoms) = static_neighbourlist(at, cutoff(at))

# --------------- data Dict handling ------------------

"""
`has_data(at::Atoms, key)`: checks whether `at` is storing a datum with key `key`; see also `set_data!`
"""
has_data(at::Atoms, key) = haskey(at.data, key)

"""
`get_data(at::Atoms, key)`: retrieves a datum;
see also `set_data!`. Unlike `get_positions`, etc, this returns the original
datum, not a copy (unless it is immutable).
"""
get_data(at::Atoms, key) = (at.data[key]).data

"""
`set_data!(a::Atoms, key, value, max_change=Inf)`:

Stores `value` under an arbitrary key `key` inside `at`. The entry is paired
with an updatemeasure `chg`. Whenever atom positions or the cell change by a
distance `d`, the counter is incremented `chg += d`. When `chg > max_change`,
the entry is deleted from the dictionary.

The two most common uses are with `max_change = Inf` (default) - in this case
the entry will never be deleted; and `max_change = 0.0` - in this case the entry
will be deleted whenever atom positions or the cell change.

The change counter is used, e.g., for neighbourlists with buffer, and for
preconditioners so that the list need not be recomputed each time the
configuration changes.
"""
function set_data!(at::Atoms, key, value, max_change=Inf)
   at.data[key] = JData(max_change, 0.0, value)
   return at
end

"""
`function set_transient!(a::Atoms, key, value, max_change=0.0)`: this is an
alias for `set_data!` but with default `max_change = 0.0` instead of `max_change =
Inf`.
"""
set_transient!(at::Atoms, key, value, max_change=0.0) =
      set_data!(at, key, value, max_change)

"""
`update_data!(at::Atoms, r::Real)`:
increment the change counter in all stored data and delete the
expired ones.
"""
function update_data!(at::Atoms, r::Real)
   for (key, t) in at.data
      t.accum_change += r
      if t.accum_change >= t.max_change
         delete!(at.data, key)
      end
   end
   return at
end

# ------------------------ workaround for JLD bug  ----------------------

_read_X(X) = X
_read_X(X::Vector{Any}) = [ JVecF(x) for x in X ]
_read_X(X::Matrix) = [ JVecF(X[:,n]) for n = 1:size(X,2) ]

Dict(at::Atoms) =
   Dict( "id"   => "JuLIP.Atoms",
         "X"    => at.X,
         "P"    => at.P,
         "M"    => at.M,
         "Z"    => at.Z,
         "cell" => at.cell,
         "pbc"  => at.pbc,
         "calc" => nothing,
         "cons" => nothing,
         "data" => nothing )
         # "data" => at.data )

Atoms(D::Dict) =
   Atoms( _read_X(D["X"]),
          _read_X(D["P"]),
          Vector{Float64}(D["M"]),
          Vector{Int}(D["Z"]),
          D["cell"],
          D["pbc"] )
          # calc = D["calc"],
          # cons = D["cons"],
          # data = Dict{Any, JData}(),
          # data = D["data"] )


text:
#=
 █████╗  ██████╗██████╗ ██╗
██╔══██╗██╔════╝██╔══██╗██║
███████║██║     ██████╔╝██║
██╔══██║██║     ██╔══██╗██║
██║  ██║╚██████╗██║  ██║███████╗
╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝╚══════╝
File:       rendezvous_review.jl
Author:     Gabriel Barsi Haberfeld, 2020. gbh2@illinois.edu
Function:   This program simulates all results in the paper "Geometry-Informed
            Minimum Risk Rendezvous Algorithm for Heterogeneous Agents in Urban
            Environments"
Instructions:   Run this file in juno with Julia 1.2.0 or later.
Requirements:   JuMP, Ipopt, Plots, LinearAlgebra, BenchmarkTools.
=#

using Plots, LinearAlgebra, Random, Statistics, ColorSchemes, LazySets
using BenchmarkTools, Dates, Measures
using GaussianProcesses
using Distributions
using Random
import Statistics
using Optim
using QuadGK
using Measures

NoiseVar = 0.05
NoiseStd = sqrt(NoiseVar)
NoiseLog = log10(NoiseVar)
NoiseStdLog = log10(NoiseStd)
# setup driver learning problem.
default(size = 0.6 .* [800, 600])
default(palette = :tol_bright)
default(dpi = 300)
default(lw = 2)
default(margin = 1mm)
FontSize = 12
default(xtickfontsize = FontSize)
default(ytickfontsize = FontSize)
default(xguidefontsize = FontSize)
default(yguidefontsize = FontSize)
default(legendfontsize = FontSize)
default(titlefontsize = FontSize)
rng = MersenneTwister(1234)

rectangle(w, h, x, y) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])

function LearnedDeviation(gp, v) #learned deviation function
    predict_y(gp, [v])
end

function LearnDeviationFunction(D, useConst = false, method = "full")
    #=This function takes in the dataset D and outputs a GP.
    #D[1,:] has historical velocities
    #D[2,:] has measured velocities
    Want to learn measured velocities as a function of velocities
    =#

    x = D[:, 1]   #predictors
    y = D[:, 2]   #regressors
    #Select mean and covariance function
    if !useConst
        mFcn = MeanZero()
    else
        mFcn = MeanConst(mean(D[:, 2]))
    end
    #kern = SE(0.0, 0.0)
    kern = Matern(3 / 2, zeros(5), log10(sqrt(5.25)))
    #lik = BernLik()
    logObsNoise = NoiseLog
    if method == "full"
        return GPE(x, y, mFcn, kern, logObsNoise)
    elseif method == "SOR"
        Xu = Matrix(
            quantile(
                x,
                [
                    0.2,
                    0.25,
                    0.3,
                    0.35,
                    0.4,
                    0.45,
                    0.5,
                    0.55,
                    0.6,
                    0.65,
                    0.7,
                    0.98,
                ],
            )',
        )
        X = Matrix(x')
        return GaussianProcesses.SoR(X, Xu, y, mFcn, kern, logObsNoise)
    elseif method == "DTC"
        Xu = Matrix(
            quantile(
                x,
                [
                    0.02,
                    0.2,
                    0.3,
                    0.35,
                    0.4,
                    0.45,
                    0.5,
                    0.55,
                    0.6,
                    0.65,
                    0.8,
                    0.98,
                ],
            )',
        )
        X = Matrix(x')
        return GaussianProcesses.DTC(X, Xu, y, mFcn, kern, logObsNoise)
    else
        error("Invalid GP method")
    end
end

max_stars_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_issues_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_forks_repo_name: gbarsih/Multi-Path-Safe-Rendezvous


text:
using DataFrames, CSV, BenchmarkTools

randnothing(pct_nothing) = let v = rand()
    v < pct_nothing ? nothing : v
end

randnothing(n, pct) = [randnothing(pct) for _ in 1:n]

df1(n) = DataFrame(x = rand(n), y = rand(n), z = rand(n))

df2(n, pct = 0.1) =
    DataFrame(x = randnothing(n, pct), y = randnothing(n, pct), z = randnothing(n, pct))

test1(df) = @benchmark CSV.write("/tmp/tom.csv", $df)

test2(df) = @benchmark CSV.write("/tmp/tom.csv", $df, transform = (col,val) -> something(val, missing))

function test(msg, f, df)
    println(msg)
    for i in 4:7
        rows = 10^i
        println(rows, " rows: ", f(df(rows)))
    end
end

length(ARGS) > 0 || (println("Usage: julia perf_write.jl [master|new]"); exit(1))
command = ARGS[1]

if command == "master"
    test("Testing without nothing values", test1, df1)
elseif command == "new"
    test("Testing without nothing values (default transform)", test1, df1)
    test("Testing without nothing values (to-missing transform)", test2, df1)
    test("Testing with nothing values (to-missing transform)", test2, df2)
else
    error("Unknown command: $command")
end


text:
module LoadingSpinner

using ..Ahorn
using Gtk, Gtk.ShortNames
using Cairo

animationPreviousTime = 0
animationDelay = 5

const baseWidth = ceil(Int, width(Ahorn.Assets.crowIdle))
const baseHeight = ceil(Int, height(Ahorn.Assets.crowIdle))

const startFrame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowIdle), baseWidth * 4, baseHeight * 4)
const idleFrame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowIdle), baseWidth * 4, baseHeight * 4)
const peck1Frame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowPeck1), baseWidth * 4, baseHeight * 4)
const peck2Frame = Ahorn.scalePixbufSimple(Ahorn.pixbufFromSurface(Ahorn.Assets.crowPeck2), baseWidth * 4, baseHeight * 4)

# Peck twice, return to idle position
const frames = Tuple{Gtk.Pixbuf, Number}[
    (peck1Frame, 0.08),
    (peck2Frame, 0.08),
    (idleFrame, 0.08),
    (peck1Frame, 0.08),
    (peck2Frame, 0.08),
    (idleFrame, 0.08),
]

function animate(dlg::Gtk.GObject)
    global animationPreviousTime

    if time() > animationPreviousTime + animationDelay
        for (frame, delay) in frames
            Ahorn.setProgressDialogPixbuf!(dlg, frame)
            sleep(delay)
        end

        animationPreviousTime = time()

        return true
    end

    return false
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn