text:
export state_space, ndofs, get_state,
	get_u, get_̇u, get_v, get_̇v, get_w, get_̇w,
	get_θx, get_θy, get_θz, get_̇θx, get_̇θy, get_̇θz,
	set_u, set_v, set_w, set_̇u, set_̇v, set_̇w,
	set_θx, set_θy, set_θz, set_̇θx, set_̇θy, set_̇θz

"""
```
	ndofs(acft::Aircraft)
```

Get number of degrees of freedom for an aircraft
"""
@inline ndofs(acft::Aircraft{Fg}) where Fg = size(acft.points, 2) * 12

"""
```
	get_u(q::AbstractVector, i::Int64)
```

Get x axis position of a point in the flexible structure
"""
@inline get_u(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 1]
"""
```
	get_v(q::AbstractVector, i::Int64)
```

Get y axis position of a point in the flexible structure
"""
@inline get_v(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 2]
"""
```
	get_w(q::AbstractVector, i::Int64)
```

Get z axis position of a point in the flexible structure
"""
@inline get_w(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 3]

"""
```
	get_θx(q::AbstractVector, i::Int64)
```

Get x axis angular defflection of a point in the flexible structure
"""
@inline get_θx(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 4]
"""
```
	get_θy(q::AbstractVector, i::Int64)
```

Get y axis angular defflection of a point in the flexible structure
"""
@inline get_θy(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 5]
"""
```
	get_θz(q::AbstractVector, i::Int64)
```

Get z axis angular defflection of a point in the flexible structure
"""
@inline get_θz(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 6]

"""
```
	get_̇u(q::AbstractVector, i::Int64)
```

Get derivative in time of x axis position of a point in the flexible structure
"""
@inline get_̇u(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 1 + length(q) ÷ 2]
"""
```
	get_̇v(q::AbstractVector, i::Int64)
```

Get derivative in time of y axis position of a point in the flexible structure
"""
@inline get_̇v(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 2 + length(q) ÷ 2]
"""
```
	get_̇w(q::AbstractVector, i::Int64)
```

Get derivative in time of z axis position of a point in the flexible structure
"""
@inline get_̇w(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 3 + length(q) ÷ 2]

"""
```
	get_̇̇θx(q::AbstractVector, i::Int64)
```

Get derivative in time of x axis angular defflection of a point in the flexible structure
"""
@inline get_̇θx(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 4 + length(q) ÷ 2]
"""
```
	get_̇̇θy(q::AbstractVector, i::Int64)
```

Get derivative in time of y axis angular defflection of a point in the flexible structure
"""
@inline get_̇θy(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 5 + length(q) ÷ 2]
"""
```
	get_̇̇θz(q::AbstractVector, i::Int64)
```

Get derivative in time of z axis angular defflection of a point in the flexible structure
"""
@inline get_̇θz(q::AbstractVector, i::Int64) = q[(i - 1) * 6 + 6 + length(q) ÷ 2]

"""
```
	set_u(q::AbstractVector, i::Int64, v)
```

Set x axis defflection of a point in the flexible structure
"""
@inline set_u(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 1] = v end
"""
```
	set_v(q::AbstractVector, i::Int64, v)
```

Set y axis defflection of a point in the flexible structure
"""
@inline set_v(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 2] = v end
"""
```
	set_w(q::AbstractVector, i::Int64, v)
```

Set z axis defflection of a point in the flexible structure
"""
@inline set_w(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 3] = v end

"""
```
	set_θx(q::AbstractVector, i::Int64, v)
```

Set x axis angular defflection of a point in the flexible structure
"""
@inline set_θx(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 4] = v end
"""
```
	set_θy(q::AbstractVector, i::Int64, v)
```

Set y axis angular defflection of a point in the flexible structure
"""
@inline set_θy(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 5] = v end
"""
```
	set_θz(q::AbstractVector, i::Int64, v)
```

Set z axis angular defflection of a point in the flexible structure
"""
@inline set_θz(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 6] = v end

"""
```
	set_̇u(q::AbstractVector, i::Int64, v)
```

Set derivative in time of x axis defflection of a point in the flexible structure
"""
@inline set_̇u(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 1 + length(q) ÷ 2] = v end
"""
```
	set_̇v(q::AbstractVector, i::Int64, v)
```

Set derivative in time of y axis defflection of a point in the flexible structure
"""
@inline set_̇v(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 2 + length(q) ÷ 2] = v end
"""
```
	set_̇w(q::AbstractVector, i::Int64, v)
```

Set derivative in time of z axis defflection of a point in the flexible structure
"""
@inline set_̇w(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 3 + length(q) ÷ 2] = v end

"""
```
	set_̇θx(q::AbstractVector, i::Int64, v)
```

Set derivative in time of x axis angular defflection of a point in the flexible structure
"""
@inline set_̇θx(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 4 + length(q) ÷ 2] = v end
"""
```
	set_̇θy(q::AbstractVector, i::Int64, v)
```

Set derivative in time of y axis angular defflection of a point in the flexible structure
"""
@inline set_̇θy(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 5 + length(q) ÷ 2] = v end
"""
```
	set_̇θz(q::AbstractVector, i::Int64, v)
```

Set derivative in time of z axis angular defflection of a point in the flexible structure
"""
@inline set_̇θz(q::AbstractVector, i::Int64, v) = begin q[(i - 1) * 6 + 6 + length(q) ÷ 2] = v end

"""
```
	get_u(q::AbstractVector)
```

Get vector of x axis defflection of all points in the flexible structure
"""
@inline get_u(q::AbstractVector) = q[1:6:(length(q) ÷ 2)]
"""
```
	get_v(q::AbstractVector)
```

Get vector of y axis defflection of all points in the flexible structure
"""
@inline get_v(q::AbstractVector) = q[2:6:(length(q) ÷ 2)]
"""
```
	get_w(q::AbstractVector)
```

Get vector of z axis defflection of all points in the flexible structure
"""
@inline get_w(q::AbstractVector) = q[3:6:(length(q) ÷ 2)]

"""
```
	get_θx(q::AbstractVector)
```

Get vector of z axis defflection of all points in the flexible structure
"""
@inline get_θx(q::AbstractVector) = q[4:6:(length(q) ÷ 2)]
"""
```
	get_θy(q::AbstractVector)
```

Get vector of y axis defflection of all points in the flexible structure
"""
@inline get_θy(q::AbstractVector) = q[5:6:(length(q) ÷ 2)]
"""
```
	get_θz(q::AbstractVector)
```

Get vector of z axis defflection of all points in the flexible structure
"""
@inline get_θz(q::AbstractVector) = q[6:6:(length(q) ÷ 2)]

"""
```
	get_̇u(q::AbstractVector)
```

Get vector of derivatives in time of x axis defflection of all points in the flexible structure
"""
@inline get_̇u(q::AbstractVector) = q[(1 + length(q) ÷ 2):6:end]
"""
```
	get_̇v(q::AbstractVector)
```

Get vector of derivatives in time of y axis defflection of all points in the flexible structure
"""
@inline get_̇v(q::AbstractVector) = q[(2 + length(q) ÷ 2):6:end]
"""
```
	get_̇w(q::AbstractVector)
```

Get vector of derivatives in time of z axis defflection of all points in the flexible structure
"""
@inline get_̇w(q::AbstractVector) = q[(3 + length(q) ÷ 2):6:end]

"""
```
	get_̇θx(q::AbstractVector)
```

Get vector of derivatives in time of x axis angular defflection of all points in the flexible structure
"""
@inline get_̇θx(q::AbstractVector) = q[(4 + length(q) ÷ 2):6:end]
"""
```
	get_̇θy(q::AbstractVector)
```

Get vector of derivatives in time of y axis angular defflection of all points in the flexible structure
"""
@inline get_̇θy(q::AbstractVector) = q[(5 + length(q) ÷ 2):6:end]
"""
```
	get_̇θz(q::AbstractVector)
```

Get vector of derivatives in time of z axis angular defflection of all points in the flexible structure
"""
@inline get_̇θz(q::AbstractVector) = q[(6 + length(q) ÷ 2):6:end]

"""
```
	set_u(q::AbstractVector, v)
```

Set value of all x axis defflections in the structure
"""
@inline set_u(q::AbstractVector, v) = begin q[1:6:(length(q) ÷ 2)] .= v end
"""
```
	set_v(q::AbstractVector, v)
```

Set value of all y axis defflections in the structure
"""
@inline set_v(q::AbstractVector, v) = begin q[2:6:(length(q) ÷ 2)] .= v end
"""
```
	set_w(q::AbstractVector, v)
```

Set value of all z axis defflections in the structure
"""
@inline set_w(q::AbstractVector, v) = begin q[3:6:(length(q) ÷ 2)] .= v end

"""
```
	set_θx(q::AbstractVector, v)
```

Set value of all x axis angular defflections in the structure
"""
@inline set_θx(q::AbstractVector, v) = begin q[4:6:(length(q) ÷ 2)] .= v end
"""
```
	set_θy(q::AbstractVector, v)
```

Set value of all y axis angular defflections in the structure
"""
@inline set_θy(q::AbstractVector, v) = begin q[5:6:(length(q) ÷ 2)] .= v end
"""
```
	set_θz(q::AbstractVector, v)
```

Set value of all z axis angular defflections in the structure
"""
@inline set_θz(q::AbstractVector, v) = begin q[6:6:(length(q) ÷ 2)] .= v end

"""
```
	set_̇u(q::AbstractVector, v)
```

Set time derivative of all x axis defflections in the structure
"""
@inline set_̇u(q::AbstractVector, v) = begin q[(1 + length(q) ÷ 2):6:end] .= v end
"""
```
	set_̇v(q::AbstractVector, v)
```

Set time derivative of all y axis defflections in the structure
"""
@inline set_̇v(q::AbstractVector, v) = begin q[(2 + length(q) ÷ 2):6:end] .= v end
"""
```
	set_̇w(q::AbstractVector, v)
```

Set time derivative of all z axis defflections in the structure
"""
@inline set_̇w(q::AbstractVector, v) = begin q[(3 + length(q) ÷ 2):6:end] .= v end

"""
```
	set_̇θx(q::AbstractVector, v)
```

Set time derivative of all x axis angular defflections in the structure
"""
@inline set_̇θx(q::AbstractVector, v) = begin q[(4 + length(q) ÷ 2):6:end] .= v end
"""
```
	set_̇θy(q::AbstractVector, v)
```

Set time derivative of all y axis angular defflections in the structure
"""
@inline set_̇θy(q::AbstractVector, v) = begin q[(5 + length(q) ÷ 2):6:end] .= v end
"""
```
	set_̇θz(q::AbstractVector, v)
```

Set time derivative of all z axis angular defflections in the structure
"""
@inline set_̇θz(q::AbstractVector, v) = begin q[(6 + length(q) ÷ 2):6:end] .= v end

"""
Function to get incidences for each wing strip
"""
get_incidences(acft::Aircraft{Fg}, q::AbstractVector) where Fg = [
	(
		strp.Mtouni[1, 1] * (get_θx(q, strp.ipt1) + get_θx(q, strp.ipt2)) / 2 +
		strp.Mtouni[2, 1] * (get_θy(q, strp.ipt1) + get_θy(q, strp.ipt2)) / 2 +
		strp.Mtouni[3, 1] * (get_θz(q, strp.ipt1) + get_θz(q, strp.ipt2)) / 2
	) for strp in acft.wing_strips
]

function get_freestream_velocities(acft::Aircraft{Fg}, q::AbstractVector) where Fg

	#=
	[
		- (
			[
				(get_̇u(q, strp.ipt1) + get_̇u(q, strp.ipt2)) / 2,
				(get_̇v(q, strp.ipt1) + get_̇v(q, strp.ipt2)) / 2,
				(get_̇w(q, strp.ipt1) + get_̇w(q, strp.ipt2)) / 2
			] .+ cross(
				[
					(get_̇θx(q, strp.ipt1) + get_̇θx(q, strp.ipt2)) / 2,
					(get_̇θy(q, strp.ipt1) + get_̇θy(q, strp.ipt2)) / 2,
					(get_̇θz(q, strp.ipt1) + get_̇θz(q, strp.ipt2)) / 2
				],
				[strp.c / 2, 0.0, 0.0]
			)
		) for strp in acft.wing_strips
	]
	=#

	[
		(
			Real[
				- (get_̇u(q, strp.ipt1) + get_̇u(q, strp.ipt2)) / 2,
				- (get_̇v(q, strp.ipt1) + get_̇v(q, strp.ipt2)) / 2,
				- (get_̇w(q, strp.ipt1) + get_̇w(q, strp.ipt2)) / 2
			]
		) for strp in acft.wing_strips
	]

end

"""
Function to get vortex kins, collocation points and normals
"""
function vortex_info(acft::Aircraft{Fg}, q::AbstractVector, θ::AbstractVector) where Fg

	kins = [
		(
			[
				get_u(q, strp.ipt1),
				get_v(q, strp.ipt1),
				get_w(q, strp.ipt1)
			],
			[
				get_u(q, strp.ipt2),
				get_v(q, strp.ipt2),
				get_w(q, strp.ipt2)
			]
		) for strp in acft.wing_strips
	]

	cols = [
		(
			(k[1] .+ k[2]) ./ 2 # .+ [strp.c / 2, 0.0, 0.0]
		) for (k, strp) in zip(kins, acft.wing_strips)
	]

	norms = [
		strp.Mtouni * [
			0.0,
			sin(θt),
			- cos(θt)
		] for (θt, strp) in zip(θ, acft.wing_strips)
	]

	kins, cols, norms

end

"""
Get fuselage section singularities
"""
function fuselage_singularities(acft::Aircraft{Fg}, q::AbstractVector) where Fg

	ls = [
		(
			acft.points[:, sect.ipt2] .- acft.points[:, sect.ipt1]
		) for sect in acft.fuselage_cuts
	]

	Ls = [
		begin
			L = norm(l)

			l ./= L

			L
		end for l in ls
	]

	σ = [
		begin
			V = [
				- (get_̇u(q, sect.ipt1) + get_̇u(q, sect.ipt2)) / 2,
				- (get_̇v(q, sect.ipt1) + get_̇v(q, sect.ipt2)) / 2,
				- (get_̇w(q, sect.ipt1) + get_̇w(q, sect.ipt2)) / 2
			]

			(
				(sect.A2 - sect.A1) * (
					ls[i] ⋅ V
				) + Ls[i] * (sqrt(sect.A2 * π) + sqrt(sect.A1 * π)) * sect.Cf * norm(V) / 2
			) / Ls[i]
		end for (i, sect) in enumerate(acft.fuselage_cuts)
	]

	μ = [
		(sect.A1 + sect.A2) .* (
			begin
				ξ = [
					(get_̇u(q, sect.ipt1) + get_̇u(q, sect.ipt2)) / 2,
					(get_̇v(q, sect.ipt1) + get_̇v(q, sect.ipt2)) / 2,
					(get_̇w(q, sect.ipt1) + get_̇w(q, sect.ipt2)) / 2
				]

				ξ .- (ξ ⋅ ls[i]) .* ls[i]
			end
		) for (i, sect) in enumerate(acft.fuselage_cuts)
	]

	σ, μ

end

"""
Find singularity values for punctual drag components
"""
function point_drag_singularities(acft::Aircraft{Fg}, q::AbstractVector) where Fg

	vs = [
		[
			get_̇u(q, pmass.ipt),
			get_̇v(q, pmass.ipt),
			get_̇w(q, pmass.ipt)
		] for pmass in acft.point_drags
	]

	σ = [
		norm(v) * pmass.A * pmass.CD for (v, pmass) in zip(vs, acft.point_drags)
	]

	μ = [
		v .* pmass.V for (v, pmass) in zip(vs, acft.point_drags)
	]

	σ, μ

end

