text:
module File

importall Images
importall ImmutableArrays
importall FileSystem
importall Textures
import GL
import BSP.Bsp
import ..Mesh
import ..Tree

import Base:read,sizeof

const Entities       = 1
const Planes         = 2
const Vertices       = 3
const Visibility     = 4
const Nodes          = 5
const TexInfos       = 6
const Faces          = 7
const Lightmaps      = 8
const Leaves         = 9
const LeafFaceTable  = 10
const LeafBrushTable = 11
const Edges          = 12
const FaceEdgeTable  = 13
const Models         = 14
const Brushes        = 15
const BrushSides     = 16
const Pop            = 17
const Areas          = 18
const AreaPortals    = 19

const Magic = "IBSP"

immutable Header
	signature::Uint32
	version::Uint32
end

immutable Lump
	offset::Uint32
	length::Uint32
end

immutable Plane
	normal::Vector3{Float32}
	distance::Float32
	type_::Uint32
end

typealias Vertex Vector3{Float32}

immutable Face
	plane::Uint16
	plane_side::Uint16
	first_edge::Uint32
	num_edges::Uint16
	tex_info::Uint16
	lightmap_styles::Uint32
	lightmap_offset::Uint32
end

typealias FaceEdge Int32

immutable Node
	plane::Uint32
	front::Int32
	back::Int32
	bbox_min::Vector3{Int16}
	bbox_max::Vector3{Int16}
	first_face::Uint16
	num_faces::Uint16
end

immutable Leaf
	brush_or::Uint32
	cluster::Uint16
	area::Uint16
	bbox_min::Vector3{Int16}
	bbox_max::Vector3{Int16}
	first_face::Uint16
	num_faces::Uint16
	first_brush::Uint16
	num_brushes::Uint16
end

typealias LeafFace Uint16

immutable Edge
	v1::Uint16
	v2::Uint16
end

immutable VisOffsets
	sight::Uint32
	sound::Uint32
end

immutable TexInfo
	u::Vector4{Float32}
	v::Vector4{Float32}
	flags::Uint32
	value::Uint32
	name::ASCIIString
	next::Uint32
end
function read(io, ::Type{TexInfo})
	u = read(io, Vector4{Float32})
	v = read(io, Vector4{Float32})
	flags = read(io, Uint32)
	value = read(io, Uint32)
	name = lowercase(rstrip(bytestring(read(io, Uint8, 32)), "\0"))
	next = read(io, Uint32)
	return TexInfo(u,v,flags,value,name,next)
end
sizeof(::Type{TexInfo}) = 76

read(s::IO, t::DataType) = t(ntuple(length(t.types), x->read(s, t.types[x]))...)

function build(idx::Integer, nodes::Vector{Node}, planes::Vector{Plane}, leaves::Vector{Tree.Leaf})
	node = nodes[idx]
	if node.front < 0
		front = leaves[-node.front]
	else
		front = build(node.front+1, nodes, planes, leaves)
	end
	if node.back < 0
		back = leaves[-node.back]
	else
		back = build(node.back+1, nodes, planes, leaves)
	end
	normal = planes[node.plane+1].normal
	distance = planes[node.plane+1].distance
	return Tree.Node(front, back, normal, distance)
end

function readlump(io::IO, lump::Lump, t::Type)
	seek(io, lump.offset)
	read(io, t, int(lump.length / sizeof(t)))
end

function readheightmap(name::String)
	img = imread(string("/home/jay/q2renew/textures/", name, ".height.png"))
	width, height = size(img)[2:3]
	gray = uint8(convert(Array, img))[:,:,1]
	handle = GL.GenTexture()

	# create normals from heights
	normals = Array(Uint8, 0)
	for w = 1:width
		wp = w-1>=1 ? w-1 : width
		wn = w+1<=width ? w+1 : 1
		for h = 1:height
			hp = h-1>=1 ? h-1 : height
			hn = h+1<=height ? h+1 : 1
			x = int(gray[w,hp] - gray[w,hn]) / 32.0
			y = int(gray[wp,h] - gray[wn,h]) / 32.0
			n = Vector3(x,y,1.0)
			n /= norm(n)
			n = uint8(n * (255/2) + (255/2))
			push!(normals, n[1])
			push!(normals, n[2])
			push!(normals, n[3])
			push!(normals, gray[w,h])
		end
	end

	# upload image data to GPU
	GL.BindTexture(GL.TEXTURE_2D, handle)
	GL.TexImage2D(GL.TEXTURE_2D, GL.RGBA, width, height, GL.RGBA, normals)
	GL.GenerateMipmap(GL.TEXTURE_2D)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MIN_FILTER, GL.LINEAR_MIPMAP_LINEAR)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MAG_FILTER, GL.LINEAR)
	GL.TexParameteri(GL.TEXTURE_2D, 0x84FE, 8) # anisotropic filtering
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_S, GL.REPEAT)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_T, GL.REPEAT)
	GL.BindTexture(GL.TEXTURE_2D, 0)
	return handle
end

function read(io::IO, ::Type{Bsp})

	###   Read in header and lump data   #######################################

	hdr   = read(io, Header)
	lumps = read(io, Lump, 19)

	bin_edges    = readlump(io, lumps[Edges],      Edge)
	bin_faces    = readlump(io, lumps[Faces],      Face)
	bin_leaves   = readlump(io, lumps[Leaves],     Leaf)
	bin_nodes    = readlump(io, lumps[Nodes],      Node)
	bin_planes   = readlump(io, lumps[Planes],     Plane)
	bin_texinfos = readlump(io, lumps[TexInfos],   TexInfo)
	bin_vertices = readlump(io, lumps[Vertices],   Vertex)
	bin_vis      = readlump(io, lumps[Visibility], Uint8)

	face2edge = readlump(io, lumps[FaceEdgeTable], FaceEdge)
	leaf2face = readlump(io, lumps[LeafFaceTable], LeafFace)

	seek(io, lumps[Visibility].offset)
	num_clusters = read(io, Uint32)
	vis_offsets = read(io, VisOffsets, num_clusters)

	###   Convert File.TexInfo to Mesh.Textures   ##############################

	default_normal_map = GL.GenTexture()
	GL.BindTexture(GL.TEXTURE_2D, default_normal_map)
	GL.TexImage2D(GL.TEXTURE_2D, GL.RGBA, 1, 1, GL.RGBA, Uint8[127,127,255,127])
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MIN_FILTER, GL.NEAREST)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MAG_FILTER, GL.NEAREST)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_S, GL.REPEAT)
	GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_T, GL.REPEAT)
	GL.BindTexture(GL.TEXTURE_2D, 0)

	textures = Dict{String,Mesh.Texture}()
	for texinfo = bin_texinfos
		if !has(textures, texinfo.name)
			if texinfo.flags > 0
				println(texinfo.flags, ", ", texinfo.value)
			end

			f = qopen(string("textures/", texinfo.name, ".wal"))
			img = imread(f, Textures.WAL)
			close(f)

			width = uint32(size(img)[2])
			height = uint32(size(img)[3])
			diffuse = GL.GenTexture()

			# upload image data to GPU
			GL.BindTexture(GL.TEXTURE_2D, diffuse)
			GL.TexImage2D(GL.TEXTURE_2D, GL.RGB, width, height, GL.RGB, uint8(img[:]))
			GL.GenerateMipmap(GL.TEXTURE_2D)
			GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MIN_FILTER, GL.LINEAR_MIPMAP_LINEAR)
			GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_MAG_FILTER, GL.LINEAR)
			GL.TexParameteri(GL.TEXTURE_2D, 0x84FE, 8) # anisotropic filtering
			GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_S, GL.REPEAT)
			GL.TexParameteri(GL.TEXTURE_2D, GL.TEXTURE_WRAP_T, GL.REPEAT)
			GL.BindTexture(GL.TEXTURE_2D, 0)

			normal = default_normal_map
			try
				normal = readheightmap(texinfo.name)
				print("[h] ")
			catch e
				print("[ ] ")
			end

			textures[texinfo.name] = Mesh.Texture(diffuse, normal, width, height, texinfo.flags)
			println(texinfo.name)
		end
	end

	###   Convert File.Faces to Mesh.Faces   ###################################

	faces = Array(Mesh.Face,0)
	for face = bin_faces

		indices = Array(Uint16,0)
		first = face.first_edge + 1
		hub = bin_edges[abs(face2edge[first])+1].v1
		last = first + face.num_edges - 1
		for idx = face2edge[first:last]
			v1 = bin_edges[abs(idx)+1].v1
			v2 = bin_edges[abs(idx)+1].v2
			if idx < 0
				v1, v2 = v2, v1
			end
			if v1 != hub && v2 != hub
				push!(indices, hub)
				push!(indices, v2)
				push!(indices, v1)
			end
		end

		ibo = GL.GenBuffer()
		GL.BindBuffer(GL.ELEMENT_ARRAY_BUFFER, ibo)
		GL.BufferData(GL.ELEMENT_ARRAY_BUFFER, indices, GL.STATIC_DRAW)
		GL.BindBuffer(GL.ELEMENT_ARRAY_BUFFER, 0)

		normal = bin_planes[face.plane+1].normal
		if face.plane_side != 0
			normal = -normal
		end

		texinfo = bin_texinfos[face.tex_info+1]
		u = texinfo.u
		v = texinfo.v
		tex = textures[texinfo.name]

		push!(faces, Mesh.Face(indices, ibo, normal, tex, u, v, Array(Mesh.Light,0)))
	end

	###   Convert File.Leaves to Tree.Leaves   #################################

	leaves_in_cluster = Dict{Uint16,Vector{Tree.Leaf}}()

	leaves = Array(Tree.Leaf,0)
	for leaf = bin_leaves
		leaf_faces = Array(Mesh.Face,0)
		first = leaf.first_face + 1
		last = first + leaf.num_faces - 1
		for i = leaf2face[first:last]
			face = faces[i+1]
			if contains(leaf_faces, face)
				warn("duplicate face in leaf")
			else
				push!(leaf_faces, face)
			end
		end
		c = leaf.cluster
		if !has(leaves_in_cluster, c)
			leaves_in_cluster[c] = Array(Tree.Leaf,0)
		end
		leaf = Tree.Leaf(leaf_faces)
		push!(leaves, leaf)
		push!(leaves_in_cluster[c], leaf)
	end

	###   Populate Tree.Leaves visibility info   ###############################

	faces_in_cluster = Dict{Uint16,Vector{Mesh.Face}}()
	for tup = leaves_in_cluster
		c = tup[1]
		cluster_leaves = tup[2]
		visible_faces = Array(Mesh.Face,0)
		for leaf = cluster_leaves
			for face = leaf.faces
				push!(visible_faces, face)
			end
		end
		faces_in_cluster[c] = unique(visible_faces)
	end

	faces_from_cluster = Dict{Uint16,Vector{Mesh.Face}}()
	for tup = faces_in_cluster
		c = tup[1]                   # this cluster index
		visible_faces = copy(tup[2])

		if c == uint16(-1)
			continue
		end

		v = vis_offsets[c+1].sight+1 # pvs buffer index
		k = 0                        # other cluster index
		while k < num_clusters
			if bin_vis[v] == 0
				v += 1
				k += 8 * bin_vis[v]
			else
				for bit = 0:7
					if bin_vis[v] & (1 << bit) != 0 && has(faces_in_cluster, k)
						for face = faces_in_cluster[k]
							push!(visible_faces, face)
						end
					end
					k += 1
				end
			end
			v += 1
		end
		faces_from_cluster[c] = unique(visible_faces)
	end

	for i = 1:length(leaves)
		leaf = leaves[i]
		c = bin_leaves[i].cluster
		if c != uint16(-1)
			leaf.faces = faces_from_cluster[c]
		else
			leaf.faces = faces
		end
	end

	###   Build BSP tree   #####################################################

	tree = build(1, bin_nodes, bin_planes, leaves)

	###   Parse entity data   ##################################################

	# Read entity lump
	seek(io, lumps[Entities].offset)
	entityLump = bytestring(read(io, Uint8, lumps[Entities].length-1))
	entities = Array(Dict{String, String},0)

	# Create a dictionary for each entity
	for entity = split(entityLump, ['{', '}'])
		entity = strip(entity)
		if length(entity) < 1 # TODO: remove this w/ regex
			continue
		end
		dict = Dict{String, String}()
		for field = split(entity, '\n')
			fieldStr = split(strip(field), "\" \"")
			name = lstrip(strip(fieldStr[1], '"'), "_")
			value = strip(fieldStr[2], '"')
			dict[name] = value
		end
		if length(dict) < 1
			continue
		elseif !has(dict, "classname")
			warn("entity has no classname")
			continue
		else
			push!(entities, dict)
		end
	end

	###   Build list of lights from entity data   ##############################

	lights = Array(Mesh.Light,0)
	for ent = entities
		if ent["classname"] != "light"
			continue
		end
		if has(ent, "origin")
			origin = split(ent["origin"])
			origin = Vector3{Float32}(float32(origin[1]), float32(origin[2]), float32(origin[3]))
		else
			warn("light has no origin")
			continue
		end
		if has(ent, "color")
			color = split(ent["color"])
			color = Vector3{Float32}(float32(color[1]), float32(color[2]), float32(color[3]))
			color /= max(color)
		else
			color = Vector3{Float32}(1.0, 1.0, 1.0)
		end
		if has(ent, "light")
			power = float32(ent["light"])
		else
			power = float32(300)
		end
		power *= 1.5
		push!(lights, Mesh.Light(origin, color, power))
	end

	###   Associate lights with visible faces   ################################

	for light = lights
		for face = search(tree, light.origin).faces
			lit = Mesh.islit(face, bin_vertices, light)
			if lit && !contains(face.lights, light)
				push!(face.lights, light)
			end
		end
	end

	###   Calculate leaf statistics   ##########################################

	leaf_stats = Array(Int,0)
	for leaf = leaves
		if length(leaf.faces) > 0
			push!(leaf_stats, length(leaf.faces))
		end
	end
	println("visible faces per non-empty leaf:")
	println("min:    ", min(leaf_stats))
	println("mean:   ", mean(leaf_stats))
	println("median: ", median(leaf_stats))
	println("max:    ", max(leaf_stats))

	###   Calculate light statistics   #########################################

	light_stats = Array(Int,0)
	for face = faces
		push!(light_stats, length(face.lights))
	end
	println("lights per face:")
	println("min:    ", min(light_stats))
	println("mean:   ", mean(light_stats))
	println("median: ", median(light_stats))
	println("max:    ", max(light_stats))

	return Bsp(tree, entities, bin_vertices, max(light_stats))
