text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 9af1fa20-07fa-11eb-2d8d-239d4cc3f3ae
using SSpline

# ╔═╡ 2de193c2-080c-11eb-0451-5d5c46e11324
using Colors, Plots

# ╔═╡ 7b664047-b676-4a20-9c8c-5e610497f9fb
begin
	#import DarkMode
    #DarkMode.enable(theme="material-palenight", cm_config=Dict("tabSize" => 4))
	#DarkMode.enable()
    #DarkMode.Toolbox(theme="default")
end

# ╔═╡ 0c4701a0-dd1c-4f4b-91eb-cb94bc9d4440
md"cell width"

# ╔═╡ 679c9e8a-ef4f-4b2b-a1ce-245b0ffa773d
html"<style> main {max-width: 950px;}"

# ╔═╡ 3eb09f10-63f2-11eb-38c5-db7e4547b561
pyplot()#gr()

# ╔═╡ 8890e470-63f0-11eb-3df9-19f6bddcf0fd
md"""
Create x and y vectors of original points:
"""

# ╔═╡ 26bca8a0-080c-11eb-33b2-7bf3e34cf10a
begin
	x=[10,20,35,50]
	y=[10,30,20,22]
	#x=[0,10,30,50,70,90,100]
	#y=[30,130,150,150,170,220,320]
	#x=[0.0,25.0,50.0,75.0,100.0]
	#y=[20.0,50.0,10.0,60.0,35.0]
	#x=[0,10,30,45,70,90,100]
	#y=[20,50,70,-70,60,50,20]
	#x=[0,10,30,50,70,80,100]
	#y=[150,200,200,200,180,100,0]
end

# ╔═╡ 9240d97e-63eb-11eb-0a53-278a141ee9a6
md"""Create nearest neighbor spline:

`nearestneighbor(x, y)`, `stepspline(x, y)`, `spline0(x,y)`

Result:

`NearestNeighborSpline(a, x)`
"""

# ╔═╡ 39f7e0b0-63ec-11eb-0a32-eb2d79a226d5
nns = nearestneighbor(x,y)

# ╔═╡ 2f3b8ef0-5920-11eb-3e40-39ee8b03512a
md"""Create linear spline:

`linearspline(x, y)`, `lspline(x, y)`, `spline1(x,y)`

Result:

`LinearSpline(b, a, x)`
"""

# ╔═╡ 61ab5d39-f5e4-46ea-8fc8-2be9bec862b3
md"### Test: `spline1` with extrapolation type `:zero`, `:boundary`, `:linear`"

# ╔═╡ 6d4484bb-7924-4a93-b225-d163484d0315
ls = spline1(x,y)

# ╔═╡ ef1bf6db-2352-4d40-9888-70d93ff471f0
begin
	plot(legend=:outertopright,size=(850,400))

	lz = extrap(ls,-10:110,:zero)
	lzi = interp(lz,-10:1:110)

	scatter!(collect(lz.x[1:end-1]),lz.a,markersize=10,color=:orange,label="knots :zero")
	plot!(lzi[1],lzi[2],line = 2, color=:orange, label="spline1, :zero")

	lb = extrap(ls,-10:110,:boundary)
	lbi = interp(lb,-10:1:110)
	scatter!(x,y,markersize=7,color=:skyblue,label="knots :boundary")
	plot!(lbi[1:end-1],lbi[2],line = 2, style = :dash, color=:skyblue, label="spline1, :boundary")

	ll = extrap(ls,-10:110,:linear)
	lli = interp(ll,-10:1:110)
	scatter!(x,y,markersize=3,color=:grey,label="knots :linear")
	plot!(lli[1:end-1],lli[2],line = 2, style=:dashdot, color=:grey, label="spline1, :linear")
end

# ╔═╡ 8bac622b-eb80-4bcd-a9e4-3cb375089b62
ls,lz

# ╔═╡ 5c543811-e06c-4ba7-aa8f-265c6602a90e
ls,lb

# ╔═╡ 8ca3a6d0-6c49-4fe8-8add-11b2b013a562
ls,ll

# ╔═╡ d13e17c8-693a-4f94-9b3f-cf4621f43484
md"### Test: `spline3` with extrapolation type `:zero`, `:boundary`, `:linear`"

# ╔═╡ a525192c-5641-4058-9e34-73c0452285b3
	cs = spline3(x,y,:natural)

# ╔═╡ 5d7b2ff7-8b94-4144-9206-5d07e5d60425
begin
	cz = extrapolate(cs,-10:110,:zero)
	czi = interp(cz,-10:1:110)
	scatter(x,y,markersize=5,color=:white,legend = :bottom,label="knots")
	plot!(czi[1:end-1],czi[2],linewidth=3, linestyle=:dot,color=:orange, label="spline3, :zero")
	cb = extrapolate(cs,-10:110,:boundary)
	cbi = interp(cb,-10:1:110)
	plot!(cbi[1:end-1],cbi[2],linewidth = 3, linestyle=:dash,color=:cyan3, label="spline3, :boundary")
	cl = extrapolate(cs,-10:110,:linear)
	cli = interp(cl,-10:1:110)
	plot!(cli[1],cli[2],linewidth = 2, linestyle=:dashdot,color=:red, label="spline3, :linear")
end

# ╔═╡ 7f6d5943-673e-4619-91b4-2ec66515cbc6
cs,cz

# ╔═╡ 58e5f56b-0230-4cab-a463-aa72e783d48b
cs,cb

# ╔═╡ 1fc24cbb-d237-469d-af7b-5a975ae584a4
cs,cl

