text:
type Header
	tracenum::Int32
	o1::Float32
	n1::Int32
	d1::Float32
	sx::Float32
	sy::Float32
	gx::Float32
	gy::Float32
	mx::Float32
	my::Float32
	hx::Float32
	hy::Float32
	h::Float32
	az::Float32
	ang::Float32
	isx::Int32
	isy::Int32
	igx::Int32
	igy::Int32
	imx::Int32
	imy::Int32
	ihx::Int32
	ihy::Int32
	ih::Int32
	iaz::Int32
	iang::Int32
	selev::Float32
	gelev::Float32
	sstat::Float32
	gstat::Float32
	trid::Int32
end

header_count = Dict{AbstractString,Int32}()
header_count["tracenum"] = 0
header_count["o1"]    = 4
header_count["n1"]    = 8
header_count["d1"]    = 12
header_count["sx"]    = 16
header_count["sy"]    = 20
header_count["gx"]    = 24
header_count["gy"]    = 28
header_count["mx"]    = 32
header_count["my"]    = 36
header_count["hx"]    = 40
header_count["hy"]    = 44
header_count["h"]     = 48
header_count["az"]    = 52
header_count["ang"]   = 56
header_count["isx"]   = 60
header_count["isy"]   = 64
header_count["igx"]   = 68
header_count["igy"]   = 72
header_count["imx"]   = 76
header_count["imy"]   = 80
header_count["ihx"]   = 84
header_count["ihy"]   = 88
header_count["ih"]    = 92
header_count["iaz"]   = 96
header_count["iang"]  = 100
header_count["selev"] = 104
header_count["gelev"] = 108
header_count["sstat"] = 112
header_count["gstat"] = 116
header_count["trid"]  = 120

function InitSeisHeader()
	h = Header(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
	0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
	0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
	0.0)
	return h
end

function GrabHeader(stream,j)
	@compat position = 4*length(fieldnames(Header))*(j-1)
	seek(stream,position)
	h = InitSeisHeader()
	h.tracenum = read(stream,typeof(h.tracenum))
	h.o1    = read(stream,typeof(h.o1))
	h.n1    = read(stream,typeof(h.n1))
	h.d1    = read(stream,typeof(h.d1))
	h.sx    = read(stream,typeof(h.sx))
	h.sy    = read(stream,typeof(h.sy))
	h.gx    = read(stream,typeof(h.gx))
	h.gy    = read(stream,typeof(h.gy))
	h.mx    = read(stream,typeof(h.mx))
	h.my    = read(stream,typeof(h.my))
	h.hx    = read(stream,typeof(h.hx))
	h.hy    = read(stream,typeof(h.hy))
	h.h     = read(stream,typeof(h.h))
	h.az    = read(stream,typeof(h.az))
	h.ang   = read(stream,typeof(h.ang))
	h.isx   = read(stream,typeof(h.isx))
	h.isy   = read(stream,typeof(h.isy))
	h.igx   = read(stream,typeof(h.igx))
	h.igy   = read(stream,typeof(h.igy))
	h.imx   = read(stream,typeof(h.imx))
	h.imy   = read(stream,typeof(h.imy))
	h.ihx   = read(stream,typeof(h.ihx))
	h.ihy   = read(stream,typeof(h.ihy))
	h.ih    = read(stream,typeof(h.ih))
	h.iaz   = read(stream,typeof(h.iaz))
	h.iang   = read(stream,typeof(h.iang))
	h.selev = read(stream,typeof(h.selev))
	h.gelev = read(stream,typeof(h.gelev))
	h.sstat = read(stream,typeof(h.sstat))
	h.gstat = read(stream,typeof(h.gstat))
	h.trid  = read(stream,typeof(h.trid))
	return h
end

function PutHeader(stream,h,j)
	@compat position = 4*length(fieldnames(Header))*(j-1)
	seek(stream,position)
	write(stream,h.tracenum)
	write(stream,h.o1)
	write(stream,h.n1)
	write(stream,h.d1)
	write(stream,h.sx)
	write(stream,h.sy)
	write(stream,h.gx)
	write(stream,h.gy)
	write(stream,h.mx)
	write(stream,h.my)
	write(stream,h.hx)
	write(stream,h.hy)
	write(stream,h.h)
	write(stream,h.az)
	write(stream,h.ang)
	write(stream,h.isx)
	write(stream,h.isy)
	write(stream,h.igx)
	write(stream,h.igy)
	write(stream,h.imx)
	write(stream,h.imy)
	write(stream,h.ihx)
	write(stream,h.ihy)
	write(stream,h.ih)
	write(stream,h.iaz)
	write(stream,h.iang)
	write(stream,h.selev)
	write(stream,h.gelev)
	write(stream,h.sstat)
	write(stream,h.gstat)
	write(stream,h.trid)
end

bitstype 32 Header32Bits

function BitsToHeader(h_in)
	h = InitSeisHeader()
	h.tracenum = reinterpret(typeof(h.tracenum),h_in[1])
	h.o1    = reinterpret(typeof(h.o1),h_in[2])
	h.n1    = reinterpret(typeof(h.n1),h_in[3])
	h.d1    = reinterpret(typeof(h.d1),h_in[4])
	h.sx    = reinterpret(typeof(h.sx),h_in[5])
	h.sy    = reinterpret(typeof(h.sy),h_in[6])
	h.gx    = reinterpret(typeof(h.gx),h_in[7])
	h.gy    = reinterpret(typeof(h.gy),h_in[8])
	h.mx    = reinterpret(typeof(h.mx),h_in[9])
	h.my    = reinterpret(typeof(h.my),h_in[10])
	h.hx    = reinterpret(typeof(h.hx),h_in[11])
	h.hy    = reinterpret(typeof(h.hy),h_in[12])
	h.h     = reinterpret(typeof(h.h),h_in[13])
	h.az    = reinterpret(typeof(h.az),h_in[14])
	h.ang   = reinterpret(typeof(h.ang),h_in[15])
	h.isx   = reinterpret(typeof(h.isx),h_in[16])
	h.isy   = reinterpret(typeof(h.isy),h_in[17])
	h.igx   = reinterpret(typeof(h.igx),h_in[18])
	h.igy   = reinterpret(typeof(h.igy),h_in[19])
	h.imx   = reinterpret(typeof(h.imx),h_in[20])
	h.imy   = reinterpret(typeof(h.imy),h_in[21])
	h.ihx   = reinterpret(typeof(h.ihx),h_in[22])
	h.ihy   = reinterpret(typeof(h.ihy),h_in[23])
	h.ih    = reinterpret(typeof(h.ih),h_in[24])
	h.iaz   = reinterpret(typeof(h.iaz),h_in[25])
	h.iang   = reinterpret(typeof(h.iang),h_in[26])
	h.selev = reinterpret(typeof(h.selev),h_in[27])
	h.gelev = reinterpret(typeof(h.gelev),h_in[28])
	h.sstat = reinterpret(typeof(h.sstat),h_in[29])
	h.gstat = reinterpret(typeof(h.gstat),h_in[30])
	h.trid  = reinterpret(typeof(h.trid),h_in[31])
	return h
end

function HeaderToBits(h);
	h_out = [reinterpret(Header32Bits,h.tracenum);
	reinterpret(Header32Bits,h.o1);
	reinterpret(Header32Bits,h.n1);
	reinterpret(Header32Bits,h.d1);
	reinterpret(Header32Bits,h.sx);
	reinterpret(Header32Bits,h.sy);
	reinterpret(Header32Bits,h.gx);
	reinterpret(Header32Bits,h.gy);
	reinterpret(Header32Bits,h.mx);
	reinterpret(Header32Bits,h.my);
	reinterpret(Header32Bits,h.hx);
	reinterpret(Header32Bits,h.hy);
	reinterpret(Header32Bits,h.h);
	reinterpret(Header32Bits,h.az);
	reinterpret(Header32Bits,h.ang);
	reinterpret(Header32Bits,h.isx);
	reinterpret(Header32Bits,h.isy);
	reinterpret(Header32Bits,h.igx);
	reinterpret(Header32Bits,h.igy);
	reinterpret(Header32Bits,h.imx);
	reinterpret(Header32Bits,h.imy);
	reinterpret(Header32Bits,h.ihx);
	reinterpret(Header32Bits,h.ihy);
	reinterpret(Header32Bits,h.ih);
	reinterpret(Header32Bits,h.iaz);
	reinterpret(Header32Bits,h.iang);
	reinterpret(Header32Bits,h.selev);
	reinterpret(Header32Bits,h.gelev);
	reinterpret(Header32Bits,h.sstat);
	reinterpret(Header32Bits,h.gstat);
	reinterpret(Header32Bits,h.trid)]
	return h_out
end

function GetNumTraces(in)

	filename_h = ParseHeaderName(in)
	@compat nhead = length(fieldnames(Header))
	stream_h = open(filename_h)
	nx = round(Int,filesize(stream_h)/(nhead*4))
	close(stream_h)
	return nx

end

function ParseHeaderName(filename::ASCIIString)

	return success(`grep "headers=" $filename`) ? chomp(readall(pipeline(`grep "headers=" $filename` ,pipeline(`tail -1` , `awk '{print substr($1,10,length($1)-10) }' `))    )) : "NULL"

end

function ParseDataName(filename::ASCIIString)

	return chomp(readall(pipeline(`grep "in=" $filename` ,pipeline(`tail -1` , `awk '{print substr($1,5,length($1)-5) }' `))    ))

end

function ParseDataFormat(filename::ASCIIString)

	return chomp(readall(pipeline(`grep "data_format" $filename` , pipeline(`tail -1` , `awk '{print substr($1,14,length($1)-14)}'`)) ))

end

function ParseDataESize(filename::ASCIIString)

	@compat return parse(Int,chomp(readall(pipeline(`grep "esize" $filename` , pipeline(`tail -1` , `awk '{print substr($1,7,length($1))}'`)) )))

end

function ExtractHeader(h::Array{Header,1},key::ASCIIString)

	keytype = eval(parse("typeof(Seismic.InitSeisHeader().$(key))"))
	out = keytype[]
	for ix = 1 : length(h)
	push!(out,getfield(h[ix],symbol(key)))
	end
	return out

end

type Extent
        n1::Int32
        n2::Int32
        n3::Int32
        n4::Int32
        n5::Int32
        o1::Float32
        o2::Float32
        o3::Float32
        o4::Float32
        o5::Float32
        d1::Float32
        d2::Float32
        d3::Float32
        d4::Float32
        d5::Float32
        label1::ASCIIString
        label2::ASCIIString
        label3::ASCIIString
        label4::ASCIIString
        label5::ASCIIString
        unit1::ASCIIString
        unit2::ASCIIString
        unit3::ASCIIString
        unit4::ASCIIString
        unit5::ASCIIString
        title::ASCIIString
end

