text:
using DataStructures
using Compat
import Base.Ordering
import Base.Forward
import Base.Reverse
import DataStructures.eq
import Base.lt
import Base.ForwardOrdering

## Function fulldump dumps the entire tree; helpful for debugging.

function fulldump(t::DataStructures.BalancedTree23)
    thislevstack = Int[]
    rl = t.rootloc
    dpth = t.depth
    push!(thislevstack, rl)
    println("  rootloc = $rl depth = $dpth")
    nextlevstack = Int[]
    levcount = 0
    for mydepth = 1 : dpth
        isleaf = (dpth == mydepth)
        sz = size(thislevstack,1)
        if sz == 0
            break
        end
        levcount += 1
        println("-------------\n------LEVEL $levcount")
        for i = 1 : sz
            ii = thislevstack[i]
            p = t.tree[ii].parent
            sk1 = t.tree[ii].splitkey1
            sk2 = t.tree[ii].splitkey2
            c1 = t.tree[ii].child1
            c2 = t.tree[ii].child2
            c3 = t.tree[ii].child3
            dt = (mydepth == dpth)? "child(data)" : "child(tree)"
            if c3 == 0
                println("ii = $ii splitkey1 = /$sk1/ $dt.1 = $c1 $dt.2 = $c2 parent = $p")
                if !isleaf
                    push!(nextlevstack,c1)
                    push!(nextlevstack,c2)
                end
            else
                println("ii = $ii splitkey1 = /$sk1/ splitkey2 = /$sk2/ $dt.1 = $c1 $dt.2 = $c2 $dt.3 =$c3 parent = $p")
                if !isleaf
                    push!(nextlevstack,c1)
                    push!(nextlevstack,c2)
                    push!(nextlevstack,c3)
                end
            end
        end
        thislevstack = nextlevstack
        nextlevstack = Int[]
    end
    println("----- FREE TREE CELLS -----")
    for i = 1 : size(t.freetreeinds,1)
        println(t.freetreeinds[i])
    end

    println("----- DATA ARRAY ---")
    for j = 1 : size(t.data,1)
        k = t.data[j].k
        d = t.data[j].d
        p = t.data[j].parent
        println("j = $j k = /$k/ d = /$d/ parent = /$p/")
    end
    println("----- FREE DATA CELLS -----")
    for i = 1 : size(t.freedatainds,1)
        println(t.freedatainds[i])
    end
end

## Function checkcorrectness checks a balanced tree for correctness.

function checkcorrectness{K,D,Ord <: Ordering}(t::DataStructures.BalancedTree23{K,D,Ord})
    dsz = size(t.data, 1)
    tsz = size(t.tree, 1)
    r = t.rootloc
    bfstreenodes = Array(Int, 0)
    tdpth = t.depth
    intree = IntSet()
    levstart = Array(Int, tdpth)
    push!(bfstreenodes, r)
    levstart[1] = 1
    push!(intree, r)
    for curdepth = 2 : tdpth
        levstart[curdepth] = size(bfstreenodes, 1) + 1
        for l = levstart[curdepth - 1] : levstart[curdepth] - 1
            anc = bfstreenodes[l]
            c1 = t.tree[anc].child1
            if in(c1, intree)
                println("anc = $anc c1 = $c1")
                error("Tree contains loops 1")
            end
            push!(bfstreenodes, c1)
            push!(intree, c1)
            c2 = t.tree[anc].child2
            if in(c2, intree)
                error("Tree contains loops 2")
            end
            push!(bfstreenodes, c2)
            push!(intree, c2)
            c3 = t.tree[anc].child3
            if c3 > 0
                if in(c3, intree)
                    error("Tree contains loops 3")
                end
                push!(bfstreenodes, c3)
                push!(intree, c3)
            end
        end
    end
    bfstreesize = size(bfstreenodes, 1)
    dataused = IntSet()
    minkeys = Array(K, bfstreesize)
    maxkeys = Array(K, bfstreesize)
    for s = levstart[tdpth] : bfstreesize
        anc = bfstreenodes[s]
        c1 = t.tree[anc].child1
        if s == levstart[tdpth]
            if c1 != 1
                error("Leftmost data descendant should be node 1")
            end
        else
            minkeys[s] = t.data[c1].k
        end
        c2 = t.tree[anc].child2
        c3 = t.tree[anc].child3
        lastchild = c3 > 0? c3 : c2
        if s == bfstreesize
            if lastchild != 2
                error("Rightmost data descendant should be node 2")
            end
        else
            maxkeys[s] = t.data[lastchild].k
        end
        if s > levstart[tdpth] && !lt(t.ord, maxkeys[s - 1], minkeys[s])
            error("Data nodes out of order")
        end
        if s < bfstreesize || c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey1, t.data[c2].k)
                error("Splitkey1 of leaf should match key of 2nd child")
            end
        end
        if s < bfstreesize && c3 > 0
            if !eq(t.ord, t.tree[anc].splitkey2, t.data[c3].k)
                error("Splitkey2 of leaf should match key of 1st child")
            end
        end
        if t.data[c1].parent != anc || t.data[c2].parent != anc ||
            (c3 > 0 && t.data[c3].parent != anc)
            println("c1 = $c1 c2 = $c2 c3 = $c3 anc = $anc")
            println("t.data[c1].parent = $(t.data[c1].parent) t.data[c2].parent = $(t.data[c2].parent)")
            if c3 > 0
                println("t.data[c3].parent = $(t.data[c3].parent)")
            end
            error("Incorrect parent node for data child")
        end
        push!(dataused, c1)
        push!(dataused, c2)
        if c3 > 0
            push!(dataused, c3)
        end
    end
    for curdepth = tdpth - 1 : -1 : 1
        cp = levstart[curdepth + 1]
        for s = levstart[curdepth] : levstart[curdepth + 1] - 1
            anc = bfstreenodes[s]
            c1 = t.tree[anc].child1
            @assert(c1 == bfstreenodes[cp])
            if s > levstart[curdepth]
                mk1 = minkeys[cp]
            end
            cp += 1
            if t.tree[c1].parent != anc
                error("Parent/child1 links do not match")
            end
            c2 = t.tree[anc].child2
            @assert(c2 == bfstreenodes[cp])
            mk2 = minkeys[cp]
            cp += 1
            if t.tree[c2].parent != anc
                error("Parent/child2 links do not match")
            end
            c3 = t.tree[anc].child3
            @assert(s == levstart[curdepth] || lt(t.ord,mk1,mk2))
            if c3 > 0
                if t.tree[c3].parent != anc
                    error("Parent/child3 links do not match")
                end
                mk3 = minkeys[cp]
                cp += 1
                @assert(lt(t.ord,mk2, mk3))
            end
            if s > levstart[curdepth]
                minkeys[s] = mk1
            end
            if !eq(t.ord, t.tree[anc].splitkey1, mk2)
                error("Minkey2 not equal to minimum key among descendants of child2")
            end
            if c3 > 0 && !eq(t.ord, t.tree[anc].splitkey2, mk3)
                error("Minkey3 not equal to minimum key among descendants of child3")
            end
        end
    end
    freedata = IntSet()
    for i = 1 : size(t.freedatainds,1)
        fdi = t.freedatainds[i]
        if in(fdi, freedata)
            error("t.freedatainds has repeated element $i")
        end
        if fdi < 1 || fdi > dsz
            error("t.freedatainds entry out of range")
        end
        push!(freedata, fdi)
    end
    if in(:useddatacells, names(DataStructures.BalancedTree23))
        if last(t.useddatacells) > dsz
            error("t.useddatacells has indices larger than t.data size")
        end
        for i = 1 : dsz
            if (in(i, dataused) && !in(i, t.useddatacells)) ||
                (!in(i,dataused) && in(i, t.useddatacells))
                error("Mismatch between actual data cells used and useddatacells array")
            end
            if (in(i, freedata) && in(i, dataused)) ||
                (!in(i,freedata) && !in(i, dataused))
                error("Mismatch between t.freedatainds and t.useddatacells")
            end
        end
    end
    freetree = IntSet()
    for i = 1 : size(t.freetreeinds,1)
        tfi = t.freetreeinds[i]
        if in(tfi, freetree)
            error("Free tree index repeated twice")
        end
        if tfi < 1 || tfi > tsz
            error("Free tree index out of range")
        end
        push!(freetree, tfi)
    end
    for i = 1 : tsz
        if (!in(i, intree) && !in(i, freetree)) ||
            (in(i, intree) && in(i, freetree))
            error("Mismatch between t.freetreeinds and actual cells used")
        end
    end
end

function test1()
    # a few basic tests to start
    m1 = SortedDict((@compat Dict{ASCIIString,ASCIIString}()), Forward)
    kdarray = ["hello", "jello", "alpha", "beta", "fortune", "random",
               "july", "wednesday"]
    checkcorrectness(m1.bt)
    for i = 1 : div(size(kdarray,1), 2)
        k = kdarray[i*2-1]
        d = kdarray[i*2]
        #println("- inserting: k = $k d = $d")
        m1[k] = d
        checkcorrectness(m1.bt)
        # fulldump(m1.bt)
    end
    i1 = startof(m1)
    count = 0
    while i1 != pastendtoken(m1)
        count += 1
        k,d = deref(i1)
        #println("+ reading: k = $k, d = $d")
        i1 = advance(i1)
        checkcorrectness(m1.bt)
    end
    @test count == 4
end