end

end


text:
include("utilities.jl")

using SEAL
using Printf

"""
    example_bfv_basics()

Perform some basic operations such as encryption/decryption, multiplication, addition etc. using the
BFV scheme. This routine is based on the file `native/examples/1_bfv_basics.cpp` of the original
SEAL library and should yield the exact same output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/1_bfv_basics.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/1_bfv_basics.cpp)

See also: [`example_ckks_basics`](@ref)
"""
function example_bfv_basics()
  print_example_banner("Example: BFV Basics")

  parms = EncryptionParameters(SchemeType.bfv)

  poly_modulus_degree = 4096
  set_poly_modulus_degree!(parms, poly_modulus_degree)

  set_coeff_modulus!(parms, coeff_modulus_bfv_default(poly_modulus_degree))

  set_plain_modulus!(parms, 1024)

  context = SEALContext(parms)

  print_line(@__LINE__)
  println("Set encryption parameters and print")
  print_parameters(context)

  println("Parameter validation (success): ", parameter_error_message(context))

  println()
  println("~~~~~~ A naive way to calculate 4(x^2+1)(x+1)^2. ~~~~~~")

  keygen = KeyGenerator(context)
  secret_key_ = secret_key(keygen)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)

  encryptor = Encryptor(context, public_key_)

  evaluator = Evaluator(context)

  decryptor = Decryptor(context, secret_key_)

  print_line(@__LINE__)
  x = 6
  x_plain = Plaintext(string(x))
  println("Express x = " * string(x) * " as a plaintext polynomial 0x" * to_string(x_plain) * ".")

  print_line(@__LINE__)
  x_encrypted = Ciphertext()
  println("Encrypt x_plain to x_encrypted.")
  encrypt!(x_encrypted, x_plain, encryptor)

  println("    + size of freshly encrypted x: ", length(x_encrypted))

  println("    + noise budget in freshly encrypted x: ",
          invariant_noise_budget(x_encrypted, decryptor),
          " bits")

  x_decrypted = Plaintext()
  decrypt!(x_decrypted, x_encrypted, decryptor)
  println("    + decryption of x_encrypted: 0x", to_string(x_decrypted), " ...... Correct.")

  print_line(@__LINE__)
  println("Compute x_sq_plus_one (x^2+1).")
  x_sq_plus_one = Ciphertext()
  square!(x_sq_plus_one, x_encrypted, evaluator)
  plain_one = Plaintext("1")
  add_plain_inplace!(x_sq_plus_one, plain_one, evaluator)

  println("    + size of x_sq_plus_one: ", length(x_sq_plus_one))
  println("    + noise budget in x_sq_plus_one: ",
          invariant_noise_budget(x_sq_plus_one, decryptor),
          " bits")

  decrypted_result = Plaintext()
  decrypt!(decrypted_result, x_sq_plus_one, decryptor)
  println("    + decryption of x_sq_plus_one: 0x", to_string(decrypted_result), " ...... Correct.")

  print_line(@__LINE__)
  println("Compute x_plus_one_sq ((x+1)^2).")
  x_plus_one_sq = Ciphertext()
  add_plain!(x_plus_one_sq, x_encrypted, plain_one, evaluator)
  square_inplace!(x_plus_one_sq, evaluator)
  println("    + size of x_plus_one_sq: ", length(x_plus_one_sq))
  println("    + noise budget in x_plus_one_sq: ",
          invariant_noise_budget(x_plus_one_sq, decryptor),
          " bits")
  decrypt!(decrypted_result, x_plus_one_sq, decryptor)
  println("    + decryption of x_plus_one_sq: 0x", to_string(decrypted_result),
          " ...... Correct.")

  print_line(@__LINE__)
  println("Compute encrypted_result (4(x^2+1)(x+1)^2).")
  encrypted_result = Ciphertext()
  plain_four = Plaintext("4")
  multiply_plain_inplace!(x_sq_plus_one, plain_four, evaluator)
  multiply!(encrypted_result, x_sq_plus_one, x_plus_one_sq, evaluator)
  println("    + size of encrypted_result: ", length(encrypted_result))
  println("    + noise budget in encrypted_result: ",
          invariant_noise_budget(encrypted_result, decryptor),
          " bits")
  println("NOTE: Decryption can be incorrect if noise budget is zero.")

  println()
  println("~~~~~~ A better way to calculate 4(x^2+1)(x+1)^2. ~~~~~~")

  print_line(@__LINE__)
  println("Generate relinearization keys.")
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)

  print_line(@__LINE__)
  println("Compute and relinearize x_squared (x^2),")
  println(" "^13, "then compute x_sq_plus_one (x^2+1)")
  x_squared = Ciphertext()
  square!(x_squared, x_encrypted, evaluator)
  println("    + size of x_squared: ", length(x_squared))
  relinearize_inplace!(x_squared, relin_keys_, evaluator)
  println("    + size of x_squared (after relinearization): ", length(x_squared))
  add_plain!(x_sq_plus_one, x_squared, plain_one, evaluator)
  println("    + noise budget in x_sq_plus_one: ",
          invariant_noise_budget(x_sq_plus_one, decryptor),
          " bits")
  decrypt!(decrypted_result, x_sq_plus_one, decryptor)
  println("    + decryption of x_sq_plus_one: 0x", to_string(decrypted_result), " ...... Correct.")

  print_line(@__LINE__)
  x_plus_one = Ciphertext()
  println("Compute x_plus_one (x+1),")
  println(" "^13, "then compute and relinearize x_plus_one_sq ((x+1)^2).")
  add_plain!(x_plus_one, x_encrypted, plain_one, evaluator)
  square!(x_plus_one_sq, x_plus_one, evaluator)
  println("    + size of x_plus_one_sq: ", length(x_plus_one_sq))
  relinearize_inplace!(x_plus_one_sq, relin_keys_, evaluator)
  println("    + noise budget in x_plus_one_sq: ",
          invariant_noise_budget(x_plus_one_sq, decryptor),
          " bits")
  decrypt!(decrypted_result, x_plus_one_sq, decryptor)
  println("    + decryption of x_plus_one_sq: 0x", to_string(decrypted_result),
          " ...... Correct.")

  print_line(@__LINE__)
  println("Compute and relinearize encrypted_result (4(x^2+1)(x+1)^2).")
  multiply_plain_inplace!(x_sq_plus_one, plain_four, evaluator)
  multiply!(encrypted_result, x_sq_plus_one, x_plus_one_sq, evaluator)
  println("    + size of encrypted_result: ", length(encrypted_result))
  relinearize_inplace!(encrypted_result, relin_keys_, evaluator)
  println("    + size of encrypted_result (after relinearization): ", length(encrypted_result))
  println("    + noise budget in encrypted_result: ",
          invariant_noise_budget(encrypted_result, decryptor),
          " bits")
  println()
  println("NOTE: Notice the increase in remaining noise budget.")

  print_line(@__LINE__)
  println("Decrypt encrypted_result (4(x^2+1)(x+1)^2).")
  decrypt!(decrypted_result, encrypted_result, decryptor)
  println("    + decryption of 4(x^2+1)(x+1)^2 = 0x", to_string(decrypted_result),
          " ...... Correct.")
  println()

  print_line(@__LINE__)
  println("An example of invalid parameters")
  set_poly_modulus_degree!(parms, 2048)
  context = SEALContext(parms)
  print_parameters(context)
  println("Parameter validation (failed): ", parameter_error_message(context))
  println()

  return
end

ext: jl


text:
module UnsplashImages

using Pkg.Artifacts
using FileIO
using ImageShow

export unsplash_image

const artifacts_toml = normpath(@__DIR__, "..", "Artifacts.toml")

"""
    unsplash_image(name; ops...)
"""
function unsplash_image(name; ops...)
    imagefile = image_path(name)
    load(imagefile; ops...)
end

function image_path(filename)
    ensure_artifact_installed("images", artifacts_toml)
    image_dir = artifact_path(artifact_hash("images", artifacts_toml))
    normpath(image_dir, filename)
end

end # module UnsplashImages


text:
using PyPlot
using PyCall
@pyimport numpy as np

y = -0.5:0.1:3
#y = -0.5:0.1:3
x = -1.5:0.1:2
(X,Y) = np.meshgrid(x,y)
xv = vec(X)
yv = vec(Y)
g = [xv yv]
rosen(x,y) = (1-x)^2 + 100(y -x^2)^2
z = mapslices(p -> rosen(p[1], p[2]), g, 2)

fig = figure()
ax = fig[:add_subplot](111, projection="3d")
ax[:plot_surface](X, Y, reshape(z, length(x), length(y)), cmap="Spectral")
#ax[:plot_surface](X, Y, reshape(z, length(x), length(y)), cmap="terrain")

#ax[:clear]()


text:
# This file was generated, do not modify it. # hide
mtm = machine(tm, Xhp, y)
fit!(mtm)
rep = report(mtm)
@show round.(rep.measurements.^2, digits=2)
@show argmin(rep.measurements)


text:
##################################################
# Parsing of the command-line arguments
##################################################
tabargs = ArgParseSettings()
@add_arg_table! tabargs begin
    "--parallel"
    help = "Parallel computation: yes/no"
    arg_type = String
    default = "no"
    "--lmax"
    help = "Maximum harmonic number"
    arg_type = Int64
    default = 10
    "--nbK"
    help = "Number of sampling points for Klnnp"
    arg_type = Int64
    default = 100
    "--nbResPoints"
    help = "Number of points on the resonance lines"
    arg_type = Int64
    default = 100
    "--mBH"
    help = "Mass of the central BH (in Msun)"
    arg_type = Float64
    default = 4300000.0
    "--rh"
    help = "Influence radius of the BH (in mpc)"
    arg_type = Float64
    default = 2000.0
    "--inputBath"
    help = "Name of the input file for the bath"
    arg_type = String
    default = "../background/inputBath.txt"
    "--inputCluster"
    help = "Name of the input file for the S-cluster"
    arg_type = String
    default = "../background/inputCluster.txt"