function ReadTextHeader(filename)
	@compat n1 = success(`grep "n1" $filename`) ? parse(Int,chomp(readall(pipeline(`grep "n1" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	@compat n2 = success(`grep "n2" $filename`) ? parse(Int,chomp(readall(pipeline(`grep "n2" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	@compat n3 = success(`grep "n3" $filename`) ? parse(Int,chomp(readall(pipeline(`grep "n3" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	@compat n4 = success(`grep "n4" $filename`) ? parse(Int,chomp(readall(pipeline(`grep "n4" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	@compat n5 = success(`grep "n5" $filename`) ? parse(Int,chomp(readall(pipeline(`grep "n5" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	o1 = success(`grep "o1" $filename`) ? float(chomp(readall(pipeline(`grep "o1" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 0
	o2 = success(`grep "o2" $filename`) ? float(chomp(readall(pipeline(`grep "o2" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 0
	o3 = success(`grep "o3" $filename`) ? float(chomp(readall(pipeline(`grep "o3" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 0
	o4 = success(`grep "o4" $filename`) ? float(chomp(readall(pipeline(`grep "o4" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 0
	o5 = success(`grep "o5" $filename`) ? float(chomp(readall(pipeline(`grep "o5" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 0
	d1 = success(`grep "d1" $filename`) ? float(chomp(readall(pipeline(`grep "d1" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	d2 = success(`grep "d2" $filename`) ? float(chomp(readall(pipeline(`grep "d2" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	d3 = success(`grep "d3" $filename`) ? float(chomp(readall(pipeline(`grep "d3" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	d4 = success(`grep "d4" $filename`) ? float(chomp(readall(pipeline(`grep "d4" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	d5 = success(`grep "d5" $filename`) ? float(chomp(readall(pipeline(`grep "d5" $filename` , pipeline(`tail -1` , `awk '{print substr($1,4,length($1))}'`)) ))) : 1
	label1 = success(`grep "label1" $filename`) ? chomp(readall(pipeline(`grep "label1" $filename` , pipeline(`tail -1` , `awk '{print substr($0,10,length($0)-10)}'`)) )) : ""
	label2 = success(`grep "label2" $filename`) ? chomp(readall(pipeline(`grep "label2" $filename` , pipeline(`tail -1` , `awk '{print substr($0,10,length($0)-10)}'`)) )) : ""
	label3 = success(`grep "label3" $filename`) ? chomp(readall(pipeline(`grep "label3" $filename` , pipeline(`tail -1` , `awk '{print substr($0,10,length($0)-10)}'`)) )) : ""
	label4 = success(`grep "label4" $filename`) ? chomp(readall(pipeline(`grep "label4" $filename` , pipeline(`tail -1` , `awk '{print substr($0,10,length($0)-10)}'`)) )) : ""
	label5 = success(`grep "label5" $filename`) ? chomp(readall(pipeline(`grep "label5" $filename` , pipeline(`tail -1` , `awk '{print substr($0,10,length($0)-10)}'`)) )) : ""
	unit1 = success(`grep "unit1" $filename`) ? chomp(readall(pipeline(`grep "unit1" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	unit2 = success(`grep "unit2" $filename`) ? chomp(readall(pipeline(`grep "unit2" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	unit3 = success(`grep "unit3" $filename`) ? chomp(readall(pipeline(`grep "unit3" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	unit4 = success(`grep "unit4" $filename`) ? chomp(readall(pipeline(`grep "unit4" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	unit5 = success(`grep "unit5" $filename`) ? chomp(readall(pipeline(`grep "unit5" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	title = success(`grep "title" $filename`) ? chomp(readall(pipeline(`grep "title" $filename` , pipeline(`tail -1` , `awk '{print substr($0,9,length($0)-9)}'`)) )) : ""
	extent = Extent(convert(Int32,n1),convert(Int32,n2),convert(Int32,n3),convert(Int32,n4),convert(Int32,n5),
		   convert(Float32,o1),convert(Float32,o2),convert(Float32,o3),convert(Float32,o4),convert(Float32,o5),
		   convert(Float32,d1),convert(Float32,d2),convert(Float32,d3),convert(Float32,d4),convert(Float32,d5),
		   label1,label2,label3,label4,label5,
		   unit1,unit2,unit3,unit4,unit5,
		   title)
	return extent
end

function WriteTextHeader(filename,extent,format,esize,filename_d,filename_h)
	# write the text header
	stream = open(filename, "w")
	write(stream,join(["	n1=",extent.n1,"\n"]))
	write(stream,join(["	n2=",extent.n2,"\n"]))
	write(stream,join(["	n3=",extent.n3,"\n"]))
	write(stream,join(["	n4=",extent.n4,"\n"]))
	write(stream,join(["	n5=",extent.n5,"\n"]))
	write(stream,join(["	o1=",extent.o1,"\n"]))
	write(stream,join(["	o2=",extent.o2,"\n"]))
	write(stream,join(["	o3=",extent.o3,"\n"]))
	write(stream,join(["	o4=",extent.o4,"\n"]))
	write(stream,join(["	o5=",extent.o5,"\n"]))
	write(stream,join(["	d1=",extent.d1,"\n"]))
	write(stream,join(["	d2=",extent.d2,"\n"]))
	write(stream,join(["	d3=",extent.d3,"\n"]))
	write(stream,join(["	d4=",extent.d4,"\n"]))
	write(stream,join(["	d5=",extent.d5,"\n"]))
	write(stream,join(["	label1=\"",extent.label1,"\"\n"]))
	write(stream,join(["	label2=\"",extent.label2,"\"\n"]))
	write(stream,join(["	label3=\"",extent.label3,"\"\n"]))
	write(stream,join(["	label4=\"",extent.label4,"\"\n"]))
	write(stream,join(["	label5=\"",extent.label5,"\"\n"]))
	write(stream,join(["	unit1=\"",extent.unit1,"\"\n"]))
	write(stream,join(["	unit2=\"",extent.unit2,"\"\n"]))
	write(stream,join(["	unit3=\"",extent.unit3,"\"\n"]))
	write(stream,join(["	unit4=\"",extent.unit4,"\"\n"]))
	write(stream,join(["	unit5=\"",extent.unit5,"\"\n"]))
	write(stream,join(["	title=\"",extent.title,"\"\n"]))
	write(stream,join(["	data_format=\"",format,"\"\n"]))
	write(stream,join(["	esize=",esize,"\n"]))
	write(stream,join(["	in=\"",filename_d,"\"\n"]))
	write(stream,join(["	headers=\"",filename_h,"\"\n"]))
	close(stream)
end


text:
## Flow field parameters

"""
    abstract type FlowFields

Data structure that provide access to flow fields (gridded as arrays) which will be
used to interpolate velocities to individual locations later on (once embedded in
an `Individuals` struct).

Following the C-grid convention also used in `MITgcm` (https://mitgcm.readthedocs.io)
flow fields are expected to be staggered as follows: grid cell i,j has its center located at i-1/2,j-1/2 while the
corresponding `u[i,j]` (resp. `v[i,j]) is located at i-1,j-1/2 (resp. i-1/2,j-1).

Also by convention, velocity fields are expected to have been normalized to grid units (e.g. 1/s rather than m/s)
before sending them to one of the supported `FlowFields` constructors (using either `Array` or `MeshArray`):

```
𝐹_Array2D (u0,u1,v0,v1,𝑇)
𝐹_Array3D (u0,u1,v0,v1,w0,w1,𝑇)
𝐹_MeshArray2D (u0,u1,v0,v1,𝑇,update_location!)
𝐹_MeshArray3D (u0,u1,v0,v1,w0,w1,𝑇,update_location!)
```

For example, constructor calls may look like

```
𝐹=𝐹_Array3D{eltype(u)}(u,u,v,v,0*w,1*w,[0.0,10.0])
or
𝐹=𝐹_MeshArray2D{eltype(u)}(u,u,v,v,[0.0,10.0],func)
```

as shown in the online documentation examples.

"""
abstract type FlowFields end

struct 𝐹_Array2D{T} <: FlowFields
    u0::Array{T,2}
    u1::Array{T,2}
    v0::Array{T,2}
    v1::Array{T,2}
    𝑇::Array{T}
end

struct 𝐹_Array3D{T} <: FlowFields
    u0::Array{T,3}
    u1::Array{T,3}
    v0::Array{T,3}
    v1::Array{T,3}
    w0::Array{T,3}
    w1::Array{T,3}
    𝑇::Array{T}
end

struct 𝐹_MeshArray2D{T} <: FlowFields
    u0::AbstractMeshArray{T,1}
    u1::AbstractMeshArray{T,1}
    v0::AbstractMeshArray{T,1}
    v1::AbstractMeshArray{T,1}
    𝑇::Array{T}
    update_location!::Function
end

struct 𝐹_MeshArray3D{T} <: FlowFields
    u0::AbstractMeshArray{T,2}
    u1::AbstractMeshArray{T,2}
    v0::AbstractMeshArray{T,2}
    v1::AbstractMeshArray{T,2}
    w0::AbstractMeshArray{T,2}
    w1::AbstractMeshArray{T,2}
    𝑇::Array{T}
    update_location!::Function
end

"""
    defaults for Individuals constructor
"""

default_solver(prob) = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
a=fill(0.0,1,1)
default_flowfields = 𝐹_Array2D{Float64}(a,a,a,a,[0. 1.])
default_recorder = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[])
default_postproc = (x->x)

"""
    struct Individuals{T,N}

- Data:           📌 (position),   🔴(record), 🆔 (ID), 𝑃 (`FlowFields`)
- Functions:      🚄 (velocity),   ∫ (integration), 🔧(postprocessing)
- NamedTuples:    𝐷 (diagnostics),      𝑀 (metadata)

The velocity function 🚄 typically computes velocity at individual positions (📌 to start) within the
specified space-time domain by interpolating gridded variables (provided via 𝑃). Individual trajectories
are computed by integrating (∫) interpolated velocities through time. Normally, integration is done by
calling ∫! which updates 📌 at the end and records results in 🔴 via 🔧. Unicode cheatsheet:

- 📌=`\\:pushpin:<tab>`,          🔴=`\\:red_circle:<tab>`, 🆔=`\\:id:<tab>`
- 🚄=`\\:bullettrain_side:<tab>`, ∫=`\\int<tab>`,          🔧=`\\:wrench:<tab>`
- 𝑃=`\\itP<tab>`,                 𝐷=`\\itD<tab>`,           𝑀=`\\itM<tab>`

Simple constructors that use `FlowFields` to choose adequate defaults:

- Individuals(𝐹::𝐹_Array2D,x,y)
- Individuals(𝐹::𝐹_Array3D,x,y,z)
- Individuals(𝐹::𝐹_MeshArray2D,x,y,fid)
- Individuals(𝐹::𝐹_MeshArray3D,x,y,z,fid)

Further customization is achievable via keyword constructors:

```
df=DataFrame( ID=[], x=[], y=[], z=[], t = [])
𝐼=Individuals{Float64,2}(📌=zeros(3,10),🆔=1:10,🔴=deepcopy(df))
𝐼=Individuals(📌=zeros(3,2),🆔=collect(1:2),🔴=deepcopy(df))
```

Or via the plain text (or no-unicode) constructors:

```
df=DataFrame( ID=[], x=[], y=[], z=[], t = [])
I=(position=zeros(3,2),ID=1:2,record=deepcopy(df))
I=Individuals(I)
```
"""
Base.@kwdef struct Individuals{T,N}
   📌  ::Array{T,N} = Array{T,N}(undef, Tuple(Int.(zeros(1,N)))) #\:pushpin:<tab>
   🔴  ::DataFrame = similar(default_recorder) #\:red_circle:<tab>
   🆔   ::Array{Int,1} = Array{Int,1}(undef, 0) #\:id:<tab>
   🚄  ::Function = dxy_dt #\:bullettrain_side:<tab>
   ∫   ::Function = default_solver #\int<tab>
   🔧  ::Function = default_postproc #\:wrench:<tab>
   𝑃   ::FlowFields = default_flowfields #\itP<tab>
   𝐷   ::NamedTuple = NamedTuple() #\itD<tab>
   𝑀   ::NamedTuple = NamedTuple() #\itM<tab>
end

function Individuals(NT::NamedTuple)

    haskey(NT,:position) ? 📌=NT.position : 📌=Array{Float64,2}(undef, Tuple(Int.(zeros(1,2))))
    haskey(NT,:record) ? 🔴=NT.record : 🔴=similar(default_recorder)
    haskey(NT,:ID) ? 🆔=NT.ID : 🆔=collect(1:size(📌,2))
    haskey(NT,:velocity) ? 🚄=NT.velocity : 🚄=dxy_dt
    haskey(NT,:integration) ? ∫=NT.integration : ∫=default_solver
    haskey(NT,:postprocessing) ? 🔧=NT.postprocessing : 🔧=default_postproc
    haskey(NT,:parameters) ? 𝑃=NT.parameters : 𝑃=default_flowfields
    haskey(NT,:diagnostics) ? 𝐷=NT.diagnostics : 𝐷=NamedTuple()
    haskey(NT,:metadata) ? 𝑀=NT.metadata : 𝑀=NamedTuple()
    isa(📌,UnitRange) ? 📌=collect(📌) : nothing
    haskey(NT,:type) ? T=NT.type : T=eltype(📌)

    Individuals{T,ndims(📌)}(📌=📌,🔴=🔴,🆔=🆔,🚄=🚄,∫=∫,🔧=🔧,𝑃=𝑃,𝐷=𝐷,𝑀=𝑀)
end

function Individuals(𝐹::𝐹_Array2D,x,y, NT::NamedTuple = NamedTuple())
    📌=permutedims([[x[i];y[i]] for i in eachindex(x)])
    length(📌)==1 ? 📌=📌[1] : nothing
    T=eltype(📌)

    🔴 = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[])
    haskey(NT,:🔴) ? 🔴=NT.🔴 : nothing

    🔧 = postprocess_xy
    haskey(NT,:🔧) ? 🔧=NT.🔧 : nothing

    🆔=collect(1:size(📌,2))
    haskey(NT,:🆔) ? 🆔=NT.🆔 : nothing

    ∫=default_solver
    haskey(NT,:∫) ? ∫=NT.∫ : nothing

    Individuals{T,ndims(📌)}(𝑃=𝐹,📌=📌,🔴=🔴,🆔=🆔,🚄=dxy_dt,∫=∫,🔧=🔧)
end

function Individuals(𝐹::𝐹_Array3D,x,y,z, NT::NamedTuple = NamedTuple())
    📌=permutedims([[x[i];y[i];z[i]] for i in eachindex(x)])
    length(📌)==1 ? 📌=📌[1] : nothing
    T=eltype(📌)

    🔴 = DataFrame(ID=Int[], x=Float64[], y=Float64[], z=Float64[], t=Float64[])
    haskey(NT,:🔴) ? 🔴=NT.🔴 : nothing

    function 🔧(sol,𝑄::FlowFields;id=missing,𝑇=missing)
        df=postprocess_xy(sol,𝐹,id=id,𝑇=𝑇)
        z=sol[3,:]
        df.z=z[:]
        return df
    end
    haskey(NT,:🔧) ? 🔧=NT.🔧 : nothing

    🆔=collect(1:size(📌,2))
    haskey(NT,:🆔) ? 🆔=NT.🆔 : nothing

    ∫=default_solver
    haskey(NT,:∫) ? ∫=NT.∫ : nothing

    Individuals{T,ndims(📌)}(𝑃=𝐹,📌=📌,🔴=🔴,🆔=🆔,🚄=dxyz_dt,∫=∫,🔧=🔧)
end

function Individuals(𝐹::𝐹_MeshArray2D,x,y,fid, NT::NamedTuple = NamedTuple())
    📌=permutedims([[x[i];y[i];fid[i]] for i in eachindex(x)])
    length(📌)==1 ? 📌=📌[1] : nothing
    T=eltype(📌)

    🔴 = DataFrame(ID=Int[], x=Float64[], y=Float64[], fid=Int64[], t=Float64[])
    haskey(NT,:🔴) ? 🔴=NT.🔴 : nothing

    🔧 = postprocess_MeshArray
    haskey(NT,:🔧) ? 🔧=NT.🔧 : nothing

    🆔=collect(1:size(📌,2))
    haskey(NT,:🆔) ? 🆔=NT.🆔 : nothing

    ∫=default_solver
    haskey(NT,:∫) ? ∫=NT.∫ : nothing

    Individuals{T,ndims(📌)}(𝑃=𝐹,📌=📌,🔴=🔴,🆔=🆔,🚄=dxy_dt!,∫=∫,🔧=🔧)
end

function Individuals(𝐹::𝐹_MeshArray3D,x,y,z,fid, NT::NamedTuple = NamedTuple())
    📌=permutedims([[x[i];y[i];z[i];fid[i]] for i in eachindex(x)])
    length(📌)==1 ? 📌=📌[1] : nothing
    T=eltype(📌)

    🔴 = DataFrame(ID=Int[], x=Float64[], y=Float64[], z=Float64[], fid=Int64[], t=Float64[])
    haskey(NT,:🔴) ? 🔴=NT.🔴 : nothing

    function 🔧(sol,𝑄::FlowFields;id=missing,𝑇=missing)
        df=postprocess_MeshArray(sol,𝐹,id=id,𝑇=𝑇)
        z=[sol[1,i,j][1] for i in 1:size(sol,2), j in 1:size(sol,3)]
        df.z=z[:]
        return df
    end
    haskey(NT,:🔧) ? 🔧=NT.🔧 : nothing

    🆔=collect(1:size(📌,2))
    haskey(NT,:🆔) ? 🆔=NT.🆔 : nothing

    ∫=default_solver
    haskey(NT,:∫) ? ∫=NT.∫ : nothing

    Individuals{T,ndims(📌)}(𝑃=𝐹,📌=📌,🔴=🔴,🆔=🆔,🚄=dxyz_dt!,∫=∫,🔧=🔧)
end

"""
    ∫!(𝐼::Individuals,𝑇::Tuple)

Displace simulated individuals continuously through space over time period 𝑇 starting from position 📌.

- This is typically achieved by computing the cumulative integral of velocity experienced by each individual along its trajectory (∫ 🚄 dt).
- The current default is `solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)` but all solver options from the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package are available.
- After this, `∫!` is also equipped to postprocess results recorded into 🔴 via the 🔧 workflow, and the last step in `∫!` consists in updating 📌 to be ready for continuing in a subsequent call to `∫!`.
"""
function ∫!(𝐼::Individuals,𝑇::Tuple)
    @unpack 🚄,📌,𝑃, 🔧, 🆔, 🔴, ∫ = 𝐼

    prob = ODEProblem(🚄,📌, 𝑇 ,𝑃)
    sol = ∫(prob)

    tmp = 🔧(sol,𝑃, id=🆔, 𝑇=𝑇)

    isempty(🔴) ? np =0 : np=length(🆔)
    append!(🔴,tmp[np+1:end,:])

    nd=length(size(sol))
    nd==3 ? 📌[:,:] = deepcopy(sol[:,:,end]) : 📌[:] = deepcopy(sol[:,end])

end

"""
    ∫!(𝐼::Individuals)

Call ∫!(𝐼::Individuals,𝐼.𝑃.𝑇)
"""
∫!(𝐼::Individuals) = ∫!(𝐼::Individuals,(𝐼.𝑃.𝑇[1],𝐼.𝑃.𝑇[2]))

## Convenience Methods (size,show,similar)

Base.size(A::Individuals) = size(A.📌)

function Base.show(io::IO, 𝐼::Individuals) where {T}
    @unpack 🚄,📌,𝑃, 𝐷, 𝑀, 🔧, 🆔, 🔴, ∫ = 𝐼
    printstyled(io, "  📌 details     = ",color=:normal)
    printstyled(io, "$(size(📌)) $(typeof(𝐼).parameters[1])\n",color=:blue)
    printstyled(io, "  🔴 details     = ",color=:normal)
    printstyled(io, "$(size(🔴)) $(names(🔴))\n",color=:blue)
    printstyled(io, "  🆔 range       = ",color=:normal)
    printstyled(io, "$(extrema(🆔))\n",color=:blue)
    printstyled(io, "  🚄 function    = ",color=:normal)
    printstyled(io, "$(🚄)\n",color=:blue)
    printstyled(io, "  ∫  function    = ",color=:normal)
    printstyled(io, "$(∫)\n",color=:blue)
    printstyled(io, "  🔧 function    = ",color=:normal)
    printstyled(io, "$(🔧)\n",color=:blue)
    printstyled(io, "  𝑃  details     = ",color=:normal)
    printstyled(io, "$(fieldnames(typeof(𝑃)))\n",color=:blue)
  return
end

function Base.similar(𝐼::Individuals)
    @unpack 🚄,📌,𝑃, 𝐷, 𝑀, 🔧, 🆔, 🔴, ∫ = 𝐼
    T = typeof(𝐼).parameters[1]
    N = ndims(𝐼.📌)
    return Individuals{T,N}(📌=similar(📌),🔴=similar(🔴),🆔=similar(🆔),
                          🚄=🚄, ∫=∫, 🔧=🔧, 𝑃=𝑃, 𝐷=𝐷, 𝑀=𝑀)
end

"""
    Base.diff(𝐼::Individuals)

Difference in grid unit coordinates (dx,dy) between final and initial positions.
"""
function Base.diff(𝐼::Individuals)
    f(x)=last(x).-first(x)
    🔴_by_ID = groupby(𝐼.🔴, :ID)
    return combine(🔴_by_ID,nrow,:x => f => :dx,:y => f => :dy)
end

"""
    gcdist(𝐼::Individuals)

Great circle distance (gcd in radians) between final and initial positions.
"""
function gcdist(𝐼::Individuals)
    🔴_by_ID = groupby(𝐼.🔴, :ID)
    tmp = combine(🔴_by_ID,
    :lon => first => :lo1,:lon => last => :lo2,
    :lat => first => :la1,:lat => last => :la2)

    gcdist(lo1,lo2,la1,la2) = acos(sind(la1)*sind(la2)+cosd(la1)*cosd(la2)*cosd(lo1-lo2))
    tmp.gcd=[gcdist(tmp.lo1[i],tmp.lo2[i],tmp.la1[i],tmp.la2[i]) for i in 1:size(tmp,1)]
    return tmp
end

ext: jl


text:
function TA_dijkstra_shortest_paths(graph, travel_time, origin, start_node, end_node)

    no_node = nv(graph)
    no_arc = ne(graph)

    distmx = Inf*ones(no_node, no_node)
	for i in 1:no_arc
		distmx[start_node[i], end_node[i]] = travel_time[i]
	end

    state = dijkstra_shortest_paths(graph, origin, distmx)
    return state
end

function create_graph(start_node, end_node)
    @assert Base.length(start_node)==Base.length(end_node)

    no_node = max(maximum(start_node), maximum(end_node))
    no_arc = Base.length(start_node)

    graph = DiGraph(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin && origin != destination && current != 0
        parent = state.parents[current]

        # println("origin=$origin, destination=$destination, parent=$parent, current=$current")

        if parent != 0
            link_idx = link_dic[parent,current]
            if link_idx != 0
                x[link_idx] = 1
            end
        end

        current = parent
    end

    return x
end

# function get_shortest_path(start_node, end_node, link_length, origin, destination)
#     @assert Base.length(start_node)==Base.length(end_node)
#     @assert Base.length(start_node)==Base.length(link_length)
#
#     graph = create_graph(start_node, end_node)
#
#     state = dijkstra_shortest_paths(graph, link_length, origin)
#
#     path = get_path(state, origin, destination)
#     x = get_vector(path, start_node, end_node)
#
#     return path, x
# end
#
# function get_vector(state, origin, destination, start_node, end_node)
#     current = destination
#     parent = -1
#     x = zeros(Int, Base.length(start_node))
#
#     while parent != origin
#         parent = state.parents[current]
#
#         for j=1:Base.length(start_node)
#             if start_node[j]==parent && end_node[j]==current
#                 x[j] = 1
#                 break
#             end
#         end
#
#         current = parent
#     end
#
#     return x
# end


text:
module ModV5Packet
export V5Packet

using UUIDs

const MAGIC = 0x2b2b3556
const REPLY_MASK = 0x1

mutable struct V5Packet
    magic::UInt32
    requestid::UUID
    flags::UInt8
    length::UInt16
    payload::Vector{UInt8}

    function V5Packet(requestid, flags, payload)
        len = length(payload)
        len > typemax(UInt16) && throw(ArgumentError("Payload too large: $len"))
        new(MAGIC, requestid, flags, len, payload)
    end

end # struct

function V5Packet(payload::Vector{UInt8})::V5Packet
    V5Packet(uuid4(), 0, payload)
end

function V5Packet(payload::Vector{UInt8}, replyto::UUID)::V5Packet
    p = V5Packet(replyto, 0, payload)
    p.REPLY = true
    return p
end

getflag(pkt, mask) = (pkt.flags & mask) != 0
setflag(pkt, mask, x) = x ? pkt.flags |= mask : pkt.flags &= ~mask

function Base.getproperty(pkt::V5Packet, sym::Symbol)
    if sym == :REPLY
        getflag(pkt, REPLY_MASK)
    else
        getfield(pkt, sym)
    end
end

function Base.setproperty!(pkt::V5Packet, sym::Symbol, x)
    if sym == :REPLY
        setflag(pkt, REPLY_MASK, x)
    else
        setfield!(pkt, sym, x)
    end
end

function readguid(io::IO)
    a = read(io, 16)
    reverse!(a, 1, 4)
    reverse!(a, 5, 6)
    reverse!(a, 7, 8)
    reinterpret(UInt128, a)[1] |> ntoh |> UUID
end

function writeguid(io::IO, x::UUID)
    a = reinterpret(UInt8, [hton(x.value)])
    reverse!(a, 1, 4)
    reverse!(a, 5, 6)
    reverse!(a, 7, 8)
    write(io, a)
end

function Base.read(io::IO, ::Type{V5Packet})::V5Packet
    magic = read(io, UInt32)
    magic == MAGIC || throw(ErrorException("Invalid magic $magic"))
    requestid = readguid(io)
    flags = read(io, UInt8)
    len = read(io, UInt16)
    payload = read(io, len)
    V5Packet(requestid, flags, payload)
end

function Base.write(io::IO, x::V5Packet)
    nb = write(io, x.magic)
    nb += writeguid(io, x.requestid)
    nb += write(io, x.flags)
    nb += write(io, x.length)
    nb += write(io, x.payload)
end

end # module


text:
export CrossEntropy, start!, step!, release!

immutable CrossEntropy{C, A} <: PopulationBasedOptAlgo{C}
    ctx::C
    populationSize::Int
    dnaSize::Int
    mutationChance::Float32
    mutationStrength::Float32
    selectionStdDev::Float32
    keepElitesRate::Float32
    comparer::AbstractComparer
    popMan::PopulationManager
    offsprings::A
end

CrossEntropy(
    ctx,
    dnaSize,
    comparer::AbstractComparer,
    decode::Function;
    populationSize = 100,
    mutationChance = 0.05f0,
    mutationStrength = 0.02f0,
    selectionStdDev = 0.3f0,
    keepElitesRate = 0.05f0) =
CrossEntropy(
    ctx,
    populationSize,
    dnaSize,
    mutationChance,
    mutationStrength,
    selectionStdDev,
    keepElitesRate,
    comparer,
    PopulationManager(
        ctx,
        populationSize,
        dnaSize,
        comparer,
        decode),
    array(ctx, Float32, dnaSize, populationSize))

release!(ce::CrossEntropy) = (release!(ce.popMan); release!(ce.offsprings))

function start!(ce::CrossEntropy)
    reset!(ce.comparer)
    ce.offsprings[:,:] = 0f0
    randomize!(ce.popMan)
end

function step!(ce::CrossEntropy)
    reset!(ce.comparer)

    for i in 1:ce.populationSize
        parentIndex = chooseParentIndex(ce.popMan, ce.selectionStdDev);
        parent = ce.popMan.population[parentIndex]
        ce.offsprings[col(ce.ctx, i - 1)] = parent.dna.array;
    end

    candidate = createCandidate(ce.popMan)
    keepElites!(ce.popMan, candidate, ce.keepElitesRate)

    _mean = mean(ce.offsprings, 1)
    stdDev = stdev(ce.offsprings, 1)
    try
        while length(candidate) < ce.populationSize
            push!(candidate, scope!(ce.ctx) do this
                rnd = randn(ce.ctx, Float32, ce.dnaSize)
                dnaArray = rnd .* stdDev .+ _mean
                childDNA = DNA(ce.ctx, dnaArray)
                this.result(childDNA.array)
                mutate!(childDNA, ce.mutationChance, ce.mutationStrength)
                eval!(childDNA)
            end)
        end
    finally
        release!(_mean)
        release!(stdDev)
    end

    set!(ce.popMan, candidate)
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
#=
Created on Monday 18th of January 2021
Last update: Tuesday 19th of January 2021

@author: Bram De Jaegher
bram.de.jaegher@gmail.com
=#

#=______________________
|                       |
|         Day 1         |
|_______________________|
=#

#= Notebook 1: basics =#
### Clipping
function clip(x)
  if x ≤ 0 && return 0.
  elseif x ≥ 1 && return 1.
  else
  return x
  end
end

### Stirling
stirling(n) = sqrt(2*π*n)*(n/exp(1))^n

### Time is relative
function since_epoch(t)
  now = t
  seconds = now % 60 # 60 seconds in a minute
  now -= seconds
  # to minutes
  now = div(now, 60)
  minutes = now % 60 # 60 minutes in an hour
  now -= minutes
  # to hours
  now = div(now, 60)
  hours = now % 24 # 24 hours in a day
  now -= hours
  # to days
  now = div(now, 24)
  days = now % 365 # ± 365 days in a year
  return days, hours, minutes, seconds
end

### Fermat
function checkfermat(a::Int, b::Int, c::Int, n::Int)
  if a^n + b^n == c^n
      println("Holy smokes, Fermat was wrong!")
  else
      println("No, that doesn’t work.")
  end
end

### Justify
rightjustify(s) = println(" "^(70-length(s)) * s)

### Grid print
function printgrid()
  dashed = " -"^4*" "
  plusdash = "+"*dashed*"+"*dashed*"+"
  verticals = ("|"*" "^9)^2*"|"
  s = (plusdash * "\n" * (verticals * "\n")^4)^2*plusdash * "\n"
  println(s)
end

### Molecular mass spectrometry
f1(s) = replace(s, "-CH3"=> "-COOH") |> s -> replace(s, "CH3-"=> "COOH-")
f2(s) = replace(s, "-CO-COOH"=> "") |> s -> replace(s, "COOH-CO-"=> "")

function f3(molecule)
  if length(molecule) < 4
    molecule = "CH4"
    return molecule
  end

  if molecule[1:4] == "CH2-"
    molecule = "CH3-" * molecule[5:end]
  end
  if molecule[end-3:end] == "-CH2"
    molecule = molecule[1:end-4] * "-CH3"
  end
  return molecule
end

f4(s) = replace(s, "-CH2-CH2-"=> "-CO-")

function f5(molecule)
  if molecule[1:3] == "CO-"
    molecule = "COOH-" * molecule[4:end]
  end

  if molecule[end-2:end] == "-CO"
    molecule = molecule[1:end-3] * "-COOH"
  end
  return molecule
end

f6(s) = replace(s, "CO-CO-CO" => "CO-CH2-CO")

function cycle(molecule)
  molecule |> f1 |> f2 |> f3 |> f4 |> f5 |> f6
end

function degrade(molecule)
  mol_prev = ""
  while mol_prev !== molecule
    mol_prev = molecule
    molecule = cycle(molecule)
  end
  return molecule
end

function molecular_weight(molecule)
	return count("CH3", molecule)* 15.0 +
			count("CH4", molecule) * 16.0 +
			count("CH2", molecule) * 14.0 +
			count("CO", molecule) * 28.0 +
			count("OH", molecule) * 17
end

function simulate_path(space, offset, down=1)
  x, y = 1, 1
  n, m = size(space)
  n_trees = 0
  for y in 1:down:n
    n_trees += space[y,x]
    x += offset
    while x > m
      x -= m
    end
  end
  return n_trees
end

# n_trees = simulate_path(map, 3)
n_trees = 254
# product_all_trees = simulate_path(space, 1) * simulate_path(space, 3) * simulate_path(space, 5) * simulate_path(space, 7) * simulate_path(space, 1, 2)
product_all_trees = 1666768320
# answers to advanced solution:
#goingdowntheslope(3, 1, slope)
#slopestyles = (
#    (1, 1),
#    (3, 1),
#    (5, 1),
#    (7, 1),
#    (1, 2),
#)
#prod([goingdowntheslope(direction..., slope) for direction in slopestyles])

# WIP bigprint

#= Notebook 2: collections =#
### Riemann sum

function riemannsum(f, a, b; n=100)
  dx = (b - a) / n
  return sum(f.(a:dx:b)) * dx
end

### Markdown table
function markdowntable(table, header)
  n, m = size(table)
  table_string = ""
  # add header
  table_string *= "| " * join(header, " | ") * " |\n"
  # horizontal lines for separating the columns
  table_string *= "| " * repeat(":--|", m) * "\n"
  for i in 1:n
      table_string *= "| " * join(table[i,:], " | ") * " |\n"
  end
  #clipboard(table_string)
  return table_string
end

### Determinant
# WIP

### Pi Time
function estimatepi(n)
  hits = 0
  for i in 1:n
      x = rand()
      y = rand()
      if x^2 + y^2 ≤ 1.0
          hits += 1
      end
  end
  return 4hits/n
end

estimatepi2(n) = 4count(x-> x ≤ 1.0, sum(rand(n, 2).^2, dims=2)) / n

### Vandermonde
vandermonde(α, n) = [αᵢ^j for αᵢ in α, j in 0:n-1]

#= Notebook 3: project 1 =#
### Mean
mean(x) = sum(x)/length(x)

weighted_mean(x, w) = sum(w.*x)

### 1D convolution
function convolve_1d(x::Vector, w::Vector)
	@assert length(w) % 2 == 1 "length of `w` has to be odd!"
	@assert length(w) < length(x) "length of `w` should be smaller than `x`"
	n = length(x)
	m = length(w) ÷ 2
	y = zeros(n)

	fill!(y, 0.0)
	for (i, xj) in enumerate(x)
		for (j, wj) in enumerate(w)
			k = clamp(i + j - m, 1, n)
			y[i] += w[j] * x[k]
		end
	end
	return y
end

### Weight vectors
function uniform_weights(m)
	return ones(2m+1) / (2m+1)
end

function triangle_weights(m)
	w = zeros(2m+1)
	for i in 1:(m+1)
		w[i] = i
		w[end-i+1] = i
	end
	w ./= sum(w)
	return w
end

function gaussian_weights(m; σ=4)
	w = exp.(-(-m:m).^2 / 2σ^2)
	return w ./ sum(w)
end

### Protein sliding window
function protein_sliding_window(sequence, m, zscales)
	n = length(sequence)
	y = zeros(n)
	for i in (m+1):(n-m)
		for k in -m:m
			y[i] += zscales[sequence[i+k]]
		end
	end
	return y / (2m + 1)
end

### Decimate
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

### 2D convolution
function convolve_2d(M::Matrix, K::Matrix)
	out = similar(M)
	n_rows, n_cols = size(M)
	fill!(out, 0.0)
	m = div(size(K, 1), 2)
	for i in 1:n_rows
		for j in 1:n_cols
			for k in -m:m
				for l in -m:m
					out[i,j] += M[clamp(i+k, 1, n_rows), clamp(j+l, 1, n_cols)] * K[k+m+1,l+m+1]
				end
			end
		end
	end
	return out
end

function gaussian_kernel(m; σ=4)
  K = [exp(-(x^2 + y^2) / 2σ^2) for x in -m:m, y in -m:m]
  K ./= sum(K)
  return K
end

function convolve_image(M::Matrix{<:AbstractRGB}, K::Matrix)
	Mred = convolve_2d(red.(M) .|> Float32, K)
	Mgreen = convolve_2d(green.(M) .|> Float32, K)
	Mblue = convolve_2d(blue.(M) .|> Float32, K)
	return RGB.(Mred, Mgreen, Mblue)
end

### Edge detection
Gx = [1 0 -1; 2 0 -2; 1 0 -1]
Gy = [1 2 1; 0 0 0; -1 -2 -1]
function edge_detection(M)
	M = M .|> Gray .|> Float64
	return sqrt.(Solutions.convolve_2d(M, Gx).^2 + Solutions.convolve_2d(M, Gy).^2) .|> Gray
end

### Get index of bitstring
getbinarydigit(rule, i) = isodd(rule >> i)

### Next state
nextstate(l::Bool, s::Bool, r::Bool, rule::Int) = nextstate(l, s, r, UInt8(rule))

function nextstate(l::Bool, s::Bool, r::Bool, rule::UInt8)
  return getbinarydigit(rule, 4l+2s+1r)
end

### Update 1 step of CA
function update1dca!(xnew, x, rule::Integer)
	n = length(x)
	xnew[1] = nextstate(x[end], x[1], x[2], rule)
	xnew[end] = nextstate(x[end-1], x[end], x[1], rule)
	for i in 2:n-1
		xnew[i] = nextstate(x[i-1], x[i], x[i+1], rule)
	end
	return xnew
end

update1dca(x, rule::Integer) = update1dca!(similar(x), x, rule)

### simulating the CA
"""
    simulate(x0, rule; nsteps=100)

Simulate `nsteps` time steps according to `rule` with `X0` as the initial condition.
Returns a matrix X, where the rows are the state vectors at different time steps.
"""
function simulate(x0, rule::UInt8; nsteps=100)
	n = length(x0)
    X = zeros(Bool, nsteps+1, n)
	X[1,:] = x0
	for t in 1:nsteps
		x = @view X[t,:]
		xnew = @view X[t+1,:]
		update1dca!(xnew, x, rule)
	end
    return X
end

### plotting CA
#plot(ca_image(X), size=(1000, 1000))

#=______________________
|                       |
|         Day 2         |
|_______________________|
=#

#= Notebook 1: types =#
### String parsing
bunchofnumbers_parser(bunchofnumbers) = parse.(Float64, split(rstrip(bunchofnumbers), ", ")) |> sum

#= Notebook 2: composite types =#
### Wizarding currency
struct WizCur
  galleons::Int
  sickles::Int
  knuts::Int
  function WizCur(galleons::Int, sickles::Int, knuts::Int)
      sickles += knuts ÷ 29
      knuts %= 29
      galleons += sickles ÷ 17
      sickles %= 17
      return new(galleons, sickles, knuts)
  end
end

galleons(money::WizCur) = money.galleons
sickles(money::WizCur) = money.sickles
knuts(money::WizCur) = money.knuts

moneyinknuts(money::WizCur) = 29*17galleons(money) + 29sickles(money) + knuts(money)

function Base.show(io::IO, money::WizCur)
  print(io, "$(galleons(money))G, $(sickles(money))S, $(knuts(money))K")
end

Base.isless(m1::WizCur, m2::WizCur) = moneyinknuts(m1) < moneyinknuts(m2)
Base.isgreater(m1::WizCur, m2::WizCur) = moneyinknuts(m1) > moneyinknuts(m2)
Base.isequal(m1::WizCur, m2::WizCur) = moneyinknuts(m1) == moneyinknuts(m2)

Base.:+(m1::WizCur, m2::WizCur) = WizCur(galleons(m1)+galleons(m2),
                                          sickles(m1)+sickles(m2),
                                          knuts(m1)+knuts(m2))

money_ron = WizCur(0, 19, 732)
money_harry = WizCur(3, 1, 7)

dungbomb_fund = money_ron + money_harry

### Vandermonde
struct Vandermonde{T,VT} <: AbstractMatrix{T}
  α::VT
  m::Int
  Vandermonde(α::AbstractVector{T}, m) where {T} = new{T,typeof(α)}(α,m)
end

Vandermonde(α::Vector{<:Number}) = Vandermonde(α, length(α))

Base.size(V::Vandermonde) = (length(V.α), V.m)
Base.getindex(V::Vandermonde, i, j) = V.α[i]^(j-1)

α = [1, 2, 3, 4]

V = Vandermonde(α, 4)

determinant(V::Vandermonde) =
  ((xi-xj) for (i,xi) in enumerate(V.α), (j, xj) in enumerate(V.α) if i < j) |> prod


text:
const Link_Collection_Type = Vector{String}

# +Links
#
# This function cleans extracted links by removing duplicates
#
# See: [[extract_links_string][]]
function remove_link_duplicate(toClean::Link_Collection_Type)::Link_Collection_Type
    # a priori sort is not necessary
    return unique(toClean)
end

# +Links L:extract_links_string
#
# This function returns the list of links found in the string. It
# returns nothing if no link is found
#
# !s="# some text [[some_target][]], [[link_target][link_name]]";
# !J4Org.extract_links(s)
#
# *Caveat:* only use links of the forme "[[something][]]", which are not
#           valid Org mode links, see [[doc_link_substitution][]]
#
# *Test link:* [[doc_link_substitution][]]
#
# *Note*: duplicates are removed
#
function extract_links(input::String)::Link_Collection_Type
    # remove blalba "toremove" kmsdqkm (do not consider quoted links)
    input = replace(input,r"(\".*\")" => "")
    v=Link_Collection_Type()

    match_link = r"\[\[(\w+)\]\[\]\]"

    if !occursin(match_link,input)
        return v
    end

    offset=1
    while (m=match(match_link,input,offset))!=nothing
        push!(v,m[1])
        offset=m.offsets[end]+length(m[1])
    end

    return remove_link_duplicate(v)
end

# +Links
#
# This function is like [[extract_links_string][]], except that it
# processes a [[Documented_Item][]]
#
# *Note*: duplicates are removed
#
function extract_links(di::Documented_Item)::Link_Collection_Type
    v=extract_links(raw_string_doc(di))
    return remove_link_duplicate(v)
end

# +Links
#
# This function is like [[extract_links_string][]], except that it
# processes an array of [[Documented_Item][]]
#
# *Note*: duplicates are removed
#
function extract_links(di_array::Array{Documented_Item,1})::Link_Collection_Type
    v=Link_Collection_Type()

    for di in di_array
        v=vcat(v,extract_links(di))
    end

    return remove_link_duplicate(v)
end

# +Links
#
# Returns the indices of [[Documented_Item][]] containing the
# link_target (have a tag line with L:link_target)
#
# *Note:* a normal situation is to have zero or one indices. Several
# indices means that we do not have a unique target.
#
function get_item_idx_from_link_target(link_target::String,di_array::Array{Documented_Item,1})::Vector{Int}
    v=Vector{Int}()
    n = length(di_array)
    for i in 1:n
        if link_target==link(di_array[i])
            push!(v,i)
        end
    end
    return v
end

# +Links
# Check for if links are well formed
#
# Two kinds of errors:
# - links without target
# - duplicate link targets
#
# *Returns*: true if ok, false if some errors are detected
#
function check_for_link_error(di_array::Array{Documented_Item,1},
                              di_array_universe::Array{Documented_Item,1})::Bool

    visited_links = Link_Collection_Type() # store visited links -> do not print error several times

    ok=true

    for di in di_array
        links_to_check = extract_links(di)

        for link_target in links_to_check

            if link_target ∉ visited_links

                push!(visited_links,link_target)

                item_idx = get_item_idx_from_link_target(link_target,di_array_universe)

                if isempty(item_idx)
                    ok=false
                    warning_message("link target $(link_target) in $(create_file_org_link(di)) not found")
                elseif length(item_idx)>1
                    ok=false
                    for idx in item_idx
                        warning_message("duplicate link target $(duplicated_link) presents in $(create_file_org_link(di_array_universe[idx]))")
                    end
                end
            end
        end
    end

    return ok
end

# +Links
#
# Creates readable part of the link.
#
# By default use item identifier.
#
# Improve it if we can, for instance
#
# - for functions: identifier -> identifier(...)
# - for enums: identifier -> @enum identifier
# - ...
#
# *Parameters*:
# - di: the item containing the link target (L:link_target)
#
# *Post-condition*:
# - returns a non-empty string
#
function create_magnified_link_name(di::Documented_Item)::String

    # start with the identifier
    identifier_as_string = identifier(di)

    # try to magnify
    # - function add identifier()
    # - strucure add struct identifier
    # - ...
    if is_documented_function(di)
        magnified_link_name=identifier_as_string*"(...)"
    elseif is_documented_structure(di)
        magnified_link_name="struct "*identifier_as_string
    elseif is_documented_abstract_type(di)
        magnified_link_name="abstract "*identifier_as_string
    elseif is_documented_enum_type(di)
        magnified_link_name="@enum "*identifier_as_string
    elseif is_documented_macro_type(di)
        magnified_link_name="@"*identifier_as_string
    elseif is_documented_variable_type(di)
        magnified_link_name="variable "*identifier_as_string
    end

    @assert !isempty(magnified_link_name)

    return magnified_link_name
end

# +Links
#
# - new_target=="": transforms [[link_target][]] -> _link_magnified_
# - otherwise: transforms [[link_target][]] -> [[link_new_target][link_magnified]]
#
function doc_link_substitution_helper(doc::String,
                                     link_target,
                                     link_new_target,
                                     link_magnified)::String
    src="[[$(link_target)][]]"
    if isempty(link_new_target)
        dest="_$(link_magnified)_"
    else
        dest="[[$(link_new_target)][$(link_magnified)]]"
    end
    return replace(doc,src => dest)
end

# +Links, TODO                                       L:doc_link_substitution
# From doc string performs links substitution
#
# *Note:* in order to do not interfere with org mode link we only process "links" of the form "[[something][]]"
#         see https://orgmode.org/manual/Link-format.html
#
# *Note:* to be able to write a "unactive" link, use C-x 8 RET 200b
#         (see: https://emacs.stackexchange.com/a/16702)
#
# - [ ] TODO Finally duplicate links are not checked... do it before
#
function doc_link_substitution(doc::String,
                               di_array::Array{Documented_Item,1},
                               di_array_universe::Array{Documented_Item,1},
                               link_prefix::String)::String
    # Extract links
    links = extract_links(doc)

    # Something to?
    if isempty(links)
        return doc
    end

    # Process each link, sequentially modify doc
    for link_target in links
        # default values (to be modified)
        # link_new_target = "" if target not found
        link_new_target = link_prefix*link_target
        # use create_magnified_link_name()
        link_magnified = link_target

        # Find item indices having link_target as target (L:link_target)
        link_target_idx = get_item_idx_from_link_target(link_target,di_array)

        # No target found
        if isempty(link_target_idx)
            link_new_target = "" # inactive link

            # Try to magnify name from di_array_universe
            if !(di_array===di_array_universe)
                # Try to find one in di_array_universe
                link_target_idx = get_item_idx_from_link_target(link_target,di_array_universe)

                if !isempty(link_target_idx)
                    link_magnified = create_magnified_link_name(di_array_universe[link_target_idx[1]])
                end
            end
        else
            link_magnified= create_magnified_link_name(di_array[link_target_idx[1]])
        end

        doc = doc_link_substitution_helper(doc,
                                           link_target,
                                           link_new_target,
                                           link_magnified)
    end

    return doc
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Ball(center, radius)

A Ball with `center` and `radius`.
"""
struct Ball{Dim,T} <: Primitive{Dim,T}
  center::Point{Dim,T}
  radius::T
end

Ball(center::Tuple, radius) = Ball(Point(center), radius)

paramdim(::Type{<:Ball{Dim}}) where {Dim} = Dim

center(b::Ball) = b.center
radius(b::Ball) = b.radius

# https://en.wikipedia.org/wiki/Volume_of_an_n-ball
function measure(b::Ball{Dim}) where {Dim}
  r, n = b.radius, Dim
  (π^(n/2) * r^n) / gamma(n/2 + 1)
end

function Base.in(p::Point, b::Ball)
  x = coordinates(p)
  c = coordinates(b.center)
  r = b.radius
  sum(abs2, x - c) ≤ r^2
end

boundary(b::Ball) = Sphere(b.center, b.radius)


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>2.0, :alg=>TMJets(maxsteps=1_000, abstol=1e-10, orderT=6, orderQ=2)),
    "X0"   => Hyperrectangle(low=fill(0.99, 7), high=fill(1.01, 7)),
    "info" => Dict(
        "name" => "biomodel7d",
        "dim" => 7,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
"""
    MarkovChains

An optimized module for efficient Markov chain simulation.
"""
module MarkovChains
using Distributions: AbstractRNG, Categorical, Discrete, Multivariate, rand, Sampleable
using LightGraphs: AbstractGraph
using StochasticProcesses
using StochasticProcesses.Samplers

export

# abstract types
AbstractMarkovChain, AbstractStochasticProcess,

# concrete types
MarkovChain,

# methods
rand

# include scripts
include("markovchains.jl")
include("simulatemarkovchains.jl")
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
using Flux: loadparams!

import ..loading_method

loading_method(::Val{:perfrnn}) = load_perfrnn

function load_perfrnn(weights, config)
    lstm_units = config["lstm_units"]
    input_dims = config["input_dims"]

    layers = []
    for i in 1:config["num_layers"]
        layer = BasicLSTM(
                     (i == 1) ? input_dims : lstm_units,
                     lstm_units
                )
        push!(layers, layer)
    end
    push!(layers, Dense(lstm_units, input_dims))

    model = Chain(layers...)

    weight_names = keys(weights)
    rnn_weights = filter(name -> occursin("rnn", name), weight_names)
    dense_weights = filter(name -> occursin("fully_connected", name), weight_names)

    for i in 1:length(rnn_weights)
        cell_weights = filter(name -> occursin("cell_$(i-1)", name), rnn_weights)
        for j in cell_weights
            if occursin("kernel", j)
                kernel = weights[j]
                if i == 1
                    Wi = @view kernel[:, 1:input_dims]
                    Wh = @view kernel[:, input_dims+1:end]
                else
                    Wi = @view kernel[:, 1:lstm_units]
                    Wh = @view kernel[:, lstm_units+1:end]
                end
                loadparams!(model[i], [Wi, Wh])
            elseif occursin("bias", j)
                loadparams!(model[i].cell.b, [weights[j]])
            end
        end
    end

    for j in dense_weights
        if occursin("weights", j)
            loadparams!(model[end].W, [weights[j]])
        elseif occursin("biases", j)
            loadparams!(model[end].b, [weights[j]])
        end
    end

    num_velocitybins = config["velocity_bins"]
    perfencoder = PerformanceOneHotEncoding(num_velocitybins=num_velocitybins)
    perfrnn = PerfRNN(model, perfencoder, num_velocitybins)
end


text:
#These test/demonstrate the basic funcionality

@diff_eq F begin
    dy[1] = -y[1]
end
@diff_eq_jac J begin
    pd[1,1] = -1.0
end

#From the book Numerical recipes
@diff_eq F2 begin
    dy[1] = 998*y[1] + 1998*y[2]
    dy[2] = -999*y[1] - 1999*y[2]
end
@diff_eq_jac J2 begin
    pd[1,1] = 998
    pd[1,2] = 1998
    pd[2,1] = -999
    pd[2,2] = -1999
end

@diff_eq F3 begin
    dy .= -y
end
@diff_eq_jac J3 begin
    pd .= -1.0*ones(1,1)
end

y0 = [1.0]
t = range(0,stop=1,length=10)

y02 = [1.0,0.0]

s_exp = exp.(-t)

s_exp2 = hcat( [2*exp.(-t) - exp.(-1000*t),
                -exp.(-t) + exp.(-1000*t)]... )

tol = 1.0e-2

s1 = hcat(ode(F,t,y0,verbose=false)...)'   #jacobian generated, stiff
@test sqrt(sum(abs2, s1 - s_exp)) < tol
s2 = hcat(ode(F,J,t,y0,verbose=false)...)'  #jacobian supplied
@test sqrt(sum(abs2, s2 - s_exp)) < tol
s3 = hcat(ode(F,t,y0,stiff=false,verbose=false)...)'  #jacobian ignored, not stiff
@test sqrt(sum(abs2, s3 - s_exp)) < tol
s4 = hcat(ode(F,J,t,y0,[1.0e-4],verbose=false)...)'  #increased absolute tolerance
@test sqrt(sum(abs2, s2 - s_exp)) > sqrt(sum(abs2, s4 - s_exp))
s5 = hcat(ode(F,J,t,y0,[1.0e-4],1.0e-4,verbose=false)...)'  #increased relative tolerance
@test sqrt(sum(abs2, s5 - s_exp)) < sqrt(sum(abs2, s4 - s_exp))

#Stiffs
s6 = hcat(ode(F2,t,y02,verbose=false)...)'   #jacobian generated, stiff
@test sqrt(sum(abs2, s6 - s_exp2)) < tol
s7 = hcat(ode(F2,J2,t,y02,verbose=false)...)'  #jacobian supplied
@test sqrt(sum(abs2, s7 - s_exp2)) < tol
s8 = hcat(ode(F2,t,y02,stiff=false,verbose=false)...)'  #jacobian ignored, not stiff
@test sqrt(sum(abs2, s8 - s_exp2)) < tol

#Vectorized
s9 = hcat(ode(F3,J3,t,y0,verbose=false)...)'  #jacobian supplied
@test all( s9 .≈ s2 )


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8, 6))
plot(y_hat, color="blue", label="Predicted")
plot(y_wind[test], color="red", label="Observed")
xlabel("Time", fontsize=14)
ylabel("Power generation", fontsize=14)
xticks([])
yticks(fontsize=12)
xlim(0, 100)
legend(fontsize=14)

savefig(joinpath(@OUTPUT, "obs_v_pred.svg")) # hide


text:
using Ipopt, StaticArrays, LinearAlgebra, BenchmarkTools, Plots
const TO = TrajectoryOptimization

# Car escape
T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-2,
    constraint_tolerance=max_con_viol,
    penalty_scaling=50.,
    penalty_initial=10.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

# AL-iLQR
prob = copy(Problems.car_3obs_static)
U0 = deepcopy(controls(prob))
al = AugmentedLagrangianSolver(prob, opts_al)
initial_controls!(al, U0)
solve!(al)
max_violation(al)

@btime begin
    initial_controls!($al, $U0)
    solve!($al)
end

p = plot(prob.model, get_trajectory(al))
plot!(prob.constraints[1].con)
display(p)

# ALTRO
altro = ALTROSolver(prob, opts_altro)
initial_controls!(altro, U0)
solve!(altro)
max_violation(altro)

@btime begin
    initial_controls!($altro, $U0)
    solve!($altro)
end
max_violation(altro)

# Ipopt
prob_ipopt = copy(Problems.car_3obs_static)
prob_ipopt = TO.change_integration(prob_ipopt, HermiteSimpson)
rollout!(prob_ipopt)
ipopt = DIRCOLSolver(prob_ipopt, opts_ipopt)
ipopt.opts.verbose = true
solve!(ipopt)
max_violation(ipopt)

ipopt.opts.verbose = false
@btime solve!($ipopt)

p = plot(prob.model, get_trajectory(ipopt))
plot!(prob.constraints[1].con)
display(p)
plot(controls(ipopt))


text:
# Dates, LeapSeconds
# TAI UTC
# GPS BEIDOU GLONASS

include("Constant.jl")

include("TimeSystem.jl")
include("TimeSystemCorr.jl")

include("Convert.jl")


text:
using BinaryBuilder

# Collection of sources required to build ZMQ
sources = [
    "https://github.com/zeromq/libzmq.git" =>
    "2cb1240db64ce1ea299e00474c646a2453a8435b", # v4.3.1
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libzmq
sh autogen.sh
./configure --prefix=$prefix --host=${target} \
    --without-docs --disable-libunwind --disable-perf --disable-Werror \
    --disable-eventfd --without-gcov --disable-curve-keygen \
    CXXFLAGS="-g -O2 -fms-extensions"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms() # build on all supported platforms

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libzmq", :libzmq),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "ZMQ", v"4.3.1", sources, script, platforms, products, dependencies)


text:
import Base.show

mutable struct Neo4jInstance
  connection::Neo4j.Connection
  graph::Neo4j.Graph
end

mutable struct CloudGraphsDFG{T <: AbstractParams} <: AbstractDFG{T}
    neo4jInstance::Neo4jInstance
    userId::String
    robotId::String
    sessionId::String
    description::String #TODO Maybe remove description
    addHistory::Vector{Symbol}
    solverParams::T # Solver parameters

    # inner constructor for all constructors in common
    function CloudGraphsDFG{T}(neo4jInstance::Neo4jInstance,
                               userId::String,
                               robotId::String,
                               sessionId::String,
                               description::String,
                               addHistory::Vector{Symbol},
                               solverParams::T;
                               createSessionNodes::Bool=true,
                               userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                               robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                               sessionData::Dict{Symbol, String} = Dict{Symbol, String}()) where T <: AbstractParams
        # Validate the userId, robotId, and sessionId
        !isValidLabel(userId) && error("'$userId' is not a valid User ID")
        !isValidLabel(robotId) && error("'$robotId' is not a valid Robot ID")
        !isValidLabel(sessionId) && error("'$sessionId' is not a valid Session ID")

        # neo4jConnection = Neo4j.Connection(host, port=port, user=dbUser, password=dbPassword);
        # graph = Neo4j.getgraph(neo4jConnection)
        # neo4jInstance = Neo4jInstance(neo4jConnection, graph)

        dfg = new{T}(neo4jInstance, userId, robotId, sessionId, description, addHistory, solverParams)
        # Create the session if it doesn't already exist
        if createSessionNodes
            createDfgSessionIfNotExist(dfg)
            setUserData!(dfg, userData)
            setRobotData!(dfg, robotData)
            setSessionData!(dfg, sessionData)
            setDescription!(dfg, description)
        end

        return dfg
    end
end

"""
    $(SIGNATURES)
Create a new CloudGraphs-based DFG factor graph using a Neo4j.Connection or by specifying the Neo4j connection information
"""
function CloudGraphsDFG{T}(neo4jConnection::Neo4j.Connection,
                           userId::String,
                           robotId::String,
                           sessionId::String,
                           description::String;
                           solverParams::T=NoSolverParams(),
                           kwargs...) where T <: AbstractParams

    graph = Neo4j.getgraph(neo4jConnection)
    neo4jInstance = Neo4jInstance(neo4jConnection, graph)

    return CloudGraphsDFG{T}(neo4jInstance, userId, robotId, sessionId, description, Symbol[], solverParams; kwargs...)

end

function CloudGraphsDFG{T}(host::String,
                           port::Int,
                           dbUser::String,
                           dbPassword::String,
                           userId::String,
                           robotId::String,
                           sessionId::String,
                           description::String;
                           kwargs...) where T <: AbstractParams
    neo4jConnection = Neo4j.Connection(host, port=port, user=dbUser, password=dbPassword)
    return CloudGraphsDFG{T}(neo4jConnection, userId, robotId, sessionId, description; kwargs...)
end

# construct using the default settings for localhost
function CloudGraphsDFG(; hostname="localhost",
                          port=7474,
                          username="neo4j",
                          password="test",
                          userId::String="DefaultUser",
                          robotId::String="DefaultRobot",
                          sessionId::String="Session_$(string(uuid4())[1:6])", #TODO randstring(['a':'z';'A':'Z'],1) ipv Session
                          description::String="CloudGraphsDFG implementation",
                          solverParams::T=NoSolverParams(),
                          kwargs...) where T <: AbstractParams

    @info "Creating $sessionId"
    return CloudGraphsDFG{T}(hostname,
                             port,
                             username,
                             password,
                             userId,
                             robotId,
                             sessionId,
                             description;
                             solverParams=solverParams,
                             kwargs...)
end

function CloudGraphsDFG(description::String,
                        userId::String,
                        robotId::String,
                        sessionId::String,
                        userData::Dict{Symbol, String},
                        robotData::Dict{Symbol, String},
                        sessionData::Dict{Symbol, String},
                        solverParams::AbstractParams;
                        host::String = "localhost",
                        port::Int = 7474,
                        dbUser::String = "neo4j",
                        dbPassword::String = "test")

    return CloudGraphsDFG{typeof(solverParams)}(host,
                                                port,
                                                dbUser,
                                                dbPassword,
                                                userId,
                                                robotId,
                                                sessionId,
                                                description;
                                                solverParams=solverParams,
                                                userData=userData,
                                                robotData=robotData,
                                                sessionData=sessionData)

end

function show(io::IO, ::MIME"text/plain", c::CloudGraphsDFG)
    println(io, "CloudGraphsDFG:")
    println(io, " - Neo4J instance: $(c.neo4jInstance.connection.host)")
    println(io, " - Session: $(c.userId):$(c.robotId):$(c.sessionId)")
    println(io, " - Description: ", c.description)
end


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_tso_out_fo_dp_limitable" begin

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     | S2: 60
    DP => 10h40         |
      S1: 55            |
      S2: 65            |
    =#

    TS = [DateTime("2015-01-01T11:00:00")]
    network = PSCOPF.Networks.Network()
    # Buses
    PSCOPF.Networks.add_new_bus!(network, "bus_1")
    # Limitables
    PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                            0., 100.,
                                            0., 1.,
                                            Dates.Second(20*60), Dates.Second(20*60))
    # Uncertainties
    uncertainties = PSCOPF.Uncertainties()
    # initial generators state : need to pay starting cost at TS[1]
    generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()
    #ManagementMode
    mode = PSCOPF.ManagementMode("test_mode", Dates.Minute(5))

    #=
    ECH1*                    ECH2              ECH3         <---FO--->TS
    |                        |
    10h                      10h40             10h45                 11h
                             <--------------------DP(wind1)----------->
    =#
    @testset "tso_out_fo_can_change_limit_by_scenario_before_dp" begin
        ech = DateTime("2015-01-01T10:00:00")
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 65.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

        # firmness
        firmness = PSCOPF.Firmness(
                SortedDict(),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                )

        context = PSCOPF.PSCOPFContext(network, TS, mode,
                                        generators_init_state,
                                        uncertainties, nothing)

        #This limit should be ignored because firmness is FREE
        OLD_LIMIT = 55.
        PSCOPF.set_limitation_value!(context.tso_actions, "wind_1_1", TS[1], OLD_LIMIT)

        context.tso_schedule = PSCOPF.Schedule(PSCOPF.TSO(), Dates.DateTime("2015-01-01T06:00:00"), SortedDict(
                                            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                                SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>55.,"S2"=>55.)))
                                                ),
                                            )
                                    )

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T10:00:00"), SortedDict(
                                        "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                            SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                            SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>57.,"S2"=>62.)))
                                            ),
                                        )
                                )

        tso = PSCOPF.TSOOutFO()

        @test firmness == PSCOPF.compute_firmness(tso,
                                                ech, DateTime("2015-01-01T10:40:00"),
                                                TS, context)

        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # Limit is by scenario before DP
        @test 1e-09 < abs( value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"])
                            - value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) )

        # Limit was changed by scenario since we are before DP
        @test 55. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]) # ≈ 55. == OLD_LIMIT #changed
        @test 60. ≈ value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) # > 55. == OLD_LIMIT #changed

        #"S1" this is NOT an active limit
        @test value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        @test value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        #"S2" this is an active limit
        @test 1. ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"])
        @test 60. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"])

        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"]) < 1e-09
    end

    #=
    ECH1                    ECH2              ECH3*         <---FO--->TS
    |                        |
    10h                      10h40             10h45                 11h
                             <--------------------DP(wind1)----------->
    =#
    @testset "tso_out_fo_can_increase_the_common_limit_after_dp" begin
        ech = DateTime("2015-01-01T10:45:00")
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 65.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

        # firmness
        firmness = PSCOPF.Firmness(
                SortedDict(),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                )

        context = PSCOPF.PSCOPFContext(network, TS, mode,
                                        generators_init_state,
                                        uncertainties, nothing)

        OLD_LIMIT = 55.
        PSCOPF.set_limitation_value!(context.tso_actions, "wind_1_1", TS[1], OLD_LIMIT)

        context.tso_schedule = PSCOPF.Schedule(PSCOPF.TSO(), Dates.DateTime("2015-01-01T10:40:00"), SortedDict(
                                            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                                SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>55.,"S2"=>55.)))
                                                ),
                                            )
                                    )

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T10:45:00"), SortedDict(
                                        "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                            SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                            SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>57.,"S2"=>62.)))
                                            ),
                                        )
                                )

        tso = PSCOPF.TSOOutFO()

        @test firmness == PSCOPF.compute_firmness(tso,
                                                ech, DateTime("2015-01-01T10:50:00"),
                                                TS, context)

        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution has slacks
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # Limit is firm after DP : ie common to all scenarios
        @test ≈(value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]),
                value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) )

        # Limit can be changed after DP :
        @test 55. == OLD_LIMIT < value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]) ≈ 60
        @test 55. == OLD_LIMIT < value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) ≈ 60

        @test value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        @test value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        @test 55. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1], "S1"])
        @test 1. ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"])
        @test 60. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"])
        @test 60. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1], "S2"])
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"]) < 1e-09
    end

    #=
    ECH1                    ECH2              ECH3*         <---FO--->TS
    |                        |
    10h                      10h40             10h45                 11h
                             <--------------------DP(wind1)----------->
    =#
    @testset "tso_out_fo_can_decrease_the_common_limit_after_dp" begin
        ech = DateTime("2015-01-01T10:45:00")
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 65.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)

        # firmness
        firmness = PSCOPF.Firmness(
                SortedDict(),
                SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.DECIDED), )
                )

        context = PSCOPF.PSCOPFContext(network, TS, mode,
                                        generators_init_state,
                                        uncertainties, nothing)

        OLD_LIMIT = 70.
        PSCOPF.set_limitation_value!(context.tso_actions, "wind_1_1", TS[1], OLD_LIMIT)

        context.tso_schedule = PSCOPF.Schedule(PSCOPF.TSO(), Dates.DateTime("2015-01-01T10:40:00"), SortedDict(
                                            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                                SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                                SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>70.,"S2"=>70.)))
                                                ),
                                            )
                                    )

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T10:45:00"), SortedDict(
                                        "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                                            SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                                            SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>57.,"S2"=>62.)))
                                            ),
                                        )
                                )

        tso = PSCOPF.TSOOutFO()

        @test firmness == PSCOPF.compute_firmness(tso,
                                                ech, DateTime("2015-01-01T10:50:00"),
                                                TS, context)

        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution has slacks
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # Limit is firm after DP : ie common to all scenarios
        @test ≈(value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]),
                value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) )

        # Limit can be changed after DP :
        @test 70. == OLD_LIMIT > value(result.limitable_model.p_limit["wind_1_1",TS[1], "S1"]) ≈ 60
        @test 70. == OLD_LIMIT > value(result.limitable_model.p_limit["wind_1_1",TS[1], "S2"]) ≈ 60

        @test value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        @test value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S1"]) < 1e-09
        @test 55. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1], "S1"])
        @test 1. ≈ value(result.limitable_model.b_is_limited["wind_1_1",TS[1], "S2"])
        @test 60. ≈ value(result.limitable_model.p_limit_x_is_limited["wind_1_1",TS[1], "S2"])
        @test 60. ≈ value(result.limitable_model.p_injected["wind_1_1",TS[1], "S2"])
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S2"]) < 1e-09
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# helper functions