# ╔═╡ cbd7235f-76b7-451b-ad2f-83f1a5f1089c
function extra(Spl,xrange,extr::Symbol)
    envx=collect(xrange)
    xs=0.0 # start of extrapolatin range
    xe=0.0 # end of extrapolation range
    envx[1] < Spl.x[1] ? xs=envx[1] : xs=Spl.x[1]
    envx[1] < Spl.x[1] ? xs=envx[1] : xs=Spl.x[1]
    envx[end] > Spl.x[end] ? xe=envx[end] : xe=Spl.x[end]

    x=copy(Spl.x)
    a=copy(Spl.a)
	nf=nextfloat(Float64(x[end]))

    if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
        b=copy(Spl.b)
    end
    if typeof(Spl) in (QuadraticSpline, CubicSpline)
        c=copy(Spl.c)
    end
    if typeof(Spl) == CubicSpline
        d=copy(Spl.d)
    end
	if extr == :zero
		if xs < x[1]
			pushfirst!(x,xs)
			pushfirst!(a, 0.0)
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				x[end]=nf
				push!(x,xe)
				push!(a,0.0)
				push!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				push!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				push!(d,0.0)
			end
		end
	elseif extr == :boundary
		de=x[end]-x[end-1]
		if xs < x[1]
			pushfirst!(x,xs)
			pushfirst!(a, a[1])
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,0.0)
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) == LinearSpline
				push!(a,Spl.b[end]*de+Spl.a[end])
				push!(b,0.0)
				x[end]=nf
			end
			if typeof(Spl) == QuadraticSpline
				push!(a,c[end]*de^2 + b[end]*de + a[end])
				push!(b,0.0)
				push!(c,0.0)
				x[end]=nf
			end
			if typeof(Spl) == CubicSpline
				push!(a,d[end]*de^3 + c[end]*de^2 + b[end]*de + a[end])
				push!(b,0.0)
				push!(c,0.0)
				push!(d,0.0)
				x[end]=nf
			end
		end
		push!(x,xe)
	elseif extr == :linear
		if xs < x[1]
			pushfirst!(x,xs)
			#
			pushfirst!(a,a[1]+b[1]*(x[1]-x[2]))
			if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				pushfirst!(b,b[1])
			end
			if typeof(Spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				pushfirst!(d,0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(Spl) == LinearSpline
				push!(a,a[end]+b[end]*(x[end]-x[end-1]))
				push!(b,b[end])
			end
			if typeof(Spl) == QuadraticSpline
				push!(c,0.0)
			end
			if typeof(Spl) == CubicSpline
				de=x[end]-x[end-1]
				push!(a,d[end]*de^3 + c[end]*de^2 + b[end]*de + a[end])
				push!(b,3*d[end]*de^2 + 2*c[end]*de + b[end])
				push!(c,0.0)
				push!(d,0.0)
				x[end]=nf
			end
		end
		push!(x,xe)
    elseif extr == :quadratic
        if xs < x[1]
            pushfirst!(x,xs)
			pushfirst!(a,a[1]/(b[1]*(x[1]-(x[1]-x[2]))))
            if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                #pushfirst!(b,x[1]*((a[3]-a[2])+x[2]*(a[1]-b[3])+x[3]*(a[2]-a[1]))/((x[1]-x[2])*(x[1]-x[3])*(x[2]-x[3])))
				pushfirst!(b,x[2]*((a[3]-a[2])+x[3]*(a[1]-b[3])+x[4]*(a[2]-a[1]))/((x[2]-x[3])*(x[2]-x[4])*(x[3]-x[4])))
            end
            if typeof(Spl) in (QuadraticSpline, CubicSpline)
                pushfirst!(c,0.0)
            end
            if typeof(Spl) == CubicSpline
                pushfirst!(d,0.0)
            end
        end
        if xe ≥ x[end]
            push!(x,xe)
            if typeof(Spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                push!(b,x[end-2]*((a[end]-a[end-1])+x[end-1]*(a[end-2]-b[end])+x[end]*(a[end-1]-a[end-2]))/((x[end-2]-x[end-1])*(x[end-2]-x[end])*(x[end-1]-x[end])))
                push!(a,a[end]/(b[end]*(x[end]-xs)))
            end
            if typeof(Spl) in (QuadraticSpline, CubicSpline)
                push!(c,0.0)
            end
            if typeof(Spl) == CubicSpline
                push!(d,0.0)
            end
        end
    else error("Extrapolation type $extr does not exist.")
    end
    if typeof(Spl) == NearestNeighborSpline
        NearestNeighborSpline(a,x)
    elseif typeof(Spl) == LinearSpline
        LinearSpline(b,a,x)
    elseif typeof(Spl) == QuadraticSpline
        QuadraticSpline(c,b,a,x)
    else
        CubicSpline(d,c,b,a,x)
    end
end

# ╔═╡ 62c1dde5-6a6f-4404-9890-e14c614d7c25

# ╔═╡ 778ac396-d942-4bf8-b6b1-2aacfb3c0d8b

# ╔═╡ 5376ead3-58be-4887-83e8-24f95802a0d9
begin
	scatter(x,y,markersize=5,color=:white,legend = :bottom,label="knots")
	cq = extrap(cs,-10:110,:linear)
	cqi = interp(cq,-10:1:110)
	plot!(cqi[1],cqi[2],linewidth = 2, linestyle=:dash,color=:red, label="spline3, :quadratic")
end

# ╔═╡ daf3d3fd-8033-401a-b106-115a65479136

# ╔═╡ 028b3dc0-63f1-11eb-2b34-ab517656c8d3
md"""Create cubic splines:

* `cubicspline(x::Vector,y::Vector,style...)`
   * shorthands for `cubicspline(...)`:
     * `cspline(...)`
     * `spline3(...)`

`style...` =

* `:natural`
  * natural spline, zero curvature at both boundaries
    * shorthand: `cubicspline(x::Vector,y::Vector)`

* `:clamped, ds, de`
  * clamped spline, arbitrary slopes at start (`ds`) and end (`de`) boundaries
* `:periodic`
   * periodic spline, identical values, slope and curvature at both boundaries.
* `:constrained` (in the works)
   * constrained spline, no overshooting of the interpolation function. zero slope at spline values if at least one slope is zero or the slope has a a sign change.

Result:

object of type `CubicSpline` with the fields `d, c, b, a, x`
"""

# ╔═╡ a520a3e0-63f1-11eb-3d5c-31ad2e72d89c
md"spline interpolations:"

# ╔═╡ 73c8b780-38a7-11eb-247b-f3e0fe47db66
function cons(x,y)
    length(x)!=length(y) ? error("Vector length mismatch. Length(x): $(length(x)) != length(y): $(length(y))") : nothing
    length(x) < 3 ? error("Too few points: $(length(x)). Cubic splines need at least 3 points") : nothing

    n = length(x)-1
    Δx= zeros(n+1)
    Δy= zeros(n+1)
    σ = zeros(n+1) #slope

    @inbounds for i in 1:n-1 # 1:n-1
        Δx[i] = x[i+1] - x[i] #h_i = (x_i+1 - x_i)
        Δy[i] = y[i+1] - y[i]
    end

    @inbounds for i in 2:n-1
        slope = Δy[i-1]*Δy[i]
        if slope > 0.0
			σ[i] = 2.0/(Δx[i+1]/Δy[i] + Δx[i-1]/Δy[i-1])
		elseif slope <= 0.0
			σ[i] = 0.0
		end
    end

    σ[1] = 3.0* Δy[1]  /(2.0* Δx[1])   - σ[2]/2.0
    σ[n] = 3.0* Δy[n-1]/(2.0* Δx[n-1]) - σ[n-1]/2.0

    a = zeros(n+1)
    b = zeros(n+1)
    c = zeros(n+1)
    d = zeros(n+1)

    @inbounds for i in 2:n
        #κa =  -2.0*(σ[i+1]+2*σ[i])/Δx[i] + 6*Δy[i]/Δx[i]^2 # f''[i-1]
        #κb =  2.0*(2.0*σ[i+1]+σ[i])/Δx[i] - 6*Δy[i]/Δx[i]^2 # f''[i]
        κa = 2.0 * (σ[i] + 2* σ[i-1]) / Δx[i-1] + 6 * Δy[i-1] / (Δx[i-1]*Δx[i-1])
        κb = 2.0* (2.0 * σ[i] + σ[i-1]) / Δx[i-1] - 6 * Δy[i-1] / (Δx[i-1] * Δx[i-1])
        d[i]=(κb-κa)/(6*Δx[i-1])
        c[i]=(x[i]*κa - x[i-1]*κb)/(2*Δx[i-1])
        b[i]=(Δy[i-1]- c[i]*(x[i]^2 - x[i-1]^2) - d[i]*(x[i]^3 - x[i-1]^3))/Δx[i-1]
        a[i] = y[i-1] - b[i]*x[i-1] - c[i]*x[i-1]^2 - d[i]*x[i-1]^3
    end
    CubicSpline(d,c,b,a,x)
end

# ╔═╡ 017de170-5d17-11eb-2964-3b6996cb7eb8
begin
	csnat=cubicspline(x,y,:natural)
	csclamp=cubicspline(x,y,:clamped,-10,10)
	cstrain=cons(x,y)
	#csconstr=conspline(x,y)
end

# ╔═╡ 6b1bc1ee-59c6-11eb-3a5c-2d0c5ab35d0c
begin
	nni=interp(nns,collect(0:1:100))
	lint=interp(ls,collect(0:100))
	cintn=interp(csnat,collect(0:100))
	cintc=interp(csclamp,collect(0:100))
	cints=interp(cstrain,collect(0:100))
end

# ╔═╡ 51f3c956-9462-490d-8375-a5001509a9d5
begin
	plot(legend=:outertopright, size=(850,450))
	plot!(nni[1],nni[2],color=:red, style=:dot, label="nearest neighbor")
	plot!(lint[1],lint[2], color=:orange1, label="linear")
	plot!(cintn[1],cintn[2],color=:green4,style=:dash, label="cubic natural")
	#plot!(cintc[1],cintc[2],color=:blue2,style=:dashdot, label="cubic clamped, slope at endpoints = 0, 0")
	plot!(cints[1],cints[2],color=:violet, label="cubic constrained")
	#plot!(cintcon[1],cintcon[2],color=:black,style=:dash,label="cubic constrained")
	scatter!(x,y,label="knots",color=:black, markersize=3)
end

# ╔═╡ 55e3217a-e35e-496e-9f8d-2e042392bb68
function ocons(x,y)
    length(x)!=length(y) ? error("Vector length mismatch. Length(x): $(length(x)) != length(y): $(length(y))") : nothing
    length(x) < 3 ? error("Too few points: $(length(x)). Cubic splines need at least 3 points") : nothing

    n = length(x)
    Δx= zeros(n+1)
    Δy= zeros(n+1)
    σ = zeros(n+1) #slope

    @inbounds for i in 1:n-1 # 1:n-1
        Δx[i] = x[i+1] - x[i] #h_i = (x_i+1 - x_i)
        Δy[i] = y[i+1] - y[i]
    end

    @inbounds for i in 2:n-1
        slope = Δy[i-1]*Δy[i]
        if slope > 0.0
			σ[i] = 2.0/(Δx[i+1]/Δy[i] + Δx[i-1]/Δy[i-1])
		elseif slope <= 0.0
			σ[i] = 0.0
		end
    end

    σ[1] = 3.0* Δy[1]  /(2.0* Δx[1])   - σ[2]/2.0
    σ[n] = 3.0* Δy[n-1]/(2.0* Δx[n-1]) - σ[n-1]/2.0

    a = zeros(n+1)
    b = zeros(n+1)
    c = zeros(n+1)
    d = zeros(n+1)

    @inbounds for i in 2:n
        #κa =  -2.0*(σ[i+1]+2*σ[i])/Δx[i] + 6*Δy[i]/Δx[i]^2 # f''[i-1]
        #κb =  2.0*(2.0*σ[i+1]+σ[i])/Δx[i] - 6*Δy[i]/Δx[i]^2 # f''[i]
        κa = -2.0 * (σ[i] + 2* σ[i-1]) / Δx[i-1] + 6 * Δy[i-1] / Δx[i-1]^2
        κb = 2.0* (2.0 * σ[i] + σ[i-1]) / Δx[i-1] - 6 * Δy[i-1] / Δx[i-1]^2
        a[i-1]=(κb-κa)/(6*Δx[i-1])
        b[i-1]=(x[i]*κa - x[i-1]*κb)/(2*Δx[i-1])
        c[i-1]=(Δy[i-1]- b[i]*(x[i]^2 - x[i-1]^2) - a[i]*(x[i]^3 - x[i-1]^3))/Δx[i-1]
        d[i-1] = y[i-1] - c[i]*x[i-1] - b[i]*x[i-1]^2 - a[i]*x[i-1]^3
    end
    CubicSpline(a,b,c,d,x)
end

# ╔═╡ 50b92ab0-5920-11eb-3f85-01c690d30a78
lspline(x,y)

# ╔═╡ c17b386b-d17f-4930-bd4f-bf5b4c0b302a
function constrained(x,y)

end

# ╔═╡ 013f0c19-c116-48b3-99af-e9ca7a7974bc
x

# ╔═╡ 6f8fadd2-9157-452f-b66e-7d798077356f
y

# ╔═╡ 7893260e-b68e-470a-b42d-6f9563bd3407
constrainedspline(x,y)

# ╔═╡ 1433e532-2f95-11eb-238b-8d3969489ce0
function conspline(x,y)
    length(x)!=length(y) ? error("Vector length mismatch. Length(x): $(length(x)) != length(y): $(length(y))") : nothing
    length(x) < 3 ? error("Too few points: $(length(x)). Spline needs at least 3 points") : nothing

    n = length(x)#-1
    global Δx= zeros(n)
    Δy= zeros(n)
    σ= zeros(n+1) #slope
    a = zeros(n+1)
    b = zeros(n+1)
    c = zeros(n+1)
    d = zeros(n+1)
    for i in 1:n-1
        Δx[i] = x[i+1] - x[i] #h_i = (x_i+1 - x_i)
        Δy[i] = y[i+1] - y[i]
    end

    @inbounds for i in 1:n-1
        slopechange = y[i]*y[i-1]
        if slopechange > 0.0
			σ[i] = 2.0/(Δx[i]/Δy[i] + Δx[i-1]/Δy[i-1])
		elseif slopechange <= 0.0
			σ[i] = 0.0
		end
    end
	σ[1] = 1.5 * (y[2]-y[1]) / (x[2]-x[1]) - σ[2] / 2.0
    σ[n] = 1.5 * (y[n]-y[n-1]) / (x[n]-x[n-1]) - σ[n-1] / 2.0

	κi_1=0.0
	κi=0.0

	global κineg=zeros(n)
	global κipos=zeros(n)
	κineg[1]=0.0
	κipos[1]=0.0
    for i in 2:n
        κi_1 = -2*(σ[i]+2*σ[i-1]) / (x[i]-x[i-1]) + 6*(y[i]-y[i-1])/(x[i]-x[i-1])^2 # f''[i-1]
		#κineg[i]=κi_1
        κi =  2*(2*σ[i]+σ[i-1])/ (x[i]-x[i-1]) - 6*(y[i]-y[i-1])/(x[i]-x[i-1])^2 # f''[i]
		#κipos[i]=κi
        #d[i]=(κi-κi_1)/(6*(x[i]-x[i-1]))
		#c[i]=(x[i]*κi_1 - x[i-1]*κi)/(2*(x[i]-x[i-1]))
		b[i]=((y[i]-y[i-1])-c[i]*(x[i]^2-x[i-1]^2)-d[i]*(x[i]^3-x[i-1]^3))/(x[i]-x[i-1])
		a[i] = y[i-1] - b[i]*x[i-1] - c[i]*x[i-1]^2 - d[i]*x[i-1]^3
    end
	CubicSpline(d[2:end],c[2:end],b[2:end],vcat(a[2:end-1],y[end]),x)
end

# ╔═╡ 5793d204-3531-43c3-aab6-67f9d68b0a1f
ocons(x,y)

# ╔═╡ bdf29d20-49c9-11eb-3993-393f8ddf9ef6

# ╔═╡ c509f3b0-38b2-11eb-0388-a30708ab7cc8
md"Splines"

# ╔═╡ ca6d58ae-38b2-11eb-3fc1-bf442497c9df
md"Interpolations"

# ╔═╡ 8322c440-38b2-11eb-0d19-27627b9a99c0
md"Graphs"

# ╔═╡ 4ce80fc2-2ff2-11eb-3ef4-8315ae577413
md"""

Spline degrees:

* linearspline, lspline, spline1

* quadraticspline, qspline, spline2

* cubicspline, cspline, spline3

* cubicsplineip

Spline styles:

* constrainedspline, naturalspline, clampedspline, periodicspline

Interpolations and extrapolations:

* interp, interptest, extrap, extrapolate

Derivatives:

* slope, deriv1

* curvature, deriv2

* curvrate, deriv3

Types:

* CubicSpline, QuadraticSpline, LinearSpline
"""

# ╔═╡ 5ceaf4c0-87e3-11eb-0821-733526b5a838
x

# ╔═╡ c2cf3461-052d-4265-a7ca-a354c921a160
y

# ╔═╡ c6a21b54-182f-4fc7-af2d-948c9b1d26f6
begin
	slope=zeros(length(x))
	n=length(x)
	for i in 2:n-1
	if (y[i]-y[i-1])*(y[i+1]-y[i]) ≤ 0.0
		slope[i]=0.0
	else
		slope[i]=2.0/((x[i+1]-x[i])/(y[i+1]-y[i]) + (x[i]-x[i-1])/(y[i]-y[i-1]))
	end
	slope[1] = 1.5*(y[2]-y[1])/(x[2]-x[1]) - slope[2]/2.0
    slope[n] = 1.5*(y[n]-y[n-1])/(x[n]-x[n-1]) - slope[n-1]/2.0
	end
	slope
end

# ╔═╡ ca80905a-5845-4ac3-bcef-068e47f366b5
begin
	curveminus=zeros(n)
	curveplus=zeros(n)
	for i in 2:n-1
		cm = 2.0*(slope[i]+2*slope[i-1])/(x[i]-x[i-1]) + 6*(y[i])-y[i-1]/(x[i]-x[i-1])^2
		cp = 2.0*(2.0*slope[i]+slope[i-1])/(x[i]-x[i-1]) - 6*(y[i])-y[i-1]/(x[i]-x[i-1])^2
		curveminus[i]=cm
		curveplus[i]=cp
	end
	curveminus,curveplus
end

# ╔═╡ 04111612-7c8f-47c8-a9a3-c20dcd693ae1
begin
	d=zeros(n)
	c=zeros(n)
	b=zeros(n)
	a=zeros(n)
	for i in 2:n
	cm = -2*(slope[i]+2*slope[i-1])/(x[i]-x[i-1]) + 6*(y[i]-y[i-1])/(x[i]-x[i-1])^2
	cp =2*(2*slope[i]+slope[i-1])/(x[i]-x[i-1]) - 6*(y[i]-y[i-1])/(x[i]-x[i-1])^2
		d[i-1]=(cp-cm)/(6*(x[i]-x[i-1]))
		c[i-1]=(x[i]*cm-x[i-1]*cp)/(2*(x[i]-x[i-1]))
		b[i-1]=((y[i]-y[i-1])-c[i]*(x[i]^2-x[i-1]^2)-d[i]*(x[i]^3-x[i-1]^3))/(x[i]-x[i-1])
		a[i-1]=y[i-1]-b[i]*x[i-1]-c[i]*x[i-1]^2-d[i]*x[i-1]^3
	end
	test=CubicSpline(d[1:n-1],c[1:n-1],b[1:n-1],a[1:n-1],x)
end

# ╔═╡ a455dd2f-bb06-44f4-9988-d0ee0b4df1f3
ctest=interp(test,0:1:100)

# ╔═╡ b25f5b8f-e460-4ab8-80e2-9f49791d6b90
begin
	#plot(ctest[1],ctest[2])
	scatter(x,y,legend=:outertopright, label="original samples")
end

# ╔═╡ Cell order:
# ╠═9af1fa20-07fa-11eb-2d8d-239d4cc3f3ae
# ╠═2de193c2-080c-11eb-0451-5d5c46e11324
# ╠═7b664047-b676-4a20-9c8c-5e610497f9fb
# ╟─0c4701a0-dd1c-4f4b-91eb-cb94bc9d4440
# ╠═679c9e8a-ef4f-4b2b-a1ce-245b0ffa773d
# ╠═3eb09f10-63f2-11eb-38c5-db7e4547b561
# ╟─8890e470-63f0-11eb-3df9-19f6bddcf0fd
# ╠═26bca8a0-080c-11eb-33b2-7bf3e34cf10a
# ╟─9240d97e-63eb-11eb-0a53-278a141ee9a6
# ╠═39f7e0b0-63ec-11eb-0a32-eb2d79a226d5
# ╟─2f3b8ef0-5920-11eb-3e40-39ee8b03512a
# ╟─61ab5d39-f5e4-46ea-8fc8-2be9bec862b3
# ╠═6d4484bb-7924-4a93-b225-d163484d0315
# ╠═ef1bf6db-2352-4d40-9888-70d93ff471f0
# ╠═8bac622b-eb80-4bcd-a9e4-3cb375089b62
# ╠═5c543811-e06c-4ba7-aa8f-265c6602a90e
# ╠═8ca3a6d0-6c49-4fe8-8add-11b2b013a562
# ╟─d13e17c8-693a-4f94-9b3f-cf4621f43484
# ╠═a525192c-5641-4058-9e34-73c0452285b3
# ╠═5d7b2ff7-8b94-4144-9206-5d07e5d60425
# ╠═7f6d5943-673e-4619-91b4-2ec66515cbc6
# ╠═58e5f56b-0230-4cab-a463-aa72e783d48b
# ╠═1fc24cbb-d237-469d-af7b-5a975ae584a4
# ╟─cbd7235f-76b7-451b-ad2f-83f1a5f1089c
# ╠═62c1dde5-6a6f-4404-9890-e14c614d7c25
# ╠═778ac396-d942-4bf8-b6b1-2aacfb3c0d8b
# ╠═5376ead3-58be-4887-83e8-24f95802a0d9
# ╠═daf3d3fd-8033-401a-b106-115a65479136
# ╟─028b3dc0-63f1-11eb-2b34-ab517656c8d3
# ╠═017de170-5d17-11eb-2964-3b6996cb7eb8
# ╠═a520a3e0-63f1-11eb-3d5c-31ad2e72d89c
# ╠═6b1bc1ee-59c6-11eb-3a5c-2d0c5ab35d0c
# ╠═51f3c956-9462-490d-8375-a5001509a9d5
# ╠═73c8b780-38a7-11eb-247b-f3e0fe47db66
# ╟─55e3217a-e35e-496e-9f8d-2e042392bb68
# ╠═50b92ab0-5920-11eb-3f85-01c690d30a78
# ╠═c17b386b-d17f-4930-bd4f-bf5b4c0b302a
# ╠═013f0c19-c116-48b3-99af-e9ca7a7974bc
# ╠═6f8fadd2-9157-452f-b66e-7d798077356f
# ╠═7893260e-b68e-470a-b42d-6f9563bd3407
# ╟─1433e532-2f95-11eb-238b-8d3969489ce0
# ╠═5793d204-3531-43c3-aab6-67f9d68b0a1f
# ╠═bdf29d20-49c9-11eb-3993-393f8ddf9ef6
# ╟─c509f3b0-38b2-11eb-0388-a30708ab7cc8
# ╟─ca6d58ae-38b2-11eb-3fc1-bf442497c9df
# ╟─8322c440-38b2-11eb-0d19-27627b9a99c0
# ╟─4ce80fc2-2ff2-11eb-3ef4-8315ae577413
# ╠═5ceaf4c0-87e3-11eb-0821-733526b5a838
# ╠═c2cf3461-052d-4265-a7ca-a354c921a160
# ╠═c6a21b54-182f-4fc7-af2d-948c9b1d26f6
# ╠═ca80905a-5845-4ac3-bcef-068e47f366b5
# ╠═04111612-7c8f-47c8-a9a3-c20dcd693ae1
# ╠═a455dd2f-bb06-44f4-9988-d0ee0b4df1f3
# ╠═b25f5b8f-e460-4ab8-80e2-9f49791d6b90


text:
# TODO test
function getFileName(cxf::CXFile)
    c = ccall( (:clang_getFileName, libclang), CXString, (CXFile,), cxf)
    !(c.data == C_NULL) ? unsafe_string(c.data) : ""
end

# TODO move this to straight call
function getIncludedFile(cu::InclusionDirective)
    getFileName(ccall( (:clang_getIncludedFile, libclang),
            CXFile,
            (CXCursor,),
            cu))
end

function getNullLocation()
    ccall( (:clang_getNullLocation, libclang), CXSourceLocation, ())
end
function equalLocations(l1::CXSourceLocation, l2::CXSourceLocation)
    ccall( (:clang_equalLocations, libclang), UInt32, (CXSourceLocation,CXSourceLocation), l1,l2 ) == 1
end
function getNullRange()
    ccall( (:clang_getNullRange, libclang), CXSourceRange, ())
end

# returns CXSourceRange
function getRange(l1::CXSourceLocation,l2::CXSourceLocation,)
    ccall( (:clang_getRange, libclang), CXSourceRange, (CXSourceLocation,CXSourceLocation), l1, l2) |>
        CXSourceRange
end
function equalRanges(r1::CXSourceRange,r2::CXSourceRange)
    ccall( (:clang_equalRanges, libclang), UInt32, (CXSourceRange, CXSourceRange), r1, r2) == 1
end
function Range_isNull(cxr::CXSourceRange,)
    ccall( (:clang_Range_isNull, libclang), Int32, (CXSourceRange,), cxr)
end
function getNullCursor()::CLCursor
    ccall( (:clang_getNullCursor, libclang), CXCursor, () )
end
function getTranslationUnitCursor(tu::CXTranslationUnit)::CLCursor
    ccall( (:clang_getTranslationUnitCursor, libclang), CXCursor, (CXTranslationUnit,), tu) |> CLCursor
end
function equalCursors(cu1,cu2)
    ccall( (:clang_equalCursors, libclang), UInt32, (CXCursor, CXCursor), cu1, cu2) == 1
end
function Cursor_isNull(cu)
    ccall( (:clang_Cursor_isNull, libclang), Int32, (CXCursor,), cu) == 1
end
function hashCursor(cu)
    ccall( (:clang_hashCursor, libclang), UInt32, (CXCursor,), cu)
end
function getCursorKind(cu)
    ccall( (:clang_getCursorKind, libclang), UInt32, (CXCursor,), cu)
end

# TODO test
function isDeclaration(k::CXCursorKind)
    ccall( (:clang_isDeclaration, libclang), UInt32, (CXCursorKind,), k) == 1
end
function isReference(k::CXCursorKind)
    ccall( (:clang_isReference, libclang), UInt32, (CXCursorKind,), k) == 1
end
#function getCursorLinkage(cu)
#    ccall( (:clang_getCursorLinkage, libclang), CXLinkageKind, (CXCursor,), cu)
#end TODO
#function getCursorAvailability(cu::CXCursor,)
#    ccall( (:clang_getCursorAvailability, libclang), CXAvailabilityKind, (CXCursor,), cu)
#end

#function getCursorLanguage(cu)
#    ccall( (:clang_getCursorLanguage, libclang), CXLanguageKind, (CXCursor,), cu)
#end
function getCursorSemanticParent(cu)::CLCursor
    ccall( (:clang_getCursorSemanticParent, libclang), CXCursor, (CXCursor,), cu)
end
function getCursorLexicalParent(cu)::CLCursor
    ccall( (:clang_getCursorLexicalParent, libclang), CXCursor, (CXCursor,), cu)
end
function getCursorType(cu)::CLType
    ccall( (:clang_getCursorType, libclang), CXType, (CXCursor,), cu)
end

function getTypedefDeclUnderlyingType(cu)::CLType
    ccall( (:clang_getTypedefDeclUnderlyingType, libclang), CXType, (CXCursor,), cu)
end
function getEnumDeclIntegerType(cu)::CLType
    ccall( (:clang_getEnumDeclIntegerType, libclang), CXType, (CXCursor,), cu)
end
function getEnumConstantDeclValue(cu)
    ccall( (:clang_getEnumConstantDeclValue, libclang), Int64, (CXCursor,), cu)
end
function getEnumConstantDeclUnsignedValue(cu)
    ccall( (:clang_getEnumConstantDeclUnsignedValue, libclang), UInt64, (CXCursor,), cu)
end
function Cursor_getNumArguments(cu)
    ccall( (:clang_Cursor_getNumArguments, libclang), Int32, (CXCursor,), cu)
end
function Cursor_getArgument(cu,arg::Int32)::CLCursor
    ccall( (:clang_Cursor_getArgument, libclang), CXCursor, (CXCursor,Int32), cu, arg)
end
function equalTypes(t1,t2)
    # TODO doc here is weird
    #   "Returns non-zero if the CXTypes represent the same type and zero otherwise."
    r = ccall( (:clang_equalTypes, libclang), UInt32, (CXType,CXType), t1, t2)
    return r != 0
end
function getCanonicalType(t)::CLType
    ccall( (:clang_getCanonicalType, libclang), CXType, (CXType,), t)
end
function isConstQualifiedType(t)
    ccall( (:clang_isConstQualifiedType, libclang), UInt32, (CXType,), t) == 1
end
function isVolatileQualifiedType(t)
    ccall( (:clang_isVolatileQualifiedType, libclang), UInt32, (CXType,), t) == 1
end
function isRestrictQualifiedType(t)
    ccall( (:clang_isRestrictQualifiedType, libclang), UInt32, (CXType,), t)
end
function getPointeeType(t)::CLType
    ccall( (:clang_getPointeeType, libclang), CXType, (CXType,), t)
end
function getTypeDeclaration(t)::CLCursor
    ccall( (:clang_getTypeDeclaration, libclang), CXCursor, (CXType,), t)
end
function getDeclObjCTypeEncoding(cu)::String
    ccall( (:clang_getDeclObjCTypeEncoding, libclang), CXString, (CXCursor,), cu)
end
function getTypeKindSpelling(tk)
    s = ccall( (:clang_getTypeKindSpelling, libclang), CXString, (CXTypeKind,), tk)
end
function getFunctionTypeCallingConv(t::CXType)::CLCallingConv
    ccall( (:clang_getFunctionTypeCallingConv, libclang), Int32, (CXType,), t)
end
function getResultType(t)::CLType
    ccall( (:clang_getResultType, libclang), CXType, (CXType,), t)
end
function getNumArgTypes(t)
    ccall( (:clang_getNumArgTypes, libclang), Int32, (CXType,), t)
end
function getArgType(t, arg::UInt32)::CLType
    ccall( (:clang_getArgType, libclang), CXType, (CXType,UInt32), t, arg)
end
function isFunctionTypeVariadic(t)
    ccall( (:clang_isFunctionTypeVariadic, libclang), UInt32, (CXType,), t) == 1
end
function getCursorResultType(cu)::CLType
    ccall( (:clang_getCursorResultType, libclang), CXType, (CXCursor,), cu)
end
function isPODType(t)
    ccall( (:clang_isPODType, libclang), UInt32, (CXType,), t) == 1
end
function getElementType(t)::CLType
    ccall( (:clang_getElementType, libclang), CXType, (CXType,), t)
end
function getNumElements(t)
    ccall( (:clang_getNumElements, libclang), Int64, (CXType,), t)
end
function getArrayElementType(t)::CLType
    ccall( (:clang_getArrayElementType, libclang), CXType, (CXType,), t)
end
function getArraySize(t)
    ccall( (:clang_getArraySize, libclang), Int64, (CXType,), t)
end
function isVirtualBase(cu)
    ccall( (:clang_isVirtualBase, libclang), UInt32, (CXCursor,), cu) == 1
end
function getCXXAccessSpecifier(cu)::CLCXXAccessSpecifier
    ccall( (:clang_getCXXAccessSpecifier, libclang), Int64, (CXCursor,), cu)
end
function getNumOverloadedDecls(cu)
    ccall( (:clang_getNumOverloadedDecls, libclang), UInt32, (CXCursor,), cu)
end
function getOverloadedDecl(cu,arg::UInt32)::CLCursor
    ccall( (:clang_getOverloadedDecl, libclang), CXCursor, (CXCursor,UInt32), cu, arg)
end
function getCursorUSR(cu)::String
    c = ccall( (:clang_getCursorUSR, libclang), CXString, (CXCursor,), cu)
end
function getCursorSpelling(cu)::String
    c = ccall( (:clang_getCursorSpelling, libclang), CXString, (CXCursor,), cu)
end
function getCursorDisplayName(cu)::String
    c = ccall( (:clang_getCursorDisplayName, libclang), CXString, (CXCursor,), cu)
end
function getCursorReferenced(cu)::CLCursor
    ccall( (:clang_getCursorReferenced, libclang), CXCursor, (CXCursor,), cu)
end
function getCursorDefinition(cu)::CLCursor
    ccall( (:clang_getCursorDefinition, libclang), CXCursor, (CXCursor,), cu)
end
function isCursorDefinition(cu)
    ccall( (:clang_isCursorDefinition, libclang), UInt32, (CXCursor,), cu) == 1
end
function getCanonicalCursor(cu)::CLCursor
    ccall( (:clang_getCanonicalCursor, libclang), CXCursor, (CXCursor,), cu)
end
function CXXMethod_isStatic(cu)
    ccall( (:clang_CXXMethod_isStatic, libclang), UInt32, (CXCursor,), cu) == 1
end
function CXXMethod_isVirtual(cu)
    ccall( (:clang_CXXMethod_isVirtual, libclang), UInt32, (CXCursor,), cu) == 1
end
function CXXMethod_isProtected(cu)
    ccall( (:clang_CXXMethod_isProtected, libclang), UInt32, (CXCursor,), cu) == 1
end
function getTemplateCursorKind(cu)::CLCursorKind
    ccall( (:clang_getTemplateCursorKind, libclang), Int32, (CXCursor,), cu)
end
function getSpecializedCursorTemplate(cu)::CLCursor
    ccall( (:clang_getSpecializedCursorTemplate, libclang), CXCursor, (CXCursor,), cu)
end
function getTokenKind(cu)
    ccall( (:clang_getTokenKind, libclang), Int32, (CXCursor,), cu)
end
function getTokenSpelling(tu,token)::String
    ccall( (:clang_getTokenSpelling, libclang), CXString, (CXTranslationUnit,CXToken), tu, token)
end
function getCursorKindSpelling(ck)::String
    ccall( (:clang_getCursorKindSpelling, libclang), CXString, (CXCursorKind,), ck)
end
function getClangVersion()
    c = ccall( (:clang_getClangVersion, libclang), CXString, ())
end
function getCursorExtent(cu)
    ccall( (:clang_getCursorExtent, libclang), CXSourceRange, (CXCursor,), cu)
end
function disposeTokens(tu::CXTranslationUnit, tokens::Ptr{CXToken}, numtokens::Cuint)
    ccall( (:clang_disposeTokens, libclang), Nothing, (Ptr{Cvoid}, Ptr{CXToken}, Cuint),
           tu, tokens, numtokens)
end
function tokenize(tu::CXTranslationUnit, sr::CXSourceRange)
    token_p = Ref{Ptr{CXToken}}()
    numtokens = Ref{Cuint}()

    ccall(  (:clang_tokenize, libclang), Cuint,
            (CXTranslationUnit, CXSourceRange, Ref{Ptr{CXToken}}, Ref{Cuint}),
            tu, sr, token_p, numtokens)
    return TokenList(token_p[], convert(Cuint, numtokens[]), tu)
end
function Cursor_getTranslationUnit(cu)
    ccall( (:clang_Cursor_getTranslationUnit, libclang), CXTranslationUnit, (CXCursor,), cu)
end

ext: jl


text:
@testset "test_discrete" begin

C_111 = ss([-5], [2], [3], [0])
C_212 = ss([-5 -3; 2 -9], [1; 2], [1 0; 0 1], [0; 0])
C_221 = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0], [0 0])
C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], [1 0; 0 1], [1 0; 0 1])

