text:
"""
```julia
malloc(size::Integer)
```
Libc `malloc` function, accessed by direct StaticCompiler-safe `llvmcall`.

Allocate `size` bytes of memory and return a pointer to that memory.

See also: `calloc`, `free`.

## Examples
```julia
julia> p = malloc(500)
Ptr{UInt8} @0x00007ff0e9e74290

julia> free(p)
0
```
"""
@inline malloc(size::Integer) = malloc(Int64(size))
@inline function malloc(size::Int64)
    Base.llvmcall(("""
    ; External declaration of the `malloc` function
    declare i8* @malloc(i64)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %size) #0 {
      %ptr = call i8* (i64) @malloc(i64 %size)
      %jlp = ptrtoint i8* %ptr to i64
      ret i64 %jlp
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Ptr{UInt8}, Tuple{Int64}, size)
end
@inline malloc(size::Unsigned) = malloc(UInt64(size))
@inline function malloc(size::UInt64)
    Base.llvmcall(("""
    ; External declaration of the `malloc` function
    declare i8* @malloc(i64)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %size) #0 {
      %ptr = call i8* (i64) @malloc(i64 %size)
      %jlp = ptrtoint i8* %ptr to i64
      ret i64 %jlp
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Ptr{UInt8}, Tuple{UInt64}, size)
end

"""
```julia
calloc([n], size::Integer)
```
Libc `calloc` function, accessed by direct StaticCompiler-safe `llvmcall`.

Allocate `size` bytes of zero-initialized memory and return a pointer to that memory.
As `malloc`, but initializes the memory to all zero.

See also: `malloc`, `free`.

## Examples
```julia
julia> p = calloc(100*sizeof(Int64))
Ptr{UInt8} @0x00007fb74ff04360

julia> MallocArray{Int64}(p, 10, 10)
10×10 MallocMatrix{Int64}:
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0

julia> free(p)
0
```
"""
@inline calloc(size::Integer) = calloc(Int64(size))
@inline calloc(size::Int64) = calloc(1, size)
@inline function calloc(n::Int64, size::Int64)
    Base.llvmcall(("""
    ; External declaration of the `calloc` function
    declare i8* @calloc(i64, i64)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %n, i64 %size) #0 {
      %ptr = call i8* (i64, i64) @calloc(i64 %n, i64 %size)
      %jlp = ptrtoint i8* %ptr to i64
      ret i64 %jlp
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Ptr{UInt8}, Tuple{Int64, Int64}, n, size)
end

"""
```julia
free(ptr::Ptr)
```
Libc `free` function, accessed by direct StaticCompiler-safe `llvmcall`.

Free memory that has been previously allocated with `malloc`.

See also: `calloc`, `malloc`.

## Examples
```julia
julia> p = malloc(500)
Ptr{UInt8} @0x00007ff0e9e74290

julia> free(p)
0
```
"""
@inline free(ptr::Ptr) = free(Ptr{UInt8}(ptr))
@inline function free(ptr::Ptr{UInt8})
    Base.llvmcall(("""
    ; External declaration of the `free` function
    declare void @free(i8*)

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %jlp) #0 {
      %ptr = inttoptr i64 %jlp to i8*
      call void (i8*) @free(i8* %ptr)
      ret i32 0
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{UInt8}}, ptr)
end

"""
```julia
memset!(a, char::Integer, nbytes::Integer=sizeof(a))
```
Libc `memset` function, accessed by direct StaticCompiler-safe `llvmcall`.

Set `nbytes` bytes of the array or memory region `a` to the `Char`/`UInt8`
conversion of the integer `char`.

## Examples
```julia
julia> a = rand(5,5)
5×5 Matrix{Float64}:
 0.252808   0.6125    0.947215   0.0966341  0.637651
 0.736149   0.527729  0.928291   0.725644   0.832734
 0.704827   0.990302  0.0380948  0.768337   0.891583
 0.0826808  0.833624  0.364925   0.230345   0.366826
 0.301975   0.113886  0.329196   0.772636   0.0156762

julia> memset!(a, 0)
0

julia> a
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
@inline memset!(a, char, nbytes=sizeof(a)) = GC.@preserve a memset!(pointer(a), char, nbytes)
@inline memset!(a::AbstractMallocdMemory, char, nbytes=sizeof(a)) = memset!(pointer(a), char, nbytes)
@inline memset!(ptr::Ptr, char::Integer, nbytes::Integer) = memset!(Ptr{UInt8}(ptr), char, nbytes)
@inline memset!(ptr::Ptr{UInt8}, char::Integer, nbytes::Integer) = memset!(ptr, Int64(char), Int64(nbytes))
@inline function memset!(ptr::Ptr{UInt8}, char::Int64, nbytes::Int64)
    Base.llvmcall(("""
    ; External declaration of the `memset` function
    ; Function Attrs: argmemonly nounwind
    declare void @memset(i8* nocapture writeonly, i64, i64) #0

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %jlp, i64 %value, i64 %n) #1 {
      %ptr = inttoptr i64 %jlp to i8*
      call void @memset(i8* %ptr, i64 %value, i64 %n)
      ret i32 0
    }

    attributes #0 = { argmemonly nounwind }
    attributes #1 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{UInt8}, Int64, Int64}, ptr, char, nbytes)
end

"""
```julia
memcpy!(a, b, n=length(b))
```
Libc `memcpy` function, accessed by direct StaticCompiler-safe `llvmcall`.

Copy `n` elements from array `b` to array `a`.

## Examples
```julia
julia> a = rand(3)
3-element Vector{Float64}:
 0.8559883493421137
 0.4203692766310769
 0.5728354965961716

julia> memcpy!(a, ones(3))
0

julia> a
3-element Vector{Float64}:
 1.0
 1.0
 1.0
```
"""
@inline memcpy!(a, b) = memcpy!(a, b, length(b))
@inline memcpy!(a, b, n::Int64) = GC.@preserve a b memcpy!(pointer(a), pointer(b), n)
@inline memcpy!(dst::Ptr, src::Ptr{T}, n::Int64) where {T} = memcpy!(Ptr{UInt8}(dst), Ptr{UInt8}(src), n*sizeof(T))
@inline function memcpy!(dst::Ptr{UInt8}, src::Ptr{UInt8}, nbytes::Int64)
    Base.llvmcall(("""
    ; External declaration of the `memcpy` function
    ; Function Attrs: argmemonly nounwind
    declare void @llvm.memcpy.p0i8.p0i8.i64(i8* nocapture writeonly, i8* nocapture readonly, i64, i1) #0

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %jldest, i64 %jlsrc, i64 %nbytes) #1 {
      %dest = inttoptr i64 %jldest to i8*
      %src = inttoptr i64 %jlsrc to i8*
      call void @llvm.memcpy.p0i8.p0i8.i64(i8* %dest, i8* %src, i64 %nbytes, i1 false)
      ret i32 0
    }

    attributes #0 = { argmemonly nounwind }
    attributes #1 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{UInt8}, Ptr{UInt8}, Int64}, dst, src, nbytes)
end

"""
```julia
memcmp(a, b, n::Int64)
```
Libc `memcmp` function, accessed by direct StaticCompiler-safe `llvmcall`.

Compare the first `n` bytes of `a` and `b`, returning
* a positive value if the first `n` bytes of `a` are greater than the first `n` bytes of `b`
* a negative value if the first `n` bytes of `a` are less than the first `n` bytes of `b`
* `0` the first `n` bytes of `a` are equal to the first `n` bytes of `b`

## Examples
```julia
julia> memcmp(c"foo", c"foo", 3)
0

julia> memcmp(c"foo", c"bar", 3)
4
```
"""
@inline memcmp(a, b, n::Int64) = GC.@preserve a b memcmp(pointer(a), pointer(b), n)
@inline memcmp(a::Ptr, b::Ptr, n::Int64) = memcpy!(Ptr{UInt8}(a), Ptr{UInt8}(b), n)
@inline function memcmp(a::Ptr{UInt8}, b::Ptr{UInt8}, nbytes::Int64)
    Base.llvmcall(("""
    ; External declaration of the `memcmp` function
    declare i32 @memcmp(i8*, i8*, i64)

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %jla, i64 %jlb, i64 %nbytes) #0 {
      %a = inttoptr i64 %jla to i8*
      %b = inttoptr i64 %jlb to i8*
      %cmp = call i32 @memcmp(i8* %a, i8* %b, i64 %nbytes)
      ret i32 %cmp
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{UInt8}, Ptr{UInt8}, Int64}, a, b, nbytes)
end

"""
```julia
time()
```
Libc `time` function, accessed by direct StaticCompiler-safe `llvmcall`.

Return, as an `Int64`, the current time in seconds since the beginning of the
current Unix epoch on 00:00:00 UTC, January 1, 1970.

## Examples
```julia
julia> StaticTools.time()
1651105298
```
"""
@inline function time()
    Base.llvmcall(("""
    ; External declaration of the `time` function
    declare i64 @time(i64*)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main() {
      %time = call i64 @time(i64* null)
      ret i64 %time
    }
    """, "main"), Int64, Tuple{}, size)
end