"""
    r_lattice(prob)

Return a 2d array of size (3, N_r) where the columns are the lattice coordinates of a uniform r grid in one unit cell with size_r(prob)[i] points in direction i.
"""
function r_lattice(prob::AbstractBSEProblem)
    N_rs = size_r(prob)
    rl_1d_1 = range(0.0, stop = 1.0, length = N_rs[1] + 1)[1:N_rs[1]]
    rl_1d_2 = range(0.0, stop = 1.0, length = N_rs[2] + 1)[1:N_rs[2]]
    rl_1d_3 = range(0.0, stop = 1.0, length = N_rs[3] + 1)[1:N_rs[3]]
    rl_1 = kron(ones(N_rs[3]), ones(N_rs[2]), rl_1d_1)
    rl_2 = kron(ones(N_rs[3]), rl_1d_2, ones(N_rs[1]))
    rl_3 = kron(rl_1d_3, ones(N_rs[2]), ones(N_rs[1]))
    rl = transpose(hcat(rl_1, rl_2, rl_3))
    return rl
end

"""
    r_cartesian(prob)

Return a 2d array of size (3, N_r) where the columns are the cartesian coordinates of the r grid generated by r_lattice.
"""
function r_cartesian(prob::AbstractBSEProblem)
    a_mat = lattice_matrix(prob)
    rl = r_lattice(prob)
    rc = similar(rl)
    for ir in 1:size(rl, 2)
        rc[:, ir] = a_mat * rl[:, ir]
    end
    return rc