"""
```
	function get_state(
		acft::Aircraft{Fg},
		U∞::Fg;
		α::Fg = 0.0,
		β::Fg = 0.0,
		p::Fg = 0.0,
		q::Fg = 0.0,
		r::Fg = 0.0
	) where Fg
```

Function to instantiate state vector with a given flight condition

* `acft`: aircraft at hand
* `U∞`: freestream velocity
* `α`: angle of attack
* `β`: sideslip angle
* `p`: x axis angular velocity (in the stability coordinate system)
* `q`: y axis angular velocity (in the stability coordinate system)
* `q`: z axis angular velocity (in the stability coordinate system)
"""
function get_state(
	acft::Aircraft{Fg},
	U∞::Fg;
	α::Fg = 0.0,
	β::Fg = 0.0,
	p::Fg = 0.0,
	q::Fg = 0.0,
	r::Fg = 0.0
) where Fg

	V = [
		- cosd(α) * cosd(β) * U∞,
		- cosd(α) * sind(β) * U∞,
		- sind(α) * U∞
	]
	ω = [
		- p,
		q,
		- r
	]

	vcat(
		[
			[
				acft.points[:, ipt];
				0.0;
				0.0;
				0.0
			] for ipt in 1:size(acft.points, 2)
		]...,
		[
			[
				V .+ ω × acft.points[:, ipt];
				ω
			] for ipt = 1:size(acft.points, 2)
		]...
	)

end

"""
Get coefficients for wing sections
"""
function get_section_forces(
	acft::Aircraft{Fg},
	Γs, U∞,
	us, u0s, ρ, kins,
	∂Γ!∂t
) where Fg

	Γr = (
		isnothing(∂Γ!∂t) ?
		Γs :
		[
			s.c * Γt / norm(u) for (i, (u, s, Γt)) in enumerate(zip(us, acft.wing_strips, ∂Γ!∂t))
		] .+ Γs
	)

	Fs = [
		cross(
			u, k[2] .- k[1]
		) .* ρ .* Γ for (u, k, Γ) in zip(us, kins, Γr)
	]

	CLs = [
		2 * Γs[i] / (U∞ * sect.c) for (i, sect) in enumerate(acft.wing_strips)
	]
	Cms = [
		sect.afl.Cm0 + sect.afl.Cmα * (CL - sect.afl.CL0) / sect.afl.CLα for (sect, CL) in zip(acft.wing_strips, CLs)
	]
	CDs = [
		sect.afl.polar[3] + sect.afl.polar[2] * CL + sect.afl.polar[1] * CL ^ 2 for (sect, CL) in zip(acft.wing_strips, CLs)
	]

	for (i, (k, u, CD, F, sect)) in enumerate(zip(kins, u0s, CDs, Fs, acft.wing_strips))
		nu = norm(u)

		q = (sect.c * sqrt((k[2][2] - k[1][2]) ^ 2 + (k[2][3] - k[1][3]) ^ 2)) * ρ * nu ^ 2 / 2

		@. F += u * CD * q / nu
		CDs[i] = F ⋅ u / (nu * q)
	end

	Ms = [
		begin
			nu = norm(u)

			q = sect.c ^ 2 * ρ * nu ^ 2 / 2

			M = @. Cm * q * (k[2] - k[1])

			M
		end for (k, u, Cm, sect) in zip(kins, us, Cms, acft.wing_strips)
	]

	Fs, Ms, CLs, CDs, Cms

end

"""
Obtain fuselage loads
"""
function fuselage_loads(acft::Aircraft{Fg}, q::AbstractVector, ρ, σs, PGβ) where Fg

	[
		begin
			V = [
				- (get_̇u(q, sect.ipt1) + get_̇u(q, sect.ipt2)) / 2,
				- (get_̇v(q, sect.ipt1) + get_̇v(q, sect.ipt2)) / 2,
				- (get_̇w(q, sect.ipt1) + get_̇w(q, sect.ipt2)) / 2
			]

			(σ * sect.L * PGβ * ρ) .* V
		end for (sect, σ) in zip(acft.fuselage_cuts, σs)
	]

end

"""
Obtain loads at point drag elements
"""
function point_loads(acft::Aircraft{Fg}, q::AbstractVector, ρ, σs, PGβ) where Fg

	[
		begin
			V = [
				- get_̇u(q, pmass.ipt),
				- get_̇v(q, pmass.ipt),
				- get_̇w(q, pmass.ipt)
			]

			(σ * PGβ * ρ) .* V
		end for (pmass, σ) in zip(acft.point_drags, σs)
	]

end