function test2()
    # test all the methods here except loops
    m0 = SortedDict(@compat Dict{Int, Float64}())
    m1 = SortedDict(@compat Dict(8=>32.0, 12=>33.1, 6=>18.2))
    expected = ([6,8,12], [18.2, 32.0, 33.1])
    checkcorrectness(m1.bt)
    ii = startof(m1)
    m2 = packdeepcopy(m1)
    m3 = packcopy(m1)
    p = first(m1)
    @test p[1] == 6 && p[2] == 18.2
    @test in((8,32.0),m3)
    @test !in((8,32.1),m3)
    for j = 1 : 3
        @test ii != pastendtoken(m1)
        pr = deref(ii)
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        checkcorrectness(m1.bt)
        oldii = ii
        ii = advance(ii)
        delete!(oldii)
    end
    checkcorrectness(m1.bt)
    checkcorrectness(m2.bt)
    @test length(m2) == 3
    ii = startof(m2)
    for j = 1 : 3
        pr = deref(ii)
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        ii = advance(ii)
    end

    checkcorrectness(m3.bt)
    @test length(m3) == 3
    ii = startof(m3)
    for j = 1 : 3
        pr = deref(ii)
        @test pr[1] == expected[1][j] && pr[2] == expected[2][j]
        ii = advance(ii)
    end

    @test isempty(m1)
    @test length(m1) == 0
    N = 5000
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        if i % 50 == 0
            checkcorrectness(m1.bt)
        end
    end
    @test !isempty(m1)
    assert(length(m1) == N - 1)
    for i = 2 : N
        d = pop!(m1, i)
        @test d == convert(Float64,i)^2
        checkcorrectness(m1.bt)
    end
    @test isempty(m1)
    @test length(m1) == 0
    for i = N : -1 : 2
        m1[i] = convert(Float64,i) ^ 2
        checkcorrectness(m1.bt)
    end
    ii = endof(m1)
    for i = 1 : N - 1
        pr = deref(ii)
        @test pr[1] == N + 1 - i && pr[2] == convert(Float64,pr[1]) ^ 2
        ii = regress(ii)
    end
    lastprime = 1
    while true
        ii = searchsortedafter(m1, lastprime)
        if ii == pastendtoken(m1)
            break
        end
        j = deref_key(ii)
        for k = j * 2 : j : N
            p = find(m1, k)
            if p != pastendtoken(m1)
                delete!(p)
                checkcorrectness(m1.bt)
            end
        end
        lastprime = j
    end
    @test ii == pastendtoken(m1)
    @test status(ii) == 3
    h = assemble(m1, semi(ii))
    @test status(h) == 3
    h2 = assemble(m1, SDSemiToken(0))
    @test status(h2) == 0
    t = 0
    u = 0.0
    for pr = m1
        t += pr[1]
        u += pr[2]
    end
    numprimes = length(m1)
    pn = primes(N)
    @test t == sum(pn)
    @test u == sum(pn.^2)
    ij = endof(m1)
    @test deref_key(ij) == last(pn) && convert(Float64, last(pn)^2) ==  deref_value(ij)
    m1[6] = 49.0
    @test length(m1) == numprimes + 1
    @test m1[6] == 49.0
    b, i6 = insert!(m1, 6, 50.0)
    @test length(m1) == numprimes + 1
    @test !b
    p = deref(i6)
    @test p[1] == 6 && p[2] == 50.0
    b2, i7 = insert!(m1, 8, 51.0)
    @test b2
    st = semi(i6)
    m1[st] = 9.0
    p = deref(i6)
    @test p[1] == 6 && p[2] == 9.0
    @test m1[st] == 9.0
    @test length(m1) == numprimes + 2
    p = deref(i7)
    @test p[1] == 8 && p[2] == 51.0
    delete!(i7)
    z = pop!(m1, 6)
    checkcorrectness(m1.bt)
    @test z == 9.0
    i8 = startof(m1)
    p = deref(i8)
    @test p[1] == 2 && p[2] == 4.0
    @test i8 != beforestarttoken(m1)
    @test status(i8) == 1
    i9 = regress(i8)
    @test i9 == beforestarttoken(m1)
    @test status(i9) == 2
    i10 = find(m1, 17)
    i11 = regress(i10)
    @test deref_key(i11) == 13
    i12 = searchsortedfirst(m1, 47)
    i13 = searchsortedfirst(m1, 48)
    @test deref_key(i12) == 47
    @test deref_key(i13) == 53
    i14 = searchsortedafter(m1, 47)
    i15 = searchsortedafter(m1, 48)
    @test deref_key(i14) == 53
    @test deref_key(i15) == 53
    i16 = searchsortedlast(m1, 47)
    i17 = searchsortedlast(m1, 48)
    @test deref_key(i16) == 47
    @test deref_key(i17) == 47
    ww = primes(N)
    cc = last(m1)
    @test cc[1] == last(ww)
    wwx = first(m1)
    @test wwx[1] == 2
    tpr = eltype(m1)
    @test tpr[1] == Int && tpr[2] == Float64
    co = orderobject(m1)
    @test co == Forward
    @test haskey(m1, 71)
    @test !haskey(m1, 77)
    @test get(m1, 70, -45.2) == -45.2
    @test get(m1, 83, -45.2) == convert(Float64,83)^2
    h = get!(m1, 5, 27.0)
    @test h == 25.0
    h = get!(m1, 6, 27.0)
    @test h == 27.0
    @test m1[6] == 27.0
    @test length(m1) == length(ww) + 1
    pop!(m1, 6)
    @test getkey(m1,7, 9) == 7
    @test getkey(m1,7.0, 9) == 7
    @test getkey(m1,12, 9) == 9
    delete!(m1, 17)
    @test length(m1) == length(ww) - 1
    @test deref_key(advance(find(m1,13))) == 19
    empty!(m1)
    checkcorrectness(m1.bt)
    @test isempty(m1)
    c1 = SortedDict(@compat Dict("Eggplants"=>3,
                        "Figs"=>9,
                        "Apples"=>7))
    c2 = SortedDict(@compat Dict("Eggplants"=>6,
                        "Honeydews"=>19,
                        "Melons"=>11))
    @test !isequal(c1,c2)
    c3 = merge(c1, c2)
    checkcorrectness(c3.bt)
    c4 = SortedDict(@compat Dict("Apples"=>7,
                        "Figs"=>9,
                        "Eggplants"=>6,
                        "Melons"=>11,
                        "Honeydews"=>19))
    @test isequal(c3,c4)
    c5 = SortedDict(@compat Dict("Apples"=>7))
    @test !isequal(c4,c5)
    merge!(c1,c2)
    checkcorrectness(c1.bt)
    @test isequal(c3,c1)
    merge!(c3,c3)
    @test isequal(c3,c1)
    checkcorrectness(c3.bt)
end

function bitreverse(i)
    zeroi = zero(i)
    onei = one(i)
    twoi = onei + onei
    r = zeroi
    for j = 1 : 32
        r *= twoi
        r += (i & onei)
        i = div(i,twoi)
    end
    r
end

## Test the loop constructs.
function test3{T}(z::T)
    zero1 = zero(z)
    one1 = one(z)
    two1 = one1 + one1
    m1 = SortedDict(@compat Dict{T,T}())
    N = 1000
    for l = 1 : N
        lUi = convert(T, l)
        m1[bitreverse(lUi)] = lUi
    end
    count = 0
    for (tok,(k,v)) in tokens(startof(m1) : endof(m1))
        for (tok2,(k2,v2)) in tokens(excludelast(startof(m1), pastendtoken(m1)))
            if isless(tok,tok2)
                @test deref_key(tok) < deref_key(tok2)
            elseif isequal(tok,tok2)
                @test deref_key(tok) == deref_key(tok2)
            else
                @test deref_key(tok) > deref_key(tok2)
            end
            count += 1
        end
    end
    @test count == N^2
    N = 10000
    sk = zero1
    sv = zero1
    skhalf = zero1
    svhalf = zero1
    for l = 1 : N
        lUi = convert(T, l)
        brl = bitreverse(lUi)
        sk += brl
        m1[brl] = lUi
        sv += lUi
        if brl < div(N, 2)
            skhalf += brl
            svhalf += lUi
        end
    end
    count = 0
    sk2 = zero1
    sv2 = zero1
    for (k,v) in m1
        sk2 += k
        sv2 += v
        count += 1
    end
    @test count == N
    @test sk2 == sk
    @test sv == sv2
    sk2 = zero1
    for k in keys(m1)
        sk2 += k
    end
    @test sk2 == sk
    sv2 = zero1
    for v in values(m1)
        sv2 += v
    end
    @test sv == sv2
    count = 0
    for (t,k) in tokens(keys(m1))
        @test deref_key(t) == k
        count += 1
    end
    @test count == N
    count = 0
    for (t,v) in tokens(values(m1))
        @test deref_value(t) == v
        count += 1
    end
    @test count == N

    pos1 = searchsortedfirst(m1, div(N,2))
    sk2 = zero1
    for k in keys(excludelast(startof(m1), pos1))
        sk2 += k
    end
    @test sk2 == skhalf
    sv2 = zero1
    for v in values(excludelast(startof(m1), pos1))
        sv2 += v
    end
    @test sv2 == svhalf
    count = 0
    for (k,v) in excludelast(pastendtoken(m1), pastendtoken(m1))
        count += 1
    end
    @test count == 0
    count = 0
    for (k,v) in startof(m1) :  beforestarttoken(m1)
        count += 1
    end
    @test count == 0