end
parsed_args = parse_args(tabargs)


text:
function AbstractNFFTs.deconvolve!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if !p.params.storeDeconvolutionIdx
    deconvolve_alloc_free!(p, f, g)
  else
    p.tmpVecHat[:] .= vec(f) .* p.windowHatInvLUT[1]
    g[p.deconvolveIdx] = p.tmpVecHat
  end
  return
end

function deconvolve_alloc_free!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}) where {D,T,U}
  if D == 1
    _deconvolve_alloc_free!(p, f, g, 1)
  else
    @cthreads for o = 1:p.N[end]
        _deconvolve_alloc_free!(p, f, g, o)
    end
  end
end

@generated function _deconvolve_alloc_free!(p::NFFTPlan{T,D,1}, f::AbstractArray{U,D}, g::StridedArray{Complex{T},D}, o) where {D,T,U}
  quote
    @nexprs 1 d -> gidx_{$D-1} = rem(o+p.n[$D] - p.N[$D]÷2 - 1, p.n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:size(f,d+1)) d-> begin
        gidx_d = rem(l_d+p.n[d+1] - p.N[d+1]÷2 - 1, p.n[d+1]) + 1
      end begin
      N2 = p.N[1]÷2
      @inbounds @simd for i = 1:N2
        v = f[i, CartesianIndex(@ntuple $(D-1) l)]
        v *= p.windowHatInvLUT[1][i]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        g[i-N2+p.n[1], CartesianIndex(@ntuple $(D-1) gidx)] = v

        v = f[i+N2, CartesianIndex(@ntuple $(D-1) l)]
        v *= p.windowHatInvLUT[1][i+N2]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        g[i, CartesianIndex(@ntuple $(D-1) gidx)] = v
      end
    end
  end
end

########## deconvolve adjoint ##########

function AbstractNFFTs.deconvolve_transpose!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,T,U}
  if !p.params.storeDeconvolutionIdx
    deconvolve_transpose_alloc_free!(p, g, f)
  else
    p.tmpVecHat[:] = g[p.deconvolveIdx]
    f[:] .= vec(p.tmpVecHat) .* p.windowHatInvLUT[1]
  end
  return
end

function deconvolve_transpose_alloc_free!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}) where {D,T,U}
  if D == 1
    _deconvolve_transpose_alloc_free!(p, g, f, 1)
  else
    @cthreads for o = 1:p.N[end]
      _deconvolve_transpose_alloc_free!(p, g, f, o)
    end
  end
end

@generated function _deconvolve_transpose_alloc_free!(p::NFFTPlan{T,D,1}, g::AbstractArray{Complex{T},D}, f::StridedArray{U,D}, o) where {D,T,U}
  quote
    @nexprs 1 d -> gidx_{$D-1} = rem(o+p.n[$D] - p.N[$D]÷2 - 1, p.n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:size(f,d+1)) d-> begin
        gidx_d = rem(l_d+p.n[d+1] - p.N[d+1]÷2 - 1, p.n[d+1]) + 1
      end begin
      N2 = p.N[1]÷2
      @inbounds @simd for i = 1:N2
        v = g[i-N2+p.n[1], CartesianIndex(@ntuple $(D-1) gidx)]
        v *= p.windowHatInvLUT[1][i]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        f[i, CartesianIndex(@ntuple $(D-1) l)] = v

        v = g[i, CartesianIndex(@ntuple $(D-1) gidx)]
        v *= p.windowHatInvLUT[1][i+N2]
        @nexprs $(D-1) d -> v *= p.windowHatInvLUT[d+1][l_d]
        f[i+N2, CartesianIndex(@ntuple $(D-1) l)] = v
      end
    end
  end
end


text:
"""    with_seed(f, seed)

Perform the zero-argument function `f` after seeding the default RNG with `seed`.

Ensures that the RNG is re-seeded randomly when `f` exits normally or abnormally
(i.e. by throwing an exception).

Note the useful syntax:

    with_seed(1234) do
        ...
    end
"""
function with_seed(f, seed)
    Random.seed!(seed)
    try
        f()
    finally
        Random.seed!()
    end
end

function _names(k, d)
    string(k)
end
function _names(k, d::Vector)
    ["$(k)[$(i)]" for i in 1:length(d)]
end
function _names(k, d::Matrix)
    nx, ny = size(d)
    reshape(["$(k)[$(i),$(j)]" for i in 1:nx, j in 1:ny], nx*ny)
end

function _3Dify(d::Array{Float64, 2})
    ns, nc = size(d)
    reshape(d, (ns, 1, nc))
end
function _3Dify(d::Array{Array{Float64, 1}, 2})
    ns, nc = size(d)
    np, = size(d[1,1])

    out = zeros(ns, np, nc)
    for k in 1:nc
        for j in 1:np
            for i in 1:ns
                out[i,j,k] = d[i,k][j]
            end
        end
    end
    out
end
function _3Dify(d::Array{Matrix{Float64}, 2})
    ns, nc = size(d)
    nx, ny = size(d[1,1])

    out = zeros(ns, nx*ny, nc)
    for k in 1:nc
        for l in 1:ny
            for j in 1:nx
                for i in 1:ns
                    out[i, (l-1)*ny + j, k] = d[i,k][j,l]
                end
            end
        end
    end
    out
end

"""
    append_generated_quantities(trace, genq)

Given a trace (`MCMCChains` object) and a 2D array of named tuples corresponding
to `generated_quantities` output, concatenate the two together into a single
`MCMCChains` object.
"""
function append_generated_quantities(trace, genq)
    ks = keys(genq[1,1])
    nms = []
    ds = []
    for k in ks
        d = map(x -> getindex(x, k), genq)
        push!(nms, _names(k, d[1,1]))
        push!(ds, _3Dify(d))
    end
    nms = vcat(nms...)
    ds = cat(ds...; dims=2)

    nmap = trace.name_map
    atrace = cat(Array(trace, keys(nmap), append_chains=false)..., dims=3)
    Chains(cat(atrace, ds, dims=2), vcat(map(string, names(trace)), nms), nmap)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "tmux"
# Upstream uses version numbers like 3.1, 3.1a, 3.1b, 3.1c, we convert the
# letter into the patch number
version = v"3.3.0"
# 3.3-rc is a Release Candidate, but for stable ones you'll need to remote the -rc
tmux_tag = "$(version.major).$(version.minor)" * (version.patch > 0 ? Char('a' - 1 + version.patch) : "") * "-rc"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/tmux/tmux/releases/download/$(tmux_tag)/tmux-$(tmux_tag).tar.gz",
                  "c2f2314feab0d10868d15c5cb2eb6f4b0c30973fa98782b678e8ae94cf6ca268"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/tmux-*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-TERM=screen --enable-utf8proc
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows)