end

"""
    k_lattice(prob[, shift])

Return a 2d array of size `(3, N_k)` where the columns are the lattice coordinates of a uniform ``k`` grid in one reciprocal unit cell with `size_k(prob)[i]` points in direction `i`. The grid can optionally be shifted by a vector `shift`.
"""
function k_lattice(prob::AbstractBSEProblem, shift = [0.0, 0.0, 0.0])
    N_ks = size_k(prob)
    kl_1d_1 = mod.(range(0.0, stop = 1.0, length = N_ks[1] + 1)[1:N_ks[1]] .+ shift[1], 1.0)
    kl_1d_2 = mod.(range(0.0, stop = 1.0, length = N_ks[2] + 1)[1:N_ks[2]] .+ shift[2], 1.0)
    kl_1d_3 = mod.(range(0.0, stop = 1.0, length = N_ks[3] + 1)[1:N_ks[3]] .+ shift[3], 1.0)
    kl_1 = kron(ones(N_ks[3]), ones(N_ks[2]), kl_1d_1)
    kl_2 = kron(ones(N_ks[3]), kl_1d_2, ones(N_ks[1]))
    kl_3 = kron(kl_1d_3, ones(N_ks[2]), ones(N_ks[1]))
    kl = transpose(hcat(kl_1, kl_2, kl_3))
    return kl