@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :zoh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test c2d_x0map(ss(4*[1 0; 0 1]), 0.5, :foh) == (ss(4*[1 0; 0 1], 0.5), zeros(0, 2))
@test_tupleapprox(c2d_x0map(C_111, 0.01, :zoh),
    ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0])
@test_tupleapprox(c2d_x0map(C_111, 0.01, :foh),
    ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01),
    [1 -0.009835396005712075])
@test_tupleapprox(c2d_x0map(C_212, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009466805409666932; 0.019219966830212765], [1 0; 0 1], [0; 0], 0.01),
    [1 0 0; 0 1 0])
@test_tupleapprox(c2d_x0map(C_212, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.008957940478201584; 0.018468989584974498], [1 0; 0 1], [0.004820885889482196;
    0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675])
@test_tupleapprox(c2d_x0map(C_221, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0])
@test_tupleapprox(c2d_x0map(C_221, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5;
    0 1 -3.219047211142736e-5 -0.009706152723187249])
@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :zoh),
    ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819],
    [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465],
    [1 0; 0 1], [1 0; 0 1], 0.01), [1 0 0 0; 0 1 0 0])
@test_tupleapprox(c2d_x0map(C_222_d, 0.01, :foh),
    ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482],
    [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387],
    [1 0; 0 1], [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01),
    [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249])