"""
```
	function state_space(
		acft::Aircraft{Fg}, x::AbstractVector, U∞::Fg;
		a::Fg = 343.0,
		ρ::Fg = 1.225,
		α::Fg = 0.0,
		β::Fg = 0.0,
		p::Fg = 0.0,
		q::Fg = 0.0,
		r::Fg = 0.0,
		Sref::Fg = 1.0,
		cref::Fg = 1.0,
		bref::Fg = 1.0,
		∂Γ!∂t::Union{AbstractVector, Nothing} = nothing,
		fixed_points::Vector{Fi} = Int64[],
		engine_loads::Matrix{Fg} = Matrix{Fg}(undef, 2, 0),
		amplification_factor::Float64 = 1.0,
		aerodynamic_forces::Bool = true,
		structural_damping::Bool = true
	) where {Fg, Fi <: Int}
```

Function to get a dynamic analysis of an aircraft, given a state vector

* `acft`: aircraft object
* `q`: variables of interest
* `U∞`: freestream velocity
* `a`: sound speed
* `ρ`: density
* `α`: angle of attack
* `β`: sideslip angle
* `p`: x axis angular velocity (in the stability coordinate system)
* `q`: y axis angular velocity (in the stability coordinate system)
* `q`: z axis angular velocity (in the stability coordinate system)
* `Sref`: reference surface
* `cref`: reference chord
* `bref`: reference span
* `∂Γ!∂t`: derivative of circulation over time, for unsteady simulations
* `structural_damping`: boolean indicating whether to consider structural damping
* `engine_loads`: forces and moments at the engines, as in:

```
engine_loads = [
	F1 F2 ...;
	M1 M2 ...
]
```

* output: derivatives of variables of interest (altered in place)
* output: named tuple as in example:

```
(
	Γ = Γs,
	fs = (
		[ # elastic forces and moments acting at each beam segment
			[fx1, fy1, fz1, mx1, my1, mz1, fx2, fy2, fz2, mx2, my2, mz2],
			[fx1, fy1, fz1, mx1, my1, mz1, fx2, fy2, fz2, mx2, my2, mz2],
			...
		]
	)
	sectional = (
		CL = CLs,
		CD = CDs,
		Cm = Cms
	), # sectional coefficients for each strip
	total = (
		CL = CL,
		CD = CD,
		Cl = Cl,
		Cm = Cm,
		Cn = Cn,
		CX = CX,
		CZ = CY,
		CZ = CZ
	)
)
```
"""
function state_space(
	acft::Aircraft{Fg}, x::AbstractVector, U∞::Fg;
	a::Fg = 343.0,
	ρ::Fg = 1.225,
	α::Fg = 0.0,
	β::Fg = 0.0,
	p::Fg = 0.0,
	q::Fg = 0.0,
	r::Fg = 0.0,
	Sref::Fg = 1.0,
	cref::Fg = 1.0,
	bref::Fg = 1.0,
	∂Γ!∂t::Union{AbstractVector, Nothing} = nothing,
	fixed_points::Vector{Fi} = Int64[],
	engine_loads::Matrix{Fg} = Matrix{Fg}(undef, 2, 0),
	amplification_factor::Float64 = 1.0,
	structural_damping::Bool = true,
	aerodynamic_forces::Bool = true
) where {Fg, Fi <: Int}

	linds = [
		wsect.ipt1 for wsect in acft.wing_strips
	]
	rinds = [
		wsect.ipt2 for wsect in acft.wing_strips
	]

	cleft, cright = let cavg = _wstrip_interp(acft, :c)
		(cavg[linds], cavg[rinds])
	end
	eleft, eright = let eavg = _wstrip_interp(acft, :e)
		(eavg[linds], eavg[rinds])
	end
	xCGleft, xCGright = let xavg = _wstrip_interp(acft, :xCG)
		(xavg[linds], xavg[rinds])
	end

	dp_l, dp_r = let app_pts = (
		acft.points[:, linds] .+ acft.points[:, rinds]
	) ./ 2
		(
			acft.points[:, linds] .- app_pts,
			acft.points[:, rinds] .- app_pts
		)
	end
	dp_l[2:3, :] .= 0.0
	dp_r[2:3, :] .= 0.0

	xhat = [
		1.0, 0.0, 0.0
	]
	θ = get_incidences(acft, x)
	kins, cols, norms = vortex_info(acft, x, θ)

	_engine_loads = hcat(
		engine_loads,
		zeros(Fg, 2, length(acft.engines) - size(engine_loads, 2))
	)

	us = get_freestream_velocities(acft, x)

	add_engine_influence!(
		acft,
		x,
		_engine_loads,
		cols,
		us,
		ρ
	)

	σ, μ = fuselage_singularities(acft, x)
	σ_point, μ_point = point_drag_singularities(acft, x)

	if aerodynamic_forces
		for (ipt, pmass) in enumerate(acft.point_drags)
			pt = [
				get_u(x, pmass.ipt),
				get_v(x, pmass.ipt),
				get_w(x, pmass.ipt)
			]

			ϵ = sqrt(pmass.A / π) / 4

			for (u, colpt) in zip(us, cols)
				u .+= σ_point[ipt] .* point_source_infl(
					colpt, pt, ϵ
				)
				u .+= begin
					nμ = norm(μ_point[ipt])
					ξ = μ_point[ipt] ./ (nμ + 1e-10)

					nμ .* point_doublet_infl(colpt, pt, ξ, ϵ)
				end
			end
		end

		for (ifus, sect) in enumerate(acft.fuselage_cuts)
			pt1 = [
				get_u(x, sect.ipt1),
				get_v(x, sect.ipt1),
				get_w(x, sect.ipt1)
			]
			pt2 = [
				get_u(x, sect.ipt2),
				get_v(x, sect.ipt2),
				get_w(x, sect.ipt2)
			]

			A = (sect.A1 + sect.A2) / 2
			ϵ = sqrt(A / π) / 2

			for (u, colpt) in zip(us, cols)
				u .+= σ[ifus] .* source_line(
					colpt, pt1, pt2, ϵ
				)
				u .+= begin
					nμ = norm(μ[ifus])
					η = μ[ifus] ./ (nμ + 1e-10)

					nμ .* doublet_line(
						colpt, pt1, pt2, η, ϵ
					)
				end
			end
		end
	end

	M∞ = U∞ / a
	PGβ = sqrt(1.0 - M∞ ^ 2)

	Γ, u0s = (
		aerodynamic_forces ?
		begin
			A, b, As = linear_problem(acft, kins, cols, norms, us, xhat, PGβ)

			Γ = (A \ b) ./ PGβ
			for i = 1:length(acft.fuselage_cuts)
				σ[i] /= PGβ
				μ[i] ./= PGβ
			end
			for i = 1:length(acft.point_drags)
				σ_point[i] /= PGβ
				μ_point[i] ./= PGβ
			end

			u0s = deepcopy(us)
			add_downwash!(us, Γ, As)

			Γ, u0s
		end :
		(zeros(eltype(x), length(acft.wing_strips)), us)
	)

	Fs, Ms, CLs, CDs, Cms = get_section_forces(acft, Γ, U∞, us, u0s, ρ, kins, ∂Γ!∂t)

	fFs = fuselage_loads(acft, x, ρ, σ, PGβ)
	pFs = point_loads(acft, x, ρ, σ_point, PGβ)

	f = zeros(Real, length(x) ÷ 2)

	if aerodynamic_forces
		for (i, (strip, F, M, cl, cr, xCGl, xCGr)) in enumerate(
			zip(
				acft.wing_strips, Fs, Ms, cleft, cright, xCGleft, xCGright
			)
		)
			f1inds = (6 * (strip.ipt1 - 1) + 1):(6 * (strip.ipt1 - 1) + 3)
			f2inds = (6 * (strip.ipt2 - 1) + 1):(6 * (strip.ipt2 - 1) + 3)

			m1inds = (6 * (strip.ipt1 - 1) + 4):(6 * strip.ipt1)
			m2inds = (6 * (strip.ipt2 - 1) + 4):(6 * strip.ipt2)

			dcg_l = dp_l[:, i] .+ [cl * xCGl, 0.0, 0.0]
			dcg_r = dp_r[:, i] .+ [cr * xCGr, 0.0, 0.0]

			Ft = F ./ 2
			Mt_l = M ./ 2 .+ cross(
				Ft,
				dcg_l
			)
			Mt_r = M ./ 2 .+ cross(
				Ft,
				dcg_r
			)

			m_sect = (acft.masses[strip.ipt1] + acft.masses[strip.ipt2]) / 2
			inv_inertia = inv(m_sect.I) .* m_sect.m

			f[f1inds] .+= Ft .+ dcg_l × (inv_inertia * Mt_l)
			f[f2inds] .+= Ft .+ dcg_r × (inv_inertia * Mt_r)

			f[m1inds] .+= Mt_l
			f[m2inds] .+= Mt_r
		end

		for (cut, F) in zip(acft.fuselage_cuts, fFs)
			f1inds = (6 * (cut.ipt1 - 1) + 1):(6 * (cut.ipt1 - 1) + 3)
			f2inds = (6 * (cut.ipt2 - 1) + 1):(6 * (cut.ipt2 - 1) + 3)

			f[f1inds] .+= F ./ 2
			f[f2inds] .+= F ./ 2
		end

		for (pmass, F) in zip(acft.point_drags, pFs)
			f[(6 * (pmass.ipt - 1) + 1):(6 * (pmass.ipt - 1) + 3)] .+= F
		end
	end

	is_calculated = zeros(Bool, length(acft.beams))
	elastic_effort = Vector{Any}(undef, length(acft.beams))

	V = [
		- cosd(α) * cosd(β) * U∞,
		- cosd(α) * sind(β) * U∞,
		- sind(α) * U∞
	]
	ω = [
		- p,
		q,
		- r
	]

	for (sect, cl, cr, xCGl, xCGr, el, er) in zip(
		acft.wing_strips, cleft, cright, xCGleft, xCGright, eleft, eright
	)
		is_calculated[sect.ibeam] = true

		g = acft.beams[sect.ibeam].g * structural_damping

		θ1 = [
			get_θx(x, sect.ipt1),
			get_θy(x, sect.ipt1),
			get_θz(x, sect.ipt1)
		] .+ (
			[
				get_̇θx(x, sect.ipt1),
				get_̇θy(x, sect.ipt1),
				get_̇θz(x, sect.ipt1)
			] .- ω
		) .* g
		θ2 = [
			get_θx(x, sect.ipt2),
			get_θy(x, sect.ipt2),
			get_θz(x, sect.ipt2)
		] .+ (
			[
				get_̇θx(x, sect.ipt2),
				get_̇θy(x, sect.ipt2),
				get_̇θz(x, sect.ipt2)
			] .- ω
		) .* g

		x1 = [
			get_u(x, sect.ipt1),
			get_v(x, sect.ipt1),
			get_w(x, sect.ipt1)
		] .- acft.points[:, sect.ipt1] .+ (
			[
				get_̇u(x, sect.ipt1),
				get_̇v(x, sect.ipt1),
				get_̇w(x, sect.ipt1)
			] .- V
		) .* g
		x2 = [
			get_u(x, sect.ipt2),
			get_v(x, sect.ipt2),
			get_w(x, sect.ipt2)
		] .- acft.points[:, sect.ipt2] .+ (
			[
				get_̇u(x, sect.ipt2),
				get_̇v(x, sect.ipt2),
				get_̇w(x, sect.ipt2)
			] .- V
		) .* g

		x1 .+= cross(
			θ1,
			[cl * el, 0.0, 0.0]
		)
		x2 .+= cross(
			θ2,
			[cr * er, 0.0, 0.0]
		)

		b = acft.beams[sect.ibeam]

		pl = b.K * [
			b.Mtosys * x1;
			b.Mtosys * θ1;
			b.Mtosys * x2;
			b.Mtosys * θ2
		]

		elastic_effort[sect.ibeam] = pl

		f1inds = (6 * (sect.ipt1 - 1) + 1):(6 * (sect.ipt1 - 1) + 3)
		f2inds = (6 * (sect.ipt2 - 1) + 1):(6 * (sect.ipt2 - 1) + 3)

		m1inds = (6 * (sect.ipt1 - 1) + 4):(6 * sect.ipt1)
		m2inds = (6 * (sect.ipt2 - 1) + 4):(6 * sect.ipt2)

		m_sect = (acft.masses[sect.ipt1] + acft.masses[sect.ipt2]) / 2
		inv_inertia = inv(m_sect.I) .* m_sect.m

		fb = b.Mtouni * pl[1:3]
		mb = cross(
			fb, [cl * (xCGl - el), 0.0, 0.0]
		) .+ b.Mtouni * pl[4:6]
		f[f1inds] .+= fb .+ [cl * xCGl, 0.0, 0.0] × (inv_inertia * mb)
		f[m1inds] .+= mb

		fb .= b.Mtouni * pl[7:9]
		mb .= cross(
			fb, [cr * (xCGr - er), 0.0, 0.0]
		) .+ b.Mtouni * pl[10:12]
		f[f2inds] .+= fb .+ [cr * xCGr, 0.0, 0.0] × (inv_inertia * mb)
		f[m2inds] .+= mb
	end

	for (ibeam, beam) in enumerate(acft.beams)
		if !is_calculated[ibeam]
			is_calculated[ibeam] = true

			g = beam.g * structural_damping

			θ1 = [
				get_θx(x, beam.ipt1),
				get_θy(x, beam.ipt1),
				get_θz(x, beam.ipt1)
			] .+ (
				[
					get_̇θx(x, beam.ipt1),
					get_̇θy(x, beam.ipt1),
					get_̇θz(x, beam.ipt1)
				] .- ω
			) .* g
			θ2 = [
				get_θx(x, beam.ipt2),
				get_θy(x, beam.ipt2),
				get_θz(x, beam.ipt2)
			] .+ (
				[
					get_̇θx(x, beam.ipt2),
					get_̇θy(x, beam.ipt2),
					get_̇θz(x, beam.ipt2)
				] .- ω
			) .* g

			x1 = [
				get_u(x, beam.ipt1),
				get_v(x, beam.ipt1),
				get_w(x, beam.ipt1)
			] .- acft.points[:, beam.ipt1] .+ (
				[
					get_̇u(x, beam.ipt1),
					get_̇v(x, beam.ipt1),
					get_̇w(x, beam.ipt1)
				] .- V
			) .* g
			x2 = [
				get_u(x, beam.ipt2),
				get_v(x, beam.ipt2),
				get_w(x, beam.ipt2)
			] .- acft.points[:, beam.ipt2] .+ (
				[
					get_̇u(x, beam.ipt2),
					get_̇v(x, beam.ipt2),
					get_̇w(x, beam.ipt2)
				] .- V
			) .* g

			pl = beam.K * [
				beam.Mtosys * x1;
				beam.Mtosys * θ1;
				beam.Mtosys * x2;
				beam.Mtosys * θ2
			]

			elastic_effort[ibeam] = pl

			f1inds = (6 * (beam.ipt1 - 1) + 1):(6 * (beam.ipt1 - 1) + 3)
			f2inds = (6 * (beam.ipt2 - 1) + 1):(6 * (beam.ipt2 - 1) + 3)

			m1inds = (6 * (beam.ipt1 - 1) + 4):(6 * beam.ipt1)
			m2inds = (6 * (beam.ipt2 - 1) + 4):(6 * beam.ipt2)

			f[f1inds] .+= beam.Mtouni * pl[1:3]
			f[f2inds] .+= beam.Mtouni * pl[7:9]

			f[m1inds] .+= beam.Mtouni * pl[4:6]
			f[m2inds] .+= beam.Mtouni * pl[10:12]
		end
	end

	for (ie, eng) in enumerate(acft.engines)
		inds = (6 * (eng.ipt - 1) + 1):(6 * eng.ipt)

		θ = [
			get_θx(x, eng.ipt),
			get_θy(x, eng.ipt),
			get_θz(x, eng.ipt)
		]
		ξ⃗ = eng.ξ⃗ .+ cross(θ, eng.ξ⃗)
		ξ⃗ ./= norm(ξ⃗)

		f[inds] .+= [
			ξ⃗ .* _engine_loads[1, ie];
			ξ⃗ .* _engine_loads[2, ie]
		]
	end

	qd = similar(f)

	ident = Matrix{Float64}(1e-5 * I, 3, 3)
	for (im, mass) in enumerate(acft.masses)
		finds = (6 * (im - 1) + 1):(6 * (im - 1) + 3)
		minds = (6 * (im - 1) + 4):(6 * im)

		rp = [
			get_u(x, im),
			get_v(x, im),
			get_w(x, im)
		]
		vp = [
			get_̇u(x, im),
			get_̇v(x, im),
			get_̇w(x, im)
		]

		qd[finds] .= (
			f[finds] .-
			mass.m .* cross(ω, cross(ω, rp)) .-
			(2 * mass.m) .* cross(ω, vp)
		) ./ (mass.m + 1e-5)
		qd[minds] .= inv(mass.I .+ ident) * f[minds]
	end

	amp = amplification_factor
	Ndof2 = ndofs(acft) ÷ 2

	for f in fixed_points
		inds = (6 * (f - 1) + 1):(6 * f)
		vinds = inds .+ 1

		x0 = [acft.points[:, f]; zeros(3)]
		v0 = [V; ω]

		qd[inds] .= let (dx, dv) = (x[inds] .- x0, x[vinds] .- v0)
			- amp .* (dx .+ dv .* structural_damping)
		end
	end

	for lnk in acft.links
		if lnk.i2 == 0
			let f = lnk.i1
				inds = (6 * (f - 1) + 1):(6 * f)
				vinds = inds .+ Ndof2

				x0 = [acft.points[:, f]; zeros(3)]
				v0 = [V; ω]

				qd[inds] .= let (dx, dv) = (x[inds] .- x0, x[vinds] .- v0)
					- amp .* (dx .+ dv .* structural_damping)
				end
			end
		else false
			i1 = lnk.i1
			i2 = lnk.i2

			inds1 = (6 * (i1 - 1) + 1):(6 * i1)
			inds2 = (6 * (i2 - 1) + 1):(6 * i2)

			dp = acft.points[:, i2] .- acft.points[:, i1]

			r1 = [
				get_u(x, i1),
				get_v(x, i1),
				get_w(x, i1)
			]
			r2 = [
				get_u(x, i2),
				get_v(x, i2),
				get_w(x, i2)
			]

			θ1 = [
				get_θx(x, i1),
				get_θy(x, i1),
				get_θz(x, i1)
			]
			θ2 = [
				get_θx(x, i2),
				get_θy(x, i2),
				get_θz(x, i2)
			]

			v1 = [
				get_̇u(x, i1),
				get_̇v(x, i1),
				get_̇w(x, i1)
			]
			v2 = [
				get_̇u(x, i2),
				get_̇v(x, i2),
				get_̇w(x, i2)
			]

			ω1 = [
				get_̇θx(x, i1),
				get_̇θy(x, i1),
				get_̇θz(x, i1)
			]
			ω2 = [
				get_̇θx(x, i2),
				get_̇θy(x, i2),
				get_̇θz(x, i2)
			]

			a1 = [
				get_u(qd, i1),
				get_v(qd, i1),
				get_w(qd, i1)
			]
			a2 = [
				get_u(qd, i2),
				get_v(qd, i2),
				get_w(qd, i2)
			]

			α1 = [
				get_θx(qd, i1),
				get_θy(qd, i1),
				get_θz(qd, i1)
			]
			α2 = [
				get_θx(qd, i2),
				get_θy(qd, i2),
				get_θz(qd, i2)
			]

			m1 = acft.masses[i1]
			m2 = acft.masses[i2]

			mtot = m1 + m2

			a_constraint = - amp .* (
				r2 .- r1 .- dp .+ (v2 .- v1 .- ω × dp) .* structural_damping
			)
			α_constraint = - amp .* (
				θ2 .- θ1 .+ (ω2 .- ω1) .* structural_damping
			)

			invItot = inv(
				mtot.I .+ ident
			)
			invmtot = 1.0 / (mtot.m + 1e-5)

			qd[inds1] .= [
				a1 .- m2.m .* a_constraint .* invmtot;
				α1 .- invItot * (m2.I * α_constraint)
			]
			qd[inds2] .= [
				a2 .+ m1.m .* a_constraint .* invmtot;
				α2 .+ invItot * (m1.I * α_constraint)
			]
		end
	end

	qdot = [
		begin
			vs = x[(length(x) ÷ 2 + 1):end]

			for i = 1:6:length(vs)
				vs[i:(i + 2)] .-= V
				vs[(i + 3):(i + 5)] .-= ω
			end

			vs
		end;
		qd
	]

	F = zeros(Fg, 3)
	M = zeros(Fg, 3)
	if length(acft.wing_strips) > 0
		F = copy(Fs[1])
		M = Ms[1] .+ cross(
			(kins[1][1] .+ kins[1][2]) ./ 2, Fs[1]
		)
		for i = 2:length(acft.wing_strips)
			F .+= Fs[i]
			M .+= Ms[i] .+ cross(
				(kins[i][1] .+ kins[i][2]) ./ 2, Fs[i]
			)
		end
	end

	for (i, pF) in enumerate(pFs)
		F .+= pF
		M .+= cross(
			begin
				ipt = acft.point_drags[i].ipt

				[
					get_u(x, ipt),
					get_v(x, ipt),
					get_w(x, ipt)
				]
			end,
			pF
		)
	end

	for (i, fF) in enumerate(fFs)
		F .+= fF
		M .+= cross(
			begin
				ipt1 = acft.fuselage_cuts[i].ipt1
				ipt2 = acft.fuselage_cuts[i].ipt2

				[
					(get_u(x, ipt1) + get_u(x, ipt2)) / 2,
					(get_v(x, ipt1) + get_v(x, ipt2)) / 2,
					(get_w(x, ipt1) + get_w(x, ipt2)) / 2
				]
			end,
			fF
		)
	end

	q = U∞ ^ 2 * ρ / 2

	CX = - F[1] / (Sref * q)
	CY = F[2] / (Sref * q)
	CZ = - F[3] / (Sref * q)

	xv = [
		cosd(α) * cosd(β),
		cosd(α) * sind(β),
		sind(α)
	]

	CF = [
		- CX,
		CY,
		- CZ
	]

	CD = CF ⋅ xv
	CL = norm(
		CF .- CD .* xv
	)

	Cl = - M[1] / (q * Sref * bref)
	Cm = M[2] / (q * Sref * cref)
	Cn = - M[3] / (q * Sref * bref)

	is_stalled = [
		cl < s.afl.CLmin || cl > s.afl.CLmax for (cl, s) in zip(CLs, acft.wing_strips)
	]

	qdot, (
		Γ = Γ,
		sectional = (
			CL = CLs,
			CD = CDs,
			Cm = Cms,
			is_stalled = is_stalled
		),
		total = (
			CX = CX,
			CY = CY,
			CZ = CZ,
			CL = CL,
			CD = CD,
			Cl = Cl,
			Cm = Cm,
			Cn = Cn
		),
		f = elastic_effort
	)

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
module MixedModels

using BlockArrays
using BlockDiagonals
using CategoricalArrays
using Distributions
using GLM
using LinearAlgebra
using NamedArrays
using NLopt
using Random
using ProgressMeter
using Showoff
using SparseArrays
using StaticArrays
using Statistics
using StatsBase
using StatsModels
using Tables
using TypedTables

using LinearAlgebra: BlasFloat, BlasReal, HermOrSym, PosDefException, copytri!
using Printf: @printf, @sprintf
using GLM: Link, canonicallink

using StatsFuns: log2π

import Base: *
import GLM: dispersion, dispersion_parameter
import NLopt: Opt
import StatsBase: fit, fit!

export @formula,
       Bernoulli,
       Binomial,
       Block,
       BlockedSparse,
       DummyCoding,
       EffectsCoding,
       Gamma,
       GeneralizedLinearMixedModel,
       HelmertCoding,
       InverseGaussian,
       InverseLink,
       LinearMixedModel,
       LogitLink,
       LogLink,
       MixedModel,
       Normal,
       OptSummary,
       Poisson,
       RaggedArray,
       RandomEffectsTerm,
       ReMat,
       UniformBlockDiagonal,
       VarCorr,
       aic,
       aicc,
       bic,
       coef,
       coefnames,
       coeftable,
       cond,
       describeblocks,
       condVar,
       deviance,
       dispersion,
       dispersion_parameter,
       dof,
       dof_residual,
       fit,
       fit!,
       fitted,
       fixef,
       fulldummy,
       fnames,
       GHnorm,
       loglikelihood,
       lowerbd,
       nblocks,
       nobs,
       objective,
       parametricbootstrap,
       pirls!,
       predict,
       pwrss,
       ranef,
       refit!,
       residuals,
       response,
       shortestCovInt,
       sdest,
       setθ!,
       simulate!,
       sparse,
       statscholesky,
       std,
       stderror,
       updateL!,
       varest,
       vcov,
       zerocorr,
       zerocorr!

import Base: ==, *

abstract type MixedModel{T} <: StatsModels.RegressionModel end # model with fixed and random effects

include("utilities.jl")
include("arraytypes.jl")
include("varcorr.jl")
include("femat.jl")
include("remat.jl")
include("optsummary.jl")
include("randomeffectsterm.jl")
include("linearmixedmodel.jl")
include("gausshermite.jl")
include("generalizedlinearmixedmodel.jl")
include("mixed.jl")
include("linalg/statschol.jl")
include("linalg/cholUnblocked.jl")
include("linalg/rankUpdate.jl")
include("linalg/logdet.jl")
include("linalg.jl")
include("simulate.jl")

end # module

ext: jl


text:
using SafeTestsets

@safetestset "Grid tests" begin include("grid.jl") end


text:
# comparing EQ, MaternP, NN kernel with energetic norm for outer loop data
using LinearAlgebra
using SARA
using SARA: iSARA, random_sampling
using CovarianceFunctions
const Kernel = CovarianceFunctions
using GaussianDistributions
using Base.Threads

# generalized coefficient of determination, where errors are weighted by their standard deviation
# McFadden (in the context of logistic regression), (log?) likelihood ratio index
# x is data, y is model prediction, σ is expected noise for each data point
function likelihood_ratio(x::AbstractVector, y::AbstractVector, σ::AbstractVector = ones(length(x)))
    δ = (x .- mean(x)) ./ σ
    tss = sum(abs2, δ)
    δ = (x .- y) ./ σ
    ess = sum(abs2, δ)
    r2 = 1 - ess / tss
end

################# setup kernels for comparison ##############
# θ are stripe-dependent rescaling parameters
# LSA-bias is encoded in last parameter tuple
# l is length scale for all kernels
get_kernels(l) = get_kernels(l...)
function get_kernels(T_length::Real, τ_length::Real)
    # list of base kernels
    k1 = Kernel.EQ()
    k2 = Kernel.MaternP(2)
    k3 = Kernel.NN(1000.) # parameter governs size of bias coefficients in neural network, put proportional to temperature domain
    base_kernels = [k1, k2, k3]
    # equip base kernels with energetic norm
    # putlength(k) = Kernel.Lengthscale(k, l)
    # A = Diagonal(inv.([T_length, τ_length]).^2)
    # putenergy(k) = Kernel.Energetic(k, A)
    # kernels = putenergy.(base_kernels)
    scaling = [T_length, τ_length]
    scale(k) = (x, y) -> k(x./scaling, y./scaling)
    kernels = scale.(base_kernels)
    return kernels