end

# test all the errors
function test4()
    m = SortedDict(@compat Dict("a" => 6, "bb" => 9))
    @test_throws KeyError println(m["b"])
    m2 = SortedDict(@compat Dict{ASCIIString, Int}())
    @test_throws BoundsError println(first(m2))
    @test_throws BoundsError println(last(m2))
    state1 = start(m2)
    @test_throws BoundsError next(m2, state1)
    @test_throws ArgumentError beforestarttoken(m) : pastendtoken(m2)
    @test_throws ArgumentError excludelast(beforestarttoken(m), pastendtoken(m2))
    @test_throws ArgumentError isless(beforestarttoken(m), pastendtoken(m2))
    @test_throws ArgumentError isequal(beforestarttoken(m), pastendtoken(m2))
    i1 = find(m,"a")
    delete!(i1)
    i2 = find(m,"bb")
    @test_throws BoundsError start(i1:i2)
    @test_throws BoundsError start(excludelast(i1,i2))
    @test_throws KeyError delete!(m,"a")
    @test_throws KeyError pop!(m,"a")
    m3 = SortedDict((@compat Dict{ASCIIString, Int}()), Reverse)
    @test_throws ErrorException isequal(m2, m3)
    i1semi = semi(i1)
    @test_throws BoundsError m[i1semi]
    @test_throws BoundsError regress(beforestarttoken(m))
    @test_throws BoundsError advance(pastendtoken(m))
end

function seekfile(fname)
    fullname = joinpath(Pkg.dir("DataStructures"), "test", fname)
end

immutable CaseInsensitive <: Ordering
end

lt(::CaseInsensitive, a, b) = isless(lowercase(a), lowercase(b))
eq(::CaseInsensitive, a, b) = isequal(lowercase(a), lowercase(b))

## Test use of alternative orderings

function test5()
    keylist = ["Apple", "aPPle", "berry", "CHerry", "Dairy", "diary"]
    vallist = [6,9,-4,2,1,8]
    m = SortedDict(@compat Dict{ASCIIString,Int}())
    for j = 1:6
        m[keylist[j]] = vallist[j]
    end
    checkcorrectness(m.bt)
    expectedord1 = [1,4,5,2,3,6]
    count = 0
    for p in m
        count += 1
        @test p[1] == keylist[expectedord1[count]] &&
                p[2] == vallist[expectedord1[count]]
    end
    @test count == 6
    m2 = SortedDict((@compat Dict{ASCIIString, Int}()), Reverse)
    for j = 1 : 6
        m2[keylist[j]] = vallist[j]
    end
    checkcorrectness(m2.bt)
    expectedord2 = [6,3,2,5,4,1]
    count = 0
    for p in m2
        count += 1
        @test p[1] == keylist[expectedord2[count]] &&
                p[2] == vallist[expectedord2[count]]
    end
    @test count == 6
    m3 = SortedDict((@compat Dict{ASCIIString, Int}()), CaseInsensitive())
    for j = 1 : 6
        m3[keylist[j]] = vallist[j]
    end
    checkcorrectness(m3.bt)
    expectedord3 = [2,3,4,5,6]
    count = 0
    for p in m3
        count += 1
        @test p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]]
    end
    @test count == 5
    m4 = SortedDict((@compat Dict{ASCIIString,Int}()), Lt((x,y) -> isless(lowercase(x),lowercase(y))))
    for j = 1 : 6
        m4[keylist[j]] = vallist[j]
    end
    checkcorrectness(m4.bt)
    count = 0
    for p in m4
        count += 1
        @test p[1] == keylist[expectedord3[count]] &&
                p[2] == vallist[expectedord3[count]]
    end
    @test count == 5
end

## test6, test6a and test6b are not run by package
## testing but are included for timing tests.
## They are identical except for their usage of ordering objects.

function test6(numtrial::Int, expectedk::ASCIIString, expectedd::ASCIIString)
    NSTRINGPAIR = 50000
    m1 = SortedDict(@compat Dict{ASCIIString,ASCIIString}())
    strlist = ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!(l)
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance(l)
        end
        k,d = deref(l)
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

function test6a(numtrial::Int, expectedk::ASCIIString, expectedd::ASCIIString)
    NSTRINGPAIR = 50000
    m1 = SortedDict((@compat Dict{ASCIIString, ASCIIString}()), Lt(isless))
    strlist = ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!(l)
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance(l)
        end
        k,d = deref(l)
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

function SDConstruct(a::Associative; lt::Function=isless, by::Function=identity)
    if by == identity
        return SortedDict(a, Lt(lt))
    elseif lt == isless
        return SortedDict(a, By(by))
    else
        error("having both by and lt not both implemented")
    end
end

function test6b(numtrial::Int, expectedk::ASCIIString, expectedd::ASCIIString)
    NSTRINGPAIR = 50000
    m1 = SDConstruct((@compat Dict{ASCIIString,ASCIIString}()), lt=isless)
    strlist = ASCIIString[]
    open(seekfile("wordsScram.txt"), "r") do inio
        for j = 1 : NSTRINGPAIR * 2
            push!(strlist, chomp(readline(inio)))
        end
    end
    for trial = 1 : numtrial
        empty!(m1)
        count = 1
        for j = 1 : NSTRINGPAIR
            k = strlist[count]
            d = strlist[count + 1]
            m1[k] = d
            count += 2
        end
        delct = div(NSTRINGPAIR, 6)
        for j = 1 : delct
            l = find(m1, strlist[j * 2 + 1])
            delete!(l)
        end
        spec = div(NSTRINGPAIR * 3,4)
        l = startof(m1)
        for j = 1 : spec
            l = advance(l)
        end
        k,d = deref(l)
        ekn = expectedk
        edn = expectedd
        @test k == ekn && d == edn
    end
end

test1()
test2()
test3(0x00000000)
test4()
test5()
#test6(2, "soothingly", "compere")


text:
export NewtonResult, newton, newton!, NewtonCache, NewtonReturnCode

module NewtonReturnCode
    @doc """
        NewtonReturnCode.codes

    The possible return codes of Newton's method

    * `NewtonReturnCode.converged`
    * `NewtonReturnCode.terminated`
    * `NewtonReturnCode.terminated_no_approximate`
    * `NewtonReturnCode.maximal_iterations`
    """
    @enum codes begin
        converged
        terminated
        terminated_no_approximate
        maximal_iterations
    end
end

"""
    NewtonResult{T}

Structure holding information about the outcome of the `newton` function. The fields are.
* `return_code::NewtonReturnCode.codes`: The return code of computation. `NewtonReturnCode.converged` means that `accuracy ≤ tol`.
* `accuracy::T`: |xᵢ-xᵢ₋₁| for i = iters and x₀,x₁,…,xᵢ₋₁,xᵢ are the Newton iterates.
* `iters::Int`: The number of iterations used.
* `digits_lost::Float64` Estimate of the (relative) lost digits in the linear algebra.
"""
struct NewtonResult{T}
    return_code::NewtonReturnCode.codes
    accuracy::T
    iters::Int
    ω₀::Float64
    ω::Float64
    norm_Δx₀::T
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", r::NewtonResult) = r
Base.show(io::IO, result::NewtonResult) = print_fieldnames(io, result)

"""
    isconverged(result::Result)

Returns whether the correction was successfull.
"""
isconverged(result::NewtonResult) = result.return_code == NewtonReturnCode.converged

abstract type AbstractNewtonCache end
"""
    NewtonCache(F::AbstractSystem, x)

Cache for the [`newton`](@ref) function.
"""
struct NewtonCache{T, SC<:AbstractSystemCache} <: AbstractNewtonCache
    Jac::Jacobian{T}
    rᵢ::Vector{T}
    Δxᵢ::Vector{T}
    system_cache::SC
end

struct NewtonCacheF64{SC<:AbstractSystemCache, SC2<:AbstractSystemCache, AV<:AbstractVector{Complex{Double64}}} <: AbstractNewtonCache
    Jac::Jacobian{ComplexF64}
    xᵢ_D64::AV
    rᵢ::Vector{ComplexF64}
    rᵢ_D64::Vector{Complex{Double64}}
    Δxᵢ::Vector{ComplexF64}
    system_cache::SC
    system_cache_D64::SC2
end

function newton_cache(F::AbstractSystem, x)
    system_cache = cache(F, x)
    Jac = Jacobian(Random.randn!(jacobian(F, x, system_cache)))
    rᵢ = evaluate(F, x, system_cache)
    Δxᵢ = zeros(eltype(x), length(x))

    NewtonCache(Jac, rᵢ, Δxᵢ, system_cache)
end

function newton_cache(F::AbstractSystem, x::AbstractVector{ComplexF64})
    system_cache = cache(F, x)
    Jac = Jacobian(Random.randn!(complex.(float.(jacobian(F, x, system_cache)))))
    rᵢ = evaluate(F, x, system_cache)
    xD64 = similar(x, Complex{Double64})
    system_cache_D64 = cache(F, xD64)
    rᵢ_D64 = similar(rᵢ, Complex{Double64})
    Δxᵢ = zeros(eltype(x), length(x))

    NewtonCacheF64(Jac, xD64, rᵢ, rᵢ_D64, Δxᵢ, system_cache, system_cache_D64)
end

