text:
#layout
#auxiliary function to binary association
""" This function is used to define a list of logical elements
for use in order to define an algebra of layout function. It is
an internal function.
"""
function tree(v,niv,pos,lista,n)
	if niv<=n
		if pos==0
		aux=copy(v);
		aux[niv]=true
		niv=niv+1
		lista=push!(lista,aux)
		lista=tree(aux,niv,0,lista,n)
		lista=tree(aux,niv,1,lista,n)
		else
		niv=niv+1
		aux=copy(v);
		lista=tree(aux,niv,0,lista,n)
		lista=tree(aux,niv,1,lista,n)
		end
	end
	return lista
end

"""
This function compose the list given in tree function. This is
an internal function.
"""
function buildtree(n)
	v=falses(n)
	lista=[v]
	for i in tree(v,1,0,[v],n)[2:length(tree(v,1,0,[v],n))]
		lista=push!(lista,i)
	end
	for i in tree(v,1,1,[v],n)[2:length(tree(v,1,1,[v],n))]
		lista=push!(lista,i)
	end
	return lista
end

"""
This is an internal function used to print in REPL an UTF8 char
when the layout pretty option is trigged.
"""
function prettylayout(p,q)
	D=['₁','₂','₃','₄','₅','₆','₇','₈','₉']
	if 1<=p<10 && 1<=q<10
		return [D[p],D[q]]
	else
		if q==0
			return [D[p]]
		end
	end
end

""" This function is used to define the algebra  to work. If you
want use Liga, is mandatory to use this function.
## Example
```julia-repl
julia> layout(3,0,"GA")
```
generates a G3 algebra  with base 1,e1,e2,e3,e12,e13,e23,e123.
If you want to redefine the space, is highly recommended to restart
Julia.
"""
function layout(p::Int,q::Int,algebra::String; custom_repl="enable")
	# changing mode prompt

	eval(:( GIdLayout=($p,$q,$(algebra),$(2^(p+q))))) #Global Indentify Layout Parameters
	eval(:( LigaPrecision=10.0^(-10)))
	dim=p+q
	bn=buildtree(dim);
	#abn=BitArray{1}(dim)
	abn=BitArray(undef,dim,1)
	abnval=0.0
	#sorting bn
	bngr=zeros(2^dim)
	abngr=0
	#sorting bn
	bnval=zeros(2^dim)
	for i=1:2^dim
		bnval[i]=0.0
		bngr[i]=0
		for j=1:length(bn[i])
			if bn[i][j]==true
				bnval[i]+=2^(j-1)
				bngr[i]+=1
			end
		end
	end
	#sorting by grade
	for i=1:2^dim #sorting by grade
		for j=i+1:2^dim
			if bngr[i]>bngr[j]
				abnval=bnval[j]
				abn=bn[j]
				abngr=bngr[j]
				bnval[j]=bnval[i]
				bn[j]=bn[i]
				bngr[j]=bngr[i]
				bnval[i]=abnval
				bn[i]=abn
				bngr[i]=abngr
			end
		end
	end
	#sorting by value
	for i=1:2^dim
		for j=i+1:2^dim
			if bnval[i]>bnval[j] && bngr[i]==bngr[j]
				abnval=bnval[j]
				abn=bn[j]
				#abngr=bngr[j]
				bnval[j]=bnval[i]
				bn[j]=bn[i]
				#bngr[j]=bngr[i]
				bnval[i]=abnval
				bn[i]=abn
				#bngr[i]=abngr
			end
		end
	end

	if algebra  == "GA"
		if custom_repl=="enable"
			#OhMyREPL.colorscheme!("JuliaDefault")
			OhMyREPL.input_prompt!("𝔾$(String(prettylayout(p,q))) >", :green)
			#OhMyREPL.output_prompt!("𝔾$(String(prettylayout(p,q))) >", :yellow)
			#Base.active_repl.interface.modes[1].prompt=" 𝔾$(String(prettylayout(p,q))) > "
		else
			Base.active_repl.interface.modes[1].prompt="julia> "
		end
		bnew=Array{Vector{Bool},1}(undef,2^dim)
		ind=1
		for v in bn
			s=findall(x->x==true,v)
			if isempty(s)
				eval(Meta.parse(" id = BasisBlade(1) ;"));
				eval(Meta.parse("export id;"));
				bnew[ind]=v
				ind+=1
			else
				conc=string(s[1])
				for k=2:length(s)
					conc=string(conc,s[k])
				end
				eval(Meta.parse(" e$(conc) = BasisBlade($(ind)) ;"));
				eval(Meta.parse("export e$(conc);"));
				bnew[ind]=v
				ind+=1
			end
		end
		eval(Meta.parse("export LogicalBasisBlade  ;"));
		eval(Meta.parse(" LogicalBasisBlade = $(bnew);"));
		#Run operation table
		eval(Meta.parse("export BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ"));
		eval(Meta.parse("BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ=OperationTable()"));

	else
		if algebra == "Projective"
			if q != 0
				error("Something wrong with the second argument of layout(arg1,arg2,arg3), probably you need to setup equal 0.")
			else
				if custom_repl=="enable"
					#OhMyREPL.colorscheme!("JuliaDefault")
					OhMyREPL.input_prompt!("𝔾$(String(prettylayout(p,q))) [Proj]>", :green)
					#OhMyREPL.output_prompt!("𝔾$(String(prettylayout(p,q))) [Proj]>", :yellow)
					#Base.active_repl.interface.modes[1].prompt=" 𝔾$(String(prettylayout(p,q))) > "
				else
					Base.active_repl.interface.modes[1].prompt="julia> "
				end
				#printstyled("\n An apropriated embedding for Projective Space was created \n",color=:blue)
				bnew=Array{Vector{Bool},1}(undef,2^dim)
				ind=1
				for v in bn
					s=findall(x->x==true,v)
					if isempty(s)
						eval(Meta.parse(" id = BasisBlade(1) ;"));
						eval(Meta.parse("export id;"));
						bnew[ind]=v
						ind+=1
					else
						conc=string(s[1])
						for k=2:length(s)
							conc=string(conc,s[k])
						end
						eval(Meta.parse(" e$(conc) = BasisBlade($(ind)) ;"));
						eval(Meta.parse("export e$(conc);"));
						bnew[ind]=v
						ind+=1
					end
				end
				eval(Meta.parse("export LogicalBasisBlade  ;"));
				eval(Meta.parse(" LogicalBasisBlade = $(bnew);"));
				#Run operation table
				eval(Meta.parse("export BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ"));
				eval(Meta.parse("BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ=OperationTable()"));
			end
		end
		if algebra =="Conformal"
			if q != 1
				error("Something wrong with the second argument of layout(arg1,arg2,arg3), probably you need to setup equal 1.")
			else
				if custom_repl=="enable"
					OhMyREPL.input_prompt!("𝔾$(String(prettylayout(p,q)))[Conf]>", :green)
			#		OhMyREPL.output_prompt!("𝔾$(String(prettylayout(p,q)))[Conf]>", :yellow)
					#Base.active_repl.interface.modes[1].prompt=" 𝔾$(String(prettylayout(p,q))) > "
				else
					Base.active_repl.interface.modes[1].prompt="julia> "
				end
				#printstyled("\n An apropriated embedding for Conformal Space was created \n",color=:blue)
				#auxmask=zeros(2^dim)
				bnew=Array{Vector{Bool},1}(undef,2^dim)
				ind=1
				for v in bn
					s=findall(x->x==true,v)
					if isempty(s)
						eval(Meta.parse(" id = BasisBlade(1) ;"));
						eval(Meta.parse("export id;"));
						bnew[ind]=v
						ind+=1
					else
						if s[1]==dim-1
							conc="p"
						#auxmask[ind]+=1
						end
						if s[1]==dim
							conc="n"
							#auxmask[ind]+=2
						end
						if s[1]<dim-1
							conc=string(s[1])
						end
						for k=2:length(s)
							if s[k]==dim-1
								conc=string(conc,"p")
								#auxmask[ind]+=1
							end
							if s[k]==dim
								conc=string(conc,"n")
								#auxmask[ind]+=2
							end
							if s[k]<dim-1
								conc=string(conc,s[k])
							end
						end
						eval(Meta.parse(" e$(conc) = BasisBlade($(ind)) ;"));
						eval(Meta.parse("export e$(conc);"));
						bnew[ind]=v
						ind+=1
					end
				end
				#eval(Meta.parse("const vmask=$(auxmask)"))
				eval(Meta.parse("export LogicalBasisBlade  ;"));
				eval(Meta.parse(" LogicalBasisBlade = $(bnew);"));
				#Run operation table
				eval(Meta.parse("export BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ"));
				eval(Meta.parse("BBgeoprodTable,BBinnerprodTable,BBouterprodTable,Γ=OperationTable()"));
				eval(Meta.parse(" e∞ = 1.0*en+1.0*ep ;"));
				eval(Meta.parse("export e∞ "));
				eval(Meta.parse(" e0 = 0.5*en-0.5*ep ;"));
				eval(Meta.parse("export e0 "));
			end
		end
		if algebra!="Projective" && algebra!="Conformal"
			error("Not indentified algebra- Necessary redefine the layout!")
		end
	end

	#eval(:(const Layout=LayoutInfo([$p,$q],$(2^(p+q)),$(algebra ))))
	eval(Meta.parse(" I=BasisBlade(2^($(dim)))"));
	return ;
end

"""
This function show basic information after run layout function. For example if you
want do see basis elements and quantities related to created space is enough to follow
## Example
```julia-repl
julia> layoutinfo()
```
"""
function layoutinfo()
	printstyled(" ( p, q ) signature  = ( $(GIdLayout[1]),$(GIdLayout[2]) )\n",color=:blue)
	printstyled(" algebra   = $(GIdLayout[3])\n",color=:blue)
	printstyled(" dimension = $(GIdLayout[4])\n",color=:blue)
	printstyled(" identification of basis blade elements \n",color=:blue)
	for i=1:2^(GIdLayout[1]+GIdLayout[2])
		printstyled(" $(i) -> $(BasisBlade(i))\n",color=:blue )
	end
end


text:
# Starting example for built in functions in Learning Julia

# numeric functions
# round, floor, ceil, abs
x = 30.5
println(round(x, RoundUp))

y = 29.5
println(floor(y))
println(ceil(y))

z = -15
println(abs(z))

# text i/o
# print, println
print("Hello")
print(", ", "world")
println("!")

# printstyled
myStr = "This is some text"
printstyled(myStr, bold = true, color = :red)
println()

# read standard input
# readline
println("What is your name?")
name = readline()
print("My name is ", name, "\n")

# "is" functions
println(isascii("abc"))
println(isascii("αβγ"))

println(isdigit('9'))
println(isdigit('a'))

println(isspace(' '))
println(isspace('\r'))
println(isspace('\n'))
println(isspace('A'))

ext: jl

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
# Read massif.out.<pid> from `valgrind --tools=massif`
infile = ARGS[1]
content = open(infile) do file
    read(file, String)
end

function parse_block(s)
    lines = split(s, "\n"; limit=8)
    [split(x, "=")[2] for x in lines[3:7]]
end

blocks = split(content, "#-----------\nsnapshot=")
names = ["time", "mem_heap_B", "mem_heap_extra_B", "mem_stacks_B", "heap_tree"]
println(join(names, "\t"))
for block in blocks[2:end]
    values = parse_block(block)
    println(join(values, "\t"))
end

max_stars_repo_name: heavywatal/scribble

max_issues_repo_name: heavywatal/scribble

max_forks_repo_name: heavywatal/scribble


text:
struct Network
    neurons::Array{Float64}
    nin::Int64
    nout::Int64
    connections::BitArray
    weights::Array{Float64}
    ge::Array{Float64}
    gi::Array{Float64}
    trace::Array{Float64}
    exc::BitArray
    inh::BitArray
    inputs::BitArray
    outputs::BitArray
    da::Array{Float64}
    history::Dict
    cfg::Dict
end

function Network(n_neurons::Int64, n_input::Int64, n_output::Int64,
                 connectivity::Float64, cfg::Dict;
                 n_exc::Int64 = Int64(round(0.8 * n_neurons)))
    # starting membrane threshold for all neurons
    neurons = zeros(n_neurons, 3)
    neurons[:, 1] = cfg["vreset"]
    neurons[:, 2] = 0.0 # homeostasis
    neurons[:, 3] = 0.0 # refractory
    # index of exc neurons
    n_inh = n_neurons - n_exc
    exc = [trues(n_exc); falses(n_inh)]
    shuffle!(exc)
    inh = .~(exc)
    # set weights
    weights = cfg["wstart"] * ones(n_neurons, n_neurons)
    if cfg["wrandom"]
        weights .*= rand(n_neurons, n_neurons)
    end
    weights[inh, :] .= cfg["winh"]
    # weights .*= (1.0 - eye(n_neurons, n_neurons))
    connections = (rand(n_neurons, n_neurons) .< connectivity)
    weights .*= connections
    # set input neurons
    inputs = falses(n_neurons)
    inds = shuffle(find(exc))
    inputs[inds[1:n_input]] = true
    # set output neurons
    outputs = falses(n_neurons)
    inds = shuffle(find(exc .& .~(inputs)))
    outputs[inds[1:n_output]] = true
    trace = zeros(n_neurons, n_neurons)
    ge = zeros(n_neurons, n_neurons)
    gi = zeros(n_neurons, n_neurons)
    da = [0.0]
    history = Dict()
    history["spikes"] = BitArray(n_neurons, 0)
    Network(neurons, n_input, n_output, connections, weights, ge, gi, trace,
            exc, inh, inputs, outputs, da, history, cfg)
end

function input!(n::Network, input_spikes::BitArray)
    inds = find(n.inputs)[find(input_spikes)]
    n.neurons[inds, 1] .+= n.cfg["vinput"]
    nothing
end