end

"""
    k_cartesian(prob[, shift])

Return a 2d array of size `(3, N_k)` where the columns are the
lattice coordinates of a uniform ``k`` grid in the first Brillouin zone
with `size_k(prob)[i]` points in direction `i`. The grid can
optionally be shifted by a vector `shift`.
The cartesian representation is chosen to lie in the first Brillouin
zone, i.e. to have minimal distance from the origin.
"""
function k_cartesian(prob::AbstractBSEProblem, shift = [0.0, 0.0, 0.0])
    a_mat = lattice_matrix(prob)
    b_mat = 2 * π * inv(a_mat)
    kl = k_lattice(prob, b_mat \ shift)
    kc = similar(kl)
    for ik in 1:size(kl, 2)
        kc[:, ik] = k_lattice_to_reduced_bz(kl[:, ik], b_mat)
    end
    return kc
end

"""
    k_lattice_to_reduced_bz(k_l, b_mat)

Helper function to compute cartesian coordinate of reciprocal point given in lattice coordinates.
"""
function k_lattice_to_reduced_bz(k_l, b_mat)
    shifted_norms = zeros(2, 2, 2)
    for i1 in 1:2, i2 in 1:2, i3 in 1:2
        shifted_norms[i1, i2, i3] = norm(b_mat * (k_l - [i1 - 1, i2 - 1, i3 - 1]))
    end
    m, i = findmin(shifted_norms)
    return b_mat * (k_l - [i[1] - 1, i[2]  - 1, i[3] - 1])