# Test c2d for transfer functions
G = tf([1, 1], [1, 3, 1])
Gd = c2d(G, 0.2)
@test Gd ≈ tf([0, 0.165883310712090, -0.135903621603238], [1.0, -1.518831946985175, 0.548811636094027], 0.2) rtol=1e-14

# Issue #391
@test c2d(tf(C_111), 0.01, :zoh) ≈ tf(c2d(C_111, 0.01, :zoh))
@test c2d(tf(C_111), 0.01, :foh) ≈ tf(c2d(C_111, 0.01, :foh))

# Issue #471 discussion, test c2d for zero row C
C_210 = ss(C_212.A, C_212.B, zeros(0, 2), zeros(0, 1))
@test c2d(C_210, 0.01).A ≈ c2d(C_212, 0.01).A

# c2d on a zpk model should arguably return a zpk model
@test_broken typeof(c2d(zpk(G), 1)) <: TransferFunction{<:ControlSystems.SisoZpk}

# ERRORS
@test_throws MethodError c2d(ss([1], [2], [3], [4], 0.01), 0.01)   # Already discrete
@test_throws MethodError c2d(ss([1], [2], [3], [4], -1), 0.01)     # Already discrete

# d2c
@static if VERSION > v"1.4" # log(matrix) is buggy on previous versions, should be fixed in 1.4 and back-ported to 1.0.6
    @test d2c(c2d(C_111, 0.01)) ≈ C_111
    @test d2c(c2d(C_212, 0.01)) ≈ C_212
    @test d2c(c2d(C_221, 0.01)) ≈ C_221
    @test d2c(c2d(C_222_d, 0.01)) ≈ C_222_d
    @test d2c(Gd) ≈ G

    sys = ss([0 1; 0 0], [0;1], [1 0], 0)
    sysd = c2d(sys, 1)
    @test d2c(sysd) ≈ sys
end

# forward euler
@test c2d(C_111, 1, :fwdeuler).A == I + C_111.A
@test d2c(c2d(C_111, 0.01, :fwdeuler), :fwdeuler) ≈ C_111
@test d2c(c2d(C_212, 0.01, :fwdeuler), :fwdeuler) ≈ C_212
@test d2c(c2d(C_221, 0.01, :fwdeuler), :fwdeuler) ≈ C_221
@test d2c(c2d(C_222_d, 0.01, :fwdeuler), :fwdeuler) ≈ C_222_d
@test d2c(c2d(G, 0.01, :fwdeuler), :fwdeuler) ≈ G

Cd = c2d(C_111, 0.001, :fwdeuler)
t = 0:0.001:2
y,_ = step(C_111, t)
yd,_ = step(Cd, t)
@test norm(yd-y) / norm(y) < 1e-3

# dab
import DSP: conv
ζ = 0.2
ω = 1
B = [1]
A   = [1, 2ζ*ω, ω^2]
P  = tf(B,A)
# Control design
ζ0 = 0.7
ω0 = 2
Am = [1, 2ζ0*ω0, ω0^2]
Ao = conv(2Am, [1/2, 1]) # Observer polynomial, add extra pole due to the integrator
AR = [1,0] # Force the controller to contain an integrator ( 1/(s+0) )

B⁺  = [1] # The process numerator polynomial can be facored as B = B⁺B⁻ where B⁻ contains the zeros we do not want to cancel (non-minimum phase and poorly damped zeros)
B⁻  = [1]
Bm  = conv(B⁺, B⁻) # In this case, keep the entire numerator polynomial of the process

R,S,T = rstc(B⁺,B⁻,A,Bm,Am,Ao,AR) # Calculate the 2-DOF controller polynomials

Gcl = tf(conv(B,T),zpconv(A,R,B,S)) # Form the closed loop polynomial from reference to output, the closed-loop characteristic polynomial is AR + BS, the function zpconv takes care of the polynomial multiplication and makes sure the coefficient vectores are of equal length

@test ControlSystems.isstable(Gcl)