end

############### generate sampling paths with different kernels
function sampling_path(G::Gaussian, x::AbstractMatrix,
                       y::AbstractMatrix,  σ²::AbstractMatrix, samplesizes)
    C = G
    nsample = length(samplesizes)
    rmse, r2 = zeros(nsample), zeros(nsample)
    indices = zeros(Int, 0)
    n, m = size(x)
    helper(x, i) = @views vec(x[:, i]) # selects the columns indexed by i, and vectorizes the result
    for (i, s) in enumerate(samplesizes) # random sampling
        indices = rand(1:m, s) # choose s random columns
        sort!(indices)
        xi, yi, σ²i = helper(x, indices), helper(y, indices), helper(σ², indices)
        C = conditional(G, xi, Gaussian(yi, σ²i), tol = 1e-6)
        μ = mean(C, vec(x))
        δ = vec(y) - μ
        rmse[i] = sqrt(mean(δ.^2 ./ vec(σ²))) # weighted root mean squared deviation
        r2[i] = likelihood_ratio(vec(y), μ, vec(σ²))
    end
    return rmse, r2
end

# using random sampling
function outer_kernel_benchmark!(rmse, r2, kernels, x, y, σ², samplesizes)
    nexp = size(rmse, 1)
    for (j, k) in enumerate(kernels)
        println("kernel $j of $(length(kernels))")
        G = Gaussian(k)
        @threads for i in 1:nexp
            d, r = sampling_path(G, x, y, σ², samplesizes)
            rmse[i, :, j] = d
            r2[i, :, j] = r
        end
    end
    return rmse
end

function energetic_kernel_benchmark(x, y, σ², nexp, nkernel, samplesizes, samplesperstripe)
    nsample = length(samplesizes) # number of different sample sizes [1, ..., nsample]
    # range of plausible length scales
    T_length = range(20, 40, step = 10)
    τ_length = range(.2, .4, step = .1)
    nlength = length(T_length) * length(τ_length)
    rmse = zeros(nexp, nsample, nkernel, length(T_length), length(τ_length))
    r2 = zeros(nexp, nsample, nkernel, length(T_length), length(τ_length))
    # first, calculate with energetic input norm
    for (i, T_l) in enumerate(T_length)
        for (j, τ_l) in enumerate(τ_length)
            println("lengthscale $(length(τ_length)*(i-1) + j) of $nlength")
            kernels = get_kernels((T_l, τ_l))
            rmse_ij = @view rmse[:, :, :, i, j]
            r2_ij = @view r2[:, :, :, i, j]
            outer_kernel_benchmark!(rmse_ij, r2_ij, kernels, x, y, σ², samplesizes)
        end
    end

    if true
        f = h5open("outer_kernel_benchmark_energetic.h5", "w")
        f["weighted rmse"] = rmse
        f["likelihood ratio"] = r2
        f["temperature lengthscales"] = collect(T_length) # in resultss
        f["dwelltime lengthscales"] = collect(τ_length)
        f["nexp"] = nexp
        f["sample sizes"] = collect(samplesizes)
        f["nlength"] = nlength
        f["nkernel"] = nkernel
        f["kernel names"] = ["EQ", "Matern_2", "NN"]
        f["dimensions"] = ["nexp", "nsample", "nkernel", "temperature lengthscales", "dwelltime lengthscales"]
        f["samples per stripe"] = samplesperstripe
        close(f)
    end
    return rmse, r2
end

function uniform_kernel_benchmark(x, y, σ², nexp, nkernel, samplesizes, samplesperstripe)
    nsample = length(samplesizes) # number of different sample sizes [1, ..., nsample]
    # range of plausible length scales
    lengthscales = [.25, .5, 1., 2., 4., 8., 16., 32.]
    rmse = zeros(nexp, nsample, nkernel, length(lengthscales))
    r2 = zeros(nexp, nsample, nkernel, length(lengthscales))
    for (i, l) in enumerate(lengthscales)
        println("lengthscale $i of $(length(lengthscales))")
        kernels = get_kernels((l, l))
        rmse_i = @view rmse[:, :, :, i]
        r2_i = @view r2[:, :, :, i]
        outer_kernel_benchmark!(rmse_i, r2_i, kernels, x, y, σ², samplesizes)
    end

    if true
        f = h5open("outer_kernel_benchmark_uniform.h5", "w")
        f["weighted rmse"] = rmse
        f["likelihood ratio"] = r2
        f["lengthscales"] = collect(lengthscales) # in resultss
        f["nexp"] = nexp
        f["sample sizes"] = collect(samplesizes)
        f["nlength"] = length(lengthscales)
        f["nkernel"] = nkernel
        f["kernel names"] = ["EQ", "Matern_2", "NN"]
        f["dimensions"] = ["nexp", "nsample", "nkernel", "lengthscales"]
        f["samples per stripe"] = samplesperstripe
        close(f)
    end
    return rmse, r2
end

function outer_kernel_benchmark_driver(x::AbstractMatrix, y::AbstractMatrix, σ²::AbstractMatrix)
    # subsample observations per
    subsample(x, k::Int) = @view x[end+1-k:end, :]
    samplesperstripe = 8
    x = subsample(x, samplesperstripe)
    y = subsample(y, samplesperstripe)
    σ² = subsample(σ², samplesperstripe)

    # standardize global data
    y = y .- mean(y)
    y_normalization = std(y)
    y /= y_normalization
    σ² /= y_normalization^2 # noise variance

    # experimental parameters
    nexp = 1 # number of independent experiments per nsample
    nkernel = 3 # number of kernels

    samplemax = 64
    samplestep = 32 # number of sample
    samplesizes = range(32, samplemax, step = samplestep)
    energetic_kernel_benchmark(x, y, σ², nexp, nkernel, samplesizes, samplesperstripe)
    uniform_kernel_benchmark(x, y, σ², nexp, nkernel, samplesizes, samplesperstripe)
end

# load outer loop data
using HDF5
path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
# path = "NatCom2020/"
# file = "Bi2O3_19F44_01_outer_loop_data.h5" # without input noise
file = "Bi2O3_19F44_01_outer_loop_input_noise_data.h5"
f = h5open(path * file, "r")
dwelltimes = read(f, "dwelltimes")
temperatures = read(f, "temperatures")
gradients = read(f, "gradients")
var_gradients = read(f, "var(gradients)")
close(f)

# pre-process data
x = tuple.(temperatures, dwelltimes)
y = gradients
σ² = var_gradients # uncertainty in y

rmse, r2 = outer_kernel_benchmark_driver(x, y, σ²)

# mrmse = reshape(mean(rmse, dims = 1), size(rmse)[2:4]..., :)
# mr2 = reshape(mean(r2, dims = 1), size(r2)[2:4]..., :)