# The products that we will ensure are always built
products = [
    ExecutableProduct("tmux", :tmux)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("libevent_jll"),
    Dependency("Ncurses_jll"),
    Dependency("utf8proc_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
"""
    initialize_solver!(p,o)

Initialize the solver to the optimization [`Problem`](@ref) by initializing
the encapsulated `options` from within the [`RecordOptions`](@ref)` o`.
"""
function initialize_solver!(p::Problem, o::RecordOptions)
    initialize_solver!(p, o.options)
    get(o.recordDictionary, :Start, RecordGroup())(p, get_options(o), 0)
    return get(o.recordDictionary, :All, RecordGroup())(p, get_options(o), 0)
end
"""
    step_solver!(p,o,iter)

Do one iteration step (the `iter`th) for [`Problem`](@ref)` p` by modifying
the values in the [`Options`](@ref)` o.options` and record the result(s).
"""
function step_solver!(p::Problem, o::RecordOptions, i)
    step_solver!(p, o.options, i)
    get(o.recordDictionary, :Step, RecordGroup())(p, get_options(o), i)
    return get(o.recordDictionary, :All, RecordGroup())(p, get_options(o), i)
end
"""
    get_solver_result(o)

Return the final result after all iterations that is stored within the
(modified during the iterations) [`Options`](@ref)` o`.
"""
function get_solver_result(o::RecordOptions)
    return get_solver_result(o.options)
end
"""
    stop_solver!(p,o,i)

determine whether the solver for [`Problem`](@ref) `p` and the
[`RecordOptions`](@ref) `o` should stop at iteration `i`.
If so, do a (final) record to `:All` and `:Stop`.
"""
function stop_solver!(p::Problem, o::RecordOptions, i::Int)
    s = stop_solver!(p, o.options, i)
    if s
        get(o.recordDictionary, :Stop, RecordGroup())(p, get_options(o), typemin(Int))
        get(o.recordDictionary, :All, RecordGroup())(p, get_options(o), typemin(Int))
    end
    return s
end


text:
using Rheologies;
const R = Rheologies;

output_path = @__DIR__
## SPATIAL DOMAIN ##
dim = 2 # spatial dimensions

Lx = 20e3
Ly = 10e3
H_bdt = 1e3
n_seeds = 0
radius_bounds = (300, 400)
#el_geom = Triangle # Quadrilateral or Triangle, equivalent to Cell{dim,nnodes,nfaces}

nx = 50
ny = 50

corner1 = Vec(0.0, 0.0)
corner2 = Vec(Lx, Ly)
el_geom = Quadrilateral # Quadrilateral or Triangle, equivalent to Cell{dim,nnodes,nfaces}

# generate the grid and sets:
grid = JuAFEM.generate_grid(el_geom, (nx, ny), corner1, corner2)
# Add facesets / nodesets / cellsets on which bc will be applied (except for the top, bottom, left and right boundary that are implemented by default)
addnodeset!(
    grid,
    "clamped_lateral",
    x -> ((x[1] in (corner1[1], corner2[1])) & (x[2] == Ly / 2)),
);
addnodeset!(grid, "clamped_up", x -> ((x[1] == Lx/2) & (x[2] == Ly)));

## VARIABLES AND INTERPOLATIONS
variables = PrimitiveVariables{1}((:u,), (2,), el_geom)
#variables = PrimitiveVariables{2}((:u,:p), (2,1), el_geom) # not good !!!!

quad_order = 2
quad_type = :legendre # or :lobatto

## BOUNDARY CONDITIONS ##
# DIRICHLET :
# KEY : set on which to apply bc
# VALUES :
# 1) constrained variable
# 2) constraint function of (x,t)
# 3) constrained components
# NEUMANN :
# KEY : set on which to apply bc
# VALUES : traction as a function of (x,t)

cos_shape(x,A,Lx) = -A*cos(2π/Lx * x[1])
year = 3600 * 24 * 365.0
v_in_max = 0.01 / year # 1cm/year
cos_shape(x) = cos_shape(x,v_in_max,Lx)

bc = BoundaryConditions(
    dirichlet = Dict(
        "top" => [:u, (x,t)-> 0.0, 2],
        "bottom" => [:u, (x, t) -> cos_shape(x) * t, 2],
        "clamped_up" => [:u, (x, t) -> 0.0, 1]
        ),
    neumann = nothing,
)

bf = BodyForces([0.0, 0.0])#-9.81*2700])

## MATERIAL RHEOLOGY :

# define a rheology distribution function ...
# function multi_seeds(x,x0,radius,val_in,val_out)
#     for i in eachindex(radius)
#         (sqrt((x[1]-x0[i][1])^2 + (x[2]-x0[i][2])^2) <= radius[i]) && return val_in
#     end
#     return val_out
# end
# multi_seeds(x,val_in,val_out) = multi_seeds(x,x0,radius,val_in,val_out)
# and use it ...

elas = Elasticity(E = 70e9, ν = 0.3)

visco = Viscosity(η = x -> (x[2] >= Ly - H_bdt) ? 1e24 : 1e17)

VEP_rheology = Rheology(viscosity = visco, elasticity = elas)

### CLOCK ###
clock = Clock(tspan = (0.0, 50year), Δt = 0.1year, Δt_max = 20year)

### SOLVER ###
linsolver = BackslashSolver()
nlsolver = NewtonRaphson(max_iter_number = 10, atol = 1e-3, linear_solver = linsolver)

### MODEL ###
VEP_model = Model(
    grid = grid,
    variables = variables,
    quad_order = quad_order,
    quad_type = quad_type,
    bc = bc,
    body_forces = bf,
    rheology = VEP_rheology,
    initial_state = nothing,
    clock = clock,
    solver = nlsolver,
)

### OUTPUT ###
output_name = "viscous_example"
VEP_outputs = Dict(
    :σxx => (r, s) -> s.σ[1, 1],
    :σyy => (r, s) -> s.σ[2, 2],
    :σzz => (r, s) -> s.σ[3, 3],
    :η => (r, s) -> r.viscosity.η,
    :p => (r, s) -> -1 / 3 * tr(s.σ),
)

VEP_ow = VTKOutputWriter(
    VEP_model,
    joinpath(output_path, output_name),
    VEP_outputs,
    interval = 1,
    force_path = true,
)

### SOLVE ###
@time modelsol, u = solve(VEP_model, output_writer = VEP_ow, log = true);


text:
@testset "linear_models" begin

    @testset "trsbox" begin

        # Model specifications
        model = LinearModel(2)
        model.c[] = 0.0
        model.g .= [1.0, 2.0]
        model.kopt[] = 0

        # Box constraints
        l = [0.0, 0.0]
        u = [4.0, 4.0]

        # Trust-region
        Δ = 10.0

        # --------------------------------------------------------------------------------
        # Test 01 - Large trust-region and full active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [0.0, 0.0]
        model.fval[1] = model.c[] + dot( model.g, model.xopt )  # Interpolation model

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [0.0, 0.0] )
        @test( d == [0.0, 0.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 02 - Large trust-region and partial active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [2.0, 0.0]
        model.fval[1] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [0.0, 0.0] )
        @test( d == [-2.0, 0.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 03 - Large trust-region and partial active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [0.0, 4.0]
        model.fval[1] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [0.0, 0.0] )
        @test( d == [0.0, - 4.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 04 - Large trust-region and empty active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [2.0, 2.0]
        model.fval[1] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [0.0, 0.0] )
        @test( d == [- 2.0, - 2.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 05 - Large trust-region and empty active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [4.0, 4.0]
        model.fval[1] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [0.0, 0.0] )
        @test( d == [- 4.0, - 4.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

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

        # Model specifications
        model = LinearModel(2)
        model.c[] = 0.0
        model.g .= [1.0, 2.0]
        model.kopt[] = 1

        # Box constraints
        l = [- 2.0, - 2.0]
        u = [6.0, 6.0]

        # Trust-region
        Δ = 1.0

        # --------------------------------------------------------------------------------
        # Test 06 - Small trust-region and empty active set (full cauchy step)
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [0.0, 0.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( isapprox(x, (- Δ / norm(model.g)) * model.g) )
        @test( isapprox(d, (- Δ / norm(model.g)) * model.g) )
        @test( active_set == [false, false] )
        @test( status == :trust_region_boundary )

        # --------------------------------------------------------------------------------
        # Test 07 - Small trust-region and full active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [- 2.0, - 2.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 2.0, - 2.0] )
        @test( d == [0.0, 0.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 08 - Small trust-region and partial active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [- 2.0, 3.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 2.0, 2.0] )
        @test( d == [0.0, - 1.0] )
        @test( active_set == [true, false] )
        @test( status == :trust_region_boundary )

        # --------------------------------------------------------------------------------
        # Test 09 - Small trust-region and partial active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [2.0, -2.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [1.0, - 2.0] )
        @test( d == [- 1.0, 0.0] )
        @test( active_set == [false, true] )
        @test( status == :trust_region_boundary )

        # --------------------------------------------------------------------------------
        # Test 10 - Small trust-region and partial active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [6.0, - 1.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( isapprox(x, model.xopt + (- Δ / norm(model.g)) * model.g) )
        @test( isapprox(d, (- Δ / norm(model.g)) * model.g) )
        @test( active_set == [false, false] )
        @test( status == :trust_region_boundary )

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

        # Model specifications
        model = LinearModel(2)
        model.c[] = 0.0
        model.g .= [1.0, 2.0]
        model.kopt[] = 1

        # Box constraints
        l = [- 5.0, - 3.0]
        u = [6.0, 10.0]

        # Trust-region
        Δ = 4.0

        # --------------------------------------------------------------------------------
        # Test 11 - Intersection between the trust-region and
        #           the lower limits of the y-axis of the box.
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [0.0, 0.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- sqrt(7.0), - 3.0] )
        @test( d == [- sqrt(7.0), - 3.0] )
        @test( active_set == [false, true] )
        @test( status == :trust_region_boundary )

        # --------------------------------------------------------------------------------
        # Test 12 - Intersection between the trust-region and
        #           the lower limits of the x and y-axis of the box.
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [- 2.0, - 2.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 5.0, - 3.0] )
        @test( d == [- 3.0, - 1.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 13 - Intersection between the trust-region and
        #           the lower limits of the x-axis of the box.
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [-4.0, 5.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 5.0, 5.0 - sqrt(15)] )
        @test( d == [- 1.0, - sqrt(15)] )
        @test( active_set == [true, false] )
        @test( status == :trust_region_boundary )

        # --------------------------------------------------------------------------------
        # Test 14 - Intersection between the trust-region and
        #           the lower limits of the x and y-axis of the box (short step).
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [- 4.0, - 2.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 5.0, - 3.0] )
        @test( d == [- 1.0, - 1.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

        # --------------------------------------------------------------------------------
        # Test 15 - Full active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [- 5.0, - 3.0]
        model.fval[2] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [- 5.0, - 3.0] )
        @test( d == [ 0.0, 0.0] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

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

        # Model specifications
        model = LinearModel(2)
        model.c[] = 970.0
        model.g .= [-66.0, 692.0]
        model.kopt[] = 0

        # Box constraints
        l = [0.0, 0.0]
        u = [5.0, 5.0]

        # Trust-region
        Δ = 1.2

        # --------------------------------------------------------------------------------
        # Test 16 - Full active set
        # --------------------------------------------------------------------------------

        # Test specifications
        model.xopt .= [5.0, 0.0]
        model.fval[model.kopt[] + 1] = model.c[] + dot( model.g, model.xopt )

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v_aux = zeros(Float64, 2)

        status = LOWDER.trsbox!(model, Δ, l, u, active_set, x, d, v_aux)
        @test( x == [ 5.0, 0.0 ] )
        @test( d == [ 0.0, 0.0 ] )
        @test( active_set == [true, true] )
        @test( status == :full_active_set )

    end

    @testset "altmov" begin

        # --------------------------------------------------------------------------------
        # Test 01
        # --------------------------------------------------------------------------------

        # Model specifications
        model = LinearModel(2)
        model.c[] = 970.0
        model.g .= [1.0, 2.0]
        model.kopt[] = 0
        model.xbase .= zeros(Float64, 2)
        model.xopt .= zeros(Float64, 2)
        model.Y .= Matrix{Float64}(I, 2, 2)
        model.dst .= [1.0, 1.0]
        model.factorsY .= Matrix{Float64}(I, 2, 2)
        qrY = qr!(model.factorsY, Val(true))
        model.τY .= qrY.τ
        model.jpvtY .= qrY.jpvt

        # Box constraints
        l = [-5.0, -10.0]
        u = [20.0, 15.0]

        # Sample set radius
        δ = 2.0

        # Test specifications
        idx_t = 1

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v = zeros(Float64, 2)
        w = zeros(Float64, 2)

        status = LOWDER.altmov!(model, idx_t, δ, l, u, x, d, v, w, active_set)
        @test( x == [2.0, 0.0] )
        @test( d == [2.0, 0.0] )
        @test( status == :usual_altmov )

        # --------------------------------------------------------------------------------
        # Test 02
        # --------------------------------------------------------------------------------

        # Model specifications
        model = LinearModel(2)
        model.c[] = 0.0
        model.g .= [-66.0, 692.0]
        model.kopt[] = 0
        model.xbase .= [5.0, 0.0]
        model.xopt .= [5.0, 0.0]
        model.Y .= [-1.0 0.0; 0.0 1.0]
        model.dst .= [1.0, 1.0]
        model.factorsY .= [-1.0 0.0; 0.0 1.0]
        qrY = qr!(model.factorsY, Val(true))
        model.τY .= qrY.τ
        model.jpvtY .= qrY.jpvt

        # Box constraints
        l = [ 0.0, 0.0 ]
        u = [ 5.0, 5.0]

        # Sample set radius
        δ = 1.0

        # Test specifications
        idx_t = 1

        active_set = zeros(Bool, 2)
        x = zeros(Float64, 2)
        d = zeros(Float64, 2)
        v = zeros(Float64, 2)
        w = zeros(Float64, 2)

        status = LOWDER.altmov!(model, idx_t, δ, l, u, x, d, v, w, active_set)
        @test( x == [4.0, 0.0] )
        @test( d == [-1.0, 0.0] )
        @test( status == :usual_altmov )

    end

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
using BayesianSumProductNetworks
using Test
using SpecialFunctions

@testset "normal inverse gamma - mllh" begin

    x = [0.1, -0.1, 0.0]

    H = NormalInverseGamma(0, 1, 2, 3)
    sstat = NormalInvGammaSufficientStats(length(x), sum(x), sum( (x .- mean(x)).^2 ))

    n = length(x)
    (μ0, ν0, a0, b0) = params(H)
    (μn, νn, an, bn) = postparams(H, sstat)

    # Compute marginal likelihood according to Murphy.
    pD = (gamma(an) / gamma(a0)) * sqrt( (ν0/νn)) * ((2*b0)^a0 / (2*bn)^an) * (1/(π^(n/2)))
    logpD = log(pD)
    logpD2 = BayesianSumProductNetworks.invgmllh(μ0, ν0, a0, b0, μn, νn, an, bn, n)

    @test logpD ≈ logpD2
end

@testset "categorical - mllh" begin

    x = [0, 1, 0]

    γ = [0.5, 0.5]
    sstat = map(v -> sum(x .== v), 0:1)
    n = length(x)

    # Compute marginal likelihood according to Stephen Tu.
    pD = (gamma(sum(γ)) / prod(gamma.(γ))) * (prod(gamma.(map(j -> sum(γ[j] + sstat[j]), 1:2))) / gamma(n + sum(γ)))
    @test log(pD) ≈ BayesianSumProductNetworks.categoricalMLL(γ, sstat)
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
module Visualize_Test

using FactCheck, JNeuron

myimport=input(string(dirname(Base.source_path()),"/../examples/data/cell2.asc"));
blank_neuron=instantiate(myimport);
set_nsegs!(blank_neuron);

xy=JNeuron.plot_arrays(blank_neuron)

facts() do

    @fact length(xy) --> 198
    @fact size(xy[1],1) --> size(blank_neuron.secs[1].pt3d,1)
    @fact xy[1][1,1] --> blank_neuron.secs[1].pt3d[1].x
    @fact xy[1][1,2] --> blank_neuron.secs[1].pt3d[1].y
end

blank_neuron2=deepcopy(blank_neuron)

JNeuron.translate3d!(blank_neuron2,1.0,1.0,1.0)

facts() do

    @fact blank_neuron.secs[1].pt3d[1].x-->roughly(blank_neuron2.secs[1].pt3d[1].x-1.0)
    @fact blank_neuron.secs[1].pt3d[1].y-->roughly(blank_neuron2.secs[1].pt3d[1].y-1.0)
    @fact blank_neuron.secs[1].pt3d[1].z-->roughly(blank_neuron2.secs[1].pt3d[1].z-1.0)
end

blank_neuronx=deepcopy(blank_neuron)

JNeuron.rotate3d!(blank_neuronx,pi/2,1)

facts() do

    @fact blank_neuronx.secs[1].pt3d[1].x-->roughly(blank_neuron.secs[1].pt3d[1].x, atol=1e-2)
    @fact blank_neuronx.secs[1].pt3d[1].y-->roughly(-blank_neuron.secs[1].pt3d[1].z, atol=1e-2)
    @fact blank_neuronx.secs[1].pt3d[1].z-->roughly(blank_neuron.secs[1].pt3d[1].y, atol=1e-2)

end

blank_neurony=deepcopy(blank_neuron)

JNeuron.rotate3d!(blank_neurony,pi/2,2)

facts() do

    @fact blank_neurony.secs[1].pt3d[1].x-->roughly(blank_neuron.secs[1].pt3d[1].z, atol=1e-2)
    @fact blank_neurony.secs[1].pt3d[1].y-->roughly(blank_neuron.secs[1].pt3d[1].y, atol=1e-2)
    @fact blank_neurony.secs[1].pt3d[1].z-->roughly(-blank_neuron.secs[1].pt3d[1].x, atol=1e-2)

end

blank_neuronz=deepcopy(blank_neuron)

JNeuron.rotate3d!(blank_neuronz,pi/2,3)

facts() do

    @fact blank_neuronz.secs[1].pt3d[1].x-->roughly(-blank_neuron.secs[1].pt3d[1].y, atol=1e-2)
    @fact blank_neuronz.secs[1].pt3d[1].y-->roughly(blank_neuron.secs[1].pt3d[1].x, atol=1e-2)
    @fact blank_neuronz.secs[1].pt3d[1].z-->roughly(blank_neuron.secs[1].pt3d[1].z, atol=1e-2)

end

blank_neuronr=deepcopy(blank_neuron)

JNeuron.randomize_shape!(blank_neuronr)

facts() do

    @fact blank_neuron.secs[10].pt3d[1].x-->not(roughly(blank_neuronr.secs[1].pt3d[1].x))
    @fact blank_neuron.secs[10].pt3d[1].y-->not(roughly(blank_neuronr.secs[1].pt3d[1].y))
    @fact blank_neuron.secs[10].pt3d[1].z-->not(roughly(blank_neuronr.secs[1].pt3d[1].z))
end

end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
include("/Users/drutje/development/julia/Binance.jl/src/Binance.jl")
using Binance

using Pkg;

# packages to install, first time can take a while downloading, please be patient ...
packages=["HTTP","JSON","Dates","DataFrames","Plots","StatPlots","PlotThemes","GR","PyPlot","PyCall","LaTeXStrings","Plotly","PlotlyJS"]

for package in packages
    if get(Pkg.installed(),package,-1) == -1
        println(" getting package : ", package)
        Pkg.add(package)
    end
end

userdataChannel = Channel(10)

apiKey = ""
apiSecret = ""

listenKey = Binance.openUserData(apiKey)
#Binance.closeUserData(apiKey)
listenKey = Binance.pingUserData(apiKey,listenKey)

wsUserData


text:
# This file was generated, do not modify it. # hide
@from_network mach begin
    mutable struct MyAverageTwo
        regressor1=model1
        regressor2=model2
    end
end


text:
@reexport module BairdCounter

export BairdCounterEnv

using ReinforcementLearningBase

const ACTIONS = (:dashed, :solid)

Base.@kwdef mutable struct BairdCounterEnv <: AbstractEnv
    current::Int = rand(1:7)
end

RLBase.get_observation_space(env::BairdCounterEnv) = DiscreteSpace(7)
RLBase.get_action_space(env::BairdCounterEnv) = DiscreteSpace(length(ACTIONS))

function (env::BairdCounterEnv)(a)
    if ACTIONS[a] == :dashed
        env.current = rand(1:6)
    else
        env.current = 7
    end
    nothing
end

RLBase.observe(env::BairdCounterEnv) =
    (reward = 0.0, terminal = false, state = env.current)

function RLBase.reset!(env::BairdCounterEnv)
    env.current = rand(1:6)
    nothing
end

end


text:
function fcn()
    x = sin(1.)
    y = sin(1.)
    z = x + y
    return z
end

# run fallowing line from REPL
# Juno.@enter fcn()


text:
#get cumulative distribution (discrete) from scaled/unnormalised probability distribution (discrete)
function pdist2cdist(pdist::AbstractVector{T}) where {T}
    cdist = Vector{T}(undef,length(pdist))
    cdist[1] = pdist[1]
    @views for i ∈ 2:length(pdist)
        cdist[i] = cdist[i-1] + pdist[i]
    end
    K = 1/cdist[end]
    @turbo for i ∈ eachindex(cdist)
        cdist[i] *= K
    end
    return cdist
end

# Find index of closest element of data array to input x (same as argmin(abs.(data.-x)) but faster)
function findclosest(x::T, data::AbstractVector{T}) where T<:Real
    d = zero(T)
    mindx::Int = 1
    mmag::T = Inf
    @inbounds @views for i ∈ eachindex(data)
        d = abs(data[i] - x)
        if(d<mmag)
            mmag = d
            mindx = i
        end
    end
    return mindx
end


text:
mutable struct driving_control_cmd_t <: LCMType
    timestamp::Int64
    steering_angle::Float64
    throttle_value::Float64
    brake_value::Float64
end

@lcmtypesetup(driving_control_cmd_t)


text:
# # Step

struct Step
    n::Int
end

# constructor:
Step(; n=5) = Step(n)

@create_docs(Step,
             header="Step(; n=1)",
             example="Step(2)",
             body="Train for `n` more iterations. "*
             "Will never trigger a stop. ")

function update!(c::Step, model, verbosity, n, state=(new_iterations = 0,))
    new_iterations = state.new_iterations
    verbosity > 1 &&
        @info "Stepping model for $(c.n) more iterations. "
    train!(model, c.n)
    state  = (new_iterations = new_iterations + c.n,)
    return state
end

@inline function takedown(c::Step, verbosity, state)
    verbosity > 1 &&
        @info "A total of $(state.new_iterations) iterations added. "
    return state
end

# # Info

struct Info{F<:Function}
    f::F
end

# constructor:
Info(; f::Function=identity) = Info(f)

@create_docs(Info,
             header="Info(f=identity)",
             example="Info(my_loss_function)",
             body="Log to `Info` the value of `f(m)`, "*
             "where `m` "*
             "is the object being iterated. If "*
             "`IterativeControl.expose(m)` has been overloaded, then "*
             "log `f(expose(m))` instead.\n\n"*
             "Can be suppressed by setting the global verbosity level "*
             "sufficiently low. \n\n"*
             "See also [`Warn`](@ref), [`Error`](@ref). ")

function update!(c::Info, model, verbosity, args...)
    verbosity > 0 && @info _log_eval(c.f, model)
    return nothing
end

# # Warn

struct Warn{P<:Function,F<:Union{Function,String}}
    predicate::P
    f::F
end

# constructor:
Warn(predicate; f="") = Warn(predicate, f)

@create_docs(Warn,
             header="Warn(predicate; f=\"\")",
             example="Warn(m -> length(m.cache) > 100, "*
             "f=\"Memory low\")",
             body="If `predicate(m)` is `true`, then "*
             "log to `Warn` the value of `f` "*
             "(or `f(IterationControl.expose(m))` if `f` is a function). "*
             "Here `m` "*
             "is the object being iterated.\n\n"*
             "Can be suppressed by setting the global verbosity level "*
             "sufficiently low.\n\n"*
             "See also [`Info`](@ref), [`Error`](@ref). ")

function update!(c::Warn, model, verbosity, n, state=(warnings=(),))
    warnings = state.warnings
    if c.predicate(model)
        warning = _log_eval(c.f, model)
        verbosity > -1 && @warn warning
        newstate = (warnings = tuple(warnings..., warning),)
    else
        newstate = state
    end
    return newstate
end

function takedown(c::Warn, verbosity, state)
    warnings = join(state.warnings, "\n")
    verbosity > 1 && !isempty(warnings) &&
        @warn "A `Warn` control issued these warnings:\n$warnings"
    return state
end

# # Error

struct Error{P<:Function,F<:Union{Function,String}}
    predicate::P
    f::F
    exception::Union{Exception, Nothing}
end

# constructor:
Error(predicate; f="", exception=nothing) = Error(predicate, f, exception)

@create_docs(Error,
             header="Error(predicate; f=\"\", exception=nothing))",
             example="Error(m -> isnan(m.bias), f=\"Bias overflow!\")",
             body="If `predicate(m)` is `true`, then "*
             "log at the `Error` level the value of `f` "*
             "(or `f(IterationControl.expose(m))` if `f` is a function) "*
             "and stop iteration at the end of the current control cycle. "*
             "Here `m` "*
             "is the object being iterated.\n\n"*
             "Specify `exception=...` to throw an immediate "*
             "execption, without "*
             "waiting to the end of the control cycle.\n\n"*
             "See also [`Info`](@ref), [`Warn`](@ref). ")

function update!(c::Error,
                 model,
                 verbosity,
                 n=1,
                 state=(done = false, error=()))
    if c.predicate(model)
        error = _log_eval(c.f, model)
        @error error
        c.exception === nothing || throw(c.exception)
        state = (done = true, error=error)
    end
    return state
end

done(c::Error, state) = state.done

takedown(c::Error, verbosity, state) = state

# # CALLBACK

struct Callback{F<:Function}
    f::F
    stop_if_true::Bool
    stop_message::Union{String,Nothing}
    raw::Bool
end

# constructor:
Callback(f::Function;
         stop_if_true=false,
         stop_message=nothing,
         raw=false) = Callback(f, stop_if_true, stop_message, raw)
Callback(; f=identity, kwargs...) = Callback(f, kwargs...)

@create_docs(Callback,
             header="Callback(f=_->nothing, stop_if_true=false, "*
             "stop_message=nothing, raw=false)",
             example="Callback(m->put!(v, my_loss_function(m))",
             body="Call `f(IterationControl.expose(m))`, where "*
             "`m` is the object being iterated, unless `raw=true`, in "*
             "which case call `f(m)` (guaranteed if `expose` has not been "*
             "overloaded.) "*
             "If `stop_if_true` is `true`, then trigger an early stop "*
             "if the value returned by `f` is `true`, logging the "*
             "`stop_message` if specified. ")

function update!(c::Callback, model, verbosity, n, state=(done=false, ))
    r = c.f(expose(model, c.raw))
    done = (c.stop_if_true && r isa Bool && r) ? true : false
    return (done=done,)
end

done(c::Callback, state) = state.done

function takedown(c::Callback, verbosity, state)
    if state.done
        message = c.stop_message === nothing ?
            "Stop triggered by a `Callback` control. " :
            c.stop_message
        verbosity > 0 && @info message
        return (done = true, log = message)
    else
        return (done = false, log = "")
    end
end

# # DATA

struct Data{S}
    data::S
    stop_when_exhausted::Bool
end
Data(; data=(), stop_when_exhausted=false) = Data(data, stop_when_exhausted)
Data(data; kwargs...) = Data(data=data; kwargs...)

Base.show(io::IO, d::Data{S}) where S =
    print(io, "Data{$S}(<omitted data>; "*
          "stop_when_exhausted=$(d.stop_when_exhausted))")

@create_docs(Data,
             header="Data(my_data; stop_when_exhausted=false)",
             example="Data(rand(100))",
             body="In each application of this control a new `item` from the "*
             "iterable, `data`, is retrieved (using `iterate`) and "*
             "`IterationControl.ingest!(m, item)` is called. Here "*
             "`m` is the object being iterated. \n\n"*
             "A control becomes passive once the `data` iterable is done. "*
             "To trigger "*
             "a stop *after one passive application of the control*, set "*
             "`stop_when_exhausted=true`. ")

const DATA_EXHAUSTED = "Data exhausted. "
const DATA_STOP ="Stop triggered because data exhausted ."

function update!(c::Data, model, verbosity, n)
    iter = c.data
    data_exhausted = true
    next = iter_state = iterate(iter)
    if next !== nothing
        data_exhausted = false
        item, iter_state = next
    end
    data_exhausted && verbosity > 0 && !c.stop_when_exhausted &&
        @info DATA_EXHAUSTED
    data_exhausted || ingest!(model, item)
    done = data_exhausted && c.stop_when_exhausted
    return (iter_state = iter_state, done = done)
end

function update!(c::Data, model, verbosity, n, state)
    iter = c.data
    iter_state = state.iter_state
    data_exhausted = true
    if iter_state !== nothing
        next = iterate(iter, iter_state)
        if next === nothing
            iter_state = nothing
        else
            data_exhausted = false
            item, iter_state = next
        end
    end
    data_exhausted && verbosity > 0 && !c.stop_when_exhausted &&
        iter_state !== nothing && @info DATA_EXHAUSTED
    data_exhausted || ingest!(model, item)
    done = data_exhausted && c.stop_when_exhausted
    return (iter_state = iter_state, done = done)
end

done(c::Data, state) = state.done

function takedown(c::Data, verbosity, state)
    if state.done
        verbosity > 0 && @info DATA_STOP
        return (done = true, log = DATA_STOP)
    else
        return (done = false, log = "")
    end
end

# # WithLossDo

struct WithLossDo{F<:Function}
    f::F
    stop_if_true::Bool
    stop_message::Union{String,Nothing}
end

# constructor:
WithLossDo(f::Function;
     stop_if_true=false,
     stop_message=nothing) = WithLossDo(f, stop_if_true, stop_message)
WithLossDo(; f=x->@info("loss: $x"), kwargs...) = WithLossDo(f, kwargs...)

@create_docs(WithLossDo,
             header="WithLossDo(f=x->@info(\"loss: \$x\"), "*
             "stop_if_true=false, "*
             "stop_message=nothing)",
             example="WithLossDo(x->put!(my_losses, x))",
             body="Call `f(loss)`, where "*
             "`loss` is current loss.\n\n"*
             "If `stop_if_true` is `true`, then trigger an early stop "*
             "if the value returned by `f` is `true`, logging the "*
             "`stop_message` if specified. ")

needs_loss(::WithLossDo) = true

function update!(c::WithLossDo,
                 model,
                 verbosity,
                 n,
                 state=(loss=nothing, done=false))
    loss = IterationControl.loss(model)
    loss === nothing && throw(ERR_NEEDS_LOSS)
    r = c.f(loss)
    done = (c.stop_if_true && r isa Bool && r) ? true : false
    return (loss=loss, done=done)
end

done(c::WithLossDo, state) = state.done

function takedown(c::WithLossDo, verbosity, state)
    verbosity > 1 && @info "final loss: $(state.loss). "
    if state.done
        message = c.stop_message === nothing ?
            "Stop triggered by a `WithLossDo` control. " :
            c.stop_message
        verbosity > 0 && @info message
        return merge(state, (log = message,))
    else
        return merge(state, (log = "",))
    end
end

# # WithTrainingLossesDo

struct WithTrainingLossesDo{F<:Function}
    f::F
    stop_if_true::Bool
    stop_message::Union{String,Nothing}
end

# constructor:
WithTrainingLossesDo(f::Function;
     stop_if_true=false,
     stop_message=nothing) = WithTrainingLossesDo(f, stop_if_true, stop_message)
WithTrainingLossesDo(; f=v->@info("training: $v"), kwargs...) =
    WithTrainingLossesDo(f, kwargs...)

@create_docs(WithTrainingLossesDo,
             header="WithTrainingLossesDo(f=v->@info(\"training: \$v\"), "*
             "stop_if_true=false, "*
             "stop_message=nothing)",
             example="WithTrainingLossesDo(v->put!(my_losses, last(v))",
             body="Call `f(training_losses)`, where "*
             "`training_losses` is the vector of most recent batch "*
             "of training losses.\n\n"*
             "If `stop_if_true` is `true`, then trigger an early stop "*
             "if the value returned by `f` is `true`, logging the "*
             "`stop_message` if specified. ")

needs_training_losses(::WithTrainingLossesDo) = true

function update!(c::WithTrainingLossesDo,
                 model,
                 verbosity,
                 n,
                 state=(latest_training_loss = nothing, done = false))
    losses = IterationControl.training_losses(model)
    losses === nothing && throw(ERR_NEEDS_TRAINING_LOSSES)
    r = c.f(losses)
    done = (c.stop_if_true && r isa Bool && r) ? true : false
    return (latest_training_loss=losses[end], done=done)
end

done(c::WithTrainingLossesDo, state) = state.done

function takedown(c::WithTrainingLossesDo, verbosity, state)
    verbosity > 1 &&
        @info "final training loss: $(state.latest_training_loss). "
    if state.done
        message = c.stop_message === nothing ?
            "Stop triggered by a `WithTrainingLossesDo` control. " :
            c.stop_message
        verbosity > 0 && @info message
        return merge(state, (log = message,))
    else
        return merge(state, (log = "",))
    end
end

# # WithNumberDo

struct WithNumberDo{F<:Function}
    f::F
    stop_if_true::Bool
    stop_message::Union{String,Nothing}
end

# constructor:
WithNumberDo(f::Function;
     stop_if_true=false,
     stop_message=nothing) = WithNumberDo(f, stop_if_true, stop_message)
WithNumberDo(; f=n->@info("number: $n"), kwargs...) = WithNumberDo(f, kwargs...)

@create_docs(WithNumberDo,
             header="WithNumberDo(f=n->@info(\"number: \$n\"), "*
             "stop_if_true=false, "*
             "stop_message=nothing)",
             example="WithNumberDo(n->put!(my_channel, n))",
             body="Call `f(n + 1)`, where "*
             "`n` is the number of complete control cycles. "*
             "of the control (so, `n = 1, 2, 3, ...`, unless control is "*
             "wrapped in a [`IterationControl.skip`](@ref))`.\n\n"*
             "If `stop_if_true` is `true`, then trigger an early stop "*
             "if the value returned by `f` is `true`, logging the "*
             "`stop_message` if specified. ")

function update!(c::WithNumberDo,
                 model,
                 verbosity,
                 n,
                 state=(done = false, n = n))
    r = c.f(n)
    done = (c.stop_if_true && r isa Bool && r) ? true : false
    return (done = done, n = n)
end

done(c::WithNumberDo, state) = state.done

function takedown(c::WithNumberDo, verbosity, state)
    verbosity > 1 && @info "final number: $(state.n). "
    if state.done
        message = c.stop_message === nothing ?
            "Stop triggered by a `WithNumberDo` control. " :
            c.stop_message
        verbosity > 0 && @info message
        return merge(state, (log = message,))
    else
        return merge(state, (log = "",))
    end
end


text:
# This file was generated, do not modify it. # hide
f = @node folds(X, 3)


text:
mutable struct StepFunction{S,F,G} <: AbstractStep
    selections::S
    params::Union{Nothing, Vector{Symbol}}
    f::F
    broadcast::Bool
    trained::Bool
    skip::Bool
    prehook::G
end

function transform!(s::StepFunction, df)
    if s.broadcast
        for k in s.params
            df[k] = s.f.(df[k])
        end
    else
        for k in s.params
            df[k] = s.f(df[k])
        end
    end
    df
end
function transform(s::StepFunction, df)
    transform!(s, copy(df))
end
function fit!(s::StepFunction, df)
    s.params = getselectionkeys(df, s.selections)
end
function step_function!(r::Recipe, s...; f=nothing, broadcast::Bool=true, skip::Bool=false, prehook=identity)
    f === nothing && throw(error("No function supplied to `step_function!`"))
    push!(r.steps,
          StepFunction([s...],
                       nothing,
                       f,
                       broadcast,
                       true,
                       skip,
                       prehook))
end


text:
module Stats
import Luna: Maths, Grid, Modes, Utils, settings, PhysData, Fields, Processing
import Luna.PhysData: wlfreq, c, ε_0
import Luna.NonlinearRHS: TransModal, TransModeAvg
import Luna.Nonlinear: PlasmaCumtrapz
import Luna.Capillary: MarcatiliMode
import FFTW
import LinearAlgebra: mul!
import Printf: @sprintf
import Logging: @warn

"""
    ω0(grid)

Create stats function to calculate the centre of mass (first moment) of the spectral power
density.
"""
function ω0(grid)
    addstat! = let ω=grid.ω
        function addstat!(d, Eω, Et, z, dz)
            d["ω0"] = squeeze(Maths.moment(ω, abs2.(Eω); dim=1))
        end
    end
    return addstat!
end

squeeze(ω0::Array{T, 1}) where T = ω0[1]
squeeze(ω0::Array{T, 2}) where T = ω0[1, :]

"""
    energy(grid, energyfun_ω)

Create stats function to calculate the total energy.
"""
function energy(grid, energyfun_ω)
    function addstat!(d, Eω, Et, z, dz)
        if ndims(Eω) > 1
            d["energy"] = [energyfun_ω(Eω[:, i]) for i=1:size(Eω, 2)]
        else
            d["energy"] = energyfun_ω(Eω)
        end
    end
    return addstat!
end

"""
    energy_λ(grid, energyfun_ω, λlims; label)

Create stats function to calculate the energy in a wavelength region given by `λlims`.
If `label` is omitted, the stats dataset is named by the wavelength limits.
"""
function energy_λ(grid, energyfun_ω, λlims; label=nothing, winwidth=0)
    λlims = collect(λlims)
    ωmin, ωmax = extrema(wlfreq.(λlims))
    window = Maths.planck_taper(grid.ω, ωmin-winwidth, ωmin, ωmax, ωmax+winwidth)
    if isnothing(label)
        λnm = 1e9.*λlims
        label = @sprintf("%.2fnm_%.2fnm", minimum(λnm), maximum(λnm))
    end
    energy_window(grid, energyfun_ω, window; label=label)
end

"""
    energy_window(grid, energyfun_ω, window; label)

Create stats function to calculate the energy filtered by a `window`. The stats dataset will
be named `energy_[label]`.
"""
function energy_window(grid, energyfun_ω, window; label)
    key = "energy_$label"
    function addstat!(d, Eω, Et, z, dz)
        if ndims(Eω) > 1
            d[key] = [energyfun_ω(Eω[:, i].*window) for i=1:size(Eω, 2)]
        else
            d[key] = energyfun_ω(Eω.*window)
        end
    end
    return addstat!
end

"""
    peakpower(grid)

Create stats function to calculate the peak power.
"""
function peakpower(grid)
    function addstat!(d, Eω, Et, z, dz)
        if ndims(Et) > 1
            d["peakpower"] = dropdims(maximum(abs2.(Et), dims=1), dims=1)
        else
            d["peakpower"] = maximum(abs2.(Et))
        end
    end
    return addstat!
end

"""
    peakpower(grid, Eω, window; label)

Create stats function to calculate the peak power within a frequency range defined by the
window function `window`. `window` must have the same length as `grid.ω`. The stats
dataset is labeled as `peakpower_[label]`.
"""
function peakpower(grid, Eω, window::Vector{<:Real}; label)
    Etbuf, analytic! = plan_analytic(grid, Eω) # output buffer and function for inverse FT
    Eωbuf = similar(Eω) # buffer for Eω with window applied
    key = "peakpower_$label"
    function addstat!(d, Eω, Et, z, dz)
        Eωbuf .= Eω .* window
        analytic!(Etbuf, Eωbuf)
        if ndims(Etbuf) > 1
            d[key] = dropdims(maximum(abs2.(Etbuf), dims=1), dims=1)
        else
            d[key] = maximum(abs2.(Etbuf))
        end
    end
end

"""
    peakpower(grid, Eω, λlims; label=nothing)

Create stats function to calculate the peak power within a frequency range defined by the
wavelength limits `λlims`. If `label` is given, the stats dataset is labeled as
`peakpower_[label]`, otherwise `label` is created automatically from `λlims`.
"""
function peakpower(grid, Eω, λlims::NTuple{2, <:Real}; label=nothing, winwidth=:auto)
    window = Processing.ωwindow_λ(grid.ω, λlims; winwidth=winwidth)
    if isnothing(label)
        λnm = 1e9.*λlims
        label = @sprintf("%.2fnm_%.2fnm", minimum(λnm), maximum(λnm))
    end
    peakpower(grid, Eω, window; label=label)
end

"""
    peakintensity(grid, aeff)

Create stats function to calculate the mode-averaged peak intensity given the effective area
`aeff(z)`.
"""
function peakintensity(grid, aeff)
    function addstat!(d, Eω, Et, z, dz)
        d["peakintensity"] = maximum(abs2.(Et))/aeff(z)
    end
end

"""
    peakintensity(grid, mode)

Create stats function to calculate the peak intensity for several modes.
"""
function peakintensity(grid, modes::Modes.ModeCollection; components=:y)
    tospace = Modes.ToSpace(modes, components=components)
    npol = tospace.npol
    Et0 = zeros(ComplexF64, (length(grid.t), npol))
    function addstat!(d, Eω, Et, z, dz)
        Modes.to_space!(Et0, Et, (0, 0), tospace; z=z)
        if npol > 1
            d["peakintensity"] = c*ε_0/2 * maximum(sum(abs2.(Et0), dims=2))
        else
            d["peakintensity"] = c*ε_0/2 * maximum(abs2.(Et0))
        end
    end
end

"""
    fwhm_t(grid)

Create stats function to calculate the temporal FWHM (pulse duration) for mode average.
"""
function fwhm_t(grid)
    function addstat!(d, Eω, Et, z, dz)
        Pt = abs2.(Et)
        if ndims(Et) > 1
            Ptsum = dropdims(sum(Pt; dims=2); dims=2)
            d["fwhm_t_min"] = [Maths.fwhm(grid.t, Pt[:, i], method=:linear)
                              for i=1:size(Et, 2)]
            d["fwhm_t_max"] = [Maths.fwhm(grid.t, Pt[:, i], method=:linear, minmax=:max)
                              for i=1:size(Et, 2)]
            d["fwhm_t_min_allmodes"] = Maths.fwhm(grid.t, Ptsum, method=:linear)
            d["fwhm_t_max_allmodes"] = Maths.fwhm(grid.t, Ptsum, method=:linear, minmax=:max)
        else
            d["fwhm_t_min"] = Maths.fwhm(grid.t, Pt, method=:linear, minmax=:min)
            d["fwhm_t_max"] = Maths.fwhm(grid.t, Pt, method=:linear, minmax=:max)
        end
    end
end

"""
    fwhm_r(grid, modes; components=:y)

Create stats function to calculate the radial FWHM (aka beam size) in a modal propagation.
"""
function fwhm_r(grid, modes; components=:y)
    tospace = Modes.ToSpace(modes, components=components)
    npol = tospace.npol
    Eω0 = zeros(ComplexF64, (length(grid.ω), npol))

    function addstat!(d, Eω, Et, z, dz)
        function f(r)
            Modes.to_space!(Eω0, Eω, (r, 0), tospace; z=z)
            sum(abs2.(Eω0))
        end
        d["fwhm_r"] = 2*Maths.hwhm(f)
    end
end

"""
    electrondensity(grid, ionrate, dfun, aeff; oversampling=1)

Create stats function to calculate the maximum electron density in mode average.

If oversampling > 1, the field is oversampled before the calculation
!!! warning
    Oversampling can lead to a significant performance hit
"""
function electrondensity(grid::Grid.RealGrid, ionrate!, dfun, aeff; oversampling=1)
    to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling)
    δt = to[2] - to[1]
    # ionfrac! stores the time-dependent ionisation fraction in out and returns the max
    # ionisation rate
    function ionfrac!(out, Et)
        ionrate!(out, Et)
        ratemax = maximum(out)
        Maths.cumtrapz!(out, δt) # in-place cumulative integration
        @. out = 1 - exp(-out)
        return ratemax
    end
    frac = similar(to)
    function addstat!(d, Eω, Et, z, dz)
        # note: oversampling returns its arguments without any work done if factor==1
        to, Eto = Maths.oversample(grid.t, Et, factor=oversampling)
        @. Eto /= sqrt(ε_0*c*aeff(z)/2)
        ratemax = ionfrac!(frac, real(Eto))
        d["electrondensity"] = frac[end]*dfun(z)
        d["peak_ionisation_rate"] = ratemax
    end
end

"""
    electrondensity(grid, ionrate, dfun, modes; oversampling=1)

Create stats function to calculate the maximum electron density for multimode simulations.

If oversampling > 1, the field is oversampled before the calculation
!!! warning
    Oversampling can lead to a significant performance hit
"""
function electrondensity(grid::Grid.RealGrid, ionrate!, dfun,
                         modes::Modes.ModeCollection;
                         components=:y, oversampling=1) where N
    to, Eto = Maths.oversample(grid.t, complex(grid.t), factor=oversampling)
    δt = to[2] - to[1]
    # ionfrac! stores the time-dependent ionisation fraction in out and returns the max
    # ionisation rate
    function ionfrac!(out, Et)
        ionrate!(out, Et)
        ratemax = maximum(out)
        Maths.cumtrapz!(out, δt) # in-place cumulative integration
        @. out = 1 - exp(-out)
        return ratemax
    end
    tospace = Modes.ToSpace(modes, components=components)
    frac = similar(to)
    npol = tospace.npol
    Et0 = zeros(ComplexF64, (length(to), npol))
    function addstat!(d, Eω, Et, z, dz)
        # note: oversampling returns its arguments without any work done if factor==1
        to, Eto = Maths.oversample(grid.t, Et, factor=oversampling)
        Modes.to_space!(Et0, Eto, (0, 0), tospace; z=z)
        if npol > 1
            ratemax = ionfrac!(frac, hypot.(real(Et0[:, 1]), real(Et0[:, 2])))
        else
            ratemax = ionfrac!(frac, real(Et0[:, 1]))
        end
        d["electrondensity"] = frac[end]*dfun(z)
        d["peak_ionisation_rate"] = ratemax
    end
end

"""
    density(dfun)

Create stats function to capture the gas density as defined by `dfun(z)`
"""
function density(dfun)
    function addstat!(d, Eω, Et, z, dz)
        d["density"] = dfun(z)
    end
end

"""
    pressure(dfun, gas)

Create stats function to capture the pressure. Like [`density`](@ref) but converts to
pressure.
"""
function pressure(dfun, gas)
    function addstat!(d, Eω, Et, z, dz)
        d["pressure"] = PhysData.pressure(gas, dfun(z))
    end
end

function pressure(dfun, gases::Tuple)
    function addstat!(d, Eω, Et, z, dz)
        dens = dfun(z)
        for (di, gi) in zip(dens, gases)
            d["pressure_$gi"] = PhysData.pressure(gi, di)
        end
    end
end

"""
    core_radius(a)

Create stats function to capture core radius as defined by `a` (either a `Number` or a
callable `a(z)`)
"""
function core_radius(a::Number)
    function addstat!(d, Eω, Et, z, dz)
        d["core_radius"] = a
    end
end

function core_radius(afun)
    function addstat!(d, Eω, Et, z, dz)
        d["core_radius"] = afun(z)
    end
end

"""
    zdw(mode)

Create stats function to capture the zero-dispersion wavelength (ZDW).

!!! warning
    Since [`Modes.zdw`](@ref) is based on root-finding of a derivative, this can be slow!
"""
function zdw(mode::Modes.AbstractMode; λmin=100e-9, λmax=3000e-9)
    λ00 = Modes.zdw(mode; λmin=λmin, λmax=λmax, z=0)
    function addstat!(d, Eω, Et, z, dz)
        d["zdw"] = missnan(Modes.zdw(mode, λ00; z=z))
        λ00 = d["zdw"]
    end
end

"""
    zdw(mode)

Create stats function to capture the zero-dispersion wavelength (ZDW).

!!! warning
    Since [`Modes.zdw`](@ref) is based on root-finding of a derivative, this can be slow!
"""
function zdw(modes; λmin=100e-9, λmax=3000e-9)
    λ00 = zeros(length(modes))
    for (ii, mode) in enumerate(modes)
        tmp = Modes.zdw(mode; λmin=λmin, λmax=λmax, z=0)
        if ismissing(tmp)
            λ00[ii] = λmin
        else
            λ00[ii] = tmp
        end
    end
    function addstat!(d, Eω, Et, z, dz)
        d["zdw"] = [missnan(Modes.zdw(modes[ii], λ00[ii]; z=z)) for ii in eachindex(modes)]
        λ00 .= d["zdw"]
    end
end

# convert missing to NaN
missnan(x) = ismissing(x) ? NaN : x

function zdz!(d, Eω, Et, z, dz)
    d["z"] = z
    d["dz"] = dz
end

"""
    plan_analytic(grid, Eω)

Plan a transform from the frequency-domain field `Eω` to the analytic time-domain field.

Returns both a buffer for the analytic field and a closure to do the transform.
"""
function plan_analytic(grid::Grid.EnvGrid, Eω)
    Eta = similar(Eω)
    Utils.loadFFTwisdom()
    iFT = FFTW.plan_ifft(copy(Eω), 1, flags=settings["fftw_flag"])
    Utils.saveFFTwisdom()
    function analytic!(Eta, Eω)
        mul!(Eta, iFT, Eω) # for envelope fields, we only need to do the inverse transform
    end
    return Eta, analytic!
end

function plan_analytic(grid::Grid.RealGrid, Eω)
    s = collect(size(Eω))
    s[1] = (length(grid.ω) - 1)*2 # e.g. for 4097 rFFT samples, we need 8192 FFT samples
    Eta = Array{ComplexF64, ndims(Eω)}(undef, Tuple(s))
    Eωa = zero(Eta)
    idxhi = CartesianIndices(size(Eω)[2:end]) # index over all other dimensions
    Utils.loadFFTwisdom()
    iFT = FFTW.plan_ifft(Eωa, 1, flags=settings["fftw_flag"])
    Utils.saveFFTwisdom()
    function analytic!(Eta, Eω)
        copyto_fft!(Eωa, Eω, idxhi) # copy across to FFT-sampled buffer
        mul!(Eta, iFT, Eωa) # now do the inverse transform
    end
    return Eta, analytic!
end

"""
    copyto_fft!(Eωa, Eω, idxhi)

Copy the rFFT-sampled field `Eω` to the FFT-sampled buffer `Eωa`, ready for inverse FFT
"""
function copyto_fft!(Eωa, Eω, idxhi)
    n = size(Eω, 1)-1 # rFFT has sample at +fs/2, but FFT does not (only at -fs/2)
    for idx in idxhi
        for i in 1:n
            Eωa[i, idx] = 2*Eω[i, idx]
        end
    end
end

"""
    collect_stats(grid, Eω, funcs...)

Create a closure which collects statistics from the individual functions in `funcs`.

Each function given will be called with the arguments `(d, Eω, Et, z, dz)`, where
- d -> dictionary to store statistics values. each `func` should **mutate** this
- Eω -> frequency-domain field
- Et -> analytic time-domain field
- z -> current propagation distance
- dz -> current stepsize
"""
function collect_stats(grid, Eω, funcs...)
    # make sure z and dz are recorded
    if !(zdz! in funcs)
        funcs = (funcs..., zdz!)
    end
    Et, analytic! = plan_analytic(grid, Eω)
    f = let funcs=funcs
        function collect_stats(Eω, z, dz)
            d = Dict{String, Any}()
            analytic!(Et, Eω)
            for func in funcs
                func(d, Eω, Et, z, dz)
            end
            return d
        end
    end
    return f
end

function default(grid, Eω, mode::Modes.AbstractMode, linop, transform;
                 windows=nothing, gas=nothing, onaxis=false, userfuns=Any[])
    _, energyfunω = Fields.energyfuncs(grid)
    funs = [ω0(grid), energy(grid, energyfunω), peakpower(grid),
            fwhm_t(grid), zdw_linop(mode, linop),
            density(transform.densityfun)]
    if !isnothing(gas)
        push!(funs, pressure(transform.densityfun, gas))
    end
    if onaxis
        push!(funs, peakintensity(grid, (mode,)))
    else
        push!(funs, peakintensity(grid, transform.aeff))
    end
    for resp in transform.resp
        if resp isa PlasmaCumtrapz
            ir = resp.ratefunc
            if onaxis
                push!(funs, electrondensity(grid, resp.ratefunc, transform.densityfun, (mode,)))
            else
                push!(funs, electrondensity(grid, resp.ratefunc, transform.densityfun, transform.aeff))
            end
        end
    end
    if !isnothing(windows)
        for win in windows
            push!(funs, energy_λ(grid, energyfunω, win))
        end
    end
    for (idx, uf) in enumerate(userfuns)
        if uf in funs
            @warn("userfun $idx is already present in the default set and will be ignored")
        else
            push!(funs, uf)
        end
    end
    collect_stats(grid, Eω, funs...)
end

function default(grid, Eω, modes::Modes.ModeCollection, linop, transform;
                 windows=nothing, gas=nothing, userfuns=Any[])
    _, energyfunω = Fields.energyfuncs(grid)
    pol = transform.ts.indices == 1:2 ? :xy : transform.ts.indices == 1 ? :x : :y
    funs = [ω0(grid), energy(grid, energyfunω), peakpower(grid),
            peakintensity(grid, modes, components=pol), fwhm_t(grid),
            zdw_linop(modes, linop), density(transform.densityfun),
            fwhm_r(grid, modes; components=pol)]
    if !isnothing(gas)
        push!(funs, pressure(transform.densityfun, gas))
    end
    for resp in transform.resp
        if resp isa PlasmaCumtrapz
            ir = resp.ratefunc
            ed = electrondensity(grid, resp.ratefunc, transform.densityfun, modes,
                                 components=pol)
            push!(funs, ed)
        end
    end
    if !isnothing(windows)
        for win in windows
            push!(funs, energy_λ(grid, energyfunω, win))
        end
    end
    for (idx, uf) in enumerate(userfuns)
        if uf in funs
            @warn("userfun $uf is already present in the default set and will be ignored")
        else
            push!(funs, uf)
        end
    end
    collect_stats(grid, Eω, funs...)
end

# For constant linop, ZDW is also constant
function zdw_linop(mode::Modes.AbstractMode, linop::AbstractArray)
    zdw = missnan(Modes.zdw(mode))
    (d, Eω, Et, z, dz) -> d["zdw"] = zdw
end

function zdw_linop(modes, linop::AbstractArray)
    zdw = [missnan(Modes.zdw(mode)) for mode in modes]
    (d, Eω, Et, z, dz) -> d["zdw"] = zdw
end

zdw_linop(mode_s, linop) = zdw(mode_s)

end


text:
using PyPlot
A = imread("sample_photo.png")
imshow(A);

println("A is of type $(typeof(A)) with dimensions $(size(A))")
println("The values of A range between $(minimum(A)) and $(maximum(A))")

Ar = copy(A); Ar[:,:,[2,3]] .= 0.0;
Ag = copy(A); Ag[:,:,[1,3]] .= 0.0;
Ab = copy(A); Ab[:,:,[1,2]] .= 0.0;
subplot(2,2,1); imshow(Ar);
subplot(2,2,2); imshow(Ag);
subplot(2,2,3); imshow(Ab);
subplot(2,2,4); imshow(A);

using Statistics
B = mean(A, dims=3)[:,:,1]

function imshow_scale(A)
    # Like imshow(A) but scales the values to [0,1] and supports grayscale

    A .-= minimum(A)            # Scale and shift to [0,1]
    A ./= maximum(A)
    if ndims(A) < 3
        A = reshape(A, size(A,1), size(A,2), 1)
    end
    if size(A,3) == 1
        A = repeat(A, 1, 1, 3)  # Set R=G=B for grayscale
    end
    imshow(A)
end

imshow_scale(B);

function image_threshold(A, th)
    return Float64.(A .> th)
end

imshow_scale(1.0 .- image_threshold(B, 0.8));


text:
using Point2DTools
using Test

# Computational Geometry, Third Edition, P.243 - not in order
const hull1 = Point2D[(0.1, 0.35), (0.13, 0.22), (0.07, 0.15), (0.12, 0.3), (0.16,0.2)]
const hull1vx = Point2D[(0.07, 0.15), (0.1, 0.35), (0.16, 0.2), (0.07, 0.15)]

const rpoints = Point2D[[408, 22], [424, 65], [484, 416], [283, 161], [98, 358], [126, 288], [309, 56], [63, 410], [221, 8] ,[71, 378] ,[73, 138], [305, 105], [397, 475], [376, 212], [292, 304], [498, 303]]

@testset "Point2DTools.jl" begin
   @test hull1[convexperm(hull1)] == hull1vx
   @test convex(hull1) == hull1vx
   @test convexperm(rpoints) == ([8, 13, 3, 16, 2, 1, 9, 11, 8])
   @test bounds([Point2D(0,0), Point2D(1,2), Point2D(2,1)]) == [Point2D(0,0), Point2D(2,2)]
   @test scalefns(Point2D(0.,0.), Point2D(10., 10.), 10, 100)[2](10) == 10.0
end


text:
export Normal

doc"""
    Normal(mean, var)

Generator of ndarray with a normal distribution.

# Arguments
* mean: Mean of the random values.
* var: Variance of the random values.
"""
struct Normal
    mean
    var

    function Normal(mean, var)
        @assert var >= 0
        new(mean, var)
    end
end

function (n::Normal)(::Type{T}, dims::Int...) where T
    r = randn(T, dims)
    r .*= sqrt(n.var)
    r .+= n.mean
    r
end


text:
import Distributions: ContinuousUnivariateDistribution,
                      cdf,
                      logcdf,
                      pdf,
                      logpdf,
                      insupport

import Statistics: quantile

export RadialUniform,
       PoissonInvariant

"""
    RadialUniform(r_in, r_out) <: ContinuousUnivariateDistribution

A radially uniform distribution from `r_in` to `r_out`.

This distribution is derived from the change of variables evaluation of

``p(r, \\theta) = p(x, y)\\left|\\nabla_{r,\\theta}(x, y) \\right|``

which given ``p(x, y)\\propto 1``, leads to ``p(r, \\theta)\\propto r``

# Form

``p(r) = \\frac{2r}{r_\\text{out}^2 - r_\\text{in}^2}\\quad\\forall{r\\in (r_\\text{in}, r_\\text{out})}``

# Supported Functions

These functions have been explicitly written for `RadialUniform` from [Distributions.jl](https://github.com/juliastats/Distributions.jl). There may be more functionality available from fallback methods, but the following are guaranteed to work.
* [`pdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.pdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
* [`logpdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.logpdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
* [`cdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.cdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
* [`logcdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.logcdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
* [`quantile`](https://juliastats.org/Distributions.jl/stable/univariate/#Statistics.quantile-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
* [`minimum`](https://juliastats.org/Distributions.jl/stable/univariate/#Base.minimum-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport})
* [`maximum`](https://juliastats.org/Distributions.jl/stable/univariate/#Base.maximum-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport})

# Examples
```jldoctest
julia> using Distributions

julia> dist = RadialUniform(0, 10)
RadialUniform{Float64}(r_in=0.0, r_out=10.0)

julia> pdf(dist, -1)
0.0

julia> pdf(dist, 3)
0.06

julia> cdf(dist, quantile(dist, 0.5))
0.5000000000000001
```
"""
struct RadialUniform{T <: AbstractFloat} <: ContinuousUnivariateDistribution
    r_in::T
    r_out::T
    function RadialUniform(r_in::S, r_out::S) where S
        T = float(S)
        new{T}(T(r_in), T(r_out))
    end
end

RadialUniform(r_in, r_out) = RadialUniform(promote(r_in, r_out)...)

Base.minimum(d::RadialUniform) = d.r_in
Base.maximum(d::RadialUniform) = d.r_out

quantile(d::RadialUniform, q::Real) = sqrt((d.r_out^2 - d.r_in^2) * q + d.r_in^2)

function cdf(d::RadialUniform{T}, x::Real) where {T}
    x < d.r_in && return zero(T)
    x > d.r_out && return one(T)
    return (x^2 - d.r_in^2) / (d.r_out^2 - d.r_in^2)
end

function logcdf(d::RadialUniform{T}, x::Real) where {T}
    x < d.r_in && return -T(Inf)
    x > d.r_out && return zero(T)
    return log(x^2 - d.r_in^2) - log(d.r_out^2 - d.r_in^2)
end

pdf(d::RadialUniform{T}, x::Real) where T = insupport(d, x) ? 2 * x / (d.r_out^2 - d.r_in^2) : zero(T)
logpdf(d::RadialUniform{T}, x::Real) where T = insupport(d, x) ? log(2x) - log(d.r_out^2 - d.r_in^2) : -T(Inf)

# """
#     PoissonInvariant(a, b) <: ContinuousUnivariateDistribution

# A Poisson invariant distributions truncated from `a` to `b`.

# This distribution is derived from the determinant of the Fisher information matrix

# ``p(x) \\propto \\sqrt{E\\left[\\left(\\frac{d\\ln{L}}{dx} \\right)^2 \\right]}``

# which gives ``p(x)\\propto 1/\\sqrt{x}``

# # Form

# ``p(x) = \\frac{1}{2(\\sqrt{b} - \\sqrt{a})\\sqrt{x}}\\quad\\forall{x\\in (a,b)}``

# # Supported Functions

# These functions have been explicitly written for `RadialUniform` from [Distributions.jl](https://github.com/juliastats/Distributions.jl). There may be more functionality available from fallback methods, but the following are guaranteed to work.
# * [`pdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.pdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
# * [`logpdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.logpdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
# * [`cdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.cdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
# * [`logcdf`](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.logcdf-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
# * [`quantile`](https://juliastats.org/Distributions.jl/stable/univariate/#Statistics.quantile-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport,Real})
# * [`minimum`](https://juliastats.org/Distributions.jl/stable/univariate/#Base.minimum-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport})
# * [`maximum`](https://juliastats.org/Distributions.jl/stable/univariate/#Base.maximum-Tuple{Distribution{Univariate,S}%20where%20S%3C:ValueSupport})

# # Examples
# ```jldoctest
# julia> using Distributions

# julia> dist = PoissonInvariant(0, 1)
# PoissonInvariant{Float64}(a=0.0, b=1.0)

# julia> pdf(dist, -1)
# 0.0

# julia> pdf(dist, 3)
# 0.0

# julia> cdf(dist, quantile(dist, 0.5))
# 0.5
# ```
# """
# struct PoissonInvariant{T <: AbstractFloat} <: ContinuousUnivariateDistribution
#     a::T
#     b::T
#     function PoissonInvariant(a::S, b::S) where S
#         T = float(S)
#         new{T}(T(a), T(b))
#     end
# end

# PoissonInvariant(a, b) = PoissonInvariant(promote(a, b)...)

# Base.minimum(d::PoissonInvariant) = d.a
# Base.maximum(d::PoissonInvariant) = d.b

# quantile(d::PoissonInvariant, q::Real) = (q * (sqrt(d.b) - sqrt(d.a)) + sqrt(d.a))^2

# function cdf(d::PoissonInvariant{T}, x::Real) where T
#     x < d.a && return zero(T)
#     x > d.b && return one(T)
#     return (sqrt(x) - sqrt(d.a)) / (sqrt(d.b) - sqrt(d.a))
# end

# function logcdf(d::PoissonInvariant{T}, x::Real) where T
#     x < d.a && return -T(Inf)
#     x > d.b && return zero(T)
#     return log(sqrt(x) - sqrt(d.a)) - log(sqrt(d.b) - sqrt(d.a))
# end

# pdf(d::PoissonInvariant{T}, x::Real) where T = insupport(d, x) ? 1 / (2sqrt(x) * (sqrt(d.b) - sqrt(d.a))) : zero(T)
# logpdf(d::PoissonInvariant{T}, x::Real) where T = insupport(d, x) ? -log(2sqrt(x) * (sqrt(d.b) - sqrt(d.a))) : -T(Inf)