p = pole(Gcl)
# Test that all desired poles are in the closed-loop system
@test norm(minimum(abs.((pole(tf(Bm,Am)) .- sort(p, by=imag)')), dims=2)) < 1e-6
# Test that the observer poles are in the closed-loop system
@test norm(minimum(abs.((pole(tf(1,Ao)) .- sort(p, by=imag)')), dims=2)) < 1e-6

end


text:
#!/usr/bin/env julia

using KeyboardDriven
using Test

@time begin
@time @testset "Bytes to key" begin include("test_bytestokey.jl") end
end


text:
using FixedSizeArrays, GeometryTypes, Colors, GLAbstraction
import GLAbstraction: opengl_prefix
abstract TranspilerBackend

immutable FuncExpr
    name::Symbol
    slots
    args
    body::Vector
    returntype::DataType
    static_parameters::Tuple
end

immutable Light{T}
    position::Vec{3,T}
    ambient::Vec{3,T}
    diffuse::Vec{3,T}
    diffuse_power::T
    specular::Vec{3,T}
    specular_power::T
end

immutable Shading{T}
    ambient::Vec{3, T}
    specular::Vec{3, T}
    shininess::T
end
function get_lambda(pass, f, types)
    lambda = pass(f, types)
    if isa(lambda, Vector)
        if isempty(lambda)
            args = join(map(t->"::$t", types), ", ")
            error("$f($args) couldn't be found")
        end
        return first(lambda)
    else
        isa(lambda, LambdaInfo) && return lambda
        error("Not sure what's up with returntype of $pass. Returned: $lambda")
    end
end

function FuncExpr(f, types)
    lam_lowered = get_lambda(code_lowered, f, types)
    lam_typed = get_lambda(code_typed, f, types)
    slotnames = Base.lambdainfo_slotnames(lam_lowered)
    slottypes = lam_typed.slottypes
    nargs = length(types)
    slots = [(SlotNumber(i), (replace(name, "@", ""), slottypes[i])) for (i, name) in enumerate(slotnames)]
    ast = Base.uncompressed_ast(lam_lowered)
    name = Symbol(f)
    static_params = tuple(lam_typed.sparam_vals...)
    FuncExpr(name, slots, view(slots, 2:(nargs+1)), ast, lam_typed.rettype, static_params)
end

immutable Declaration
    dependencies::Vector{Declaration}
    source::String
end

type Backend{Identifier} <: TranspilerBackend
    #declarations::Dict{String, Declaration}
    current_slots::Vector{Dict{Slot, Tuple{String, DataType}}}
    type_map::Dict{DataType, String}
    function_map::Dict{Tuple{String, Vector{DataType}}, Tuple{String, DataType}}
    indentation_level::Int
end
type GTBlock
    entry
    block
    out
end
GTBlock(entry) = GTBlock(entry, [], Nullable{GotoNode}())
"""
Transpiled Expression
"""
immutable TExpr
    source::String
    typ::DataType
end
returntype(x::TExpr) = x.typ
source(x::TExpr) = x.source

typealias GL Backend{:GL}

function Backend(sym::Symbol)
    Backend{sym}(
        Dict{Slot, Tuple{String, DataType}}[],
        Dict{DataType, String}(),
        Dict{Tuple{String, Tuple}, Tuple{String, DataType}}(),
        1
    )
end

function insert_function!(backend::GL, jlfun, glsl_fun, types_ret::Pair)
    types, returntype = types_ret
    backend.function_map[(string(jlfun), types)] = (string(glsl_fun), returntype)
end

function insert_buildins!(backend::GL)
    # TODO do this in a cool way!
    # e.g define Module buildins with glsl inline and other tricks
    types = (Float32, Float64, Int32, UInt32)

    type_map = Dict(
        Float32 => "float",
        Float64 => "double",
        Int32   => "int",
        UInt32  => "uint"
    )
    for i=1:4
        for T in (Vec, Point)
            for ET in types
                VT = T{i, ET}
                if i==1
                    type_map[VT] = type_map[ET]
                else
                    type_map[VT] = string(opengl_prefix(ET), "vec", i)
                end
                for j=1:4
                    type_map[Mat{i,j,ET}] = string(
                        opengl_prefix(ET), "mat", i==j ? "$i" : "$(i)x$(j)"
                    )
                end
            end
        end
    end
    type_map[RGB{Float32}] = "vec3"
    type_map[RGBA{Float32}] = "vec4"
    merge!(backend.type_map, type_map)
    buildins = Dict(
        3 => Symbol[:round,:ceil,:trunc,:floor,:fma,:any,:all,:clamp],
        2 => Symbol[:floor,:max,:any,:ldexp,:ceil,:normalize,:all,:log,:mod,:cross,:trunc,:round,:min],
        1 => Symbol[:atanh,:tanh,:exp2,:floor,:max,:any,:sin,:tan,:modf,:step,:log2,:acosh,:ceil,:normalize,:length,:atan,:all,:log,:isnan,:trunc,:asinh,:transpose,:abs,:frexp,:sqrt,:round,:sinh,:sign,:acos,:asin,:cosh,:exp,:min,:isinf,:cos]
    )
    ops = (
        ".*" => "*",
        ".+" => "+",
        "./" => "/",
        "*" => "*",
        "+" => "+",
        "/" => "/",
        # "div_float" => "/",
        # "sub_float" => "-",
        # "add_float" => "+"
    )
    boolfuns = [:all, :any]
    # TODO actually figure out the signature of the GLSL build-in functions
    glsl_types = unique(collect(keys(backend.type_map)))
    for T in glsl_types
        if !isa(T, Texture) && !isa(T, Mat)
            for i=1:3
                for fun in buildins[i]
                    insert_function!(backend, string(fun), string(fun), [T for _=1:i] => T)
                end
            end
            if T <: FixedVector
                insert_function!(backend, "dot", "dot", [T, T] => eltype(T))
            end
            for (jl_op, glsl_op) in ops
                insert_function!(backend, jl_op, glsl_op, [T, T] => T)
            end
            for jl_op in ("*", "-", "+")
                insert_function!(backend, jl_op, jl_op, [T] => T)
            end
        end
    end
    # cases not yet covered generically
    insert_function!(backend, "Float32", "float", [Float32] => Float32)
    insert_function!(backend, "Float32", "float", [Int] => Float32)
    insert_function!(backend, "Float32", "float", [Int32] => Float32)
    insert_function!(backend, "Float32", "float", [UInt] => Float32)
    insert_function!(backend, "Float32", "float", [UInt32] => Float32)

    insert_function!(backend, "norm", "normalize", [Vec3f0] => Float32)
    insert_function!(backend, "^", "pow", [Float32, Float32] => Float32)
    insert_function!(backend, "^", "pow", [Int32, Int32] => Float32)

    insert_function!(backend, Symbol(Vec3f0), "vec3", [Float32, Float32, Float32] => Vec3f0)
    insert_function!(backend, Symbol(Vec3f0), "vec3", [Vec3f0] => Vec3f0)
    insert_function!(backend, "<=", "<=", [Float32, Float32] => Bool)
    insert_function!(backend, ".*", "*", [Vec3f0, Float32] => Vec3f0)
    insert_function!(backend, ".*", "*", [Float32, Vec3f0] => Vec3f0)
    insert_function!(backend, "*", "*", [Vec3f0, Float32] => Vec3f0)
    insert_function!(backend, "*", "*", [Float32, Vec3f0] => Vec3f0)
    insert_function!(GLBackend, "+", "+", [Vec3f0, Vec3f0, Vec3f0] => Vec3f0)

end

const GLBackend = Backend(:GL)
insert_buildins!(GLBackend)

isbox(x) = (x == GlobalRef(Base, :box))
extract_sym(x) = error("$x $(typeof(x)) cant be converted to Symbol")
extract_sym(x::Symbol) = x
extract_sym(x::QuoteNode) = x.value
get_types(x::AbstractVector) = map(slot_name_t -> slot_name_t[2][2], x)

"""
Type parameter are not supported, so we'll put them in the name
"""
function to_type_name(T)
    join((T.name.name, T.parameters...), "_")
end
function to_type_name(backend::GL, T::DataType)
    get!(backend.type_map, T) do
        declare_composite_type(backend, T)
    end
end

function to_type_assert(backend::GL, name::String, typ::String)
    "$typ $name"
end
function to_type_assert(backend::GL, name::String, typ::DataType)
    tn = to_type_name(backend, typ)
    to_type_assert(backend, name, tn)
end

function declare_slots(backend::GL, slots)
    join(map(slots) do slot_name_t
        slot, (name, t) = slot_name_t
        slot_declare = to_type_assert(backend, name, t)
        declare_block(backend, slot_declare)
    end, "")
end

function declare_block(backend::GL, source::String, level=backend.indentation_level)
    repeat("    ", level)*source*";\n"
end
function declare_composite_type(backend::GL, T::DataType)
    name = to_type_name(T)
    get!(backend.type_map, T) do
        fields = [to_type_assert(backend, string(n), fieldtype(T, n)) for n in fieldnames(T)]
        #TODO actualy insert this into declarations
        println("""
struct $(name){
$(join(fields, ";\n"));
};
        """)
        name
    end
end

function declare_function(backend::GL, func::FuncExpr)
    types = get_types(func.args)
    name, rettype = get!(backend.function_map, (string(func.name), types)) do
        arg_decl = map(arg->to_type_assert(backend, arg[2]...), func.args)
        rettype = to_type_name(backend, func.returntype)
        transpilat = transpile(backend, func.body)
        src = source(transpilat)
        arg_end = last(func.args.indexes[1])
        var_decl = declare_slots(backend::GL, func.slots[(arg_end+1):end])
        #@assert returntype(transpilat) == func.returntype
        name = func.name
        #TODO actualy insert this into declarations
        println("""
$rettype $name($(join(arg_decl, ", "))){
$var_decl
$src
}
""")
        string(name), func.returntype
    end
    name, rettype
end

function transpile(backend::GL, f::Symbol, types)
    func = getfield(current_module(), f)
    transpile(backend, func, types)
end
function transpile(backend::GL, f::GlobalRef, types)
    func = getfield(f.mod, f.name)
    transpile(backend, func, types)
end
function transpile(backend::GL, f::Union{Function, Type}, types::AbstractVector)
    transpile(backend, f, tuple(types...))
end

function remove_static_param!(expr, static_parameters)
    expr # leave any other untouched
end
function remove_static_param!(expr::Vector, static_parameters)
    map!(e->remove_static_param!(e, static_parameters), expr)
end
function remove_static_param!(expr::Expr, static_parameters)
    if expr.head == :static_parameter
        idx = expr.args[1]
        return static_parameters[idx]
    else
        remove_static_param!(expr.args, static_parameters)
        expr
    end
end
function transpile(backend::GL, f::Union{Function, Type}, types::Tuple)
    name, rettype = get!(backend.function_map, (string(Symbol(f)), [types...])) do
        func = FuncExpr(f, types)
        # convert slot types to glsl slot types
        map!(func.slots) do slot_name_t
            s, (n, t) = slot_name_t
            if applicable(GLAbstraction.gl_promote, t)
                t = GLAbstraction.gl_promote(t)
            end
            s, (n, t)
        end
        remove_static_param!(func.body, func.static_parameters)
        push!(backend.current_slots, Dict(func.slots))
        decl = declare_function(backend, func)
        pop!(backend.current_slots)
        decl
    end
    name, rettype
end

function transpile(backend::GL, block::Expr)
    if block.head == :block || block.head == :body
        transpile(backend, block.args)
    else
        # I knoow, Val no good here... But it makes the function overwritable and composable, which might be nice
        # in order to have unspecialized fall backs which can be specialized
        #println("will transpile: $(block.head) ", block)
        return transpile(backend, block, Val{block.head}())
    end
end

function is_if(block::GTBlock)
    (
        isa(block.entry, Expr) &&
        block.entry.head == :gotoifnot
    )
end
function is_elseif(block::GTBlock)
    (
        isa(block.entry, Tuple) &&
        isa(first(block.entry), LabelNode) &&
        isa(last(block.entry), Expr) &&
        last(block.entry).head == :gotoifnot
    )
end
function is_else(block::GTBlock)
    isa(block.entry, LabelNode) && isnull(block.out)
end

# (╯°□°）╯︵ ┻━┻
function deal_with_goto(backend, elem, state, list, block_source)
    labels = Set(Int[])
    blocks = GTBlock[]
    while true
        if isa(elem, Expr) && elem.head == :gotoifnot
            condition, label = elem.args
            push!(labels, label)
            push!(blocks, GTBlock(elem))
        elseif (isa(elem, LabelNode) && elem.label in labels)
            setdiff!(labels, elem.label)
            isempty(labels) && break # found last label
            elem0, state0 = elem, state
            elem, state = next(list, state)
            while isa(elem, LineNumberNode)
                elem, state = next(list, state)
            end
            if !(isa(elem, Expr) && elem.head == :gotoifnot)
                elem, state = elem0, state0 # reset peek
                push!(blocks, GTBlock(elem))
            else
                push!(labels, elem.args[2])
                push!(blocks, GTBlock((elem0, elem)))
            end
        elseif isa(elem, GotoNode)
            last(blocks).out = Nullable(elem) #close block
            push!(labels, elem.label)
        elseif isa(elem, Expr) && elem.head == :gotoifnot
            condition, label = elem.args
            push!(labels, label)
        else
            push!(last(blocks).block, elem) # insert statements of block
        end
        done(list, state) && break
        elem, state = next(list, state)
    end
    if done(list, state) && !isempty(labels)
        error("Oh my... Couldn't find all goto labels")
    end
    indents = backend.indentation_level
    indspace = repeat("    ", indents)
    backend.indentation_level += 1
    for block in blocks
        bodysrc = source(transpile(backend, block.block))
        if is_if(block)
            cond = block.entry.args[1]
            condsrc = source(transpile(backend, cond))
            block_source *= """
$(indspace)if($condsrc){
$bodysrc
$(indspace)}
"""
        elseif is_elseif(block)
            cond = block.entry[2].args[1]
            condsrc = source(transpile(backend, cond))
            block_source *= """
$(indspace)else if($condsrc){
$bodysrc
$(indspace)}
"""
        elseif is_else(block)
            block_source *= """
$(indspace)else{
$bodysrc
$(indspace)}
"""
        else
            error("Block $block not supported yet")
        end
    end
    backend.indentation_level = indents
    elem, state = next(list, state) # move one after last labelnode
    elem, state, block_source
end

function transpile(backend::GL, list::AbstractVector)
    last_type = Any
    block_source = ""
    state = start(list)
    while !done(list, state)
        elem, state = next(list, state)
        if isa(elem, Expr) && (elem.head == :gotoifnot || isa(elem, GotoNode))
            elem, state, block_source = deal_with_goto(backend, elem, state, list, block_source)
        end
        line = transpile(backend, elem)
        line == nothing && continue
        block_source *= declare_block(backend, source(line))
        last_type = returntype(line)
    end
    return TExpr(block_source, last_type)
end

transpile(backend::GL, expr::Void) = nothing

function transpile(backend::GL, expr::LineNumberNode)
    TExpr("//julia source line: $(expr.line)", Void)
end
"""
concrete values
"""
function transpile(backend::GL, value)
    cnv_value = gl_convert(value)
    TExpr(string(cnv_value), typeof(cnv_value))
end

function transpile(backend::GL, expr::Slot)
    src_t = last(backend.current_slots)[expr]
    TExpr(src_t...)
end
function transpile(backend::GL, expr::Expr, ::Val{:(=)})
    lh, rh = map(x->transpile(backend, x), expr.args)
    TExpr("$(source(lh)) = $(source(rh))", Void)
end

function transpile(backend::GL, expr::Expr, ::Val{:call})
    fun = shift!(expr.args)

    if isbox(fun) # ignore boxes, LOL \( ﾟ◡ﾟ)/
        @assert length(expr.args) == 2
        ret_type, inner_expr = expr.args
        t = transpile(backend, inner_expr)
        @assert returntype(t) == ret_type
        return t
    end
    if fun == GlobalRef(Core, :getfield)
        @assert length(expr.args) == 2
        expr, field = expr.args
        expr_trp = transpile(backend, expr)
        field_s = extract_sym(field)
        typ = returntype(expr_trp)
        return TExpr("$(source(expr_trp)).$(field_s)", fieldtype(typ, field_s))
    end
    args = map(x->transpile(backend, x), expr.args)

    fun_name, return_type = transpile(backend, fun, returntype.(args))
    transpilat = source.(args)
    call_expr = if Base.isoperator(Symbol(fun_name)) # most Julia operators should also be operators in GLSL
        if length(transpilat) == 1
            "$(fun_name)$(transpilat[1])"
        else
            foldl(transpilat) do v0, v1
                "$v0 $fun_name $v1"
            end
        end
    else
        "$fun_name($(join(transpilat, ", ")))"
    end
    TExpr(call_expr, return_type)
end

function transpile(backend::GL, expr::Expr, ::Val{:return})
    src = transpile(backend, expr.args[1])
    TExpr("return $(source(src))", returntype(src))
end

# Test

function blinnphong4{NV, T}(light, L::Vec{NV, T}, N, V, color, shading)
    diff_coeff = max(dot(L,N), T(0.0))
    # specular coefficient
    H = normalize(L+V)

    spec_coeff = max(dot(H,N), 0.0)^(shading.shininess)
    if diff_coeff <= 0.0
        spec_coeff = 0.0
    elseif diff_coeff <= 0.2
        # some nonesense to test elseif
        spec_coeff *= 2.0;
        spec_coeff += 1.0;
    else
        spec_coeff = spec_coeff;
        return L
    end
    # final lighting model
    return Vec3f0(
        light.ambient .* shading.ambient +
        light.diffuse .* light.diffuse_power .* color * diff_coeff +
        light.specular .* light.specular_power .* shading.specular * spec_coeff
    )
end

l = Light(Vec3f0(0), Vec3f0(0), Vec3f0(0), 1f0, Vec3f0(0), 1f0)
s = Shading(Vec3f0(0), Vec3f0(0), 1f0)
blinnphong3(l, Vec3f0(0), Vec3f0(0), Vec3f0(0), Vec3f0(0), s)
fe = FuncExpr(blinnphong3, (Light{Float32},Vec3f0, Vec3f0, Vec3f0, Vec3f0,Shading{Float32}))
fe.slots
transpile(GLBackend, blinnphong4, (Light{Float32},Vec3f0, Vec3f0, Vec3f0, Vec3f0,Shading{Float32}))
#delete!(GLBackend.function_map, ("blinnphong3", [Light{Float32},Vec3f0, Vec3f0, Vec3f0, Vec3f0,Shading{Float32}]))


text:
"""
StanDump --- a package for writing data in the CmdStan dump data format.

The single exported function is [`stan_dump`](@ref).
"""
module StanDump

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES

export stan_dump

####
#### Representation of formatting
####

"""
Wrapper for an IO stream for writing data for use by Stan. See constructor for documentation
of slots.

Part of the API, but not exported.
"""
struct StanDumpIO{T <: IO}
    io::T
    def_arrow::Bool
    def_newline::Bool
    compact::Bool
end

"""
    StanDumpIO(io; def_arrow = false, def_newline = false, compact = false)

Wrap an IO stream `io` for writing data to be read by Stan.

# Arguments

- `def_arrow::Bool`: when `true` use `<-`, otherwise `=` for variable definitions.

- `def_newline::Bool`: when `true`, each `=` or `<-` is followed by a newline.

- `compact::Bool`: when `true`, drop whitespace when possible.
"""
function StanDumpIO(io; def_arrow = true, def_newline = false, compact = false)
    StanDumpIO(io, def_arrow, def_newline, compact)
end

####
#### Internals
####

"""
    dump(sd, xs...)

Write arguments `xs...` as data for Stan into `sd`. For internal use.

NOTE: Define methods only for valid Stan objects, using `_dump` for everything else.
"""
dump(sd::StanDumpIO, x) = throw(ArgumentError("Can't represent $(x) as data for Stan."))

dump(sd::StanDumpIO, x::Float64) = print(sd.io, x)

dump(sd::StanDumpIO, x::Real) = dump(sd, Float64(x))

function dump(sd::StanDumpIO, x::Integer)
    if typemin(Int32) ≤ x ≤ typemax(Int32)
        print(sd.io, x)
    elseif typemin(Int64) ≤ x ≤ typemax(Int64)
        print(sd.io, x, "L")
    else
        throw(ArgumentError("Integer too large to represent in Stan."))
    end
end

"""
$(SIGNATURES)

Write arguments arguments as data for Stan into `sd`, passing through strings and
characters, and allowing other special objects which are not valid data.

For internal implementation. also defined for objects which are not valid in Stan.
"""
_dump(sd::StanDumpIO, xs...) = for x in xs _dump(sd, x) end

_dump(sd::StanDumpIO, x) = dump(sd, x)

_dump(sd::StanDumpIO, x::Union{Char,String}) = print(sd.io, x)

"Write a space unless output is requested to be compact."
struct Space end

const SPACE = Space()

_dump(sd::StanDumpIO, ::Space) = sd.compact || print(sd.io, " ")

"""
$(SIGNATURES)

Test if the argument is valid as a Stan variable name.

NOTE: only basic checks, does not test conflicts with reserved names.
"""
function is_valid_varname(name::String)
    isvalid(c) = isascii(c) && (isdigit(c) || isletter(c) || c == '_')
    all(isvalid, name) && isletter(name[1]) && !endswith(name, "__")
end

function _dump(sd::StanDumpIO, x::Symbol)
    v = string(x)
    @argcheck is_valid_stan_varname(v) "Invalid variable name $(v)."
    print(sd.io, v)
end

function _dump(sd::StanDumpIO, x::Pair)
    name, value = x
    @argcheck name isa Union{AbstractString, Symbol} "Use symbols or strings for variable names."
    varname = string(name)
    @argcheck is_valid_varname(varname)
    _dump(sd, varname, SPACE, sd.def_arrow ? "<-" : '=',
          sd.def_newline ? "\n" : SPACE, value, "\n")
end

"""
$(SIGNATURES)

Dump elements of a vector (or iterable, considered as a vector).

When the element type is `<:Integer`, dump as `Int`s, otherwise as `Float64`; all values are
converted for consistency.
"""
function _dump_vector(sd::StanDumpIO, x)
    S = eltype(x) <: Integer ? Int : Float64
    if isempty(x)
        _dump(sd, S ≡ Int ? "integer" : "double", "(0)")
    else
        _dump(sd, "c(")
        for (i, x) in enumerate(x)
            i > 1 && _dump(sd, ",", SPACE)
            dump(sd, S(x))      # convert value for consistency
        end
        _dump(sd, ")")
    end
end

dump(sd::StanDumpIO, x::AbstractVector) = _dump_vector(sd, x)

function dump(sd::StanDumpIO, r::UnitRange{<: Integer})
    isempty(r) ? _dump_vector(sd, r) : _dump(sd, minimum(r), ":", maximum(r))
end

function dump(sd::StanDumpIO, A::AbstractArray)
    _dump(sd, "structure(", view(A, :), ",", SPACE, ".Dim", SPACE, "=",
          SPACE, collect(size(A)), ")")
end

####
#### Interface
####

"""
    stan_dump(filename, data; force = false, kwargs...)
    stan_dump(io, data; kwargs...)
    stan_dump(StanDump.StanDumpIO(io; kwargs...), data)

Write `data`, which is a value that supports `pairs` (eg a `NamedTuple` or a `Dict`) to
`filename` or `io`.

Using a `filename`, it will not be overwritten unless `force = true` is specified.

Keyword arguments are passed to `StanDumpIO` to govern formatting (most users should not
care about this, except for debugging purposes).

# Example

```jldoctest
julia> stan_dump(stdout, (N = 1, a = 1:5, b = ones(2, 2)))
N <- 1
a <- 1:5
b <- structure(c(1.0, 1.0, 1.0, 1.0), .Dim = c(2, 2))
```
"""
stan_dump(sd::StanDumpIO, data) = foreach(p -> _dump(sd, p), pairs(data))

stan_dump(io::IO, data; kwargs...) = stan_dump(StanDumpIO(io; kwargs...), data)

function stan_dump(filename::AbstractString, data; force::Bool = false, kwargs...)
    @argcheck force || !ispath(filename) "$(filename) already exists"
    open(io -> stan_dump(StanDumpIO(io; kwargs...), data), filename, "w")
end

end # module


text:
"""
Bridges `CP.ValuePrecedence` to reification.
"""
struct ValuePrecedence2ReificationBridge{T} <: MOIBC.AbstractBridge
    vars_reif_value::Vector{MOI.VariableIndex}
    vars_reif_precv::Vector{MOI.VariableIndex}
    vars_reif_value_bin::Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    vars_reif_precv_bin::Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_reif_value::Vector{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}} # Compare values from 2 to end.
    cons_reif_precv::Vector{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}} # Compare values from 1 to end-1.
    cons_implication::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{ValuePrecedence2ReificationBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.ValuePrecedence{T},
) where {T}
    return MOIBC.bridge_constraint(
        ValuePrecedence2ReificationBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{ValuePrecedence2ReificationBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.ValuePrecedence{T},
) where {T}
    f_scalars = MOIU.scalarize(f)

    vars_reif_value, vars_reif_value_bin = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:(s.dimension - 1)])
    vars_reif_precv, vars_reif_precv_bin = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:(s.dimension - 1)])

    cons_reif_value = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}[
        MOI.add_constraint(
            model,
            MOIU.vectorize(
                MOI.ScalarAffineFunction{T}[
                    one(T) * vars_reif_value[i - 1],
                    f_scalars[i] - T(s.value),
                ]
            ),
            CP.Reification(MOI.EqualTo(zero(T)))
        )
        for i in 2:s.dimension
    ]

    cons_reif_precv = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}}[
        MOI.add_constraint(
            model,
            MOIU.vectorize(
                MOI.ScalarAffineFunction{T}[
                    one(T) * vars_reif_precv[i],
                    f_scalars[i] - T(s.before)
                ]
            ),
            CP.Reification(MOI.EqualTo(zero(T)))
        )
        for i in 1:(s.dimension - 1)
    ]

    # If there is one `value`, then there must be at least one `before` before.
    cons_implication = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}[
        MOI.add_constraint(
            model,
            vars_reif_value[i - 1] - sum(one(T) .* vars_reif_precv[1:(i - 1)]),
            MOI.LessThan(zero(T))
        )
        for i in 2:s.dimension
    ]

    return ValuePrecedence2ReificationBridge(vars_reif_value, vars_reif_precv, vars_reif_value_bin, vars_reif_precv_bin, cons_reif_value, cons_reif_precv, cons_implication)