"""
    newton(F::AbstractSystem, x₀, norm=euclidean_norm, cache=NewtonCache(F, x₀); tol=1e-6, miniters=1, maxiters=3, simplified_last_step=true)

An ordinary Newton's method. If `simplified_last_step` is `true`, then for the last iteration
the previously Jacobian will be used. This uses an LU-factorization for square systems
and a QR-factorization for overdetermined.
"""
function newton(F::AbstractSystem, x₀, norm=euclidean_norm, cache=NewtonCache(F, x₀); kwargs...)
    x = copy(x₀)
    x, newton!(x, F, x₀, norm, cache; kwargs...)
end

"""
    newton!(out, F::AbstractSystem, x₀, norm, cache::AbstractNewtonCache; tol=1e-6, miniters=1, maxiters=3, simplified_last_step=true)

In-place version of [`newton`](@ref). Needs a [`NewtonCache`](@ref) and `norm` as input.
"""
@inline function newton!(out, F::AbstractSystem, x₀, norm, cache::AbstractNewtonCache;
                 tol::Float64=1e-6, miniters::Int=1, maxiters::Int=3,
                 simplified_last_step::Bool=true,
                 update_jacobian_infos::Bool=false,
                 use_qr::Bool=false,
                 ω::Float64=0.0, expected_Δx₀::Float64=0.0,
                 precision::PrecisionOption=PRECISION_ADAPTIVE)
    newton!(out, F, x₀, norm, cache, cache.Jac, tol,
            miniters, maxiters, simplified_last_step,
            update_jacobian_infos, use_qr, ω, expected_Δx₀, precision)
end

function newton!(out, F::AbstractSystem, x₀, norm, cache::AbstractNewtonCache, Jac::Jacobian,
                 tol::Float64, miniters::Int, maxiters::Int, simplified_last_step::Bool,
                 update_jacobian_infos::Bool, use_qr::Bool, ω::Float64,
                 expected_Δx₀::Float64, precision::PrecisionOption)
    rᵢ, Δxᵢ = cache.rᵢ, cache.Δxᵢ
    Jᵢ = Jac.J
    copyto!(out, x₀)
    xᵢ₊₁ = xᵢ = out # just alias to make logic easier
    rᵢ₊₁ = rᵢ

    # initialize variables
    T = real(eltype(xᵢ))
    Θ₀ = Θᵢ₋₁ = norm_Δxᵢ₋₁ = norm_Δxᵢ = norm_Δx₀ = zero(T)
    accuracy = T(Inf)
    ω₀ = ω
    for i ∈ 0:maxiters
        curr_accuracy = i == 0 ? expected_Δx₀ : accuracy
        newton_step!(cache, Jac, F, xᵢ;
                        tol=tol,
                        curr_accuracy=curr_accuracy,
                        ω = ω,
                        simple_step=(i == maxiters && simplified_last_step),
                        use_qr=use_qr,
                        update_infos=(update_jacobian_infos && i == 0),
                        precision=precision)
        norm_Δxᵢ₋₁ = norm_Δxᵢ
        norm_Δxᵢ = Float64(norm(Δxᵢ))
        @inbounds for k in eachindex(xᵢ)
            xᵢ₊₁[k] = xᵢ[k] - Δxᵢ[k]
        end

        if i == 0
            accuracy = norm_Δx₀ = norm_Δxᵢ₋₁ = norm_Δxᵢ
            if norm_Δx₀ ≤ tol && i + 1 ≥ miniters
                return NewtonResult(NewtonReturnCode.converged, norm_Δx₀, i + 1, ω₀, ω, norm_Δx₀)
            end

        else
            Θᵢ₋₁ = norm_Δxᵢ / norm_Δxᵢ₋₁
            ωᵢ₋₁ = 2Θᵢ₋₁ / norm_Δxᵢ₋₁
            if i == 1
                ω = ω₀ = ωᵢ₋₁
            else
                ω = @fastmath max(ω, ωᵢ₋₁)
            end

            if Θᵢ₋₁ > 0.5
                return NewtonResult(NewtonReturnCode.terminated, accuracy, i + 1, ω₀, ω, norm_Δx₀)
            end

            accuracy = norm_Δxᵢ / (1 - 2Θᵢ₋₁^2)
            if accuracy ≤ tol && i + 1 ≥ miniters
                return NewtonResult(NewtonReturnCode.converged, accuracy, i + 1, ω₀, ω, norm_Δx₀)
            end
        end
    end

    return NewtonResult(NewtonReturnCode.maximal_iterations, accuracy, maxiters+1, ω₀, ω, norm_Δx₀)
end

@inline function newton_step!(cache::NewtonCache, Jac::Jacobian, F, xᵢ;
                            tol::Float64=1e-7, ω::Float64=NaN, curr_accuracy::Float64=NaN,
                            simple_step::Bool=false, update_infos::Bool=true,
                            use_qr::Bool=false,
                            precision::PrecisionOption=PRECISION_ADAPTIVE)
    rᵢ, Δxᵢ = cache.rᵢ, cache.Δxᵢ
    if simple_step
        evaluate!(rᵢ, F, xᵢ, cache.system_cache)
    else
        evaluate_and_jacobian!(rᵢ, Jᵢ, F, xᵢ, cache.system_cache)
        updated_jacobian!(Jac, update_infos=(update_infos || use_qr))
    end
    solve!(Δxᵢ, Jac, rᵢ, update_digits_lost=update_infos)
end

@inline function newton_step!(cache::NewtonCacheF64, Jac::Jacobian, F, xᵢ;
                    tol::Float64=1e-7, ω::Float64=NaN, curr_accuracy::Float64=NaN,
                    simple_step::Bool=false, update_infos::Bool=true,
                    use_qr::Bool=false,
                    precision::PrecisionOption=PRECISION_ADAPTIVE)
    Jᵢ = Jac.J
    @unpack rᵢ, rᵢ_D64, Δxᵢ, xᵢ_D64 = cache
    if simple_step
        if higher_precision(Jac, tol, curr_accuracy, ω, precision)
            xᵢ_D64 .= xᵢ
            evaluate!(rᵢ, F, xᵢ_D64, cache.system_cache_D64)
        else
            evaluate!(rᵢ, F, xᵢ, cache.system_cache)
        end
    else
        if higher_precision(Jac, tol, curr_accuracy, ω, precision)
            xᵢ_D64 .= xᵢ
            evaluate!(rᵢ, F, xᵢ_D64, cache.system_cache_D64)
            jacobian!(Jᵢ, F, xᵢ, cache.system_cache)
        else
            evaluate_and_jacobian!(rᵢ, Jᵢ, F, xᵢ, cache.system_cache)
        end
        updated_jacobian!(Jac, update_infos=(update_infos || use_qr))
    end
    solve!(Δxᵢ, Jac, rᵢ, update_digits_lost=update_infos)
end

@inline function higher_precision(Jac::Jacobian{ComplexF64}, tol::Float64,
                            curr_accuracy::Float64, ω::Float64,
                            precision::PrecisionOption)
    if precision == PRECISION_FIXED_64
        return false
    elseif precision == PRECISION_FIXED_128
        return true
    end
    digits_lost = unpack(Jac.digits_lost, 0.0)
    if !isfinite(curr_accuracy) || !isfinite(ω) || iszero(ω)
        digits_lost > 12 + log₁₀(tol)
    else
        digits_lost > 10 + log₁₀(max(0.5ω * curr_accuracy^2, tol))
    end
end

ext: jl


text:
export
AVLData,
BalancedTree,
BalanceData,
AVLTree,
AVLData,
avltree,
avlitree,
avldata,
weigh,
height,
balance,
maxheight
#build_wrapped_node

abstract type BalanceTrait end
struct Unbalanced <: BalanceTrait end
struct AVLBalanced <: BalanceTrait end
struct WAVLBalanced <: BalanceTrait end
struct RAVLBalanced <: BalanceTrait end

abstract type BalanceData{T} <: WrappedData{T} end
struct AVLData{T} <: BalanceData{T}
  kv::T
  rank::Int64
end

wrapped_data(n::AVLData{T}) where T = n.kv

AVLData(kv::T, h::S) where {T,S<:Integer} = AVLData(kv,Int64(h))
AVLData(kv::T) where T = AVLData(kv,Int64(0))
AVLData{T}(kv::T) where T = AVLData(kv,Int64(0))

const AVLNode{T} = Union{Node{AVLData{T}}, IndexedBinaryTree{AVLData{T}}}

balance_trait(::Type) = Unbalanced()
balance_trait(::Type{EmptyNode}) = Unbalanced()
balance_trait(::Type{<:AVLNode{T}}) where T = AVLBalanced()

avltree( d::T ) where T = tree( AVLData(d) )
avltree( dc::Vector{T} ) where T = tree( AVLData.(dc) )
avldata( kv::T ) where T = AVLData(kv)
avlitree( d::T ) where T = itree( AVLData(d) )
avlitree( dc::Vector{T} ) where T = itree( AVLData.(dc) )

""" maximum height of avltree with n leaves """
maxheight(n) = 1.440*log(2,n+1.065) - 0.328

function rotate_right( c::N ) where {T, N<:AVLNode{T}}
  b = left(c)
  indata = rebuild_parent_data( Val(:either), keyval(right(b)), keyval(right(c)), keyval(c) )
  cdata = AVLData(indata, height(right(b),right(c)))
  c = rebuild_node(Val(:left), right(b), right(c), c, cdata)

  indata = rebuild_parent_data( Val(:either), keyval(left(b)), keyval(c), keyval(b) )
  rebuild_node(Val(:right), left(b), c, b, AVLData(indata, height(left(b),c)))
end