function spike!(n::Network)
    # calculate spikes, reset spiked pre-synaptic membranes
    spikes = n.neurons[:, 1] .>= (n.cfg["vthresh"] .+ n.neurons[:, 2])
    spikes[n.exc] .&= (n.neurons[n.exc, 3] .> n.cfg["refrac_e"])
    spikes[n.inh] .&= (n.neurons[n.inh, 3] .> n.cfg["refrac_i"])
    # calculate conductance
    espikes = spikes .& n.exc
    n.ge[espikes, :] .+= n.weights[espikes, :]
    ispikes = spikes .& n.inh
    n.gi[ispikes, :] .+= n.weights[ispikes, :]
    # calculate post-synaptic membrane delta
    v = n.neurons[:, 1]
    dv = ((n.cfg["eleak"] - v) + (sum(n.ge, 1)' .* (n.cfg["eexc"] - v))
          + (sum(n.gi, 1)' .* (n.cfg["einh"] - v))
          + n.cfg["noise"] * rand(size(v)))
    dv[n.exc] /= n.cfg["tme"]
    dv[n.inh] /= n.cfg["tmi"]
    # update neuron state values
    nspikes = .~(spikes)
    n.neurons[spikes, 1] .= n.cfg["vreset"]
    n.neurons[nspikes, 1] .+= dv[nspikes]
    n.neurons[spikes, 2] .+= n.cfg["theta_plus"]
    n.neurons[nspikes, 2] .-= n.neurons[nspikes, 2] / n.cfg["ttheta"]
    n.neurons[spikes, 3] .= 0.0
    n.neurons[nspikes, 3] .+= 1.0
    n.trace[spikes, :] .+= 1
    n.trace .-= (n.trace / n.cfg["ttrace"])
    n.ge .-= n.ge / n.cfg["te"]
    n.gi .-= n.gi / n.cfg["ti"]
    # n.history["spikes"] = [n.history["spikes"] spikes]
    spikes
end

function learn!(n::Network, spikes::BitArray)
    if n.da[1] > 0.0
        dw = (exp(-n.cfg["lr"]) * n.da[1] * (n.trace[n.exc, spikes] - n.cfg["target"]) .*
              (n.cfg["wmax"] - n.weights[n.exc, spikes]).^n.cfg["mu"])
        n.weights[n.exc, spikes] .+= dw
        n.weights[n.weights .> n.cfg["wmax"]] .= n.cfg["wmax"]
        n.weights[n.weights .< 0.0] .= 0.0
        n.weights .*= n.connections
        n.da .-= n.da / n.cfg["tdop"]
    end
    nothing
end

function step!(n::Network, train::Bool=true)
    spikes = spike!(n)
    if train
        learn!(n, spikes)
    end
    spikes[n.outputs]
end

function step!(n::Network, input_spikes::BitArray, train::Bool=true)
    input!(n, input_spikes)
    step!(n, train)
end


text:
function something()::String
    return "test"
end

function lookup_and_write(values::Vector)
    output = nothing
    if length(values) == 1
        output = values[1]
    elseif length(values) == 2
        output = values[2]
    elseif length(values) == 3
        output = values[3]
    else
        output = values
    end
    return output
end

function lookup_and_write_without_else(values::Vector)
    output = nothing
    if length(values) == 1
        output = values[1]
    elseif length(values) == 2
        output = values[2]
    elseif length(values) == 3
        output = values[3]
    end
    return output
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(lookup_and_write([]) == [])
    @assert(lookup_and_write(convert(Vector, [1])) == 1)
    @assert(lookup_and_write(convert(Vector, [1, 2])) == 2)
    @assert(lookup_and_write(convert(Vector, [1, 2, 3])) == 3)
    @assert(lookup_and_write(convert(Vector, [1, 2, 3, 4])) == [1, 2, 3, 4])
    @assert(lookup_and_write_without_else([]) === nothing)
    @assert(lookup_and_write_without_else(convert(Vector, [1])) == 1)
    @assert(lookup_and_write_without_else(convert(Vector, [1, 2])) == 2)
    @assert(lookup_and_write_without_else(convert(Vector, [1, 2, 3])) == 3)
    @assert(lookup_and_write_without_else(convert(Vector, [1, 2, 3, 4])) === nothing)
end


text:
import Base: setindex!, reduce, select
import DataValues: dropna
export NextTable, table, colnames, pkeynames, columns, pkeys, reindex, dropna

"""
A permutation

# Fields:

- `columns`: The columns being indexed as a vector of integers (column numbers)
- `perm`: the permutation - an array or iterator which has the sorted permutation
"""
struct Perm{X}
    columns::Vector{Int}
    perm::X
end

abstract type AbstractIndexedTable end

"""
A table.
"""
struct NextTable{C<:Columns} <: AbstractIndexedTable
    # `Columns` object which iterates to give an array of rows
    columns::C
    # columns that are primary keys (Vector{Int})
    pkey::Vector{Int}
    # Cache permutations by various subsets of columns
    perms::Vector{Perm}
    # store what percent of the data in each column is unique
    cardinality::Vector{Nullable{Float64}}

    columns_buffer::Any
end

"""
`table(cols::AbstractVector...; names, <options>)`

Create a table with columns given by `cols`.
```jldoctest table
julia> a = table([1,2,3], [4,5,6])
Table with 3 rows, 2 columns:
1  2
────
1  4
2  5
3  6
```

`names` specify names for columns. If specified, the table will be an iterator of named tuples.

```jldoctest table
julia> b = table([1,2,3], [4,5,6], names=[:x, :y])
Table with 3 rows, 2 columns:
x  y
────
1  4
2  5
3  6

```

`table(cols::Union{Tuple, NamedTuple}; <options>)`

Convert a struct of columns to a table of structs.

```jldoctest table
julia> table(([1,2,3], [4,5,6])) == a
true

julia> table(@NT(x=[1,2,3], y=[4,5,6])) == b
true
```

`table(cols::Columns; <options>)`

Construct a table from a vector of tuples. See [`rows`](@ref).

```jldoctest table
julia> table(Columns([1,2,3], [4,5,6])) == a
true

julia> table(Columns(x=[1,2,3], y=[4,5,6])) == b
true
```

`table(t::Union{Table, NDSparse}; <options>)`

Copy a Table or NDSparse to create a new table. The same primary keys as the input are used.

```jldoctest table
julia> b == table(b)
true
```

`table(iter; <options>)`

Construct a table from an iterable table.

# Options:

- `pkey`: select columns to act as the primary key. By default, no columns are used as primary key.
- `presorted`: is the data pre-sorted by primary key columns? If so, skip sorting. `false` by default. Irrelevant if `chunks` is specified.
- `copy`: creates a copy of the input vectors if `true`. `true` by default. Irrelavant if `chunks` is specified.
- `chunks`: distribute the table into `chunks` (Integer) chunks (a safe bet is nworkers()). Table is not distributed by default. See [Distributed](@distributed) docs.

# Examples:

Specifying `pkey` will cause the table to be sorted by the columns named in pkey:

```jldoctest table
julia> b = table([2,3,1], [4,5,6], names=[:x, :y], pkey=:x)
Table with 3 rows, 2 columns:
x  y
────
1  6
2  4
3  5

julia> b = table([2,1,2,1],[2,3,1,3],[4,5,6,7],
                 names=[:x, :y, :z], pkey=(:x,:y))
Table with 4 rows, 3 columns:
x  y  z
───────
1  3  5
1  3  7
2  1  6
2  2  4
```
Note that the keys do not have to be unique.

`chunks` option creates a distributed table.

`chunks` can be:

1. An integer -- number of chunks to create
2. An vector of `k` integers -- number of elements in each of the `k` chunks.
3. The distribution of another array. i.e. `vec.subdomains` where `vec` is a distributed array.

```jldoctest table
julia> t = table([2,3,1,4], [4,5,6,7],
                  names=[:x, :y], pkey=:x, chunks=2)
Distributed Table with 4 rows in 2 chunks:
x  y
────
1  6
2  4
3  5
4  7
```

A distributed table will be constructed if one of the arrays passed into `table` constructor is a distributed
array. A distributed Array can be constructed using `distribute`:

```jldoctest table

julia> x = distribute([1,2,3,4], 2);

julia> t = table(x, [5,6,7,8], names=[:x,:y])
Distributed Table with 4 rows in 2 chunks:
x  y
────
1  5
2  6
3  7
4  8

julia> table(columns(t)..., [9,10,11,12],
             names=[:x,:y,:z])
Distributed Table with 4 rows in 2 chunks:
x  y  z
────────
1  5  9
2  6  10
3  7  11
4  8  12

```

Distribution is done to match the first distributed column from left to right. Specify `chunks` to override this.
"""
function table end

function table(::Val{:serial}, cols::Tup;
               pkey=Int[],
               chunks=nothing, # unused here
               perms=Perm[],
               presorted=false,
               copy=true,
               cardinality=fill(Nullable{Float64}(), length(cols)))

    cs = rows(cols)

    if isa(pkey, Union{Int, Symbol})
        pkey = [pkey]
    elseif isa(pkey, Tuple)
        pkey = [pkey...]
    end

    if !presorted && !isempty(pkey)
        pkeys = rows(cs, (pkey...))
        if !issorted(pkeys)
            perm = sortperm(pkeys)
            if copy
                cs = cs[perm]
            else
                cs = permute!(cs, perm)
            end
        elseif copy
            cs = Base.copy(cs)
        end
    elseif copy
        cs = Base.copy(cs)
    end

    intpkey = map(k->colindex(cs, k), pkey)

    NextTable{typeof(cs)}(cs,
           intpkey,
           perms,
           cardinality,
           similar(cs, 0))
end

function table(::Val{impl}, cols; kwargs...) where impl
    if impl == :distributed && isa(cols, Tup)
        error("""You requested to create a distributed table.
                 Distributed table is implemented by JuliaDB.
                 run `using JuliaDB` and try again.""")
    else
        error("unknown table implementation invoked")
    end
end

# detect if a distributed table has to be constructed.
_impl(impl::Val) = impl
_impl(impl::Val, x::AbstractArray, z...) = _impl(impl, z...)
_impl(x::AbstractArray...) = _impl(Val{:serial}(), x...)

function table(cs::Tup; chunks=nothing, kwargs...)
    if chunks !== nothing
        impl = Val{:distributed}()
    else
        impl = _impl(astuple(cs)...)
    end
    table(impl, cs; chunks=chunks, kwargs...)
end

table(cs::Columns; kwargs...) = table(columns(cs); kwargs...)
table(c::Columns{<:Pair}; kwargs...) = convert(NextTable, c.columns.first, c.columns.second; kwargs...)

function table(cols::AbstractArray...; names=nothing, kwargs...)
    if isa(names, AbstractArray) && all(x->isa(x, Symbol), names)
        cs = namedtuple(names...)(cols...)
    else
        cs = cols
    end
    table(cs; kwargs...)
end

# Easy constructor to create a derivative table
function table(t::NextTable;
               columns=t.columns,
               chunks=nothing,
               pkey=t.pkey,
               perms=t.perms,
               cardinality=t.cardinality,
               presorted=false,
               copy=true)

    table(columns,
          pkey=pkey,
          perms=perms,
          chunks=chunks,
          cardinality=cardinality,
          presorted=presorted,
          copy=copy)
end

Base.@pure colnames(t::AbstractIndexedTable) = fieldnames(eltype(t))
columns(t::NextTable) = columns(t.columns)
# throw a better error message when a custom array
# of different size is used
function column(t::NextTable, a::AbstractArray)
    if length(t) != length(a)
        throw(ArgumentError(
                "vector provided must have the same length as table"
             )
        )
    end
    column(rows(t), a)
end

Base.eltype(::Type{NextTable{C}}) where {C} = eltype(C)
Base.eltype(t::NextTable) = eltype(typeof(t))
Base.copy(t::NextTable) = table(t, copy=true)
function Base.empty!(t::NextTable)
    empty!(t.perms)
    empty!(rows(t))
    t
end
Base.:(==)(a::NextTable, b::NextTable) = rows(a) == rows(b)

Base.getindex(t::NextTable, i::Integer) = getindex(t.columns, i)
Base.getindex(t::NextTable, i::Colon) = copy(t)
Base.endof(t::NextTable) = length(t)

function Base.view(t::NextTable, I)
    t.pkey == Int64[] || eltype(I) == Bool || issorted(I) ||
        throw(ArgumentError("`view` called with unsorted index."))
    table(
        view(t.columns, I),
        pkey = t.pkey,
        copy = false,
        presorted = true)
end

Base.length(t::NextTable) = length(t.columns)
Base.start(t::NextTable) = start(t.columns)
Base.next(t::NextTable, i) = next(t.columns, i)
Base.done(t::NextTable, i) = done(t.columns, i)
function getindex(t::NextTable, idxs::AbstractVector{<:Integer})
    if t.pkey == Int64[] || eltype(idxs) == Bool || issorted(idxs)
       #perms = map(t.perms) do p
       #    # TODO: make the ranks continuous
       #    Perm(p.columns, p.perm[idxs])
       #end
        perms = Perm[]
        table(t, columns=t.columns[idxs], perms=perms, copy=false, presorted=true)
    else
        # this is for gracefully allowing this later
        throw(ArgumentError("`getindex` called with unsorted index. This is not allowed at this time."))
    end
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable}, key::Tuple)
    t = d[]
    idx = [colindex(t, k) for k in key]
    pkey = Int[]
    for (i, pk) in enumerate(t.pkey)
        j = findfirst(idx, pk)
        if j > 0
            push!(pkey, j)
        end
    end
    table(d.src, columns=columns(t, key), pkey=pkey)
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable}, key::SpecialSelector)
    getindex(d, lowerselection(d[], key))
end

function ColDict(t::AbstractIndexedTable; copy=nothing)
    ColDict(Base.copy(t.pkey), t,
            Base.copy(colnames(t)), Any[columns(t)...], copy)
end

function Base.getindex(d::ColDict{<:AbstractIndexedTable})
    table(d.columns...;
          names=d.names,
          copy=d.copy === nothing ? false : d.copy,
          pkey=d.pkey)
end

function subtable(t::Union{Columns, NextTable}, idxs; presorted=true)
    table(t, columns=rows(t)[idxs], perms=t.perms, copy=false, presorted=presorted)
end

function primaryperm(t::NextTable)
    Perm(t.pkey, Base.OneTo(length(t)))
end

permcache(t::NextTable) = [primaryperm(t), t.perms;]
cacheperm!(t::NextTable, p) = push!(t.perms, p)

"""
`pkeynames(t::Table)`

Names of the primary key columns in `t`.

# Example

```jldoctest pkeynames

julia> t = table([1,2], [3,4]);

julia> pkeynames(t)
()

julia> t = table([1,2], [3,4], pkey=1);

julia> pkeynames(t)
(1,)

julia> t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2));

julia> pkeys(t)
2-element IndexedTables.Columns{NamedTuples._NT_x_y{Int64,Int64},NamedTuples._NT_x_y{Array{Int64,1},Array{Int64,1}}}:
 (x = 1, y = 3)
 (x = 2, y = 1)

```
"""
function pkeynames(t::AbstractIndexedTable)
    if eltype(t) <: NamedTuple
        (colnames(t)[t.pkey]...)
    else
        (t.pkey...)
    end