end

function MOI.supports_constraint(
    ::Type{ValuePrecedence2ReificationBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.ValuePrecedence{T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{ValuePrecedence2ReificationBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{ValuePrecedence2ReificationBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
    ]
end

function MOI.get(b::ValuePrecedence2ReificationBridge, ::MOI.NumberOfVariables)
    return length(b.vars_reif_precv) + length(b.vars_reif_value)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_reif_precv_bin) + length(b.vars_reif_value_bin)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_reif_value) + length(b.cons_reif_precv)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T}
    },
) where {T}
    return length(b.cons_implication)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return vcat(b.vars_reif_precv, b.vars_reif_value)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return vcat(b.vars_reif_precv_bin, b.vars_reif_value_bin)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Reification{MOI.EqualTo{T}},
    },
) where {T}
    return vcat(b.cons_reif_value, b.cons_reif_precv)
end

function MOI.get(
    b::ValuePrecedence2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T}
    },
) where {T}
    return copy(b.cons_implication)
end


text:
using Documenter, CircGeometry

push!(LOAD_PATH,"../src/")

makedocs(;
    modules=[CircGeometry],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/pseastham/CircGeometry.jl/blob/{commit}{path}#L{line}",
    sitename="CircGeometry.jl",
    authors="Patrick Eastham",
    assets=String[],
)

deploydocs(;
    repo="github.com/pseastham/CircGeometry.jl",
)

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
using PyPlot, Images
using PyCall
@pyimport numpy

function makeDICimage(img_name="../images/mandrill.png";sigma = 0.1)

  raw_img = PyPlot.imread(img_name);
  img = convert(Image{Images.Gray},raw_img[:,:,1:3]);
  M = convert(Array,img);
  M = float(M);

  gX, gY = numpy.gradient(M);
  gY = gY + sigma*rand(size(gY));

  maxi = maximum(gY)
  mini = minimum(gY)
  clamped = 1 +  (-gY' + maxi) * 1 / ( mini-maxi);
  out_name = "../DICimages/DIC" * basename(img_name)
  imwrite(grayim(clamped),out_name)

  @printf("File outputted to ../DICimages/")
  return

end

max_stars_repo_name: bcollyer/DICkit

max_issues_repo_name: bcollyer/DICkit

max_forks_repo_name: bcollyer/DICkit


text:
"""
_column_to_pointer{T}(p::Pointer)

construct JSONPointer.Pointer with specified type
"""
function _column_to_pointer(token_string::AbstractString)::Pointer
    if !startswith(token_string, JSONPointer.TOKEN_PREFIX)
        token_string = "/" * token_string
    end
    if endswith(token_string, "}")
        x = split(token_string, "{")
        p = Pointer(x[1])
        T = jsontype_to_juliatype(x[2][1:end-1])

        return Pointer{Array{T, 1}}(p.tokens)
    else
        return Pointer(token_string)
    end
end
_column_to_pointer(token) = _column_to_pointer(string(token))

function jsontype_to_juliatype(t)
    if t == "string"
        return String
    elseif t == "number"
        return Float64
    # JSON does not have distinct types for integers and floating-point values
    # but Excel does, and distinguishing integer is useful for many things.
    elseif t == "integer"
        return Int
    elseif t == "object"
        return OrderedCollections.OrderedDict{String,Any}
    elseif t == "array"
        return Vector{Any}
    elseif t == "boolean"
        return Bool
    elseif t == "null"
        return Missing
    else
        error(
            "You specified a type that JSON doesn't recognize! Instead of " *
            "`::$t`, you must use one of `::string`, `::number`, " *
            "`::object`, `::array`, `::boolean`, or `::null`."
        )
    end
end

function pointer_to_colname(p::Pointer{T})::String where T
    col = "/" * join(p.tokens, "/")
    t = juliatype_to_jsontype(T)
    if t == "array"
        col *= "::$t"
        t2 = juliatype_to_jsontype(eltype(T))
        if !isempty(t2)
            col *= "{$t2}"
        end
    elseif !isempty(t)
        col *= "::$t"
    end
    return col
end

function juliatype_to_jsontype(T)
    if T <: OrderedDict
        t = "object"
    elseif T <: Array
        t = "array"
    elseif T == String
        t = "string"
    elseif T == Float64
        t = "number"
    elseif T == Int
        t = "integer"
    elseif T == Bool
        t = "boolean"
    elseif T == Missing
        t = "null"
    elseif T == Nothing
        t = "null"
    elseif T == Any
        t = ""
    else
        @warn("cannot find jsontype from $T, returning empty string")
        t = ""
    end
    return t
end


text:
"""
	kml2gmt(cmd0::String="", arg1=[], kwargs...)

kml2gmt - Extract GMT table data from Google Earth KML files

Full option list at [`kml2gmt`](http://gmt.soest.hawaii.edu/doc/latest/kml2gmt.html)

Parameters
----------

- **F** : **select** : -- Str --        Flags = s|l|p

    Specify a particular feature type to output. Choose from points (s), line (l), or polygon (p).
    By default we output all geometries.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/kml2gmt.html#f)
- **Z** : **altitudes** : -- Bool or [] --

    Output the altitude coordinates as GMT z coordinates [Default will output just longitude and latitude].
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/kml2gmt.html#z)
- $(GMT.opt_V)
- $(GMT.opt_bo)
- $(GMT.opt_do)
- $(GMT.opt_swap_xy)
"""
function kml2gmt(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("kml2gmt", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_bo(cmd, d)
	cmd, = parse_do(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'F', d, [:F :select])
	cmd = add_opt(cmd, 'Z', d, [:Z :altitudes])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "kml2gmt", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
kml2gmt(arg1=[], cmd0::String=""; kw...) = kml2gmt(cmd0, arg1; kw...)


text:
using TerminalCalendar
using Test
using Dates

@testset "TerminalCalendar.jl" begin
    feb2021 = Date(2021, 2)

    feb2021_auto = TerminalCalendar.calendarize(feb2021)
    # feb2021_manual = Matrix{Union{Missing, Date}}()
    feb_range = Date(2021, 2, 1):Day(1):Date(2021, 2, 28)
    feb_matrix = permutedims(reshape(feb_range, 7, :))

    @test feb2021_auto == feb_matrix
end


text:
using Test
import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 20.0)
        MOIU.mock_optimize!(mock, [4, 5, 1])
    end,
)
MOIT.int1test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 2]),
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 2]),
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
        mock,
        [0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 3.0, 12.0],
    ),
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
        mock,
        [0.0, 0.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0, 6.0, 24.0],
    ),
)
MOIT.int2test(mock, config)
# FIXME [1, 0...] is not the correct optimal solution but it passes the test
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0; zeros(10)]),
)
MOIT.int3test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 0, 0, 1, 1]),
)
MOIT.knapsacktest(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) ->
        MOIU.mock_optimize!(mock, [1.25, 8.75, 0.0, 1.0]),
)
MOIT.indicator1_test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) ->
        MOIU.mock_optimize!(mock, [2.0, 8.0, 1.0, 0.0]),
)
MOIT.indicator2_test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) ->
        MOIU.mock_optimize!(mock, [1.25, 8.75, 1.0, 1.0]),
)
MOIT.indicator3_test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) ->
        MOIU.mock_optimize!(mock, [1.25, 8.75, 0.0, 1.0]),
)
MOIT.indicator4_test(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 0.0)
        MOIU.mock_optimize!(mock, [0.0, 0.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 2.0)
        MOIU.mock_optimize!(mock, [2.0, 1.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 2.0)
        MOIU.mock_optimize!(mock, [2.0, 2.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 2.5)
        MOIU.mock_optimize!(mock, [2.5, 2.5])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 3.0)
        MOIU.mock_optimize!(mock, [3.0, 3.0])
    end,
    (mock::MOIU.MockOptimizer) ->
        MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
)
MOIT.semiconttest(mock, config)
MOIU.set_mock_optimize!(
    mock,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 0.0)
        MOIU.mock_optimize!(mock, [0.0, 0.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 2.0)
        MOIU.mock_optimize!(mock, [2.0, 1.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 2.0)
        MOIU.mock_optimize!(mock, [2.0, 2.0])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 3.0)
        MOIU.mock_optimize!(mock, [3.0, 2.5])
    end,
    (mock::MOIU.MockOptimizer) -> begin
        MOI.set(mock, MOI.ObjectiveBound(), 3.0)
        MOIU.mock_optimize!(mock, [3.0, 3.0])
    end,
    (mock::MOIU.MockOptimizer) ->
        MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
)
MOIT.semiinttest(mock, config)


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
const AFFEQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}