end

"""
    size_q(prob)

Number of differences of k-points in the different directions.
"""
function size_q(prob::AbstractBSEProblem)
    size_q(size_k(prob))
end
function size_q(N_ks::Tuple{Int,Int,Int})
    return (n -> n == 1 ? 1 : 2 * n).(N_ks) # double size unless size is 1
end

"""
    q_lattice(prob)

Return a 2d array of size `(3, N_q)` where the columns are the lattice coordinates of the uniform ``q`` grid generated by a uniform ``k`` grid in one reciprocal unit cell with `size_k(prob)[i]` points in direction `i`.
"""

function q_lattice(prob::AbstractBSEProblem)
    return q_lattice(size_k(prob))
end
function q_lattice(N_ks::Tuple{Int,Int,Int})
    N_qs = size_q(N_ks)
    ql_1d_1 = vcat(range(0, stop = 1, length = N_ks[1] + 1)[1:(end - 1)], range(-1, stop = 0, length = N_ks[1] + 1)[(end - (N_qs[1] - N_ks[1])):(end - 1)])
    ql_1d_2 = vcat(range(0, stop = 1, length = N_ks[2] + 1)[1:(end - 1)], range(-1, stop = 0, length = N_ks[2] + 1)[(end - (N_qs[2] - N_ks[2])):(end - 1)])
    ql_1d_3 = vcat(range(0, stop = 1, length = N_ks[3] + 1)[1:(end - 1)], range(-1, stop = 0, length = N_ks[3] + 1)[(end - (N_qs[3] - N_ks[3])):(end - 1)])
    ql_1 = kron(ones(N_qs[3]), ones(N_qs[2]), ql_1d_1)
    ql_2 = kron(ones(N_qs[3]), ql_1d_2, ones(N_qs[1]))
    ql_3 = kron(ql_1d_3, ones(N_qs[2]), ones(N_qs[1]))
    ql = transpose(hcat(ql_1, ql_2, ql_3))
    return ql