end

# for a table, selecting the "value" means selecting all fields
valuenames(t::AbstractIndexedTable) = (colnames(t)...)

"""
    pkeys(itr::Table)

Primary keys of the table. If Table doesn't have any designated
primary key columns (constructed without `pkey` argument) then
a default key of tuples `(1,):(n,)` is generated.

# Example

```jldoctest pkeys

julia> a = table(["a","b"], [3,4]) # no pkey
Table with 2 rows, 2 columns:
1    2
──────
"a"  3
"b"  4

julia> pkeys(a)
2-element Columns{Tuple{Int64}}:
 (1,)
 (2,)

julia> a = table(["a","b"], [3,4], pkey=1)
Table with 2 rows, 2 columns:
1    2
──────
"a"  3
"b"  4

julia> pkeys(a)
2-element Columns{Tuple{String}}:
 ("a",)
 ("b",)
```

"""
function pkeys(t::NextTable)
    if isempty(t.pkey)
        Columns(Base.OneTo(length(t)))
    else
        rows(t, pkeynames(t))
    end
end

Base.values(t::NextTable) = rows(t)

"""
`sort(t[, by::Selection]; select::Selection, kwargs...)`

Sort rows by `by`. All of `Base.sort` keywords can be used.

# Examples

```jldoctest sort
julia> t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6],
names=[:x,:y,:z]);

julia> sort(t, :z; select = (:y, :z), rev = true)
Table with 6 rows, 2 columns:
y  z
────
1  6
1  5
2  4
2  3
1  2
1  1
```
"""
sort(t::NextTable, by...; select = valuenames(t), kwargs...) =
    table(rows(t, select)[sortperm(rows(t, by...); kwargs...)], copy = false)

"""
`sort!(t[, by::Selection]; kwargs...)`

Sort rows by `by` in place. All of `Base.sort` keywords can be used.

# Examples

```jldoctest sort!
julia> t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6],
names=[:x,:y,:z]);

julia> sort!(t, :z, rev = true);

julia> t
Table with 6 rows, 3 columns:
x  y  z
───────
2  1  6
2  1  5
2  2  4
1  2  3
1  1  2
1  1  1
"""
function sort!(t::NextTable, by...; kwargs...)
    isempty(t.pkey) || error("Tables with primary keys can't be sorted in place")
    permute!(rows(t), sortperm(rows(t, by...); kwargs...))
    t
end

"""
    excludecols(itr, cols)

Names of all columns in `itr` except `cols`. `itr` can be any of
`Table`, `NDSparse`, `Columns`, or `AbstractVector`

```jldoctest excludecols
julia> t = table([2,1],[1,3],[4,5], names=[:x,:y,:z], pkey=(1,2))
Table with 2 rows, 3 columns:
x  y  z
───────
1  3  5
2  1  4

julia> excludecols(t, (:x,))
(:y, :z)

julia> excludecols(t, (2,))
(:x, :z)

julia> excludecols(t, pkeynames(t))
(:z,)

julia> excludecols([1,2,3], (1,))
()

```
"""
function excludecols(t, cols)
    if cols isa SpecialSelector
        return excludecols(t, lowerselection(t, cols))
    end
    if !isa(cols, Tuple)
        return excludecols(t, (cols,))
    end
    ns = colnames(t)
    mask = ones(Bool, length(ns))
    for c in cols
        i = colindex(t, c)
        if i !== 0
            mask[i] = false
        end
    end
    ((1:length(ns))[mask]...)
end

"""
    convert(NextTable, pkeys, vals; kwargs...)

Construct a table with `pkeys` as primary keys and `vals` as corresponding non-indexed items.
keyword arguments will be forwarded to [`table`](@ref) constructor.

# Example

```jldoctest
julia> convert(NextTable, Columns(x=[1,2],y=[3,4]), Columns(z=[1,2]), presorted=true)
Table with 2 rows, 3 columns:
x  y  z
───────
1  3  1
2  4  2
```
"""
function convert(::Type{NextTable}, key, val; kwargs...)
    cs = concat_cols(key, val)
    table(cs, pkey=[1:ncols(key);]; kwargs...)
end

convert(T::Type{NextTable}, c::Columns{<:Pair}; kwargs...) = convert(T, c.columns.first, c.columns.second; kwargs...)
# showing

import Base.Markdown.with_output_format

global show_compact_when_wide = true
function set_show_compact!(flag=true)
    global show_compact_when_wide
    show_compact_when_wide = flag
end

function showtable(io::IO, t; header=nothing, cnames=colnames(t), divider=nothing, cstyle=[], full=false, ellipsis=:middle, compact=show_compact_when_wide)
    height, width = displaysize(io)
    showrows = height-5 - (header !== nothing)
    n = length(t)
    header !== nothing && println(io, header)
    if full
        rows = [1:n;]
        showrows = n
    else
        if ellipsis == :middle
            lastfew = div(showrows, 2) - 1
            firstfew = showrows - lastfew - 1
            rows = n > showrows ? [1:firstfew; (n-lastfew+1):n] : [1:n;]
        elseif ellipsis == :end
            lst = n == showrows ?
                showrows : showrows-1 # make space for ellipse
            rows = [1:min(length(t), showrows);]
        else
            error("ellipsis must be either :middle or :end")
        end
    end
    nc = length(columns(t))
    reprs  = [ sprint(io->showcompact(io,columns(t)[j][i])) for i in rows, j in 1:nc ]
    strcnames = map(string, cnames)
    widths  = [ max(strwidth(get(strcnames, c, "")), isempty(reprs) ? 0 : maximum(map(strwidth, reprs[:,c]))) for c in 1:nc ]
    if compact && !isempty(widths) && sum(widths) + 2*nc > width
        return showmeta(io, t, cnames)
    end
    for c in 1:nc
        nm = get(strcnames, c, "")
        style = get(cstyle, c, nothing)
        txt = c==nc && divider!=nc ? nm : rpad(nm, widths[c]+(c==divider ? 1 : 2), " ")
        if style == nothing
            print(io, txt)
        else
            with_output_format(style, print, io, txt)
        end
        if c == divider
            print(io, "│")
            length(cnames) > divider && print(io, " ")
        end
    end
    println(io)
    if divider !== nothing
        print(io, "─"^(sum(widths[1:divider])+2*divider-1), "┼")
        if !isempty(widths[divider+1:end])
            print(io, "─"^(sum(widths[divider+1:end])+2*(nc-divider)-1))
        end
    else
        if !isempty(widths)
            print(io, "─"^(sum(widths)+2*nc-2))
        end
    end
    for r in 1:size(reprs,1)
        println(io)
        for c in 1:nc
            print(io, c==nc && nc!=divider ? reprs[r,c] : rpad(reprs[r,c], widths[c]+(c==divider ? 1 : 2), " "))
            if c == divider
                print(io, "│ ")
            end
        end
        if n > showrows && ((ellipsis == :middle && r == firstfew) || (ellipsis == :end && r == size(reprs, 1)))
            if divider === nothing
                println(io)
                print(io, "⋮")
            else
                println(io)
                print(io, " "^(sum(widths[1:divider]) + 2*divider-1), "⋮")
            end
        end
    end
end

function showmeta(io, t, cnames)
    nc = length(columns(t))
    println(io, "Columns:")
    metat = Columns(([1:nc;], [Text(string(get(cnames, i, "<noname>"))) for i in 1:nc],
                     [map(eltype, columns(t))...]))
    showtable(io, metat, cnames=["#", "colname", "type"], cstyle=fill(:bold, nc), full=true, compact=false)
end

function subscriptprint(x::Integer)
    s = string(x)
    cs = Char[]
    lookup = ["₀₁₂₃₄₅₆₇₈₉"...]
    join([lookup[parse(Int, c)+1] for c in s],"")
end

function show(io::IO, t::NextTable{T}) where {T}
    header = "Table with $(length(t)) rows, $(length(columns(t))) columns:"
    cstyle = Dict([i=>:bold for i in t.pkey])
    cnames = string.(colnames(t))
    for (i, k) in enumerate(t.pkey)
        cstyle[k] = :bold
        #cnames[k] = cnames[k] * "$(subscriptprint(i))"
    end
    showtable(io, t, header=header, cnames=cnames, cstyle=cstyle)
end


text:
using DReal
using Compat
import DReal: opensmt_context, opensmt_set_verbosity, opensmt_assert, opensmt_check, opensmt_print_expr, opensmt_mk_context,
              opensmt_mk_bool_var, opensmt_mk_int_var, opensmt_mk_minus, opensmt_mk_leq, opensmt_mk_num_from_string,
              opensmt_push, opensmt_pop, opensmt_reset, opensmt_del_context, opensmt_set_precision, opensmt_mk_real_var

import DReal: opensmt_mk_geq,opensmt_mk_leq

i = 0;
ctx = opensmt_mk_context(@compat UInt32(1));
opensmt_set_precision (ctx, 0.0000001);
x = opensmt_mk_real_var(ctx, "xleaks4" , 0.0, 1.0);
y = opensmt_mk_real_var(ctx, "yleaks4" , 0.0, 1.0);

for i=1:1000
    opensmt_push(ctx);
    lb = opensmt_mk_num_from_string(ctx, "$(rand())");
    ub = opensmt_mk_num_from_string(ctx, "$(rand())");

    if rand() > 0.5
    # print(lb,ub)
        geq = opensmt_mk_geq(ctx, x, y);
        opensmt_assert(ctx, geq);
    else
        leq = opensmt_mk_leq(ctx, x, y);
        opensmt_assert(ctx, leq);
    end
    res = opensmt_check( ctx );
    opensmt_pop(ctx);
end
opensmt_del_context(ctx);


text:
# ground state tests

import FourierGPE:V
V(x,t) = 0.5*x^2

# ==== set simulation parameters ====
L = (40.0,)
N = (512,)
μ = 25.0

# === collect parameters ===========
sim = Sim(L,N)
@pack! sim = μ
@unpack_Sim sim
# ===================================

# useful TF state
ψ0(x,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,0.0)/μ,0.0)+im*0.0)

#evaluate on grid
x = X[1]
ψi = ψ0.(x,μ,g)

#make initial state
ϕi = kspace(ψi,sim)
@pack! sim = ϕi

# ====== Evolve in k space without error ==========

runsim(sim;info=false)

sol,err = testsim(sim)
@test err == false

psitest = xspace(sol[end],sim)
npeak = abs2.(psitest[256])
@test isapprox(g*npeak,μ,rtol=1e-2)


text:
"""
    PolyhedraToLPBridge{T}

The `PolyhedraToLPBridge` converts a constraint `VF`-in-`PolyhedraOptSet`
into the constraints `F`-in-`EqualTo` for the hyperplanes and `F`-to-`LessThan`
for halfspaces.
"""
struct PolyhedraToLPBridge{T, F} <: MOI.Bridges.Constraint.AbstractBridge
    hyperplanes::Vector{MOI.ConstraintIndex{F, MOI.EqualTo{T}}}
    halfspaces::Vector{MOI.ConstraintIndex{F, MOI.LessThan{T}}}
end
function MOI.Bridges.Constraint.bridge_constraint(
    ::Type{PolyhedraToLPBridge{T, F}}, model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction, p::PolyhedraOptSet) where {T, F}

    vf = MOIU.eachscalar(f)
    hps = [
        MOIU.normalize_and_add_constraint(model, func(T, h.a, vf, F), MOI.EqualTo(convert(T, h.β)))
        for h in hyperplanes(p.rep)
    ]
    hss = [
        MOIU.normalize_and_add_constraint(model, func(T, h.a, vf, F), MOI.LessThan(convert(T, h.β)))
        for h in halfspaces(p.rep)
    ]
    return PolyhedraToLPBridge{T, F}(hps, hss)
end

function func(T::Type, a::AbstractVector, vf, F::Type)
    func = zero(F)
    for (α, f) in zip(a, vf)
        MOIU.operate!(+, T, func, MOIU.operate!(*, T, f, convert(T, α)))
    end
    return func
end

MOI.supports_constraint(::Type{PolyhedraToLPBridge{T}},
                        ::Type{<:MOI.AbstractVectorFunction},
                        ::Type{<:PolyhedraOptSet}) where {T} = true
function MOI.Bridges.added_constrained_variable_types(::Type{<:PolyhedraToLPBridge})
    return Tuple{DataType}[]
end
function MOI.Bridges.added_constraint_types(::Type{PolyhedraToLPBridge{T, F}}) where {T, F}
    return [(F, MOI.EqualTo{T}), (F, MOI.LessThan{T})]
end
function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{<:PolyhedraToLPBridge{T}},
    VF::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:PolyhedraOptSet}) where T

    SF = MOIU.scalar_type(VF)
    TermType = MOIU.promote_operation(*, T, T, SF)
    F = MOIU.promote_operation(+, T, TermType, TermType)
    return PolyhedraToLPBridge{T, F}
end

# Attributes, Bridge acting as an model
function MOI.get(b::PolyhedraToLPBridge{T, F}, ::MOI.NumberOfConstraints{F, MOI.EqualTo{T}}) where {T, F}
    return length(b.hyperplanes)
end
function MOI.get(b::PolyhedraToLPBridge{T, F}, ::MOI.ListOfConstraintIndices{F, MOI.EqualTo{T}}) where {T, F}
    return b.hyperplanes
end
function MOI.get(b::PolyhedraToLPBridge{T, F}, ::MOI.NumberOfConstraints{F, MOI.LessThan{T}}) where {T, F}
    return length(b.halfspaces)
end
function MOI.get(b::PolyhedraToLPBridge{T, F}, ::MOI.ListOfConstraintIndices{F, MOI.LessThan{T}}) where {T, F}
    return b.halfspaces
end

# Indices
function MOI.delete(model::MOI.ModelLike, b::PolyhedraToLPBridge)
    for h in b.hyperplanes
        MOI.delete(model, h)
    end
    for h in b.halfspaces
        MOI.delete(model, h)
    end
end


text:
# Normalizing the data after pre-processing and before post-processing.

abstract type Scaling end

# scales the data so that it ranges from 0 to 1.
struct min_max_scaling <: Scaling
    Δv::Float64 # max - min
    v_min::Float64
end

# scale: scales the data so that it ranges from 0 to 1.
scale(x, scaling::min_max_scaling) = (x .- scaling.v_min) ./ scaling.Δv

# unscale: undo scale
unscale(x, scaling::min_max_scaling) = (x .* scaling.Δv) .+ scaling.v_min