mutable struct Optimizer <: MOI.AbstractOptimizer
    objconstant::Cdouble
    objsign::Int
    blockdims::Vector{CSDP_INT}
    varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
    num_entries::Dict{Tuple{Int, Int}, Int}
    b::Vector{Cdouble}
    C::blockmatrix
    problem::Union{Nothing, LoadingProblem}
    X::blockmatrix
    y::Union{Nothing, Vector{Cdouble}}
    Z::blockmatrix
    status::CSDP_INT
    pobj::Cdouble
    dobj::Cdouble
    solve_time::Float64
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
        optimizer = new(
            zero(Cdouble), 1, CSDP_INT[], Tuple{Int, Int, Int}[],
            Dict{Tuple{Int, Int}, Int}(), Cdouble[],
            blockmatrix(), nothing, blockmatrix(), nothing, blockmatrix(),
            -1, NaN, NaN, NaN, false, Dict{Symbol, Any}())
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawParameter(String(key)), value)
        end
        # May need to call `free_loaded_prob` and `free_loading_prob`.
        finalizer(MOI.empty!, optimizer)
        return optimizer
    end
end

varmap(optimizer::Optimizer, vi::MOI.VariableIndex) = optimizer.varmap[vi.value]

function MOI.supports(optimizer::Optimizer, param::MOI.RawParameter)
    return Symbol(param.name) in ALLOWED_OPTIONS
end
function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
    if !(param.name isa String)
        Base.depwarn(
            "passing `$(param.name)` to `MOI.RawParameter` as type " *
            "`$(typeof(param.name))` is deprecated. Use a string instead.",
            Symbol("MOI.set")
        )
    end
    if !MOI.supports(optimizer, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    optimizer.options[Symbol(param.name)] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
    # TODO: This gives a poor error message if the name of the parameter is invalid.
    if !(param.name isa String)
        Base.depwarn(
            "passing `$(param.name)` to `MOI.RawParameter` as type " *
            "`$(typeof(param.name))` is deprecated. Use a string instead.",
            Symbol("MOI.set")
        )
    end
    return optimizer.options[Symbol(param.name)]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.get(::Optimizer, ::MOI.SolverName) = "CSDP"
# See table "Return codes for easy_sdp() and CSDP" in `doc/csdpuser.pdf`.
const RAW_STATUS = [
    "Problem solved to optimality.",
    "Problem is primal infeasible.",
    "Problem is dual infeasible.",
    "Problem solved to near optimality.",
    "Maximum iterations reached.",
    "Stuck at edge of primal feasibility.",
    "Stuck at edge of dual feasibility.",
    "Lack of progress.",
    "X, Z, or O is singular.",
    "NaN or Inf values encountered.",
    "Program stopped by signal (SIXCPU, SIGTERM, etc.)"]

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return RAW_STATUS[optimizer.status + 1]
end
function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
    return optimizer.solve_time
end

function MOI.is_empty(optimizer::Optimizer)
    return iszero(optimizer.objconstant) &&
        isone(optimizer.objsign) &&
        isempty(optimizer.blockdims) &&
        isempty(optimizer.varmap) &&
        isempty(optimizer.num_entries) &&
        isempty(optimizer.b) &&
        iszero(optimizer.C.nblocks) &&
        optimizer.C.blocks == C_NULL &&
        optimizer.problem === nothing
end

function MOI.empty!(optimizer::Optimizer)
    optimizer.objconstant = zero(Cdouble)
    optimizer.objsign = 1
    empty!(optimizer.blockdims)
    empty!(optimizer.varmap)
    empty!(optimizer.num_entries)
    empty!(optimizer.b)
    if optimizer.problem !== nothing
        if optimizer.y !== nothing
            free_loaded_prob(optimizer.problem, optimizer.X, optimizer.y, optimizer.Z)
        end
        free_loading_prob(optimizer.problem)
    end
    optimizer.problem = nothing
    optimizer.C.nblocks = 0
    optimizer.C.blocks = C_NULL
    optimizer.X.nblocks = 0
    optimizer.X.blocks = C_NULL
    optimizer.y = nothing
    optimizer.Z.nblocks = 0
    optimizer.Z.blocks = C_NULL
    optimizer.status = -1
    optimizer.pobj = 0.0
    optimizer.dobj = 0.0
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Cdouble}}})
    return true
end

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{MOI.Reals})
    return false
end
const SupportedSets = Union{MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables},
    ::Type{<:SupportedSets})
    return true
end
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Cdouble}},
    ::Type{MOI.EqualTo{Cdouble}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end
MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    # To be sure that it is done before load(optimizer, ::ObjectiveFunction, ...), we do it in allocate
    optimizer.objsign = sense == MOI.MIN_SENSE ? -1 : 1
end
function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction, ::MOI.ScalarAffineFunction) end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
# Loads objective coefficient α * vi
function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    # in SDP format, it is max and in MPB Conic format it is min
    coef = optimizer.objsign * α
    if i != j
        coef /= 2
    end
    addentry(optimizer.problem, 0, blk, i, j, coef, true)
end
function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    obj = MOIU.canonical(f)
    optimizer.objconstant = f.constant
    for t in obj.terms
        if !iszero(t.coefficient)
            load_objective_term!(optimizer, t.coefficient, t.variable_index)
        end
    end
end

function new_block(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.blockdims, -MOI.dimension(set))
    blk = length(optimizer.blockdims)
    for i in 1:MOI.dimension(set)
        push!(optimizer.varmap, (blk, i, i))
    end
end

function new_block(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.blockdims, set.side_dimension)
    blk = length(optimizer.blockdims)
    for i in 1:set.side_dimension
        for j in 1:i
            push!(optimizer.varmap, (blk, i, j))
        end
    end
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer,
                                             set::SupportedSets)
    offset = length(optimizer.varmap)
    new_block(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.load_constrained_variables(
    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    set::SupportedSets)
end

function MOIU.load_variables(optimizer::Optimizer, nvars)
    @assert nvars == length(optimizer.varmap)
    dummy = isempty(optimizer.b)
    if dummy
        # See https://github.com/coin-or/Csdp/issues/2
        optimizer.b = [one(Cdouble)]
        optimizer.blockdims = [optimizer.blockdims; CSDP_INT(-1)]
        count_entry(optimizer, 1, length(optimizer.blockdims))
    end
    optimizer.C.nblocks = length(optimizer.blockdims)
    num_entries = zeros(CSDP_INT, length(optimizer.b), length(optimizer.blockdims))
    for (key, value) in optimizer.num_entries
        num_entries[key...] = value
    end
    optimizer.problem = allocate_loading_prob(Ref(optimizer.C), optimizer.blockdims, length(optimizer.b), num_entries, 3)
    if dummy
        # See https://github.com/coin-or/Csdp/issues/2
        duplicate = addentry(optimizer.problem, 1, length(optimizer.blockdims), 1, 1, 1.0, true)
        @assert !duplicate
    end

end

function count_entry(optimizer::Optimizer, con_idx::Integer, blk::Integer)
    key = (con_idx, blk)
    optimizer.num_entries[key] = get(optimizer.num_entries, key, 0) + 1
end

function MOIU.allocate_constraint(optimizer::Optimizer,
                                  func::MOI.ScalarAffineFunction{Cdouble},
                                  set::MOI.EqualTo{Cdouble})
    if !iszero(MOI.constant(func))
        throw(MOI.ScalarFunctionConstantNotZero{
            Cdouble, MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}(
                MOI.constant(func)))
    end
    push!(optimizer.b, MOI.constant(set))
    func = MOIU.canonical(func) # sum terms with same variables and same output_index
    for t in func.terms
        if !iszero(t.coefficient)
            blk, i, j = varmap(optimizer, t.variable_index)
            count_entry(optimizer, length(optimizer.b), blk)
        end
    end
    return AFFEQ(length(optimizer.b))
end

function MOIU.load_constraint(optimizer::Optimizer, ci::AFFEQ,
                              f::MOI.ScalarAffineFunction, s::MOI.EqualTo)
    if !iszero(MOI.constant(f))
        throw(MOI.ScalarFunctionConstantNotZero{
            Cdouble, MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}(
                MOI.constant(f)))
    end
    setconstant(optimizer.problem, ci.value, MOI.constant(s))
    f = MOIU.canonical(f) # sum terms with same variables and same output_index
    if isempty(f.terms)
        throw(ArgumentError("Empty constraint $ci: $f-in-$s. Not supported by CSDP."))
    end
    for t in f.terms
        if !iszero(t.coefficient)
            blk, i, j = varmap(optimizer, t.variable_index)
            coef = t.coefficient
            if i != j
                coef /= 2
            end
            duplicate = addentry(optimizer.problem, ci.value, blk, i, j, coef, true)
            @assert !duplicate
        end
    end
end

function MOI.optimize!(optimizer::Optimizer)
    write_prob(optimizer)

    start_time = time()
    optimizer.y = loaded_initsoln(optimizer.problem, length(optimizer.b), Ref(optimizer.X), Ref(optimizer.Z))

    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:printlevel] = 0
    end

    optimizer.status, optimizer.pobj, optimizer.dobj = loaded_sdp(
        optimizer.problem, Ref(optimizer.X), optimizer.y,
        Ref(optimizer.Z), options)
    optimizer.solve_time = time() - start_time
end

function MOI.get(m::Optimizer, ::MOI.TerminationStatus)
    status = m.status
    if status == -1
        return MOI.OPTIMIZE_NOT_CALLED
    elseif status == 0
        return MOI.OPTIMAL
    elseif status == 1
        return MOI.INFEASIBLE
    elseif status == 2
        return MOI.DUAL_INFEASIBLE
    elseif status == 3
        return MOI.ALMOST_OPTIMAL
    elseif status == 4
        return MOI.ITERATION_LIMIT
    elseif 5 <= status <= 7
        return MOI.SLOW_PROGRESS
    elseif 8 <= status <= 9
        return MOI.NUMERICAL_ERROR
    else
        error("Internal library error: status=$status")
    end
end

function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = m.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 1
        return MOI.INFEASIBLE_POINT
    elseif status == 2
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == 3
        return MOI.NEARLY_FEASIBLE_POINT
    elseif 4 <= status <= 9
        return MOI.UNKNOWN_RESULT_STATUS
    else
        error("Internal library error: status=$status")
    end
end

function MOI.get(m::Optimizer, attr::MOI.DualStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = m.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 1
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == 2
        return MOI.INFEASIBLE_POINT
    elseif status == 3
        return MOI.NEARLY_FEASIBLE_POINT
    elseif 4 <= status <= 9
        return MOI.UNKNOWN_RESULT_STATUS
    else
        error("Internal library error: status=$status")
    end
end

MOI.get(m::Optimizer, ::MOI.ResultCount) = 1
function MOI.get(m::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objsign * m.pobj + m.objconstant
end
function MOI.get(m::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objsign * m.dobj + m.objconstant
end
struct PrimalSolutionMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::PrimalSolutionMatrix) = true
MOI.get(optimizer::Optimizer, ::PrimalSolutionMatrix) = optimizer.X

struct DualSolutionVector <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSolutionVector) = true
MOI.get(optimizer::Optimizer, ::DualSolutionVector) = optimizer.y

struct DualSlackMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSlackMatrix) = true
MOI.get(optimizer::Optimizer, ::DualSlackMatrix) = optimizer.Z

function block(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    return optimizer.varmap[ci.value][1]
end
function dimension(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    blockdim = optimizer.blockdims[block(optimizer, ci)]
    if blockdim < 0
        return -blockdim
    else
        return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(blockdim))
    end
end
function vectorize_block(M, blk::Integer, s::Type{MOI.Nonnegatives})
    return diag(block(M, blk))
end
function vectorize_block(M::AbstractMatrix{Cdouble}, blk::Integer, s::Type{MOI.PositiveSemidefiniteConeTriangle}) where T
    B = block(M, blk)
    d = LinearAlgebra.checksquare(B)
    n = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))
    v = Vector{Cdouble}(undef, n)
    k = 0
    for j in 1:d
        for i in 1:j
            k += 1
            v[k] = B[i, j]
        end
    end
    @assert k == n
    return v
end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    blk, i, j = varmap(optimizer, vi)
    return block(MOI.get(optimizer, PrimalSolutionMatrix()), blk)[i, j]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, PrimalSolutionMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.b[ci.value]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, DualSlackMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return -MOI.get(optimizer, DualSolutionVector())[ci.value]
end


text:
using ConstrainedDynamics
using ConstrainedDynamics: posargsk, fullargssol, Lmat, VLᵀmat, Lᵀmat, ωbar
using Rotations
using Rotations: params
using StaticArrays

function getxqkvector(state)
    x, q = posargsk(state)
    q = params(q)

    return [x;q]
end

function getsolestimate(state)
    x, v, q, ω = fullargssol(state)
    q = params(q)

    return [x;v;q;ω]
end

getxk(state) = posargsk(state)[1]
getqk(state) = posargsk(state)[2]

function transfuncpos(vars)
    xa = vars[1:3]
    qa = UnitQuaternion(SA[vars[4:7]...])
    xb = vars[8:10]
    qb = UnitQuaternion(SA[vars[11:14]...])

    pa = vars[15:17]
    pb = vars[18:20]

    return vrotate(xb + vrotate(pb, qb) - (xa + vrotate(pa, qa)), inv(qa))