text:
#function bodies
@createOpenGLFun wglUseFontBitmapsA(hDC::HDC, first::DWORD, count::DWORD, listBase::DWORD)::Cint
@createOpenGLFun wglGetPbufferDCARB(hPbuffer::HPBUFFERARB)::Ptr{Void}
@createOpenGLFun wglGetPixelFormatAttribfvEXT(hdc::HDC, iPixelFormat::Cint, iLayerPlane::Cint, nAttributes::Cuint, piAttributes::Ptr{Cint}, pfValues::Ptr{Cfloat})::Cint
@createOpenGLFun wglReleaseVideoDeviceNV(hVideoDevice::HPVIDEODEV)::Cint
@createOpenGLFun wglBindVideoCaptureDeviceNV(uVideoSlot::Cuint, hDevice::HVIDEOINPUTDEVICENV)::Cint
@createOpenGLFun wglGetGenlockSourceI3D(hDC::HDC, uSource::Ptr{Cuint})::Cint
@createOpenGLFun wglGetFrameUsageI3D(pUsage::Ptr{Cfloat})::Cint
@createOpenGLFun wglRealizeLayerPalette(hdc::HDC, iLayerPlane::Cint, bRealize::Bool)::Cint
@createOpenGLFun wglEnumGpusNV(iGpuIndex::Cuint, phGpu::HGPUNV)::Cint
@createOpenGLFun wglEnableGenlockI3D(hDC::HDC)::Cint
@createOpenGLFun wglGetDigitalVideoParametersI3D(hDC::HDC, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglCreateContextAttribsARB(hDC::HDC, hShareContext::HGLRC, attribList::Ptr{Cint})::Ptr{Void}
@createOpenGLFun wglReleaseVideoCaptureDeviceNV(hDc::HDC, hDevice::HVIDEOINPUTDEVICENV)::Cint
@createOpenGLFun wglEnumGpusFromAffinityDCNV(hAffinityDC::HDC, iGpuIndex::Cuint, hGpu::HGPUNV)::Cint
@createOpenGLFun wglQueryPbufferEXT(hPbuffer::HPBUFFEREXT, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglQueryVideoCaptureDeviceNV(hDc::HDC, hDevice::HVIDEOINPUTDEVICENV, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglSendPbufferToVideoNV(hPbuffer::HPBUFFERARB, iBufferType::Cint, pulCounterPbuffer::Ptr{Culong}, bBlock::Bool)::Cint
@createOpenGLFun wglDXObjectAccessNV(hObject::HANDLE, access::GLenum)::Cint
@createOpenGLFun wglEnumGpuDevicesNV(hGpu::HGPUNV, iDeviceIndex::Cuint, lpGpuDevice::PGPU_DEVICE)::Cint
@createOpenGLFun wglGetGPUIDsAMD(maxCount::Cuint, ids::Ptr{Cuint})::Cuint
@createOpenGLFun wglCopyImageSubDataNV(hSrcRC::HGLRC, srcName::GLuint, srcTarget::GLenum, srcLevel::GLint, srcX::GLint, srcY::GLint, srcZ::GLint, hDstRC::HGLRC, dstName::GLuint, dstTarget::GLenum, dstLevel::GLint, dstX::GLint, dstY::GLint, dstZ::GLint, width::GLsizei, height::GLsizei, depth::GLsizei)::Cint
@createOpenGLFun wglGetPixelFormatAttribivARB(hdc::HDC, iPixelFormat::Cint, iLayerPlane::Cint, nAttributes::Cuint, piAttributes::Ptr{Cint}, piValues::Ptr{Cint})::Cint
@createOpenGLFun wglGetCurrentReadDCARB()::Ptr{Void}
@createOpenGLFun wglIsEnabledFrameLockI3D(pFlag::Ptr{Bool})::Cint
@createOpenGLFun wglDeleteAssociatedContextAMD(hglrc::HGLRC)::Cint
@createOpenGLFun wglDXUnlockObjectsNV(hDevice::HANDLE, count::GLint, hObjects::HANDLE)::Cint
@createOpenGLFun wglUseFontBitmapsW(hDC::HDC, first::DWORD, count::DWORD, listBase::DWORD)::Cint
@createOpenGLFun wglBindSwapBarrierNV(group::GLuint, barrier::GLuint)::Cint
@createOpenGLFun wglGetCurrentReadDCEXT()::Ptr{Void}
@createOpenGLFun wglSetPbufferAttribARB(hPbuffer::HPBUFFERARB, piAttribList::Ptr{Cint})::Cint
@createOpenGLFun wglQueryFrameTrackingI3D(pFrameCount::Ptr{DWORD}, pMissedFrames::Ptr{DWORD}, pLastMissedUsage::Ptr{Cfloat})::Cint
@createOpenGLFun wglCreateContext(hDc::HDC)::Ptr{Void}
@createOpenGLFun wglWaitForMscOML(hdc::HDC, target_msc::Clonglong, divisor::Clonglong, remainder::Clonglong, ust::Ptr{Clonglong}, msc::Ptr{Clonglong}, sbc::Ptr{Clonglong})::Cint
@createOpenGLFun wglSetGammaTableParametersI3D(hDC::HDC, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglIsEnabledGenlockI3D(hDC::HDC, pFlag::Ptr{Bool})::Cint
@createOpenGLFun wglLoadDisplayColorTableEXT(table::Ptr{GLushort}, length::GLuint)::Bool
@createOpenGLFun wglSwapBuffers(hdc::HDC)::Cint
@createOpenGLFun wglBeginFrameTrackingI3D()::Cint
@createOpenGLFun wglGetPixelFormatAttribivEXT(hdc::HDC, iPixelFormat::Cint, iLayerPlane::Cint, nAttributes::Cuint, piAttributes::Ptr{Cint}, piValues::Ptr{Cint})::Cint
@createOpenGLFun wglGetDefaultProcAddress(lpszProc::LPCSTR)::Ptr{Void}
@createOpenGLFun wglCreateAssociatedContextAttribsAMD(id::Cuint, hShareContext::HGLRC, attribList::Ptr{Cint})::Ptr{Void}
@createOpenGLFun wglMakeContextCurrentEXT(hDrawDC::HDC, hReadDC::HDC, hglrc::HGLRC)::Cint
@createOpenGLFun wglGetSwapIntervalEXT()::Cint
@createOpenGLFun wglMakeAssociatedContextCurrentAMD(hglrc::HGLRC)::Cint
@createOpenGLFun wglDestroyPbufferEXT(hPbuffer::HPBUFFEREXT)::Cint
@createOpenGLFun wglGetVideoDeviceNV(hDC::HDC, numDevices::Cint, hVideoDevice::HPVIDEODEV)::Cint
@createOpenGLFun wglDXOpenDeviceNV(dxDevice::Ptr{Void})::Ptr{Void}
@createOpenGLFun wglShareLists(hrcSrvShare::HGLRC, hrcSrvSource::HGLRC)::Cint
@createOpenGLFun wglWaitForSbcOML(hdc::HDC, target_sbc::Clonglong, ust::Ptr{Clonglong}, msc::Ptr{Clonglong}, sbc::Ptr{Clonglong})::Cint
@createOpenGLFun wglRestoreBufferRegionARB(hRegion::HANDLE, x::Cint, y::Cint, width::Cint, height::Cint, xSrc::Cint, ySrc::Cint)::Cint
@createOpenGLFun wglGetGenlockSourceEdgeI3D(hDC::HDC, uEdge::Ptr{Cuint})::Cint
@createOpenGLFun wglCreateDisplayColorTableEXT(id::GLushort)::Bool
@createOpenGLFun wglGetMscRateOML(hdc::HDC, numerator::Ptr{Cint}, denominator::Ptr{Cint})::Cint
@createOpenGLFun wglGenlockSourceDelayI3D(hDC::HDC, uDelay::Cuint)::Cint
@createOpenGLFun wglSwapIntervalEXT(interval::Cint)::Cint
@createOpenGLFun wglCreateBufferRegionARB(hDC::HDC, iLayerPlane::Cint, uType::Cuint)::Ptr{Void}
@createOpenGLFun wglAssociateImageBufferEventsI3D(hDC::HDC, pEvent::HANDLE, pAddress::LPVOID, pSize::Ptr{DWORD}, count::Cuint)::Cint
@createOpenGLFun wglGetGenlockSampleRateI3D(hDC::HDC, uRate::Ptr{Cuint})::Cint
@createOpenGLFun wglGetPixelFormat(hdc::HDC)::Cint
@createOpenGLFun wglSetGammaTableI3D(hDC::HDC, iEntries::Cint, puRed::Ptr{Cushort}, puGreen::Ptr{Cushort}, puBlue::Ptr{Cushort})::Cint
@createOpenGLFun wglDestroyPbufferARB(hPbuffer::HPBUFFERARB)::Cint
@createOpenGLFun wglGetGammaTableParametersI3D(hDC::HDC, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglDXCloseDeviceNV(hDevice::HANDLE)::Cint
@createOpenGLFun wglGetCurrentContext()::Ptr{Void}
@createOpenGLFun wglQuerySwapGroupNV(hDC::HDC, group::Ptr{GLuint}, barrier::Ptr{GLuint})::Cint
@createOpenGLFun wglGetPixelFormatAttribfvARB(hdc::HDC, iPixelFormat::Cint, iLayerPlane::Cint, nAttributes::Cuint, piAttributes::Ptr{Cint}, pfValues::Ptr{Cfloat})::Cint
@createOpenGLFun wglDXUnregisterObjectNV(hDevice::HANDLE, hObject::HANDLE)::Cint
@createOpenGLFun wglChoosePixelFormatARB(hdc::HDC, piAttribIList::Ptr{Cint}, pfAttribFList::Ptr{Cfloat}, nMaxFormats::Cuint, piFormats::Ptr{Cint}, nNumFormats::Ptr{Cuint})::Cint
@createOpenGLFun wglGenlockSourceI3D(hDC::HDC, uSource::Cuint)::Cint
@createOpenGLFun wglCreateAffinityDCNV(phGpuList::HGPUNV)::Ptr{Void}
@createOpenGLFun wglJoinSwapGroupNV(hDC::HDC, group::GLuint)::Cint
@createOpenGLFun wglMakeCurrent(hDc::HDC, newContext::HGLRC)::Cint
@createOpenGLFun wglGetGPUInfoAMD(id::Cuint, property::Cint, dataType::GLenum, size::Cuint, data::Ptr{Void})::Cint
@createOpenGLFun wglCreateAssociatedContextAMD(id::Cuint)::Ptr{Void}
@createOpenGLFun wglDestroyImageBufferI3D(hDC::HDC, pAddress::LPVOID)::Cint
@createOpenGLFun wglDXLockObjectsNV(hDevice::HANDLE, count::GLint, hObjects::HANDLE)::Cint
@createOpenGLFun wglEndFrameTrackingI3D()::Cint
@createOpenGLFun wglQueryFrameCountNV(hDC::HDC, count::Ptr{GLuint})::Cint
@createOpenGLFun wglGetSyncValuesOML(hdc::HDC, ust::Ptr{Clonglong}, msc::Ptr{Clonglong}, sbc::Ptr{Clonglong})::Cint
@createOpenGLFun wglGetCurrentAssociatedContextAMD()::Ptr{Void}
@createOpenGLFun wglDisableFrameLockI3D()::Cint
@createOpenGLFun wglLockVideoCaptureDeviceNV(hDc::HDC, hDevice::HVIDEOINPUTDEVICENV)::Cint
@createOpenGLFun wglCopyContext(hglrcSrc::HGLRC, hglrcDst::HGLRC, mask::Cuint)::Cint
@createOpenGLFun wglBindVideoDeviceNV(hDC::HDC, uVideoSlot::Cuint, hVideoDevice::HVIDEOOUTPUTDEVICENV, piAttribList::Ptr{Cint})::Cint
@createOpenGLFun wglSetPixelFormat(hdc::HDC, ipfd::Cint, ppfd::PIXELFORMATDESCRIPTOR)::Cint
@createOpenGLFun wglReleaseTexImageARB(hPbuffer::HPBUFFERARB, iBuffer::Cint)::Cint
@createOpenGLFun wglEnableFrameLockI3D()::Cint
@createOpenGLFun wglSetDigitalVideoParametersI3D(hDC::HDC, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglChoosePixelFormat(hDc::HDC, pPfd::PIXELFORMATDESCRIPTOR)::Cint
@createOpenGLFun wglGetPbufferDCEXT(hPbuffer::HPBUFFEREXT)::Ptr{Void}
@createOpenGLFun wglReleasePbufferDCEXT(hPbuffer::HPBUFFEREXT, hDC::HDC)::Cint
@createOpenGLFun wglChoosePixelFormatEXT(hdc::HDC, piAttribIList::Ptr{Cint}, pfAttribFList::Ptr{Cfloat}, nMaxFormats::Cuint, piFormats::Ptr{Cint}, nNumFormats::Ptr{Cuint})::Cint
@createOpenGLFun wglGenlockSourceEdgeI3D(hDC::HDC, uEdge::Cuint)::Cint
@createOpenGLFun wglDXRegisterObjectNV(hDevice::HANDLE, dxObject::Ptr{Void}, name::GLuint, type_::GLenum, access::GLenum)::Ptr{Void}
@createOpenGLFun wglSaveBufferRegionARB(hRegion::HANDLE, x::Cint, y::Cint, width::Cint, height::Cint)::Cint
@createOpenGLFun wglAllocateMemoryNV(size::GLsizei, readfreq::GLfloat, writefreq::GLfloat, priority::GLfloat)::Ptr{Void}
@createOpenGLFun wglGetCurrentDC()::Ptr{Void}
@createOpenGLFun wglGetProcAddress(lpszProc::LPCSTR)::Ptr{Void}
@createOpenGLFun wglBindTexImageARB(hPbuffer::HPBUFFERARB, iBuffer::Cint)::Cint
@createOpenGLFun wglReleaseImageBufferEventsI3D(hDC::HDC, pAddress::LPVOID, count::Cuint)::Cint
@createOpenGLFun wglSwapLayerBuffers(hdc::HDC, fuFlags::Cuint)::Cint
@createOpenGLFun wglGetExtensionsStringARB(hdc::HDC)::Ptr{Cchar}
@createOpenGLFun wglGetExtensionsStringEXT()::Ptr{Cchar}
@createOpenGLFun wglCreateLayerContext(hDc::HDC, level::Cint)::Ptr{Void}
@createOpenGLFun wglDXSetResourceShareHandleNV(dxObject::Ptr{Void}, shareHandle::HANDLE)::Cint
@createOpenGLFun wglBindVideoImageNV(hVideoDevice::HPVIDEODEV, hPbuffer::HPBUFFERARB, iVideoBuffer::Cint)::Cint
@createOpenGLFun wglQueryCurrentContextNV(iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglSetStereoEmitterState3DL(hDC::HDC, uState::Cuint)::Cint
@createOpenGLFun wglReleaseVideoImageNV(hPbuffer::HPBUFFERARB, iVideoBuffer::Cint)::Cint
@createOpenGLFun wglGetVideoInfoNV(hpVideoDevice::HPVIDEODEV, pulCounterOutputPbuffer::Ptr{Culong}, pulCounterOutputVideo::Ptr{Culong})::Cint
@createOpenGLFun wglQueryPbufferARB(hPbuffer::HPBUFFERARB, iAttribute::Cint, piValue::Ptr{Cint})::Cint
@createOpenGLFun wglDisableGenlockI3D(hDC::HDC)::Cint
@createOpenGLFun wglCreateImageBufferI3D(hDC::HDC, dwSize::DWORD, uFlags::Cuint)::Ptr{Void}
@createOpenGLFun wglQueryMaxSwapGroupsNV(hDC::HDC, maxGroups::Ptr{GLuint}, maxBarriers::Ptr{GLuint})::Cint
@createOpenGLFun wglBindDisplayColorTableEXT(id::GLushort)::Bool
@createOpenGLFun wglGenlockSampleRateI3D(hDC::HDC, uRate::Cuint)::Cint
@createOpenGLFun wglGetContextGPUIDAMD(hglrc::HGLRC)::Cuint
@createOpenGLFun wglCreatePbufferARB(hDC::HDC, iPixelFormat::Cint, iWidth::Cint, iHeight::Cint, piAttribList::Ptr{Cint})::Ptr{Void}
@createOpenGLFun wglMakeContextCurrentARB(hDrawDC::HDC, hReadDC::HDC, hglrc::HGLRC)::Cint
@createOpenGLFun wglSwapLayerBuffersMscOML(hdc::HDC, fuPlanes::Cint, target_msc::Clonglong, divisor::Clonglong, remainder::Clonglong)::CClonglong
@createOpenGLFun wglSwapBuffersMscOML(hdc::HDC, target_msc::Clonglong, divisor::Clonglong, remainder::Clonglong)::CClonglong
@createOpenGLFun wglDeleteContext(oldContext::HGLRC)::Cint
@createOpenGLFun wglQueryFrameLockMasterI3D(pFlag::Ptr{Bool})::Cint
@createOpenGLFun wglDescribeLayerPlane(hDc::HDC, pixelFormat::Cint, layerPlane::Cint, nBytes::Cuint, plpd::LAYERPLANEDESCRIPTOR)::Cint
@createOpenGLFun wglGetGenlockSourceDelayI3D(hDC::HDC, uDelay::Ptr{Cuint})::Cint
@createOpenGLFun wglQueryGenlockMaxSourceDelayI3D(hDC::HDC, uMaxLineDelay::Ptr{Cuint}, uMaxPixelDelay::Ptr{Cuint})::Cint
@createOpenGLFun wglGetGammaTableI3D(hDC::HDC, iEntries::Cint, puRed::Ptr{Cushort}, puGreen::Ptr{Cushort}, puBlue::Ptr{Cushort})::Cint
@createOpenGLFun wglEnumerateVideoCaptureDevicesNV(hDc::HDC, phDeviceList::HVIDEOINPUTDEVICENV)::Cuint
@createOpenGLFun wglReleasePbufferDCARB(hPbuffer::HPBUFFERARB, hDC::HDC)::Cint
@createOpenGLFun wglEnumerateVideoDevicesNV(hDC::HDC, phDeviceList::HVIDEOOUTPUTDEVICENV)::Cint
@createOpenGLFun wglGetLayerPaletteEntries(hdc::HDC, iLayerPlane::Cint, iStart::Cint, cEntries::Cint, pcr::COLORREF)::Cint
@createOpenGLFun wglDeleteDCNV(hdc::HDC)::Cint
@createOpenGLFun wglCreatePbufferEXT(hDC::HDC, iPixelFormat::Cint, iWidth::Cint, iHeight::Cint, piAttribList::Ptr{Cint})::Ptr{Void}
@createOpenGLFun wglSetLayerPaletteEntries(hdc::HDC, iLayerPlane::Cint, iStart::Cint, cEntries::Cint, pcr::COLORREF)::Cint
@createOpenGLFun wglResetFrameCountNV(hDC::HDC)::Cint
#constants
const GL_WGL_MIPMAP_LEVEL_ARB                                         = 0x207B
const GL_WGL_GENLOCK_SOURCE_MULTIVIEW_I3D                             = 0x2044
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_RGB_NV                         = 0x20A0
const GL_WGL_DIGITAL_VIDEO_GAMMA_CORRECTED_I3D                        = 0x2053
const GL_WGL_AUX8_ARB                                                 = 0x208F
const GL_WGL_SWAP_LAYER_BUFFERS_EXT                                   = 0x2006
const GL_WGL_BIND_TO_TEXTURE_RGBA_ARB                                 = 0x2071
const GL_WGL_BIND_TO_VIDEO_RGB_AND_DEPTH_NV                           = 0x20C2
const GL_WGL_TYPE_RGBA_EXT                                            = 0x202B
const GL_WGL_AUX2_ARB                                                 = 0x2089
const GL_WGL_GENLOCK_SOURCE_DIGITAL_FIELD_I3D                         = 0x2049
const GL_WGL_AUX5_ARB                                                 = 0x208C
const GL_WGL_SAMPLE_BUFFERS_3DFX                                      = 0x2060
const GL_WGL_BLUE_SHIFT_ARB                                           = 0x201A
const GL_WGL_LOSE_CONTEXT_ON_RESET_ARB                                = 0x8252
const GL_WGL_STENCIL_BITS_EXT                                         = 0x2023
const GL_WGL_MAX_PBUFFER_WIDTH_ARB                                    = 0x202F
const GL_WGL_TRANSPARENT_EXT                                          = 0x200A
const GL_WGL_CONTEXT_RESET_NOTIFICATION_STRATEGY_ARB                  = 0x8256
const GL_WGL_GPU_NUM_PIPES_AMD                                        = 0x21A5
const GL_WGL_SWAP_EXCHANGE_ARB                                        = 0x2028
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_R_NV                     = 0x20B1
const GL_WGL_GPU_RENDERER_STRING_AMD                                  = 0x1F01
const GL_WGL_COVERAGE_SAMPLES_NV                                      = 0x2042
const GL_WGL_ACCUM_GREEN_BITS_ARB                                     = 0x201F
const GL_WGL_CONTEXT_ROBUST_ACCESS_BIT_ARB                            = 0x00000004
const GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_Y_ARB                          = 0x207F
const GL_WGL_GENLOCK_SOURCE_EXTERNAL_TTL_I3D                          = 0x2047
const GL_WGL_BIND_TO_VIDEO_RGB_NV                                     = 0x20C0
const GL_WGL_SHARE_DEPTH_ARB                                          = 0x200C
const GL_WGL_PBUFFER_LOST_ARB                                         = 0x2036
const GL_WGL_TYPE_COLORINDEX_ARB                                      = 0x202C
const GL_WGL_SAMPLE_BUFFERS_EXT                                       = 0x2041
const GL_WGL_GENLOCK_SOURCE_EDGE_FALLING_I3D                          = 0x204A
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_RGBA_NV                        = 0x20A1
const GL_WGL_ACCUM_ALPHA_BITS_EXT                                     = 0x2021
const GL_WGL_NUMBER_UNDERLAYS_EXT                                     = 0x2009
const GL_WGL_VIDEO_OUT_STACKED_FIELDS_2_1                             = 0x20CC
const GL_WGL_SAMPLES_3DFX                                             = 0x2061
const GL_WGL_GREEN_BITS_ARB                                           = 0x2017
const GL_WGL_AUX7_ARB                                                 = 0x208E
const GL_WGL_DEPTH_BITS_EXT                                           = 0x2022
const GL_WGL_DRAW_TO_BITMAP_EXT                                       = 0x2002
const GL_WGL_CONTEXT_MINOR_VERSION_ARB                                = 0x2092
const GL_WGL_SUPPORT_OPENGL_EXT                                       = 0x2010
const GL_WGL_TYPE_RGBA_ARB                                            = 0x202B
const GL_WGL_AUX9_ARB                                                 = 0x2090
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RGBA_NV                  = 0x20B4
const GL_WGL_RED_BITS_EXT                                             = 0x2015
const GL_WGL_SUPPORT_GDI_EXT                                          = 0x200F
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RGB_NV                   = 0x20B3
const GL_WGL_GAMMA_EXCLUDE_DESKTOP_I3D                                = 0x204F
const GL_WGL_BACK_COLOR_BUFFER_BIT_ARB                                = 0x00000002
const GL_WGL_TEXTURE_FLOAT_R_NV                                       = 0x20B5
const GL_WGL_GENERIC_ACCELERATION_EXT                                 = 0x2026
const GL_WGL_IMAGE_BUFFER_LOCK_I3D                                    = 0x00000002
const GL_WGL_NUMBER_PIXEL_FORMATS_ARB                                 = 0x2000
const GL_WGL_ALPHA_SHIFT_EXT                                          = 0x201C
const GL_WGL_RED_BITS_ARB                                             = 0x2015
const GL_WGL_SAMPLES_EXT                                              = 0x2042
const GL_WGL_GREEN_SHIFT_EXT                                          = 0x2018
const GL_WGL_ACCESS_READ_ONLY_NV                                      = 0x00000000
const GL_WGL_TEXTURE_RGBA_ARB                                         = 0x2076
const GL_WGL_CUBE_MAP_FACE_ARB                                        = 0x207C
const GL_WGL_NO_RESET_NOTIFICATION_ARB                                = 0x8261
const GL_WGL_AUX1_ARB                                                 = 0x2088
const GL_WGL_ALPHA_BITS_ARB                                           = 0x201B
const GL_WGL_STEREO_EMITTER_ENABLE_3DL                                = 0x2055
const GL_WGL_ACCUM_BITS_ARB                                           = 0x201D
const GL_WGL_STEREO_POLARITY_NORMAL_3DL                               = 0x2057
const GL_WGL_CONTEXT_MAJOR_VERSION_ARB                                = 0x2091
const GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_Z_ARB                          = 0x2081
const GL_WGL_COLOR_SAMPLES_NV                                         = 0x20B9
const GL_WGL_STENCIL_BUFFER_BIT_ARB                                   = 0x00000008
const GL_WGL_TRANSPARENT_VALUE_EXT                                    = 0x200B
const GL_WGL_GREEN_SHIFT_ARB                                          = 0x2018
const GL_WGL_SWAP_UNDEFINED_ARB                                       = 0x202A
const GL_WGL_DEPTH_TEXTURE_FORMAT_NV                                  = 0x20A5
const GL_WGL_NEED_PALETTE_EXT                                         = 0x2004
const GL_WGL_DRAW_TO_WINDOW_ARB                                       = 0x2001
const GL_WGL_AUX3_ARB                                                 = 0x208A
const GL_WGL_FULL_ACCELERATION_EXT                                    = 0x2027
const GL_WGL_ERROR_INCOMPATIBLE_AFFINITY_MASKS_NV                     = 0x20D0
const GL_WGL_RED_SHIFT_ARB                                            = 0x2016
const GL_WGL_MAX_PBUFFER_HEIGHT_EXT                                   = 0x2030
const GL_WGL_COLOR_BITS_EXT                                           = 0x2014
const GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_X_ARB                          = 0x207D
const GL_WGL_NUM_VIDEO_CAPTURE_SLOTS_NV                               = 0x20CF
const GL_WGL_GENERIC_ACCELERATION_ARB                                 = 0x2026
const GL_WGL_BACK_RIGHT_ARB                                           = 0x2086
const GL_WGL_DEPTH_COMPONENT_NV                                       = 0x20A7
const GL_WGL_OPTIMAL_PBUFFER_WIDTH_EXT                                = 0x2031
const GL_WGL_ACCUM_ALPHA_BITS_ARB                                     = 0x2021
const GL_WGL_VIDEO_OUT_FIELD_2                                        = 0x20CA
const GL_ERROR_INCOMPATIBLE_DEVICE_CONTEXTS_ARB                       = 0x2054
const GL_WGL_DIGITAL_VIDEO_CURSOR_INCLUDED_I3D                        = 0x2052
const GL_WGL_AUX_BUFFERS_EXT                                          = 0x2024
const GL_WGL_DRAW_TO_WINDOW_EXT                                       = 0x2001
const GL_WGL_FLOAT_COMPONENTS_NV                                      = 0x20B0
const GL_WGL_ACCUM_BLUE_BITS_ARB                                      = 0x2020
const GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_Z_ARB                          = 0x2082
const GL_WGL_BACK_LEFT_ARB                                            = 0x2085
const GL_WGL_GREEN_BITS_EXT                                           = 0x2017
const GL_WGL_SWAP_LAYER_BUFFERS_ARB                                   = 0x2006
const GL_WGL_VIDEO_OUT_COLOR_NV                                       = 0x20C3
const GL_WGL_CONTEXT_FORWARD_COMPATIBLE_BIT_ARB                       = 0x00000002
const GL_WGL_STEREO_EXT                                               = 0x2012
const GL_WGL_PBUFFER_WIDTH_EXT                                        = 0x2034
const GL_WGL_DRAW_TO_PBUFFER_ARB                                      = 0x202D
const GL_WGL_DIGITAL_VIDEO_CURSOR_ALPHA_VALUE_I3D                     = 0x2051
const GL_WGL_CONTEXT_LAYER_PLANE_ARB                                  = 0x2093
const GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_X_ARB                          = 0x207E
const GL_WGL_TEXTURE_DEPTH_COMPONENT_NV                               = 0x20A6
const GL_WGL_NUMBER_PIXEL_FORMATS_EXT                                 = 0x2000
const GL_WGL_SWAP_COPY_EXT                                            = 0x2029
const GL_WGL_DEPTH_FLOAT_EXT                                          = 0x2040
const GL_WGL_GENLOCK_SOURCE_EDGE_RISING_I3D                           = 0x204B
const GL_WGL_TEXTURE_FORMAT_ARB                                       = 0x2072
const GL_WGL_VIDEO_OUT_COLOR_AND_ALPHA_NV                             = 0x20C6
const GL_WGL_BLUE_BITS_EXT                                            = 0x2019
const GL_WGL_MIPMAP_TEXTURE_ARB                                       = 0x2074
const GL_WGL_VIDEO_OUT_STACKED_FIELDS_1_2                             = 0x20CB
const GL_WGL_GPU_OPENGL_VERSION_STRING_AMD                            = 0x1F02
const GL_WGL_TEXTURE_1D_ARB                                           = 0x2079
const GL_WGL_FRONT_LEFT_ARB                                           = 0x2083
const GL_WGL_NUMBER_OVERLAYS_ARB                                      = 0x2008
const GL_WGL_DOUBLE_BUFFER_EXT                                        = 0x2011
const GL_WGL_DRAW_TO_BITMAP_ARB                                       = 0x2002
const GL_WGL_SHARE_STENCIL_ARB                                        = 0x200D
const GL_WGL_BLUE_BITS_ARB                                            = 0x2019
const GL_WGL_AUX4_ARB                                                 = 0x208B
const GL_WGL_BIND_TO_TEXTURE_RGB_ARB                                  = 0x2070
const GL_WGL_VIDEO_OUT_FIELD_1                                        = 0x20C9
const GL_WGL_MAX_PBUFFER_PIXELS_ARB                                   = 0x202E
const GL_WGL_MAX_PBUFFER_HEIGHT_ARB                                   = 0x2030
const GL_WGL_TYPE_COLORINDEX_EXT                                      = 0x202C
const GL_WGL_NO_ACCELERATION_ARB                                      = 0x2025
const GL_WGL_NUM_VIDEO_SLOTS_NV                                       = 0x20F0
const GL_WGL_STEREO_POLARITY_INVERT_3DL                               = 0x2058
const GL_WGL_NUMBER_OVERLAYS_EXT                                      = 0x2008
const GL_WGL_NEED_PALETTE_ARB                                         = 0x2004
const GL_WGL_STEREO_ARB                                               = 0x2012
const GL_WGL_SWAP_EXCHANGE_EXT                                        = 0x2028
const GL_WGL_SHARE_ACCUM_EXT                                          = 0x200E
const GL_WGL_FULL_ACCELERATION_ARB                                    = 0x2027
const GL_WGL_ACCELERATION_EXT                                         = 0x2003
const GL_WGL_SWAP_COPY_ARB                                            = 0x2029
const GL_WGL_NO_ACCELERATION_EXT                                      = 0x2025
const GL_WGL_OPTIMAL_PBUFFER_HEIGHT_EXT                               = 0x2032
const GL_WGL_PBUFFER_HEIGHT_ARB                                       = 0x2035
const GL_WGL_PBUFFER_WIDTH_ARB                                        = 0x2034
const GL_WGL_GPU_NUM_SIMD_AMD                                         = 0x21A6
const GL_WGL_PIXEL_TYPE_ARB                                           = 0x2013
const GL_WGL_BIND_TO_VIDEO_RGBA_NV                                    = 0x20C1
const GL_WGL_MAX_PBUFFER_PIXELS_EXT                                   = 0x202E
const GL_WGL_GPU_NUM_RB_AMD                                           = 0x21A7
const GL_WGL_SWAP_UNDEFINED_EXT                                       = 0x202A
const GL_WGL_GENLOCK_SOURCE_EXTERNAL_SYNC_I3D                         = 0x2045
const GL_WGL_MAX_PBUFFER_WIDTH_EXT                                    = 0x202F
const GL_WGL_TRANSPARENT_ARB                                          = 0x200A
const GL_WGL_CONTEXT_ES2_PROFILE_BIT_EXT                              = 0x00000004
const GL_WGL_ACCUM_RED_BITS_ARB                                       = 0x201E
const GL_WGL_VIDEO_OUT_ALPHA_NV                                       = 0x20C4
const GL_WGL_NEED_SYSTEM_PALETTE_EXT                                  = 0x2005
const GL_WGL_GENLOCK_SOURCE_DIGITAL_SYNC_I3D                          = 0x2048
const GL_WGL_FRONT_RIGHT_ARB                                          = 0x2084
const GL_WGL_AUX_BUFFERS_ARB                                          = 0x2024
const GL_WGL_ALPHA_BITS_EXT                                           = 0x201B
const GL_WGL_NEED_SYSTEM_PALETTE_ARB                                  = 0x2005
const GL_WGL_SWAP_METHOD_ARB                                          = 0x2007
const GL_ERROR_INVALID_PIXEL_TYPE_ARB                                 = 0x2043
const GL_WGL_SHARE_STENCIL_EXT                                        = 0x200D
const GL_WGL_DEPTH_BITS_ARB                                           = 0x2022
const GL_WGL_NO_TEXTURE_ARB                                           = 0x2077
const GL_WGL_BIND_TO_TEXTURE_DEPTH_NV                                 = 0x20A3
const GL_WGL_SAMPLES_ARB                                              = 0x2042
const GL_WGL_TEXTURE_TARGET_ARB                                       = 0x2073
const GL_ERROR_INVALID_PROFILE_ARB                                    = 0x2096
const GL_WGL_GPU_RAM_AMD                                              = 0x21A3
const GL_WGL_CONTEXT_PROFILE_MASK_ARB                                 = 0x9126
const GL_WGL_TEXTURE_RECTANGLE_NV                                     = 0x20A2
const GL_WGL_SAMPLE_BUFFERS_ARB                                       = 0x2041
const GL_WGL_ERROR_MISSING_AFFINITY_MASK_NV                           = 0x20D1
const GL_WGL_GENLOCK_SOURCE_EDGE_BOTH_I3D                             = 0x204C
const GL_WGL_VIDEO_OUT_COLOR_AND_DEPTH_NV                             = 0x20C7
const GL_WGL_GPU_CLOCK_AMD                                            = 0x21A4
const GL_WGL_STEREO_EMITTER_DISABLE_3DL                               = 0x2056
const GL_WGL_ALPHA_SHIFT_ARB                                          = 0x201C
const GL_WGL_FRAMEBUFFER_SRGB_CAPABLE_EXT                             = 0x20A9
const GL_WGL_PBUFFER_HEIGHT_EXT                                       = 0x2035
const GL_WGL_CONTEXT_DEBUG_BIT_ARB                                    = 0x00000001
const GL_WGL_ACCESS_READ_WRITE_NV                                     = 0x00000001
const GL_WGL_TYPE_RGBA_UNSIGNED_FLOAT_EXT                             = 0x20A8
const GL_WGL_BLUE_SHIFT_EXT                                           = 0x201A
const GL_WGL_VIDEO_OUT_DEPTH_NV                                       = 0x20C5
const GL_WGL_ACCELERATION_ARB                                         = 0x2003
const GL_WGL_ACCESS_WRITE_DISCARD_NV                                  = 0x00000002
const GL_WGL_CONTEXT_COMPATIBILITY_PROFILE_BIT_ARB                    = 0x00000002
const GL_WGL_SUPPORT_GDI_ARB                                          = 0x200F
const GL_WGL_ACCUM_GREEN_BITS_EXT                                     = 0x201F
const GL_WGL_TEXTURE_FLOAT_RGB_NV                                     = 0x20B7
const GL_WGL_TEXTURE_FLOAT_RG_NV                                      = 0x20B6
const GL_WGL_FRONT_COLOR_BUFFER_BIT_ARB                               = 0x00000001
const GL_WGL_DRAW_TO_PBUFFER_EXT                                      = 0x202D
const GL_WGL_NUMBER_UNDERLAYS_ARB                                     = 0x2009
const GL_WGL_FRAMEBUFFER_SRGB_CAPABLE_ARB                             = 0x20A9
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RG_NV                    = 0x20B2
const GL_WGL_TYPE_RGBA_FLOAT_ARB                                      = 0x21A0
const GL_WGL_ACCUM_BITS_EXT                                           = 0x201D
const GL_WGL_BIND_TO_TEXTURE_RECTANGLE_DEPTH_NV                       = 0x20A4
const GL_WGL_DEPTH_BUFFER_BIT_ARB                                     = 0x00000004
const GL_WGL_IMAGE_BUFFER_MIN_ACCESS_I3D                              = 0x00000001
const GL_WGL_CONTEXT_CORE_PROFILE_BIT_ARB                             = 0x00000001
const GL_WGL_SHARE_DEPTH_EXT                                          = 0x200C
const GL_WGL_PBUFFER_LARGEST_ARB                                      = 0x2033
const GL_WGL_AUX6_ARB                                                 = 0x208D
const GL_WGL_UNIQUE_ID_NV                                             = 0x20CE
const GL_WGL_TYPE_RGBA_FLOAT_ATI                                      = 0x21A0
const GL_WGL_GPU_NUM_SPI_AMD                                          = 0x21A8
const GL_WGL_ACCUM_RED_BITS_EXT                                       = 0x201E
const GL_WGL_COLOR_BITS_ARB                                           = 0x2014
const GL_WGL_TEXTURE_CUBE_MAP_ARB                                     = 0x2078
const GL_WGL_TEXTURE_RGB_ARB                                          = 0x2075
const GL_WGL_GENLOCK_SOURCE_EXTERNAL_FIELD_I3D                        = 0x2046
const GL_WGL_GPU_VENDOR_AMD                                           = 0x1F00
const GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_Y_ARB                          = 0x2080
const GL_WGL_TEXTURE_FLOAT_RGBA_NV                                    = 0x20B8
const GL_WGL_SHARE_ACCUM_ARB                                          = 0x200E
const GL_WGL_PBUFFER_LARGEST_EXT                                      = 0x2033
const GL_WGL_GPU_FASTEST_TARGET_GPUS_AMD                              = 0x21A2
const GL_ERROR_INVALID_PIXEL_TYPE_EXT                                 = 0x2043
const GL_WGL_DOUBLE_BUFFER_ARB                                        = 0x2011
const GL_WGL_RED_SHIFT_EXT                                            = 0x2016
const GL_WGL_PIXEL_TYPE_EXT                                           = 0x2013
const GL_WGL_TEXTURE_2D_ARB                                           = 0x207A
const GL_WGL_ACCUM_BLUE_BITS_EXT                                      = 0x2020
const GL_WGL_DIGITAL_VIDEO_CURSOR_ALPHA_FRAMEBUFFER_I3D               = 0x2050
const GL_WGL_VIDEO_OUT_FRAME                                          = 0x20C8
const GL_WGL_AUX0_ARB                                                 = 0x2087
const GL_WGL_SWAP_METHOD_EXT                                          = 0x2007
const GL_WGL_CONTEXT_FLAGS_ARB                                        = 0x2094
const GL_WGL_STENCIL_BITS_ARB                                         = 0x2023
const GL_WGL_SUPPORT_OPENGL_ARB                                       = 0x2010
const GL_WGL_GAMMA_TABLE_SIZE_I3D                                     = 0x204E
const GL_ERROR_INVALID_VERSION_ARB                                    = 0x2095
# Export everything!
export wglUseFontBitmapsA, wglGetPbufferDCARB, wglGetPixelFormatAttribfvEXT, wglReleaseVideoDeviceNV, wglBindVideoCaptureDeviceNV, wglGetGenlockSourceI3D, wglGetFrameUsageI3D, wglRealizeLayerPalette, wglEnumGpusNV, wglEnableGenlockI3D, wglGetDigitalVideoParametersI3D, wglCreateContextAttribsARB, wglReleaseVideoCaptureDeviceNV, wglEnumGpusFromAffinityDCNV, wglQueryPbufferEXT, wglQueryVideoCaptureDeviceNV, wglSendPbufferToVideoNV, wglDXObjectAccessNV, wglEnumGpuDevicesNV, wglGetGPUIDsAMD, wglCopyImageSubDataNV, wglGetPixelFormatAttribivARB, wglGetCurrentReadDCARB, wglIsEnabledFrameLockI3D, wglDeleteAssociatedContextAMD, wglDXUnlockObjectsNV, wglUseFontBitmapsW, wglBindSwapBarrierNV, wglGetCurrentReadDCEXT, wglSetPbufferAttribARB, wglQueryFrameTrackingI3D, wglCreateContext, wglWaitForMscOML, wglSetGammaTableParametersI3D, wglIsEnabledGenlockI3D, wglLoadDisplayColorTableEXT, wglSwapBuffers, wglBeginFrameTrackingI3D, wglGetPixelFormatAttribivEXT, wglGetDefaultProcAddress, wglCreateAssociatedContextAttribsAMD, wglMakeContextCurrentEXT, wglGetSwapIntervalEXT, wglMakeAssociatedContextCurrentAMD, wglDestroyPbufferEXT, wglGetVideoDeviceNV, wglDXOpenDeviceNV, wglShareLists, wglWaitForSbcOML, wglRestoreBufferRegionARB, wglGetGenlockSourceEdgeI3D, wglCreateDisplayColorTableEXT, wglGetMscRateOML, wglGenlockSourceDelayI3D, wglSwapIntervalEXT, wglCreateBufferRegionARB, wglAssociateImageBufferEventsI3D, wglGetGenlockSampleRateI3D, wglGetPixelFormat, wglSetGammaTableI3D, wglDestroyPbufferARB, wglGetGammaTableParametersI3D, wglDXCloseDeviceNV, wglGetCurrentContext, wglQuerySwapGroupNV, wglGetPixelFormatAttribfvARB, wglDXUnregisterObjectNV, wglChoosePixelFormatARB, wglGenlockSourceI3D, wglCreateAffinityDCNV, wglJoinSwapGroupNV, wglMakeCurrent, wglGetGPUInfoAMD, wglCreateAssociatedContextAMD, wglDestroyImageBufferI3D, wglDXLockObjectsNV, wglEndFrameTrackingI3D, wglQueryFrameCountNV, wglGetSyncValuesOML, wglGetCurrentAssociatedContextAMD, wglDisableFrameLockI3D, wglLockVideoCaptureDeviceNV, wglCopyContext, wglBindVideoDeviceNV, wglSetPixelFormat, wglReleaseTexImageARB, wglEnableFrameLockI3D, wglSetDigitalVideoParametersI3D, wglChoosePixelFormat, wglGetPbufferDCEXT, wglReleasePbufferDCEXT, wglChoosePixelFormatEXT, wglGenlockSourceEdgeI3D, wglDXRegisterObjectNV, wglSaveBufferRegionARB, wglAllocateMemoryNV, wglGetCurrentDC, wglGetProcAddress, wglBindTexImageARB, wglReleaseImageBufferEventsI3D, wglSwapLayerBuffers, wglGetExtensionsStringARB, wglGetExtensionsStringEXT, wglCreateLayerContext, wglDXSetResourceShareHandleNV, wglBindVideoImageNV, wglQueryCurrentContextNV, wglSetStereoEmitterState3DL, wglReleaseVideoImageNV, wglGetVideoInfoNV, wglQueryPbufferARB, wglDisableGenlockI3D, wglCreateImageBufferI3D, wglQueryMaxSwapGroupsNV, wglBindDisplayColorTableEXT, wglGenlockSampleRateI3D, wglGetContextGPUIDAMD, wglCreatePbufferARB, wglMakeContextCurrentARB, wglSwapLayerBuffersMscOML, wglSwapBuffersMscOML, wglDeleteContext, wglQueryFrameLockMasterI3D, wglDescribeLayerPlane, wglGetGenlockSourceDelayI3D, wglQueryGenlockMaxSourceDelayI3D, wglGetGammaTableI3D, wglEnumerateVideoCaptureDevicesNV, wglReleasePbufferDCARB, wglEnumerateVideoDevicesNV, wglGetLayerPaletteEntries, wglDeleteDCNV, wglCreatePbufferEXT, wglSetLayerPaletteEntries, wglResetFrameCountNV, GL_WGL_MIPMAP_LEVEL_ARB, GL_WGL_GENLOCK_SOURCE_MULTIVIEW_I3D, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_RGB_NV, GL_WGL_DIGITAL_VIDEO_GAMMA_CORRECTED_I3D, GL_WGL_AUX8_ARB, GL_WGL_SWAP_LAYER_BUFFERS_EXT, GL_WGL_BIND_TO_TEXTURE_RGBA_ARB, GL_WGL_BIND_TO_VIDEO_RGB_AND_DEPTH_NV, GL_WGL_TYPE_RGBA_EXT, GL_WGL_AUX2_ARB, GL_WGL_GENLOCK_SOURCE_DIGITAL_FIELD_I3D, GL_WGL_AUX5_ARB, GL_WGL_SAMPLE_BUFFERS_3DFX, GL_WGL_BLUE_SHIFT_ARB, GL_WGL_LOSE_CONTEXT_ON_RESET_ARB, GL_WGL_STENCIL_BITS_EXT, GL_WGL_MAX_PBUFFER_WIDTH_ARB, GL_WGL_TRANSPARENT_EXT, GL_WGL_CONTEXT_RESET_NOTIFICATION_STRATEGY_ARB, GL_WGL_GPU_NUM_PIPES_AMD, GL_WGL_SWAP_EXCHANGE_ARB, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_R_NV, GL_WGL_GPU_RENDERER_STRING_AMD, GL_WGL_COVERAGE_SAMPLES_NV, GL_WGL_ACCUM_GREEN_BITS_ARB, GL_WGL_CONTEXT_ROBUST_ACCESS_BIT_ARB, GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_Y_ARB, GL_WGL_GENLOCK_SOURCE_EXTERNAL_TTL_I3D, GL_WGL_BIND_TO_VIDEO_RGB_NV, GL_WGL_SHARE_DEPTH_ARB, GL_WGL_PBUFFER_LOST_ARB, GL_WGL_TYPE_COLORINDEX_ARB, GL_WGL_SAMPLE_BUFFERS_EXT, GL_WGL_GENLOCK_SOURCE_EDGE_FALLING_I3D, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_RGBA_NV, GL_WGL_ACCUM_ALPHA_BITS_EXT, GL_WGL_NUMBER_UNDERLAYS_EXT, GL_WGL_VIDEO_OUT_STACKED_FIELDS_2_1, GL_WGL_SAMPLES_3DFX, GL_WGL_GREEN_BITS_ARB, GL_WGL_AUX7_ARB, GL_WGL_DEPTH_BITS_EXT, GL_WGL_DRAW_TO_BITMAP_EXT, GL_WGL_CONTEXT_MINOR_VERSION_ARB, GL_WGL_SUPPORT_OPENGL_EXT, GL_WGL_TYPE_RGBA_ARB, GL_WGL_AUX9_ARB, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RGBA_NV, GL_WGL_RED_BITS_EXT, GL_WGL_SUPPORT_GDI_EXT, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RGB_NV, GL_WGL_GAMMA_EXCLUDE_DESKTOP_I3D, GL_WGL_BACK_COLOR_BUFFER_BIT_ARB, GL_WGL_TEXTURE_FLOAT_R_NV, GL_WGL_GENERIC_ACCELERATION_EXT, GL_WGL_IMAGE_BUFFER_LOCK_I3D, GL_WGL_NUMBER_PIXEL_FORMATS_ARB, GL_WGL_ALPHA_SHIFT_EXT, GL_WGL_RED_BITS_ARB, GL_WGL_SAMPLES_EXT, GL_WGL_GREEN_SHIFT_EXT, GL_WGL_ACCESS_READ_ONLY_NV, GL_WGL_TEXTURE_RGBA_ARB, GL_WGL_CUBE_MAP_FACE_ARB, GL_WGL_NO_RESET_NOTIFICATION_ARB, GL_WGL_AUX1_ARB, GL_WGL_ALPHA_BITS_ARB, GL_WGL_STEREO_EMITTER_ENABLE_3DL, GL_WGL_ACCUM_BITS_ARB, GL_WGL_STEREO_POLARITY_NORMAL_3DL, GL_WGL_CONTEXT_MAJOR_VERSION_ARB, GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_Z_ARB, GL_WGL_COLOR_SAMPLES_NV, GL_WGL_STENCIL_BUFFER_BIT_ARB, GL_WGL_TRANSPARENT_VALUE_EXT, GL_WGL_GREEN_SHIFT_ARB, GL_WGL_SWAP_UNDEFINED_ARB, GL_WGL_DEPTH_TEXTURE_FORMAT_NV, GL_WGL_NEED_PALETTE_EXT, GL_WGL_DRAW_TO_WINDOW_ARB, GL_WGL_AUX3_ARB, GL_WGL_FULL_ACCELERATION_EXT, GL_WGL_ERROR_INCOMPATIBLE_AFFINITY_MASKS_NV, GL_WGL_RED_SHIFT_ARB, GL_WGL_MAX_PBUFFER_HEIGHT_EXT, GL_WGL_COLOR_BITS_EXT, GL_WGL_TEXTURE_CUBE_MAP_POSITIVE_X_ARB, GL_WGL_NUM_VIDEO_CAPTURE_SLOTS_NV, GL_WGL_GENERIC_ACCELERATION_ARB, GL_WGL_BACK_RIGHT_ARB, GL_WGL_DEPTH_COMPONENT_NV, GL_WGL_OPTIMAL_PBUFFER_WIDTH_EXT, GL_WGL_ACCUM_ALPHA_BITS_ARB, GL_WGL_VIDEO_OUT_FIELD_2, GL_ERROR_INCOMPATIBLE_DEVICE_CONTEXTS_ARB, GL_WGL_DIGITAL_VIDEO_CURSOR_INCLUDED_I3D, GL_WGL_AUX_BUFFERS_EXT, GL_WGL_DRAW_TO_WINDOW_EXT, GL_WGL_FLOAT_COMPONENTS_NV, GL_WGL_ACCUM_BLUE_BITS_ARB, GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_Z_ARB, GL_WGL_BACK_LEFT_ARB, GL_WGL_GREEN_BITS_EXT, GL_WGL_SWAP_LAYER_BUFFERS_ARB, GL_WGL_VIDEO_OUT_COLOR_NV, GL_WGL_CONTEXT_FORWARD_COMPATIBLE_BIT_ARB, GL_WGL_STEREO_EXT, GL_WGL_PBUFFER_WIDTH_EXT, GL_WGL_DRAW_TO_PBUFFER_ARB, GL_WGL_DIGITAL_VIDEO_CURSOR_ALPHA_VALUE_I3D, GL_WGL_CONTEXT_LAYER_PLANE_ARB, GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_X_ARB, GL_WGL_TEXTURE_DEPTH_COMPONENT_NV, GL_WGL_NUMBER_PIXEL_FORMATS_EXT, GL_WGL_SWAP_COPY_EXT, GL_WGL_DEPTH_FLOAT_EXT, GL_WGL_GENLOCK_SOURCE_EDGE_RISING_I3D, GL_WGL_TEXTURE_FORMAT_ARB, GL_WGL_VIDEO_OUT_COLOR_AND_ALPHA_NV, GL_WGL_BLUE_BITS_EXT, GL_WGL_MIPMAP_TEXTURE_ARB, GL_WGL_VIDEO_OUT_STACKED_FIELDS_1_2, GL_WGL_GPU_OPENGL_VERSION_STRING_AMD, GL_WGL_TEXTURE_1D_ARB, GL_WGL_FRONT_LEFT_ARB, GL_WGL_NUMBER_OVERLAYS_ARB, GL_WGL_DOUBLE_BUFFER_EXT, GL_WGL_DRAW_TO_BITMAP_ARB, GL_WGL_SHARE_STENCIL_ARB, GL_WGL_BLUE_BITS_ARB, GL_WGL_AUX4_ARB, GL_WGL_BIND_TO_TEXTURE_RGB_ARB, GL_WGL_VIDEO_OUT_FIELD_1, GL_WGL_MAX_PBUFFER_PIXELS_ARB, GL_WGL_MAX_PBUFFER_HEIGHT_ARB, GL_WGL_TYPE_COLORINDEX_EXT, GL_WGL_NO_ACCELERATION_ARB, GL_WGL_NUM_VIDEO_SLOTS_NV, GL_WGL_STEREO_POLARITY_INVERT_3DL, GL_WGL_NUMBER_OVERLAYS_EXT, GL_WGL_NEED_PALETTE_ARB, GL_WGL_STEREO_ARB, GL_WGL_SWAP_EXCHANGE_EXT, GL_WGL_SHARE_ACCUM_EXT, GL_WGL_FULL_ACCELERATION_ARB, GL_WGL_ACCELERATION_EXT, GL_WGL_SWAP_COPY_ARB, GL_WGL_NO_ACCELERATION_EXT, GL_WGL_OPTIMAL_PBUFFER_HEIGHT_EXT, GL_WGL_PBUFFER_HEIGHT_ARB, GL_WGL_PBUFFER_WIDTH_ARB, GL_WGL_GPU_NUM_SIMD_AMD, GL_WGL_PIXEL_TYPE_ARB, GL_WGL_BIND_TO_VIDEO_RGBA_NV, GL_WGL_MAX_PBUFFER_PIXELS_EXT, GL_WGL_GPU_NUM_RB_AMD, GL_WGL_SWAP_UNDEFINED_EXT, GL_WGL_GENLOCK_SOURCE_EXTERNAL_SYNC_I3D, GL_WGL_MAX_PBUFFER_WIDTH_EXT, GL_WGL_TRANSPARENT_ARB, GL_WGL_CONTEXT_ES2_PROFILE_BIT_EXT, GL_WGL_ACCUM_RED_BITS_ARB, GL_WGL_VIDEO_OUT_ALPHA_NV, GL_WGL_NEED_SYSTEM_PALETTE_EXT, GL_WGL_GENLOCK_SOURCE_DIGITAL_SYNC_I3D, GL_WGL_FRONT_RIGHT_ARB, GL_WGL_AUX_BUFFERS_ARB, GL_WGL_ALPHA_BITS_EXT, GL_WGL_NEED_SYSTEM_PALETTE_ARB, GL_WGL_SWAP_METHOD_ARB, GL_ERROR_INVALID_PIXEL_TYPE_ARB, GL_WGL_SHARE_STENCIL_EXT, GL_WGL_DEPTH_BITS_ARB, GL_WGL_NO_TEXTURE_ARB, GL_WGL_BIND_TO_TEXTURE_DEPTH_NV, GL_WGL_SAMPLES_ARB, GL_WGL_TEXTURE_TARGET_ARB, GL_ERROR_INVALID_PROFILE_ARB, GL_WGL_GPU_RAM_AMD, GL_WGL_CONTEXT_PROFILE_MASK_ARB, GL_WGL_TEXTURE_RECTANGLE_NV, GL_WGL_SAMPLE_BUFFERS_ARB, GL_WGL_ERROR_MISSING_AFFINITY_MASK_NV, GL_WGL_GENLOCK_SOURCE_EDGE_BOTH_I3D, GL_WGL_VIDEO_OUT_COLOR_AND_DEPTH_NV, GL_WGL_GPU_CLOCK_AMD, GL_WGL_STEREO_EMITTER_DISABLE_3DL, GL_WGL_ALPHA_SHIFT_ARB, GL_WGL_FRAMEBUFFER_SRGB_CAPABLE_EXT, GL_WGL_PBUFFER_HEIGHT_EXT, GL_WGL_CONTEXT_DEBUG_BIT_ARB, GL_WGL_ACCESS_READ_WRITE_NV, GL_WGL_TYPE_RGBA_UNSIGNED_FLOAT_EXT, GL_WGL_BLUE_SHIFT_EXT, GL_WGL_VIDEO_OUT_DEPTH_NV, GL_WGL_ACCELERATION_ARB, GL_WGL_ACCESS_WRITE_DISCARD_NV, GL_WGL_CONTEXT_COMPATIBILITY_PROFILE_BIT_ARB, GL_WGL_SUPPORT_GDI_ARB, GL_WGL_ACCUM_GREEN_BITS_EXT, GL_WGL_TEXTURE_FLOAT_RGB_NV, GL_WGL_TEXTURE_FLOAT_RG_NV, GL_WGL_FRONT_COLOR_BUFFER_BIT_ARB, GL_WGL_DRAW_TO_PBUFFER_EXT, GL_WGL_NUMBER_UNDERLAYS_ARB, GL_WGL_FRAMEBUFFER_SRGB_CAPABLE_ARB, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_FLOAT_RG_NV, GL_WGL_TYPE_RGBA_FLOAT_ARB, GL_WGL_ACCUM_BITS_EXT, GL_WGL_BIND_TO_TEXTURE_RECTANGLE_DEPTH_NV, GL_WGL_DEPTH_BUFFER_BIT_ARB, GL_WGL_IMAGE_BUFFER_MIN_ACCESS_I3D, GL_WGL_CONTEXT_CORE_PROFILE_BIT_ARB, GL_WGL_SHARE_DEPTH_EXT, GL_WGL_PBUFFER_LARGEST_ARB, GL_WGL_AUX6_ARB, GL_WGL_UNIQUE_ID_NV, GL_WGL_TYPE_RGBA_FLOAT_ATI, GL_WGL_GPU_NUM_SPI_AMD, GL_WGL_ACCUM_RED_BITS_EXT, GL_WGL_COLOR_BITS_ARB, GL_WGL_TEXTURE_CUBE_MAP_ARB, GL_WGL_TEXTURE_RGB_ARB, GL_WGL_GENLOCK_SOURCE_EXTERNAL_FIELD_I3D, GL_WGL_GPU_VENDOR_AMD, GL_WGL_TEXTURE_CUBE_MAP_NEGATIVE_Y_ARB, GL_WGL_TEXTURE_FLOAT_RGBA_NV, GL_WGL_SHARE_ACCUM_ARB, GL_WGL_PBUFFER_LARGEST_EXT, GL_WGL_GPU_FASTEST_TARGET_GPUS_AMD, GL_ERROR_INVALID_PIXEL_TYPE_EXT, GL_WGL_DOUBLE_BUFFER_ARB, GL_WGL_RED_SHIFT_EXT, GL_WGL_PIXEL_TYPE_EXT, GL_WGL_TEXTURE_2D_ARB, GL_WGL_ACCUM_BLUE_BITS_EXT, GL_WGL_DIGITAL_VIDEO_CURSOR_ALPHA_FRAMEBUFFER_I3D, GL_WGL_VIDEO_OUT_FRAME, GL_WGL_AUX0_ARB, GL_WGL_SWAP_METHOD_EXT, GL_WGL_CONTEXT_FLAGS_ARB, GL_WGL_STENCIL_BITS_ARB, GL_WGL_SUPPORT_OPENGL_ARB, GL_WGL_GAMMA_TABLE_SIZE_I3D, GL_ERROR_INVALID_VERSION_ARB


text:
function double_integrator_dynamics!(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T
    ẋ[1] = x[2]
    ẋ[2] = u[1]
end

n = 2
m = 1

doubleintegrator = Model(double_integrator_dynamics!,n,m)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using Singular

include("DifferentialPolynomials.jl")

# Should not be here?
AbstractAlgebra.needs_parentheses(x::Singular.n_unknown{Nemo.fmpq}) = false
AbstractAlgebra.needs_parentheses(x::Singular.n_unknown{AbstractAlgebra.Generic.Frac{Nemo.fmpq_mpoly}}) = true

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

function lex_ranking(varnames::Array{String, 1})
    return (x, y) -> begin
        vx, ordx = decompose_derivative(x)
        vy, ordy = decompose_derivative(y)
        indx = findfirst(a -> a == vx, varnames)
        indy = findfirst(a -> a == vy, varnames)
        if indx != indy
            return indx < indy
        end
        return ordx < ordy
    end
end

function lex_ranking(R::DiffPolyRing)
    return lex_ranking(R.varnames)
end

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

function get_singular_polys(diff_polys::Array{DiffPoly, 1}, ordering=:lex, ranking=nothing)
    R = parent(diff_polys[1])
    if ranking === nothing
        ranking = lex_ranking(R)
    end

    orders = Array{Integer, 1}()
    for v in gens(R)
        push!(orders, max([order(p, v) for p in diff_polys]...))
    end
    newvarnames = Array{String, 1}()
    for (i, name) in enumerate(R.varnames)
        for ord in 0:orders[i]
            push!(newvarnames, form_derivative(name, ord))
        end
    end
    sort!(newvarnames, lt=ranking)
    bring = base_ring(R.poly_ring)
    Rsing, varsing = Singular.PolynomialRing(bring, newvarnames, ordering=ordering)
    return Rsing, [parent_ring_change(algdata(p), Rsing) for p in diff_polys]
end

function from_singular(p, R::DiffPolyRing)
    return DiffPoly(R, parent_ring_change(p, R.poly_ring))
end

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

max_stars_repo_name: pogudingleb/DifferentialAlgebra

max_issues_repo_name: pogudingleb/DifferentialAlgebra

max_forks_repo_name: pogudingleb/DifferentialAlgebra


text:
# This implements the common parts of Evolutionary Programming, i.e. the
# population-based "culling" between generations. Different variants of EP
# can then be implemented by changing the mutation-operator etc.

PopulationBasedEvolutionaryProgrammingDefaultParameters = {
  :PopulationSize => 50,
  :InitialStandardDeviation => 3.0,  # Initial standard deviation
}

# The specific functions of a particular EP strategy is attached to a strategy
# datum while the generic population-based EP is kept in a single type. This
# way we can easily change the strategy without having to rewrite the common parts
# throughout.
abstract EvolutionaryProgrammingStrategy

type PopulationBasedEvolutionaryProgramming <: PopulationOptimizer
  parameters::Parameters

  generation::Int64
  dimension::Int64
  population_size::Int64

  # Pre-calc the tau values so we do not have to recalc it again and again during evolution.
  tau::Float64
  tauprim::Float64

  # Two populations that we switch between as being the current population.
  # The generation mod(number, 2) determines which one is current. This way
  # we can avoid continues memory allocation and freeing of the population
  # structures.
  pop0::Array{Float64, 2}
  pop1::Array{Float64, 2}

  # Child population used in creating the next population
  children::Array{Float64, 2}

  # And the sigma values for each individual are kept in their separate matrices.
  # Exact same structure as for the populations and children matrices.
  sigmas0::Array{Float64, 2}
  sigmas1::Array{Float64, 2}
  childsigmas::Array{Float64, 2}

  PopulationBasedEvolutionaryProgramming(parameters) = begin

    dim = parameters[:NumDimensions]
    popsize = parameters[:PopulationSize]
    initial_population = rand_individuals_lhs(parameters[:SearchSpace], popsize)
    initial_sigmas = parameters[:InitialStandardDeviation] * ones(dims, popsize)

    new(parameters, 0, dims, popsize,
      1.0 / sqrt(2 * sqrt(popsize)), 1.0 / sqrt(2 * popsize),
      initial_population, zeros(Float64, dims, popsize), zeros(Float64, dims, popsize),
      initial_sigmas, zeros(Float64, dims, popsize), zeros(Float64, dims, popsize))

  end
end

function get_current_population(pep::PopulationBasedEvolutionaryProgramming)
  if mod(pep.generation, 2) == 0
    return pep.pop0, pep.sigmas0
  else
    return pep.pop1, pep.sigmas1
  end
end

# The generic EP strategy functions have the abstract type and thus specific
# strategies can fall back on them simply by not "overriding" that function.
# This implements a classic Gaussian EP with gaussian mutations (GEP).
has_ask_tell_interface(ep::EvolutionaryProgrammingStrategy) = false

function step(ep::EvolutionaryProgrammingStrategy)

  pep = ep.pep
  children = pep.children
  childsigmas = pep.childsigmas
  currentpop, currentsigmas, nextpop, nextsigmas = get_current_and_next_population_and_sigma(pep)

  # Mutate each parent to create a child.
  for i in 1:pep.population_size
    children[:,i], childsigmas[:,i] = mutate_individual(ep, currentpop, currentsigmas,
      i, pep.dimension, pep.tauprim, pep.tau)
  end

  # Select the individuals to survive and copy them to the next population.
  select_winning_individuals(ep, currentpop, currentsigmas, children, childsigmas, nextpop, nextsigmas)

end

function select_winning_individuals(ep::EvolutionaryProgrammingStrategy,
  currentpop::Array{Float64, 2}, currentsigmas::Array{Float64, 2},
  children::Array{Float64, 2}, childsigmas::Array{Float64, 2},
  nextpop::Array{Float64, 2}, nextsigmas::Array{Float64, 2})

end

function mutate_individual(ep::EvolutionaryProgrammingStrategy,
  population::Array{Float64, 2}, sigmas::Array{Float64, 2},
  i::Int64, dims::Int64, tauprim::Float64, tau::Float64)

  tauprimconst = tauprim * sample_tauprimconst_distribution(ep)

  childsigmas = sigmas[:,i] .* exp( tauprimconst + tau * sample_sigma_distribution(ep, dims)' )
  child = population[:,i] .+ childsigmas .* sample_mutation_distribution(ep, dims)'

  return child, childsigmas

end

function sample_tauprimconst_distribution(ep::EvolutionaryProgrammingStrategy)
  randn() # Standard is to use a N(0, 1) gaussian...
end

# The standard distribution is a Gaussian(0, 1)
function sample_mutation_distribution(ep::EvolutionaryProgrammingStrategy, numsamples::Int64)
  randn(numsamples) # Standard is to use a N(0, 1) gaussians...
end

# The standard distribution for mutating the sigmas is a Gaussian(0, 1) for each dimension.
function sample_sigma_distribution(ep::EvolutionaryProgrammingStrategy, numsamples::Int64)
  randn(numsamples) # Standard is to use a N(0, 1) gaussians...
end

# Specific EP variants are implemented in their own types and have the PBEP
# structure in them

type GaussianEvolutionaryProgramming <: EvolutionaryProgrammingStrategy
  pep::PopulationBasedEvolutionaryProgramming
end


text:
function Δt(ts::AbstractVector{<:Date})

    @assert length(ts) >= 3 "time series must have a length >= 3"

    ts1 = ts[1]
    ts2 = ts[2]
    ts3 = ts[3]

    dy = year(ts[2])-year(ts[1])
    dm = month(ts[2])-month(ts[1])
    dd = day(ts[2])-day(ts[1])

    same_day = day(ts1) == day(ts2) == day(ts3)
    same_month = month(ts1) == month(ts2) == month(ts3)
    same_year = year(ts1) == year(ts2) == year(ts3)

    !same_day && dd % 7 == 0 && return Week(dd÷7)
    !same_day && return Day(dd)
    same_day && !same_month && dm % 4 == 0 && return Quarter(dm÷4)
    same_day && !same_month && return Month(dm)
    same_day && !same_year && return Year(dy)
    return ts[2]-ts[1]

end

function Δt(ts::AbstractVector{<:DateTime})

    @assert length(ts) >= 3 "time series must have a length >= 3"

    ts1 = ts[1]
    ts2 = ts[2]
    ts3 = ts[3]

    dt = Dates.value(ts2-ts1) # Milliseconds

    dt % (24*3_600_000) == 0 && return Δt(convert(AbstractVector{Date},ts))
    dt % 3_600_000 == 0 && return Hour(dt/3_600_000)
    dt % 60_000 == 0 && return Minute(dt/60_000)
    dt % 1_000 == 0 && return Second(dt/1_000)
    return ts[2]-ts[1]

end

function Δt(ts::AbstractVector{<:Time})

    @assert length(ts) >= 3 "time series must have a length >= 3"

    ts1 = ts[1]
    ts2 = ts[2]
    ts3 = ts[3]

    dt = Dates.value(ts2-ts1) # Nanoseconds

    dt % 3_600_000_000_000 == 0 && return Hour(dt/3_600_000_000_000)
    dt % 60_000_000_000 == 0 && return Minute(dt/60_000_000_000)
    dt % 1_000_000_000 == 0 && return Second(dt/1_000_000_000)
    dt % 1_000_000 == 0 && return Millisecond(dt/1_000_000)
    dt % 1_000 == 0  && return Microsecond(dt/1_000)
    return Nanosecond(dt)

end

Δt(ts::AbstractVector) = ts[2]-ts[1]

# Next n increments for a timestamp vector
function nΔt(ts::AbstractVector, n::Integer)
    dt = Δt(ts)
    ts1 = n > 0 ? ts[end] + dt : ts[1] - dt
    n > 0 && return reshape(collect(ts1 : dt : ts1 + (n-1)*dt),:,1)
    n < 0 && return reshape(collect(ts1 + (n+1)*dt: dt : ts1),:,1)
    return reshape(ts,:,1)
end

# DataFrame to canon Time Series Dataframe
function tots(df::DataFrame, interval = Day; start::Integer = 1)::DataFrame

    x = df[:,eltype.(eachcol(df)) .<: Union{Missing,Real}]
    names_x = names(x)
    ts = df[:,eltype.(eachcol(df)) .<: Union{Date, DateTime, Time,
                                             Year, Quarter, Month, Day,
                                             Hour, Second, Millisecond,
                                             Microsecond, Nanosecond}]
    if size(ts) != (0,0)
        names_x = vcat(names(ts)[:,1],names_x)
        return df[!,names_x]
    else
        ts = reshape(collect(interval(start):interval(1):interval(start-1+size(x,1))),:,1)
        tdf = DataFrame(ts,["Time"])
        return hcat(tdf,df[!,names_x])
    end

end

function tots(ar::Matrix, interval = Day; start::Integer = 1)::DataFrame
    df = DataFrame(ar,:auto)
    tots(df,interval; start)
end

function tots(ar::Vector, interval = Day; start::Integer = 1)::DataFrame
    df = DataFrame(reshape(ar,:,1),:auto)
    tots(df,interval; start)
end


text:
using Spotify
using Spotify.JSON3
using Dates: now
!endswith(pwd(), "test") && cd("test")
"""
- 1-2   Albums
- 3-6   Artists
- 8-13  Browse
- 14    Episodes
- 15-21 Follow
- 22-33 Library
- 34-35 Personalization
- 36-37 Player
- 38-40 Playlists
- Search
- Shows
- Tracks
- UsersProfile
- Objects
"""
verb_vec = include("verb_vec.jl")
tuples_vec = include("parameter_tuples_vec.jl")
param_dic = include("parameter_default_dic.jl")
get_defaults_tuple(tu) = map(x -> param_dic[x], tu)
get_defaults_tuple(tu::Symbol) = (param_dic[tu],)

function get_default_call_expression(n)
    verb = verb_vec[n]
    tup = tuples_vec[n]
    defatup = get_defaults_tuple(tup)
    Expr(:call, verb, defatup...)
end

"""
    make_default_call(n)

Make Spotify API call number n, ref. 'verb_vec'.
"""
make_default_call(n) = eval(get_default_call_expression(n))

"""
    show_default_call(n)
    show_default_call(io::IO, n)

Show Spotify API call number n, ref. 'verb_vec'.
"""
function show_default_call(n)
    show(stderr, get_default_call_expression(n))
    print(stderr, "\n")
end
show_default_call(io::IO, n) = show(io, get_default_call_expression(n))