# struct min_max_scaling <: Scaling
#     ΔT::Float64 # range of Temp values in initial profile
#     T_max::Float64 # max value of the initial temperature profile (before cooling)
# end
#
# #scale: normalize the data so that it ranges from [-1,0]
# scale(x, scaling::min_max_scaling) = (x .- scaling.T_max) ./ scaling.ΔT
# #unscale: undo normalization
# unscale(x, scaling::min_max_scaling) = x .* scaling.ΔT .+ scaling.T_max

# struct wTscaling <: Scaling end
#
# #scale: divide all wT values by the maximum value across all timesteps and gridpoints
# scale(x, scaling::wTscaling) = (x .+ 1e-5) ./ 4e-5
# #scale: undo scale
# unscale(x, scaling::wTscaling) = x .* 4e-5 .- 1e-5

# struct wTscaling <: Scaling
#     nc::Float64 #normalization constant
# end
#
# #scale: divide all wT values by the maximum value across all timesteps and gridpoints
# scale(x, scaling::wTscaling) = x ./ scaling.nc
# #scale: undo scale
# unscale(x, scaling::wTscaling) = x .* scaling.nc

# function get_scaling(V_name, vavg)
#     if V_name=="T"; return Tscaling(vavg[1][end]-vavg[1][1], maximum(vavg[1])) # Tscaling(ΔT, T_max)
#     elseif V_name=="wT"; return wTscaling( maximum(maximum, vavg) ) # wTscaling(nc)
#     else; throw(error()) end
# end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
abstract type AbstractReca <: AbstractReservoirComputer end

struct RECA{S,R,E,T,Q} <: AbstractReca
    #res_size::I
    train_data::S
    automata::R
    input_encoding::E
    nla_type::ReservoirComputing.NonLinearAlgorithm
    states::T
    states_type::Q
end

"""
    RECA(train_data,
        automata;
        generations = 8,
        input_encoding=RandomMapping(),
        nla_type = NLADefault(),
        states_type = StandardStates())

[1] Yilmaz, Ozgur. “_Reservoir computing using cellular automata._”
arXiv preprint arXiv:1410.0162 (2014).

[2] Nichele, Stefano, and Andreas Molund. “_Deep reservoir computing using cellular
automata._” arXiv preprint arXiv:1703.02806 (2017).
"""
function RECA(
    train_data,
    automata;
    generations = 8,
    input_encoding=RandomMapping(),
    nla_type = NLADefault(),
    states_type = StandardStates()
)
    in_size = size(train_data, 1)
    #res_size = obtain_res_size(input_encoding, generations)
    state_encoding = create_encoding(input_encoding, train_data, generations)
    states = reca_create_states(state_encoding, automata, train_data)

    return RECA(train_data, automata, state_encoding, nla_type, states, states_type)
end

#training dispatch
function train(reca::AbstractReca, target_data, training_method=StandardRidge(0.0))
    states_new = reca.states_type(reca.nla_type, reca.states, reca.train_data)
    return _train(states_new, target_data, training_method)
end

#predict dispatch
function (reca::RECA)(
    prediction,
    output_layer::AbstractOutputLayer,
    initial_conditions=output_layer.last_value,
    last_state=zeros(reca.input_encoding.ca_size)
)
    return obtain_prediction(
        reca, prediction, last_state, output_layer; initial_conditions=initial_conditions
    )
end

function next_state_prediction!(reca::RECA, x, out, i, args...)
    rm = reca.input_encoding
    x = encoding(rm, out, x)
    ca = CellularAutomaton(reca.automata, x, rm.generations+1)
    ca_states = ca.evolution[2:end ,:]
    x_new = reshape(transpose(ca_states), rm.states_size)
    x = ca.evolution[end, :]
    return x, x_new
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
module StochasticIntegrals

using DataFrames
using DataStructures: OrderedDict
using Dates
using Distributions
using LinearAlgebra: Hermitian, cholesky, inv, det, LowerTriangular, diag
using UnivariateFunctions
using Random
using Sobol: SobolSeq, next!
using FixedPointAcceleration

include("4_number_generators.jl")
export NumberGenerator, Mersenne, SobolGen, next!
include("1_main_functions.jl")
export ItoIntegral, volatility, variance, covariance, correlation, ItoSet, ForwardCovariance
export get_draws, get_zero_draws, pdf, make_covariance_matrix, log_likelihood, brownians_in_use
export StochasticIntegralsCovariance, SimpleCovariance, update!
include("2_data_conversions.jl")
export to_draws, to_dataframe, to_array
include("3_getConfidenceHypercube.jl")
export get_confidence_hypercube
include("5_conditional_distributions.jl")
export generate_conditioned_distribution
include("6_ItoProcesses.jl")
export ItoProcess, evolve!, evolve_covar_and_ito_processes!, make_ito_process_syncronous_time_series, make_ito_process_non_syncronous_time_series

end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["Build"] = BenchmarkGroup()

# ==============================================================================
# Decomposition-based approach
# ==============================================================================
include("building_GLGM06.jl")

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# # tune parameters
# tune!(SUITE)
#
# # run the benchmarks
# results = run(SUITE, verbose=false)
#
# # return the sample with the smallest time value in each test
# println("minimum: ", minimum(results))
#
# # return the median for each test
# println("median: ", median(results))

# ==============================================================================
# Create plots
# ==============================================================================

function plot_dashed_line!(v, time_horizon=20.)
    plot!(x->x, x->v, 0.0, time_horizon, line=2, color="red", linestyle=:dash,
          legend=nothing)
end

plot_vars = (0, 25)

res = solve(build_TV, alg=algo_dense, T=time_horizon)