end

"""
    Ω0_volume(prob)

Return the volume of the unit cell.
"""
function Ω0_volume(prob::AbstractBSEProblem)
    return abs(det(lattice_matrix(prob)))
end

"""
    BZ_volume(prob)

Return the volume of the first Brillouin zone. The convention is to use
the normalization b_mat^* * a_mat = 2π with the lattice vectors a_mat
and the reciprocal vectors b_mat.
"""
function BZ_volume(prob::AbstractBSEProblem)
    return (2 * π)^3 / Ω0_volume(prob)
end

"""
    compute_v_hat(prob[, cutoff])

Compute the Fourier coefficients of the Coulomb potential. The frequencies are stored in the order compatible with the 3d FFT.
"""
function compute_v_hat(prob::AbstractBSEProblem, cutoff = Inf)
    N_rs = size_r(prob)
    b_mat = 2 * π * inv(lattice_matrix(prob))
    v_hat = 4 * pi * vec(mapslices(x -> norm(x) < 1e-10 || norm(b_mat * x) > cutoff ? 0.0 : 1 / norm(b_mat * x)^2, fftfreq(N_rs...); dims = 1))
    return v_hat
end

"""
    fftfreq(n1, n2, n3)

Generate an 3d array of Fourier frequencies in the order that they appear in the 3d FFT.
"""
function fftfreq(n1, n2, n3)
    f1 = vcat(collect(0:(div(n1, 2) - (1 - n1 % 2))), collect(-div(n1, 2):-1))
    f2 = vcat(collect(0:(div(n2, 2) - (1 - n2 % 2))), collect(-div(n2, 2):-1))
    f3 = vcat(collect(0:(div(n3, 2) - (1 - n3 % 2))), collect(-div(n3, 2):-1))

    F1 = kron(ones(n3), ones(n2), f1)
    F2 = kron(ones(n3), f2, ones(n1))
    F3 = kron(f3, ones(n2), ones(n1))

    return vcat(F1', F2', F3')
end

"""
    G_vector_to_index(G, N_rs)

Return the index `i` such that `fftfreq(N_rs...)[:, i] == G`.
"""
function G_vector_to_index(G, N_rs) # TODO: make fast for all dimensions
    if length(G) == 3
        ind = mod1(G[1] + 1, N_rs[1]) + N_rs[1] * (mod1(G[2] + 1, N_rs[2]) - 1) + N_rs[1] * N_rs[2] * (mod1(G[3] + 1, N_rs[3]) - 1)
    else
        ind = mod1(G[1] + 1, N_rs[1])
        for id in 2:length(G)
            ind += prod(N_rs[1:(id - 1)]) * (mod1(G[id] + 1, N_rs[id]) - 1)
        end
    end
    return ind
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.1.3
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

# +
using Loess
using Random
using Statistics

Random.seed!(100)
xs = 10 .* rand(100)
ys = sin.(xs) .+ 0.5 * rand(100)

model = loess(xs, ys)

us = collect(minimum(xs):0.1:maximum(xs))
vs = predict(model, us);

# +
using Plots

plot( us, vs)
scatter!( xs, ys, markersize=2)
# -