function rotate_left( c::N ) where {T, N<:AVLNode{T}}
  b = right(c)
  indata = rebuild_parent_data( Val(:either), keyval(left(c)), keyval(left(b)), keyval(c) )
  c = rebuild_node(Val(:right), left(c), left(b), c, AVLData(indata, height(left(c),left(b))))

  indata = rebuild_parent_data( Val(:either), keyval(c), keyval(right(b)), keyval(b) )
  rebuild_node(Val(:left), c, right(b), b, AVLData(indata, height(c,right(b))))
end

balance( n::T ) where T = balance( n, balance_trait(T) )

balance( n::T, ::Unbalanced ) where T = n

function balance( st::NODE, trait::AVLBalanced ) where {T, NODE<:AVLNode{T}}
  scale = weigh(st)
  if scale > 1
    # left heavy
    subscale = weigh(left(st))
    if subscale > 0
      # left child heavy
      rotate_right(st)
    else
      # right child heavy
      rotate_right(rebuild_node(:left, rotate_left(left(st)), right(st), st))
    end
  elseif scale < -1
    # right heavy
    subscale = weigh(right(st))
    if subscale > 0
      # left child heavy
      rotate_left(rebuild_node(:right, left(st), rotate_right(right(st)), st))
    else
      # right child heavy
      rotate_left(st)
    end
  else
    st
  end
end

function rebuild_node(::Val{:left},
                      child::Branch{AVLData{T}},
                      _::Branch{AVLData{T}},
                      n::NODE, d::T ) where {T, NODE<:AVLNode{T}}
  avldata = AVLData{T}(d, max(height(right(n)),height(child))+1)
  balance(NODE(avldata, child, right(n)))
end

function rebuild_node(::Val{:right},
                      _::Branch{AVLData{T}},
                      child::Branch{AVLData{T}},
                      n::NODE, d::T ) where {T, NODE<:AVLNode{T}}
  avldata = AVLData{T}(d, max(height(left(n)),height(child))+1)
  balance(NODE(avldata,left(n),child))
end

function rebuild_node(::Val{:both},
                      lc::Branch{AVLData{T}},
                      rc::Branch{AVLData{T}},
                      n::NODE, d::T ) where {T, NODE<:AVLNode{T}}
  avldata = AVLData{T}(d, max(height(lc), height(rc))+1)
  balance(NODE(avldata,lc,rc))
end

function rebuild_node(dir::Val{:left},
                      lc::Branch{AVLData{T}},
                      rc::Branch{AVLData{T}},
                      n::IndexedBinaryTree{AVLData{T}}, d::T ) where {T}
  avldata = AVLData{T}(d, max(height(right(n)),height(lc))+1)
  balance(rebuild_node!(dir, lc, right(n), n, avldata))
end

function rebuild_node(dir::Val{:right},
                      lc::Branch{AVLData{T}},
                      rc::Branch{AVLData{T}},
                      n::IndexedBinaryTree{AVLData{T}}, d::T ) where {T}
  avldata = AVLData{T}(d, max(height(left(n)),height(rc))+1)
  balance(rebuild_node!(dir, left(n), rc, n, avldata))
end

function rebuild_node(dir::Val{:both},
                      lc::Branch{AVLData{T}},
                      rc::Branch{AVLData{T}},
                      n::IndexedBinaryTree{AVLData{T}}, d::T ) where {T}
  avldata = AVLData{T}(d, max(height(lc), height(rc))+1)
  balance(rebuild_node!(dir, lc, rc, n, avldata))
end

function ileaves(t::IndexedBinaryTree{AVLData{T}}) where T
  L = Vector{IndexedNode{AVLData{T}}}(undef, 1<<height(t))
  resize!(L, 0);
  ileaves(t.nodes, t.node, L);
  return L
end

update(n::Branch{AVLData{T}}) where T = balance(n)

weigh(n::NODE) where {T, NODE<:AVLNode{T}} = height(left(n))-height(right(n))
weigh(n::Nothing) = 0

key(d::AVLData{T}) where T = key(d.kv)
value(d::AVLData{T}) where T = value(d.kv)

keyval(n::NODE) where {T, NODE<:AVLNode{T}} = data(n).kv

height(n::NODE) where {T, NODE<:AVLNode{T}} = data(n).rank
height(n::EmptyNode) = -1

height(l::Branch{T}, r::Branch{T}) where T = max(height(l),height(r)) + 1

#!me why do I need to wrap? is that a key or a value?  should delete work on key(n)?
#delete(n::AVLNode{T}, k::T) where T = delete(n, avldata(k))

import Base.show
show(io::IO, n::AVLData{T}) where T = print(string(n.kv,"↑",n.rank ))


text:
## generate data from a PH model with constant baseline hazard log(2)/12
## return object of type Ph.Data, PhSpline.Data or WeibullPh.Data
function ph_exp_1(par, fpar, data_type=WeibullPHData)
    n = fpar.n

    z = zeros(Float64, n, 1)
    z[:,1] = rand(Binomial(1, 0.5), n) .- 0.5

    ## 10 patients per time unit
    max_entry = n / 15 ##fpar.max_entry ##/ 10

    entry = rand(Uniform(0, max_entry), n)
    U = -log.(1 .-rand(Uniform(0, 1), n))

    ## 12 time units median survival time
    T = (U./(par.scale .* exp.(z*par.beta))).^(1/par.shape)
    ##T = U.*exp.(-par.beta*z)./(log(2)/12)

    ## approx: 10% admin cens, 10% drop-out
    C = rand(Exponential(45/log(2)), n)
    ##C = min.(rand(Uniform(0, 3*36), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    data_type(id, entry, time, status, z)
end

function ph_biweibull_1(par, fpar, data_type=BiPHWeibullData)
    n = fpar.n

    z = rand(Binomial(1, 0.5), n) .- 0.5
    b = 1 .+ rand(Binomial(1, par.prev), n)

    ## 25 patients per month
    max_entry = n / 15 ## fpar.max_entry
    entry = rand(Uniform(0, max_entry), n)

    U = -log.(1 .-rand(Uniform(0, 1), n))
    T = zeros(Float64, n)

    for i in 1:n
        T[i] = (U[i] / (exp(z[i]*par.beta[b[i]])*par.scale[b[i]]))^(1/par.shape[b[i]])
    end

    C = rand(Exponential(45/log(2)), n)
    ##C = min.(rand(Uniform(0, 3*36), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    data_type(id, entry, time, status, b, z)
end