"""
```julia
usleep(μsec::Integer)
```
Libc `usleep` function, accessed by direct StaticCompiler-safe `llvmcall`.

Suspend execution of the calling thread for (at least) μsec microseconds.

## Examples
```julia
julia> usleep(1000000)
0
```
"""
@inline usleep(μsec::Integer) = malloc(Int64(μsec))
@inline function usleep(μsec::Int64)
    Base.llvmcall(("""
    ; External declaration of the `usleep` function
    declare i32 @usleep(i64)

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %usec) #0 {
      %status = call i32 (i64) @usleep(i64 %usec)
      ret i32 %status
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Int64}, μsec)
end

"""
```julia
system(s)
```
Libc `system` function, accessed by direct StaticCompiler-safe `llvmcall`.

Pass the null-terminated string (or pointer thereto) `s` to the libc `system`
function for evaluation.

Returns `0` on success.

## Examples
```julia
julia> StaticTools.system(c"time echo hello")
hello

real    0m0.001s
user    0m0.000s
sys 0m0.000s
0
```
"""
@inline system(s::AbstractMallocdMemory) = system(pointer(s))
@inline system(s) = GC.@preserve s system(pointer(s))
@inline function system(s::Ptr{UInt8})
    Base.llvmcall(("""
    ; External declaration of the `system` function
    declare i32 @system(...)

    ; Function Attrs: nounwind ssp uwtable
    define i32 @main(i64 %jlstr) #0 {
      %str = inttoptr i64 %jlstr to i8*
      %status = call i32 (i8*, ...) bitcast (i32 (...)* @system to i32 (i8*, ...)*)(i8* %str)
      ret i32 %status
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int32, Tuple{Ptr{UInt8}}, s)
end

"""
```julia
strlen(s)
```
Libc `strlen` function, accessed by direct StaticCompiler-safe `llvmcall`.

Returns the length in bytes of the null-terminated string `s`, not counting the
terminating null character.

## Examples
```julia
julia> strlen("foo") # Not documented, but Julia strings are null-terminated in practice every time I've checked
3

julia> strlen(c"foo")
3
```
"""
@inline strlen(s::AbstractMallocdMemory) = strlen(pointer(s))
@inline strlen(s) = GC.@preserve s strlen(pointer(s))
@inline function strlen(s::Ptr{UInt8})
    Base.llvmcall(("""
    ; External declaration of the `strlen` function
    declare i64 @strlen(i8*)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %jlstr) #0 {
      %str = inttoptr i64 %jlstr to i8*
      %li = call i64 (i8*) @strlen (i8* %str)
      ret i64 %li
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int64, Tuple{Ptr{UInt8}}, s)
end

"""
```julia
strtod(s)
```
Libc `strtod` function, accessed by direct StaticCompiler-safe `llvmcall`.

Returns a `Float64` ("double") containing the number written out in decimal form
in null-terminated string `s`.

## Examples
```julia
julia> num, pbuf = StaticTools.strtod(c"3.1415")
(3.1415, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000010aeee946,)))

julia> num, pbuf = StaticTools.strtod(c"5")
(5.0, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000010d8f2bb1,)))
```
"""
@inline strtod(s::AbstractMallocdMemory) = strtod(pointer(s))
@inline strtod(s) = GC.@preserve s strtod(pointer(s))
@inline function strtod(p::Ptr{UInt8})
    pbuf = MemoryBuffer{1,Ptr{UInt8}}(undef)
    num = GC.@preserve pbuf strtod(p, pointer(pbuf))
    return num, pbuf
end
@inline function strtod(s::Ptr{UInt8}, p::Ptr{Ptr{UInt8}})
    Base.llvmcall(("""
    ; External declaration of the `strtod` function
    declare double @strtod(i8*, i8**)

    ; Function Attrs: nounwind ssp uwtable
    define double @main(i64 %jlstr, i64 %jlp) #0 {
      %str = inttoptr i64 %jlstr to i8*
      %ptr = inttoptr i64 %jlp to i8**
      %d = call double (i8*, i8**) @strtod (i8* %str, i8** %ptr)
      ret double %d
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Float64, Tuple{Ptr{UInt8}, Ptr{Ptr{UInt8}}}, s, p)
end

"""
```julia
strtol(s)
```
Libc `strtol` function, accessed by direct StaticCompiler-safe `llvmcall`.

Returns an `Int64` ("long") containing the number written out in decimal form in
null-terminated string `s`.

## Examples
```julia
julia> num, pbuf = StaticTools.strtol(c"3.1415")
(3, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000010dd827f1,)))

julia> num, pbuf = StaticTools.strtol(c"5")
(5, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000015dbdda41,)))
```
"""
@inline strtol(s::AbstractMallocdMemory) = strtol(pointer(s))
@inline strtol(s) = GC.@preserve s strtol(pointer(s))
@inline function strtol(p::Ptr{UInt8}, base::Int32=Int32(10))
    pbuf = MemoryBuffer{1,Ptr{UInt8}}(undef)
    num = GC.@preserve pbuf strtol(p, pointer(pbuf), base)
    return num, pbuf
end
@inline function strtol(s::Ptr{UInt8}, p::Ptr{Ptr{UInt8}}, base::Int32)
    Base.llvmcall(("""
    ; External declaration of the `strtol` function
    declare i64 @strtol(i8*, i8**, i32)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %jlstr, i64 %jlp, i32 %base) #0 {
      %str = inttoptr i64 %jlstr to i8*
      %ptr = inttoptr i64 %jlp to i8**
      %li = call i64 (i8*, i8**, i32) @strtol (i8* %str, i8** %ptr, i32 %base)
      ret i64 %li
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), Int64, Tuple{Ptr{UInt8}, Ptr{Ptr{UInt8}}, Int32}, s, p, base)
end

"""
```julia
strtoul(s)
```
Libc `strtol` function, accessed by direct StaticCompiler-safe `llvmcall`.

Returns an `UInt64` ("unsigned long") containing the number written out in decimal
form in null-terminated string `s`.

## Examples
```julia
julia> num, pbuf = StaticTools.strtoul(c"3.1415")
(0x0000000000000003, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000010d6976a1,)))

julia> num, pbuf = StaticTools.strtoul(c"5")
(0x0000000000000005, ManualMemory.MemoryBuffer{1, Ptr{UInt8}}((Ptr{UInt8} @0x000000015ed45d11,)))
```
"""
@inline strtoul(s::AbstractMallocdMemory) = strtoul(pointer(s))
@inline strtoul(s) = GC.@preserve s strtoul(pointer(s))
@inline function strtoul(p::Ptr{UInt8}, base::Int32=Int32(10))
    pbuf = MemoryBuffer{1,Ptr{UInt8}}(undef)
    num = GC.@preserve pbuf strtoul(p, pointer(pbuf), base)
    return num, pbuf
end
@inline function strtoul(s::Ptr{UInt8}, p::Ptr{Ptr{UInt8}}, base::Int32)
    Base.llvmcall(("""
    ; External declaration of the `strtoul` function
    declare i64 @strtoul(i8*, i8**, i32)

    ; Function Attrs: nounwind ssp uwtable
    define i64 @main(i64 %jlstr, i64 %jlp, i32 %base) #0 {
      %str = inttoptr i64 %jlstr to i8*
      %ptr = inttoptr i64 %jlp to i8**
      %li = call i64 (i8*, i8**, i32) @strtoul (i8* %str, i8** %ptr, i32 %base)
      ret i64 %li
    }

    attributes #0 = { alwaysinline nounwind ssp uwtable }
    """, "main"), UInt64, Tuple{Ptr{UInt8}, Ptr{Ptr{UInt8}}, Int32}, s, p, base)
end

"""
```julia
parse(::Type{T}, s::Union{StaticString, MallocString})
```
Parse a number from a `StaticString` or `MallocString` `s`.

## Examples
```julia
julia> parse(Float64, c"3.141592")
3.141592

julia> parse(Int64, c"3.141592")
3
```
"""
@inline function Base.parse(::Type{Float64}, s::Union{StaticString, MallocString})
    num, pbuf = strtod(s)
    load(pointer(pbuf)) == pointer(s) && return NaN
    return num
end
@inline function Base.parse(::Type{Float64}, s::Ptr{UInt8})
    num, pbuf = strtod(s)
    load(pointer(pbuf)) == s && return NaN
    return num
end
@inline Base.parse(::Type{T}, s::Union{StaticString, MallocString, Ptr{UInt8}}) where {T <: AbstractFloat} = T(parse(Float64, s))

@inline function Base.parse(::Type{Int64}, s::Union{StaticString, MallocString, Ptr{UInt8}})
    num, pbuf = strtol(s)
    return num
end
@inline Base.parse(::Type{T}, s::Union{StaticString, MallocString, Ptr{UInt8}}) where {T <: Integer} = T(parse(Int64, s))

@inline function Base.parse(::Type{UInt64}, s::Union{StaticString, MallocString, Ptr{UInt8}})
    num, pbuf = strtoul(s)
    return num
end
@inline Base.parse(::Type{T}, s::Union{StaticString, MallocString, Ptr{UInt8}}) where {T <: Unsigned} = T(parse(UInt64, s))

# Convenient parsing for argv (slight type piracy)
@inline Base.parse(::Type{T}, argv::Ptr{Ptr{UInt8}}, n::Integer) where {T} = parse(T, MallocString(argv, n))

## ---


text:
################################################################################
# Impulse Transform
################################################################################
function impulse_transform(relative::Symbol, joint::Joint, xa::AbstractVector,
        qa::Quaternion, xb::AbstractVector, qb::Quaternion) where {T}
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=true)
    Diagonal([sones(3); 0.5 * sones(3)]) * transpose([X Q]) #TODO: 0.5 Q
end

################################################################################
# Derivatives
################################################################################
function impulse_map_jacobian(relative::Symbol, jacobian::Symbol, joint::Joint, pbody::Node{T}, cbody::Node{T}, λ) where T
    # ∂(G*λ)/∂(x,q)
    p = impulse_projector(joint) * λ
    impulse_transform_jacobian(relative, jacobian,
        joint,
        current_configuration(pbody.state)...,
        current_configuration(cbody.state)...,
        p)
end

ext: jl


text:
using Genie.Router
using CalculatorController
import Genie.Renderer.Json: json
using Genie.Requests

route("/") do
  serve_static_file("welcome.html")
end

route("/beep") do
  print(println("\007"))
end

route("/params", method = POST) do
  messages = jsonpayload()
  out = Dict()
  for mk in keys(messages)
    message = messages[mk]
    out[mk] = [CalculatorController.params_json([point[k] for k in CalculatorController.KEY_ORDER]) for point in message]
  end
  return out |> json
end

route("/opt", method = POST) do
  messages = jsonpayload()
  out = Dict()
  for mk in keys(messages)
    message = messages[mk]
    last_result = nothing
    last_weight = nothing
    out[mk] = []
    for point in message
      p = [point[k] for k in CalculatorController.PARAM_ORDER]
      if CalculatorController.KEY_ORDER[1] in keys(point)
        V0 = [convert(Float64, point[v]) for v in CalculatorController.KEY_ORDER]
      elseif last_result !== nothing
        V0 = last_result
      else
        V0 = nothing
      end
      if "weights" in keys(point)
        w = [convert(Float64, point["weights"][k]) for k in CalculatorController.PARAM_ORDER]
      elseif last_weight !== nothing
        w = last_weight
      else
        w = CalculatorController.weights
      end
      result = CalculatorController.opt_json(p; V0=V0, w=w)
      last_result = [convert(Float64, result[v]) for v in CalculatorController.KEY_ORDER]
      last_weight = w
      append!(out[mk], result)
    end
  end
  return out |> json
end


text:
abstract Solver

# A Newton Raphsoin solver
immutable NRSolver <: Solver
    rel_tol::Float64
    abs_tol::Float64
    alt_tol::Float64
    max_iters::Int
    err_on_nonconv::Bool
end

NRSolver(;rel_tol=0.0::Float64, abs_tol=0.0::Float64,
          max_iters=10::Int, err_on_nonconv=true::Bool, alt_tol = 0.0) =
    NRSolver(rel_tol, abs_tol, alt_tol, max_iters, err_on_nonconv)

# Solves a nonlinear problem by iterating until convergence
function solve(solver::NRSolver, fp::FEProblem, exporter::AbstractDataExporter)
    println("Starting Newton-Raphson solver..")
    load = extload(fp)
    force_imbalance = similar(load)
    iteration = 0
    tstep = 0
    n_print = 0

    K = create_sparse_structure(fp::FEProblem)
    colptrs = get_colptrs(K, fp::FEProblem)

    for t in 0.0 #TODO Implement proper time stepping
        println("Current time $t")
        iteration = 0
        tstep += 1
        updatebcs!(fp, t)
        while true
            iteration += 1
            if iteration > solver.max_iters
                warn("Max iterations hit.")
                if solver.err_on_nonconv
                    error("Did not converge")
                end
                break
            end

            int_f = assemble_intf(fp)
            @devec force_imbalance[:] = load .- int_f

            abs_res = norm(force_imbalance)
            rel_res = norm(force_imbalance) / norm(load)

            @printf("\tIteration %d, relative residual %1.4e, absolute residual %1.4e", iteration, rel_res, abs_res)

            print("\n\tDofType residuals:\n \t\t")
            for (dof_type, eqs) in fp.doftype_eqs
                @devec r = sum(sqr(int_f[eqs]))
                @printf("%s: %1.2e   \t", dof_type, sqrt(r))
            end

            if  abs_res < solver.abs_tol
                println("\nConverged with absolute tolerance!")
                break
            end

            if rel_res < solver.rel_tol
                println("\nConverged with relative tolerance!")
                break
            end

            assembleK!(K, fp, colptrs, fp.dof_vals)

            #du = cholfact(Symmetric(K, :L)) \ force_imbalance
            du = K \ force_imbalance

            print("\n\tDofType updates:\n \t\t")
            for (dof_type, eqs) in fp.doftype_eqs
                @devec r = sum(sqr(du[eqs]))
                @printf("%s: %1.2e   \t", dof_type, sqrt(r))
            end
            updatedofs!(fp, du)

            if norm(du) < solver.alt_tol
                     println("\n Converged with alternative tolerance!")
                break
            end

            print("\n\n")
        end
    update_feproblem(fp)

    n_print += 1
    write_data(fp, exporter, n_print)
    end

end


text:
# sanity checks

using NeuralVerification
using Test

function printtest(solver, problem_sat, problem_unsat)
    println(typeof(solver))
    res_sat = solve(solver, problem_sat)
    res_unsat = solve(solver, problem_unsat)

    function col(s, rev = false)
        cols = [:green, :red]
        rev && reverse!(cols)
        s == :SAT   && return cols[1]
        s == :UNSAT && return cols[2]
        return (:yellow) #else
    end

    print("\tSAT test.   Result: "); printstyled("$(res_sat.status)\n", color = col(res_sat.status))
    print("\tUNSAT test. Result: "); printstyled("$(res_unsat.status)\n", color = col(res_unsat.status, true))
    println("_"^70, "\n")
end

printtest(solvers::Vector, p1, p2) = ([printtest(s, p1, p2) for s in solvers]; nothing)

at = @__DIR__
small_nnet = read_nnet("$at/../examples/networks/small_nnet.nnet")

# The input set is always [-1:1]
input_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])
input_hpoly  = HPolytope(input_hyper)

out_hyper_30_80 = Hyperrectangle(low = [20.0], high = [90.0])
out_hyper_50    = Hyperrectangle(low = [-1.0], high = [50.0]) # includes points in the output region ie y > 30.5

problem_sat_hyper_hyper           = Problem(small_nnet, input_hyper, out_hyper_30_80)                      # 40.0 < y < 60.0
problem_unsat_hyper_hyper         = Problem(small_nnet, input_hyper, out_hyper_50)                         # -1.0 < y < 50.0
problem_sat_hpoly_hpoly_bounded   = Problem(small_nnet, input_hpoly, HPolytope(out_hyper_30_80))
problem_unsat_hpoly_hpoly_bounded = Problem(small_nnet, input_hpoly, HPolytope(out_hyper_50))
# halfspace constraints:
problem_sat_hyper_hs              = Problem(small_nnet, input_hyper, HPolytope([HalfSpace([1.], 100.)]))     # y < 100.0
problem_unsat_hyper_hs            = Problem(small_nnet, input_hyper, HPolytope([HalfSpace([1.], 10.)]))      # y < 10.0

# GROUP 1           # Input: HPolytope, Output: HPolytope
# group1 = [MaxSens(), ExactReach(), Ai2()]
group1 = [MaxSens(resolution = 0.6), ExactReach()] # Ai2 is 100% broken right now so dropping it
for solver in group1
    printtest(solver, problem_sat_hpoly_hpoly_bounded, problem_unsat_hpoly_hpoly_bounded)
    sat   = solve(solver, problem_sat_hpoly_hpoly_bounded)
    unsat = solve(solver, problem_unsat_hpoly_hpoly_bounded)

    @test sat.status ∈ (:SAT, :Unknown)
    @test unsat.status ∈ (:UNSAT, :Unknown)
end

# GROUP 2, 3, 4     # Input: HPolytope, Output: HPolytope
glpk = GLPKSolverMIP()
group2 = [S(optimizer = glpk) for S in (NSVerify, MIPVerify, ILP)]
group3 = [ConvDual(), Duality(optimizer = glpk)]
group4 = [FastLin(), FastLip()]
group6 = [Reluplex(), Planet(GLPKSolverMIP())]

for solver in [group2; group3; group4; group6]
    printtest(solver, problem_sat_hyper_hs, problem_unsat_hyper_hs)
    sat   = solve(solver, problem_sat_hyper_hs)
    unsat = solve(solver, problem_unsat_hyper_hs)

    @test sat.status ∈ (:SAT, :Unknown)
    @test unsat.status ∈ (:UNSAT, :Unknown)
end

# GROUP 5, 6        # Input: Hyperrectangle, Output: Hyperrectangle
group5 = [ReluVal(max_iter = 10), DLV(), Sherlock(glpk, 1.0), BaB(optimizer = glpk)]

for solver in [group5;]
    printtest(solver, problem_sat_hyper_hyper, problem_unsat_hyper_hyper)
    sat   = solve(solver, problem_sat_hyper_hyper)
    unsat = solve(solver, problem_unsat_hyper_hyper)

    @test sat.status ∈ (:SAT, :Unknown)
    @test unsat.status ∈ (:UNSAT, :Unknown)
end

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

### Certify - only works for single hidden layer
tiny_nnet = read_nnet("$at/../examples/networks/tiny_nnet.nnet")
solver_certify = Certify()
inputSet  = Hyperrectangle([2.0], [.5])
outputSet = HPolytope(ones(1,1), [2.5])
problem_certify = Problem(tiny_nnet, inputSet, outputSet)
@test @no_error solve(solver_certify, problem_certify)


text:
npsList = DataFrame(NPS = [0.125,
	0.25,
	0.375,
	0.5,
	0.75,
	1,
	1.25,
	1.5,
	2,
	2.5,
	3,
	3.5,
	4,
	5,
	6,
	8,
	10,
	12,
	14,
	16,
	18,
	20,
	22,
	24,
	26,
	28,
	30,
	32,
	36,
	42,
	48] )


text:
module Koala

# new:
export @more, @dbg, @colon, keys_ordered_by_values, params, partition
export bootstrap_resample_of_mean, bootstrap_resample_of_median
export load_boston, load_ames, load_iris, datanow
export hasmissing, countmissing, ismissingtype, purify
export fit!, predict, rms, rmsl, rmslp1, rmsp, err, transform, inverse_transform
export ConstantRegressor
export IdentityTransformer, FeatureSelector
export default_transformer_X, default_transformer_y, clean!
export Machine
export learning_curve, cv, compete, @colon, @curve, @pcurve
export split_seen_unseen
export bootstrap_histogram_of_mean, bootstrap_histogram_of_mean!, PlotableDict
export bootstrap_histogram_of_median, bootstrap_histogram_of_median!

# for use in this module:
using Statistics
import DataFrames: DataFrame, AbstractDataFrame, names, eltypes, describe
import CSV
import StatsBase: sample, countmap
import HypothesisTests: OneSampleTTest, pvalue
# import ROC
using  RecipesBase # for plotting recipes
import Distributed: @distributed, nworkers

# reexport:
export describe, countmap, names, eltypes

# extended in this module:
import Base: show, isempty, split

# constants:
const COLUMN_WIDTH = 24           # for displaying dictionaries with `show`
const DEFAULT_SHOW_DEPTH = 1      # how deep to display fields of `BaseType` objects
const srcdir = dirname(@__FILE__) # the full path for this file

# functions to be extended (provided methods) in dependent packages:
function fit end
function predict end
function setup end
function transform end
function inverse_transform end
function default_transformer_X end
function default_transformer_y end

## ON TERMINOLOGY AND THE MAIN TYPES

# 1. *Transformers and schemes*. A subtype of `Transformer` is a
# `struct` identifying a type of data transformation whose objects
# (*transformers*) store associated hyperparemeters. For example, a
# `ToIntTransformer` object (defined in KoalaTransformer.jl) relabels
# non-numeric vectors with integer vectors, the first integer being
# given by the value of the hyperparameter `initial_label` (0 or 1,
# say). A transformer does *not* store the actual data encoding the
# transformation (the relabelling in our example). This information
# (called a *scheme*) is part of a transformer *machine* (see below).

# 2. *Models and predictors.* A subtype of `Model` is a `struct`
# associated with some learning algorithm whose objects (*models*)
# store associated hyperparameters. For example, a `TreeRegressor`
# object (defined in KoalaTrees.jl) has fields describing the growth
# of a decision tree, such as the minimim number of training patterns
# triggering a node split. Such parameters are distinct from variables
# that are *learned* by the model, such as a decision tree's actual
# structure, or the coefficients (weights) in a linear model. So a
# model, in our sense, is distinct from a *predictor*, which is part
# of a superivised learning *machine* (see below).

# 3. A *transformer machine* is an object used for actually
# transforming data. It is fit once (to some possibly different data,
# usually the training data) during instantiation. This fitting wraps
# the provided transformer with the scheme encoding the
# transformation.

# 4. A *supervised learning machine* is an object eventually used to
# make predictions but is not pre-trained. Rather, during
# instantiation, it pre-transforms (once) the data to be used for
# training *and* testing, into a form suitable for the algorithm
# associated with model supplied (eg, by standardizing or one-hot
# encoding inputs). The hyperparmaters stored in the model may be
# changed after instantiation and retraining the machine retrains
# using the new hyperparameters..

# 5. *Summary*. Tranformers and models just store
# hyperparemters. Machines do the actual transforming and
# learning/predicting. So, a supervised learning machine wraps a model
# in data and, after training, also a predictor; a transformer machine
# wraps a transformer with the scheme encoding the actual
# transformation.

# All Koala types and types created in its associated libraries
# subtype the following:

abstract type BaseType end

## HELPERS AND CONVENIENCE FUNCTIONS

# extension of `show` method for `BaseType` objects:
include("show.jl")

ismissingtype(T::Type) = isa(T, Union) && T.a == Missing

principaltype(T::Union) = T.b::Type # for extacting pure type from missing type

hasmissing(v::AbstractVector) = findfirst(ismissing, v) != nothing

countmissing(v::AbstractVector) = count(ismissing, v)

""" convert a vector of eltype `Union{T, Missing}` to one of eltype `T`"""
function purify(v::AbstractVector)
    T = eltype(v)
    if ismissingtype(T)
        !hasmissing(v) || @error "Can't purify a vector with missing values."
        return convert(Array{principaltype(T)}, v)
    else
        return v
    end
end

# """convenience macro for printing variable values (eg for degugging)"""
# macro dbg(v)
#     esc(quote
#         print((@colon $v), "=", $v, " ")
#     end)
# end

# macro dbg(v1, v2)
#     esc(quote
#         print((@colon $v1), "=", $v1, " ")
#         println((@colon $v2), "=", $v2, " ")
#     end)
# end

# macro dbg(v1, v2, v3)
#     esc(quote
#         print((@colon $v1), "=", $v1, " ")
#         print((@colon $v2), "=", $v2, " ")
#         println((@colon $v3), "=", $v3, " ")
#     end)
# end

""" a version of warn that only warns given a *non-empty* string."""
function softwarn(str)
    if !isempty(str)
        @warn str
    end
end

"""Load a well-known public regression dataset with nominal features."""
function load_boston()
    df = CSV.read(joinpath(srcdir, "data", "Boston.csv"),
                  categorical=false, allowmissing=:none)
    features = filter(names(df)) do f
        f != :MedV
    end
    X = df[features]
    for j in 1:size(X, 2)
        X[j] = purify(X[j])
    end
    y = purify(df[:MedV])
    return X, y
end

"""Load a reduced version of the well-known Ames Housing dataset,
having six numerical and six categorical features."""
function load_ames()
    df = CSV.read(joinpath(srcdir, "data", "reduced_ames.csv"),
                  categorical=false, allowmissing=:none)
    features = filter(names(df)) do f
        f != :target
    end
    X = df[features]
    for j in 1:size(X,2)
        X[j] = purify(X[j])
    end
    y = purify(exp.(df[:target]))
    return X, y
end
datanow=load_ames

"""Load a well-known public classification dataset with nominal features."""
function load_iris()
    df = CSV.read(joinpath(srcdir, "data", "iris.csv"),
                  categorical=false, allowmissing=:none)
    features = filter(names(df)) do f
        f != :target
    end
    X = df[features]
    for j in 1:size(X, 2)
        X[j] = purify(X[j])
    end
    y = purify(df[:target])
    return X, y
end

function keys_ordered_by_values(d::Dict{T,S}) where {T, S<:Real}

    items = collect(d) # 1d array containing the (key, value) pairs
    sort!(items, by=pair->pair[2], alg=QuickSort)

    return T[pair[1] for pair in items]

end

"""
    bootstrap_resample_of_mean(v; n=10^6)

Returns a vector of `n` estimates of the of the distribution
generating the samples in `v` from `n` bootstrap resamples of `v`.

"""
function bootstrap_resample_of_mean(v; n=10^6)

    n_samples = length(v)
    mu = mean(v)

    simulated_means = Array{Float64}(undef, n)

    for i in 1:n
        pseudo_sample = sample(v, n_samples, replace=true)
        simulated_means[i]=mean(pseudo_sample)
    end
    return simulated_means
end

"""
    bootstrap_resample_of_median(v; n=10^6)

Returns a vector of `n` estimates of the median of the distribution
generating the samples in `v` from `n` bootstrap resamples of `v`.

"""
function bootstrap_resample_of_median(v; n=10^6)

    n_samples = length(v)
    mu = median(v)

    simulated_medians = Array{Float64}(n)

    for i in 1:n
        pseudo_sample = sample(v, n_samples, replace=true)
        simulated_medians[i]=median(pseudo_sample)
    end
    return simulated_medians
end

# """
#     get_meta(df)

# Return very short description of an `AbstractDataFrame`, `df`.

# """
# function get_meta(df::AbstractDataFrame)
#     return DataFrame(feature=names(df),
#                      eltype=eltypes(df),
#                      n_missing=[countmissing(df[j]) for j in 1:size(df, 2)])
# end

""" Return a dictionary of values keyed on the fields of specified
object, suitable for display"""
function params(object::BaseType)
    value_given_parameter  = Dict{Symbol,Any}()
    for name in fieldnames(typeof(object))
        if isdefined(object, name)
            value_given_parameter[name] = getfield(object,name)
        else
            value_given_parameter[name] = "<undefined>"
        end
    end
    return value_given_parameter
end

# function name(T::Type)
#     if isa(T, Union)
#         types = [getfield(T, f) for f in fieldnames(T)]
#         if isempty(types)
#             return ""
#         else
#             return string("Union{", types[1], types[2:end]..., "}")
#         end
#     else
#         return T.name.name
#     end
# end

# """ Extract type parameters of the type of an object."""
# type_parameters(object) = typeof(object).parameters

## ABSTRACT MODEL TYPES

abstract type Model <: BaseType end

# An abstract type for transformation models:
abstract type Transformer <: Model end

# Supervised model types are collected together according to their
# corresponding predictor type, `P`. For example, the model type
# `ConstantRegressor` (for predicting the historical target average on
# all input patterns) has predictor type `Float64`. One reason for
# this is that ensemble methods involve ensembles of predictors whose
# type should be known ahead of time.
abstract type SupervisedModel{P} <: Model end
abstract type Regressor{P} <: SupervisedModel{P} end
abstract type Classifier{P} <: SupervisedModel{P} end

## LOSS AND LOW-LEVEL ERROR FUNCTIONS

function rms(y, yhat, rows)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in rows
        dev = y[i] - yhat[i]
        ret += dev*dev
    end
    return sqrt(ret/length(rows))
end

function rms(y, yhat)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in eachindex(y)
        dev = y[i] - yhat[i]
        ret += dev*dev
    end
    return sqrt(ret/length(y))
end

function rmsl(y, yhat, rows)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in rows
        dev = log(y[i]) - log(yhat[i])
        ret += dev*dev
    end
    return sqrt(ret/length(rows))
end

function rmsl(y, yhat)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in eachindex(y)
        dev = log(y[i]) - log(yhat[i])
        ret += dev*dev
    end
    return sqrt(ret/length(y))
end

function rmslp1(y, yhat, rows)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in rows
        dev = log(y[i] + 1) - log(yhat[i] + 1)
        ret += dev*dev
    end
    return sqrt(ret/length(rows))
end

function rmslp1(y, yhat)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    for i in eachindex(y)
        dev = log(y[i] + 1) - log(yhat[i] + 1)
        ret += dev*dev
    end
    return sqrt(ret/length(y))
end

""" Root mean squared percentage loss """
function rmsp(y, yhat, rows)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    count = 0
    for i in rows
        if y[i] != 0.0
            dev = (y[i] - yhat[i])/y[i]
            ret += dev*dev
            count += 1
        end
    end
    return sqrt(ret/count)
end

function rmsp(y, yhat)
    length(y) == length(yhat) || throw(DimensionMismatch())
    ret = 0.0
    count = 0
    for i in eachindex(y)
        if y[i] != 0.0
            dev = (y[i] - yhat[i])/y[i]
            ret += dev*dev
            count += 1
        end
    end
    return sqrt(ret/count)
end

# function auc(truelabel::L) where L
#     _auc(y::AbstractVector{L}, yhat::AbstractVector{T}) where T<:Real =
#         ROC.AUC(ROC.roc(yhat, y, truelabel))
#     return _auc
# end

function err(rgs::SupervisedModel, predictor, X, y, rows,
             parallel, verbosity, loss::Function=rms)
    return loss(view(y, rows), predict(rgs, predictor, X, rows, parallel, verbosity))
end

function err(rgs::SupervisedModel, predictor, X, y,
             parallel, verbosity, loss::Function=rms)
    return loss(y, predict(rgs, predictor, X, parallel, verbosity))
end

# TODO: More classifier loss and metrics

## MACHINES

abstract type Machine <: BaseType end

## MACHINES FOR TRANSFORMING

struct TransformerMachine{T<:Transformer} <: Machine
    transformer::T
    scheme

    function TransformerMachine{T}(transformer::T, scheme, ::Colon) where T <: Transformer
        return new{T}(transformer, scheme)
    end
end

function TransformerMachine(transformer::T, scheme, ::Colon) where T <: Transformer
    return TransformerMachine{T}(transformer, scheme, :)
end

# usual contructor:
function TransformerMachine(transformer::T, X;
                               parallel::Bool=true,
                               verbosity::Int=1, args...) where T <: Transformer
    scheme = fit(transformer, X, parallel, verbosity; args...)
    return TransformerMachine(transformer, scheme, :)
end

Machine(transformer::Transformer, X; args...) =
    TransformerMachine(transformer, X; args...)

function transform(mach::TransformerMachine, X)
    return transform(mach.transformer, mach.scheme, X)
end

function inverse_transform(mach::TransformerMachine, X; parallel=true, verbosity=1)
    return inverse_transform(mach.transformer, mach.scheme, X)
end

## COMMONLY USED TRANSFORMERS

# a null transformer, useful as default keyword argument
# (has no fit, transform or inverse_transform methods):
struct EmptyTransformer <: Transformer end

# no transformer is empty except the null transformer:
Base.isempty(transformer::Transformer) = false
Base.isempty(transformer::EmptyTransformer) = true

# For (a) remembering the features used in `fit` (calibration), in
# order presented, and selecting only those on tranforming new data
# frames; or (b) selecting only feature labels specified in the
# tranformer:
mutable struct FeatureSelector <: Transformer
    features::Vector{Symbol} # empty means do (a) above
end
FeatureSelector(;features=Symbol[]) = FeatureSelector(features)
function fit(transformer::FeatureSelector, X, parallel, verbosity)
    namesX = names(X)
    issubset(Set(transformer.features), Set(namesX)) ||
        @error "Attempting to select non-existent feature(s)."
    if isempty(transformer.features)
        return namesX
    else
        return transformer.features
    end
end
function transform(transformer::FeatureSelector, features, X)
    issubset(Set(features), Set(names(X))) ||
        @error "Supplied frame does not admit previously selected features."
    return X[features]
end

# identity transformations:
struct IdentityTransformer <: Transformer end
fit(transformer::IdentityTransformer, y, parallel, verbosity) = nothing
transform(transformer::IdentityTransformer, scheme, y) = y
inverse_transform(transformer::IdentityTransformer, scheme, y) = y

"""
    struct RowsTransformer <: Transformer

A special transformer mapping `AbstractVector{Int}`s to
`AbstractVector{Int}`s, to deal with dropping rows from an input
`DataFrame`.

Suppose `df` is a `DataFrame` and `bad` is a vector of indices for
rows that have been dropped from `df` to form `df_safe`. Given a
vector of indices `rows`, I want `df[good_rows,:]` where `good_rows`
is `rows` with any elements of `bad` removed. However, I only have
access to `df_safe` (although I know `size(df, 1)`).  An appropriate
transformation of `rows` gives a vector of indices to use on
`df_safe`:

    mach = Machine(RowTransformer(size(df,1)), bad)
    df[good_rows,:] == df_safe[transform(mach, rows),:] # true

Analogous statements hold for array-like collections apart from
`DataFrames`.

"""
struct RowsTransformer <: Transformer
    original_num_rows::Int
end

function fit(transformer::RowsTransformer, missing_indices, parallel, verbosity)
    scheme = Array{Union{Int,Missing}}(undef, transformer.original_num_rows)
    counter = 1
    for i in 1:transformer.original_num_rows
        if i in missing_indices
            scheme[i] = missing
        else
            scheme[i] = counter
            counter += 1
        end
    end
    return scheme
end

transform(transformer::RowsTransformer, scheme, rows) =
    collect(skipmissing(view(scheme, rows)))

## MACHINES FOR SUPERVISED LEARNING

"""
## `mutable struct SupervisedMachine`

### Keyword options

    drop_unseen=false

All calls to `fit!`, `err`, `cv` and `learning_curve` on the machine
will ignore rows with categorical features taking on values not seen
in training. Note: The `predict` method cannot be safely called on
data containing such rows.

    features = Symbol[]

Only specified features are used. Empty means *all* features are
used. The features used can also be controlled by `transformer_X`.

"""
mutable struct SupervisedMachine{P, M <: SupervisedModel{P}} <: Machine

    model::M
    transformer_X::Transformer
    transformer_y::Transformer
    scheme_X
    scheme_y
    rows_with_unseen::AbstractVector{Int}
    rows_transformer_machine::TransformerMachine{RowsTransformer}
    n_iter::Int
    Xt
    yt
    predictor::P
    report
    cache

    function SupervisedMachine{P, M}(
        model::M,
        X::AbstractDataFrame,
        y::AbstractVector,
        train_rows::AbstractVector{Int}; # for defining data transformations
        features=Symbol[],
        transformer_X=EmptyTransformer(),
        transformer_y=EmptyTransformer(),
        drop_unseen::Bool=false,
        verbosity::Int=1) where {P, M <: SupervisedModel{P}}

        # check dimension match:
        size(X,1) == length(y) || throw(DimensionMismatch())

        # check valid `features`; if none provided, take to be all
        if isempty(features)
            features = names(X)
        end
        allunique(features) || @error "Duplicate features."
        issubset(Set(features), Set(names(X))) || @error "Invalid feature vector."

        # bind `X` to view of `X` containing only the features
        # specified, in the order specified (does not generate copy):
        X = X[features]

        # check for missing data and report eltypes:
        if verbosity > 0
            kind_given_feature = Dict{Symbol,String}()
            @info "Element types of input features before transformation:"
        end
        for field in names(X)
            T = eltype(X[field])
            real_categorical = (T <: AbstractFloat ? "continuous" : "categorical")
            verbosity > 0 ? kind_given_feature[field] = real_categorical : nothing
            if ismissingtype(T) && verbosity > -1
                    @warn ":$field has a missing-element type. Treated as categorical."
            end
        end
        if verbosity > 0
            show(stdout, MIME("text/plain"), kind_given_feature)
            println()
        end

        # report size of data used for transformations
        percent_train = round(Int, 1000*length(train_rows)/length(y))/10
        verbosity < 0 || @info "$percent_train% of data used to compute "*
                              "transformation parameters."

        # assign transformers if not provided:
        if isempty(transformer_X)
            transformer_X = default_transformer_X(model)
        end
        if isempty(transformer_y)
            transformer_y = default_transformer_y(model)
        end

        # if necessary, determine rows with categorical levels not
        # seen in X[train_rows, :] and build corresponding rows
        # transformer:
        if drop_unseen

            # get test_rows:
            test_rows = filter(eachindex(y)) do i
                !(i in train_rows)
            end

            seen, unseen = split_seen_unseen(X, train_rows, test_rows)
            if verbosity > -1
                if length(unseen) == 0
                    @warn "All remaining rows of the inputs data have categorical "*
                         "features taking values not seen during data transformation."
                else
                    bad_percentage = round(Int, 1000*length(unseen)/length(test_rows))/10
                    @warn "$bad_percentage% of the remaining rows of input data "*
                         "(recorded in the attribute `rows_with_unseen`) "*
                         "contain "*
                         "patterns for which some categorical feature "*
                         "takes on values not seen during data "*
                         "transformation. These will be ignored in "*
                         "calls to err(), cv() and  learning_curve(). "*
                         "However, you will be unable to safely "*
                         "call predict() on such data."
                end
            end
        else
            unseen = Int[]
        end

        rows_transformer_machine = Machine(RowsTransformer(length(y)), unseen)

        mach = new{P, M}(model::M)
        mach.transformer_X = transformer_X
        mach.transformer_y = transformer_y
        mach.scheme_X = fit(mach.transformer_X, view(X, train_rows),
                            true, verbosity - 1)
        mach.scheme_y = fit(mach.transformer_y, y[train_rows], verbosity - 1, 1)
        mach.rows_with_unseen = unseen
        mach.rows_transformer_machine = rows_transformer_machine

        # transform the data, dumping any untransformable data if
        # drop_unseen = true, while throwing an exception otherwise:
        if isempty(unseen)
            try
                mach.Xt = transform(mach.transformer_X, mach.scheme_X, X)
                mach.yt = transform(mach.transformer_y, mach.scheme_y, y)
            catch exception
                if isa(exception, KeyError)
                    @error "KeyError: key $(exception.key) not found. Problably "*
                        "a categorical feature takes on values "*
                        "not encountered in the rows provided for "*
                        "computing data transformations. Try calling "*
                        "machine constructor with `drop_unseen=true`."
                else
                    throw(exception)
                end
            end
        else
            all_good_rows = filter(eachindex(y)) do i
                !(i in unseen)
            end
            mach.Xt = transform(mach.transformer_X, mach.scheme_X, view(X, all_good_rows))
            mach.yt = transform(mach.transformer_y, mach.scheme_y, view(y, all_good_rows))
        end

        mach.n_iter = 0
        mach.report = Dict{Symbol,Any}()

        return mach
    end

end

function SupervisedMachine(model::M, X, y,
                           train_rows; args...) where {P, M <: SupervisedModel{P}}
    return SupervisedMachine{P, M}(model, X, y, train_rows; args...)
end

Machine(model::SupervisedModel, X, y, train_rows; args...) =
    SupervisedMachine(model, X, y, train_rows; args...)

function fit!(mach::SupervisedMachine, rows;
              add=false, verbosity=1, parallel=true, args...)
    verbosity < 0 || softwarn(clean!(mach.model))

    # transform `rows` to account for rows that may have been dropped
    # during transformation:
    if !isempty(mach.rows_with_unseen)
        rows = transform(mach.rows_transformer_machine, rows)
    end

    if !add
        mach.n_iter = 0
    end
    if  mach.n_iter == 0
        mach.cache = setup(mach.model, mach.Xt, mach.yt, rows, mach.scheme_X,
                           parallel, verbosity)
    end
    mach.predictor, report, mach.cache =
        fit(mach.model, mach.cache, add, parallel, verbosity; args...)
    merge!(mach.report, report)
    if isdefined(mach.model, :n)
        mach.n_iter += mach.model.n
    else
        mach.n_iter = 1
    end
    return mach
end

# When rows are not specified, the cache is not recalculated. Ie,
# `setup` is skipped:
function fit!(mach::SupervisedMachine;
              add=false, verbosity=1, parallel=true, args...)
    verbosity < 0 || softwarn(clean!(mach.model))
    if !isdefined(mach, :cache)
        @error "You must specify training rows in the first call to fit!\n"*
              "E.g., fit!(mach, train_rows)."
    end
    if !add
        mach.n_iter = 0
    end
    mach.predictor, report, mach.cache =
        fit(mach.model, mach.cache, add, parallel, verbosity; args...)
    merge!(mach.report, report)
    if isdefined(mach.model, :n)
        mach.n_iter += mach.model.n
    else
        mach.n_iter = 1
    end
    return mach
end

function predict(mach::SupervisedMachine, X, rows; parallel=true, verbosity=1)
    mach.n_iter > 0 || @error string(mach, " has not been fitted.")
    Xt = transform(mach.transformer_X, mach.scheme_X, view(X, rows))
    yt = predict(mach.model, mach.predictor, Xt, parallel, verbosity)
    return inverse_transform(mach.transformer_y, mach.scheme_y, yt)
end

function predict(mach::SupervisedMachine, X; parallel=true, verbosity=1)
    mach.n_iter > 0 || @error string(mach, " has not been fitted.")
    Xt = transform(mach.transformer_X, mach.scheme_X, X)
    yt = predict(mach.model, mach.predictor, Xt, parallel, verbosity)
    return inverse_transform(mach.transformer_y, mach.scheme_y, yt)
end

function err(mach::SupervisedMachine, test_rows;
             loss=rms, parallel=false, verbosity=1, raw=false)

    mach.n_iter > 0 || @error "Attempting to predict using untrained machine."

    !raw || verbosity < 0 || @warn "Reporting errors for *transformed* target. "*
                                  "Use `raw=false` to report true errors."

    # transform `test_rows` to account for rows that may have been
    # dropped during transformation:
    if !isempty(mach.rows_with_unseen)
        verbosity < 1 ||
            @info "Ignoring rows with categorical data values "*
                 "not seen in transformation."
        test_rows = transform(mach.rows_transformer_machine, test_rows)
    end

    # transformed version of target predictions:
    raw_predictions = predict(mach.model, mach.predictor, mach.Xt, test_rows,
                            parallel, verbosity)

    if raw # return error on *transformed* target, which is faster
        return loss(view(mach.yt, test_rows), raw_predictions)
    else  # return error for untransformed target
        return loss(inverse_transform(mach.transformer_y, mach.scheme_y,
                                      view(mach.yt, test_rows)),
                    inverse_transform(mach.transformer_y, mach.scheme_y,
                                      raw_predictions))
    end
end

## `SupervisedModel`  fall-back methods

# By default all supervised model types can be ensembled.
supports_ensembling(model::SupervisedModel) = true

# default default_transformers:
default_transformer_X(model::SupervisedModel) = FeatureSelector()
default_transformer_y(model::SupervisedModel) = IdentityTransformer()

# for enforcing model parameter invariants:
clean!(model::SupervisedModel) = "" # no checks, emtpy message

# To allow for an extra `rows` argument in `setup` and `predict`:

# fallbacks:
setup(model::SupervisedModel, Xt, yt, rows, scheme_X, parallel, verbosity) =
    setup(model, Xt[rows,:], view(yt, rows), scheme_X, parallel, verbosity)
predict(model::SupervisedModel, predictor, Xt, rows, parallel, verbosity) =
    predict(model, predictor, Xt[rows,:], parallel, verbosity)

## CONSTANT-PREDICTING REGRESSOR (A `REGRESSOR` EXAMPLE)

# Note: To test iterative methods, we give the following simple regressor
# model a "bogus" field `n` for counting the number of iterations (which
# make no difference to predictions):

mutable struct ConstantRegressor <: Regressor{Float64}
    n::Int
end

ConstantRegressor() = ConstantRegressor(1)

function setup(rgs::ConstantRegressor, X, y, scheme_X, parallel, verbosity)
    return mean(y)
end

function fit(rgs::ConstantRegressor, cache, add, parallel, verbosity)
    predictor = cache
    report = Dict{Symbol, Any}()
    report[:mean] = predictor
    return predictor, report, cache
end

function predict(rgs::ConstantRegressor, predictor, X, parallel, verbosity)
    return  Float64[predictor for i in 1:size(X,1)]
end

## VALIDATION TOOLS

"""
    partition(rows::AbstractVector{Int}, fractions...)

Splits the vector `rows` into a tuple of vectors whose lengths are
given by the corresponding `fractions` of `length(rows)`. The last
fraction is not provided, as it is inferred from the preceding
ones. So, for example,

    julia> partition(1:1000, 0.2, 0.7)
    (1:200, 201:900, 901:1000)

"""
function partition(rows::AbstractVector{Int}, fractions...)
    rows = collect(rows)
    rowss = []
    if sum(fractions) >= 1
        throw(DomainError)
    end
    n_patterns = length(rows)
    first = 1
    for p in fractions
        n = round(Int, p*n_patterns)
        n == 0 ? (@warn "A split has only one element"; n = 1) : nothing
        push!(rowss, rows[first:(first + n - 1)])
        first = first + n
    end
    if first > n_patterns
        @warn "Last vector in the split has only one element."
        first = n_patterns
    end
    push!(rowss, rows[first:n_patterns])
    return tuple(rowss...)
end

"""
## `function learning_curve(mach::SupervisedMachine, train_rows, test_rows,
##                      range; restart=true, loss=rms, raw=false, parallel=true,
##                      verbosity=1, fit_args...)`

    u,v = learning_curve(mach, test_rows, 1:10:200)
    plot(u, v)

Assming, say, `Plots` is installed, the above produces a plot of the
RMS error for the machine `mach`, on the test data with rows
`test_rows`, against the number of iterations `n` of the algorithm it
implements (assumed to be iterative). Here `n` ranges over `1:10:200`
and training is performed using `train_rows`. For parallization, the
value of the optional keyword `parallel` is passed to each call to
`fit`, along with any other keyword arguments `fit_args` that `fit`
supports.

"""
function learning_curve(mach::SupervisedMachine, train_rows, test_rows,
                        range; restart=true, loss=rms, raw=false, parallel=true,
                        verbosity=1, fit_args...)

    isdefined(mach.model, :n) || @error "$(mach.model) does not support iteration."

    # save to be reset at end:
    old_n = mach.model.n

    !raw || verbosity < 0 ||
        @warn "Reporting errors for *transformed* target. Use `raw=false` "*
             " to report true errors."

    # transform `train_rows` and `test_rows` to account for rows that
    # may have been dropped during transformation:
    if !isempty(mach.rows_with_unseen)
        verbosity < 1 ||
            @info "Ignoring rows with categorical data values "*
                 "not seen in transformation."
        train_rows = transform(mach.rows_transformer_machine, train_rows)
    end

    range = collect(range)
    sort!(range)

    if restart
        mach.n_iter = 0
        mach.cache = setup(mach.model, mach.Xt, mach.yt, train_rows, mach.scheme_X,
                           parallel, verbosity - 1)
    end

    n_iter_list = Float64[]
    errors = Float64[]

    filter!(range) do x
        x > mach.n_iter
    end

    report = Dict{Symbol,Any}()
    while !isempty(range)
        verbosity < 1 || print("\rNext iteration number: ", range[1])
        # set number of iterations for `fit` call:
        mach.model.n = range[1] - mach.n_iter
        if mach.n_iter == 0 # then use add=false
            mach.predictor, report, mach.cache =
                fit(mach.model, mach.cache, false, parallel, verbosity - 1; fit_args...)
        else # use add=true
            mach.predictor, report, mach.cache =
                fit(mach.model, mach.cache, true, parallel, verbosity - 1; fit_args...)
        end
        mach.n_iter += mach.model.n
        push!(n_iter_list, mach.n_iter)
        push!(errors, err(mach, test_rows, raw=raw, loss=loss, verbosity=verbosity - 1))
        filter!(range) do x
            x > mach.n_iter
        end
    end

    verbosity < 1 || println("\nLearning curve added to machine report.")

    report[:learning_curve] = (n_iter_list, errors)
    merge!(mach.report, report)

    mach.model.n = old_n

    return n_iter_list, errors

end

"""
    cv(mach::SupervisedMachine, rows;
       n_folds=9, loss=rms, parallel=true, verbosity=1, raw=false, randomize=false)

Return a list of cross-validation losses for patterns with row indices
in `rows`, an iterator that is initially randomized when an optional
parameter `randomize` is set to `true`.

"""
function cv(mach::SupervisedMachine, rows; n_folds=9, loss=rms,
             parallel=true, verbosity=1, raw=false, randomize=false)

    !raw || verbosity < 0 ||
        @warn "Reporting errors for *transformed* target. Use `raw=false` "*
             " to report true errors."

    # transform `rows` to account for rows that may have been dropped
    # during transformation:
    if !isempty(mach.rows_with_unseen)
        verbosity < 1 ||
            @info "Ignoring rows with categorical data values "*
                 "not seen in transformation."
        rows = transform(mach.rows_transformer_machine, rows)
    end

    n_samples = length(rows)
    if randomize
         rows = sample(collect(rows), n_samples, replace=false)
    end
    k = floor(Int,n_samples/n_folds)

    # function to return the error for the fold `rows[f:s]`:
    function get_error(f, s)
        test_rows = rows[f:s]
        train_rows = vcat(rows[1:(f - 1)], rows[(s + 1):end])
        fit!(mach, train_rows; parallel=false, verbosity=0)
        return err(mach, test_rows; loss=loss,
                   parallel=false, verbosity=verbosity - 1,
                   raw=raw)
    end

    firsts = 1:k:((n_folds - 1)*k + 1) # itr of first test_rows index
    seconds = k:k:(n_folds*k)          # itr of ending test_rows index

    if parallel && nworkers() > 1
        if verbosity > 0
            println("Distributing cross-validation computation "*
                    "among $(nworkers()) workers.")
        end
        return @distributed vcat for n in 1:n_folds
            Float64[get_error(firsts[n], seconds[n])]
        end
    else
        errors = Array{Float64}(undef, n_folds)
        for n in 1:n_folds
            verbosity < 1 || print("\rfold: $n")
            errors[n] = get_error(firsts[n], seconds[n])
        end
        verbosity < 1 || println()
        return errors
    end

end

macro colon(p)
    Expr(:quote, p)
end

"""
## `@curve`

The code,

    @curve var range code

evaluates `code`, replacing appearances of `var` therein with each
value in `range`. The range and corresponding evaluations are returned
as a tuple of arrays. For example,

    @curve  x 1:3 (x^2 + 1)

evaluates to

    ([1,2,3], [2, 5, 10])

This is convenient for plotting functions using, eg, the `Plots` package:

    plot(@curve x 1:3 (x^2 + 1))

A macro `@pcurve` parallelizes the same behaviour.  A two-variable
implementation is also available, operating as in the following
example:

    julia> @curve x [1,2,3] y [7,8] (x + y)
    ([1,2,3],[7 8],[8.0 9.0; 9.0 10.0; 10.0 11.0])

    julia> ans[3]
    3×2 Array{Float64,2}:
      8.0   9.0
      9.0  10.0
     10.0  11.0

N.B. The second range is returned as a *row* vector for consistency
with the output matrix. This is also helpful when plotting, as in:

    julia> u1, u2, A = @curve x range(0, stop=1, length=100) α [1,2,3] x^α
    julia> u2 = map(u2) do α "α = "*string(α) end
    julia> plot(u1, A, label=u2)

which generates three superimposed plots - of the functions x, x^2 and x^3 - each
labels with the exponents α = 1, 2, 3 in the legend.

"""
macro curve(var1, range, code)
    quote
        output = []
        N = length($(esc(range)))
        for i in eachindex($(esc(range)))
            local $(esc(var1)) = $(esc(range))[i]
            print((@colon $var1), "=", $(esc(var1)), "                   \r")
            flush(stdout)
            # print(i,"\r"); flush(stdout)
            push!(output, $(esc(code)))
        end
        collect($(esc(range))), [x for x in output]
    end
end

macro curve(var1, range1, var2, range2, code)
    quote
        output = Array{Float64}(undef, length($(esc(range1))), length($(esc(range2))))
        for i1 in eachindex($(esc(range1)))
            local $(esc(var1)) = $(esc(range1))[i1]
            for i2 in eachindex($(esc(range2)))
                local $(esc(var2)) = $(esc(range2))[i2]
                # @dbg $(esc(var1)) $(esc(var2))
                print((@colon $var1), "=", $(esc(var1)), " ")
                print((@colon $var2), "=", $(esc(var2)), "                    \r")
                flush(stdout)
                output[i1,i2] = $(esc(code))
            end
        end
        collect($(esc(range1))), permutedims(collect($(esc(range2)))), output
    end
end

macro pcurve(var1, range, code)
    quote
        N = length($(esc(range)))
        pairs = @distributed vcat for i in eachindex($(esc(range)))
            local $(esc(var1)) = $(esc(range))[i]
#            print((@colon $var1), "=", $(esc(var1)), "                    \r")
#            flush(stdout)
#            print(i,"\r"); flush(stdout)
            [( $(esc(range))[i], $(esc(code)) )]
        end
        sort!(pairs, by=first)
        collect(map(first,pairs)), collect(map(last, pairs))
    end
end

"""
## `function compete(e0, e1; alpha=0.05)`

Given paired samples `e0` and `e1`, we test the
null-hypothesis that the underlying distributions have the same mean,
using the significance level `alpha`. Normality of the underlying
distributions is assumed and a two-sided t-test applied.

### Return value

`(c, p)` where `p` is the p-value for the test and `c` is the competition outcome:

- '0' if the null is rejected and `e0` has the smaller mean (M1 "wins")
- '1' if the null is rejected and `e1` has the smaller mean (M0 "wins")
- 'D' if the null is accepted ("draw" because `p > alpha`)

"""
function compete(e0, e1; alpha=0.05)
    t = OneSampleTTest(e1-e0)
    p = pvalue(t)
    if p > alpha
        return 'D', p
    end
    if mean(e0) < mean(e1)
        return '0', p
    end
    return '1', p
end

"""
## `function split_seen_unseen`

    split_seen_unseen(v::Vector, train_indices, test_indices)`

Detects on which indices in `test_indices` the `AbstractVector`, `v`,
takes on values not occuring in the `train_indices` part of the
vector. Returns a pair of integer vectors `(S, U)` where `S` are the
seen indices and `U` the unseen ones.

    split_seen_unseen(df::AbstractDataFrame, train_rows, test_rows)

Detects the rows in `test_rows` for which any of the categorical
features of `df` takes on a value not observed in
`train_rows`. Returns a pair of integer vectors `(S, U)` where `S` are
the unseen rows and `U` the unseen ones.

"""
function split_seen_unseen(v::AbstractVector, train_rows, test_rows)
    train_values = Set(v[train_rows])
    test_values  = Set(v[test_rows])
    unseen_values = setdiff(test_values, train_values)

    is_unseen = map(test_rows) do i
        v[i] in unseen_values
    end

    return test_rows[.!is_unseen], test_rows[is_unseen]

end

function split_seen_unseen(df::AbstractDataFrame, train_rows, test_rows)

    # collect all rows with unseen categorical values together:
    unseen_rows = Int[]
    for j in 1:size(df, 2)
        if !(eltype(df[j]) <: AbstractFloat)
            seen, unseen = split_seen_unseen(df[j], train_rows, test_rows)
            append!(unseen_rows, unseen)
        end
    end
    unseen_rows = unique(unseen_rows)

    ordered_unseen = Int[]
    ordered_seen = Int[]

    for i in test_rows
        if i in unseen_rows
            push!(ordered_unseen, i)
        else
            push!(ordered_seen, i)
        end
    end

    return ordered_seen, ordered_unseen

end

## RECIPES FOR USE WITH Plots.jl

"""
    PlotableDict(d)

Returns a wrapped version of the dictionary `d` (assumed to have
`Real` values) to enable plotting with calls to `Plot.plot`, etc
(provided by a `Plots` recipe). In particular,

    plot(PlotableDict(d), ordered_by_keys=true)

plots the bars in the lexographic order of the keys. Otherwise, the
bars are ordered by the values corresponding to to each key.

"""
struct PlotableDict{KeyType,ValueType<:Real}
    dict::Dict{KeyType,ValueType}
end

@recipe function dummy(::Type{PlotableDict{KeyType,ValueType}},
                       pd::PlotableDict{KeyType,ValueType};
                       ordered_by_keys=false) where {KeyType,ValueType<:Real}

    seriestype := :bar

    dict_to_be_plotted = pd.dict
    x = String[]
    y = ValueType[]

    if ordered_by_keys
        kys = sort(collect(keys(dict_to_be_plotted)))
    else
        kys = keys_ordered_by_values(dict_to_be_plotted)
    end

    for k in kys
        push!(x, string(k))
        push!(y, dict_to_be_plotted[k])
    end
    offset = 0.05*maximum(abs.(y))
    xticks := Float64[]
    annotations := [(Float64(i), offset, x[i]) for i in eachindex(x)]
    y

end

mutable struct BootstrapHistogramOfMean
        args
end

"""
    bootstrap_histogram_of_mean(v; n=1e6, plotting_kws...)

Create a bootstrap sample of size `n` from `v` and plot the
corresponding histogram.

"""
bootstrap_histogram_of_mean(args...; kw...) = begin
            RecipesBase.plot(BootstrapHistogramOfMean(args); kw...)
end

bootstrap_histogram_of_mean!(args...; kw...) = begin
    RecipesBase.plot!(BootstrapHistogramOfMean(args); kw...)
end

@recipe function dummy(h::BootstrapHistogramOfMean; n=10^6)
    length(h.args) == 1 || typeof(h.args) <: AbstractVector ||
        @error "A BootstrapHistogramOfMean should be given one vector. "*
    "Got: $(typeof(h.args))"
    v = h.args[1]
    bootstrap = bootstrap_resample_of_mean(v; n=n)
    @series begin
        seriestype := :histogram
        alpha --> 0.5
        normalized := true
        bins --> 50
        bootstrap
    end
    if @isdefined StatPlots
        @series begin
            seriestype := :density
            label := ""
            linewidth --> 2
            color --> :black
            bootstrap
        end
    else
        @info "For denisty approximation in a bootstrap_histogram_of_mean, "*
        "import StatPlots."
    end
end

mutable struct BootstrapHistogramOfMedian
        args
end

"""
    bootstrap_histogram_of_median(v; n=1e6, plotting_kws...)

Create a bootstrap sample of size `n` from `v` and plot the
corresponding histogram.

"""
bootstrap_histogram_of_median(args...; kw...) = begin
            RecipesBase.plot(BootstrapHistogramOfMean(args); kw...)
end

bootstrap_histogram_of_median!(args...; kw...) = begin
    RecipesBase.plot!(BootstrapHistogramOfMean(args); kw...)
end

@recipe function dummy(h::BootstrapHistogramOfMedian; n=10^6)
    length(h.args) == 1 || typeof(h.args) <: AbstractVector ||
        @error "A BootstrapHistogramOfMedian should be given one vector. "*
    "Got: $(typeof(h.args))"
    v = h.args[1]
    bootstrap = bootstrap_resample_of_median(v; n=n)
    @series begin
        seriestype := :histogram
        alpha --> 0.5
        normalized := true
        bins --> 50
        bootstrap
    end
    if @isdefined StatPlots
        @series begin
            seriestype := :density
            label := ""
            linewidth --> 2
            color --> :black
            bootstrap
        end
    else
        @info "For denisty approximation in a bootstrap_histogram_of_median, "*
        "import StatPlots."
    end
end

end # module


text:
using ValidatedNumerics
using FactCheck

facts("Interval loop tests") do
    i = 1

    @fact Interval(i,i).lo --> 1
    @fact @interval(i).lo --> 1

    for i in 1:10
        a = @interval(i)
        @fact a.lo --> i
    end

end

## Calculate pi by summing 1/i^2 to give pi^2/6:

setprecision(Interval, 53)

function calc_pi1(N)
    S1 = @interval(0)

    for i in 1:N
        S1 += @interval(1/i^2)
    end
    S1 += @interval(1/(N+1), 1/N)

    sqrt(6*S1)
end

function calc_pi2(N)
    S2 = @interval(0)

    for i in 1:N
        S2 += 1/i^2
    end
    S2 += @interval(1/(N+1), 1/N)

    sqrt(6*S2)
end

function calc_pi3(N)
    S3 = @floatinterval(0)

    for i in 1:N
        S3 += 1/i^2
    end
    S3 += @floatinterval(1/(N+1), 1/N)

    sqrt(6*S3)
end

function calc_pi4(N)
    S4 = @floatinterval(0)
    II = @floatinterval(1)

    for i in N:-1:1
        S4 += II / (i^2)
    end
    S4 += II / @floatinterval(N, N+1)

    sqrt(6*S4)
end

function calc_pi5(N)
    S5 = @floatinterval(0)
    II = @floatinterval(1)

    for i in N:-1:1
        S5 += 1 // (i^2)
    end
    S5 += 1 / @floatinterval(N, N+1)

    sqrt(6*S5)
end

facts("Pi tests") do

    big_pi = setprecision(256) do
        big(pi)
    end

    N = 10000
    pi1 = calc_pi1(N)
    pi2 = calc_pi2(N)
    pi3 = calc_pi3(N)
    pi4 = calc_pi4(N)
    pi5 = calc_pi5(N)

    @fact big_pi ∈ pi1 --> true
    @fact big_pi ∈ pi2 --> true
    @fact big_pi ∈ pi3 --> true
    @fact big_pi ∈ pi4 --> true
    @fact big_pi ∈ pi5 --> true

    @pending pi1 == pi2 --> true
    @pending pi2 == pi3 --> true

end


text:
using Winston

## We create spline bases on for the function space [5,15] -> Real with evenly-spaced nodes, and
## apply them to the grid x
x = linspace(5,15,101)
X1 = linearSplineBFE(x, (5,15), 11)
X2 = quadraticSplineBFE(x, (5,15), 11)
X3 = cubicSplineBFE(x, (5,15), 11)

p = FramedPlot()
map(i -> add(p, Curve(x, X1[:,i], "color", "black")), 1:3)
Winston.display(p)

p = FramedPlot()
map(i -> add(p, Curve(x, X2[:,i], "color", "red")), 1:5)
Winston.display(p)

p = FramedPlot()
map(i -> add(p, Curve(x, X3[:,i], "color", "blue")), 1:7)
Winston.display(p)

## Now we use a pre-specified set of nodes
u = [8.05,11.34,12.31,14.53, 14.75]
X1 = splineBFE(1, x, u)
X2 = splineBFE(2, x, u)

p = FramedPlot()
map(i -> add(p, Curve(x, X1[:,i], "color", "black")), 1:3)
Winston.display(p)

p = FramedPlot()
map(i -> add(p, Curve(x, X2[:,i], "color", "black")), 1:2)
Winston.display(p)


text:
@inline BioSequences.encoded_data_eltype(::Type{<:Kmer}) = UInt64

@inline function BioSequences.extract_encoded_element(seq::Kmer, i::Integer)
    bi = BioSequences.bitindex(seq, i % UInt)
    return BioSequences.extract_encoded_element(bi, seq.data)
end

@inline Base.copy(seq::Kmer) = typeof(seq)(seq.data)

@inline encoded_data(x::Kmer) = x.data

@inline BioSequences.bitindex(seq::Kmer, i::Integer) = BioSequences.bitindex(BioSequences.BitsPerSymbol(seq), BioSequences.encoded_data_eltype(typeof(seq)), i + n_unused(seq))


text:
using Test, Random, FillArrays
using LuxurySparse
using SparseArrays
using LinearAlgebra

Random.seed!(2)

p1 = IMatrix{4}()
sp = sprand(ComplexF64, 4, 4, 0.5)
ds = rand(ComplexF64, 4, 4)
pm = PermMatrix([2, 3, 4, 1], randn(4))
v = [0.5, 0.3im, 0.2, 1.0]
dv = Diagonal(v)

@testset "basic" begin
    @test p1 == copy(p1)
    @test eltype(p1) == Bool
    @test size(p1) == (4, 4)
    @test size(p1, 1) == size(p1, 2) == 4
    @test Matrix(p1) == [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]

    p0 = similar(p1, ComplexF64)
    @test p0 !== p1
end

@testset "conversion" begin
    for mat in [p1, pm, dv]
        @test mat == SparseMatrixCSC(mat)
        @test mat == Matrix(mat)
    end
    for mat in [p1, pm, dv]
        @test mat == PermMatrix(mat)
    end
    @test Diagonal(p1) == p1

    @test SparseMatrixCSC(Diagonal(Fill(2, 4))) ≈ Diagonal(Fill(2, 4))
end

@testset "sparse" begin
    @test LuxurySparse.nnz(p1) == 4
    @test LuxurySparse.nonzeros(p1) == ones(4)
end

@testset "linalg" begin
    for op in [conj, real, transpose, copy, inv]
        @test op(p1) == Matrix(I, 4, 4)
        @test typeof(op(p1)) == typeof(p1)
    end
    @test imag(p1) == zeros(4, 4)
    @test p1' == Matrix(I, 4, 4)

    # This will be lazy evaluated in 0.7+
    @static if VERSION < v"0.7-"
        @test typeof(p1') == typeof(p1)
    end

    @test ishermitian(p1)
end

@testset "elementary" begin
    @test all(isapprox.(conj(p1), conj(Matrix(p1))))
    @test all(isapprox.(real(p1), real(Matrix(p1))))
    @test all(isapprox.(imag(p1), imag(Matrix(p1))))
end

@testset "basicmath" begin
    @test p1 * 2im == Matrix(p1) * 2im
    @test p1 / 2.0 == Matrix(p1) / 2.0
end


text:
include("MolecularDynamicsTrial.jl")            # type def
include("mdenergies.jl")                        # energies and index
include("mdplace.jl")                           # initial placement
include("mdverlet.jl")                          # Verlet algorithm
include("mdmetropolis.jl")                      # Metropolis algorithm
include("mdio.jl")                              # save and load
include("mdrelax.jl")                           # relax temperature
include("mdthermalize.jl")                      # thermalize, relax, repeat
include("extract_tenth_t_pe_vir_v.jl")          # extract data


text:
@testset "Isaac" begin

println("Basic correctness tests.")

# TODO: add tests determining what happens on failure!!!!

@testset "dirder" begin
   @info("testing `dirder`")
   for (f, df, init, randinit, name) in tests
      for i = 1:10
         x = randinit()
         w = rand(length(x)) .- 0.5; w /= norm(w)
         err = norm( df(x) * w - Isaac.dirder(x, w, f, f(x)) )
         print_tf(@test err < 1e-6)
      end
   end
   println()
end

@testset "dgmres" begin
   @info("testing dgmres")
   for (f, df, init, randinit, name) in tests
      x = randinit()
      fx = f(x)
      for tol in (1e-1, 1e-2, 1e-3)
         maxiter = 100
         u, error, total_iters = dgmres(fx, f, x, tol, maxiter)
         err = norm(Isaac.dirder(x, u, f, fx) - (-fx))
         print_tf(@test err / norm(fx) < tol)
         # TODO: This next test seems to fail quite often
         # @test abs(error[end] - err) < 1e-7
         print_tf(@test total_iters < maxiter)
      end
   end
   println()
end

@testset "nsoli" begin
   println("testing nsoli")
   for (f, df, init, randinit, name) in tests
      for i = 1:10
         x0 = randinit()
         atol = 1e-5
         maxit = 100
         x, it_hist, ierr, x_hist = nsoli(x0, f;
                                       atol = atol, rtol = 0.0, maxit = maxit)
         if norm(f(x)) > atol
            print_tf(@test ierr != 0)
            print_tf(@test length(x_hist) >= maxit)
         else
            print_tf(@test norm(f(x)) <= atol)
         end
      end
   end
   println()
end

end


text:
function determine_closed(closed_order, id_str, client, status, config; weightedAvgTranslator = 0)
    tmp = deepcopy(closed_order)
    sell_date, sell_time  = split(tmp["SELL"]["time"], "T")
    yr, M, d   = split(sell_date, "-")
    new_closed = Dict()
    new_closed["orders"] = tmp
    @info(id_str)
    sleep(1)
    tmp_vol    = 0
    tmp_pro    = 0
    timestamps = []
    order      = Dict()
    new_closed["buys_exec"] = 0
    new_closed["commission"] = 0
    new_closed["symbol"] = ""
    for (k2, v2) in tmp
        if haskey(v2, "time")
            push!(timestamps, datetime2unix(DateTime(v2["time"])))
        elseif haskey(v2, "updateTime")
            push!(timestamps, v2["updateTime"]/1000)
        else
            @info("no valid timestamp")
            push!(timestamps, datetime2unix(now()))
        end
        buy_sign = 1
        if v2["side"] == "BUY"
            new_closed["buys_exec"] += 1
            buy_sign = -1
        end
        if !haskey(v2, "trades")
            sleep(0.3)
            trades = client.get_my_trades(symbol=v2["symbol"], orderId=v2["orderId"])
            v2["trades"] = trades
            v2 = add_fees_to_closed_order(v2, config, "", client)
        else
            trades = v2["trades"]
        end

        for t in v2["trades"]
            new_closed["symbol"] = t["symbol"]
            #
            ## should be deprecated
#             commission = typeof(t["commission"]) == Float64 ? t["commission"] : parse(Float64, t["commission"])
#             if !haskey(t, "com")
#                 if commission != 0
#                     if t["commissionAsset"] == config["farm_coin"]
#                         commission *= parse(Float64, t["price"])
#                     elseif t["commissionAsset"] != config["base_coin"]
#                         translator = string(t["commissionAsset"], config["base_coin"])
#                         #
#                         # Switch to historical data !!!
#                         #
#         #                                 commission *= parse(Float64, client.get_ticker(symbol=translator)["weightedAvgPrice"])
#                         weightedAvgPrice = try client.get_ticker(symbol=translator)["weightedAvgPrice"]
#                         catch
#                             string(1/parse(Float64, client.get_ticker(symbol=string(config["farm_coin"], t["commissionAsset"]))["weightedAvgPrice"]))
#                         end
#                         if typeof(weightedAvgPrice) == String
#                             weightedAvgPrice = parse(Float64, weightedAvgPrice)
#                         end
#                         println(weightedAvgPrice, " -- ", typeof(weightedAvgPrice))
#                         commission *= weightedAvgPrice
#                     end
#                 end
#             else

#             end
            ##
            #
            commission = typeof(t["com"]) == Float64 ? t["com"] : parse(Float64, t["com"])
            tmp_pro += buy_sign * parse(Float64, t["quoteQty"])
            tmp_pro -= commission
            new_closed["commission"] += commission
            tmp_vol += parse(Float64, t["qty"])
        end
    end

    new_closed["duration"]       = datetime2unix(unix2datetime(maximum(timestamps))) - datetime2unix(DateTime(split(id_str, "--")[end]))
    new_closed["trigger"]        = split(id_str, "--")[1]
    new_closed["profit"]         = tmp_pro
    new_closed["profit_percent"] = tmp_pro * 100 / status["volume"]
    new_closed["profit_thes"]    = tmp_pro * config["thes_factor"]
    new_closed["sell_time"]      = string(unix2datetime(maximum(timestamps)))

    #
    ## If base coin is NOT EUR
    if config["base_coin"] != "EUR"
        translator = string(config["base_coin"], "EUR")
        @info(translator)
        if weightedAvgTranslator == 0
            weightedAvgPrice = client.get_ticker(symbol=translator)["weightedAvgPrice"]
            if typeof(weightedAvgPrice) == String
                weightedAvgPrice = parse(Float64, weightedAvgPrice)
            end
        else
             weightedAvgPrice = weightedAvgTranslator
        end
        @info(typeof(weightedAvgPrice))
        @info(weightedAvgPrice)
        @info(new_closed["profit"])
        new_closed["weightedAvgPrice"] = weightedAvgPrice
        new_closed["profit_EUR"] = new_closed["profit"] * weightedAvgPrice
    else
        new_closed["profit_EUR"] = new_closed["profit"]
    end
    return new_closed, sell_date, yr, M, d, tmp_vol, tmp_pro
end
export determine_closed


text:
md(s) = string(AsMarkdown(s))

function hasmethod_invariant(fn, args...; title = _title_hasmethod(fn, args), kwargs...)

    return invariant(title; kwargs...) do inputs
        if !(_validate_hasmethod(args)(inputs))
            return "Got invalid inputs $inputs"
        end
        argnames = map(arg -> arg isa Symbol ? arg : arg[1], args)
        argvalues = map(arg -> arg isa Symbol ? inputs[arg] : arg[2], args)
        try
            fn(argvalues...)
            return nothing
        catch e
            sig = _signature(fn, argnames, argvalues)
            if e isa MethodError && e.f == fn
                return md("""When calling `$fn`, got a `MethodError`. This means that there
                is no method implemented for the given arguments. To fix this, please
                implement the following method:
                """) * "\n\n    " * sig
            else
                return (md("When calling `$fn`, got an unexpected error:") * "\n\n" *
                    (sprint(Base.showerror, e; context = (:color => false,)) |> indent |> faint) *
                    "\n\n" * md("""This means that there is a method matching the given arguments,
                    but calling it throws an error. To fix this, please debug the following
                    method:""") * "\n\n    " * sig)
            end
        end
    end
end

indent(s, n = 4) = wrap(s; initial_indent=repeat(" ", n), subsequent_indent=repeat(" ", n))
faint(s) = "\e[2m$s\e[22m"

function _title_hasmethod(fn, args)
    buf = IOBuffer()
    print(buf, "Method `$(nameof(parentmodule(fn))).$(nameof(fn))(")
    for (i, arg) in enumerate(args)
        name = arg isa Symbol ? arg : first(arg)
        print(buf, name)
        i != length(args) && print(buf, ", ")
    end
    print(buf, ")` implemented")
    return String(take!(buf))
end

function _validate_hasmethod(args)
    return function (inputs)
        for arg in args
            inputs isa NamedTuple || return false
            if arg isa Symbol
                haskey(inputs, arg) || return false
            end
        end
        return true
    end
end

function _signature(fn, argnames, argvalues)
    buf = IOBuffer()
    bold(s) = "\e[1m$s\e[22m"
    print(buf, parentmodule(fn), ".", nameof(fn), faint("("))
    for (i, (name, val)) in enumerate(zip(argnames, argvalues))
        print(buf, name, faint("::"), bold(nameof(typeof(val))))
        i != length(argnames) && print(buf, faint(", "))
    end
    print(buf, faint(")"))
    return String(take!(buf))
end


text:
@testset "Solvers" begin
  @testset "CookieCutter" begin
    problem = SimulationProblem(RegularGrid{Float64}(100,100), (:facies => Int, :property => Float64), 3)
    solver = CookieCutter(DummySimSolver(:facies => NamedTuple()),
                          [0 => DummySimSolver(), 1 => DummySimSolver()])

    @test sprint(show, solver) == "CookieCutter"
    @test sprint(show, MIME"text/plain"(), solver) == "CookieCutter\n  └─facies ⇨ DummySimSolver\n    └─0 ⇨ DummySimSolver\n    └─1 ⇨ DummySimSolver\n"

    Random.seed!(1234)
    solution = solve(problem, solver)

    if visualtests
      gr(size=(800,600))
      @plottest plot(solution) joinpath(datadir,"CookieCutter.png") !istravis
    end
  end

  @testset "SeqSim" begin
    sdomain = RegularGrid{Float64}(100, 100)
    problem = SimulationProblem(sdomain, :var => Float64, 3)
    solver = SeqSim(:var => (estimator=DummyEstimator(),
                             neighborhood=BallNeighborhood{2}(10.),
                             minneighbors=1, maxneighbors=10,
                             marginal=Normal(), path=LinearPath()))

    Random.seed!(1234)
    solution = solve(problem, solver)

    if visualtests
      gr(size=(900,300))
      @plottest plot(solution) joinpath(datadir,"SeqSim.png") !istravis
    end
  end
end


text:
module ACRLab

using LinearAlgebra
using StaticArrays
using Rotations
using Rotations: rotation_error
using ConstrainedDynamics
using ConstrainedDynamicsVis
using ConstrainedControl
using Interact
using Blink

export main,
    Controller1,
    Controller2,
    Controller3,
    Controller4

include("e1.jl")
include("e2.jl")
include("e3.jl")
include("e4.jl")
include("functions.jl")
include("gui.jl")

end


text:
using TestPackage
using Test

@testset "TestPackage.jl" begin
    @test 8 == testFunction(4)
end


text:
using StaticArrays
using MieSeries
using Plots

n_terms = 80   # number of terms in Mie series
f = 300e6      # frequency
λ = 3e8 / f    # wavelength
k = 2π / λ     # wavenumber
a = λ          # sphere radius

## Vertical x-polarized
TOL = 1e-5   # tolerance
θlist = range(0+TOL, π-TOL, length=100)
mesh = [SVector(sin(θ), 0., cos(θ)) for θ in θlist]

E_far = MieSeries.mieseries(mesh; k, a, n_terms, farfield=true)
E_far_norm2 = [sum(abs2.(E)) for E in E_far]
σ = 4*π*E_far_norm2
σ_dBsm = 10*log10.(σ/(π*a^2))

# plot
plot(θlist, σ_dBsm, legend=false)
title!("Bistatic RCS vertical plane x-polarized")
xlabel!("θ (rad)")
ylabel!("bistatic RCS (dBsm)")

## Vertical y-polarized
θlist = range(0+TOL, π-TOL, length=100)
mesh = [SVector(0., -sin(θ), cos(θ)) for θ in θlist]

E_far = MieSeries.mieseries(mesh; k, a, n_terms, farfield=true)
E_far_norm2 = [sum(abs2.(E)) for E in E_far]
σ = 4*π*E_far_norm2
σ_dBsm = 10*log10.(σ/(π*a^2))

# plot
plot(θlist, σ_dBsm, legend=false)
title!("Bistatic RCS vertical plane y-polarized")
xlabel!("θ (rad)")
ylabel!("bistatic RCS (dBsm)")

## Horizontal
ϕlist = range(0, π, length=100)
mesh = [SVector(cos(ϕ), sin(ϕ), 0.) for ϕ in ϕlist]

E_far = MieSeries.mieseries(mesh; k, a, n_terms, farfield=true)
E_far_norm2 = [sum(abs2.(E)) for E in E_far]
σ = 4*π*E_far_norm2
σ_dBsm = 10*log10.(σ/(π*a^2))

# plot
plot(ϕlist, σ_dBsm, legend=false)
title!("Bistatic RCS horizontal plane")
xlabel!("ϕ (rad)")
ylabel!("bistatic RCS (dBsm)")

max_stars_repo_name: Riarrieta/MieSeries

max_issues_repo_name: Riarrieta/MieSeries

max_forks_repo_name: Riarrieta/MieSeries


text:
using RoME
using TransformUtils
using Statistics
# , IncrementalInference, TransformUtils
# using Distributions
# using KernelDensityEstimate
using Test

##

@testset "test 3D convolutions and products" begin

##

tf = SE3([0.0;0.0;0.0], TU.AngleAxis(pi/4,[0;0;1.0]))# Euler(pi/4,0.0,0.0) )

N = 1
initCov = 0.0001*Matrix{Float64}(LinearAlgebra.I, 6,6)
[initCov[i,i] = 0.000001 for i in 4:6];
odoCov = deepcopy(initCov)
odo = Pose3Pose3( MvNormal(veeEuler(tf), odoCov) )

X = [0.01*randn(5,N); (0*pi/4 .+ 0.01*randn(1,N))]

fg = initfg()
addVariable!(fg, :x0, Pose3)
X_ = manikde!(X, ones(6), Pose3)
initManual!(fg, :x0, X_)

addVariable!(fg, :x1, Pose3)
addFactor!(fg, [:x0;:x1], odo)

Y = approxConv(fg, :x0x1f1, :x1, N=1) # X ⊕ odo

@error "TODO fix approxConv(.. N=1) case"
# @test norm(Y[1:3]-zeros(3)) < 1.0
# @test norm(Y[4:6]-[zeros(2);pi/4]) < 0.15

##

tf = SE3([0.0;0.0;0.0], TU.AngleAxis(pi/4,[0;0;1.0]))# Euler(pi/4,0.0,0.0) )

N = 1
initCov = 0.01*Matrix{Float64}(LinearAlgebra.I, 6,6)
[initCov[i,i] = 0.001 for i in 4:6];
odoCov = deepcopy(initCov)
odo = Pose3Pose3(  MvNormal(veeEuler(tf), odoCov) )

X = [0.01*randn(5,N); (0*pi/4.0 .+ 0.01*randn(1,N))]

fg = initfg()
addVariable!(fg, :x0, Pose3)
X_ = manikde!(X, ones(6), Pose3)
initManual!(fg, :x0, X_)

addVariable!(fg, :x1, Pose3)
addFactor!(fg, [:x0;:x1], odo)

Y = approxConv(fg, :x0x1f1, :x1, N=1) # X ⊕ odo

@error "TODO fix approxConv(.. N=1) case"
# @test norm(Y[1:3]-zeros(3)) < 1.0
# @test norm(Y[4:6]-[zeros(2);pi/4]) < 0.2

##

end

##

global N = 100
global fg = initfg()

global initCov = Matrix{Float64}(LinearAlgebra.I, 6,6)
[initCov[i,i] = 0.01 for i in 4:6];
global odoCov = deepcopy(initCov)

@testset "Adding PriorPose3 to graph..." begin

global v1 = addVariable!(fg, :x1, Pose3,  N=N) # 0.1*randn(6,N)
global initPosePrior = PriorPose3( MvNormal(zeros(6), initCov) )
global f1  = addFactor!(fg,[:x1;], initPosePrior)
# @test !isInitialized(fg, :x1)

end

##

@testset "Ensure vertex initialized properly" begin
  # start with initialization
  ensureAllInitialized!(fg)
  @test isInitialized(fg, :x1)
  @show muX1 = Statistics.mean(getVal(fg,:x1),dims=2)
  @show stdX1 = Statistics.std(getVal(fg,:x1),dims=2)
  @test sum(map(Int,abs.(muX1[1:3]) .< 0.5)) == 3
  @test sum(map(Int,abs.(muX1[4:6]) .< 0.05)) == 3
  @test sum(map(Int, 0.5 .< stdX1[1:3] .< 1.5)) == 3
  @test sum(map(Int, 0.05 .< stdX1[4:6] .< 0.15)) == 3
end

@testset "Testing PriorPose3 evaluation..." begin

global priorpts = approxConv(fg, :x1f1, :x1) # fg.g.vertices[2], 1
global means = Statistics.mean(priorpts,dims=2)
@test sum(map(Int,abs.(means[1:3]) .> 0.5)) == 0
@test sum(map(Int,abs.(means[4:6]) .> 0.05)) == 0

end

##

@testset "Adding Pose3Pose3 to graph..." begin

##

global N
global odo = SE3([10;0;0], Quaternion(0))
# global pts0X2 = projectParticles(getVal(fg,:x1), MvNormal(veeEuler(odo), odoCov) )
global odoconstr = Pose3Pose3( MvNormal(veeEuler(odo), odoCov) )
global v2 = addVariable!(fg,:x2, Pose3, N=N) # pts0X2
addFactor!(fg,[:x1;:x2],odoconstr, inflation=0.1)
# @test !isInitialized(fg, :x2)

## test opposites

# should be all zero
res = calcFactorResidual(fg, :x1x2f1, [10;0;0;0;0;0.0], zeros(6), [10;0;0;0;0;0])
@test norm(res) < 1e-10

# trivial fail case
# res = calcFactorResidual(fg, :x1x2f1, [10;0;0;0;0;0.0], zeros(6), [10;0;0;pi;pi;pi])
@warn "suppressing trivial Pose3 fail case until RoME.jl #244 has been completed."

## test following introduction of inflation, see IIF #1051

# force the inflation trivial error, https://github.com/JuliaRobotics/RoME.jl/issues/380#issuecomment-778795848
# IIF._getCCW(fg,:x1x2f1).inflation = 10.0

pts = approxConv(fg, :x1x2f1, :x2)
# X2 = manikde!(pts, Pose3)

# test translations through convolution
@test 0.8N < sum( 5 .< pts[1,:] .<15 )
@test 0.7N < sum( -5 .< pts[2,:] .< 5 )
@test 0.7N < sum( -5 .< pts[3,:] .< 5 )

# test rotations through convolution
@test 0.95N < sum( -2 .< pts[4,:] .< 2 )
@test 0.95N < sum( -2 .< pts[5,:] .< 2 )
@test 0.95N < sum( -2 .< pts[6,:] .< 2 )

##

end

##

# Noticed a DomainError on convolutions here after mutlithreading upgrade.  Previously used fill(PP3REUSE, Threads.nthreads())
@testset "Testing Pose3Pose3 evaluation..." begin

ensureAllInitialized!(fg)
@test isInitialized(fg, :x2)
global X1pts = approxConv(fg, :x1x2f1, :x1)
# X1pts = evalFactor(fg, fg.g.vertices[4], 1)
global X2pts = approxConv(fg, :x1x2f1, :x2)
# X2pts = evalFactor(fg, fg.g.vertices[4], 3)
global X2ptsMean = Statistics.mean(X2pts,dims=2)
global X1ptsMean = Statistics.mean(X1pts,dims=2)
@show X1ptsMean
@test  sum(map(Int, abs.(X1ptsMean) .< 1.25 )) == 6
@test  sum(map(Int, abs.(X2ptsMean .- [10.0;0;0;0;0;0]) .< 1.25 )) == 6

end

@testset "Construct Bayes tree and perform inference..." begin
  tree, smt, hist = solveTree!(fg)
  @test true
end

##

@testset "Ensure basic parameters on x1,x2 after inference..." begin

##

# check mean and covariances after one up and down pass over the tree
global muX1 = getPPE(fg, :x1).suggested # Statistics.mean(getVal(fg,:x1),dims=2)
@test sum(map(Int,abs.(muX1[1:3]) .< 1.0)) == 3

# sidestep trivial case, #412
# case where points land on trivial rotation error [pi;pi;pi]
@warn "simplify mean and std testing on Pose3 after #244, see #412"
pts = getVal(fg,:x1)
mask = (2.7 .< abs.(pts[4,:])) .& (2.7 .< abs.(pts[5,:])) .& (2.7 .< abs.(pts[6,:]))
mask .= xor.(mask,1)
mu1tmp = Statistics.mean(pts[:,mask],dims=2)
@test sum(map(Int, abs.(mu1tmp[4:6]) .< 0.1)) == 3
global stdX1 = Statistics.std(pts[:,mask],dims=2)

@test sum(map(Int, 0.4 .< stdX1[1:3] .< 1.6)) == 3 # had a 2==3 failure here
@show stdX1[4:6]
@test sum(map(Int, 0.02 .< stdX1[4:6] .< 0.5)) == 3
global muX2 = getPPE(fg, :x2).suggested # Statistics.mean(getVal(fg,:x2),dims=2)
@show muX2[1:3]-[10.0;0;0]
@test sum(map(Int, abs.(muX2[1:3]-[10.0;0;0]) .< 1.5)) == 3

# case where points land on trivial rotation error [pi;pi;pi]
@warn "simplify mean and std testing on Pose3 after #244, see #412"
pts = getVal(fg,:x2)
mask = (2.7 .< abs.(pts[4,:])) .& (2.7 .< abs.(pts[5,:])) .& (2.7 .< abs.(pts[6,:]))
mask .= xor.(mask,1)
mu2tmp = Statistics.mean(pts[:,mask],dims=2)
@test sum(map(Int, abs.(mu2tmp[4:6]) .< 0.1)) == 3
global stdX2 = Statistics.std(pts[:,mask],dims=2)
# @test sum(map(Int, abs.(muX2[4:6]) .< 0.1)) == 3
# global stdX2 = Statistics.std(getVal(fg,:x2),dims=2)

println("previous test failure 0.75 .< $(round.(stdX2[1:3],digits=2)) .< 2.25")
@test sum(map(Int, 0.75 .< stdX2[1:3] .< 2.5)) == 3
println("previous test failure 0.05 .< $(round.(stdX2[4:6],digits=2)) .< 0.35")
@test sum(map(Int, 0.05 .< stdX2[4:6] .< 0.5)) == 3

##

end

# println("Plot marginals to see what is happening")
# plotKDE(marginal(getBelief(fg,:x1),[1]))
# plotKDE(marginal(getBelief(fg,:x2),[1]))

#


text:
function metric(J::AbstractArray{T,2}) where {T}
    g = J'*J
    return g
end

function covariant_basis(c::Coordinate{T}) where {T}
    return ForwardDiff.jacobian(c.R,c.u)
end

function covariant_metric(c::Coordinate)
    J = covariant_basis(c)
    return metric(J)
end

function contravariant_basis(c::Coordinate)
    Jco = covariant_basis(c)
    gco = metric(Jco)
    g = inv(gco)
    J = Jco*g
    return J
end

function contravariant_metric(c::Coordinate)
    J = contravariant_basis(c)
    g = metric(J)
end

abstract type CoordinateVector{T} <: AbstractArray{T,1} end

struct ContravariantVector{T,F} <: CoordinateVector{T}
    data::Vector{T}       #Contravariant Components of the Vector
    coord::Coordinate{T,F}  #Vector Coordinate
    J::Matrix{T}
    g::Matrix{T}
    h::Vector{T}
end

function ContravariantVector(data::AbstractVector{T}, coord::Coordinate{T,F}; unit_basis = false) where {T,F}
    J = covariant_basis(coord)
    g = metric(J)
    h = sqrt.(diag(g))
    if unit_basis
        data = data./h
    end
    ContravariantVector{T,F}(data, coord, J, g, h)
end

struct CovariantVector{T,F} <: CoordinateVector{T}
    data::Vector{T}       #Covariant Components of the Vector
    coord::Coordinate{T,F}  #Vector Coordinate
    J::Matrix{T}
    g::Matrix{T}
    h::Vector{T}
end

function CovariantVector(data::AbstractVector{T}, coord::Coordinate{T,F}; unit_basis = false) where {T,F}
    J = contravariant_basis(coord)
    g = metric(J)
    h = sqrt.(diag(g))
    if unit_basis
        data = data./h
    end
    CovariantVector{T,F}(data, coord, J, g, h)
end

# Array Interface for CoordinateVector
parent(A::T) where {T<:CoordinateVector} = A.data
size(A::T) where {T<:CoordinateVector} = size(A.data)
axes(A::T) where {T<:CoordinateVector} = axes(A.data)
parenttype(::Type{CoordinateVector{T}}) where {T} = Vector{T}
IndexStyle(::Type{T}) where {T<:CoordinateVector} = IndexStyle(parenttype(T))

@propagate_inbounds getindex(A::T, i::Int) where {T<:CoordinateVector} = A.data[i]
@propagate_inbounds setindex!(A::CoordinateVector{T}, v::T, i::Int) where {T} = A.data[i] = v

function convert(::Type{T}, x::S) where {T<:CoordinateVector,S<:CoordinateVector}
    data = x.g*x.data
    g = inv(x.g)
    J = x.J*g
    h = sqrt.(diag(g))
    T(data, x.coord, J, g, h)
end

convert(::Type{T}, x::T) where {T<:CoordinateVector} = x
convert(x::CovariantVector) = convert(ContravariantVector,x)
convert(x::ContravariantVector) = convert(CovariantVector,x)

function unit_basis_components(x::T) where {T<:CoordinateVector}
    return x .* x.h
end


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"12.0.1-iains")


text:
using Test, MPI
include(joinpath("..", "..", "testhelpers.jl"))

include("iterativesolvers.jl")

@testset "Linear Solvers Poisson" begin
    runmpi(joinpath(@__DIR__, "columnwiselu.jl"))
    runmpi(joinpath(@__DIR__, "poisson.jl"))
    runmpi(joinpath(@__DIR__, "bandedsystem.jl"))
    runmpi(joinpath(@__DIR__, "cg.jl"))
    runmpi(joinpath(@__DIR__, "bgmres.jl"))
end