plot(res,
     vars=plot_vars,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t",
     ylab=L"x_{25}",
     xtick=[0., 0.5, 1.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
     xlims=(0., 1.), ylims=(-9e-3, 6e-3),
     bottom_margin=0mm, left_margin=1mm, right_margin=6mm, top_margin=3mm,
     size=(1000, 1000))
plot_dashed_line!(4e-3)
plot_dashed_line!(5.1e-3)
savefig(@relpath "BLDF01_time_horizon_1.png")

# plot(res,
#      vars=plot_vars,
#      tickfont=font(30, "Times"), guidefontsize=45,
#      xlab=L"t",
#      ylab=L"x_{25}",
#      xtick=[0., 10., 20.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
#      bottom_margin=6mm, left_margin=2mm, top_margin=3mm,
#      ylims=(-Inf, 6e-3),
#      size=(1000, 1000))
# plot_dashed_line!(4e-3)
# plot_dashed_line!(5.1e-3)
# savefig(@relpath "BLDF01_time_horizon_20.png")
#
# res = solve(build_CONST, algo=algo_dense, T=time_horizon)
#
# plot(res,
#      vars=plot_vars,
#      tickfont=font(30, "Times"), guidefontsize=45,
#      xlab=L"t",
#      ylab=L"x_{25}",
#      xtick=[0., 0.5, 1.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
#      xlims=(0., 1.), ylims=(-Inf, 6e-3),
#      bottom_margin=0mm, left_margin=1mm, right_margin=6mm, top_margin=3mm,
#      size=(1000, 1000))
# plot_dashed_line!(4e-3)
# plot_dashed_line!(5.1e-3)
# savefig(@relpath "BLDC01_time_horizon_1.png")
#
# plot(res,
#      vars=plot_vars,
#      tickfont=font(30, "Times"), guidefontsize=45,
#      xlab=L"t",
#      ylab=L"x_{25}",
#      xtick=[0., 10., 20.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
#      bottom_margin=6mm, left_margin=2mm, top_margin=3mm,
#      ylims=(-Inf, 6e-3),
#      size=(1000, 1000))
# plot_dashed_line!(4e-3)
# plot_dashed_line!(5.1e-3)
# savefig(@relpath "BLDC01_time_horizon_20.png")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
"""
`module TriangularLattice`

The triangular lattice primarily serves as a useful toy model.
"""
module TriangularLattice

using JuLIP

function bulktri(; a0 = 1.0, L = 1)
   X = a0 * [ JVec([0.0, 0.0, 0.0]), JVec([0.5, √3/2, 0.0]) ]
   F = [ 1.0 0.0 0.0; 0.0 √3 0.0; 0.0 0.0 1.0 ]
   at = ASEAtoms("Cu2")
   set_positions!(at, X)
   set_defm!(at, F)

   if L == 1
      set_pbc!(at, (true, true, false))
   else
      Lx = L
      Ly = ceil(Int, L / √3)
      at = at * (Lx, Ly, 1)
   end
   return at
end

function edge_geom(; a0 = 1.0, L = 7, b = 1.0, ν = 0.25,
                     xicorr = true, edgevacancy = true,
                     calc = nothing )
   if calc != nothing
      a0 = get_a0(calc)
   end
   # create a homogeneous lattice
   at = bulktri(a0=a0, L=L)
   X = positions(at)
   # find center-atom
   F = defm(at)
   x0 = JVec( [0.5 * diag(F)[1:2]; 0.0] )
   _ , I0 = findmin( [norm(x - x0) for x in X] )
   # dislocation core
   tcore = a0 * JVec([0.5, √3/6, 0.0])
   xc = X[I0] + tcore
   # shift configuration to move core to 0
   X = [x - xc for x in positions(at)]
   set_positions!(at, X)
   # remove the center-atom
   if edgevacancy
      deleteat!(at, I0)
   end
   # apply dislocation FF predictor
   edge_predictor!(at; b=b, xicorr=xicorr, ν=ν)
   if edgevacancy
      X = [x + tcore for x in positions(at)]
      set_positions!(at, X)
   end
   # TODO: compute ν properly
   return at
end

function edge_cluster(R::Number; kwargs...)
   at = edge_geom(;L = ceil(Int, R) * 2 + 3, kwargs...)
   X = positions(at)
   X = X[find(norm.(X) .<= R)]
   cl = ASEAtoms("Al$(length(X))")
   set_positions!(cl, X)
   set_defm!(cl, defm(at))
   set_pbc!(cl, (false, false, false))
   return cl
end

function cluster(R::Number; a0=1.0)
   at = bulktri(a0=a0, L = ceil(Int, R) * 2 + 3)
   X = positions(at)
   xc = mean(X)
   _, I0 = findmin( [norm(x - xc) for x in X] )
   x0 = X[I0]
   X = [x - x0 for x in X]
   X = X[find(norm.(X) .<= R)]
   _, I0 = findmin( norm.(X) )
   X[I0], X[1] = X[1], X[I0]
   cl = ASEAtoms("Al$(length(X))")
   set_positions!(cl, X)
   set_defm!(cl, defm(at))
   set_pbc!(cl, (false, false, false))
   return cl
end

"""
standard isotropic CLE edge dislocation solution
"""
function ulin_edge_isotropic(X, b, ν)
    x, y = X[1,:], X[2,:]
    r² = x.^2 + y.^2
    ux = b/(2*π) * ( angle(x + im*y) + (x .* y) ./ (2*(1-ν) * r²) )
    uy = -b/(2*π) * ( (1-2*ν)/(4*(1-ν)) * log(r²) + - 2 * y.^2 ./ (4*(1-ν) * r²) )
    return [ux'; uy']
end

"""
lattice corrector to CLE edge solution; cf EOS paper
"""
function xi_solver(Y::Vector, b; TOL = 1e-10, maxnit = 5)
   ξ1(x::Real, y::Real, b) = x - b * angle(x + im * y) / (2*π)
   dξ1(x::Real, y::Real, b) = 1 + b * y / (x^2 + y^2) / (2*π)
    y = Y[2]
    x = y
    for n = 1:maxnit
        f = ξ1(x, y, b) - Y[1]
        if abs(f) <= TOL; break; end
        x = x - f / dξ1(x, y, b)
    end
    if abs(ξ1(x, y, b) - Y[1]) > TOL
        warn("newton solver did not converge at Y = $Y; returning input")
        return Y
    end
    return [x, y]
end

"""
EOSShapeev edge dislocation solution
"""
function ulin_edge_eos(X, b, ν)
    Xmod = zeros(2, size(X, 2))
    for n = 1:size(X,2)
        Xmod[:, n] = xi_solver(X[1:2,n], b)
    end
    return ulin_edge_isotropic(Xmod, b, ν)
end

function edge_predictor!(at::AbstractAtoms; b = 1.0, xicorr = true, ν = 0.25)
   X = positions(at) |> mat
   if xicorr
      X[1:2,:] += ulin_edge_eos(X, b, ν)
   else
      X[1:2,:] += ulin_edge_isotropic(X, b, ν)
   end
   set_positions!(at, X)
   return at
end

function get_a0(calc::AbstractCalculator)
   mystress(a) = trace(stress(calc, bulktri(a0 = a)))
   a0, a1 = 1.0, 1.02
   s0, s1 = mystress(a0), mystress(a1)
   while abs(a1 - a0) > 1e-8 && abs(s1) > 1e-8
      anew = (a0 * s1 - a1 * s0) / (s1 - s0)
      a0, a1 = a1, anew
      s0, s1 = s1, mystress(a1)
   end
   return a1
end

end


text:
function has_opf_tables(db::SQLite.DB)
    tables = ["instances", "decompositions", "merges", "combinations"]
    return all(map(in(SQLite.tables(db)[:name]), tables))
end

function count_missing_columns(db::SQLite.DB, table, columns)
    query = """
    SELECT $(join(columns, ',')) FROM $table
    WHERE $(join(columns, " IS NULL OR ")) IS NULL
    """
    results = DBInterface.execute(db, query) |> DataFrame
    states = Dict(zip(names(results), map(x -> sum(ismissing.(x)), eachcol(results))))
    states["total"] = nrow(results)
    return states
end

function table_count(db::SQLite.DB, table)
    query = "SELECT COUNT(*) FROM $table"
    results = DBInterface.execute(db, query) |> DataFrame
    return results[1, Symbol("COUNT(*)")][1]
end

nb_instances(db::SQLite.DB) = table_count(db, "instances")
nb_decompositions(db::SQLite.DB) = table_count(db, "decompositions")

function get_table_ids(db::SQLite.DB, table)
    query = "SELECT id FROM $table"
    results = DBInterface.execute(db, query) |> DataFrame
    return results[!, :id]
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
    struct CXXScopeSpec <: Any
Hold a pointer to a `clang::CXXScopeSpec` object.
"""
struct CXXScopeSpec
    ptr::CXCXXScopeSpec
end

Base.unsafe_convert(::Type{CXCXXScopeSpec}, x::CXXScopeSpec) = x.ptr
Base.cconvert(::Type{CXCXXScopeSpec}, x::CXXScopeSpec) = x


text:
#------------------------------
# functions to define and manipulate ODE-VAE model
#------------------------------

#------------------------------
# ODE systems
#------------------------------

# Linear
function linear_2d_system(du,u,p,t)
    a11, a12, a21, a22 = p
    z1,z2 = u
    du[1] = dz1 = a11 * z1 + a12 * z2
    du[2] = dz2 = a21 * z1 + a22 * z2
end

# Lotka-Volterra
function lotka_volterra(du,u,p,t)
    α, β, γ, δ = p
    z1,z2 = u
    du[1] = dz1 = α*z1 - β*z1*z2 #prey; α: growth, β: meeting rate of prey&predator=predation rate
    du[2] = dz2 = γ*z1*z2 - δ*z2 #predator; γ: growth of predator, δ: loss rate of predator
end

#------------------------------
# define and initialize model
#------------------------------

mutable struct odevae
    p
    q
    zdim
    paramNN
    encoder
    encodedμ
    encodedlogσ
    decoder
    decodedμ
    decodedlogσ
    ODEprob
end

downscaled_glorot_uniform(dims...) = (rand(Float32, dims...) .- 0.5f0) .* sqrt(0.01f0/sum(dims)) # smaller weights initialisation

# initialise model
function init_vae(p, q, zdim, nODEparams, ODEprob; seed::Int=1234)
    myinit = ODEprob.f.f == lotka_volterra ? downscaled_glorot_uniform : Flux.glorot_uniform
    shift(arg) = ODEprob.f.f == lotka_volterra ? tanh.(arg) .+ 1 : sigmoid(arg).-0.5
    # seed
    Random.seed!(seed)
    # parameter network
    paramNN = Chain(Dense(q,q,tanh, init=myinit),
                Dense(q, nODEparams, arg ->(shift(arg)), init=myinit),
                #Dense(q, nODEparams, arg -> (2 ./(1 .+ exp.(-arg))) .- 1.0),
                Flux.Diagonal(nODEparams))
    #   VAE encoder
    Dz, Dh = zdim, p
    encoder, encodedμ, encodedlogσ = Dense(p, Dh, arg ->(tanh.(arg) .+ 1)), Dense(Dh, Dz), Chain(Dense(Dh, Dz, arg -> -relu(arg)), Flux.Diagonal(Dz))
    # VAE decoder
    decoder, decodedμ, decodedlogσ = Dense(Dz, Dh, tanh), Dense(Dh, p), Dense(Dh, p)

    model = odevae(p, q, zdim, paramNN, encoder, encodedμ, encodedlogσ, decoder, decodedμ, decodedlogσ, ODEprob)
    return model
end

#------------------------------
# define model loss
#------------------------------

latentz(μ, logσ) = μ .+ sqrt.(exp.(logσ)) .* randn(Float32,size(μ)...) # sample latent z,

kl_q_p(μ, logσ) = 0.5 .* sum(exp.(logσ) + μ.^2 .- 1 .- (logσ),dims=1)

logp_x_z(m::odevae, x, z) = sum(logpdf.(Normal.(m.decodedμ(m.decoder(z)), sqrt.(exp.(m.decodedlogσ(m.decoder(z))))), x),dims=1) # get reconstruction error

sqnorm(x) = sum(abs2, x)
reg(m::odevae) = sum(sqnorm, Flux.params(m.decoder,m.decodedμ,m.decodedlogσ)) # regularisation term in loss

getparams(m::odevae) = Flux.params(m.encoder, m.encodedμ, m.encodedlogσ, m.decoder, m.decodedμ, m.decodedlogσ, m.paramNN) # get parameters of VAE model

function loss(X, Y, t, m::odevae)
    latentμ, latentlogσ = m.encodedμ(m.encoder(X)), m.encodedlogσ(m.encoder(X))
    learnedparams = m.paramNN(Y)
    curparams = [learnedparams[1], m.ODEprob.p[2], m.ODEprob.p[3], learnedparams[2]]
    curts = Int.(floor.(t .*(1 ./dt) .+1))
    smoothμ = Array(solve(m.ODEprob, Tsit5(), u0 = [latentμ[1,1], latentμ[2,1]], p=curparams, saveat=dt))[:,curts]
    combinedμ = hcat(latentμ[:,1],smoothμ)
    combinedz = latentz.(combinedμ, latentlogσ)
    ELBO = 1.0 .* logp_x_z(m, X, combinedz) .- 0.5 .* kl_q_p(combinedμ, latentlogσ)
    lossval = sum(-ELBO) + 0.01*reg(m) + 1.0*sum((smoothμ .- latentμ[:,2:end]).^2)
    return lossval
end

function loss_wrapper(m::odevae)
    return function (X, Y, t) loss(X, Y, t, m) end
end

# to experiment:
# xs, tvals, group1, group2 = generate_xs(n, p, true_u0, sol_group1, sol_group2);
# trainingdata = zip(xs, x_baseline, tvals);
# X, Y, t = first(trainingdata)


text:
using OMEinsum
using OMEinsum.ContractionOrder
using OMEinsum.ContractionOrder: analyze_contraction, contract_pair!, evaluate_costs, contract_tree!, log2sumexp2
using TropicalNumbers

using Test, Random
@testset "analyze contraction" begin
    incidence_list = IncidenceList(Dict('A' => ['a', 'b', 'k', 'o', 'f'], 'B'=>['a', 'c', 'd', 'm', 'f'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']), openedges=['c', 'f', 'o'])
    info = analyze_contraction(incidence_list, 'A', 'B')
    @test Set(info.l1) == Set(['k'])
    @test Set(info.l2) == Set(['m'])
    @test Set(info.l12) == Set(['a'])
    @test Set(info.l01) == Set(['b','o'])
    @test Set(info.l02) == Set(['c', 'd'])
    @test Set(info.l012) == Set(['f'])
end

@testset "tree greedy" begin
    Random.seed!(2)
    incidence_list = IncidenceList(Dict('A' => ['a', 'b'], 'B'=>['a', 'c', 'd'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']))
    log2_edge_sizes = Dict([c=>i for (i,c) in enumerate(['a', 'b', 'c', 'd', 'e', 'f'])]...)
    edge_sizes = Dict([c=>(1<<i) for (i,c) in enumerate(['a', 'b', 'c', 'd', 'e', 'f'])]...)
    il = copy(incidence_list)
    contract_pair!(il, 'A', 'B', log2_edge_sizes)
    target = IncidenceList(Dict('A' => ['b', 'c', 'd'], 'C'=>['b', 'c', 'e', 'f'], 'D'=>['e'], 'E'=>['d', 'f']))
    @test il.v2e == target.v2e
    @test length(target.e2v) == length(il.e2v)
    for (k,v) in il.e2v
        @test sort(target.e2v[k]) == sort(v)
    end
    costs = evaluate_costs(MinSpaceOut(), incidence_list, log2_edge_sizes)
    @test costs == Dict(('A', 'B')=>9, ('A', 'C')=>15, ('B','C')=>18, ('B','E')=>10, ('C','D')=>11, ('C', 'E')=>14)
    tree, log2_tcs, log2_scs = tree_greedy(incidence_list, log2_edge_sizes)
    tcs_, scs_ = [], []
    contract_tree!(copy(incidence_list), tree, log2_edge_sizes, tcs_, scs_)
    @test all((log2sumexp2(tcs_), maximum(scs_)) .<= (log2(exp2(10)+exp2(16)+exp2(15)+exp2(9)), 11))
    vertices = ['A', 'B', 'C', 'D', 'E']
    optcode1 = parse_eincode(incidence_list, tree, vertices=vertices)
    @test optcode1 isa OMEinsum.NestedEinsum
    tree2 = OMEinsum.parse_tree(optcode1, vertices)
    @test tree2 == tree

    eincode = ein"ab,acd,bcef,e,df->"
    size_dict = Dict([c=>(1<<i) for (i,c) in enumerate(['a', 'b', 'c', 'd', 'e', 'f'])]...)
    Random.seed!(2)
    optcode2 = optimize_greedy(eincode, size_dict)
    tc, sc = timespace_complexity(optcode2, edge_sizes)
    @test 16 <= tc <= log2(exp2(10)+exp2(16)+exp2(15)+exp2(9))
    @test sc == 11
    @test optcode1 == optcode2
    eincode3 = ein"(ab,acd),bcef,e,df->"
    Random.seed!(2)
    optcode3 = optimize_greedy(eincode3, size_dict)
    tc, sc = timespace_complexity(optcode3, edge_sizes)
    @test 16 <= tc <= log2(exp2(10)+exp2(16)+exp2(15)+exp2(9)+1e-8)
end

@testset "fullerene" begin
    function fullerene()
        φ = (1+√5)/2
        res = NTuple{3,Float64}[]
        for (x, y, z) in ((0.0, 1.0, 3φ), (1.0, 2 + φ, 2φ), (φ, 2.0, 2φ + 1.0))
            for (α, β, γ) in ((x,y,z), (y,z,x), (z,x,y))
                for loc in ((α,β,γ), (α,β,-γ), (α,-β,γ), (α,-β,-γ), (-α,β,γ), (-α,β,-γ), (-α,-β,γ), (-α,-β,-γ))
                    if loc ∉ res
                        push!(res, loc)
                    end
                end
            end
        end
        return res
    end

    c60_xy = fullerene()
    c60_edges = [(i,j) for (i,(i2,j2,k2)) in enumerate(c60_xy), (j,(i1,j1,k1)) in enumerate(c60_xy) if i<j && (i2-i1)^2+(j2-j1)^2+(k2-k1)^2 < 5.0]
    code = EinCode((c60_edges..., [(i,) for i=1:60]...), ())
    size_dict = Dict([i=>2 for i in 1:60])
    log2_edge_sizes = Dict([i=>1 for i in 1:60])
    edge_sizes = Dict([i=>2 for i in 1:60])
    tc, sc = timespace_complexity(code, edge_sizes)
    @test tc == 60
    @test sc == 0
    optcode = optimize_greedy(code, size_dict)
    tc2, sc2 = timespace_complexity(optcode, edge_sizes)
    @test sc2 == 10
    xs = vcat([TropicalF64.([-1 1; 1 -1]) for i=1:90], [TropicalF64.([0, 0]) for i=1:60])
    @test OMEinsum.flatten(optcode) == code
    @test OMEinsum.flatten(code) == code
    @test optcode(xs...)[].n == 66
end

@testset "regression test" begin
    code = ein"i->"
    optcode = optimize_greedy(code, Dict('i'=>3))
    @test optcode isa NestedEinsum
    x = randn(3)
    @test optcode(x) ≈ code(x)

    code = ein"i,j->"
    optcode = optimize_greedy(code, Dict('i'=>3, 'j'=>3))
    @test optcode isa NestedEinsum
    x = randn(3)
    y = randn(3)
    @test optcode(x, y) ≈ code(x, y)

    code = ein"ij,jk,kl->ijl"
    optcode = optimize_greedy(code, Dict('i'=>3, 'j'=>3, 'k'=>3, 'l'=>3))
    @test optcode isa NestedEinsum
    a, b, c = [rand(3,3) for i=1:3]
    @test optcode(a, b, c) ≈ code(a, b, c)
end


text:
## namedarraytypes.jl.
## (c) 2013 David A. van Leeuwen

## Julia type that implements a drop-in replacement of Array with named dimensions.

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## DT is a tuple of Dicts, characterized by the types of the keys.
## This way NamedArray is dependent on the dictionary type of each dimensions.
## The inner constructor checks for consistency, the values must all be 1:d
if !isdefined(:NamedArray)

type NamedArray{T,N,AT,DT} <: AbstractArray{T,N}
    array::AT
    dicts::DT
    dimnames::NTuple{N}
    function NamedArray(array::AbstractArray{T,N}, dicts::NTuple{N,OrderedDict}, dimnames::NTuple{N})
        size(array) == map(length, dicts) || error("Inconsistent dictionary sizes")
#        for (d,dict) in zip(size(array),dicts)
#            Set(values(dict)) == Set(1:d) || error("Inconsistent values in dict")
#        end
        new(array, dicts, dimnames)
    end
end

## a type that negates any index
immutable Not{T}
    index::T
end

typealias NamedVector{T} NamedArray{T,1}
typealias NamedMatrix{T} NamedArray{T,2}
@compat typealias NamedVecOrMat{T} Union{NamedVector{T},NamedMatrix{T}}
@compat typealias ArrayOrNamed{T,N} Union{Array{T,N}, NamedArray{T,N}}

end


text:
# Main GaussianProcess type

mutable struct GPA{X<:AbstractMatrix,Y<:AbstractVector{<:Real},M<:Mean,K<:Kernel,L<:Likelihood,
                    CS<:CovarianceStrategy, D<:KernelData} <: GPBase
    # Observation data
    "Input observations"
    x::X
    "Output observations"
    y::Y

    # Model
    "Mean object"
    mean::M
    "Kernel object"
    kernel::K
    "Likelihood"
    lik::L
    "Strategy for computing or approximating covariance matrices"
    covstrat::CS

    # Auxiliary data
    "Dimension of inputs"
    dim::Int
    "Number of observations"
    nobs::Int
    "Auxiliary observation data (to speed up calculations)"
    data::D
    "Latent (whitened) variables - N(0,1)"
    v::Vector{Float64}
    "Mean values"
    μ::Vector{Float64}
    "`(k + exp(2*obsNoise))`"
    cK::PDMat{Float64,Matrix{Float64}}
    "Log-likelihood"
    ll::Float64
    "Gradient of log-likelihood"
    dll::Vector{Float64}
    "Log-target (marginal log-likelihood + log priors)"
    target::Float64
    "Gradient of log-target (gradient of marginal log-likelihood + gradient of log priors)"
    dtarget::Vector{Float64}

    function GPA{X,Y,M,K,L,CS,D}(x::X, y::Y, mean::M, kernel::K, lik::L, covstrat::CS, data::D) where {X,Y,M,K,L,CS,D}
        dim, nobs = size(x)
        length(y) == nobs || throw(ArgumentError("Input and output observations must have consistent dimensions."))
        gp = new{X,Y,M,K,L,CS,D}(x, y, mean, kernel, lik, covstrat, dim, nobs,
                                 data, zeros(nobs))
        initialise_target!(gp)
    end
end
@deprecate GPMC GPA

function GPA(x::AbstractMatrix, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel, lik::Likelihood, covstrat::CovarianceStrategy)
    data = KernelData(kernel, x, x, covstrat)
    return GPA{typeof(x),typeof(y),typeof(mean),typeof(kernel),typeof(lik),typeof(covstrat),typeof(data)}(
                x, y, mean, kernel, lik, covstrat, data)
end

"""
    GPA(x, y, mean, kernel, lik)

Fit a Gaussian process to a set of training points. The Gaussian process with
non-Gaussian observations is defined in terms of its user-defined likelihood function,
mean and covaiance (kernel) functions.

The non-Gaussian likelihood is handled by an approximate method (e.g. Monte Carlo). The latent function
values are represented by centered (whitened) variables ``f(x) = m(x) + Lv`` where
``v ∼ N(0, I)`` and ``LLᵀ = K_θ``.

# Arguments:
- `x::AbstractVecOrMat{Float64}`: Input observations
- `y::AbstractVector{<:Real}`: Output observations
- `mean::Mean`: Mean function
- `kernel::Kernel`: Covariance function
- `lik::Likelihood`: Likelihood function
"""
function GPA(x::AbstractMatrix, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel, lik::Likelihood)
    covstrat = FullCovariance()
    return GPA(x, y, mean, kernel, lik, covstrat)
end

GPA(x::AbstractVector, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel, lik::Likelihood) =
    GPA(x', y, mean, kernel, lik)

"""
    GP(x, y, mean::Mean, kernel::Kernel, lik::Likelihood)

Fit a Gaussian process that is defined by its `mean`, its `kernel`, and its likelihood
function `lik` to a set of training points `x` and `y`.

See also: [`GPA`](@ref)
"""
GP(x::AbstractVecOrMat{Float64}, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel,
   lik::Likelihood) = GPA(x, y, mean, kernel, lik)

"""
    initialise_ll!(gp::GPA)

Initialise the log-likelihood of Gaussian process `gp`.
"""
function initialise_ll!(gp::GPA)
    # log p(Y|v,θ)
    gp.μ = mean(gp.mean,gp.x)
    Σ = cov(gp.kernel, gp.x, gp.x, gp.data)
    gp.cK = PDMat(Σ + 1e-6*I)
    F = unwhiten(gp.cK,gp.v) + gp.μ
    gp.ll = sum(log_dens(gp.lik,F,gp.y)) #Log-likelihood
    gp
end

"""
    update_cK!(gp::GPA)

Update the covariance matrix and its Cholesky decomposition of Gaussian process `gp`.
"""
function update_cK!(gp::GPA)
    old_cK = gp.cK
    Σbuffer = old_cK.mat
    cov!(Σbuffer, gp.kernel, gp.x, gp.x, gp.data)
    for i in 1:gp.nobs
        Σbuffer[i,i] += 1e-6 # no logNoise for GPA
    end
    chol_buffer = old_cK.chol.factors
    copyto!(chol_buffer, Σbuffer)
    chol = cholesky!(Symmetric(chol_buffer))
    gp.cK = PDMat(Σbuffer, chol)
    gp
end

# modification of initialise_ll! that reuses existing matrices to avoid
# unnecessary memory allocations, which speeds things up significantly
function update_ll!(gp::GPA; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    if kern
        # only need to update the covariance matrix
        # if the covariance parameters have changed
        update_cK!(gp)
    end
    gp.μ = mean(gp.mean,gp.x)
    F = unwhiten(gp.cK,gp.v) + gp.μ
    gp.ll = sum(log_dens(gp.lik,F,gp.y)) #Log-likelihood
    gp
end

function get_L_bar!(L_bar::AbstractMatrix, dl_df::AbstractVector, v::AbstractVector, cK::PDMat)
    fill!(L_bar, 0.0)
    BLAS.ger!(1.0, dl_df, v, L_bar)
    tril!(L_bar)
    # ToDo:
    # the following two steps allocates memory
    # and are fickle, reaching into the internal
    # implementation of the cholesky decomposition
    L = cK.chol.L.data
    tril!(L)
    #
    chol_unblocked_rev!(L, L_bar)
    return L_bar
end

struct FullCovMCMCPrecompute <: AbstractGradientPrecompute
    L_bar::Matrix{Float64}
    dl_df::Vector{Float64}
    f::Vector{Float64}
end
function FullCovMCMCPrecompute(nobs::Int)
    buffer1 = Matrix{Float64}(undef, nobs, nobs)
    buffer2 = Vector{Float64}(undef, nobs)
    buffer3 = Vector{Float64}(undef, nobs)
    return FullCovMCMCPrecompute(buffer1, buffer2, buffer3)
end
init_precompute(gp::GPA) = FullCovMCMCPrecompute(gp.nobs)

function precompute!(precomp::FullCovMCMCPrecompute, gp::GPBase)
    f = unwhiten(gp.cK, gp.v)  + gp.μ
    dl_df = dlog_dens_df(gp.lik, f, gp.y)
    precomp.dl_df[:] = dl_df
    precomp.f[:] = f
end
function dll_kern!(dll::AbstractVector, gp::GPBase, precomp::FullCovMCMCPrecompute, covstrat::CovarianceStrategy)
    L_bar = precomp.L_bar
    get_L_bar!(L_bar, precomp.dl_df, gp.v, gp.cK)
    nobs = gp.nobs
    @inbounds for i in 1:nobs
        L_bar[i,i] *= 2
    end
    # in GPA, L_bar plays the role of ααinvcKI
    return dmll_kern!(dll, gp.kernel, gp.x, gp.data, L_bar, covstrat)
end
function dll_mean!(dll::AbstractVector, gp::GPBase, precomp::FullCovMCMCPrecompute)
    dmll_mean!(dll, gp.mean, gp.x, precomp.dl_df)
end

"""
     update_dll!(gp::GPA, ...)

Update the gradient of the log-likelihood of Gaussian process `gp`.
"""
function update_dll!(gp::GPA, precomp::AbstractGradientPrecompute;
    process::Bool=true, # include gradient components for the process itself
    lik::Bool=true,  # include gradient components for the likelihood parameters
    domean::Bool=true, # include gradient components for the mean parameters
    kern::Bool=true, # include gradient components for the spatial kernel parameters
    )

    n_lik_params = num_params(gp.lik)
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)

    gp.dll = Array{Float64}(undef, process * gp.nobs + lik * n_lik_params +
                            domean * n_mean_params + kern * n_kern_params)
    precompute!(precomp, gp)

    i = 1
    if process
        mul!(view(gp.dll, i:i+gp.nobs-1), gp.cK.chol.U, precomp.dl_df)
        i += gp.nobs
    end
    if lik && n_lik_params > 0
        Lgrads = dlog_dens_dθ(gp.lik, precomp.f, gp.y)
        for j in 1:n_lik_params
            gp.dll[i] = sum(Lgrads[:,j])
            i += 1
        end
    end
    # if domean && n_mean_params > 0
        # Mgrads = grad_stack(gp.mean, gp.x)
        # for j in 1:n_mean_params
            # gp.dll[i] = dot(precomp.dl_df,Mgrads[:,j])
            # i += 1
        # end
    # end
    if domean && n_mean_params>0
        dll_m = @view(gp.dll[i:i+n_mean_params-1])
        dll_mean!(dll_m, gp, precomp)
        i += n_mean_params
    end
    if kern
        dll_k = @view(gp.dll[i:end])
        dll_kern!(dll_k, gp, precomp, gp.covstrat)
    end

    gp
end

function update_ll_and_dll!(gp::GPA, precomp::AbstractGradientPrecompute; kwargs...)
    update_ll!(gp; kwargs...)
    update_dll!(gp, precomp; kwargs...)
end

"""
    initialise_target!(gp::GPA)

Initialise the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp`.
"""
function initialise_target!(gp::GPA)
    initialise_ll!(gp)
    gp.target = gp.ll - (sum(abs2, gp.v) + log2π * gp.nobs) / 2 +
        prior_logpdf(gp.lik) + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

"""
    update_target!(gp::GPA, ...)

Update the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp`.
"""
function update_target!(gp::GPA; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    update_ll!(gp; process=process, lik=lik, domean=domean, kern=kern)
    gp.target = gp.ll - (sum(abs2, gp.v) + log2π * gp.nobs) / 2 +
        prior_logpdf(gp.lik) + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

function update_dtarget!(gp::GPA, precomp::AbstractGradientPrecompute; kwargs...)
    update_dll!(gp, precomp; kwargs...)
    gp.dtarget = gp.dll + prior_gradlogpdf(gp; kwargs...)
    gp
end

"""
    update_target_and_dtarget!(gp::GPA, ...)

Update the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp` and its derivative.
"""
function update_target_and_dtarget!(gp::GPA, precomp::AbstractGradientPrecompute; kwargs...)
    update_target!(gp; kwargs...)
    update_dtarget!(gp, precomp; kwargs...)
end

function update_target_and_dtarget!(gp::GPA; kwargs...)
    precomp = init_precompute(gp)
    update_target_and_dtarget!(gp, precomp; kwargs...)
end

predict_full(gp::GPA, xpred::AbstractMatrix) = predictMVN(xpred, gp.x, gp.y, gp.kernel, gp.mean, gp.cK\unwhiten(gp.cK,gp.v), gp.covstrat, gp.cK)
"""
    predict_y(gp::GPA, x::Union{Vector{Float64},Matrix{Float64}}[; full_cov::Bool=false])

Return the predictive mean and variance of Gaussian Process `gp` at specfic points which
are given as columns of matrix `x`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""

function predict_y(gp::GPA, x::AbstractMatrix; full_cov::Bool=false)
    μ, σ2 = predict_f(gp, x; full_cov=full_cov)
    return predict_obs(gp.lik, μ, σ2)
end

appendlikbounds!(lb, ub, gp::GPA, bounds) = appendbounds!(lb, ub, num_params(gp.lik), bounds)

#—————————————————————————————————————————————————————–
# Function for sampling from the prior of the GP object hyperparameters.

function sample_params(gp::GPA; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    samples = Float64[]
    if lik && num_params(gp.lik)>0
        like_priors = get_priors(gp.lik)
        @assert !isempty(like_priors) "prior distributions of likelihood hyperparameters should be set"
        like_sample = rand(Product(like_priors))
        append!(samples, like_sample)
    end
    if domean && num_params(gp.mean)>0
        mean_priors = get_priors(gp.mean)
        @assert !isempty(mean_priors) "prior distributions of mean hyperparameters should be set"
        mean_sample = rand(Product(mean_priors))
        append!(samples, mean_sample)
    end
    if kern && num_params(gp.kernel)>0
        kernel_priors = get_priors(gp.kernel)
        @assert !isempty(kernel_priors) "prior distributions of kernel hyperparameters should be set"
        kernel_sample = rand(Product(kernel_priors))
        append!(samples, kernel_sample)
    end
    return samples
end

function sample_params(rng::AbstractRNG, gp::GPA; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    samples = Float64[]
    if lik && num_params(gp.lik)>0
        like_priors = get_priors(gp.lik)
        @assert !isempty(like_priors) "prior distributions of likelihood hyperparameters should be set"
        like_sample = rand(rng, Product(like_priors))
        append!(samples, like_sample)
    end
    if domean && num_params(gp.mean)>0
        mean_priors = get_priors(gp.mean)
        @assert !isempty(mean_priors) "prior distributions of mean hyperparameters should be set"
        mean_sample = rand(rng, Product(mean_priors))
        append!(samples, mean_sample)
    end
    if kern && num_params(gp.kernel)>0
        kernel_priors = get_priors(gp.kernel)
        @assert !isempty(kernel_priors) "prior distributions of kernel hyperparameters should be set"
        kernel_sample = rand(rng, Product(kernel_priors))
        append!(samples, kernel_sample)
    end
    return samples
end

function get_params(gp::GPA; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    params = Float64[]
    append!(params, gp.v)
    if lik  && num_params(gp.lik)>0
        append!(params, get_params(gp.lik))
    end
    if domean && num_params(gp.mean)>0
        append!(params, get_params(gp.mean))
    end
    if kern
        append!(params, get_params(gp.kernel))
    end
    return params
end

function num_params(gp::GPA; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    n = length(gp.v)
    lik && (n += num_params(gp.lik))
    domean && (n += num_params(gp.mean))
    kern && (n += num_params(gp.kernel))
    n
end

function set_params!(gp::GPA, hyp::AbstractVector; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    n_lik_params = num_params(gp.lik)
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)

    i = 1
    if process
        gp.v = hyp[1:gp.nobs]
        i += gp.nobs
    end
    if lik  && n_lik_params>0
        set_params!(gp.lik, hyp[i:i+n_lik_params-1]);
        i += n_lik_params
    end
    if domean && n_mean_params>0
        set_params!(gp.mean, hyp[i:i+n_mean_params-1])
        i += n_mean_params
    end
    if kern
        set_params!(gp.kernel, hyp[i:i+n_kern_params-1])
        i += n_kern_params
    end
end

function prior_gradlogpdf(gp::GPA; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    if process
        grad = -gp.v
    else
        grad = Float64[]
    end
    if lik
        append!(grad, prior_gradlogpdf(gp.lik))
    end
    if domean
        append!(grad, prior_gradlogpdf(gp.mean))
    end
    if kern
        append!(grad, prior_gradlogpdf(gp.kernel))
    end
    return grad
end

function Base.show(io::IO, gp::GPA)
    println(io, "GP Approximate object:")
    println(io, "  Dim = ", gp.dim)
    println(io, "  Number of observations = ", gp.nobs)
    println(io, "  Mean function:")
    show(io, gp.mean, 2)
    println(io, "\n  Kernel:")
    show(io, gp.kernel, 2)
    println(io, "\n  Likelihood:")
    show(io, gp.lik, 2)
    if gp.nobs == 0
        println("\n  No observation data")
    else
        println(io, "\n  Input observations = ")
        show(io, gp.x)
        print(io,"\n  Output observations = ")
        show(io, gp.y)
        print(io,"\n  Log-posterior = ")
        show(io, round(gp.target; digits = 3))
    end
end


text:
using ReservoirComputing

#padding
test_array = [1, 2, 3, 4, 5, 6, 7, 8, 9]
standard_array = zeros(length(test_array), 1)
extension = [0, 0, 0]
padded_array = zeros(length(test_array)+1, 1)
extended_array = zeros(length(test_array)+length(extension), 1)
padded_extended_array = zeros(length(test_array)+length(extension)+1, 1)
padding = 10.0

#testing non linear algos
nla_array = ReservoirComputing.nla(NLADefault(), test_array)
@test nla_array == test_array

nla_array = ReservoirComputing.nla(NLAT1(), test_array)
@test nla_array == [1, 2, 9, 4, 25, 6, 49, 8, 81]

nla_array = ReservoirComputing.nla(NLAT2(), test_array)
@test nla_array == [1, 2, 2, 4, 12, 6, 30, 8, 9]

nla_array = ReservoirComputing.nla(NLAT3(), test_array)
@test nla_array == [1, 2, 8, 4, 24, 6, 48, 8, 9]

#testing padding and extension
states_type = StandardStates()
standard_array = states_type(NLADefault(), test_array, extension)
@test standard_array == test_array

states_type = PaddedStates(padding = padding)
padded_array = states_type(NLADefault(), test_array, extension)
@test padded_array == reshape(vcat(padding, test_array), length(test_array)+1, 1)

states_type = PaddedStates(padding)
padded_array = states_type(NLADefault(), test_array, extension)
@test padded_array == reshape(vcat(padding, test_array), length(test_array)+1, 1)

states_type = PaddedExtendedStates(padding = padding)
padded_extended_array = states_type(NLADefault(), test_array, extension)
@test padded_extended_array == reshape(vcat(padding, extension, test_array), length(test_array)+length(extension)+1, 1)

states_type = PaddedExtendedStates(padding)
padded_extended_array = states_type(NLADefault(), test_array, extension)
@test padded_extended_array == reshape(vcat(padding, extension, test_array), length(test_array)+length(extension)+1, 1)

states_type = ExtendedStates()
extended_array = states_type(NLADefault(), test_array, extension)
@test extended_array == vcat(extension, test_array)
#reshape(vcat(extension, test_array), length(test_array)+length(extension), 1)


text:
function RBMs.initialize!(
    rbm::DenseConvRBM, data_dense::AbstractArray, data_conv::AbstractArray;
    ϵ::Real = 1e-6, λ::Real = 0.1
)
    RBMs.initialize!(rbm.visible_dense, data_dense; ϵ)
    RBMs.initialize!(rbm.visible_conv, data_conv; ϵ)
    RBMs.initialize!(hidden(rbm))
    RBMs.initialize_w!(parts(rbm).dense, data_dense; λ = λ/2)
    RBMs.initialize_w!(parts(rbm).conv, data_conv; λ = λ/2)
    RBMs.zerosum!(rbm)
    return rbm
end


text:
# only for internal use to avoid allocations
mutable struct InternalRay{T<:Real}
    x::T
    y::T
    z::T
    cx::T
    cy::T
    cz::T
    _e1::T
end

"""
A Ray described by its (x, y, z) coordinates and the
direction cosines (cx, cy, cz). The constructor ensures, that the
relation cx^2 + cy^2 + cz^2 == 1 holds.
"""
struct Ray{T<:Real}
    x::T
    y::T
    z::T
    cx::T
    cy::T
    cz::T
    function Ray(a, b, c, d, e)
        (x, y, z, l, m) = promote(a, b, c, d, e)
        return new{typeof(x)}(x, y, z, l, m, sqrt(1 - l^2 - m^2))
    end
end

function to_internal(t, ray::Ray)
    InternalRay(
        convert(t, ray.x),
        convert(t, ray.y),
        convert(t, ray.z),
        convert(t, ray.cx),
        convert(t, ray.cy),
        convert(t, ray.cz),
        zero(t),
    )
end

"""
    ray_from_NA(y, NA)

Create a Ray in the tangential plane with ray height y and direction cosine in y given by NA
"""
function ray_from_NA(y, NA)
    return Ray(zero(y), y, zero(y), zero(y), NA)
end

# extends ≈ operator
function Base.isapprox(ray1::Ray, ray2::Ray; kwargs...)
    isapprox(ray1.x, ray2.x; kwargs...) &&
        isapprox(ray1.y, ray2.y; kwargs...) &&
        isapprox(ray1.z, ray2.z; kwargs...) &&
        isapprox(ray1.cx, ray2.cx; kwargs...) &&
        isapprox(ray1.cy, ray2.cy; kwargs...) &&
        isapprox(ray1.cz, ray2.cz; kwargs...)
end


text:
using ContinuousTransformations
using Test

using Distributions: logpdf, Normal, LogNormal, MvNormal, MvLogNormal
import ForwardDiff: derivative, jacobian
using Parameters
using LinearAlgebra, Statistics, Markdown

# utilities

"Test singleton type."
ContinuousTransformations.@define_singleton TestSingleton <: Real

@testset "define singleton" begin
    @test TestSingleton <: Real
    @test isa(TESTSINGLETON, TestSingleton)
    @test repr(@doc(TESTSINGLETON)) == repr(@doc(TestSingleton)) ==
        repr(md"Test singleton type.")
end

"""
    rand_in(x)

Random number in `x` (eg a finite or infinite interval).
"""
function rand_in(segment::Segment)
    @unpack left, right = segment
    rand()*(right-left) + left
end

rand_in(ray::PositiveRay) = ray.left + randn()^2

rand_in(ray::NegativeRay) = ray.right - randn()^2

rand_in(::RealLine) = randn()

"Elements below the diagonal as a vector (by row)."
lower_to_vec(L) = vcat((L[i, 1:(i-1)] for i in 1:size(L, 1))...)

"Elements below and including the diagonal as a vector (by row)."
lowerdiag_to_vec(L) = vcat((L[i, 1:i] for i in 1:size(L, 1))...)

"""
Reconstruct lower triangular matrix (including diagonal) from a vector of
elemenst (by row).
"""
function vec_to_lowerdiag(l::AbstractVector{T}, n::Int) where T
    A = zeros(T, n, n)
    cumulative_index = 0
    for i in 1:n
        A[i, 1:i] .= l[cumulative_index .+ (1:i)]
        cumulative_index += i
    end
    LowerTriangular(A)
end

@testset "vec lower utilities" begin
    l = Float64.(1:6)
    L = vec_to_lowerdiag(l, 3)
    @test size(L) == (3, 3)
    @test eltype(L) == eltype(l)
    @test L == [1.0 0.0 0.0;
                2.0 3.0 0.0;
                4.0 5.0 6.0]
    @test lowerdiag_to_vec(L) == l
    @test lower_to_vec(L) == [2.0, 4.0, 5.0]
end

# test: intervals

@testset "interval constructors" begin
    @test_throws ArgumentError Segment(NaN, NaN)
    @test_throws ArgumentError Segment(-Inf, Inf)
    @test_throws ArgumentError Segment(2, -1)
    @test_throws ArgumentError Segment(2, 2)
    @test_throws ArgumentError PositiveRay(-Inf)
    @test_throws ArgumentError PositiveRay(NaN)
    @test_throws ArgumentError NegativeRay(-Inf)
    @test_throws ArgumentError NegativeRay(NaN)
end

@testset "interval equality" begin
    @test PositiveRay(1.0) == PositiveRay(1)
    @test NegativeRay(2.0) == NegativeRay(2)
    @test RealLine() == ℝ
    @test isa(Segment(1,2.0), Segment{Float64})
    @test Segment(1, 2.0) == Segment(1, 2) == Segment(1.0, 2.0)
end

@testset "named interval constants" begin
    @test ℝ == RealLine()
    @test ℝ⁺ == PositiveRay(0.0)
    @test ℝ⁻ == NegativeRay(0.0)
end

@testset "interval isapprox" begin
    @test ℝ ≈ ℝ
    @test !(ℝ ≈ Segment(1,2))
    @test !(ℝ ≈ ℝ⁺)
    @test !(ℝ ≈ ℝ⁻)
    @test PositiveRay(1) ≈ PositiveRay(1+eps())
    @test NegativeRay(-2) ≈ NegativeRay(-2+eps())
    @test Segment(1, 2) ≈ Segment(1+eps(), 2+eps())
end

@testset "intervals ∈, width, extrema, finiteness" begin
    seg = Segment(1.0, 2.0)
    posray = PositiveRay(0.0)
    negray = NegativeRay(1.5)
    # numbers in seg
    @test 1.0 ∈ seg
    @test 1.5 ∈ seg
    @test 2.0 ∈ seg
    @test 0.0 ∉ seg
    @test Inf ∉ seg
    @test_throws MethodError "string" ∈ seg
    # methods of seg
    @test width(seg) == 1.0
    @test middle(seg) == 1.5
    # numbers in posray
    @test 1.0 ∈ posray
    @test Inf ∈ posray
    @test 0 ∈ posray
    @test -1 ∉ posray
    @test -Inf ∉ posray
    @test_throws MethodError "string" ∈ posray
    # numbers in negray
    @test -Inf ∈ negray
    @test 0 ∈ negray
    @test 1.5 ∈ negray
    @test 2 ∉ negray
    @test Inf ∉ negray
    @test_throws MethodError "string" ∈ negray
    # numbers in the real line
    @test 1 ∈ ℝ
    @test Inf ∈ ℝ
    @test -Inf ∈ ℝ
    @test_throws MethodError "string" ∈ ℝ
    # finiteness
    @test isfinite(seg) && !isinf(seg)
    @test !isfinite(posray) && isinf(posray)
    @test !isfinite(negray) && isinf(negray)
    @test !isfinite(posray) && isinf(posray)
    @test !isfinite(ℝ) && isinf(ℝ)
end

@testset "intervals intersections" begin
    seg = Segment(1.0, 2.0)
    posray = PositiveRay(0.0)
    negray = NegativeRay(1.5)
    # intersections with ℝ
    @test seg ∩ ℝ == seg
    @test ℝ ∩ seg == seg
    @test posray ∩ ℝ == posray
    @test ℝ ∩ posray == posray
    @test negray ∩ ℝ == negray
    @test ℝ ∩ negray == negray
    @test ℝ ∩ ℝ == ℝ
    # empty intersections
    @test_throws Exception 𝕀∩ seg
    @test_throws Exception posray ∩ ℝ⁻
    @test_throws Exception ℝ⁻ ∩ seg
    # non-empty intersections
    let seg2 = Segment(1.5, 3.0)
        @test seg ∩ seg2 == seg2 ∩ seg == Segment(1.5, 2.0)
    end
    @test seg ∩ posray == posray ∩ seg == seg
    @test seg ∩ negray == negray ∩ seg == Segment(1.0, 1.5)
    @test negray ∩ posray == posray ∩ negray == Segment(0.0, 1.5)
    @test posray ∩ PositiveRay(2) == PositiveRay(2)
    @test posray ∩ PositiveRay(-2) == posray
    @test negray ∩ NegativeRay(-7) == NegativeRay(-7)
    @test negray ∩ NegativeRay(7) == negray
end

@testset "interval printing" begin
    @test repr(ℝ) == "ℝ"
end

# test: univariate transformation

@testset "type hierarchy" begin
    @test ContinuousTransformation <: Function
    @test UnivariateTransformation <: ContinuousTransformation
    @test GroupedTransformation <: ContinuousTransformation
    @test CorrelationCholeskyFactor <: ContinuousTransformation
    @test UnitVector <: ContinuousTransformation
    @test TransformationWrapper <: Function
end

@testset "univariate transformation basics" begin
    @test_throws DomainError Affine(0, 1)
    @test_throws DomainError Affine(-1, 2.0)
    @test Affine(1, 2.0) == Affine(1.0, 2.0)
    a = Affine(1,2)
    @test domain(a) == ℝ
    @test image(a) == ℝ
    @test isincreasing(a)

    @test domain(NEGATION) == ℝ
    @test image(NEGATION) == ℝ
    @test !isincreasing(NEGATION)
    @test inverse(NEGATION) == NEGATION

    @test domain(LOGISTIC) == ℝ
    @test image(LOGISTIC) == Segment(0, 1)
    @test isincreasing(LOGISTIC)
    @test inverse(LOGISTIC) == LOGIT

    @test domain(LOGIT) == Segment(0, 1)
    @test image(LOGIT) == ℝ
    @test isincreasing(LOGIT)
    @test inverse(LOGIT) == LOGISTIC

    @test domain(REALCIRCLE) == ℝ
    @test image(REALCIRCLE) == Segment(-1, 1)
    @test isincreasing(REALCIRCLE)
    @test inverse(REALCIRCLE) == INVREALCIRCLE

    @test domain(INVREALCIRCLE) == Segment(-1, 1)
    @test image(INVREALCIRCLE) == ℝ
    @test isincreasing(INVREALCIRCLE)
    @test inverse(INVREALCIRCLE) == REALCIRCLE

    @test domain(EXP) == ℝ
    @test image(EXP) == ℝ⁺
    @test isincreasing(EXP)
    @test inverse(EXP) == LOG

    @test domain(LOG) == ℝ⁺
    @test image(LOG) == ℝ
    @test isincreasing(LOG)
    @test inverse(LOG) == EXP
end

"""
    test_univariate(f, x; [AD_exceptions::Dict])

Test univariate transformation `f` with `x`. Tests for:

1. type stability,

2. transformed value being in the image,

3. inverse,

4. log Jacobian determinant (using automatic differentiation).

`forwarddiff_exceptions` is a dictionary handling exceptions that ForwardDiff
cannot cope with at the moment. See [this
discussion](https://github.com/JuliaDiff/ForwardDiff.jl/issues/209).
"""
function test_univariate(t::UnivariateTransformation, x; AD_exceptions = Dict())
    @test length(t) == 1
    @test size(t) == ()
    y = @inferred t(x)
    @test y ∈ image(t)
    @test @inferred inverse(t, y) ≈ x
    lj = @inferred logjac(t, x)
    deriv = get(AD_exceptions, x, derivative(t, x))
    @test lj ≈ log(abs(deriv)) rtol=1e-5
end

function test_univariate_random(t::UnivariateTransformation; N=500,
                                AD_exceptions = Dict())
    for _ in 1:N
        test_univariate(t, rand_in(domain(t)), AD_exceptions = AD_exceptions)
    end
end

const logistic_AD_exceptions = Dict(-Inf => 0.0)

@testset "basic univariate transformations" begin
    test_univariate_random(Affine(1,2))
    test_univariate_random(NEGATION)
    test_univariate_random(LOGISTIC; AD_exceptions = logistic_AD_exceptions)
    test_univariate_random(LOGIT)
    test_univariate_random(REALCIRCLE)
    test_univariate_random(INVREALCIRCLE)
    test_univariate_random(EXP)
    test_univariate_random(LOG)
end

random_segment() = Segment(sort(randn(2))...)

function test_affine_bridge(x, y)
    x1, x2 = extrema(x)
    y1, y2 = extrema(y)
    a = affine_bridge(x, y)
    if !isincreasing(a)
        y2, y1 = y1, y2
    end
    @test a(x1) ≈ y1
    @test a(x2) ≈ y2
    @test a(x) ≈ y
end

@testset "affine composition" begin
    a = Affine(2, 3) ∘ Affine(1, 9)
    test_univariate_random(a)
    @test domain(a) == ℝ
    @test image(a) == ℝ
end

@testset "non-RR stable composition" begin
    c = EXP ∘ INVREALCIRCLE     # ℝ ↦ ℝ⁺ and (-1,1) ↦ ℝ, not RR stable
    @test ContinuousTransformations.RR_stability(c) == ContinuousTransformations.NotRRStable()
    @test_throws MethodError domain(c)
    @test_throws MethodError image(c)
end

@testset "affine bridge" begin
    test_affine_bridge(ℝ, ℝ)
    for _ in 1:100
        test_affine_bridge(random_segment(), random_segment())
        test_affine_bridge(PositiveRay(randn()), PositiveRay(randn()))
        test_affine_bridge(NegativeRay(randn()), NegativeRay(randn()))
        test_affine_bridge(PositiveRay(randn()), NegativeRay(randn()))
        test_affine_bridge(NegativeRay(randn()), PositiveRay(randn()))
    end
end

@testset "interval negation" begin
    @test NEGATION(Segment(1.0, 2.0)) == Segment(-2.0, -1.0)
    @test NEGATION(PositiveRay(9.0)) == NegativeRay(-9.0)
    @test NEGATION(NegativeRay(7.0)) == PositiveRay(-7.0)
    @test NEGATION(ℝ) == ℝ
end

function test_bridge(dom, img)
    t = @inferred bridge(dom, img)
    @test image(t) == img
    @test domain(t) == dom
    test_univariate_random(t)
end

@testset "bridge tests" begin
    test_bridge(ℝ, Segment(1, 2))
    test_bridge(ℝ, PositiveRay(9.0))
    test_bridge(ℝ, NegativeRay(-7.0))
    test_bridge(ℝ, ℝ)
    test_bridge(Segment(1,2), ℝ)
    test_bridge(PositiveRay(9.0), ℝ)
    test_bridge(NegativeRay(-7.0), ℝ)
end

@testset "show" begin
    @test repr(EXP) == "x ↦ exp(x)"
    @test repr(LOG) == "x ↦ log(x)"
    @test repr(REALCIRCLE) == "x ↦ realcircle(x)"
    @test repr(INVREALCIRCLE) == "x ↦ realcircle⁻¹(x)"
    @test repr(LOGISTIC) == "x ↦ logistic(x)"
    @test repr(LOGIT) == "x ↦ logit(x)"
    @test repr(NEGATION) == "x ↦ -x"
    @test repr(Affine(1,0)) == "x ↦ x"
    @test repr(Affine(2,0)) == "x ↦ 2.0⋅x"
    @test repr(Affine(2,3)) == "x ↦ 2.0⋅x + 3.0"
end

@testset "show MIME method" begin
    io = IOBuffer()
    t = EXP
    show(io, MIME"text/plain"(), t)
    s = String(take!(io))
    @test s == repr(t) * "\n"
end

# unit vector transformation

@testset "UnitVector calculations" begin
    for _ in 1:1000
        u = UnitVector(rand(2:10))
        x = randn(length(u))
        y, lj = transform_and_logjac(u, x)
        lj2 = logjac(u, x)
        @test lj == lj2
        x2 = inverse(u, y)
        @test norm(y, 2) ≈ 1.0
        @test logdet(jacobian(x -> transform(u, x)[1:(end-1)], x)) ≈ lj
        @test x ≈ x2
    end
end

@testset "UnitVector misc" begin
    @test repr(UnitVector(4)) == "x ↦ UnitVector(4)(x)"
    @test_throws ArgumentError UnitVector(-7)
    @test_throws ArgumentError UnitVector(0)
    @test transform(UnitVector(1), Float64[]) == [1.0]
    @test_throws ArgumentError transform(UnitVector(4), ones(4))
end

@testset "CorrelationCholeskyFactor calculations, full logjac" begin
    for _ in 1:1000
        n = rand(2:10)
        c = CorrelationCholeskyFactor(n)
        z = randn(length(c))
        L, lj = transform_and_logjac(c, z)
        lj2 =logjac(c, z)
        @test lj == lj2
        z2 = inverse(c, L)
        @test L isa LowerTriangular
        @test size(L) == (n, n)
        Σ = L*L'
        @test all(diag(Σ) .≈ 1)         # unit diagonal
        @test all(eigvals(Σ) .> 0)      # PD
        @test z ≈ z2
        # test Jacobian for this transformation
        J = jacobian(z) do z
            L = transform(c, z)
            lower_to_vec(L)
        end
        @test logdet(J) ≈ lj
        # test Jacobian for reconstruction of a full matrix
        # NOTE: test is here because it is related and we already generated L
        z_full = lowerdiag_to_vec(L)
        J_full = jacobian(z_full) do z
            L = vec_to_lowerdiag(z, n)
            Ω = L*L'
            lowerdiag_to_vec(Ω)
        end
        @test logdet(J_full) ≈ lkj_correlation_cholesky_logpdf(L, 1.0)
    end
end

@testset "CorrelationCholeskyFactor misc" begin
    @test repr(CorrelationCholeskyFactor(4)) ==
        "x ↦ CorrelationCholeskyFactor(4)(x)"
    @test_throws ArgumentError CorrelationCholeskyFactor(-7)
    @test_throws ArgumentError CorrelationCholeskyFactor(0)
    @test transform(CorrelationCholeskyFactor(1), Float64[]) ==
        LowerTriangular(fill(1.0, 1, 1))
    @test_throws ArgumentError transform(CorrelationCholeskyFactor(4), ones(4))
end

# array transformations

"""
    rand_Inf!(x, [p])

Replace each element of `x` with Inf or -Inf (equal probability), total with IID
probability `p`.
"""
function rand_Inf!(x, p = 0.02)
    for i in 1:length(x)
        a = rand()
        a < p && (x[i] = Inf)
        rand() < 0.5 && (x[i] *= -1)
    end
end

"""
Test array transformations (method consistency).
"""
function test_array_transformation(t, dims; N = 500)
    at = ArrayTransformation(t, dims)
    @test image(at) == image(t)
    @test domain(at) == domain(t)
    @test isincreasing(at) == isincreasing(t)
    @test length(at) == prod(dims)
    @test size(at) == dims
    @test_throws DimensionMismatch at(ones(dims .+ 1))
    @test_throws DimensionMismatch logjac(at, ones(dims .+ 1))
    for _ in 1:N
        x = randn(dims)
        @test @inferred logjac(at, x) == sum(logjac.(t, x))
        ## log jacobian may be meaningless at Inf, introduce Inf's after
        rand_Inf!(x)
        y = t.(x)
        @test @inferred at(x) == y
        @test @inferred inverse(at, y) == inverse.(t, y)
    end
end

@testset "array transformations" begin
    test_array_transformation(bridge(ℝ, Segment(1,2)), (2,3))
    test_array_transformation(bridge(ℝ, ℝ), (4,5))
    test_array_transformation(EXP, (3,7,2))
    test_array_transformation(REALCIRCLE, (3,2))
    @test_throws ArgumentError ArrayTransformation(EXP, -1, 2)
    @test_throws MethodError ArrayTransformation(EXP, "a fish")
    @test repr(ArrayTransformation(EXP, 2, 3)) ==
        repr(EXP) * " for (2, 3) elements"
    @test repr(ArrayTransformation(EXP, 2)) == repr(EXP) * " for 2 elements"
end

# transformation tuple

@testset "transformation tuple univariate" begin
    ts = bridge.(Ref(ℝ),
                 (PositiveRay(1.0), NegativeRay(1.0), NegativeRay(1.0),
                  ℝ, Segment(0.0,1.0)))
    tt = TransformationTuple(ts)
    @test length(tt) == sum(length, ts)
    @test image(tt) == image.(ts)
    @test domain(tt) == domain.(ts)
    @test repr(tt) == """
TransformationTuple
    x[1] ↦ exp(x[1]) + 1.0
    x[2] ↦ -exp(x[2]) + 1.0
    x[3] ↦ -exp(x[3]) + 1.0
    x[4] ↦ x[4]
    x[5] ↦ logistic(x[5])"""
    for i in 1:length(ts)
        @test tt[i] == ts[i]
    end
    x = randn(length(tt))
    y = @inferred tt(x)
    @test y == map((t,x) -> t(x), ts, tuple(x...))
    @test @inferred logjac(tt, x) == sum(map(logjac, ts, tuple(x...)))
    @test @inferred inverse(tt, y) == [map(inverse, ts, y)...]
end

@testset "transformation tuple mixed" begin
    ts = (EXP, ArrayTransformation(REALCIRCLE, 2))
    tt = TransformationTuple(ts)
    @test length(tt) == sum(length, ts) == 3
    @test image(tt) == image.(ts)
    @test repr(tt) == """
TransformationTuple
    x[1] ↦ exp(x[1])
    x[2:3] ↦ realcircle(x[2:3]) for 2 elements"""
    for i in 1:length(ts)
        @test tt[i] == ts[i]
    end
    x = randn(length(tt))
    y = @inferred tt(x)
    @test y == (ts[1](x[1]), ts[2](x[2:3]))
    @test @inferred logjac(tt, x) == logjac(ts[1], x[1]) + logjac(ts[2], x[2:3])
    @test @inferred inverse(tt, y) == vcat(inverse.(ts, y)...)
end

@testset "transformation tuple inference" begin
    t = TransformationTuple(bridge(ℝ, Segment(0.0,10.0)),
                            ArrayTransformation(Affine(1,0), 2))
    @inferred t(ones(3))
    @inferred logjac(t, ones(3))
    @inferred inverse(t, (1.0, ones(2)))
end

# log likelihood transform

@testset "log likelihood transformation" begin
    ℓ1(x) = 0.3*log(x) + 0.6*log(1-x) # unnormalized Beta, on (0, 1)
    ℓ2(x) = 2*log(x) - 0.3*x          # unnormalized Γ, on (0, ∞)
    ℓ(x) = ℓ1(x[1]) + ℓ2(x[2])

    t1 = bridge(ℝ, Segment(0,1))
    t2 = bridge(ℝ, PositiveRay(0))
    tℓ = TransformLogLikelihood(ℓ, t1, t2)

    @test get_transformation(tℓ) == TransformationTuple(t1, t2)
    @test get_loglikelihood(tℓ) ≡ ℓ

    # other constructor

    tℓ2 = TransformLogLikelihood(ℓ, (t1, t2))
    @test get_transformation(tℓ) == get_transformation(tℓ2)
    @test get_loglikelihood(tℓ) ≡ get_loglikelihood(tℓ2)

    for _ in 1:100
        x = randn(length(tℓ))
        @test tℓ(x) ≈ ℓ((t1(x[1]), t2(x[2]))) + logjac(t1, x[1]) + logjac(t2, x[2])
    end

    @test repr(tℓ) == """
TransformLogLikelihood of length 2, with TransformationTuple
    x[1] ↦ logistic(x[1])
    x[2] ↦ exp(x[2])"""
end

# transforming distributions

@testset "transform distribution with array transformation" begin
    μ = [-0.117965, -0.263465, -0.932187]
    A = [-0.15368 1.12831 0.364249;
         1.63777 1.5392 0.101908;
         -1.22376 -1.11266 0.246365]
    Σ = A'*A                   # positive semidefinite, positive definite w.p. 1
    Dx = MvNormal(μ, Σ)
    t = ArrayTransformation(EXP, 3)
    Dy = TransformDistribution(Dx, t)
    Dz = MvLogNormal(μ, Σ)
    @test length(Dy) == length(Dx)
    @test get_transformation(Dy) ≡ t
    @test get_distribution(Dy) ≡ Dx
    for _ in 1:1000
        x = rand(Dx)
        y = Dy(x)
        l = logpdf(Dz, y)       # true logpdf from distribution
        @test logpdf_in_domain(Dy, x) ≈ l
        @test logpdf_in_image(Dy, y) ≈ l
    end
    # test below is somewhat weak, but acceptable; large variance
    mean_sim = mean(rand(Dy) for _ in 1:100000)
    @test maximum(abs.(mean_sim .- mean(Dz))) ≤ 1
end

@testset "transform univariate distribution" begin
    μ = 2.7
    σ = 1.3
    Dx = Normal(μ, σ)
    t = EXP
    Dy = TransformDistribution(Dx, t)
    Dz = LogNormal(μ, σ)
    @test length(Dy) == length(Dx)
    @test get_transformation(Dy) ≡ t
    @test get_distribution(Dy) ≡ Dx
    for _ in 1:1000
        x = rand(Dx)
        y = Dy(x)
        l = logpdf(Dz, y)       # true logpdf from distribution
        @test logpdf_in_domain(Dy, x) ≈ l
        @test logpdf_in_image(Dy, y) ≈ l
    end
    # test below is somewhat weak, but acceptable; large variance
    mean_sim = mean(rand(Dy) for _ in 1:100000)
    @test maximum(abs.(mean_sim .- mean(Dz))) ≤ 1
end

# misc utilities

@testset "ungrouping" begin
    a = bridge(ℝ, Segment(0, 5))
    b = ArrayTransformation(bridge(ℝ, PositiveRay(7.0)), 2)
    t = TransformationTuple((a, b))
    N = 100
    A = [randn(3) for _ in 1:N]
    B = ungrouping_map(Vector, t, A)
    C = ungrouping_map(Array, t, A)
    for i in 1:N
        z = t(A[i])
        @test B[1][i] == z[1]
        @test B[2][1][i] == z[2][1]
        @test B[2][2][i] == z[2][2]
        @test C[1][i] == z[1]
        @test C[2][i, :] == z[2]
    end
end

@testset "all_finite" begin
    # scalars
    @test all_finite(9)
    @test all_finite(1.0)
    @test !all_finite(Inf)
    @test !all_finite(NaN)
    @test !all_finite(-Inf)
    # arrays
    @test all_finite(fill(9.0, 2, 2))
    @test !all_finite([1.0 Inf;
                       0.0 2.0])
    @test !all_finite([NaN])
    # tuples
    @test all_finite((1.0, [2.0, 3.0]))
    @test !all_finite((1.0, [2.0, NaN]))
    # iterables
    @test all_finite((Float64(i)^2 for i in 1:10))
    let k = 0
        @test !all_finite(((k += 1; i ≤ 5 ? i : Inf) for i in 1:10))
        @test k == 6            # test early termination
    end
end