text:
using HDF5, JLD
using Random
using Profile
using LinearAlgebra
using Test
# Define variables of different types
x = 3.7
A = reshape(collect(1:15), 3, 5)
Aarray = Vector{Float64}[[1.2,1.3],[2.2,2.3,2.4]]
str = "Hello"
stringsA = String["It", "was", "a", "dark", "and", "stormy", "night"]
stringsU = String["It", "was", "a", "dark", "and", "stormy", "night"]
empty_string = ""
empty_string_array = String[]
empty_array_of_strings = String[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = ComplexF32(3,7)
cint = 1+im  # issue 108
C = reinterpret(ComplexF64, vec(B))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
simple_bool = true
try
    global MyStruct
    mutable struct MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
symhard = Symbol("troublesome \"symbol\"")
d = Dict([(syms[1],"aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
# anonymous function
module FunConstMod
    a = 2
    module Sub
        a = 4
    end
end
funconst = 11
fun = (x, y) -> funconst * x + y
function_referencing_module = (x, y) -> FunConstMod.a * x + FunConstMod.Sub.a * y
T = UInt8
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]]  # issue #123
typevar = Array{Int}[[1]]
eval(Meta.parse("typevar_lb = (Vector{U} where U<:Integer)[[1]]"))
eval(Meta.parse("typevar_ub = (Vector{U} where Int<:U<:Any)[[1]]"))
eval(Meta.parse("typevar_lb_ub = (Vector{U} where Int<:U<:Real)[[1]]"))
# Unexported type:
cpus = Base.Sys.cpu_info()
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
struct ObjWithPointer
    a::Ptr{Nothing}
end
objwithpointer = ObjWithPointer(0)
# Custom PrimitiveType (#99)
primitive type MyBT 64 end
bt = reinterpret(MyBT, Int64(55))
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = view([1:5;], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = (Tuple{})[]
struct EmptyImmutable end
emptyimmutable = EmptyImmutable()
arr_emptyimmutable = [emptyimmutable]
mutable struct EmptyType end
emptytype = EmptyType()
arr_emptytype = [emptytype]
struct EmptyII
    x::EmptyImmutable
end
emptyii = EmptyII(EmptyImmutable())
struct EmptyIT
    x::EmptyType
end
emptyit = EmptyIT(EmptyType())
mutable struct EmptyTI
    x::EmptyImmutable
end
emptyti = EmptyTI(EmptyImmutable())
mutable struct EmptyTT
    x::EmptyType
end
emptytt = EmptyTT(EmptyType())
struct EmptyIIOtherField
    x::EmptyImmutable
    y::Float64
end
emptyiiotherfield = EmptyIIOtherField(EmptyImmutable(), 5.0)

# Unicode type field names (#118)
mutable struct MyUnicodeStruct☺{τ}
    α::τ
    ∂ₓα::τ
    MyUnicodeStruct☺{τ}(α::τ, ∂ₓα::τ) where {τ} = new{τ}(α, ∂ₓα)
end
unicodestruct☺ = MyUnicodeStruct☺{Float64}(1.0, -1.0)
# Arrays of matrices (#131)
array_of_matrices = Matrix{Int}[[1 2; 3 4], [5 6; 7 8]]
# Tuple of arrays and bitstype
tup = (1, 2, [1, 2], [1 2; 3 4], bt)
# Empty tuple
empty_tup = ()
# Non-pointer-free immutable
struct MyImmutable{T}
    x::Int
    y::Vector{T}
    z::Bool
end
nonpointerfree_immutable_1 = MyImmutable(1, [1., 2., 3.], false)
nonpointerfree_immutable_2 = MyImmutable(2, Any[3., 4., 5.], true)
struct MyImmutable2
    x::Vector{Int}
    MyImmutable2() = new()
end
nonpointerfree_immutable_3 = MyImmutable2()
# Immutable with a non-concrete datatype (issue #143)
struct Vague
    x::Integer
end
vague = Vague(7)
# Immutable with a union of BitsTypes
struct BitsUnion
    x::Union{Int64, Float64}
end
bitsunion = BitsUnion(5.0)
# Immutable with a union of Types
let UT = eval(Meta.parse("Type{T} where T <: Union{Int64, Float64}"))
    @eval struct TypeUnionField
        x::$UT
    end
end
typeunionfield = TypeUnionField(Int64)
# Generic union type field
struct GenericUnionField
    x::Union{Vector{Int},Int}
end
genericunionfield = GenericUnionField(1)
# Array references
arr_contained = [1, 2, 3]
arr_ref = typeof(arr_contained)[]
push!(arr_ref, arr_contained, arr_contained)
# Object references
mutable struct ObjRefType
    x::ObjRefType
    y::ObjRefType
    ObjRefType() = new()
    ObjRefType(x, y) = new(x, y)
end
ref1 = ObjRefType()
obj_ref = ObjRefType(ObjRefType(ref1, ref1), ObjRefType(ref1, ref1))
# Immutable that requires padding between elements in array
struct PaddingTest
    x::Int64
    y::Int8
end
padding_test = PaddingTest[PaddingTest(i, i) for i = 1:8]
# Empty arrays of various types and sizes
empty_arr_1 = Int[]
empty_arr_2 = Matrix{Int}(undef, 56, 0)
empty_arr_3 = Any[]
empty_arr_4 = Matrix{Any}(undef, 0, 97)
# Moderately big dataset (which will be mmapped)
bigdata = [1:10000;]
# BigFloats and BigInts
bigints = big(3).^(1:100)
bigfloats = big(3.2).^(1:100)
# None
none = Union{}
nonearr = Vector{Union{}}(undef, 5)
# nothing
scalar_nothing = nothing
vector_nothing = Union{Int,Nothing}[1,nothing]

# some data big enough to ensure that compression is used:
Abig = kron(Matrix(1.0I, 10, 10), rand(20,20))
Bbig = Any[i for i=1:3000]
Sbig = "A test string "^1000

# Bits type type parameters
mutable struct BitsParams{x}; end
bitsparamfloat  = BitsParams{1.0}()
bitsparambool   = BitsParams{true}()
bitsparamsymbol = BitsParams{:x}()
bitsparamint    = BitsParams{1}()
bitsparamuint   = BitsParams{0x01}()
bitsparamint16  = BitsParams{Int16(1)}()

# Tuple of tuples
tuple_of_tuples = (1, 2, (3, 4, [5, 6]), [7, 8])

# SimpleVector
simplevec = Core.svec(1, 2, Int64, "foo")
iseq(x::Core.SimpleVector, y::Core.SimpleVector) = collect(x) == collect(y)

namedtupl = (x = 1, y = 2, z = "hi")

# Issue #243
# Type that overloads != so that it is not boolean
struct NALikeType; end
Base.:(!=)(::NALikeType, ::NALikeType) = NALikeType()
Base.:(!=)(::NALikeType, ::Nothing) = NALikeType()
Base.:(!=)(::Nothing, ::NALikeType) = NALikeType()
natyperef = Any[NALikeType(), NALikeType()]

# Issue #110
ver = v"0.1.2"

iseq(x,y) = isequal(x,y)
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(x::MyImmutable, y::MyImmutable) = (isequal(x.x, y.x) && isequal(x.y, y.y) && isequal(x.z, y.z))
iseq(x::EmptyType, y::EmptyType) = true
iseq(x::EmptyIT, y::EmptyIT) = true
iseq(x::Array{EmptyType}, y::Array{EmptyType}) = size(x) == size(y)
iseq(x::BitsParams{T}, y::BitsParams{T}) where {T} = true
iseq(x::BitsParams, y::BitsParams) = false
iseq(x::Union{EmptyTI, EmptyTT}, y::Union{EmptyTI, EmptyTT}) = iseq(x.x, y.x)
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in fieldnames(Base.Sys.CPUinfo)
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
iseq(x::MyUnicodeStruct☺, y::MyUnicodeStruct☺) = (x.α == y.α && x.∂ₓα == y.∂ₓα)
iseq(x::Array{Union{}}, y::Array{Union{}}) = size(x) == size(y)
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                @warn string("Error reading ", $(string(sym)))
                rethrow(e)
            end
            if !iseq(tmp, $sym)
                written = $sym
                error("For ", $(string(sym)), ", read value $tmp does not agree with written value $written")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    # try to make test failure backtraces easier to follow
    Base.remove_linenums!(ex)
    pushfirst!(ex.args, __source__)
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @test a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if (isa(a.args[i], Expr) && a.args[i].head == :line) || isa(a.args[i], LineNumberNode)
            i += 1
            continue
        end
        if (isa(b.args[j], Expr) && b.args[j].head == :line) || isa(b.args[j], LineNumberNode)
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

# check for equality of function asts
# skip line numbers in function body,
# as well as lines asserting a variable is type Any
function isAssertAny(line::Expr)
    headcheck = (line.head == :(=))
    if !headcheck
        return false
    end
    assertcheck = (length(line.args) >= 2 && isa(line.args[2], Expr) && length(line.args[2].args) >= 1 &&
        isa(line.args[2].args[1], TopNode) && line.args[2].args[1].name == :typeassert)
    if !assertcheck
        return false
    end
    anycheck = (length(line.args[2].args) >= 3 && (eval(line.args[2].args[3]) == Base.Any))
    return anycheck
end
function checkfuns(f, g)
    f_ast = Base.uncompressed_ast(f.code)
    g_ast = Base.uncompressed_ast(g.code)
    @assert f_ast.args[1] == g_ast.args[1]

    f_body = f_ast.args[3]
    g_body = g_ast.args[3]
    checkfunexpr(f_body, g_body)

    @assert f.code.module == g.code.module
end
checkfunexpr(a, b) = @test a == b
function checkfunexpr(f_body::Expr, g_body::Expr)
    i = 1
    j = 1
    while i <= length(f_body.args) && j <= length(g_body.args)
        if (isa(f_body.args[i], Expr) && ((f_body.args[i].head == :line) || isAssertAny(f_body.args[i]))) || isa(f_body.args[i], LineNumberNode)
            i += 1
            continue
        end
        if (isa(g_body.args[j], Expr) && ((g_body.args[j].head == :line) || isAssertAny(g_body.args[j]))) || isa(g_body.args[i], LineNumberNode)
            j += 1
            continue
        end
        checkfunexpr(f_body.args[i], g_body.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(f_body.args) && j >= length(g_body.args)
end

fn = joinpath(tempdir(),"test-$(randstring()).jld")

# Issue #106
module Mod106
primitive type Typ{T} 64 end
typ(x::Int64, ::Type{T}) where {T} = reinterpret(Typ{T}, x)
abstract type UnexportedT end
end

# test mmapping of small arrays (Issue #192)
let fid = jldopen(fn, "w", mmaparrays = true)
    write(fid, "a", [1:3;])
    @test ismmappable(fid["a"])
    close(fid)
    rm(fn)
end

let fid = jldopen(fn, "w", mmaparrays=false)
    write(fid, "a", [1:3;]; mmap = true)
    @test ismmappable(fid["a"])
    close(fid)
    rm(fn)
end

let fid = jldopen(fn, "w", compress = true)
    write(fid, "a", [1:3;])
    @test ismmappable(fid["a"]) == false
    close(fid)
    rm(fn)
end

let fid = jldopen(fn, "w", compatible = true, compress = true)
    write(fid, "a", [1:3;])
    @test ismmappable(fid["a"]) == false
    close(fid)
    rm(fn)
end

# Hyperslab
for compatible in (false, true), compress in (false, true)
    jldopen(fn, "w", compatible=compatible, compress=compress) do fid
        write(fid, "a", [1:3;])
        aset = fid["a"]
        aset[1:2] = [5,7]
        b = read(fid, "a")
        @test b == [5,7,3]
        @test aset[2:3] == [7,3]
        aset[3] = 27
        @test aset[1:3] == [5,7,27]
        Arnd = rand(5,3)
        write(fid, "A", Arnd)
        Aset = fid["A"]
        Aset[:,2] = 15   # FIXME: broadcasting with .= doesn't work for JldDataset
        Arnd[:,2] .= 15
        @test read(fid, "A") == Arnd
    end
end

for compatible in (false, true), compress in (false, true)
    undefv = Vector{Any}(undef, 1)
    undefm = Matrix{Any}(undef, 2, 2)
    ms_undef = MyStruct(0)
    fid = jldopen(fn, "w", compatible=compatible, compress=compress)
    @write fid x
    @write fid A
    @write fid Aarray
    @write fid str
    @write fid stringsA
    @write fid stringsU
    @write fid empty_string
    @write fid empty_string_array
    @write fid empty_array_of_strings
    @write fid tf
    @write fid TF
    @write fid AB
    @write fid t
    @write fid c
    @write fid cint
    @write fid C
    @write fid emptyA
    @write fid emptyB
    @write fid ms
    @write fid msempty
    @write fid sym
    @write fid syms
    @write fid symhard
    @write fid d
    @write fid ex
    @write fid T
    @write fid char
    @write fid unicode_char
    @write fid α
    @write fid β
    @write fid vv
    @write fid cpus
    @write fid rng
    @write fid typevar
    @write fid typevar_lb
    @write fid typevar_ub
    @write fid typevar_lb_ub
    @write fid undefv
    @write fid undefm
    @write fid ms_undef
    @write fid objwithpointer
    @write fid bt
    @write fid sa_asc
    @write fid sa_utf8
    @write fid subarray
    @write fid arr_empty_tuple
    @write fid emptyimmutable
    @write fid emptytype
    @write fid arr_emptyimmutable
    @write fid arr_emptytype
    @write fid emptyii
    @write fid emptyit
    @write fid emptyti
    @write fid emptytt
    @write fid emptyiiotherfield
    @write fid unicodestruct☺
    @write fid array_of_matrices
    @write fid tup
    @write fid empty_tup
    @write fid nonpointerfree_immutable_1
    @write fid nonpointerfree_immutable_2
    @write fid nonpointerfree_immutable_3
    @write fid vague
    @write fid bitsunion
    @write fid typeunionfield
    @write fid genericunionfield
    @write fid arr_ref
    @write fid obj_ref
    @write fid padding_test
    @write fid empty_arr_1
    @write fid empty_arr_2
    @write fid empty_arr_3
    @write fid empty_arr_4
    @write fid bigdata
    @write fid bigfloats
    @write fid bigints
    @write fid none
    @write fid nonearr
    @write fid scalar_nothing
    @write fid vector_nothing
    @write fid Abig
    @write fid Bbig
    @write fid Sbig
    @test_throws ErrorException @write fid bitsparamint16
    @write fid bitsparamfloat
    @write fid bitsparambool
    @write fid bitsparamsymbol
    @write fid bitsparamint
    @write fid bitsparamuint
    @write fid tuple_of_tuples
    @write fid simplevec
    @write fid natyperef
    @write fid ver
    @write fid simple_bool
    @write fid namedtupl

    # Make sure we can create groups (i.e., use HDF5 features)
    g = create_group(fid, "mygroup")
    i = 7
    @write g i
    write(fid, "group1/x", Any[1])
    write(fid, "group2/x", Any[2])
    close(fid)

    # mmapping currently fails on Windows; re-enable once it can work
    for mmap = (@static Sys.iswindows() ? false : (false, true))
        fidr = jldopen(fn, "r", mmaparrays=mmap)
        @test creator(fidr, "VERSION") == VERSION
        @test creator(fidr, "WORD_SIZE") == Sys.WORD_SIZE
        @test creator(fidr, "ENDIAN_BOM") == ENDIAN_BOM
        @check fidr x
        @check fidr A
        dsetA = fidr["A"]
        @test ndims(dsetA) == ndims(A)
        @test size(dsetA) == size(A)
        @test size(dsetA, 1) == size(A, 1)
        @test size(dsetA, 2) == size(A, 2)
        @test size(dsetA, 3) == size(A, 3)
        @test eltype(dsetA) == Int
        @check fidr Aarray
        @check fidr str
        @check fidr stringsA
        @check fidr stringsU
        @check fidr empty_string
        @check fidr empty_string_array
        @check fidr empty_array_of_strings
        @check fidr tf
        @check fidr TF
        @check fidr AB
        @check fidr t
        @check fidr c
        @check fidr cint
        @check fidr C
        @check fidr emptyA
        @check fidr emptyB
        @check fidr ms
        @check fidr msempty
        @check fidr sym
        @check fidr syms
        @check fidr symhard
        @check fidr d
        exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
        checkexpr(ex, exr)
        @check fidr T
        @check fidr char
        @check fidr unicode_char
        @check fidr α
        @check fidr β
        @check fidr vv
        @check fidr cpus
        @check fidr rng
        @check fidr typevar
        @check fidr typevar_lb
        @check fidr typevar_ub
        @check fidr typevar_lb_ub
        @check fidr simple_bool

        # Special cases for reading undefs
        undefv = read(fidr, "undefv")
        if !isa(undefv, Array{Any, 1}) || length(undefv) != 1 || isassigned(undefv, 1)
            error("For undefv, read value does not agree with written value")
        end
        undefm = read(fidr, "undefm")
        if !isa(undefm, Array{Any, 2}) || length(undefm) != 4 || any(map(i->isassigned(undefm, i), 1:4))
            error("For undefm, read value does not agree with written value")
        end
        ms_undef = read(fidr, "ms_undef")
        if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
            error("For ms_undef, read value does not agree with written value")
        end

        @check fidr bt
        @check fidr sa_asc
        @check fidr sa_utf8
        @check fidr subarray
        @check fidr arr_empty_tuple
        @check fidr emptyimmutable
        @check fidr emptytype
        @check fidr arr_emptyimmutable
        @check fidr arr_emptytype
        @check fidr emptyii
        @check fidr emptyit
        @check fidr emptyti
        @check fidr emptytt
        @check fidr emptyiiotherfield
        @check fidr unicodestruct☺
        @check fidr array_of_matrices
        @check fidr tup
        @check fidr empty_tup
        @check fidr nonpointerfree_immutable_1
        @check fidr nonpointerfree_immutable_2
        @check fidr nonpointerfree_immutable_3
        vaguer = read(fidr, "vague")
        @test typeof(vaguer) == typeof(vague) && vaguer.x == vague.x
        # @check fidr bitsunion # FIXME: fails on 0.7 with message: "reference encountered in pointerfree immutable; this is a bug"
        @check fidr typeunionfield
        @check fidr genericunionfield

        arr = read(fidr, "arr_ref")
        @test arr == arr_ref
        @test arr[1] === arr[2]

        obj = read(fidr, "obj_ref")
        @test obj.x.x === obj.x.y == obj.y.x === obj.y.y
        @test obj.x !== obj.y

        @check fidr padding_test
        @check fidr empty_arr_1
        @check fidr empty_arr_2
        @check fidr empty_arr_3
        @check fidr empty_arr_4
        !mmap && @check fidr bigdata # FIXME: fails on 0.7 due to data alignment issues
        @check fidr bigfloats
        @check fidr bigints
        @check fidr none
        @check fidr nonearr
        @check fidr scalar_nothing
        @check fidr vector_nothing
        !mmap && @check fidr Abig # FIXME: fails on 0.7 due to data alignment issues
        @check fidr Bbig
        @check fidr Sbig
        @check fidr bitsparamfloat
        @check fidr bitsparambool
        @check fidr bitsparamsymbol
        @check fidr bitsparamint
        @check fidr bitsparamuint
        @check fidr tuple_of_tuples
        @check fidr simplevec
        @check fidr natyperef
        @check fidr ver
        @check fidr namedtupl

        x1 = read(fidr, "group1/x")
        @test x1 == Any[1]
        x2 = read(fidr, "group2/x")
        @test x2 == Any[2]

        close(fidr)
    end
    # HDF5 1.10 gives an error when locking the file otherwise.
    rm(fn)
end # compress in (true,false)

for compatible in (false, true), compress in (false, true)
    # object references in a write session
    r = ObjRefType()
    a = [r, r]
    b = [r, r]
    @save fn a b
    jldopen(fn, "r") do fid
        a = read(fid, "a")
        b = read(fid, "b")
        @test a[1] === a[2] === b[2] === a[1]

        # Let gc get rid of a and b
        a = nothing
        b = nothing
        GC.gc()

        a = read(fid, "a")
        b = read(fid, "b")
        @test typeof(a[1]) == ObjRefType
        @test a[1] === a[2] === b[2] === a[1]
    end

    # do syntax
    jldopen(fn, "w", compatible=compatible, compress=compress) do fid
        create_group(fid, "mygroup") do g
            write(g, "x", 3.2)
        end
    end
    fid = jldopen(fn, "r")
    @test keys(fid) == String["mygroup"]
    g = fid["mygroup"]
    @test keys(g) == String["x"]
    @test read(g, "x") == 3.2
    close(g)
    close(fid)

    # Function load() and save() syntax
    d1 = Dict([("x",3.2), ("β",β), ("A",A)])
    save(fn, d1, compatible=compatible, compress=compress)
    d2 = load(fn)
    @test d1 == d2
    β2 = load(fn, "β")
    @test β == β2
    β2, A2 = load(fn, "β", "A")
    @test β == β2
    @test A == A2

    save(fn, "x", 3.2, "β", β, "A", A, compatible=compatible, compress=compress)
    d3 = load(fn)
    @test d1 == d3

    # #71
    jldopen(fn, "w", compatible=compatible, compress=compress) do file
        file["a"] = 1
    end
    jldopen(fn, "r") do file
        @test read(file, "a") == 1
    end

    # Issue #106
    save(fn, "i106", Mod106.typ(Int64(1), Mod106.UnexportedT), compress=compress)
    i106 = load(fn, "i106")
    @test i106 == Mod106.typ(Int64(1), Mod106.UnexportedT)

    # bracket syntax for datasets
    jldopen(fn, "w", compatible=compatible, compress=compress) do file
        file["a"] = [1:100;]
        file["b"] = [x*y for x=1:10,y=1:10]
        file["c"] = Any[1, 2, 3]
        file["d"] = [1//2, 1//4, 1//8]
    end
    jldopen(fn, "r+", compatible=compatible, compress=compress) do file
        @test(file["a"][1:50] == [1:50;])
        file["a"][1:50] = 1:2:100
        @test(file["a"][1:50] == [1:2:100;])
        @test(file["b"][5,6][1]==5*6)
        @test(file["c"][1:2] == [1, 2])
        file["c"][2:3] = [5, 7]
        @test(read(file, "c") == [1, 5, 7])
        @test(file["d"][2:3] == [1//4, 1//8])
        file["d"][1:1] = [9]
        @test(read(file, "d") == [9, 1//4, 1//8])
    end

    # bracket syntax when created by HDF5
    println("The following unrecognized JLD file warning is a sign of normal operation.")
    if compress
        h5open(fn, "w") do file
            file["a", blosc = 5] = [1:100;]
            file["a"][51:100] = [1:50;]
            file["b", blosc = 5] = [x*y for x=1:10,y=1:10]
        end
    else
        h5open(fn, "w") do file
            file["a"] = [1:100;]
            file["a"][51:100] = [1:50;]
            file["b"] = [x*y for x=1:10,y=1:10]
        end
    end
    jldopen(fn, "r") do file
        @test(file["a"][1:50] == [1:50;])
        @test(file["a"][:] == [1:50;1:50])
        @test(file["b"][5,6][1]==5*6)
    end

    # delete!
    jldopen(fn, "w", compatible=compatible, compress=compress) do file
        file["ms"] = ms
        delete!(file, "ms")
        file["ms"] = β
        g = create_group(file,"g")
        file["g/ms"] = ms
        @test_throws ErrorException delete!(file, "_refs/g/ms")
        delete!(file, "g/ms")
        file["g/ms"] = ms
        delete!(file, "/g/ms")
        g["ms"] = ms
        delete!(g,"ms")
        g["ms"] = ms
        delete!(g["ms"])
        g["ms"] = ms
        delete!(g)
        g = create_group(file,"g")
        g["ms"] = ms
        delete!(g)
    end
    jldopen(fn, "r") do file
        @test(read(file["ms"]) == β)
        @test(!haskey(file, "g/ms"))
        @test(!haskey(file, "g"))
    end

end # compress in (false,true)

# mismatched and missing types
module JLDTemp1
using JLD
import ..fn

mutable struct TestType1
    x::Int
end
mutable struct TestType2
    x::Int
end
struct TestType3
    x::TestType2
end

mutable struct TestType4
    x::Int
end
mutable struct TestType5
    x::TestType4
end
mutable struct TestType6 end
primitive type TestType7 8 end
struct TestType8
    a::TestType4
    b::TestType5
    c::TestType6
    d::TestType7
end

jldopen(fn, "w") do file
    truncate_module_path(file, JLDTemp1)
    write(file, "x1", TestType1(1))
    write(file, "x2", TestType3(TestType2(1)))
    write(file, "x3", TestType4(1))
    write(file, "x4", TestType5(TestType4(2)))
    write(file, "x5", [TestType5(TestType4(i)) for i = 1:5])
    write(file, "x6", TestType6())
    write(file, "x7", reinterpret(TestType7, 0x77))
    write(file, "x8", TestType8(TestType4(2), TestType5(TestType4(3)),
                                TestType6(), reinterpret(TestType7, 0x12)))
    write(file, "x9", (TestType4(1),
                       (TestType5(TestType4(2)),
                        [TestType5(TestType4(i)) for i = 1:5]),
                       TestType6()))
end
end

mutable struct TestType1
    x::Float64
end
mutable struct TestType2
    x::Int
end
struct TestType3
    x::TestType1
end

jldopen(fn, "r") do file
    @test_throws JLD.TypeMismatchException read(file, "x1")
    @test_throws MethodError read(file, "x2")
    println("The following missing type warnings are a sign of normal operation.")
    @test read(file, "x3").x == 1
    @test read(file, "x4").x.x == 2

    x5 = read(file, "x5")
    for i = 1:5
        @test x5[i].x.x == i
    end
    @test isempty(fieldnames(typeof(read(file, "x6"))))
    @test reinterpret(UInt8, read(file, "x7")) == 0x77

    x8 = read(file, "x8")
    @test x8.a.x == 2
    @test x8.b.x.x == 3
    @test isempty(fieldnames(typeof(x8.c)))
    @test reinterpret(UInt8, x8.d) == 0x12

    x9 = read(file, "x9")
    @test isa(x9, Tuple)
    @test length(x9) == 3
    @test x9[1].x == 1
    @test isa(x9[2], Tuple)
    @test length(x9[2]) == 2
    @test x9[2][1].x.x == 2
    for i = 1:5
        @test x9[2][2][i].x.x == i
    end
    @test isempty(fieldnames(typeof(x9[3])))
end

# Issue #176
exx = quote
    function incrementby1(x::Int)
        x+1
    end
end
for i = 1:2
    fid = jldopen(fn, "w")
    @write fid exx
    close(fid)
end

# Issue #226
t1 = mtime(fn)
sleep(3)
jldopen(fn, "r") do file
end
t2 = mtime(fn)
@test t2 == t1

jldopen(fn , "w") do file
    g = create_group(file, "autocorr") # create a group
    g["countsmapdistfreq"] = "abracadabra"
    g["xvect"] = 1:20
    file["stray"] = "cat"
end

a = load(fn)
@test a["autocorr"]["xvect"] == 1:20
@test a["autocorr"]["countsmapdistfreq"] == "abracadabra"
@test a["stray"] == "cat"

# Test whether we get a stack overflow with FileIO interface (#46)
@test_throws Union{ErrorException,CapturedException} save(fn, 7)

# Test load macros
jldopen(fn, "w") do file
    write(file, "loadmacrotestvar1", ['a', 'b', 'c'])
    write(file, "loadmacrotestvar2", 1)
end

@eval begin # wrapped in eval since @load with no args needs file at compile time
function f1()
    @load $fn
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end
end

f1()

function f2()
    @load fn loadmacrotestvar1 loadmacrotestvar2
    @test loadmacrotestvar1 == ['a', 'b', 'c']
    @test loadmacrotestvar2 == 1
end

f2()

@test !@isdefined loadmacrotestvar1 # should not be in global scope
@test (@eval @load $fn) == [:loadmacrotestvar1, :loadmacrotestvar2]
@test loadmacrotestvar1 == ['a', 'b', 'c']
@test loadmacrotestvar2 == 1

# Test StackFrame by saving profile output
@profile eigvals(randn(3,3))
li, lidict = Profile.retrieve()
f = tempname()*".jld"
@save f li lidict
@test isa(JLD.load(f)["lidict"], Dict{UInt64,Array{Base.StackFrame,1}})
rm(f)

# Issue #173
struct TestBenchmarkGroup
    tags::Vector{Any}
    data::Dict{Any,Any}
end
mktempdir() do d
    file = joinpath(d, "test173.jld")
    JLD.save(file, "t173", TestBenchmarkGroup(Any[], Dict{Any,Any}()))
    t = JLD.load(file)["t173"]
    @test t isa TestBenchmarkGroup
    @test isempty(t.tags) && eltype(t.tags) == Any
    @test isempty(t.data) && eltype(t.data) == Pair{Any,Any}
    rm(file)
end

# Issues #249 and #251
# TODO: debug mmap problems on Windows
if !Sys.iswindows()
    let fid = jldopen(fn, "w", mmaparrays = true)
        write(fid, "a", true)
        @test read(fid, "a") == true

        write(fid, "b", false)
        @test read(fid, "b") == false

        write(fid, "c", Bool[true, true, false, false, true])
        @test read(fid, "c") == Bool[true, true, false, false, true]

        local t = Dict{String, Any}("a" => true, "b" => false, "c" => "xyz", "d" => 13)
        write(fid, "d", t)
        @test read(fid, "d") == t

        close(fid)
        rm(fn)
    end
end

# Issues #294, #297
struct T294{T}
    a::T
end
let fid = jldopen(fn, "w")
    a = T294((1, 2))
    b = T294(a)
    c = T294(b)
    write(fid, "a", a)
    write(fid, "b", b)
    write(fid, "c", c)
    close(fid)

    @test load(fn, "a") == a
    @test load(fn, "b") == b
    @test load(fn, "c") == c
end


text:
using Test, EvalMetrics, Random
using EvalMetrics.Encodings

include("encodings.jl")
include("utilities.jl")
include("confusion_matrix.jl")
include("metrics.jl")
include("thresholds.jl")
include("curves.jl")

ext: jl


text:
"""
Add values to extras
"""
function addextras(data, result)
    data.extras[ModelSelection.generate_extra_key(CROSSVALIDATION_EXTRAKEY, data.extras)] = Dict(
        :ttest => result.ttest,
        :kfolds => result.k,
        :tsetsize => result.s,
        :panel => data.panel,
        :time => data.time,
        :datanames => result.datanames,
        :median => result.median_data,
        :average => result.average_data
    )
    return data
end


text:
"""
Estimate the covariance derivatives with forward differences
"""
function est_dΣdθ(prob_def::GLO, kernel_hyperparameters::Vector{T}; return_est::Bool=true, return_anal::Bool=false, return_dif::Bool=false, return_bool::Bool=false, dif::Real=1e-6, print_stuff::Bool=true) where {T<:Real}

    total_hyperparameters = append!(collect(Iterators.flatten(prob_def.a0)), kernel_hyperparameters)

    if print_stuff
        println()
        println("est_dΣdθ: Check that our analytical dΣdθ are close to numerical estimates")
        println("hyperparameters: ", total_hyperparameters)
    end

    x = prob_def.x_obs
    return_vec = []
    coeff_orders = coefficient_orders(prob_def.n_out, prob_def.n_dif, a=prob_def.a0)

    # construct estimated dΣdθs
    if return_est || return_dif || return_bool
        val = covariance(prob_def, total_hyperparameters)
        est_dΣdθs = zeros(length(total_hyperparameters), prob_def.n_out * length(x), prob_def.n_out * length(x))
        for i in 1:length(total_hyperparameters)
            if total_hyperparameters[i]!=0
                hold = copy(total_hyperparameters)
                hold[i] += dif
                est_dΣdθs[i, :, :] =  (covariance(prob_def, hold) - val) / dif
            end
        end
        if return_est; append!(return_vec, [est_dΣdθs]) end
    end

    # construct analytical dΣdθs
    if return_anal || return_dif || return_bool
        anal_dΣdθs = zeros(length(total_hyperparameters), prob_def.n_out * length(x), prob_def.n_out * length(x))
        for i in 1:length(total_hyperparameters)
            anal_dΣdθs[i, :, :] =  covariance(prob_def, total_hyperparameters; dΣdθs_total=[i])
        end
        if return_anal; append!(return_vec, [anal_dΣdθs]) end
    end

    if return_dif || return_bool
        difs = est_dΣdθs - anal_dΣdθs
        append!(return_vec, [difs])
    end

    # find whether the analytical and estimated dΣdθs are approximately the same
    if return_bool
        no_differences = true
        min_thres = 5e-3
        max_val = 0
        for ind in 1:length(total_hyperparameters)
            est = est_dΣdθs[ind,:,:]
            est[est .== 0] .= 1e-8
            val = mean(abs.(difs[ind,:,:] ./ est))
            max_val = maximum([max_val, val])
            if val>min_thres
                no_differences = false
                if print_stuff; println("dΣ/dθ$ind has a average ratioed difference of: ", val) end
            end
        end
        if print_stuff; println("Maximum average dΣ/dθ ratioed difference of: ", max_val) end
        return no_differences
    end

    return return_vec

end

function est_∇(f::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)
    # original value
    val = f(inputs)

    #estimate gradient
    j = 1
    grad = zeros(length(remove_zeros(inputs)))
    for i in 1:length(inputs)
        if !ignore_0_inputs || inputs[i]!=0
            hold = copy(inputs)
            hold[i] += dif
            grad[j] =  (f(hold) - val) / dif
            j += 1
        end
    end
    return grad
end

"""
Estimate the gradient of nlogL_GLOM() with forward differences

Parameters:

prob_def (GLO): A structure that holds all of the relevant
    information for constructing the model used in the GLOM et al. 2017+ paper
total_hyperparameters (vector): The hyperparameters for the GP model, including
    both the coefficient hyperparameters and the kernel hyperparameters
dif (float): How much to perturb the hyperparameters

Returns:
vector: An estimate of the gradient

"""
function est_∇nlogL_GLOM(prob_def::GLO, total_hyperparameters::Vector{T}; dif::Real=1e-7) where {T<:Real}

    f(inputs) = nlogL_GLOM(prob_def, inputs)
    return est_∇(f, total_hyperparameters; dif=dif, ignore_0_inputs=true)

end

function test_∇(est_G::Vector{T}, G::Vector{T}; print_stuff::Bool=true, function_name::String="function") where {T<:Real}

    if print_stuff
        println()
        println("test_∇: Check that our analytical $function_name is close to numerical estimates")
        # println("only values for non-zero hyperparameters are shown!")
        # println("hyperparameters: ", total_hyperparameters)
        println("analytical: ", G)
        println("numerical : ", est_G)
    end

    no_mismatch = true
    for i in 1:length(G)
        if !isapprox(G[i], est_G[i], rtol=5e-2)
            println("mismatch df/dθ" * string(i))
            no_mismatch = false
        end
    end

    return no_mismatch
end

"""
Test that the analytical and numerically estimated ∇nlogL_GLOM() are approximately the same

Parameters:

prob_def (GLO): A structure that holds all of the relevant
    information for constructing the model used in the GLOM et al. 2017+ paper
kernel_hyperparameters (vector): The kernel hyperparameters for the GP model
dif (float): How much to perturb the hyperparameters
print_stuff (bool): if true, prints extra information about the output

Returns:
bool: Whether or not every part of the analytical and numerical gradients match

"""
function test_∇nlogL_GLOM(prob_def::GLO, kernel_hyperparameters::Vector{T}; dif::Real=1e-4, print_stuff::Bool=true) where {T<:Real}

    total_hyperparameters = append!(collect(Iterators.flatten(prob_def.a0)), kernel_hyperparameters)
    G = ∇nlogL_GLOM(prob_def, total_hyperparameters)
    est_G = est_∇nlogL_GLOM(prob_def, total_hyperparameters; dif=dif)

    return test_∇(est_G, G; print_stuff=print_stuff, function_name="∇nlogL_GLOM")

end

function est_∇∇(g::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)

    val = g(inputs)

    #estimate hessian
    j = 1
    dim = length(remove_zeros(inputs))
    hess = zeros(dim, dim)
    for i in 1:length(inputs)
        if !ignore_0_inputs || inputs[i]!=0
            hold = copy(inputs)
            hold[i] += dif
            hess[j, :] =  (g(hold) - val) / dif
            j += 1
        end
    end
    return symmetric_A(hess)

end

function est_∇∇_from_f(f::Function, inputs::Vector{<:Real}; dif::Real=1e-7, ignore_0_inputs::Bool=false)

    val = est_∇(f, inputs; dif=dif, ignore_0_inputs=ignore_0_inputs)

    #estimate hessian
    j = 1
    hess = zeros(length(inputs), length(inputs))
    for i in 1:length(inputs)
        if !ignore_0_inputs || inputs[i]!=0
            hold = copy(inputs)
            hold[i] += dif
            hess[j, :] =  (est_∇(f, hold; dif=dif, ignore_0_inputs=ignore_0_inputs) - val) / dif
            j += 1
        end
    end
    return symmetric_A(hess)

end

"""
Estimate the Hessian of nlogL_GLOM() with forward differences

Parameters:

prob_def (GLO): A structure that holds all of the relevant
    information for constructing the model used in the GLOM et al. 2017+ paper
total_hyperparameters (vector): The hyperparameters for the GP model, including
    both the coefficient hyperparameters and the kernel hyperparameters
dif (float): How much to perturb the hyperparameters

Returns:
matrix: An estimate of the hessian
est_∇∇nlogL_GLOM
"""
function est_∇∇nlogL_GLOM(prob_def::GLO, total_hyperparameters::Vector{T}; dif::Real=1e-4) where {T<:Real}

    g(inputs) = ∇nlogL_GLOM(prob_def, inputs)
    return est_∇∇(g, total_hyperparameters; ignore_0_inputs=true)

end

function test_∇∇(est_H::Union{Symmetric{T,Matrix{T}},Matrix{T}}, H::Union{Symmetric{T,Matrix{T}},Matrix{T}}; print_stuff::Bool=true, function_name::String="function", rtol::Real=1e-3) where {T<:Real}

    if print_stuff
        println()
        println("test_∇∇: Check that our analytical $function_name is close to numerical estimates")
        # println("only values for non-zero hyperparameters are shown!")
        # println("hyperparameters: ", total_hyperparameters)
        println("analytical:")
        for i in 1:size(H, 1)
            println(H[i, :])
        end
        println("numerical:")
        for i in 1:size(est_H, 1)
            println(est_H[i, :])
        end
    end

    no_mismatch = true
    matches = fill(1, size(H))
    for i in 1:size(H, 1)
        for j in 1:size(H, 2)
            if !(isapprox(H[i, j], est_H[i, j], rtol=rtol))
                # println("mismatch d2nlogL/dθ" * string(i) * "dθ" * string(j))
                matches[i, j] = 0
                no_mismatch = false
            end
        end
    end

    if !no_mismatch
        println("mismatches at 0s")
        for i in 1:size(H, 1)
            println(matches[i, :])
        end
    end

    return no_mismatch
end

"""
Test that the analytical and numerically estimated ∇∇nlogL_GLOM() are approximately the same

Parameters:

prob_def (GLO): A structure that holds all of the relevant
    information for constructing the model used in the GLOM et al. 2017+ paper
kernel_hyperparameters (vector): The kernel hyperparameters for the GP model
dif (float): How much to perturb the hyperparameters
print_stuff (bool): if true, prints extra information about the output

Returns:
bool: Whether or not every part of the analytical and numerical Hessians match

"""
function test_∇∇nlogL_GLOM(prob_def::GLO, kernel_hyperparameters::Vector{T}; dif::Real=1e-4, print_stuff::Bool=true) where {T<:Real}

    total_hyperparameters = append!(collect(Iterators.flatten(prob_def.a0)), kernel_hyperparameters)
    H = ∇∇nlogL_GLOM(prob_def, total_hyperparameters)
    est_H = est_∇∇nlogL_GLOM(prob_def, total_hyperparameters; dif=dif)

    return test_∇∇(est_H, H; print_stuff=print_stuff, function_name="∇∇nlogL_GLOM")

end


text:
"""
The inverse logit function, which is also the CDF of the Logistic(0, 1)
distribution. Note that this function will generate an exact value of 0 when
z <= -710 and will generate an exact value of 1 when z >= 37. This introduces
numeric problems if you provide inputs outside of that range.
"""
@inline invlogit(z::Real) = 1 / (1 + exp(-z))

"""
The logit function, which is also the quantile function of the Logistic(0, 1)
distribution.
"""
@inline logit(p::Real) = log(p / (1 - p))


text:
using DSP
using WAV
function AD_chebyshev(order,x)

    return-(x*cos(order*acos(x))+order*sqrt(1-x^2)*sin(order*acos(x)))/(-1+order^2)

end
#=
this code is based on this post
https://jatinchowdhury18.medium.com/practical-considerations-for-antiderivative-anti-aliasing-d5847167f510

this is very similar reworking of his code for my purposes
=#

#anti derivative anti aliasing order 1
function processAD1(signal::Vector{Float64},f,f1, tolerance = 1.0e-5)

    x1 = 0.0
    returnSignal = Float64[]
    for s in signal
        if abs(s - x1) < tolerance
            push!(returnSignal, f((s + x1) / 2))
        else
            push!(returnSignal, (f1(s) - f1(x1)) / (s-x1))
            x1 = s
        end
    end
    return returnSignal
end
#anti derivative anti aliasing order 2
function processAD2(signal, f, f1, f2, ε=1.0e-5)

    returnSignal = Float64[]

    calcD(x0, x1) = abs(x0-x1) < ε ? f1((x0+x1)/2) : (f2(x0) - f2(x1)) / (x0-x1)

    x₁ = 0.0
    x₂ = 0.0
    for x₀ in signal
        if abs(x₀ - x₁) < ε
            x̄ = (x₀ + x₂)/2
            Δ = x̄ - x₀
            if Δ < ε
                push!(returnSignal,  f((x̄  + x₀) / 2.0))
            else
                push!(returnSignal,  (2/Δ) * (f1(x̄) + (f2(x₀) - f2(x̄)) / Δ))
            end
        else
           push!(returnSignal, (2.0 / (x₀ - x₂)) * (calcD(x₀, x₁) - calcD(x₁, x₂)) )
           x₂ = x₁
           x₁ = x₀
        end
    end

    return returnSignal
end

#hardclip first antiderivatie
function hc1(x)
    return x*sign(x) - .5
end
#hardclip second order derivative
function hc2(x)
    return (x^2/2 + 1/6) * sign(x) - x/2
end

function sc(x::Float64)::Float64
    return 2*atan(x)/pi
end
#softClip first order anti derivative
function sc1(x::Float64)::Float64
    return (2x*atan(x)-log(x^2+1))/pi
end
#softClip second order anti derivatie
function sc2(x::Float64)::Float64
    return (x*(-log(x^2+1))+(x^2-1)*(atan(x))+x)/pi
end
function hardC(signal::Vector{Float64}, gain)::Vector{Float64}
    c(x,y) = (abs(x*y) < 1 && return x*y)  || return sign(x)
    return c.(signal,gain)
end
function hardClipAD1(signal::Vector{Float64}, gain)::Vector{Float64}
    f(x) = (abs(x) < 1 ? x : sign(x) )
    f1(x) = (abs(x) < 1 ? x^2/2 : x*sign(x) - .5)
    return processAD1(gain*signal, f,f1)
end

function hardClipAD2(signal::Vector{Float64}, gain)::Vector{Float64}
    f(x) = (abs(x) < 1 ? x : sign(x) )
    f1(x) = (abs(x) < 1 ? x^2/2 : x*sign(x) - .5)
    f2(x) = (abs(x) < 1 ? x^3/6 : (x^2/2 + 1/6) * sign(x) - x/2)
    return processAD2(gain*signal, f,f1,f2)
end
function softClipAD2(signal::Vector{Float64}, gain)::Vector{Float64}
    return processAD2(gain*signal, sc,sc1,sc2)
end

max_stars_repo_name: justinbroce/julia-audio-stuffs

max_issues_repo_name: justinbroce/julia-audio-stuffs

max_forks_repo_name: justinbroce/julia-audio-stuffs


text:
using WaveletMatrices

@inline function raw_getindex(seq::BioSequence{A}, i::Integer) where A
    j = BioSequences.bitindex(seq, i)
    @inbounds return convert(UInt8, (seq.data[BioSequences.index(j)] >> BioSequences.offset(j)) & BioSequences.mask(A))
end

## This function is necessary as a result of changes made in Bio.jl as of Julia v0.5
## Since all nucleotides are converted to UInt8 as one hot 4-bit encodings by default
## we can override FMIndexes.sa_range function for those sequences specifically
## to re-encode as DNAAlphabet{2}
@inline function FMIndexes.sa_range(query::BioSequence{A}, index::FMIndex, init_range::UnitRange{Int}) where A <: BioSequences.Alphabet
    sp, ep = init_range.start, init_range.stop
    # backward search
    i = length(query)
    @inbounds while sp ≤ ep && i ≥ 1
        char = raw_getindex(query, i)
        c = index.count[char+1]
        sp = c + WaveletMatrices.rank(char, index.bwt, (sp > index.sentinel ? sp - 1 : sp) - 1) + 1
        ep = c + WaveletMatrices.rank(char, index.bwt, (ep > index.sentinel ? ep - 1 : ep))
        i -= 1
    end
    return sp:ep
end


text:
#=
testfftbig:
- Julia version:
- Author: ymocquar
- Date: 2019-11-18
=#

using Test
using FFTW
using Random
using LinearAlgebra
using HOODESolver

function getfalse( tab )
    for i=1:size(tab,1), j=1:size(tab,2)
        if !tab[i,j]
            return i,j
        end
    end
    return 0, 0
end

function testfftbig( s, T::DataType, seed_val )

    Random.seed!(seed_val)
    tab = zeros(Complex{T}, 1, s)
    tab .= rand(T, 1, s)
    if T == Float64
        tabfftref = fft(tab,2)
    else
        tab2 = zeros(Complex{Float64}, 1, s)
        tab2 .= tab
        tabfftref = fft(tab2,2)
    end

    tab_test = copy(tab)

    p = HOODESolver.PrepareFftBig(s, real(tab[1, 1]))

    HOODESolver.fftbig!(p, tab_test)

    @test isapprox(tabfftref, tab_test, atol=1e-15, rtol=1e-15)

    tol = (T == BigFloat) ? 1e-50 : 1e-15

    HOODESolver.fftbig!(p, tab_test, flag_inv=true)

    @test getfalse(isapprox.(tab, tab_test, atol=tol, rtol=tol)) == (0, 0)

    @test isapprox(tab, tab_test, atol=tol, rtol=tol)

end
function testfftbig2( s, T::DataType, seed_val, nb_v )

    Random.seed!(seed_val)
    tab = zeros(Complex{T}, nb_v, s)
    tab .= rand(T, nb_v, s)
    if T == Float64
        tabfftref = fft(tab,2)
    else
        tab2 = zeros(Complex{Float64}, nb_v, s)
        tab2 .= tab
        tabfftref = fft(tab2,2)
    end

    tab_test = copy(tab)

    p = HOODESolver.PrepareFftBig(s, one(T) )

    tab_test2 = HOODESolver.fftbig(p, tab_test)

    @test isapprox(tabfftref, tab_test2, atol=1e-15, rtol=1e-15)

    HOODESolver.fftbig!(p, tab_test)

    @test isapprox(tabfftref, tab_test, atol=1e-15, rtol=1e-15)

    tol = (T == BigFloat) ? 1e-50 : 1e-15

    tab_test3 = HOODESolver.fftbig(p, tab_test, flag_inv=true)
    @test isapprox(tab, tab_test3, atol=tol, rtol=tol)

    HOODESolver.fftbig!(p, tab_test, flag_inv=true)

    @test isapprox(tab, tab_test, atol=tol, rtol=tol)

end

function testfftbigprec(s, seed_v)

    @time @testset "test fftbig changing precision s=$s seed=$seed_v" begin
        Random.seed!(seed_v)
        for k=1:10
            setprecision(k*256) do
                p = HOODESolver.PrepareFftBig(s)
                tab = rand(BigFloat, 5, s)
                tabRef = Complex.(tab)
                tabRes = HOODESolver.fftgen(p,tab)
                tabRes2 = HOODESolver.ifftgen(p, tabRes)
                println("k=$k norm=$(norm(tabRef-tabRes2))")
#                @test isapprox(tabRef, tabRes2, atol=1e-78^k)
#                @test isapprox(real.(tabRef), real.(tabRes2), atol=1e-78^k)
                @test isapprox(tabRef, tabRes2, atol=1e-75^k)
                @test isapprox(real.(tabRef), real.(tabRes2), atol=1e-75^k)
            end
        end
    end

end

tab_decl = [[8, Float64, 12345678], [8, BigFloat, 9876],[2^8, Float64, 1928], [2^10, BigFloat, 5656]]

for t in tab_decl
    s = t[1]
    type = t[2]
    seed_v = t[3]
    @time @testset "test fftbig for value size=$s type=$type" begin testfftbig(s, type, seed_v) end
end
tab_decl2 = [[8, Float64, 4556789, 4], [8, BigFloat, 4563, 4],[2^10, Float64, 9900771, 4], [2^13, BigFloat, 125609, 4]]

for t in tab_decl2
    s = t[1]
    type = t[2]
    seed_v = t[3]
    nb_v = t[4]
    @time @testset "test fftbig for vector size=$s type=$type nb_v=$nb_v" begin testfftbig2(s, type, seed_v, nb_v) end
end

tab_decl3 =[ 8965, 1919191, 56188827, 9017555]
for sd in tab_decl3
    testfftbigprec(32,sd)
end


text:
#=
    Contents: Utils for handling authentication.
=#

"""
Returns true if username and password are valid.

INPUT:
- acdata: Either the store of Access Control data, or a connection to it.
          For example, an intra-process Dict or a connection to a database from a connection pool.
"""
function login_credentials_are_valid(username::AbstractString, password::AbstractString, acdata)
    salt, hashed_pwd = get_salt_hashedpwd(acdata, username)
    password_is_valid(password, salt, hashed_pwd)
end

"Returns: username and password extracted from login POST request."
function extract_username_password(req)
   qry = bytestring(req.data)        # query = "username=xxx&password=yyy"
   dct = JSON.parse(qry)             # Dict("username" => xxx, "password" => yyy)
   dct["username"], dct["password"]
end

"""
Returns: username, current password and new password extracted from process_pwdreset POST request.

Notes:
1. New password is supplied twice, via new_pwd and new_pwd2.
"""
function extract_currpwd_newpwd_newpwd2(req)
    qry = bytestring(req.data)
    dct = JSON.parse(qry)    # Dict("form_id" => "pwdreset", "nonce" => nonce, "current_pwd" => yyy, "new_pwd" => zzz, "new_pwd2" => www)
    dct["form_id"], dct["nonce"], dct["current_pwd"], dct["new_pwd"], dct["new_pwd2"]
end

"Returns true if req has a valid 'sessionid' cookie."
function is_logged_in(req::Request)
    result     = false
    cookiename = config[:session][:cookiename]
    sessions   = config[:session][:datastore]
    if sessions == :cookie
	result = read_session(req, cookiename) != ""
    else
	session_id = read_sessionid(req, cookiename)           # Read session_id from cookie
	result     = session_is_valid(sessions, session_id)    # Check session_id against server-side session data
    end
    result
end

is_not_logged_in(req) = !is_logged_in(req)

"Returns: username, derived from request regardless of whether the session is client-side or server-side."
function get_username(req::Request)
    username = ""
    if config[:session][:datastore] == :cookie
	session  = session_read(req)
	username = session["username"]
    else
	session_id = read_sessionid(req)
	username   = session_get(session_id, "username")
    end
    username
end

"""
Returns: Redirect location according to config.

The config either specifies:
- A location. Example: "/home", or
- A function of the user: Example: f(username) = '/members/get_role(username)'
"""
function get_redirect_location(username::AbstractString, config_key1::Symbol, config_key2::Symbol)
    loc = ""
    _redirect = config[config_key1][config_key2]
    if typeof(_redirect) == Function
	loc = _redirect(username)
    else
	loc = _redirect
    end
    loc
end

# EOF


text:
projectdir() = "C:\\Users\\nritchie\\Documents\\DrWatson\\K-Ratio Project"

using NeXLDatabase
using DataFrames
using Gadfly
using Weave

krdPath = joinpath(projectdir(), "data","exp_raw")
krdb = joinpath(projectdir(), "Data", "exp_pro", "kratio.db")
db = openNeXLDatabase(krdb)
if false
  dale = read(db,DBPerson,"dale.newbury@nist.gov")
  nicholas = read(db, DBPerson, "nicholas.ritchie@nist.gov")
  analyst=dale
  krp = read(db, DBProject, 0, "K-ratio Project")
  lab = findall(db,DBLaboratory, analyst)[1]
  jeol = read(db, DBInstrument, write(db, DBInstrument, lab, "JEOL", "JXA-8500F", "217 D113"))
  tescan = read(db, DBInstrument, write(db, DBInstrument, lab, "TESCAN", "MIRA-3", "217 F101"))

  brukersingle = read(db, DBDetector, write(db, DBDetector, jeol, "Bruker", "XFlash 6|30", "Bruker Single", 128.9, 121, -483.9, 5.0163, z(n"Be"), z(n"Sc"), z(n"Sn"), z(n"Pu")))
  brukerquad = read(db, DBDetector, write(db, DBDetector, jeol, "Bruker", "XFlash QUAD 5040", "Bruker Quad", 131.0, 121, -478.0, 5.003, 4, 21, 55, 94))
  for ii in 0:3
      write(db, DBDetector, tescan, "Pulsetor", "Torrent", "Torrent $ii", 131.0, 20, 0.0, 10.0, 4, 21, 55, 94)
  end
  write(db, DBDetector, tescan, "Pulsetor", "Torrent", "Torrent (combined)", 131.0, 20, 0.0, 10.0, 4, 21, 55, 94)
  pulsetor = findall(db, DBDetector, tescan, vendor="Pulsetor",model="Torrent")

  blockC = read(db, DBSample, write(db, DBSample, lab, "Standard Mount C", "NIST's Standard Mount C - Mostly pure elements"))
  blockElm = read(db, DBSample, write(db, DBSample, lab, "Element Mount I", "NIST's Element Mount I - Mostly pure elements"))
  blockElm2 = read(db, DBSample, write(db, DBSample, lab, "Element Mount II", "NIST's Element Mount I - Mostly pure elements"))
  blockHiTC = read(db, DBSample, write(db, DBSample, lab, "HighTC Block", "NIST's High Temperature Superconductor Block"))
  blockChMixed = read(db, DBSample, write(db, DBSample, lab, "Chuck's Mixed", "NIST's Chuck's Mixed Standards Block"))
  blockGM3 = read(db, DBSample, write(db, DBSample, lab, "Glass Mount III", "NIST's Glass Mount III"))
  blockUnk = read(db, DBSample, write(db, DBSample, lab, "Unknown", "Unknown"))
  blockFe3Si = read(db, DBSample, write(db, DBSample, lab, "??Fe3Si??", "SRM-483 block in ???"))
  blockAuAgCu = read(db, DBSample, write(db, DBSample, lab, "SRM-481 & 482 block", "SRM-481 & 482 block in ???"))
  alphaAesar = read(db, DBSample, write(db, DBSample, lab, "Alfa-Aesar Sample", "Sourced from Alfa-Aesar"))
  maurine = read(db, DBSample, write(db, DBSample, lab, "Maureen William's Sample", "Co alloys"))
  colby = read(db, DBSample, write(db, DBSample, lab, "Colby's Fe-Cr alloys", "Fe-Cr alloys"))
  geller583 = read(db, DBSample, write(db, DBSample, lab, "Geller 583","Geller Low-Z APA Mount"))
  geller588 = read(db, DBSample, write(db, DBSample, lab, "Geller 588","Geller Higher-Z APA Mount"))
  mmm1 = read(db, DBSample, write(db, DBSample, lab, "Mengason Mineral Mount 1", "Mengason Mineral Mount 1"))
  mmm2 = read(db, DBSample, write(db, DBSample, lab, "Mengason Mineral Mount 2", "Mengason Mineral Mount 2"))

  carbonCoating, noCoating= [ n"C" ], Element[]

  sulfides10=joinpath(krdPath,"Sulfides_NewBruker","10_keV_Sulfides")
  project = read(db, DBProject, write(db, DBProject, "Sulfides", "Sulfides - CuS, FeS, FeS2, CdS, PbS, ZnS", analyst, krp))
  sample = blockC
  e0 = 10.0e3

  unks = [
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_1_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_2_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_3_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_4_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_5_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_6_10kV10nA130kHz_62kHz12DT_100s.msa"),
    joinpath(sulfides10,"EDITED_FeS_StdC_10keV","FeS_StdC_7_10kV10nA130kHz_62kHz12DT_100s.msa"),
  ]

  stds = [
    #(blockChMixed, mat"CdSe", joinpath(sulfides10,"EDITED_Standards_10keV","CdSe_ChMxd_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockChMixed, mat"CdS", joinpath(sulfides10,"EDITED_Standards_10keV","CdS_ChMxd_10kV10nA130kHz_700s.msa"), e0, noCoating ),
    #(blockChMixed, mat"Cd", joinpath(sulfides10,"EDITED_Standards_10keV","Cd_ChMxd_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockC, mat"Cd", joinpath(sulfides10,"EDITED_Standards_10keV","Cd_StdC_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockChMixed, mat"CuS", joinpath(sulfides10,"EDITED_Standards_10keV","CuS_ChMxd_10kV10nA130kHz_700s.msa"), e0, noCoating ),
    #(blockC, mat"Cu", joinpath(sulfides10,"EDITED_Standards_10keV","Cu_StdC_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockC, mat"C", joinpath(sulfides10,"EDITED_Standards_10keV","C_StdC_10kV10nA130kHz_53kHz10DT_200s.msa"), e0, noCoating ),
    #(blockChMixed, mat"CdS", joinpath(sulfides10,"EDITED_Standards_10keV","EDITED_CdS_ChMxd_1_10kV10nA130kHz_78kHz15DT_100s.msa"), e0, noCoating ),
    (blockChMixed, mat"FeS2", joinpath(sulfides10,"EDITED_Standards_10keV","FeS2_ChMxd_10kV10nA130kHz_700s.msa"), e0, noCoating ),
    (blockChMixed, mat"Fe", joinpath(sulfides10,"EDITED_Standards_10keV","Fe_ChMxd_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockC, mat"Fe," joinpath(sulfides10,"EDITED_Standards_10keV","Fe_StdC_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockC, mat"MgO", joinpath(sulfides10,"EDITED_Standards_10keV","MgO_StdC_10kV10nA130kHz_81kHz15DT_100s.msa"), e0, noCoating ),
    #(blockC, mat"PbSe", joinpath(sulfides10,"EDITED_Standards_10keV","PbSe_StdC_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockChMixed, mat"PbTe", joinpath(sulfides10,"EDITED_Standards_10keV","PbTe_ChMxd_10kV10nA130kHz_300s.msa"), e0, noCoating ),
    #(blockC, mat"Se", joinpath(sulfides10,"EDITED_Standards_10keV","Se_StdC_10kV10nA130kHz_107kHz20DT_100s.msa"), e0, noCoating ),
    #(blockChMixed, mat"Te", joinpath(sulfides10,"EDITED_Standards_10keV","Te_ChMxd_10kV10nA130kHz_74kHz14DT_100s.msa"), e0, noCoating ),
    #(blockC, mat"Zn", joinpath(sulfides10,"EDITED_Standards_10keV","Zn_StdC_10kV10nA130kHz_300s.msa"), e0, noCoating ),
  ]

  cfs = write(db, DBCampaign, project, sample, unkComp, brukersingle, analyst, e0, unks, stds, noCoating)
else
  cfs = 85
end

unkComp = parse(Material, "FeS", pedigree="Stoichiometric")
ffrs = fit_spectrum(db, DBCampaign, cfs, unkComp)
display(asa(DataFrame, ffrs))
display(plot(ffrs[1]))
#display(plot(loadspectrum.(unks)...,klms=[ n"C", keys(unkComp)... ], xmax=10.0e3))


text:
config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 0.0],
    (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0]],
    (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0]])
for mock in mocks(optimize!)
    Tests.sos_term_fixed_test(mock, config)
    Tests.sdsos_term_fixed_test(mock, config)
    Tests.dsos_term_fixed_test(mock, config)
end


text:
module Kwesi
export *, VarNode, OpsNode, runNode
#need to define new method of * to handle single element arrays of Real
function *(x::Array{Real}, y::Array{Real})
  x2 = (x[1])::Real
  y2 = (y[1])::Real
  x2 * y2
end
#base type of node
abstract BaseNode
#VarNode defines variables that simply return their state value
type VarNode <: BaseNode
  name::String
  state::Array{Real, 1}
  function VarNode(name="VarNode", state=[0])
    this = new()
    this.name = name
    this.state = state
    this
  end
end
#OpsNodes return the result from their operation (a defined function)
type OpsNode <: BaseNode
  operation::Function
  state::Array{Real, 1}
  inputs::Array{BaseNode}
  function OpsNode{T<:BaseNode}(operation, state, inputs::Array{T, 1})
    this = new()
    this.operation = operation
    this.state = state
    this.inputs = inputs
    this
  end
end

function runNode(node::VarNode)
  return node.state
end

function runNode(node::OpsNode)
  inputs = [runNode(input) for input in node.inputs]
  state::Array{Real, 1}
  state = Real[]
  result = apply(node.operation, inputs)
  if (typeof(result) <: Array)
    state = result
  elseif (typeof(result) <: Real)
    push!(state, result)
  end
  node.state = state
end
end

max_stars_repo_name: outlace/Kwesi

max_issues_repo_name: outlace/Kwesi

max_forks_repo_name: outlace/Kwesi


text:
using Plots:plot, heatmap
using WAV:wavread

fnlist = [
barthann,
bartlett,
blackman,
blackmanharris,
bohman,
flattop,
hamming,
hanning,
nuttall,
parzen,
rectangular,
triangular]

plts = []
data, Fs = wavread("AMonoAudioFile.wav");
for winfn in fnlist
    fbank = MelFilters(
          fs = floor(Int,Fs),
       alpha = 0.5,
      winlen = 1024,
      stride = 128,
      nbanks = 128,
      nffts  = 1024,
     epsilon = 1e-6,
     winfunc = winfn)

    feat = fbank(copy(data));
    temp = heatmap(feat, legend=nothing, framestyle=:box, ticks=nothing, titlefont=(20, "times"))
    title!("$winfn")
    push!(plts, temp)
end

plot(plts...,layout=(4,3),size=(1600,1200))


text:
#testbug_GB18
include("GB19.jl")
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV, Random, DataFrames, Gadfly

rng = Random.seed!(2019)

gbdt = GB18.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                 # Learning rate
  num_iterations = 100,                # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB18.GBLearner(gbdt, :regression)

x, y = GB18.test_data(500; p = 10, sigma = sqrt(1))
x_test, y_test = GB18.test_data(500; p = 10, sigma = sqrt(1))
mu = GB18.fried(x)
mu_test = GB18.fried(x_test)

# initialize the hypers
hypers = GB18.Hypers(x, y, 1, 1, 0.8, 2)
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
hypers.α = 0.95
hypers.γ = 2.0
hypers.s =  ones(size(x, 2)) / size(x, 2)
hypers

# initialize opt
opts = GB18.Opts()
opts.num_burn = 5

rmse2, Store_s = GB18.TBD_predict_and_rmse!(
    x,
    y,
    x_test,
    y_test,
    mu_test,
    gbdt,
    hypers,
    opts,
    rng)

hypers
sum(hypers.varcount) # 2765(burn 460)
sum(hypers.varcount[1:5]) # 1218(burn 460)
sum(hypers.s[1:5])

hypers.s

for k in 1:500
    GB18.UpdateAlpha!(hypers)
end
hypers
alpha_current = hypers.α
hypers.s .!= 0.0

current_s     = hypers.s[hypers.s .!= 0.0]

function alpha_to_rho(alpha::Float64, scale::Float64)
    return alpha / (alpha + scale)
end

function rho_to_alpha(rho::Float64, scale::Float64)
    return scale * rho / (1 - rho)
end
hypers.λₐ
rho_current = alpha_to_rho(alpha_current, hypers.λₐ)

current_s = hypers.s[hypers.s .> 0.]

psi = mean([log(current_s[i]) for i in 1:length(current_s)])
p = 1.0 * length(current_s)

loglik_current = alpha_current * psi + lgamma(alpha_current) -
    p * lgamma(alpha_current / p) +
    logpdf(Beta(hypers.αₐ, hypers.αᵦ), rho_current)

rho_propose = rand(Beta(hypers.αₐ, hypers.αᵦ))
alpha_propose = rho_to_alpha(rho_propose, hypers.α)

loglik_propose = alpha_propose * psi + lgamma(alpha_propose) -
    p * lgamma(alpha_propose / p) +
        logpdf(Beta(hypers.αₐ, hypers.αᵦ), rho_propose)

log(rand()) < loglik_propose - loglik_current

alpha_current = alpha_propose
rho_current = rho_propose
loglik_current = loglik_propose

hypers.α = alpha_current

hypers

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text: println("Julia 👩🦰 sagt: Hallo, Jörg!")

max_stars_repo_name: kantel/julia

max_issues_repo_name: kantel/julia

max_forks_repo_name: kantel/julia


text:
# Copyright 2020 Digital Domain 3.0
#
# Licensed under the Apache License, Version 2.0 (the "Apache License")
# with the following modification; you may not use this file except in
# compliance with the Apache License and the following modification to it:
# Section 6. Trademarks. is deleted and replaced with:
#
# 6. Trademarks. This License does not grant permission to use the trade
#    names, trademarks, service marks, or product names of the Licensor
#    and its affiliates, except as required to comply with Section 4(c) of
#    the License and to reproduce the content of the NOTICE file.
#
# You may obtain a copy of the Apache License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the Apache License with the above modification is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the Apache License for the specific
# language governing permissions and limitations under the Apache License.

"""
Simple k-vectors formed from the wedge product starting with orthogonal 1-vector basis.
Simple bivectors, trivectors, etc.

Types for blades are generated as a unique type for each different k-vector with acending indices
i.e. e₁, e₂, e₁₂

Equipped with a metric such that eᵢ⋅eᵢ = {1,-1,0} depending on if the metric is +,-, or 0

product of parallel basis
eᵢ*eⱼ == eᵢ⋅eⱼ when i=j

product of orthogonal basis
eᵢ*eⱼ == eᵢ∧eⱼ == -eⱼ*eᵢ when i≠j

operators:
* geometric product
∧ wedge/exterior product
-+ subraction, addition between blades of same type (span the same subspace )
!,⟂ orthogonal complement
~ reverse
⋆ Hodge star
"""
module Blades

export
@generate_basis,
e₀, e₊, e₋, e₊₋₀, e⁺, e⁻, e⁰, eₓʸ, Blade,
basis_1blades,
basis_kblades,
pseudoscalar,
∧,
quadmetric,
reverse,
grade,
⟂,
orthogonal_complement,
dual,
nonmetric_dual,
subspace,
lower,
raise,
alle,
alld,
show_basis,
kform,
magnitude,
scalar,
⋆,
factor,
untype

using Combinatorics
using Base.Iterators
using FunctionWrappers: FunctionWrapper
using LinearAlgebra

partial = (f::Function,y...)->(z...)->f(y...,z...)
swap(f::Function) = (a,b...)->f(b...,a)
second(c) = c[2]

gauss_sum(first, last) = ((last-first+1)/2)*(first+last)
ishetro(c) = !mapreduce(ab->ab[1]==ab[2],&,zip(drop(c,1),take(c,length(c)-1)))

mirror(g) = vcat(g,reverse(g,dims=1))
prefix(g) = hcat(vcat(zeros(Int,size(g)[1]>>1), ones(Int,size(g)[1]>>1)), g)
graycodes( g = [0, 1] ) = prefix(mirror(g))
"""
  gray codes of length n are the rows of the returned matrix
"""
graycodes(n::Int) = n == 1 ? [0,1] : mapreduce(i->graycodes, (gc,gcf)->gcf(gc), 1:n-2; init = graycodes())

abstract type Blade{T,K} <: Number end
abstract type e₊₋₀{T} <: Blade{T,1} end
abstract type e⁺⁻⁰{T} <: Blade{T,1} end
abstract type e₊{T} <: e₊₋₀{T} end
abstract type e₋{T} <: e₊₋₀{T} end
abstract type e₀{T} <: e₊₋₀{T} end
abstract type e⁺{T} <: e⁺⁻⁰{T} end
abstract type e⁻{T} <: e⁺⁻⁰{T} end
abstract type e⁰{T} <: e⁺⁻⁰{T} end
abstract type eₓʸ{T} <: Blade{T,1} end

quadmetric(::Type{E}) where {T, E<:e₊{T}} = one(T)
quadmetric(::Type{E}) where {T, E<:e₋{T}} = -one(T)
quadmetric(::Type{E}) where {T, E<:e₀{T}} = zero(T)
quadmetric(::Type{E}) where {T, E<:e⁺{T}} = one(T)
quadmetric(::Type{E}) where {T, E<:e⁻{T}} = -one(T)
quadmetric(::Type{E}) where {T, E<:e⁰{T}} = zero(T)

tosub(i) = [ '₁','₂','₃','₄','₅','₆','₇','₈','₉' ][i]
tosup(i) = [ '¹','²','³','⁴','⁵','⁶','⁷','⁸','⁹' ][i]

export ZForm

"""
    ZForm{D}(f,df)

Differential 0-form
Represents a scalar function to be evaluated at a point.  Takes a vector returns a scalar.
df is a list of partial
The index is that of the coordinate basis 1-form
D is the dimensionality of the coordinate basis 1-forms

Example:
julia> dist(x,y) = ZForm(()->(sqrt(x^2+y^2)), (i)->[x/sqrt(x^2+y^2), y/sqrt(x^2+y^2)][i])

julia> dist(1.0,-1.0).df(2)
-0.7071067811865475

julia> dist(1.0,-1.0).df(1)
0.7071067811865475

julia> dist(1.0,-1.0).f()
1.4142135623730951
"""
struct ZForm{D} <: Number  # calling a form a number is a stretch, but it works nicely inside a blade
  f # ::Function
  gradient # ::Vector{ZForm}
  hessian # ::Array{ZForm,2}
end

const KForm{D} = Blade{ZForm{D}}

"AlwaysZero(): object that always returns zero when indexed"
struct AlwaysZero
end

import Base.getindex

Base.length(x::AT) where {AT<:AlwaysZero} = 0
getindex(x::AT, i...) where {AT<:AlwaysZero} = 0
isazero(a) = false
isazero(a::T) where T<:Number = iszero(a)
isazero(b::B) where {T, B<:Blade{T}} = isazero(b.x)
isazero(z::Z) where {D,Z<:ZForm{D}} = isazero(z.f)
isazero(a::A) where {A<:AlwaysZero} = true
Base.similar(a::A) where {A<:AlwaysZero} = A

constantfun(x) = (y...)->x
zerofun = (y,z...)->zero(y)
ZForm{D}() where D = ZForm{D}((y,z...)->one(y), AlwaysZero(), AlwaysZero())
ZForm{D}(x::F) where {D,F<:Function} = ZForm{D}(x, AlwaysZero(), AlwaysZero())
ZForm{D}(f::F, df) where {D,F<:Function} = ZForm{D}(f, df, AlwaysZero())
ZForm{D}(x::T) where {D,T<:Real} = x

Base.length(zf::ZF) where {D,ZF<:ZForm{D}} = D

Base.zero(z::ZForm{D}) where D = 0.0
Base.zero(z::Type{ZForm{D}}) where D = 0.0
Base.iszero(z::ZForm{D}) where D = typeof(z.f) <: AlwaysZero
Base.iszero(b::B) where {D, F<:ZForm{D}, B<:Blade{F}} = typeof(b.x.f) <: AlwaysZero
Base.abs(z::K) where {D, K<:ZForm{D}} = z

mapij(f, A::AT) where {AT<:AlwaysZero} = A
mapij(f, A::AT, B) where {AT<:AlwaysZero} = A

function mapij(f, A)
  Af = similar(A)
  for ij in CartesianIndices(A)
    Af[ij] = f(Tuple(ij), A[ij])
  end
  Af
end

function mapij(f, A, B)
  ABf = similar(A)
  for ij in CartesianIndices(A)
    ABf[ij] = f(Tuple(ij), A[ij], B[ij])
  end
  ABf
end

# with second derivatives
function Base.:*(a::ZA, b::ZB) where {D,ZA<:ZForm{D}, ZB<:ZForm}
  ZForm{D}((k...)->a.f(k...)*b.f(k...),
        map(i->(k...)->a.gradient[i](k...)*b.f(k...) +
                       a.f(k...)*b.gradient[i](k...),
            1:D),
        mapij((ij, addfᵢⱼ, bddfᵢⱼ)->(k...)->
                        addfᵢⱼ(k...)*b.f(k...) +
                        a.gradient[ij[1]](k...)*b.gradient[ij[2]](k...) +
                        a.gradient[ij[2]](k...)*b.gradient[ij[1]](k...) +
                        a.f(k...)*bddfᵢⱼ(k...),
            a.hessian, b.hessian))
end

function Base.:*(s::R, b::ZB) where {R<:Real, D, ZB<:ZForm{D}}
  if iszero(s)
    s
  else
    ZForm{D}((k...)->s*b.f(k...),
          isazero(b.gradient) ? b.gradient :
                         map(dfᵢ->(k...)->s*dfᵢ(k...), b.gradient),
          isazero(b.hessian) ? b.hessian :
                         map(bddfᵢⱼ->(k...)->s*bddfᵢⱼ(k...), b.hessian))
  end
end

Base.:*(a::ZA, s::R) where {ZA<:ZForm, R<:Real} = s*a

function Base.:+(a::ZA, b::ZB) where {D, ZA<:ZForm{D}, ZB<:ZForm{D}}
  ZForm{D}((u...)->a.f(u...) + b.f(u...),
           map(i->(u...)->a.gradient[i](u...)+b.gradient[i](u...), 1:D),
           mapij((ij, addfᵢⱼ, bddfᵢⱼ)->(u...)->addfᵢⱼ(u...)+bddfᵢⱼ(u...), a.hessian, b.hessian))
end

Base.:+(a::ZA, s::R) where {D, ZA<:ZForm{D}, R<:Real} = iszero(s) ? a : ZA((u...)->a.f(u...)+s, a.gradient, a.hessian)

Base.:*(a::ZA, b::B) where {D, ZA<:ZForm{D}, B<:Blade} = untype(b)(b.x*a)
#!me Base.:*(a::BZ, b::B) where {Z<:ZForm{D}, BZ<:Blade{Z}, T<:Real, B<:Blade{T}} = (raise(untype(a)(one(T)))*b)
∧(a::ZA, b::B) where {D, ZA<:ZForm{D}, B<:Blade} = a*b
∧(b::B, a::ZA) where {D, B<:Blade, ZA<:ZForm{D}} = a*b
∧(a::ZA, b::ZB) where {D, ZA<:ZForm{D}, ZB<:ZForm{D}} = a*b
∧(s::R, a::ZA) where {R<:Real, D, ZA<:ZForm{D}} = s*a
∧(a::ZA, s::R) where {D, ZA<:ZForm{D}, R<:Real} = s*a
∧(a::R, b::S) where {R<:Real, S<:Real} = a*b
⋆(a::B) where {D, ZA<:ZForm{D}, B<:Blade{ZA}} = a.x*(⋆(one(untype(a))))
Base.:*(a::B, s::R) where {D, ZA<:ZForm{D}, B<:Blade{ZA}, R<:Real} = B(a.x*s)
Base.:*(s::R, a::B) where {D, ZA<:ZForm{D}, R<:Real, B<:Blade{ZA}} = B(s*a.x)
Base.sign(a::K) where{D, K<:ZForm{D}} = 1

Base.zero(::Type{B}) where {T,B<:Blade{T}} = zero(T)

"""
    pseudoscalar(A)

The pseudoscalar associated with algebra of blade A.
The pseudoscalar is the singleton blade that spans the entire space.
The highest grade blade in the algebra.
"""
pseudoscalar(a::K) where {T<:Number,K<:Blade{T}} = pseudoscalar(K)(one(T))

"""
    dual(k)

Map from the subspace k to the complementary dual subspace n-k where n is the highest grade of the algebra.
"""
#dual(E::Type{KT}) where {T<:Real, K, KT<:Blade{T,K}} = pseudoscalar(KT)
dual(R::Type{RT}) where {RT<:Real} = typeof(dual(one(R)))

"""
    ⟂(b)

The orthogonal complement of this blade against the pseudoscalar I.

defined as k*I⁻¹

If the metric is degenerate, the non-metric duality mapping dual(b) should be used instead.
"""
⟂(k::K) where {T<:Number, K<:Blade{T}} = k*inv(pseudoscalar(k))
orthogonal_complement(k) = ⟂(k)

function nonmetric_dual(k::K) where {T<:Number, K<:Blade{T}}
  𝐼 = raise(pseudoscalar(k))
  lower(k*inv(𝐼))
end

Blade{T,0}(s::T) where T = s
Base.:+(b::KV, c::KV) where {T<:Number, K, KV<:Blade{T,K}} = KV(b.x+c.x)
Base.:-(b::KV, c::KV) where {T<:Number, K, KV<:Blade{T,K}} = KV(b.x-c.x)

"The magnitude of a blade.  Oriented."
magnitude(a::B) where {B<:Blade} = a.x
const scalar = magnitude

Base.show(io::IO, k::K) where {K<:Blade} = print(io,k.x,split(string(K), "{")[1])

strnot0(z,s) = isazero(z) ? "" : s

function Base.show(io::IO, z::Z) where {N, Z<:ZForm{N}}
  subs = ["ₓ", "y", "z", "w", "ᵢ", "ᵢ", "ᵢ", "ᵢ", "ᵢ", "ᵢ", "ᵢ", "ᵢ"]

  gstring = mapreduce(i->strnot0(z.gradient[i],", ϕ′"*subs[i]), *, 1:length(z.gradient); init = "")
  hstring = mapreduce(i->strnot0(z.hessian[i],", ϕ′′"*subs[i]), *, 1:length(z.hessian); init = "")

  print("ZForm["*strnot0(z.f,"ϕ")*gstring*hstring*"]")
end

struct BasisVectorInfo
  name::Symbol
  dualname::Symbol
  eᵢ::Vector{Int}
  eemetric::Float64 # maybe should change to Int, so as not to force promotion
  edemetric::Float64 # maybe should change to Int, so as not to force promotion
  n₋::Int64
end

function extract_basis(metric::String)
  swaps(n) = n*(n-1)/2
  countm(s, m::Char) = (length∘findall)(c->c==m, mapreduce( i->metric[i], *, s ))

  subspacei = combinations( 1:length(metric) )
  names = map( c->mapreduce( i->tosub(i), *, c), subspacei )
  dualnames = map( c->mapreduce( i->tosup(i), *, c), subspacei )
  name = map( s->Symbol("e"*reduce(*,s)), names )
  dualname = map( s->Symbol("e"*reduce(*,s)), dualnames )
  eemetric = map( s->((-1)^swaps(length(s)))*
                      ((-1)^countm(s, '-'))*
                      (0^countm(s, '0')), subspacei )
  edemetric = map( s->((-1)^swaps(length(s)))*
                      ((-1)^countm(s, '-'))*
                      (1^countm(s, '0')), subspacei )
  n₋ = map( s->countm(s,'-'), subspacei)
  map(zip(subspacei, name, dualname, eemetric, edemetric, n₋)) do (s, n, d, m, md,nneg)
    BasisVectorInfo(n,d,s,m,md,nneg)
  end
end

"""
    alle(m,s)

An Array with all k-blade types for a given metric string s
"""
alle(gamodule, metrics) = map(bi->getproperty(gamodule,bi.name),extract_basis(metrics))

"""
    alld(m,s)

An Array with all covector types for a given metric string s
"""
alld(gamodule, metrics) = map(bi->getproperty(gamodule,bi.dualname),extract_basis(metrics))

"""
    alle(modulename, n)

Array with all k-blade types where k <= n"
"""
alle(gamodule, n::T) where {T<:Integer} = alle(gamodule, repeat("+",n))

"""
    alld(modulename, n)

Array with all dual covector types where k <= n"
"""
alld(gamodule, n::T) where {T<:Integer} = alld(gamodule, repeat("+",n))

"convenience function to show graded basis at top-level module scope"
show_basis() = alle(Main,grade(dual(1))) |> show

"implementation specific helper.  since Symbol('Main.foo') != :Main.foo need to do some name mangling"
tosym(ns,vare) = Symbol(replace(string(ns)*"_"*string(vare), "."=>"_"))

function basis_expr(name, super)
  Expr(:struct, false, Expr(:<:, Expr(:curly, name, :T), super),
      Expr(:block,
           Expr(:(::), :x, :T)))
end

function quadmetric_expr(ns, namea, nameb, abmetric)
    Expr(:(=),
         Expr(:call,
              Expr(:., :Blades, Expr(:quote, #QuoteNode
                                     :quadmetric,
                                    )),
              Expr(:(::), Expr(:curly, :Val, Expr(:quote, #QuoteNode
                                                  tosym(ns,namea)))),
              Expr(:(::), Expr(:curly, :Val, Expr(:quote, #QuoteNode
                                                  tosym(ns,nameb))))
             ),
         Expr(:block, abmetric ))
end

function dual_expr( ns, ename, dual_fcn )
  Expr(:(=),
       Expr(:where,
            Expr(:call,
                 Expr(:., :Blades, Expr(:quote, #QuoteNode
                                        :dual
                                       )),
                 Expr(:(::), :k, :T)),
            Expr(:<:, :T, ename)), Expr(:block,
                                         Expr(:call, dual_fcn, :k)
                                        ))
end

"""
    @generate_basis(s, [export_types, generate_reciprocals, zero_mixed_index])

generate types for basis Blades and algebra given a string encoding the desired metric.
example metric strings :
"+++" for 3D space
"-+++" for minkowski spacetime
"0+++" 3D projective space ( homegeneous coords )

Optional boolean parameters:
export_types -       export the basis types into calling module scope

generate_reciprocals - generate reciprocal basis types with raised indices for all grades

zero_mixed_index -   any operator that results in a blade with mixed dual and primal indices will be set to zero.  You want this true when working with differential k-forms where eᵢ*eʲ = δᵢʲ
"""
macro generate_basis(metric, export_types=false, generate_reciprocals=false, zero_mixed_index=false)
  ns = __module__
  degenerate_metric = false
  e = extract_basis(metric)
  dualbasis = Vector{Expr}()
  types = Vector{Symbol}()
  basis = map(e) do eᵢ
    if length(eᵢ.eᵢ) == 1
      if isone(eᵢ.eemetric)
        et = :e₊
      elseif isone(-eᵢ.eemetric)
        et = :e₋
      else
        et = :e₀
        degenerate_metric = true
        if zero_mixed_index
          @warn "degenerate metric '0' and zeroing mixed index basis can cause dual to stay degenerate"
        end
      end
      super = Expr(:curly, et, :T)
    else
      super = Expr(:curly, :Blade, :T, length(eᵢ.eᵢ))
    end
    push!(dualbasis,
          basis_expr(eᵢ.dualname, Expr(:curly, :Blade, :T, length(eᵢ.eᵢ)))) #Expr(:curly, :e⁺, :T)))
    push!(types, eᵢ.dualname)
    ne = length(eᵢ.eᵢ)
    if generate_reciprocals && ne > 1
      # use graycodes to get all combinations of sub/super-script
      gc = Iterators.filter(ishetro, eachrow(graycodes(ne)))
      mixede = map(gc) do (gcᵢ)
        choosescript(j) = gcᵢ[j]==1 ? tosub(eᵢ.eᵢ[j]) : tosup(eᵢ.eᵢ[j]);
        mixedname = mapreduce(j->choosescript(j), *, 1:ne; init = "e")
        push!(types, Symbol(mixedname))
        basis_expr(Symbol(mixedname), Expr(:curly, :eₓʸ, :T))
      end
      dualbasis = vcat(dualbasis,mixede)
    end
    push!(types, eᵢ.name)
    basis_expr(eᵢ.name, super)
  end

  exporteᵢ = Vector{Expr}()
  if export_types
    exporteᵢ = map(types) do typeᵢ
      Expr(:export, typeᵢ)
    end
  end

  ek1 = filter(partial((==), 1)∘length∘partial(swap(getproperty),:eᵢ),e)
  quadmetric_fcn = map(ek1) do eᵢ
    quadmetric_expr(ns, eᵢ.name, eᵢ.name, eᵢ.eemetric)
  end

  if generate_reciprocals || degenerate_metric
    quadmetric_fcn =
      vcat(quadmetric_fcn, map(ek1) do eᵢ
            quadmetric_expr(ns, eᵢ.name, eᵢ.dualname, eᵢ.edemetric)
           end)
    quadmetric_fcn =
      vcat(quadmetric_fcn, map(ek1) do eᵢ
            quadmetric_expr(ns, eᵢ.dualname, eᵢ.name, eᵢ.edemetric)
           end)
    quadmetric_fcn =
      vcat(quadmetric_fcn, map(ek1) do eᵢ
            quadmetric_expr(ns, eᵢ.dualname, eᵢ.dualname, eᵢ.edemetric)
           end)
  end

  if zero_mixed_index
    push!(quadmetric_fcn, quadmetric_expr( ns, :eₓʸ, :eₓʸ, 0.0 ))
  else
    push!(quadmetric_fcn, quadmetric_expr( ns, :eₓʸ, :eₓʸ, 1.0 ))
  end

  tehdual = :(⟂)
  if degenerate_metric
    tehdual = :nonmetric_dual
  end
  dual_fcn = map(types) do bᵢ
    dual_expr(ns, bᵢ, tehdual)
  end

  n = length(e[end].eᵢ)
  I₀ = e[end].name
  I⁰ = e[end].dualname
  invsign = (-1)^(n*(n-1)/2 + e[end].n₋) < 0 ? -1 : 1

  # dual of the pseudoscalar
  #==
  push!(dual_fcn,
        Expr(:(=),
             Expr(:call, Expr(:., :Blades, Expr(:quote, #QuoteNode
                                                      :dual
                                                     )), Expr(:(::), :k, I₀)),
              Expr(:block,
                        Expr(:., :k,
                             Expr(:quote, #QuoteNode
                                  :x)))))
                                  ==#

  # dual of scalar.
  push!(dual_fcn,
        Expr(:(=), Expr(:where,
                        Expr(:call, Expr(:., :Blades, Expr(:quote, #QuoteNode
                                                                 :dual
                                                                )),
                             Expr(:(::), :s, :T)),
                        Expr(:<:, :T, :Real)),
             Expr(:block,
                  Expr(:call, I₀, Expr(:call, :*, invsign, :s)))))

  # need to tie pseudoscalar to each blade type in the module.  dual(1.0) may end up in wrong etype system
  # actually what if you take the dual of an operation that returns a scalar?  could get wrong dual fcn
  # really need to have a type e₀ for scalars.  we depend on type choosing correct module.
  #!me add as a known bug
  pseudoscalar_fcn = map(e) do eᵢ
    eᵢname = eᵢ.name
    [:(Blades.pseudoscalar(::Type{K}) where {K<:$eᵢname} = $I₀),
     :(Blades.pseudoscalar(::Type{K}) where {T,K<:$eᵢname{T}} = $I₀{T})]
  end |> flatten
  # also pseudoscalar for dual space
  dpseudoscalar_fcn = map(e) do eᵢ
    eᵢname = eᵢ.dualname
    [:(Blades.pseudoscalar(::Type{K}) where {K<:$eᵢname} = $I⁰),
     :(Blades.pseudoscalar(::Type{K}) where {T,K<:$eᵢname{T}} = $I⁰{T})]
  end |> flatten
  # and for ZForm
  #!me only works for types we create via macro ( unique within module ).
  # zpseudoscalar_fcn = :(Blades.pseudoscalar(::Type{K}) where {K<:ZForm} = $I₀)

  return Expr(:escape, Expr(:block,
              exporteᵢ...,
              basis...,
              dualbasis...,
              quadmetric_fcn...,
              pseudoscalar_fcn..., dpseudoscalar_fcn..., # zpseudoscalar_fcn,
              dual_fcn...,
             ))
end

#!me convert to generated, see the Blade{T,1} version below
"Extract indices for basis vectors spaning the subspace of the given Blade"
function subspace( b::Type{K} ) where {K<:Blade}
  # we could probably encode all the basis information in the Blade type
  # then we wouldn't need this string based malarky.
  # i.e. Blade{T,Tuple{Ea,Eb,Ec}} where Ea,Eb and Ec are the eᵢ indices
  fromsub = Dict('₀'=>0, '₁'=>1, '₂'=>2, '₃'=>3, '₄'=>4, '₅'=>5, '₆'=>6, '₇'=>7, '₈'=>8, '₉'=>9,
                 '⁰'=>0, '¹'=>1, '²'=>2, '³'=>3, '⁴'=>4, '⁵'=>5, '⁶'=>6, '⁷'=>7, '⁸'=>8, '⁹'=>9)
  s = split(split(string(b), ".")[end],"{")[1]
  map(sᵢ->fromsub[sᵢ], collect(s)[2:end])
end

"""
    isdualsubspace(b)

array where value is 1 if the corresponding subspace basis vector is from dualspace of the given Blade
length of array is k where k is the grade of the Blade
"""
function isdualsubspace( b::Type{K} ) where {K<:Blade}
  isdual = Dict('₀'=>0, '₁'=>0, '₂'=>0, '₃'=>0, '₄'=>0, '₅'=>0, '₆'=>0, '₇'=>0, '₈'=>0, '₉'=>0,
                '⁰'=>1, '¹'=>1, '²'=>1, '³'=>1, '⁴'=>1, '⁵'=>1, '⁶'=>1, '⁷'=>1, '⁸'=>1, '⁹'=>1)

  s = split(split(string(b), ".")[end],"{")[1]
  map(sᵢ->isdual[sᵢ], collect(s)[2:end])
end

subspace( b::K ) where {K<:Blade} = subspace(K)

@generated function subspace( b::K ) where {T,K<:Blade{T,1}}
  si = subspace(b)[1]
  :($si)
end

Base.:*(s::T, a::K) where {T<:Real, KT<:Real, K<:Blade{KT}} = K(s*a.x)
Base.:*(a::K, s::T) where {T<:Real, KT<:Real, K<:Blade{KT}} = K(s*a.x)
Base.:/(a::T, s::Real) where {R,T<:Blade{R}} = R(one(R)/s)*a
Base.:/(s::Real, a::T) where {R,T<:Blade{R}} = s*inv(a)
Base.:/(a::Blade, b::Blade) = a*inv(b)
Base.:-(a::T) where {R,T<:Blade{R}} = T(-a.x)
Base.sign(a::K) where{K<:Blade} = sign(a.x)

function swap_parity(eab)
  m = 1
  # trusty old bubble sort.
  bubn = length(eab)-1
  swapped = true
  while swapped
    swapped = false
    for i in 1:bubn
      if eab[i] > eab[i+1]
        temp = eab[i]
        eab[i] = eab[i+1]
        eab[i+1] = temp
        m = -m
        swapped = true
      end
    end
  end
  (eab,m)
end

@generated function Base.:*(a::T, b::U) where {S<:Number,T<:Blade{S}, R<:Number,U<:Blade{R}}
  ns = parentmodule(a)
  @assert ns == parentmodule(b)

  ea = zip(subspace(a), isdualsubspace(a)) |> collect
  eb = zip(subspace(b), isdualsubspace(b)) |> collect
  toscript(eᵢ) = (second(eᵢ) == 0 ? tosub(first(eᵢ)) : tosup(first(eᵢ)))
  tosyme(eᵢ) = tosym(ns,"e"*toscript(eᵢ))

  # effectively implement the geometric algebra's "multiplication table" at compile time
  allzero = false
  m = one(S)
  nea = length(ea)
  eab = Vector{Tuple{Int,Int}}()
  for i in 1:nea
    eaᵢ = ea[i]
    j = findfirst(x->first(x)==first(eaᵢ),eb)
    if j == nothing
      push!(eab,eaᵢ)
    else
      aswaps = nea-i
      totalswaps = ((j-1)+aswaps)

      m = m*(-one(S))^totalswaps
      m = m*(FunctionWrapper{S,
                             Tuple{Val{tosyme(eaᵢ)},
                                   Val{tosyme(eb[j])}}}(quadmetric))(Val(tosyme(eaᵢ)),
                                                                     Val(tosyme(eb[j])))
      splice!(eb, j)
    end
  end

  eab = vcat(eab,eb)
  eab,swp = swap_parity(eab)
  m = m*swp

  hasmixedindex(es) = begin l = mapreduce(second,+,es); l < length(es) && l > 0 end
  if length(eab) > 1 && hasmixedindex(eab)
    m = m*(FunctionWrapper{S,
                           Tuple{Val{tosym(ns, :eₓʸ)},
                                 Val{tosym(ns, :eₓʸ)}}}(quadmetric))(Val(tosym(ns, :eₓʸ)),
                                                                     Val(tosym(ns, :eₓʸ)))
  end

  if isempty(eab)
    :(a.x*b.x*$m)
  else
    if iszero(m)
      m
    else
      TU = Symbol("e"*mapreduce(toscript, *, eab))
      :($ns.$TU(a.x*b.x*$m))
    end
  end
end

"""convenience operator to build higher grade Blade Types.  indices must be in acending order"""
function Base.:*(::Type{A}, ::Type{B}) where {A<:Blade, B<:Blade}
  subs = vcat(subspace(A), subspace(B))
  ascend = mapreduce(((a,b),)->a<b, (a,b)->a&&b, zip(subs[1:end-1],subs[2:end]); init = true)
  if ascend
    untype(1A*1B)
  else
    Nothing
  end
end

"""convenience operator to build higher grade Blade Types.  indices must be in acending order"""
∧(::Type{A}, ::Type{B}) where {A<:Blade, B<:Blade} = A*B

"""
    reverse(b)

Blade constructed with reversed order of basis vectors.
i.e. reverse(e₁∧e₂∧e₃) = e₃∧e₂∧e₁
"""
@generated function Base.reverse( a::T ) where {S<:Number,K,T<:Blade{S,K}}
  s = S((-1)^(K*(K-1)/2))
  # convert each 1-blade basis to it's dual (n-1)-blade basis
  # multiply them together to get the dual for this k-blade
  :($s*a)
end

"""
    lower(k)

direct map taking all upper indices ( dual basis ) to lower indices ( standard basis )
e.g. lower(1e¹²₃) == 1e₁₂₃
"""
@generated function lower( k::Type{K} ) where {K<:Blade}
  ns = parentmodule(K)
  TU = Symbol("e"*mapreduce(tosub, *, subspace(K)))
  :($ns.$TU)
end

lower( k::K ) where {K<:Blade} = lower(K)(k.x)

"""
    raise(k)

direct map taking all lower indices ( standard basis ) to upper indices ( dual basis )
e.g. raise(1e¹²₃) == 1e¹²³
assumes you have generated a dual basis, i.e. @generate_basis("...", _, true)
"""
@generated function raise( k::Type{K} ) where {K<:Blade}
  ns = parentmodule(K)
  TU = Symbol("e"*mapreduce(tosup, *, subspace(K)))
  :($ns.$TU)
end

raise( k::K ) where {K<:Blade} = raise(K)(k.x)

"""
    grade(b)

The dimension of the subspace spanned by the blade
"""
grade(s::T) where {T<:Number} = 0
grade(a::T) where {S<:Number, K, T<:Blade{S,K}} = K
grade(::Type{B}) where {B<:Blade} = grade(B{Int64})
grade(::Type{B}) where {S, K, B<:Blade{S,K}} = K

Base.:*(s::T, e::Type{K}) where {T<:Real, K<:Blade} = K(s)

"""
    inv(k)

Inverse of a blade when it exists.  Left inverse.  k*inv(k) = 1
"""
Base.inv(k::K) where {K<:Blade} = reverse(k)/(k*reverse(k))

Base.:~(k::K) where {K<:Blade} = reverse(k)
Base.:!(k::K) where {K<:Blade} = dual(k)

"""
    ∧(a,b)

Wedge product between two blades
"""
@generated function ∧(a::A, b::B) where {TA,TB,A<:Blade{TA}, B<:Blade{TB}}
  sa = subspace(a)
  sb = subspace(b)
  if isempty(sa ∩ sb)
    :(a*b)
  else
    zero(promote_type(TA,TB))
  end
end

∧(s::T,b::B) where {T<:Real, B<:Blade} = s*b
∧(a::B,s::T) where {B<:Blade, T<:Real} = a*s

∧(a::T,::Type{B}) where {T<:Number, B<:Blade} = B(a)

import LinearAlgebra: ⋅
@generated function ⋅(a::A, b::B) where {TA,TB,A<:Blade{TA}, B<:Blade{TB}}
  sa = subspace(a)
  sb = subspace(b)
  if (sa ∩ sb) == sa
    :(a*b)
  else
    zero(promote_type(TA,TB))
  end
end

⋅(s::T,b::B) where {T<:Real, B<:Blade} = s*b
⋅(a::B,s::T) where {TB,B<:Blade{TB}, T<:Real} = zero(promote_type(T,TB))

Base.:(==)(a::B, b::B) where {B<:Blade} = a.x == b.x
Base.:(==)(a::A, b::B) where {A<:Blade, B<:Blade} = iszero(a.x) && iszero(b.x)
Base.:(==)(a::A, s::T) where {A<:Blade, T<:Real} = iszero(a.x) && iszero(s)
Base.:(==)(s::T, a::A) where {A<:Blade, T<:Real} = a==s
Base.abs(b::B) where {B<:Blade} = B(abs(b.x))

"""
    conj(b)

The complex conjugate of the blade b.
"""
Base.conj(b::B) where {B<:Blade} = reverse(b)

"""
    kform(α, u)

Apply 1-form α to u.
"""
kform( α::K, u::V ) where { TF,T,N, K<:Blade{TF,N}, V<:Blade{T,N} } = α⋅reverse(u)

import LinearAlgebra: det

det( a::B ) where { B<:Blade } = (a⋅a)/magnitude(a)

"""
  det( a, b )

The metric scalar product between two Blades of same grade.
As defined in Geometric Algebra for Computer Science.
"""
@generated function det( a::A, b::B ) where { T, N, A<:Blade{T,N}, B<:Blade{T,N} }
  ea = basis_kblades(a,1)
  eb = basis_kblades(b,1)

  sa = subspace(a)
  sb = subspace(b)
  #==
  aᵢ = one(T).*ea[sa]
  bⱼ = one(T).*eb[sb]
  gramdet = det(aᵢ .⋅ reverse(bⱼ)')
  scalar(a)*scalar(b)*gramdet
  ==#
  aᵢ = ea[sa]
  bⱼ = eb[sb]
  :(a.x*b.x*det(one($T).*$aᵢ .⋅ reverse(one($T).*$bⱼ)'))
end

det( a::A, b::B ) where {T, N, M, A<:Blade{T,N}, B<:Blade{T,M}} = zero(T)

"""
    ⋆(k, 𝑖)

Hodge star operator mapping k to it's Hodge dual relative to a given psuedovector 𝑖.
Defined by k∧⋆(k) == (k⋅k)*𝑖 where k is generated from orthonormal 1-vectors.
𝑖 is the unit psuedoscalar for a subspace containing k.
"""
@generated function ⋆(b::B, i::BI) where {T, B<:Blade{T}, BI<:Blade}
  sb = subspace(b)
  n = reduce(max,sb)
  # shuffle all missing generating 1-blades from hodge dual into k to fill holes
  # track sign flips
  holes = setdiff(1:n,sb)
  # 'filling holes' is so we can make k∧⋆(k) == 𝐼. i.e. ⋆k = holes∧𝐼[n:end] * signflip
  subpseudo,sflip = swap_parity(vcat(sb,holes))
  # then adjust with sign of k⋅k
  sflip = sflip*swap_parity(vcat(sb,sb))[2]
  :(T($sflip)*sign(b)*abs(b*i))
end

⋆(b::B, i::Type{BI}) where {T, B<:Blade{T}, BI<:Blade} = ⋆(b, i(one(T)))

"""
    ⋆(k)

Hodge star operator mapping k to it's Hodge dual.
Defined by k∧⋆(k) == 𝐼 where k is generated from orthonormal 1-vectors and
𝐼 is the psuedoscalar for the generating vector space.
"""
⋆(b::B) where {T,B<:Blade{T}} = ⋆(b, pseudoscalar(b))

⋆(s::Number) = untype(⋆dual(1)*dual(1))(s)
⋆(s::Number, i::BI) where {BI<:Blade} = i
⋆(s::Number, i::Type{BI}) where {BI<:Blade} = ⋆(s, i(one(s)))

import LinearAlgebra: norm,norm_sqr,normalize

norm_sqr(b::K) where {K<:Blade} = magnitude(b)*magnitude(b)
norm(b::K) where {K<:Blade} = abs(magnitude(b))
normalize(b::K) where {T, K<:Blade{T}} = sign(magnitude(b))*one(T)∧untype(b)

"""
    basis_kblades(k, n)

a list of all basis blades of degree n for the algebra k belongs to.

i.e. basis_kblades(e₁, 2) will depend the dimension of the algebra e₁ was generated in
"""
@generated function basis_kblades(::Type{K}, d::Val{N}) where {N,K}
  ns = parentmodule(K)
  n = grade((FunctionWrapper{Any,
                             Tuple{Type{K},}}(pseudoscalar))(K))
  kone = alle(ns,n)[1:n]
  b = filter(e->grade(e)==N, alle(ns, n))
  :($b)
end

basis_kblades(k::K, d::Val{N}) where {N,K} = basis_kblades(K, d)
basis_kblades(u,k::Int) = basis_kblades(u,Val(k))

"""
    basis_1blades(k)

a list of all basis 1-blades for the algebra k belongs to.

i.e. basis_1blades(e₁) will depend the dimension of the algebra e₁ was generated in
"""
basis_1blades( k::B ) where {B<:Blade} = basis_kblades(k,1)
basis_1blades( k::Type{B} ) where {B<:Blade} = basis_kblades(k,1)

"""
    factor(b)

all the smallest factors, that when wedged together create b.
"""
@generated function factor(b::B) where {B<:Blade}
  bs = subspace(b)

  :((b.x,basis_1blades(b)[$bs]...))
end

"Convert blade value to their type with parameter removed"
@generated function untype(b::B) where {T,B<:Blade{T}}
  ns = parentmodule(b)
  n = grade((FunctionWrapper{Any,
                             Tuple{Type{b},}}(pseudoscalar))(b))
  es = alle(ns,n)
  es = vcat(es,raise.(es))
  une = es[findfirst(eᵢ->eᵢ{T} == b, es)]
  :($une)
end

end # module


text:
macro implement_iteration(t)
quote
function Base.iterate(a::$t)
    cc = eachchunk(a)
    ii = iterate(cc)
    ii === nothing && return nothing
    cifirst, cistate = ii
    datacur = (cistate,a[toRanges(cifirst)...])
    innerinds = eachindex(datacur[2])
    ii = iterate(innerinds)
    ii === nothing && return nothing
    innerfirst, innerstate = ii
    return datacur[2][innerfirst],(datacur,cc,cistate,innerinds,innerstate)
end
function Base.iterate(a::$t,i)
    datacur,cc,cistate,innerinds,innerstate = i
    ii = iterate(innerinds, innerstate)
    if ii===nothing
        cii = iterate(cc,cistate)
        cii === nothing && return nothing
        cinow, cistate = cii
        datacur = (cistate,a[toRanges(cinow)...])
        innerinds = eachindex(datacur[2])
        ii = iterate(innerinds)
        ii === nothing && return nothing
    end
    innerfirst, innerstate = ii
    datacur[1] == cistate || error("Something has messed up this iterator")

    return datacur[2][innerfirst],(datacur,cc,cistate,innerinds,innerstate)
end
end
end


text:
@inline metric(κ::Kernel) = κ.metric

## Allows to iterate over kernels
Base.length(::Kernel) = 1
Base.iterate(k::Kernel) = (k,nothing)
Base.iterate(k::Kernel, ::Any) = nothing

### Syntactic sugar for creating matrices and using kernel functions
for k in [:ExponentialKernel,:SqExponentialKernel,:GammaExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel,:LinearKernel,:PolynomialKernel,:ExponentiatedKernel,:ZeroKernel,:WhiteKernel,:ConstantKernel,:RationalQuadraticKernel,:GammaRationalQuadraticKernel]
    @eval begin
        @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
        @inline (κ::$k)(x::AbstractVector{<:Real},y::AbstractVector{<:Real}) = kappa(κ,evaluate(κ.metric,transform(κ,x),transform(κ,y)))
        @inline (κ::$k)(X::AbstractMatrix{T},Y::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,Y,obsdim=obsdim)
        @inline (κ::$k)(X::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,obsdim=obsdim)
    end
end

### Transform generics
@inline transform(κ::Kernel) = κ.transform
@inline transform(κ::Kernel,x::AbstractVecOrMat) = transform(κ.transform,x)
@inline transform(κ::Kernel,x::AbstractVecOrMat,obsdim::Int) = transform(κ.transform,x,obsdim)

## Constructors for kernels without parameters
for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52Kernel,:ExponentiatedKernel]
    @eval begin
        $kernel(ρ::T=1.0) where {T<:Real} =   $kernel{T,ScaleTransform{Base.RefValue{T}}}(ScaleTransform(ρ))
        $kernel(ρ::A) where {A<:AbstractVector{<:Real}} = $kernel{eltype(A),ScaleTransform{A}}(ScaleTransform(ρ))
        $kernel(t::Tr) where {Tr<:Transform} = $kernel{eltype(t),Tr}(t)
    end
end

function set_params!(k::Kernel,x)
    @error "Setting parameters to this kernel is either not possible or has not been implemented"
end

params(k::Kernel) = (params(k.transform),)


text:
# Original file is a part of CoolProp project http://www.coolprop.org/coolprop/wrappers/Julia/index.html#julia
# 73+13+14=100 == 106-1(common)-4(solver)-1(print)
# 12+13+39+8+1=73 == 74-1(common) => doubles OK
# 12+1=13 == 17-4(solver) => SteamStates OK
# 6+7+1=14 == 15-1(print) => Int OK
module FreeSteam

immutable SteamState
    region::Char
    x1::Float64
    x2::Float64
end

immutable XYZ
  x::Char
  y::Char
  z::Char
end

#all 12 SteamState functions (double, double) exported (part 1)
  const set_xy = [:freesteam_set_Ts, :freesteam_set_Tx, :freesteam_set_pT, :freesteam_set_ph,
                    :freesteam_set_ps, :freesteam_set_pu, :freesteam_set_pv, :freesteam_set_uv]
@eval export($(set_xy...))

#all 12 SteamState functions (double, double) exported (part 2)
const reg_set_xy = [:freesteam_region1_set_pT, :freesteam_region2_set_pT, :freesteam_region3_set_rhoT, :freesteam_region4_set_Tx]
@eval export($(reg_set_xy...))
for f in [set_xy; reg_set_xy]
  @eval function $f(x::Number, y::Number)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), SteamState, (Float64, Float64), x, y)
  end
end

#all 13 double functions (SteamState) exported
const _S = [:freesteam_k,  :freesteam_mu, :freesteam_x,   :freesteam_w,
             :freesteam_cv, :freesteam_cp, :freesteam_s, :freesteam_h,
             :freesteam_u,  :freesteam_v,  :freesteam_rho, :freesteam_p,
             :freesteam_T]
@eval export($(_S...))
for f in _S
  @eval function $f(S::SteamState)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Float64, (SteamState, ), S)
  end
end

#all 39 double functions (double,double) exported
const z_for_xy = [:freesteam_region4_v_Tx,    :freesteam_region4_u_Tx,    :freesteam_region4_h_Tx,   :freesteam_region4_s_Tx,   :freesteam_region4_cp_Tx,
                  :freesteam_region4_cv_Tx,   :freesteam_region3_p_rhoT,  :freesteam_region3_u_rhoT, :freesteam_region3_s_rhoT, :freesteam_region3_h_rhoT,
                  :freesteam_region3_cp_rhoT, :freesteam_region3_cv_rhoT, :freesteam_region3_w_rhoT, :freesteam_region3_T_ph,   :freesteam_region3_v_ph,
                  :freesteam_region3_T_ps,    :freesteam_region3_v_ps,    :freesteam_region2_v_pT,   :freesteam_region2_u_pT,   :freesteam_region2_s_pT,
                  :freesteam_region2_h_pT,    :freesteam_region2_cp_pT,   :freesteam_region2_cv_pT,  :freesteam_region2_w_pT,   :freesteam_region2_a_pT,
                  :freesteam_region2_g_pT,    :freesteam_region2_T_ph,    :freesteam_region1_u_pT,   :freesteam_region1_v_pT,   :freesteam_region1_s_pT,
                  :freesteam_region1_h_pT,    :freesteam_region1_cp_pT,   :freesteam_region1_cv_pT,  :freesteam_region1_w_pT,   :freesteam_region1_a_pT,
                  :freesteam_region1_T_ph,    :freesteam_k_rhoT,          :freesteam_mu_rhoT,        :freesteam_region1_g_pT]
@eval export($(z_for_xy...))
for f in z_for_xy
  @eval function $f(x::Float64, y::Float64)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Float64, (Float64, Float64), x, y)
  end
end

#all 8 double functions (char,SteamState) exported
const d_for_cS = [:freesteam_region3_dAdvT, :freesteam_region3_dAdTv, :freesteam_region1_dAdTp, :freesteam_region1_dAdpT,
                   :freesteam_region2_dAdTp, :freesteam_region2_dAdpT, :freesteam_region4_dAdTx, :freesteam_region4_dAdxT]
@eval export($(d_for_cS...))
for f in d_for_cS
  @eval function $f(i::Char, S::SteamState)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Float64, (Char, SteamState), c, S)
  end
end

#1 double function (SteamState,xyz) exported
export freesteam_deriv
function freesteam_deriv(S::SteamState, xyz::XYZ)
  return ccall( (:freesteam_deriv, FreeSteamLib), Float64, (SteamState, XYZ), S, xyz)
end

#all 12 double functions (double) exported
const _x = [:freesteam_region4_Tsat_p, :freesteam_region4_rhof_T, :freesteam_region4_psat_T, :freesteam_region4_rhog_T, :freesteam_region4_dpsatdT_T,
             :freesteam_b23_p_T,        :freesteam_b23_T_p,        :freesteam_region3_psat_h, :freesteam_region3_psat_s, :freesteam_drhofdT_T,
             :freesteam_drhogdT_T,      :freesteam_surftens_T]
@eval export($(_x...))
for f in _x
  @eval function $f(x::Float64)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Float64, (Float64,), x)
  end
end

#1 SteamState function (double) exported
export freesteam_bound_pmax_T
function freesteam_bound_pmax_T(x::Float64)
  return ccall( (:freesteam_bound_pmax_T, FreeSteamLib), SteamState, (Float64, ), x)
end

#all 6 int functions (double,double,int) exported
const i_for_xyi = [:freesteam_bounds_ph, :freesteam_bounds_ps, :freesteam_bounds_pv, :freesteam_bounds_Ts, :freesteam_bounds_Tx, :freesteam_bounds_uv]
@eval export($(i_for_xyi...))
for f in i_for_xyi
  @eval function $f(x::Float64, y::Float64, i::Cint)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Cint, (Float64, Float64, Cint), x, y, i)
  end
end

#all 7 int functions (double,double) exported
const i_for_xy = [:freesteam_region_ph, :freesteam_region_ps, :freesteam_region_pu, :freesteam_region_pv, :freesteam_region_Ts, :freesteam_region_Tx, :freesteam_region_uv]
@eval export($(i_for_xy...))
for f in i_for_xy
  @eval function $f(x::Float64, y::Float64)
    return ccall( ($(QuoteNode(f)), FreeSteamLib), Cint, (Float64, Float64), x, y)
  end
end

#1 int function (SteamState) exported
export freesteam_region
function freesteam_region(S::SteamState)
  return ccall( (:freesteam_region, FreeSteamLib), Cint, (SteamState, ), S)
end

@windows_only const FreeSteamLib = abspath(joinpath(@__FILE__,"..","..","lib","freesteam.dll"));
@linux_only const FreeSteamLib = abspath(joinpath(@__FILE__,"..","..","lib","libfreesteam.so"));

end #module