function ph_billog_1(par, fpar, data_type=BiPHSplineData)
    n = fpar.n
    z = rand(Binomial(1, 0.5), n) .- 0.5
    b = 1 .+ rand(Binomial(1, par.prev), n)

    ## 25 patients per month
    max_entry = n / 15 ## fpar.max_entry
    entry = rand(Uniform(0, max_entry), n)

    T = zeros(Float64, n)
    U = 1 .- rand(Uniform(0, 1), n)

    for i in 1:n
        T[i] = exp(quantile(Logistic(log(12), 1/6), 1 - exp(log(U[i])*exp(-par.beta[b[i]]*z[i]))))
    end

    ##C = min.(rand(Exponential(fpar.tmax), n), fpar.tmax)
    C = min.(rand(Uniform(0, 3*36), n), fpar.tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:fpar.n

    data_type(id, entry, time, status, b, z)
end

function getWeibullPHModel(max_n, delta0=[log(1.0), log(0.6)], delta=log(0.6); gen_data=ph_exp_1, tmax=Inf, emax=Inf)
    dshape = Uniform(0.25, 10)
    dscale = Gamma(1, 1)

    p = WeibullPHPrior(dshape, dscale, MvNormal(diagm(0 => 100*ones(eltype(delta), 1))))
    p0 = WeibullPHPrior(dshape, dscale, MultivariateDiracPM([delta0[1]]))
    p1 = WeibullPHPrior(dshape, dscale, MultivariateDiracPM([delta0[2]]))

    m = WeibullPHModel(WeibullPHParam(1.0, log(2)/5, [delta]),
                        WeibullPHFixedParam(max_n, tmax, emax),
                        (x,y) -> gen_data(x, y, WeibullPHData))

    m, p, p0, p1
end

function getBiPHWeibullModel(max_n, delta0=[log(1.0), log(0.6)], delta=(log(0.8), log(0.6)); gen_data=ph_biweibull_1, tmax=Inf, emax=Inf)##36.0, emax=48.0)
    dshape = Uniform(0.25, 10)

    p = BiPHWeibullPrior((dshape, dshape), (Gamma(1,1), Gamma(1,1)), MvNormal(10*diagm(0 => ones(Float64, 2))), Uniform(0.1, 0.9))
    p0 = BiPHWeibullPrior((dshape, dshape), (Gamma(1,1), Gamma(1,1)), MultivariateDiracPM([0.0, 0.0]), Uniform(0.1, 0.9))
    p1 = BiPHWeibullPrior((dshape, dshape), (Gamma(1,1), Gamma(1,1)), MultivariateDiracPM([log(0.8), log(0.6)]), Uniform(0.1, 0.9))

    m = Model(BiPHWeibullParam((1.0, 1.0), (log(2)/5, log(2)/5), delta, 1/3),
            BiPHWeibullFixedParam(max_n, tmax, emax), (x,y) -> gen_data(x, y, BiPHWeibullData))

    m, p, p0, p1
end

function getBiPHSplineModel(max_n, delta0=[log(1.0), log(0.6)], delta=(log(0.8), log(0.6)); gen_data=ph_billog_1, tmax=36.0, emax=Inf)

    sp = SplineHazard.Prior(Poisson(4), InverseGamma(0.01, 0.01), (theta, v) -> sum(theta.^2)/(2*v), 4, 100,
                                       Vector(0.0:tmax/501:tmax)[2:(end-1)], (0.0, tmax))

    p = BiPHSplinePrior((sp, sp), MvNormal(10*diagm(0 => ones(Float64, 2))), Uniform(0.1, 0.9))
    ##p0 = BiPHSpline.Prior((sp, sp), Dirac.MultivariateDiracPM([delta0[1], delta0[1]]), Uniform(0.1, 0.9))
    ##p1 = BiPHSpline.Prior((sp, sp), Dirac.MultivariateDiracPM([delta0[2], delta0[2]]), Uniform(0.1, 0.9))

    p0 = BiPHSplinePrior((sp, sp), MultivariateDiracPM([0.0, 0.0]), Uniform(0.1, 0.9))
    p1 = BiPHSplinePrior((sp, sp), MultivariateDiracPM([log(0.8), log(0.6)]), Uniform(0.1, 0.9))

    m = BiPHSplineModel(BiPHSplineParam((SplineHazard.Spline.const_spline(1.0, (0.0, tmax)),
                                        SplineHazard.Spline.const_spline(1.0, (0.0, tmax))),
                                                    delta, 1/3),
                                   BiPHSplineFixedParam(max_n, tmax, emax), (x,y) -> gen_data(x, y, BiPHSplineData))

    m, p, p0, p1
end

struct SimResult
    df::DataFrame
    nev::Float64
    asn::Float64
    ndistr::Vector{Float64}
    power_Sc::Float64
    power_S::Float64
    power_both::Float64
    power_any::Float64
    cov1::Float64
    cov2::Float64
    delta1::Float64
    delta2::Float64
end

function melt(x::Vector{Vector{SimResult}})
    R = length(x[1])
    res = Vector{Vector{SimResult}}(undef, R)

    for r in 1:R
        res[r] = [xx[r] for xx in x]
    end

    res
end

collapse(x::DataFrame) = [reduce(g) for g in groupby(x, :rule)]

function reduce(x::SubDataFrame)
    df = by(x, :ia) do zz
        [mean(zz[:nevents]) mean(zz[:n]) mean(zz[:time]) mean((zz[:rejectS] .== 1) .& (zz[:rejectSc] .== 1)) mean((zz[:rejectSc] .== 1) .& (zz[:rejectS] .== 0)) mean((zz[:rejectSc] .== 0) .& (zz[:rejectS] .== 1)) mean(zz[:beta1]) mean(zz[:beta2]) mean(zz[:sd1]) mean(zz[:sd2]) mean(zz[:pd1]) mean(zz[:pd2]) mean(zz[:cilo1]) mean(zz[:cihi1]) mean(zz[:cilo2]) mean(zz[:cihi2]) mean(zz[:cov1]) mean(zz[:cov2]) mean(zz[:u1]) mean(zz[:u2])]
    end

    names!(df, [:ia, :nevents, :n, :time, :reject_both, :reject_Sc_only, :reject_S_only, :beta1, :beta2, :sd1, :sd2, :pd1, :pd2, :cilo1, :cihi1, :cilo2, :cihi2, :cov1, :cov2, :u1, :u2])
    df.rule = repeat([x.rule[1]], inner=nrow(df))

    aux = by(x, :trial_id) do z
        [maximum(z.nevents) maximum(z.n) z.cov1[end] z.cov2[end] any(z.rejectSc .== 1) any(z.rejectS .== 1) (any(z.rejectSc .== 1) & any(z.rejectS .== 1)) (any(z.rejectSc .== 1) | any(z.rejectS .== 1))]
    end
    names!(aux, [:trial_id, :nevents, :n, :cov1, :cov2, :power_Sc, :power_S, :power_both, :power_any])

    SimResult(df, mean(aux.nevents), mean(aux.n), quantile(aux.n), mean(aux.power_Sc), mean(aux.power_S), mean(aux.power_both),
              mean(aux.power_any), mean(aux.cov1), mean(aux.cov2), x[:delta1][1], x[:delta2][1])
end

function summary(rule::String, x::Vector{SimResult})
    asn = [xx.asn for xx in x]
    power = [xx.power for xx in x]
    delta = [xx.delta for xx in x]
    DataFrame(rule=repeat([rule], inner=length(delta)), delta=delta, asn=asn, power=power)
end

summary(x::Dict{String, Vector{SimResult}}) = vcat([summary(v[1], v[2]) for v in x]...)

function get_fsim(B::Int, getmodel::Function, nevents=[[100, 50], [200, 100]]; max_n=400, M=5000,
                  delta0=[log(1.0), log(0.6)], delta=(log(0.8), log(0.6)), tau=0.0,
                  gen_data=ph_biweibull_1, warmup=Int(floor(M/2)), type1=0.05, type2=0.2)

    m, p, p0, p1 = getmodel(max_n, delta0, delta; gen_data=gen_data)
    tpar = GSDTrialParam(2, M, warmup, nevents, max_n, tau, DataFrame(a1=1.0, a2=1.0, b1=1.0, b2=1.0), p)
    get_initial_bounds(B, type1, type2, m, (p, p0, p1), tpar; fsim_only=true)
end

function max_power_fsim(fs)
    c = quantile(maximum(fs[2].post_pd[:,2,:], dims=2)[:,1], 0.95)
    mean(maximum(fs[3].post_pd[:,2,:], dims=2)[:,1] .> c)
end

function get_bounds(fsim, type1=0.05, type2=0.2; ngrid=10, h=0.001, fun=find_gsd_bounds, single_stage=false)
    objs = Dict("gsd" => gsd_utility, "sponsor" => sponsor_utility, "public" => public_utility)

    opt_bounds = @distributed vcat for obj in collect(values(objs))
        fun(type1, type2, fsim[1], fsim[2], fsim[3]; utility=obj, ngrid=ngrid, h=h, single_stage=single_stage)
    end

    DataFrame(rule=["OPT_$(k)" for k in keys(objs)], obj=map(x -> x[2][1], opt_bounds),
              a1=map(x -> x[1][1,1], opt_bounds), a2=map(x -> x[1][2,1], opt_bounds),
              b1=map(x -> x[1][1,2], opt_bounds), b2=map(x -> x[1][2,2], opt_bounds),
              t1=map(x -> x[3], opt_bounds), t2=map(x -> x[4], opt_bounds),
              t1_Sc=map(x -> x[5], opt_bounds), t2_Sc=map(x -> x[6], opt_bounds),
              t1_S=map(x -> x[7], opt_bounds), t2_S=map(x -> x[8], opt_bounds))
end

## run GSD simulation for each rule (each row of bounds DataFrame)
function gsd_sim(iter, bounds::DataFrame;
                 max_d=[[100, 50], [200, 100]], max_n=400, M=5000, delta0=[log(1.0), log(0.6)], delta=(0.0, 0.0),
                 model_fun=getBiPHWeibullModel, gen_data=ph_biweibull_1, warmup=Int(floor(M/2)), tau=0.0)

    rules = unique(bounds[:rule])
    res = Vector{Vector{SimResult}}(undef, 0)

    m, p, p0, p1 = model_fun(max_n, delta0, delta; gen_data=gen_data)
    tpar = GSDTrialParam(2, M, warmup, max_d, max_n, tau, bounds, p)
    tl = GSDTrial(m, tpar)

    trials = @distributed vcat for i in 1:iter
        df = run_trial(tl, i)
        df[:delta1] = delta[1]
        df[:delta2] = delta[2]
        df
    end

    push!(res, collapse(trials))
    dd = Dict(zip(rules, melt(res)))
    dd
end

function weibull_aed_sim(iter, bounds_df::DataFrameRow; delta=(log(1.0), log(0.6)), h=[(0.00025, true), (0.000250001, false)], utility=public_utility, B=10000, M=20000)
    bounds = [bounds_df.a1 bounds_df.b1; bounds_df.a2 bounds_df.b2]
    s1m, s1p, s1p0, s1p1 = getBiPHWeibullModel(400, [log(1.0), log(0.6)], delta)

    s2m_x, s2p_x, s2p0_x, s2p1_x = getWeibullPHModel(400, [log(1.0), log(0.8)], delta[2]; gen_data=ph_exp_1)
    s2m_y, s2p_y, s2p0_y, s2p1_y = getWeibullPHModel(400, [log(1.0), log(0.6)], delta[2]; gen_data=ph_exp_1)

    s2_dbl_m = Model(DoubleWeibullParam(s2m_x.param, s2m_y.param),
                    DoubleWeibullFixedParam(s2m_x.fparam, s2m_y.fparam),
                    s2m_x.generate_data)

    s2_dbl_p = DoubleWeibullPrior(s2p_x, s2p_y)
    s2_dbl_p0 = DoubleWeibullPrior(s2p0_x, s2p0_y)
    s2_dbl_p1 = DoubleWeibullPrior(s2p1_x, s2p1_y)

    s1tpar = AEDTrialParam(2, M, B, [[100, 50], [200, 100]], 0.0, bounds; h=h, utility=utility)
    s2tpar = AEDTrialParam(1, M, B, [[200, 100]], 0.0, bounds; h=h, utility=utility)

    t = AEDTrial(B,
                  s1m, Dict("full" => s1m, "double" => s2_dbl_m),
                  (s1p, s1p0, s1p1), Dict("full" => (s1p, s1p0, s1p1), "double" => (s2_dbl_p, s2_dbl_p0, s2_dbl_p1)),
                  s1tpar, Dict("full" => s2tpar, "double" => s2tpar))

    ##run_with_progress(iter, vcat, it -> AED.run_trial(t, it))
    @sync @distributed vcat for it in 1:iter
        @info "Iteration $(it)/$(iter)"
        run_trial(t, it)
    end
end


text:
mutable struct force_torque_t <: LCMType
    l_foot_force_z::Float32
    l_foot_torque_x::Float32
    l_foot_torque_y::Float32

    r_foot_force_z::Float32
    r_foot_torque_x::Float32
    r_foot_torque_y::Float32

    l_hand_force::SVector{3, Float32}
    l_hand_torque::SVector{3, Float32}

    r_hand_force::SVector{3, Float32}
    r_hand_torque::SVector{3, Float32}
end

@lcmtypesetup(force_torque_t)


text:
## Word selection

@assert extend_word("rand(", 2, true) == ("rand", 1, 4)
@assert extend_word("σ rand(", 4, true) == ("rand", 3, 6)

@assert extend_word("GtkIDE.Refactoring", 4, true) == ("GtkIDE", 1, 6)
@assert extend_word("GtkIDE.Refactoring", 4, false) == ("GtkIDE.Refactoring", 1, 18)
@assert extend_word("GtkIDE.Refactoring", 18, false) == ("GtkIDE.Refactoring", 1, 18)

@assert extend_word("sort!", 3, false) == ("sort!", 1, 5)
@assert extend_word(" @assert", 3, false) == ("@assert", 2, 8)

##
# Refactoring.arguments

args(x) = GtkIDE.Refactoring.arguments(x)

@assert args(Expr(:block,:x)) == Symbol[:x]

@assert args(:(x = 2)) == Symbol[]

@assert args(:(x = Int)) == Symbol[]

@assert args(:(x,y = 1,2)) == Symbol[]
@assert args(:(x,y,z,k = 1,2,3,4)) == Symbol[]

@assert args(:(x,y[1],z,k = 1,2,3,4)) == Symbol[]
@assert args(:(x,y[i],z,k = 1,2,3,4)) == Symbol[:i]

@assert args(:(x,y[i],z,k = 1,a,3,4)) == Symbol[:i,:a]

#@assert args(:(x,y[i[k]],z,k = 1,2,3,4)) == Symbol[:i,:k] #fail

@assert args(:(x = y)) == Symbol[:y]

#@assert args(:(f(x,y) =  x + y*a) ) == [:a]

@assert args(:(f(a,b,c))) == Symbol[:a,:b,:c]

@assert args(quote
    x,y,z,k = 1,2,3,4
    g = x*y
end) == Symbol[]

@assert args(quote
    x,y,z,k = 1,c,3,4
    g = x*d
end) == Symbol[:c,:d]

ex = quote
    using Gadfly
    mu = K1
    const alpha = 1
    const beta = (alpha+1)*mu

    for i=1:10
        x = rand(100,x,f(K2))
    end

    plot(x=K3,y=pdf(InverseGamma(K4,beta),K5),Geom.line)

    ind, K7 = K6.asd
    x = K8[K7]
end
@assert args(ex) == [:K1,:K2,:K3,:K4,:K5,:K6,:K8]

ex = quote
    for i=1:10
        x[i] = j
    end
end
@assert args(ex) == [:j]

ex = quote
    f(ex,b) = a
    ex = f(1,2)
end
@assert args(ex) == [:a]

# that's too advanced...
#ex = quote
#    y = [x for x = 1:10]
#end
#@assert args(ex) == Symbol[]

##


text:
using CFITSIO
using ImageTransformations

println(libcfitsio_version())

# src = homedir() * "/Downloads/ALMA00000006.fits"
# dst = homedir() * "/upsampled.fits"

src = "/Volumes/OWC/fits_web_ql/FITSCACHE/ALMA01567567.fits"
dst = "/Volumes/OWC/I_AM_BIG.fits"

f = fits_open_diskfile(src)

num = fits_get_num_hdus(f)
println("Number of HDUs in the file: ", num)

for i = 1:num
    hdu_type = fits_movabs_hdu(f, i)
    println(i, ") hdu_type = ", hdu_type)
end

bitpix, = fits_read_keyword(f, "BITPIX")
println("BITPIX = ", bitpix)
bitpix = parse(Int64, bitpix)

if bitpix != -32
    println("bitpix must be == -32")
    close(f)
    exit()
end

width, = fits_read_keyword(f, "NAXIS1")
height, = fits_read_keyword(f, "NAXIS2")
depth, = fits_read_keyword(f, "NAXIS3")

width = parse(Int64, width)
height = parse(Int64, height)
depth = parse(Int64, depth)

println("width: $(width), height: $(height), depth: $(depth)")

# new dimensions
new_width = 2 * width
new_height = 2 * height

println("new width: $new_width, new_height: $new_height")

keysexist, morekeys = fits_get_hdrspace(f)

out = fits_clobber_file(dst)

for i = 1:keysexist+1 # the extra one is the <END> keyword
    rec = fits_read_record(f, i)
    name, value, comment = fits_read_keyn(f, i)

    # println(rec)
    # println(name, "|", value, "|", comment)

    if name == "NAXIS1"
        fits_write_key(out, "NAXIS1", new_width, "")
        continue
    end

    if name == "NAXIS2"
        fits_write_key(out, "NAXIS2", new_height, "")
        continue
    end

    fits_write_record(out, rec)
end

println("FITS image dimensions: ", fits_get_img_size(f))

data = Array{Float32}(undef, width, height)
new_data = Array{Float32}(undef, new_width, new_height)

nelements = width * height

for frame = 1:depth
    global new_data

    fpixel = [1, 1, frame, 1]
    fits_read_pix(f, fpixel, nelements, data)

    println("HDU $(frame): ", size(data))

    new_data = Float32.(imresize(data, (new_width, new_height)))

    fits_write_pix(out, fpixel, new_width * new_height, new_data)
end

println("upsampled size:", size(new_data))

fits_close_file(f)
fits_close_file(out)


text:
function  mw_lfst_pv(x_t, itrend, param_m, param_t, nsim)
# mwlfst_pv - rewrite of Muller and Watson (2015) GAUSS code. Calculates LFST
#   test pvalues.
# Syntax:  pv::Float64 = mwlfst_pv(x_t, itrend, param_m, param_t, nsim)
#
# Inputs:
#   x_t::Vector{Float64} - Discrete cosine transforms of data.
#   itrend::Int64 - 1 if x_t is demeaned, 2 if x_t is detrended.
#   param_m::Float64 - "g" parameter for demeaned alternative.
#   param_t::Float64 - "g" parameter for detrended alternative.
#   nsim::Int64 - Number of simulations to compute p-value
# Outputs:
#   pv::Float64 - p-value
#------------- BEGIN CODE --------------
t = 1000 # Sample Size Used for Asymptotic Approximation
q = size(x, 1)
# Compute Null and Alternative Covariance Matrices
if itrend == 1
  psi = psi_compute(t, q)
  cov_h0 = psi'psi
  g_t = param_m/t
  gt2 = g_t.^2
  tmp = mwrw_cov(t)
  tmp = gt2*tmp
  tmp1 = diag(tmp) + ones(t, 1)
  cov_u = diagrv(tmp, tmp1)
  cov_ha = psi'cov_u*psi
end
# Compute Test Statistic @
  chol_cov_h0 = chol(cov_h0)
  a = inv(chol_cov_h0')
  x_t_trans = a*x_t # @ X transformed under Ho @
  cov_ha_trans = a*cov_ha*a' # @ Covariance under Ha Transformed @
#   @ Note ... LF under Ho is unity @
  lr = mw_marg_lf(x_t_trans, cov_ha_trans)
# @ Simulations for Pvalue @
  xsim = rand(Normal(), size(x,1), nsim)
  lr_sim = mw_marg_lf(xsim, cov_ha_trans)
  pv = mean(lr_sim .>(lr - 0.0001))
return lr, pv
end

function diagrv(x, v)
  y = x - diagm(diag(x)) + diagm(v[:])
  return y
end


text:
module LightDFGs

using LightGraphs
using DocStringExtensions
using UUIDs

using ...DistributedFactorGraphs

# import ...DistributedFactorGraphs: AbstractDFG, DFGNode, AbstractDFGVariable, AbstractDFGFactor, DFGSummary, AbstractParams, NoSolverParams, DFGVariable, DFGFactor

# import DFG functions to extend
import ...DistributedFactorGraphs:  setSolverParams!,
                                    getFactor,
                                    setDescription!,
                                    # getLabelDict,
                                    getUserData,
                                    setUserData!,
                                    getRobotData,
                                    setRobotData!,
                                    getSessionData,
                                    setSessionData!,
                                    addVariable!,
                                    getVariable,
                                    getAddHistory,
                                    addFactor!,
                                    getSolverParams,
                                    exists,
                                    isVariable,
                                    isFactor,
                                    getDescription,
                                    updateVariable!,
                                    updateFactor!,
                                    deleteVariable!,
                                    deleteFactor!,
                                    getVariables,
                                    listVariables,
                                    ls,
                                    getFactors,
                                    listFactors,
                                    lsf,
                                    isConnected,
                                    getNeighbors,
                                    buildSubgraph,
                                    copyGraph!,
                                    getBiadjacencyMatrix,
                                    _getDuplicatedEmptyDFG,
                                    toDot,
                                    toDotFile,
                                    findShortestPathDijkstra

include("FactorGraphs/FactorGraphs.jl")
using .FactorGraphs

# export SymbolEdge, is_directed, has_edge
# Imports
include("entities/LightDFG.jl")
include("services/LightDFG.jl")

# Exports
export LightDFG

end


text:
###
### A generic iterator over sequence positions that satisfy some function
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

struct ConditionIterator{F<:Function,S<:BioSequence}
    f::F
    seq::S
end

"""
    each(f::Function, seq::BioSequence) = ConditionIterator(f, seq)

Return an iterator over each element in `seq` that satisfies some function `f`.

This is similar to the `Iterators.filter` iterator, except that every
iteration yields a tuple of the position of the residue that satisfied `f`, and
the residue itself.
"""
each(f::Function, seq::BioSequence) = ConditionIterator(f, seq)

Base.eltype(::Type{ConditionIterator{<:Function,S}}) where {S<:BioSequence} = Tuple{Int,eltype(S)}
Base.IteratorSize(::ConditionIterator) = Base.SizeUnknown()

###
### Generic iteration for any sequence
###

function Base.iterate(it::ConditionIterator)
    return iterate(it, find_next_site(it.f, it.seq, firstindex(it.seq)))
end
function Base.iterate(it::ConditionIterator, nextpos::Int)
    if nextpos == 0
        return nothing
    else
        return (nextpos, inbounds_getindex(it.seq, nextpos)), find_next_site(it.f, it.seq, nextpos + 1)
    end
end

function find_next_site(f::Function, seq::BioSequence, start::Integer)
    for i in max(1, start):lastindex(seq)
        if f(inbounds_getindex(seq, i))
            return i
        end
    end
    return 0
end


text:
import GPLinearODEMaker.powers_of_negative_one

"""
m52_kernel function created by kernel_coder(). Requires 1 hyperparameters. Likely created using m52_kernel_base() as an input.
Use with include("src/kernels/m52_kernel.jl").
hyperparameters == ["λ"]
"""
function m52_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==1+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 3 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:2 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ = hyperparameters_view[1]

    dabs_δ = powers_of_negative_one(δ < 0)  # store derivative of abs()
    abs_δ = abs(δ)

    if dorder_view==[6, 2]
        func = 52500.0*exp(-2.23606797749979*abs_δ/λ)/λ^8 + 6250.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2/λ^10 - 39131.1896062463*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^9 + 5250.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^8 - 1500.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^7 - 10062.3058987491*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^7 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ^6 + 3354.10196624969*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^6 - 335.410196624969*exp(-2.23606797749979*abs_δ/λ)*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3)/λ^5 + 625.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^10 - 3913.11896062463*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^9 - 1677.05098312484*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^9 + 559.016994374948*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^8 + 9000.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^8 - 1500.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^7
    end

    if dorder_view==[5, 2]
        func = -11180.339887499*exp(-2.23606797749979*abs_δ/λ)/λ^7 - 1863.38998124983*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2/λ^9 + 10000.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^8 - 1677.05098312484*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^7 + 559.016994374948*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^6 + 2500.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^6 - 55.9016994374948*exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ^5 - 1000.0*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^5 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3)/λ^4 - 279.508497187474*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^9 + 1500.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^8 + 625.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^8 - 250.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^7 - 2795.08497187474*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^7 + 559.016994374948*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^6
    end

    if dorder_view==[4, 2]
        func = 2000.0*exp(-2.23606797749979*abs_δ/λ)/λ^6 + 500.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2/λ^8 - 2236.06797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^7 + 500.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 - 200.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^5 - 536.65631459995*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5 + 25.0*exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ^4 + 268.328157299975*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^4 - 44.7213595499958*exp(-2.23606797749979*abs_δ/λ)*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3)/λ^3 + 125.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^8 - 559.016994374948*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^7 - 223.606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^7 + 111.80339887499*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^6 + 800.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^6 - 200.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^5
    end

    if dorder_view==[3, 2]
        func = -268.328157299975*exp(-2.23606797749979*abs_δ/λ)/λ^5 - 111.80339887499*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2/λ^7 + 400.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^6 - 134.164078649987*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 67.0820393249937*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^4 + 90.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4 - 11.180339887499*exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ^3 - 60.0*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^3 + 15.0*exp(-2.23606797749979*abs_δ/λ)*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3)/λ^2 - 55.9016994374948*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^7 + 200.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 + 75.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^6 - 50.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^5 - 201.246117974981*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5 + 67.0820393249937*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^4
    end

    if dorder_view==[2, 2]
        func = 20.0*exp(-2.23606797749979*abs_δ/λ)/λ^4 + 16.6666666666667*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2/λ^6 - 44.7213595499958*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^5 + 30.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 - 20.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^3 + 5.0*exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ^2 + 25.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 - 67.0820393249937*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 22.3606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^4 - 8.94427190999916*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^3 + 8.94427190999916*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^2 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3)/λ - 22.3606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5 + 40.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4 - 20.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^3
    end

    if dorder_view==[1, 2]
        func = exp(-2.23606797749979*abs_δ/λ)*(16.6666666666667*abs_δ*dabs_δ/λ^4 + 4.47213595499958*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^3) + 5.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^3 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^3 + 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^2 + 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^2 - 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3)/λ - 11.180339887499*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 20.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 - 10.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^3
    end

    if dorder_view==[0, 2]
        func = exp(-2.23606797749979*abs_δ/λ)*(6.66666666666667*abs_δ^2/λ^4 + 4.47213595499958*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^3) + 5.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ^2*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^3 + 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^2
    end

    if dorder_view==[6, 1]
        func = -7500.0*exp(-2.23606797749979*abs_δ/λ)/λ^7 + 2795.08497187474*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^8 - 750.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^7 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^6 + 1677.05098312484*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^6 - 335.410196624969*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^5 + 279.508497187474*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^8 - 750.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^7
    end

    if dorder_view==[5, 1]
        func = 1863.38998124983*exp(-2.23606797749979*abs_δ/λ)/λ^6 - 833.333333333334*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^7 + 279.508497187474*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 - 55.9016994374948*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^5 - 500.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^4 - 125.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^7 + 279.508497187474*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^6
    end

    if dorder_view==[4, 1]
        func = -400.0*exp(-2.23606797749979*abs_δ/λ)/λ^5 + 223.606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^6 - 100.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 25.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^4 + 134.164078649987*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4 - 44.7213595499958*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^3 + 55.9016994374948*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 - 100.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5
    end

    if dorder_view==[3, 1]
        func = 67.0820393249937*exp(-2.23606797749979*abs_δ/λ)/λ^4 - 50.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^5 + 33.5410196624969*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 - 11.180339887499*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^3 - 30.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^3 + 15.0*exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ^2 - 25.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 33.5410196624969*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4
    end

    if dorder_view==[2, 1]
        func = -6.66666666666667*exp(-2.23606797749979*abs_δ/λ)/λ^3 + 7.4535599249993*exp(-2.23606797749979*abs_δ/λ)*abs_δ/λ^4 - 10.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^3 + 5.0*exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ^2 + 11.180339887499*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 + 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^2 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2)/λ - 10.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^3
    end

    if dorder_view==[1, 1]
        func = exp(-2.23606797749979*abs_δ/λ)*(-5.0*abs_δ*dabs_δ/λ^3 - 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ^2) + 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^2 + 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^2 - 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2)/λ - 5.0*exp(-2.23606797749979*abs_δ/λ)*abs_δ*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^3
    end

    if dorder_view==[0, 1]
        func = exp(-2.23606797749979*abs_δ/λ)*(-1.66666666666667*abs_δ^2/λ^3 - 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ^2) + 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*abs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^2
    end

    if dorder_view==[6, 0]
        func = 1250.0*exp(-2.23606797749979*abs_δ/λ)/λ^6 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^6 - 335.410196624969*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^5
    end

    if dorder_view==[5, 0]
        func = -372.677996249965*exp(-2.23606797749979*abs_δ/λ)/λ^5 - 55.9016994374948*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^5 + 125.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^4
    end

    if dorder_view==[4, 0]
        func = 100.0*exp(-2.23606797749979*abs_δ/λ)/λ^4 + 25.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^4 - 44.7213595499958*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^3
    end

    if dorder_view==[3, 0]
        func = -22.3606797749979*exp(-2.23606797749979*abs_δ/λ)/λ^3 - 11.180339887499*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^3 + 15.0*exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ^2
    end

    if dorder_view==[2, 0]
        func = 3.33333333333333*exp(-2.23606797749979*abs_δ/λ)/λ^2 + 5.0*exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ^2 - 4.47213595499958*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ)/λ
    end

    if dorder_view==[1, 0]
        func = exp(-2.23606797749979*abs_δ/λ)*(1.66666666666667*abs_δ*dabs_δ/λ^2 + 2.23606797749979*(1 + 0.74535599249993*abs_δ/λ)*dabs_δ/λ) - 2.23606797749979*exp(-2.23606797749979*abs_δ/λ)*dabs_δ*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)/λ
    end

    if dorder_view==[0, 0]
        func = exp(-2.23606797749979*abs_δ/λ)*(1 + 2.23606797749979*abs_δ*(1 + 0.74535599249993*abs_δ/λ)/λ)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return m52_kernel, 1  # the function handle and the number of kernel hyperparameters