end

function transfunc3pos(vars)
    g = transfuncpos(vars)

    return g
end

function transfunc2pos(vars)
    g = transfuncpos(vars)
    V1 = vars[21:23]
    V2 = vars[24:26]
    V12 = [V1';V2']

    return V12 * g
end

function transfunc1pos(vars)
    g = transfuncpos(vars)
    V3 = vars[21:23]

    return V3' * g
end

function transfuncvel(vars)
    Δt = 0.01

    xa1 = vars[1:3]
    va1 = vars[4:6]
    qa1 = UnitQuaternion(SA[vars[7:10]...])
    ωa1 = vars[11:13]

    xb1 = vars[14:16]
    vb1 = vars[17:19]
    qb1 = UnitQuaternion(SA[vars[20:23]...])
    ωb1 = vars[24:26]

    xa2 = xa1 + va1 * Δt
    qa2 = qa1 * ωbar(ωa1, Δt) * Δt / 2

    xb2 = xb1 + vb1 * Δt
    qb2 = qb1 * ωbar(ωb1, Δt) * Δt / 2

    pa = vars[27:29]
    pb = vars[30:32]

    return vrotate(xb2 + vrotate(pb, qb2) - (xa2 + vrotate(pa, qa2)), inv(qa2))
end

function transfunc3vel(vars)
    g = transfuncvel(vars)

    return g
end

function transfunc2vel(vars)
    g = transfuncvel(vars)
    V1 = vars[33:35]
    V2 = vars[36:38]
    V12 = [V1';V2']

    return V12 * g
end

function transfunc1vel(vars)
    g = transfuncvel(vars)
    V3 = vars[33:35]

    return V3' * g
end

function rotfuncpos(vars)
    xa = vars[1:3]
    qa = UnitQuaternion(SA[vars[4:7]...], false)
    xb = vars[8:10]
    qb = UnitQuaternion(SA[vars[11:14]...], false)

    offset = UnitQuaternion(SA[vars[15:18]...], false)

    return VLᵀmat(qa) * Lmat(qb) * params(inv(offset))
end

function rotfunc3pos(vars)
    g = rotfuncpos(vars)

    return g
end

function rotfunc2pos(vars)
    g = rotfuncpos(vars)
    V1 = vars[19:21]
    V2 = vars[22:24]
    V12 = [V1';V2']

    return V12 * g
end

function rotfunc1pos(vars)
    g = rotfuncpos(vars)
    V3 = vars[19:21]

    return V3' * g
end

function rotfuncvel(vars)
    Δt = 0.01

    xa1 = vars[1:3]
    va1 = vars[4:6]
    qa1 = UnitQuaternion(SA[vars[7:10]...], false)
    ωa1 = vars[11:13]

    xb1 = vars[14:16]
    vb1 = vars[17:19]
    qb1 = UnitQuaternion(SA[vars[20:23]...], false)
    ωb1 = vars[24:26]

    xa2 = xa1 + va1 * Δt
    qa2 = qa1 * ωbar(ωa1, Δt) * Δt / 2

    xb2 = xb1 + vb1 * Δt
    qb2 = qb1 * ωbar(ωb1, Δt) * Δt / 2

    offset = UnitQuaternion(SA[vars[27:30]...], false)

    return VLᵀmat(qa2) * Lmat(qb2) * params(inv(offset))
end

function rotfunc3vel(vars)
    g = rotfuncvel(vars)

    return g
end

function rotfunc2vel(vars)
    g = rotfuncvel(vars)
    V1 = vars[31:33]
    V2 = vars[34:36]
    V12 = [V1';V2']

    return V12 * g
end

function rotfunc1vel(vars)
    g = rotfuncvel(vars)
    V3 = vars[31:33]

    return V3' * g
end


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

function mark_kross(r)
    side = Up
    for i in 0:1:3
        steps = 0
        while isborder(r,side)==false
            move!(r,side)
            putmarker!(r)
            steps += 1
        end
        movements!(r, inverse(side), steps)
        side = left(side)
    end
    putmarker!(r)
end


text:
module TestSetOprations

using FactCheck
using DataCubes

facts("Set Operations tests") do
  context("unique test") do
    @fact unique(nalift([1,2,3,3,2,5]),1) --> nalift([1,2,3,5])
    @fact unique(darr(a=[1,2,3,3,2,5]),1) --> darr(a=[1,2,3,5])
    @fact unique(larr(a=[1,2,3,3,2,5]),1) --> larr(a=[1,2,3,5])
    @fact unique(nalift([1 1 2 1;1 2 4 1;1 1 2 1]),1,2) --> @nalift [1 NA NA;NA 2 4]
    @fact unique(nalift([1 1 2 1;1 2 4 1;1 1 2 1]),2,1) --> @nalift [1 2;NA 4]
    @fact unique(nalift([1 1 2 1;1 2 4 1;1 1 2 1]),1) --> @nalift [1 1 2 1;1 2 4 1]
    @fact unique(nalift([1 1 2 1;1 2 4 1;1 1 2 1]),2) --> @nalift [1 1 2;1 2 4;1 1 2]
    @fact unique(@darr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),1,2) --> @darr(a=[1 NA NA;NA 2 4])
    @fact unique(@darr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),2,1) --> @darr(a=[1 2;NA 4])
    @fact unique(@darr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),1) --> @darr(a=[1 1 2 1;1 2 4 1])
    @fact unique(@darr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),2) --> @darr(a=[1 1 2;1 2 4;1 1 2])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),1,2) --> @larr(a=[1 NA NA;NA 2 4])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),2,1) --> @larr(a=[1 2;NA 4])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),1) --> @larr(a=[1 1 2 1;1 2 4 1])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1]),2) --> @larr(a=[1 1 2;1 2 4;1 1 2])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1],axis1[k=[:a,:b,:c]]),1,2) --> @larr(a=[1 NA NA;NA 2 4],axis1[k=[:a,:b]])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1],axis1[k=[:a,:b,:c]]),2,1) --> @larr(a=[1 2;NA 4],axis1[k=[:a,:b]])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1],axis1[k=[:a,:b,:c]]),1) --> @larr(a=[1 1 2 1;1 2 4 1],axis1[k=[:a,:b]])
    @fact unique(@larr(a=[1 1 2 1;1 2 4 1;1 1 2 1],axis1[k=[:a,:b,:c]]),2) --> @larr(a=[1 1 2;1 2 4;1 1 2],axis1[k=[:a,:b,:c]])
  end
  context("union test") do
    @fact union(1, nalift([1,2,3]), nalift([3,4,5])) --> nalift([1,2,3,4,5])
    @fact union(1, darr(a=[1,2,3]), darr(a=[3,4,5])) --> darr(a=[1,2,3,4,5])
    @fact union(1, larr(a=[1,2,3]), larr(a=[3,4,5])) --> larr(a=[1,2,3,4,5])
    @fact union(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[:sym1 :sym2 :sym3])) --> @larr(a=[1 2 3;4 5 6;:sym1 :sym2 :sym3])
    @fact union(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2 3])) --> @larr(a=[1 2 3;4 5 6])
    @fact union(2, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2;4 7])) --> @larr(a=[1 2 3 2;4 5 6 7])
    @fact union(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 7],b=[11 12;14 15])) --> @larr(a=[1 2 3 2;4 5 6 7], b=[11 12 13 12;14 15 16 15])
    @fact union(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 5],b=[11 12;14 15])) --> @larr(a=[1 2 3;4 5 6], b=[11 12 13;14 15 16])
    @fact union(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r5,:r2]])) --> @larr(a=[1 2 3 1 2;4 5 6 4 6],b=[11 12 13 11 12;14 15 16 14 15],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3,:r5,:r2]])
    @fact union(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r1,:r2]])) --> @larr(a=[1 2 3 2;4 5 6 6],b=[11 12 13 12;14 15 16 15],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3,:r2]])

    @fact union(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[:sym1 :sym2 :sym3])) --> @darr(a=[1 2 3;4 5 6;:sym1 :sym2 :sym3])
    @fact union(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2 3])) --> @darr(a=[1 2 3;4 5 6])
    @fact union(2, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2;4 7])) --> @darr(a=[1 2 3 2;4 5 6 7])
    @fact union(2, @darr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @darr(a=[1 2;4 7],b=[11 12;14 15])) --> @darr(a=[1 2 3 2;4 5 6 7], b=[11 12 13 12;14 15 16 15])
    @fact union(1, @nalift([1 2 3;4 5 6]), @nalift([:sym1 :sym2 :sym3])) --> @nalift([1 2 3;4 5 6;:sym1 :sym2 :sym3])
    @fact union(1, @nalift([1 2 3;4 5 6]), @nalift([1 2 3])) --> @nalift([1 2 3;4 5 6])
    @fact union(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 7])) --> @nalift([1 2 3 2;4 5 6 7])
    @fact union(2, @nalift([1 2 3;4 5 6]), @nalift([11 12;14 15])) --> @nalift([1 2 3 11 12;4 5 6 14 15])
    @fact union(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 5])) --> @nalift([1 2 3;4 5 6])
  end
  context("intersect test") do
    @fact intersect(1, larr(a=[1,2,3]), larr(a=[3,4,5])) --> larr(a=[3])
    @fact intersect(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[:sym1 :sym2 :sym3])) --> isempty
    @fact intersect(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2 3])) --> @larr(a=[1 2 3])
    @fact intersect(2, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2;4 7])) --> @larr(a=[1 4]')
    @fact intersect(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 7],b=[11 12;14 15])) --> @larr(a=[1 4]',b=[11 14]')
    @fact intersect(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 5],b=[11 12;14 15])) --> @larr(a=[1 2;4 5],b=[11 12;14 15])
    @fact intersect(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r5,:r2]])) --> isempty
    @fact intersect(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r1,:r2]])) --> @larr(a=[1 4]',b=[11 14]',axis1[k=['a','b']],axis2[r=[:r1]])

    @fact intersect(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[:sym1 :sym2 :sym3])) --> isempty
    @fact intersect(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2 3])) --> @darr(a=[1 2 3])
    @fact intersect(2, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2;4 7])) --> @darr(a=[1 4]')
    @fact intersect(2, @darr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @darr(a=[1 2;4 7],b=[11 12;14 15])) --> @darr(a=[1 4]', b=[11 14]')
    @fact intersect(1, @nalift([1 2 3;4 5 6]), @nalift([:sym1 :sym2 :sym3])) --> isempty
    @fact intersect(1, @nalift([1 2 3;4 5 6]), @nalift([1 2 3])) --> @nalift([1 2 3])
    @fact intersect(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 7])) --> @nalift([1 4]')
    @fact intersect(2, @nalift([1 2 3;4 5 6]), @nalift([11 12;14 15])) --> isempty
    @fact intersect(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 5])) --> @nalift([1 2;4 5])
  end
  context("setdiff test") do
    @fact setdiff(1, nalift([1,2,3]), nalift([3,4,5])) --> nalift([1,2])
    @fact setdiff(1, darr(a=[1,2,3]), darr(a=[3,4,5])) --> darr(a=[1,2])
    @fact setdiff(1, larr(a=[1,2,3]), larr(a=[3,4,5])) --> larr(a=[1,2])
    @fact setdiff(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[:sym1 :sym2 :sym3])) --> setdiff(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[:sym1 :sym2 :sym3]))
    @fact setdiff(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2 3;4 5 6])) --> isempty
    @fact setdiff(1, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2 3])) --> @larr(a=[4 5 6])
    @fact setdiff(2, @larr(a=[1 2 3;4 5 6]), @larr(a=[1 2;4 7])) --> @larr(a=[2 3;5 6])
    @fact setdiff(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 7],b=[11 12;14 15])) --> @larr(a=[2 3;5 6],b=[12 13;15 16])
    @fact setdiff(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @larr(a=[1 2;4 5],b=[11 12;14 15])) --> @larr(a=[3 6]',b=[13 16]')
    @fact setdiff(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r5,:r2]])) --> @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]])
    @fact setdiff(2, @larr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16],axis1[k=['a','b']],axis2[r=[:r1,:r2,:r3]]), @larr(a=[1 2;4 6],b=[11 12;14 15],axis1[k=['a','b']],axis2[r=[:r1,:r2]])) --> @larr(a=[2 3;5 6],b=[12 13;15 16],axis1[k=['a','b']],axis2[r=[:r2,:r3]])

    @fact setdiff(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[:sym1 :sym2 :sym3])) --> @darr(a=[1 2 3;4 5 6])
    @fact setdiff(1, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2 3])) --> @darr(a=[4 5 6])
    @fact setdiff(2, @darr(a=[1 2 3;4 5 6]), @darr(a=[1 2;4 7])) --> @darr(a=[2 3;5 6])
    @fact setdiff(2, @darr(a=[1 2 3;4 5 6],b=[11 12 13;14 15 16]), @darr(a=[1 2;4 7],b=[11 12;14 15])) --> @darr(a=[2 3;5 6],b=[12 13;15 16])
    @fact setdiff(1, @nalift([1 2 3;4 5 6]), @nalift([:sym1 :sym2 :sym3])) --> @nalift([1 2 3;4 5 6])
    @fact setdiff(1, @nalift([1 2 3;4 5 6]), @nalift([1 2 3])) --> @nalift([4 5 6])
    @fact setdiff(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 7])) --> @nalift([2 3;5 6])
    @fact setdiff(2, @nalift([1 2 3;4 5 6]), @nalift([11 12;14 15])) --> @nalift([1 2 3;4 5 6])
    @fact setdiff(2, @nalift([1 2 3;4 5 6]), @nalift([1 2;4 5])) --> @nalift([3 6]')
  end
end

end