text:
import Base.show

function hline(io,width::Vector)
    print(io,"├")
    for i=1:length(width)
        print(io,repeat("─",width[i]+2))
        if i < length(width)
            print(io,"┼")
        end
    end
    print(io,"┤\n")
end

function cpad(str::String,w::Int)
    len = length(str)
    if len > w
        return str[1:w]
    end

    d = w - len
    lsp = rsp = floor(Int,d/2)
    if d % 2 != 0
        rsp += 1
    end
    return string(repeat(" ",lsp),str,repeat(" ",rsp))
end

function strpr(n::Real,d::Int)
    if d == 2
        return @sprintf("%.1f",n)
    elseif d == 3
        return @sprintf("%.2f",n)
    elseif d == 4
        return @sprintf("%.3f",n)
    elseif d == 5
        return @sprintf("%.4f",n)
    elseif d == 6
        return @sprintf("%.5f",n)
    elseif d == 7
        return @sprintf("%.6f",n)
    end
end

function _getwidths(hdr,varray)
    width = length.(hdr)
    width[1] = max(width[1],maximum(length.(string.(varray[1]))),8) # 8 is for "combined"
    width[2] = max(width[2],maximum(length(string.(varray[2]))),6) # 6 is minimum width for any column
    mlen = max(length(string(maximum(maximum.(round.(Int,varray[3:end]))))),
        length(string(minimum(minimum.(round.(Int,varray[3:end]))))))
    digits = mlen > 5 ? 3 : mlen > 4 ? 4 : mlen > 3 ? 4 : 5 # minimum two decimal places
    for i=3:length(varray) # varray needs to be vector of arrays
        mlen = maximum(length.(string.(round.(Int,varray[i]))))
        width[i] = max(width[i],8,mlen + digits)
    end
    return (width,digits)
end

Base.show(io::IO,::MIME"text/plain",t::TTReturn) = show(io,t)
function show(io::IO,t::TTReturn)

    # compute widths for each column
    width,digits = _getwidths(t.colnms,t.array)

    tlen = sum(width) + 3*length(width) + 1 # 2 spaces before and after `|` + 1 boundary char

    # title
    println(io,t.title)

    # header row
    print(io,"| ",lpad(t.colnms[1],width[1]))
    for i=2:length(t.colnms)
        print(io," | ",lpad(t.colnms[i],width[i]))
    end
    print(io," |\n")

    # line
    hline(io,width)

    # Array values
    alen = t.title == "One-sample t test" ? 1 : 2
    for j=1:alen
        print(io,"| ",lpad(t.array[1][j],width[1]))
        print(io," | ",lpad(t.array[2][j],width[2]))
        for i=3:length(t.array)
            print(io," | ",lpad(strpr(t.array[i][j],digits), width[i]))
        end
        print(io," |\n")
    end

    # line
    if t.title != "One-sample t test"
        hline(io,width)
    end

    # combined - do not print if paired t test
    if t.paired == false && t.title != "One-sample t test"
        print(io,"| ",lpad("combined",width[1]))
        print(io," | ",lpad(t.array[2][3],width[2]))
        for i=3:length(t.array)
            print(io," | ",lpad(strpr(t.array[i][3],digits),width[i]))
        end
        print(io," |\n")

        # line
        hline(io,width)
    end

    # diff
    if t.title != "One-sample t test"
        print(io,"| ",lpad("diff",width[1]))
        for i=2:length(t.array)
            if t.paired == false && i in [2,4]
                print(io," | ",repeat(" ",width[i]))
            else
                print(io," | ",lpad(strpr(t.array[i][4],digits),width[i]))
            end
        end
        print(io," |\n\n")
    end

    # print ttest information
    if t.title == "One-sample t test"
        println(io,"diff = mean(",t.array[1][1],")")
        println(io,"H₀: diff = ",t.μ0)
    else
        println(io,"diff = mean(",t.array[1][1],") - mean(",t.array[1][2],")")
        println(io,"H₀: diff = 0")
    end
    if t.welch == false
        println(io,"t = ",strpr(t.t,5)," (df = ",round(Int,t.dof),")")
    else
        println(io,"t = ",strpr(t.t,5)," (df = ",strpr(t.dof,5),")")
    end
    print(io,"\n")

    tlen3rd = floor(Int,tlen/3)
    addsp = tlen % 3

    d = t.title == "One-sample t test" ? string(t.μ0) : "0"
    println(io,rpad(string("Hₐ: diff < ",d),tlen3rd),
        cpad(string("Hₐ: diff != ",d),tlen3rd),
        lpad(string("Hₐ: diff > ",d),tlen3rd+addsp)
    )
    print(io,rpad(string("Pr(T < t) = ",strpr(t.p_left,5)),tlen3rd),
        cpad(string("Pr(|T| > |t|) = ",strpr(t.p_both,5)),tlen3rd),
        lpad(string("Pr(T > t) = ",strpr(t.p_right,5)),tlen3rd+addsp)
    )
    print(io,"\n\n")
end

Base.show(io::IO,::MIME"text/plain",x::XsqResult) = show(io,x)
function show(io::IO,x::XsqResult)
    println(io,"Pearson χ-square (",x.dof,") = ",@sprintf("%.5f",x.chisq),"   P-value = ",@sprintf("%.5f",x.p))
end

#
# Base.show(io::IO, ::MIME"text/plain", tr::Stella.tab_return) = show(io,tr)
#
# function show(io::IO,tr::Stella.tab_return)
#
#     if length(tr.dimnames) == 1
#         show_oneway(io,NamedArray(tr.array,tr.dicts,tr.dimnames), total = total, precision = precision)
#         exit()
#     elseif length(tr.dimnames) > 2
#         error("Only up to two dimensional arrays are currently supported")
#     end
#
#     dimnames = string(tr.dimnames[1]) * " \\ " * string(tr.dimnames[2])
#     print(io,dimnames,"\n")
#
#     # temporary
#     row = col = cell = false
#
#     # total is true when row, col, or cell is true
#     if row || col || cell
#       total = true
#     end
#
#     # row names
#     rownames = collect(keys(tr.dicts[1]))
#
#     maxrowname = 5
#     for i = 1:length(rownames)
#         maxrowname = max(maxrowname,length(rownames[i]))
#     end
#
#     # column names
#     colnames = collect(keys(tr.dicts[20]))
#
#     maxcolname = 3
#     for i = 1:length(colnames)
#         maxcolname = max(maxcolname,length(colnames[i]))
#     end
#
#     # width of data columns - the same as the length of tht grand total
#     tot = sum(tr.array) # grand total
#     colwidth = length(digits(Int(floor(tot))))
#
#     # number of columns
#     ncols = length(colnames)
#
#     # number of rows
#     nrows = length(rownames)
#
#     # floating point numbers with three digits after decimal point
#     if eltype(tr.array) <: AbstractFloat
#       colwidth += 3
#     end
#
#     # determine column widths
#     colwidth = max(maxcolname,colwidth)
#
#     #---------------------------------------------------
#     # print header
#     print(io,repeat(" ",maxrowname)," |")
#
#     for i = 1:length(colnames)
#         print(io," ",lpad(string(colnames[i]),colwidth))
#     end
#
#     if total
#       print(io," | ",lpad("Total",colwidth))
#     end
#     print(io,"\n")
#
#     print(io,repeat("-",maxrowname),"-+-",repeat("-",(colwidth+1)*(length(colnames))-1))
#
#     if total
#       print(io,"-+-",repeat("-",colwidth))
#     end
#     print(io,"\n")
#     #---------------------------------------------------
#
#     # column totals
#     colsum = sum(tr.array,1)
#
#     # row totals
#     rowsum = sum(tr.array,2)
#
#     #----------------------------------------------------
#     # print values
#     for i = 1:nrows
#
#         # row name
#         print(io,rpad(string(rownames[i]),maxrowname)," |")
#
#         for j = 1:ncols
#             val = strval(tr.array[i,j])
#             print(io," ",lpad(val,colwidth))
#         end
#
#         # row total
#         if total
#           print(io," |")
#
#           val = strval(rowsum[i])
#           print(io," ",lpad(val,colwidth))
#         end
#         print(io,"\n")
#
#         # row percentages
#         if row
#           print(io,repeat(" ",maxrowname)," |")
#           for j = 1:ncols
#               val = strval(100 * tr.array[i,j] / rowsum[i],precision)
#               print(io," ",lpad(val,colwidth," "))
#           end
#
#           # row percentage
#           print(io," |")
#
#           val = strval(100.0,precision)
#           print(io," ",lpad(val,colwidth),"\n")
#         end
#
#         # column percentages
#         if col
#           print(io,repeat(" ",maxrowname)," |")
#           for j = 1:ncols
#               val = strval(100 * tr.array[i,j] / colsum[j],precision)
#               print(io," ",lpad(val,colwidth))
#           end
#
#           # column percent
#           print(" |")
#
#           val = strval(100 * rowsum[i] / tot,precision)
#           print(io," ",lpad(val,colwidth),"\n")
#         end
#
#         # column percentages
#         if cell
#           print(io,repeat(" ",maxrowname)," |")
#           for j = 1:ncols
#               val = strval(100 * tr.array[i,j] / tot,precision)
#               print(io," ",lpad(val,colwidth))
#           end
#
#           # column percent
#           print(io," |")
#
#           val = strval(100 * rowsum[i] / tot,precision)
#           print(io," ",lpad(val,colwidth),"\n")
#         end
#
#         if row || col || cell
#           print(io,repeat("-",maxrowname),"-+=",repeat("-",(colwidth+1)*(length(colnames))),"+-",repeat("-",colwidth),"\n")
#         end
#
#     end
#
#     #----------------------------------------------------
#     # Total
#     if total
#       if !(row || col || cell)
#         print(io,repeat("-",maxrowname+1),"+",repeat("-",(colwidth+1)*(length(colnames))),"-+-",repeat("-",colwidth),"\n")
#       end
#       print(rpad("Total",maxrowname," ")," |")
#
#       for i = 1:length(colnames)
#           val=strval(colsum[i])
#           print(io," ",lpad(val,colwidth," "))
#       end
#
#       # Grand total
#       val = strval(tot)
#       print(io," | ",lpad(val,colwidth," "),"\n")
#     end
#     #----------------------------------------------------
#
#     # row percentages
#     if row
#       print(io,repeat(" ",maxrowname)," |")
#       for j = 1:ncols
#           val = strval(100 * colsum[j] / tot,precision)
#           print(io," ",lpad(val,colwidth," "))
#       end
#
#       # column percent
#       val = strval(100.0,precision)
#       print(io," | ",lpad(val,colwidth," "),"\n")
#     end
#
#     # column percentages
#     if col
#       print(io,repeat(" ",maxrowname)," |")
#       for j = 1:ncols
#           val = strval(100.,precision)
#           print(io," ",lpad(val,colwidth," "))
#       end
#
#       # column percent
#       val = strval(100.0,precision)
#       print(io," | ",lpad(val,colwidth," "),"\n")
#
#     end
#
#     # column percentages
#     if cell
#       print(io,repeat(" ",maxrowname)," |")
#       for j = 1:ncols
#           val = strval(100 * colsum[j] / tot,precision)
#           print(io," ",lpad(val,colwidth," "))
#       end
#
#       # column percent
#       val = strval(100.0,precision)
#       print(io," | ",lpad(val,colwidth," "),"\n")
#     end
#
#     # p-value
#     println(io,"\nPearson χ² (",tr.dof,") = ",@sprintf("%.3f",tr.chisq)," Pr = ",@sprintf("%.3f",tr.p))
# end
#
# function show_oneway(io::IO,na::NamedArrays.NamedArray; total = false, precision::Int8 = 2)
#
#     rownames = names(na,1)
#
#     maxrowname = max(5,maximum(length.(rownames)))
#
#     # width of data columns - the same as the length of tht grand total
#     tot = sum(na) # grand total
#     colwidth = max(7,length(digits(Int(floor(tot)))))
#
#     # print header
#     print(io,rpad("Value",maxrowname," "))
#     print(io," | ",lpad("Count",colwidth," ")," ",lpad("Percent",colwidth," "),
#         " ",lpad("Cum_pct",colwidth," "),"\n")
#     print(io,repeat("-",maxrowname),"-+-",repeat("-",colwidth),"-",repeat("-",colwidth),"-",repeat("-",colwidth),"\n")
#
#     # values
#     cumtot = cumsum(na.array)
#     for i = 1:size(na.array,1)
#         str = strval(na.array[i])
#         pct = strval(100*na.array[i] / tot, precision)
#         cumpct = strval(100*cumtot[i] / tot, precision)
#         print(io,rpad(string(rownames[i]),maxrowname," "),
#             " | ",lpad(str,colwidth," "),
#             " ",lpad(pct,colwidth," "),
#             " ",lpad(cumpct,colwidth," "),"\n")
#     end
#
#     # total
#     if total
#         print(io,repeat("-",maxrowname),"-+-",repeat("-",colwidth),"-",repeat("-",colwidth),"-",repeat("-",colwidth),"\n")
#         str = strval(tot)
#         pct = strval(100.00, precision)
#         cumpct = strval(100.00, precision)
#         print(io,rpad("Total",maxrowname," "),
#             " | ",lpad(str,colwidth," "),
#             " ",lpad(pct,colwidth," "),
#             " ",lpad(cumpct,colwidth," "),"\n")
#     end
# end
#
# Base.show(io::IO,::MIME"text/plain",pr::pwcorr_return) = show(io,pr)
#
# function show(io::IO,pr::pwcorr_return; width::Int8 = 9, precision::Int8 = 3, p = false, N = false)
#
#     ncol = size(pr.N,1)
#
#     for i = 1:ncol
#         print(io,prepend_spaces(pr.colnames[i],width))
#         if i < ncol
#             print(io," ")
#         end
#     end
#     print(io,"\n")
#
#     for i = 1:ncol
#         print(io,prepend_spaces(pr.colnames[i],width)," ")
#         for j=1:i
#             printf(io,"%9.4f ",bc[i,j])
#             if i == j
#                 print(io,"\n")
#             end
#         end
#         if p == true
#             print(io,repeat(" ",width+1))
#             for j=1:i
#                 printf(io,"%9.4f ",bp[i,j])
#                 if i == j
#                     print(io,"\n")
#                 end
#             end
#         end
#     end
# end
# #
# # function show(io::IO,tr::Stella.tab_return; row=false, col=false, cell=false, total=false, precision::Int8 = 2)
# #
# #     if length(tr.dimnames) == 1
# #         show_oneway(io,NamedArray(tr.array,tr.dicts,tr.dimnames), total = total, precision = precision)
# #         exit()
# #     elseif length(tr.dimnames) > 2
# #         error("Only up to two dimensional arrays are currently supported")
# #     end
# #
# #     dimnames = string(tr.dimnames[1]) * " \\ " * string(tr.dimnames[2])
# #     print(io,dimnames,"\n")
# #
# #     # total is true when row, col, or cell is true
# #     if row || col || cell
# #       total = true
# #     end
# #
# #     # row names
# #     rownames = collect(keys(tr.dicts[1]))
# #
# #     maxrowname = 5
# #     for i = 1:length(rownames)
# #         maxrowname = max(maxrowname,length(rownames[i]))
# #     end
# #
# #     # column names
# #     colnames = collect(keys(tr.dicts[20]))
# #
# #     maxcolname = 3
# #     for i = 1:length(colnames)
# #         maxcolname = max(maxcolname,length(colnames[i]))
# #     end
# #
# #     # width of data columns - the same as the length of tht grand total
# #     tot = sum(tr.array) # grand total
# #     colwidth = length(digits(Int(floor(tot))))
# #
# #     # number of columns
# #     ncols = length(colnames)
# #
# #     # number of rows
# #     nrows = length(rownames)
# #
# #     # floating point numbers with three digits after decimal point
# #     if eltype(tr.array) <: AbstractFloat
# #       colwidth += 3
# #     end
# #
# #     # determine column widths
# #     colwidth = max(maxcolname,colwidth)
# #
# #     #---------------------------------------------------
# #     # print header
# #     print(io,repeat(" ",maxrowname)," |")
# #
# #     for i = 1:length(colnames)
# #         print(io," ",lpad(string(colnames[i]),colwidth))
# #     end
# #
# #     if total
# #       print(io," | ",lpad("Total",colwidth))
# #     end
# #     print(io,"\n")
# #
# #     print(io,repeat("-",maxrowname),"-+-",repeat("-",(colwidth+1)*(length(colnames))-1))
# #
# #     if total
# #       print(io,"-+-",repeat("-",colwidth))
# #     end
# #     print(io,"\n")
# #     #---------------------------------------------------
# #
# #     # column totals
# #     colsum = sum(tr.array,1)
# #
# #     # row totals
# #     rowsum = sum(tr.array,2)
# #
# #     #----------------------------------------------------
# #     # print values
# #     for i = 1:nrows
# #
# #         # row name
# #         print(io,rpad(string(rownames[i]),maxrowname)," |")
# #
# #         for j = 1:ncols
# #             val = strval(tr.array[i,j])
# #             print(io," ",lpad(val,colwidth))
# #         end
# #
# #         # row total
# #         if total
# #           print(io," |")
# #
# #           val = strval(rowsum[i])
# #           print(io," ",lpad(val,colwidth))
# #         end
# #         print(io,"\n")
# #
# #         # row percentages
# #         if row
# #           print(io,repeat(" ",maxrowname)," |")
# #           for j = 1:ncols
# #               val = strval(100 * tr.array[i,j] / rowsum[i],precision)
# #               print(io," ",lpad(val,colwidth," "))
# #           end
# #
# #           # row percentage
# #           print(io," |")
# #
# #           val = strval(100.0,precision)
# #           print(io," ",lpad(val,colwidth),"\n")
# #         end
# #
# #         # column percentages
# #         if col
# #           print(io,repeat(" ",maxrowname)," |")
# #           for j = 1:ncols
# #               val = strval(100 * tr.array[i,j] / colsum[j],precision)
# #               print(io," ",lpad(val,colwidth))
# #           end
# #
# #           # column percent
# #           print(" |")
# #
# #           val = strval(100 * rowsum[i] / tot,precision)
# #           print(io," ",lpad(val,colwidth),"\n")
# #         end
# #
# #         # column percentages
# #         if cell
# #           print(io,repeat(" ",maxrowname)," |")
# #           for j = 1:ncols
# #               val = strval(100 * tr.array[i,j] / tot,precision)
# #               print(io," ",lpad(val,colwidth))
# #           end
# #
# #           # column percent
# #           print(io," |")
# #
# #           val = strval(100 * rowsum[i] / tot,precision)
# #           print(io," ",lpad(val,colwidth),"\n")
# #         end
# #
# #         if row || col || cell
# #           print(io,repeat("-",maxrowname),"-+=",repeat("-",(colwidth+1)*(length(colnames))),"+-",repeat("-",colwidth),"\n")
# #         end
# #
# #     end
# #
# #     #----------------------------------------------------
# #     # Total
# #     if total
# #       if !(row || col || cell)
# #         print(io,repeat("-",maxrowname+1),"+",repeat("-",(colwidth+1)*(length(colnames))),"-+-",repeat("-",colwidth),"\n")
# #       end
# #       print(rpad("Total",maxrowname," ")," |")
# #
# #       for i = 1:length(colnames)
# #           val=strval(colsum[i])
# #           print(io," ",lpad(val,colwidth," "))
# #       end
# #
# #       # Grand total
# #       val = strval(tot)
# #       print(io," | ",lpad(val,colwidth," "),"\n")
# #     end
# #     #----------------------------------------------------
# #
# #     # row percentages
# #     if row
# #       print(io,repeat(" ",maxrowname)," |")
# #       for j = 1:ncols
# #           val = strval(100 * colsum[j] / tot,precision)
# #           print(io," ",lpad(val,colwidth," "))
# #       end
# #
# #       # column percent
# #       val = strval(100.0,precision)
# #       print(io," | ",lpad(val,colwidth," "),"\n")
# #     end
# #
# #     # column percentages
# #     if col
# #       print(io,repeat(" ",maxrowname)," |")
# #       for j = 1:ncols
# #           val = strval(100.,precision)
# #           print(io," ",lpad(val,colwidth," "))
# #       end
# #
# #       # column percent
# #       val = strval(100.0,precision)
# #       print(io," | ",lpad(val,colwidth," "),"\n")
# #
# #     end
# #
# #     # column percentages
# #     if cell
# #       print(io,repeat(" ",maxrowname)," |")
# #       for j = 1:ncols
# #           val = strval(100 * colsum[j] / tot,precision)
# #           print(io," ",lpad(val,colwidth," "))
# #       end
# #
# #       # column percent
# #       val = strval(100.0,precision)
# #       print(io," | ",lpad(val,colwidth," "),"\n")
# #     end
# #
# #     # p-value
# #     println(io,"\nPearson χ² (",tr.dof,") = ",@sprintf("%.3f",tr.chisq)," Pr = ",@sprintf("%.3f",tr.p))
# # end
# #
# # function show_oneway(io::IO,na::NamedArrays.NamedArray; total = false, precision::Int8 = 2)
# #
# #     rownames = names(na,1)
# #
# #     maxrowname = max(5,maximum(length.(rownames)))
# #
# #     # width of data columns - the same as the length of tht grand total
# #     tot = sum(na) # grand total
# #     colwidth = max(7,length(digits(Int(floor(tot)))))
# #
# #     # print header
# #     print(io,rpad("Value",maxrowname," "))
# #     print(io," | ",lpad("Count",colwidth," ")," ",lpad("Percent",colwidth," "),
# #         " ",lpad("Cum_pct",colwidth," "),"\n")
# #     print(io,repeat("-",maxrowname),"-+-",repeat("-",colwidth),"-",repeat("-",colwidth),"-",repeat("-",colwidth),"\n")
# #
# #     # values
# #     cumtot = cumsum(na.array)
# #     for i = 1:size(na.array,1)
# #         str = strval(na.array[i])
# #         pct = strval(100*na.array[i] / tot, precision)
# #         cumpct = strval(100*cumtot[i] / tot, precision)
# #         print(io,rpad(string(rownames[i]),maxrowname," "),
# #             " | ",lpad(str,colwidth," "),
# #             " ",lpad(pct,colwidth," "),
# #             " ",lpad(cumpct,colwidth," "),"\n")
# #     end
# #
# #     # total
# #     if total
# #         print(io,repeat("-",maxrowname),"-+-",repeat("-",colwidth),"-",repeat("-",colwidth),"-",repeat("-",colwidth),"\n")
# #         str = strval(tot)
# #         pct = strval(100.00, precision)
# #         cumpct = strval(100.00, precision)
# #         print(io,rpad("Total",maxrowname," "),
# #             " | ",lpad(str,colwidth," "),
# #             " ",lpad(pct,colwidth," "),
# #             " ",lpad(cumpct,colwidth," "),"\n")
# #     end
# # end
# #
# # Base.show(io::IO,::MIME"text/plain",pr::pwcorr_return) = show(io,pr)
# #
# # function show(io::IO,pr::pwcorr_return; width::Int8 = 9, precision::Int8 = 3, p = false, N = false)
# #
# #     ncol = size(pr.N,1)
# #
# #     for i = 1:ncol
# #         print(io,prepend_spaces(pr.colnames[i],width))
# #         if i < ncol
# #             print(io," ")
# #         end
# #     end
# #     print(io,"\n")
# #
# #     for i = 1:ncol
# #         print(io,prepend_spaces(pr.colnames[i],width)," ")
# #         for j=1:i
# #             printf(io,"%9.4f ",bc[i,j])
# #             if i == j
# #                 print(io,"\n")
# #             end
# #         end
# #         if p == true
# #             print(io,repeat(" ",width+1))
# #             for j=1:i
# #                 printf(io,"%9.4f ",bp[i,j])
# #                 if i == j
# #                     print(io,"\n")
# #                 end
# #             end
# #         end
# #     end
# # end


text:
# text

# first, the 'toy' API... "Any serious application should avoid them." :)

"""
    text(str)
    text(str, pos)
    text(str, pos, angle=pi/2)
    text(str, x, y)
    text(str, pos, halign=:left)
    text(str, valign=:baseline)
    text(str, valign=:baseline, halign=:left)
    text(str, pos, valign=:baseline, halign=:left)

Draw the text in the string `str` at `x`/`y` or `pt`, placing the start of the
string at the point. If you omit the point, it's placed at the current `0/0`.

`angle` specifies the rotation of the text relative to the current x-axis.

Horizontal alignment `halign` can be `:left`, `:center`, (also `:centre`) or
`:right`.  Vertical alignment `valign` can be `:baseline`, `:top`, `:middle`, or
`:bottom`.

The default alignment is `:left`, `:baseline`.

This uses `textextents()` to query the dimensions of the text. This returns
values of the built in to the font. You can't find

This uses Cairo's Toy text API.
"""
function text(t, pt::Point;
        halign=:left,
        valign=:baseline,
        angle=0.0)
    #= text can aligned by one of the following points
        top/left       top/center       top/right
        middle/left    middle/center    middle/right
        baseline/left  baseline/center  baseline/right
        bottom/left    bottom/center    bottom/right

    # left center right x coords are
    # [:xbearing, :width-:xbearing/2, :width]
    # top middle baseline bottom y coords are
    # [:ybearing, :ybearing/2, 0, :height + :ybearing]
    =#

    xbearing, ybearing, textwidth, textheight, xadvance, yadvance = textextents(t)
    halignment = findfirst(isequal(halign), [:left, :center, :right, :centre])

    # if unspecified or wrong, default to left, also treat UK spelling centre as center
    if halignment == nothing
        halignment = 1
    elseif halignment == 4
        halignment = 2
    end

    # was textpointx = pt.x - [0, textwidth/2, textwidth][halignment]
    textpointx = pt.x - [0, xadvance/2, textwidth + xbearing][halignment]

    valignment = findfirst(isequal(valign), [:top, :middle, :baseline, :bottom])

    # if unspecified or wrong, default to baseline
    if valignment == nothing
        valignment = 3
    end

    textpointy = pt.y - [ybearing, ybearing/2, 0, textheight + ybearing][valignment]

    # need to adjust for any rotation now
    # rotate around original point
    finalpt = rotate_point_around_point(Point(textpointx, textpointy), pt, angle)

    gsave()
        translate(finalpt)
        rotate(angle)
        newpath()
        Cairo.show_text(get_current_cr(), t)
    grestore()

    return Point(textpointx, textpointy)
end

text(t; kwargs...) = text(t, O; kwargs...)
text(t, xpos, ypos; kwargs...) = text(t, Point(xpos, ypos); kwargs...)

# deprecated probably
textcentered(t, x=0, y=0) = text(t, x, y, halign=:center)
textcentered(t, pt::Point) = textcentered(t, pt.x, pt.y)
textcentred = textcentered
textright(t, x=0, y=0) = text(t, x, y, halign=:right)
textright(t, pt::Point) = textright(t, pt.x, pt.y)

"""
    fontface(fontname)

Select a font to use. (Toy API)
"""
fontface(f) =
    Cairo.select_font_face(get_current_cr(), f,
                           Cairo.FONT_SLANT_NORMAL,
                           Cairo.FONT_WEIGHT_NORMAL)

"""
    fontsize(n)

Set the font size to `n` units. The default size is 10 units. (Toy API)
"""
fontsize(n) = Cairo.set_font_size(get_current_cr(), n)

"""
    get_fontsize()

Return the font size set by `fontsize` or. more precisely. the y-scale of the Cairo font matrix
if `Cairo.set_font_matrix` is used directly. (Toy API)

> This only works if Cairo is at least at v1.0.5.
"""
function get_fontsize()
    if @isdefined get_font_matrix
        m = get_font_matrix(get_current_cr())
        font_size = sign(m.yy)*sqrt(m.yx^2+m.yy^2)
        return font_size
    else
        throw(MethodError(get_fontsize, "Please use Cairo v1.0.5 or later to use this feature."))
    end
end

"""
    textextents(str)

Return an array of six Float64s containing the measurements of the string `str`
when set using the current font settings (Toy API):

1 x_bearing

2 y_bearing

3 width

4 height

5 x_advance

6 y_advance

The x and y bearings are the displacement from the reference point to the
upper-left corner of the bounding box. It is often zero or a small positive
value for x displacement, but can be negative x for characters like "j"; it's
almost always a negative value for y displacement.

The width and height then describe the size of the bounding box. The advance
takes you to the suggested reference point for the next letter. Note that
bounding boxes for subsequent blocks of text can overlap if the bearing is
negative, or the advance is smaller than the width would suggest.

Example:

    textextents("R")

returns

    [1.18652; -9.68335; 8.04199; 9.68335; 9.74927; 0.0]
"""
textextents(str) = Cairo.text_extents(get_current_cr(), str)

"""
    textpath(t)

Convert the text in string `t` and adds closed paths to the current path, for subsequent filling/stroking etc...

Typically you'll have to use `pathtopoly()` or `getpath()` or `getpathflat()` then
work through the one or more path(s). Or use `textoutlines()`.
"""
function textpath(t)
    Cairo.text_path(get_current_cr(), t)
end

"""
    textoutlines(s::AbstractString, pos::Point=O, action::Symbol=:none;
        halign=:left,
        valign=:baseline,
        startnewpath=true)

Convert text to a graphic path and apply `action`.

By default this function discards any current path, unless you use `startnewpath=false`

See also `textpath()`.
"""
function textoutlines(s::AbstractString, pos::Point=O, action::Symbol=:none;
    halign=:left,
    valign=:baseline,
    startnewpath=true)

    # TODO this duplicates text() too much; re-factor needed
    xbearing, ybearing, textwidth, textheight, xadvance, yadvance = textextents(s)
    halignment = findfirst(isequal(halign), [:left, :center, :right, :centre])
    if halignment == nothing
        halignment = 1
    elseif halignment == 4
        halignment = 2
    end
    textpointx = pos.x - [0, xadvance/2, textwidth + xbearing][halignment]
    valignment = findfirst(isequal(valign), [:top, :middle, :baseline, :bottom])
    if valignment == nothing
        valignment = 3
    end
    textpointy = pos.y - [ybearing, ybearing/2, 0, textheight + ybearing][valignment]
    @layer begin
        translate(Point(textpointx, textpointy))
        if startnewpath
           newpath() # forget any current path, start a new one
        end
        te = textextents(s)
        textpath(s)
        tp = pathtopoly()
        if length(tp) == 1
            poly.(tp, :path, close=true) # don't clip yet
        else
            newpath()
            for path in tp
                poly(path, :path, close=true)
                newsubpath()
            end
            closepath()
        end
    end
    do_action(action)
end

"""
    textcurve(the_text, start_angle, start_radius, x_pos = 0, y_pos = 0;
          # optional keyword arguments:
          spiral_ring_step = 0,    # step out or in by this amount
          letter_spacing = 0,      # tracking/space between chars, tighter is (-), looser is (+)
          spiral_in_out_shift = 0, # + values go outwards, - values spiral inwards
          clockwise = true
          )

Place a string of text on a curve. It can spiral in or out.

`start_angle` is relative to +ve x-axis, arc/circle is centered on `(x_pos,y_pos)` with
radius `start_radius`.
"""
function textcurve(the_text, start_angle, start_radius, x_pos=0, y_pos=0;
    spiral_ring_step = 0,
    letter_spacing = 0,
    spiral_in_out_shift = 0,
    clockwise = true
    )
    # TODO this is all very hacky and out of date...
    refangle = start_angle
    current_radius = start_radius
    spiral_space_step = 0
    xx = 0
    yy = 0
    angle_step = 0
    radius_step = 0
    for i in the_text
        glyph = string(i)
        glyph_x_bearing, glyph_y_bearing, glyph_width,
        glyph_height, glyph_x_advance, glyph_y_advance = textextents(glyph)
        spiral_space_step = glyph_x_advance + letter_spacing
        cnter = (2pi * current_radius) / spiral_space_step
        radius_step = (spiral_ring_step + spiral_in_out_shift) / cnter
        current_radius += radius_step
        angle_step += (glyph_x_advance / 2.0) + letter_spacing/2.0
        if clockwise
            refangle += angle_step / current_radius
        else
            refangle -= angle_step / current_radius
        end
        angle_step = (glyph_x_advance / 2.0) + letter_spacing/2.0
        xx = cos(refangle) * current_radius + x_pos
        yy = sin(refangle) * current_radius + y_pos
        gsave()
        translate(xx, yy)
        if clockwise
            rotate(pi/2 + refangle)
        else
            rotate(-pi/2 + refangle)
        end
        text(glyph, O, halign=:center)
        grestore()
        current_radius < 10 && break
    end
end

textcurve(the_text, start_angle, start_radius, center::Point; kwargs...) =
    textcurve(the_text, start_angle, start_radius, center.x, center.y; kwargs...)

"""
    textcurvecentered(the_text, the_angle, the_radius, center::Point;
          clockwise = true,
          letter_spacing = 0,
          baselineshift = 0

This version of the `textcurve()` function is designed for shorter text strings
that need positioning around a circle. (A cheesy effect much beloved of hipster
brands and retronauts.)

`letter_spacing` adjusts the tracking/space between chars, tighter is (-),
looser is (+)).  `baselineshift` moves the text up or down away from the
baseline.

textcurvecentred (UK spelling) is a synonym
"""
function textcurvecentered(the_text, the_angle, the_radius, center::Point;
                           clockwise = true,
                           letter_spacing = 0,
                           baselineshift = 0
      )
    atextbox = textextents(the_text)
    atextwidth = atextbox[3]                         # width of text
    if clockwise
        baselineradius = the_radius + baselineshift  # could be adjusted if we knew font height
    else
        baselineradius = the_radius - baselineshift  # could be adjusted if we knew font height
    end

    # TODO
    # hack to adjust starting angle for the letter spacing
    # to do it accurately would take lots more code and brain cells
    lspaced = length(the_text) * letter_spacing
    lspacedangle = atan(lspaced, baselineradius)

    theta = atextwidth/baselineradius               # find angle
    if clockwise
        starttextangle = the_angle - (theta/2) - lspacedangle/2
    else
        starttextangle = the_angle + (theta/2) + lspacedangle/2
    end
    starttextxpos = baselineradius * cos(starttextangle)
    starttextypos = baselineradius * sin(starttextangle)
    textcurve(the_text, starttextangle, baselineradius, center,
              clockwise=clockwise, letter_spacing=letter_spacing)
end

textcurvecentred = textcurvecentered

# the professional interface

"""
    setfont(family, fontsize)

Select a font and specify the size.

Example:

    setfont("Helvetica", 24)
    settext("Hello in Helvetica 24 using the Pro API", Point(0, 10))
"""
function setfont(family::AbstractString, fontsize)
    # output is set relative to 96dpi
    # so it needs to be rescaled
    fsize = fontsize * 72/96
    set_font_face(get_current_cr(), string(family, " ", fsize))
end

"""
    settext(text, pos;
        halign = "left",
        valign = "bottom",
        angle  = 0, # degrees!
        markup = false)

    settext(text;
        kwargs)

Draw the `text` at `pos` (if omitted defaults to `0/0`). If no font is
specified, on macOS the default font is Times Roman.

To align the text, use `halign`, one of "left", "center", or "right", and
`valign`, one of "top", "center", or "bottom".

`angle` is the rotation - in counterclockwise degrees, rather than Luxor's
default clockwise (+x-axis to +y-axis) radians.

If `markup` is `true`, then the string can contain some HTML-style markup.
Supported tags include:

    , , <s>, , , , <big>, , <tt>, and

The `<span>` tag can contains things like this:

    unreadable text
"""
settext(text::AbstractString, pos::Point; kwargs...) =
    Cairo.text(get_current_cr(), pos.x, pos.y, text; kwargs...)

settext(text; kwargs...) = settext(text, O; kwargs...)

"""
    label(txt::AbstractString, alignment::Symbol=:N, pos::Point=O;
        offset=5,
        leader=false,
        leaderoffsets=[0.0, 1.0])

Add a text label at a point, positioned relative to that point, for example, `:N` signifies
North and places the text directly above that point.

Use one of `:N`, `:S`, `:E`, `:W`, `:NE`, `:SE`, `:SW`, `:NW` to position the label
relative to that point.

    label("text")          # positions text at North (above), relative to the origin
    label("text", :S)      # positions text at South (below), relative to the origin
    label("text", :S, pt)  # positions text South of pt
    label("text", :N, pt, offset=20)  # positions text North of pt, offset by 20

The default offset is 5 units.

If `leader` is true, draw a line as well.

`leaderoffsts` uses normalized fractions (see `between()`)
to specify the gap between the designated points and the
start and end of the lines.

TODO: Negative offsets don't give good results.
"""
function label(txt::AbstractString, alignment::Symbol=:N, pos::Point=O;
        offset=5,
        leader=false,
        leaderoffsets=[0.0, 1.0])
    # alignment one of :N, :S, :E, ;W, ;NE; :SE, :SW, :NW
    if in(alignment, [:N, :n])
        pt = Point(pos.x, pos.y - offset)
        text(txt, pt, halign = :center, valign = :bottom)
    elseif in(alignment, [:E, :e])
        pt = Point(pos.x + offset, pos.y)
        text(txt, pt, halign = :left, valign = :middle)
    elseif in(alignment, [:S, :s])
        pt = Point(pos.x, pos.y + offset)
        text(txt, pt, halign = :center, valign = :top)
    elseif in(alignment, [:W, :w])
        pt = Point(pos.x - offset, pos.y)
        text(txt, pt, halign = :right, valign = :middle)
    elseif in(alignment, [:NE, :ne])
        pt = Point(pos.x + (offset * cos(pi/4)), pos.y - (offset * sin(pi/4)))
        text(txt, pt, halign = :left, valign = :bottom)
    elseif in(alignment, [:SE, :se])
        pt = Point(pos.x + (offset * cos(pi/4)), pos.y + (offset * sin(pi/4)))
        text(txt, pt, halign = :left, valign = :top)
    elseif in(alignment, [:SW, :sw])
        pt = Point(pos.x - (offset * cos(pi/4)), pos.y + (offset * sin(pi/4)))
        text(txt, pt, halign = :right, valign = :top)
    elseif in(alignment, [:NW, :nw])
        pt = Point(pos.x - (offset * cos(pi/4)), pos.y - (offset * sin(pi/4)))
        text(txt, pt, halign = :right, valign = :bottom)
    end
    if leader
        line(between(pos, pt, leaderoffsets[1]), between(pos, pt, leaderoffsets[2]), :stroke)
    end
end

"""
    label(txt::AbstractString, rotation::Float64, pos::Point=O;
        offset=5,
        leader=false,
        leaderoffsets=[0.0, 1.0])

Add a text label at a point, positioned relative to that point, for example,
`0.0` is East, `pi` is West.

    label("text", pi)          # positions text to the left of the origin
"""
function label(txt::AbstractString, rotation::Real, pos::Point=O;
        offset=5,
        leader=false,
        leaderoffsets=[0.0, 1.0])
    if 0 < rotation <= pi/4
        vertalign  = :middle
        horizalign = :left
    elseif pi/4 < rotation <= 3pi/4
        vertalign  = :top
        horizalign = :center
    elseif 3pi/4 < rotation <= 5pi/4
        vertalign  = :middle
        horizalign = :right
    elseif 5pi/4 < rotation <= 7pi/4
        vertalign  = :bottom
        horizalign = :center
    else
        vertalign  = :middle
        horizalign = :left
    end
    pt = pos + polar(offset, rotation)
    if leader
        line(between(pos, pt, leaderoffsets[1]), between(pos, pt, leaderoffsets[2]), :stroke)
    end
    text(txt, pt, halign = horizalign, valign=vertalign)
end

"""
    splittext(s)

Split the text in string `s` into an array, but keep all the separators
attached to the preceding word.
"""
function splittext(s)
    # split text into array, keeping all separators
    # hyphens stay with first word
    result = Array{String, 1}()
    iobuffer = IOBuffer()
    for c in s
        if isspace(c)
            push!(result, String(take!(iobuffer)))
            iobuffer = IOBuffer()
        elseif c == '-' # hyphen splits words but needs keeping
            print(iobuffer, c)
            push!(result, String(take!(iobuffer)))
            iobuffer = IOBuffer()
        else
            print(iobuffer, c)
        end
    end
    push!(result, String(take!(iobuffer)))
    return result
end

"""
    textlines(s::AbstractString, width::Real;
         rightgutter=5)

Split the text in `s` into lines up to `width` units wide (in the current font).

Returns an array of strings. Use `textwrap` to draw an array of strings.

TODO: A `rightgutter` optional keyword adds some padding to the right hand side
of the column. This appears to be needed sometimes -— perhaps the algorithm
needs improving to take account of the interaction of `textextents` and spaces?
"""
function textlines(s::T where T<:AbstractString, width::Real; rightgutter=5)
    result = String[]
    fields = splittext(s)
    textwidth = width - rightgutter
    # how to get the width of a space?
    spacewidth = textextents("m m")[3] - textextents("mm")[3]
    spaceleft = textwidth
    currentline = String[]
    for word in fields
        # hyphenation can leave an empty string field
        if word == ""
            word = " "
            push!(currentline, " ")
        end
        wordextents =  textextents(word)
        widthofword = wordextents[3] + wordextents[6]
        isapprox(widthofword, 0.0, atol=0.1) && continue
        if (widthofword + spacewidth) > spaceleft
            # start a new line
            push!(result, strip(join(currentline)))
            currentline=[]
            spaceleft = textwidth
        end
        if endswith(word, "-")
            push!(currentline, word)
        else
            push!(currentline, word * " ")
        end
        spaceleft -= (widthofword + spacewidth)
    end
    push!(result, strip(join(currentline)))

    # strip trailing spaces
    for i in eachindex(result)
       result[i] = strip(result[i])
    end
    return result
end

"""
    textbox(lines::Array, pos::Point=O;
        leading = 12,
        linefunc::Function = (linenumber, linetext, startpos, height) -> (),
        alignment=:left)

Draw the strings in the array `lines` vertically downwards. `leading` controls
the spacing between each line (default 12), and `alignment` determines the
horizontal alignment (default `:left`).

Optionally, before each line, execute the function
`linefunc(linenumber, linetext, startpos, height)`.

Returns the position of what would have been the next line.

See also `textwrap()`, which modifies the text so that the lines fit into a
specified width.
"""
function textbox(lines::Array, pos::Point=O;
        leading = 0,
        linefunc::Function = (linenumber, linetext, startpos, leading) -> (),
        alignment=:left)

    # find height of first non-empty line
    firstrealline = filter(!isempty, lines)[1]
    te = textextents(firstrealline)
    if leading == 0
        leading = te[4]
    end
    startpos = Point(pos.x, pos.y + leading)
    for (linenumber, linetext) in enumerate(lines)
        linefunc(linenumber, linetext, startpos, leading)
        text(linetext, startpos, halign=alignment)
        startpos = Point(startpos.x, startpos.y + leading)
    end
    return startpos
end
"""
    textbox(s::AbstractString, pos::Point=O;
        leading = 12,
        linefunc::Function = (linenumber, linetext, startpos, height) -> (),
        alignment=:left)
"""
textbox(s::AbstractString, pos::Point=O; kwargs...) = textbox([s], pos; kwargs...)

"""
    textwrap(s::T where T<:AbstractString, width::Real, pos::Point;
        rightgutter=5,
        leading=0)
    textwrap(s::T where T<:AbstractString, width::Real, pos::Point, linefunc::Function;
        rightgutter=5,
        leading=0)

Draw the string in `s` by splitting it at whitespace characters into lines, so that each
line is no longer than `width` units. The text starts at `pos` such that the first line of
text is drawn entirely below a line drawn horizontally through that position. Each line is
aligned on the left side, below `pos`.

See also `textbox()`.

Optionally, before each line, execute the function `linefunc(linenumber, linetext, startpos, leading)`.

If you don't supply a value for `leading`, the font's built-in extents are used.

Text with no whitespace characters won't wrap. You can write a simple chunking function
to split a string or array into chunks:

```
chunk(x, n) = [x[i:min(i+n-1,length(x))] for i in 1:n:length(x)]
```

For example:

```
textwrap(the_text, 300, boxtopleft(BoundingBox()) + 20,
    (ln, lt, sp, ht) -> begin
        c = count(t -> occursin(r"[[:punct:]]", t), split(lt, ""))
        @layer begin
            fontface("Menlo")
            sethue("darkred")
            text(string("[", c, "]"), sp + (310, 0))
        end
    end)
```

puts a count of the number of punctuation characters in each
line at the end of the line.

Returns the position of what would have been the next line.
"""
function textwrap(s::T where T<:AbstractString, width::Real, pos::Point, linefunc::Function;
        rightgutter=5,
        leading=0)
    lines = textlines(s, width; rightgutter=rightgutter)
    textbox(lines, pos, linefunc=linefunc, leading=leading)
end

textwrap(s::T where T<:AbstractString, width::Real, pos::Point; kwargs...) =
    textwrap(s, width, pos, (linenumber, linetext, startpos, height) -> ();
             kwargs...)

"""
    texttrack(txt, pos, tracking, fontsize=12;
        action=:fill,
        halign=:left,
        valign=:baseline,
        startnewpath=true)

Place the text in `txt` at `pos`, left-justified, and letter
space ('track') the text using the value in `tracking`.

The tracking units depend on the current font size! 1 is
1/1000 em. In a 6‑point font, 1 em equals 6 points; in a
10‑point font, 1 em equals 10 points.

A value of -50 would tighten the letter spacing noticeably.
A value of 50 would make the text more open.

The text drawing action applied to each character defaults
to `textoutlines(... :fill)`.
"""
function texttrack(txt, pos, tracking, fontsize=12;
            action=:fill,
            halign=:left,
            valign=:baseline,
            startnewpath=true)
    te = textextents(txt)
    for i in txt
        glyph = string(i)
        glyph_x_bearing, glyph_y_bearing, glyph_width, glyph_height, glyph_x_advance, glyph_y_advance = textextents(glyph)
        textoutlines(glyph, pos, action, halign=halign, valign=valign, startnewpath=startnewpath)
        x = glyph_x_advance + (tracking/1000) * fontsize
        pos += (x, 0)
    end
end

ext: jl


text:
using SparseArrays, LinearAlgebra, BenchmarkTools, MKL

A = rand(200,200)
B = rand(200,200)

C = zeros(200,200)
C[3, 6] = 1
C[4, 6] = 1
C[31, 16] = 1
C[32, 56] = 1
C[43, 46] = 1
C[109, 89] = 1

ind = findall(x -> isone(x), C)

### BLAS
function genmul(A, C, B)
    D = A*C*B
    return D
end

@benchmark genmul(A, C, B)

C = sparse(C)

### SPARSE
@benchmark genmul(A, C, B)

### LOOP, column layout, without boradcasting, with @inbounds
B_2 = permutedims(B)

function outer_mul_4(A, B, ind)
    C = zeros(size(A))
    @inbounds for i in 1:length(ind)
        C += A[:, ind[i][1]]*B[:, ind[i][2]]'
    end
    return C
end

@benchmark outer_mul_4(A,B_2, ind)

### call BLAS directly
function outer_mul_5(A, B, ind)
    C = zeros(size(A))
    for i in 1:length(ind)
        BLAS.ger!(1.0, A[:, ind[i][1]], B[:, ind[i][2]], C)
    end
    return C
end

@benchmark outer_mul_5(A,B_2, ind)

### call BLAS directly, witout column mayor --> WINNER!!!
function outer_mul_7(A, B, ind)
    C = zeros(size(A))
    for i in 1:length(ind)
        BLAS.ger!(1.0, A[:, ind[i][1]], B[ind[i][2], :], C)
    end
    return C
end

@benchmark outer_mul_7(A,B, ind)

### with mul!
C = Matrix(C)

function outer_mul_6(A, B, C)
    D = similar(A)
    E = similar(C)
    mul!(D, A, C)
    mul!(E, D, B)
    return E
end

@benchmark outer_mul_6(A,B, C)

###################################
using SparseArrays, LinearAlgebra, BenchmarkTools, MKL

A = rand(100,100)
m = zeros(100)
m[10] = 1.0
m[56] = 1.0
m[84] = 1.0

A*m

C = similar(m)

ind = findall(x -> isone(x), m)

function sparse_outer_mul!(A, C,ind) #computes A*S*B -> C, where ind gives the entries of S that are 1
    C = zeros(size(A, 1))
    @views @inbounds for i in 1:length(ind)
        C += A[:, ind[i]]
    end
    return C
end

isapprox(A*m, sparse_outer_mul!(A, ind))

@benchmark A*m

@benchmark sparse_outer_mul!(A, ind)

function sparse_outer_mul_2!(A, ind) #computes FB*∇m, where ind gives the entries of S that are 1
    @views C = sum(A[:, ind], dims = 2)
    return C
end

isapprox(A*m, sparse_outer_mul_2!(A, ind))

@benchmark @views sum(A[:, ind], dims = 2)

@benchmark sparse_outer_mul_2!(A, ind)
@benchmark sparse_outer_mul_2!(C, A, ind)

###################################
using SparseArrays, LinearAlgebra, BenchmarkTools, MKL

A = rand(100,120)
S = zeros(120,120)
B = rand(120)

S[3, 6] = 1
S[4, 6] = 1
S[31, 16] = 1
S[32, 56] = 1
S[43, 46] = 1
S[109, 89] = 1

A*S*B

C = Vector{Float64}(undef, 100)

ind = findall(x -> isone(x), S)

function sparse_outer_mul!(C, A, B, ind) #computes A*S*B -> C, where ind gives the entries of S that are 1
    fill!(C, 0.0)
    @views @inbounds for i in 1:length(ind)
        BLAS.axpy!(B[ind[i][2]], A[:, ind[i][1]], C)
    end
end

sparse_outer_mul!(C, A, B, ind)

isapprox(A*S*B, C)

@benchmark genmul($A,$S,$B)

@benchmark sparse_outer_mul!(C, A, B, ind)

function sparse_outer_mul_2!(C, A, B::Vector, ind) #computes A*S*B -> C, where ind gives the entries of S that are 1
    fill!(C, 0.0)
    @views @inbounds for i in 1:length(ind)
        C .+= B[ind[i][2]].*A[:, ind[i][1]]
    end
end

sparse_outer_mul_2!(C, A, B, ind)

isapprox(A*S*B, C)

@benchmark sparse_outer_mul_2!(C, A, B, ind)


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
import StatsBase: fit!, predict
using StatsBase
using BSON

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
        required = true
        arg_type = Int
        help = "seed"
    "dataset"
        required = true
        arg_type = String
        help = "dataset"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################
modelname = "loda_opt"
function fit(data)
	# fit train data
	try
		global model, fit_t, _, _, _ = @timed GenerativeAD.Models.Loda(data[1][1])
	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models

	training_info = (
		fit_t = fit_t,
		n_histograms = length(model.histograms),
		n_bins = [length(x.p) for x in model.histograms],
		model = nothing
		)

	# now return the different scoring functions
	training_info, [(x->model(x), (n_histograms=training_info.n_histograms,))]
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# set a maximum for parameter sampling retries
try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
for seed in 1:max_seed
	savepath = datadir("experiments/tabular_clean_val_default/$(modelname)/$(dataset)/seed=$(seed)")
	mkpath(savepath)

	# get data
	data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

	@info "Trying to fit $modelname on $dataset..."

	# fit
	training_info, results = fit(data)

	# here define what additional info should be saved together with parameters, scores, labels and predict times
	save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))

	# now loop over all anomaly score funs
	for result in results
		GenerativeAD.experiment(result..., data, savepath; save_entries...)
	end
end


text:
using RoME
using Statistics
# , Distributions

using Test

@testset "basic Point2Point2 test" begin

global fg = initfg()

addVariable!(fg, :x0, Point2)
addFactor!(fg, [:x0], PriorPoint2(MvNormal(zeros(2), Matrix{Float64}(LinearAlgebra.I, 2,2))))

addVariable!(fg, :x1, Point2)
addFactor!(fg, [:x0;:x1], Point2Point2(MvNormal([10;0.0], Matrix{Float64}(LinearAlgebra.I, 2,2))))

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)

@test sum( abs.(Statistics.mean(getVal(fg, :x0),dims=2) .- [0.0;0]) .< [1.0;1.0]) == 2
@test sum( abs.(Statistics.mean(getVal(fg, :x1),dims=2) .- [10.0;0]) .< [1.0;1.0]) == 2

end

# drawGraph(fg, show=true)
#
# using RoMEPlotting, Gadfly, Fontconfig, Cairo
#
# drawTree(tree, show=true)
#
# stuff = treeProductUp(fg, tree, :x0, :l1)
# stuff = treeProductUp(fg, tree, :l1, :l1)
#
# plotKDE(kde!(stuff[1]), levels=3)
# # drawLandms(fg) #, regexLandmark=r"x") |> PDF("/tmp/test.pdf")
# plotKDE(fg, ls(fg))

@testset "test Point2Point2Range{T}..." begin

global N=100 # return to 200
global fg = initfg()

addVariable!(fg, :x0, Point2, N=N)
addFactor!(fg, [:x0], PriorPoint2(MvNormal([100.0;0], Matrix{Float64}(LinearAlgebra.I, 2,2))))

addVariable!(fg, :x1, Point2, N=N)
addFactor!(fg, [:x1], PriorPoint2(MvNormal([0.0;100.0], Matrix{Float64}(LinearAlgebra.I, 2,2))))

addVariable!(fg, :l1, Point2, N=N)
addFactor!(fg, [:x0;:l1], Point2Point2Range(Normal(100.0, 1.0)))
addFactor!(fg, [:x1;:l1], Point2Point2Range(Normal(100.0, 1.0)))

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

@test 0.1*N < sum( 90 .< getVal(fg, :l1)[1,:] .< 110 )
@test 0.1*N < sum( -10 .< getVal(fg, :l1)[2,:] .< 10 )

@test 0.1*N < sum( -10 .< getVal(fg, :l1)[1,:] .< 10 )
@test 0.1*N < sum( 90 .< getVal(fg, :l1)[2,:] .< 110 )

global voidsel1 =  10.0 .< getVal(fg, :l1)[1,:]
@test sum( getVal(fg, :l1)[2,voidsel1] .< 70 ) < 0.3*N

global voidsel2 =  10.0 .< getVal(fg, :l1)[2,:]
@test sum( getVal(fg, :l1)[1,voidsel2] .< 70 ) < 0.3*N

@test sum( 120 .< abs.(getVal(fg, :l1)[1,:]) ) < 0.3*N
@test sum( 120 .< abs.(getVal(fg, :l1)[2,:]) ) < 0.3*N

end

#
#
# # debugging
#

# setVal!(fg, :x0, zeros(2,1))
# setVal!(fg, :l1, zeros(2,1))
# #
# using RoMEPlotting, KernelDensityEstimatePlotting
#
# plotKDE(fg, :x0)
#
#
# pts1 = IIF.approxConv(fg, :x0l1f1, :l1, N=N)
#
# plotKDE(KDE.kde!(pts1))
#
#
# pts2 = IIF.approxConv(fg, :x1l1f1, :l1, N=N)
#
# plotKDE([KDE.kde!(pts1); KDE.kde!(pts2)])

# #
# #
# #
#
# #

# using RoMEPlotting, KernelDensityEstimatePlotting
#
# plotKDE(fg, [:x0; :x1; :l1], levels=3)
#
# stuff = IIF.localProduct(fg, :l1)
# stuff
#
# plotKDE(stuff[2], levels=3, c=["red";"cyan"])

##


text:
# Wisdom & Hernandez version of Kepler solver, but with quartic convergence.

using ForwardDiff

include("g3.jl")

#function calc_ds_opt(y::T,yp::T,ypp::T,yppp::T) where {T <: Real}
## Computes quartic Newton's update to equation y=0 using first through 3rd derivatives.
## Uses techniques outlined in Murray & Dermott for Kepler solver.
## Rearrange to reduce number of divisions:
#num = y*yp
#den1 = yp*yp-y*ypp*.5
#den12 = den1*den1
#den2 = yp*den12-num*.5*(ypp*den1-third*num*yppp)
#return -y*den12/den2
#end

function solve_kepler_drift!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,
   r0::T,s0::T,state::Array{T,1},drift_first::Bool) where {T <: Real}
# Solves elliptic Kepler's equation for both elliptic and hyperbolic cases,
# along with a drift before or after the kepler step.
# Initial guess (if s0 = 0):
r0inv = inv(r0)
beta0inv = inv(beta0)
signb = sign(beta0)
if s0 == zero(T)
  s = h*r0inv
else
  s = copy(s0)
end
s0 = copy(s)
sqb = sqrt(signb*beta0)
y = zero(T); yp = one(T)
iter = 0
ds = Inf
zeta = k-r0*beta0
if drift_first
  eta = dot(x0-h*v0,v0)
else
  eta = dot(x0,v0)
end
KEPLER_TOL = sqrt(eps(h))
while iter == 0 || (abs(ds) > KEPLER_TOL && iter < 10)
  xx = sqb*s
  if beta0 > 0
    sx = sin(xx); cx = cos(xx)
  else
    cx = cosh(xx); sx = exp(xx)-cx
  end
  sx *= sqb
# Third derivative:
  yppp = zeta*cx - signb*eta*sx
# First derivative:
  yp = (-yppp+ k)*beta0inv
# Second derivative:
  ypp = signb*zeta*beta0inv*sx + eta*cx
  y  = (-ypp + eta +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
  ds = calc_ds_opt(y,yp,ypp,yppp)
  s += ds
  iter +=1
end
# Since we updated s, need to recompute:
xx = 0.5*sqb*s
if beta0 > 0
  sx = sin(xx); cx = cos(xx)
else
  cx = cosh(xx); sx = exp(xx)-cx
end
# Now, compute final values.  Compute Wisdom/Hernandez G_i^\beta(s) functions:
g1bs = 2.0*sx*cx/sqb
g2bs = 2.0*signb*sx^2*beta0inv
g0bs = 1.0-beta0*g2bs
# This should be computed to prevent roundoff error. [ ]
#g3bs = (1.0-g1bs)*beta0inv
g3bs = G3(s*sqb,beta0)
#if typeof(g1bs) == Float64
#  println("g1: ",g1bs," g2: ",g2bs," g3: ",g3bs)
#end
# Compute Gauss' Kepler functions:
f = one(T) - k*r0inv*g2bs # eqn (25)
g = r0*g1bs + eta*g2bs # eqn (27)
if drift_first
  r = norm(f*(x0-h*v0)+g*v0)
else
  r = norm(f*x0+g*v0)
end
#if typeof(r) == Float64
#  println("r: ",r)
#end
rinv = inv(r)
dfdt = -k*g1bs*rinv*r0inv
if drift_first
  # Drift backwards before Kepler step: (1/22/2018)
  fm1 = -k*r0inv*g2bs
  # This is given in 2/7/2018 notes:
  gmh = k*r0inv*(r0*(g1bs*g2bs-g3bs)+eta*g2bs^2+k*g3bs*g2bs)
else
  # Drift backwards after Kepler step: (1/24/2018)
  fm1 =  k*rinv*(g2bs-k*r0inv*H1(s*sqb,beta0))
  # This is g-h*dgdt
  gmh = k*rinv*(r0*H2(s*sqb,beta0)+eta*H1(s*sqb,beta0))
end
# Compute velocity component functions:
if drift_first
  # 2/1/18 notes:
  dgdtm1 = k*r0inv*rinv*(r0*g0bs*g2bs+eta*g1bs*g2bs+k*g1bs*g3bs)
else
  # 1/22/18 notes:
  dgdtm1 = -k*rinv*g2bs
end
#if typeof(fm1) == Float64
#  println("fm1: ",fm1," dfdt: ",dfdt," gmh: ",gmh," dgdt-1: ",dgdtm1)
#end
for j=1:3
# Compute difference vectors (finish - start) of step:
  state[1+j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
  state[4+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
end
return s,f,g,dfdt,dgdtm1,cx,sx,g1bs,g2bs,g3bs,r,rinv,ds,iter
end

function jac_delxv(x0::Array{T,1},v0::Array{T,1},k::T,s::T,beta0::T,h::T,drift_first::Bool) where {T <: Real}
# Using autodiff, computes Jacobian of delx & delv with respect to x0, v0, k, s, beta0 & h.

# Autodiff requires a single-vector input, so create an array to hold the independent variables:
  input = zeros(T,10)
  input[1:3]=x0; input[4:6]=v0; input[7]=k; input[8]=s; input[9]=beta0; input[10]=h

# Create a closure so that the function knows value of drift_first:

  function delx_delv(input::Array{T2,1}) where {T2 <: Real} # input = x0,v0,k,s,beta0,h,drift_first
  # Compute delx and delv from h, s, k, beta0, x0 and v0:
  x0 = input[1:3]; v0 = input[4:6]; k = input[7]
  s = input[8]; beta0=input[9]; h = input[10]
  # Set up a single output array for delx and delv:
  delxv = zeros(T2,6)
  # Compute square root of beta0:
  signb = sign(beta0)
  sqb = sqrt(abs(beta0))
  beta0inv=inv(beta0)
  if drift_first
    r0 = norm(x0-h*v0)
    eta = dot(x0-h*v0,v0)
  else
    r0 = norm(x0)
    eta = dot(x0,v0)
  end
  r0inv = inv(r0)
  # Since we updated s, need to recompute:
  xx = 0.5*sqb*s
  if beta0 > 0
    sx = sin(xx); cx = cos(xx)
  else
    cx = cosh(xx); sx = exp(xx)-cx
  end
  # Now, compute final values.  Compute Wisdom/Hernandez G_i^\beta(s) functions:
  g1bs = 2.0*sx*cx/sqb
  g2bs = 2.0*signb*sx^2*beta0inv
  g0bs = 1.0-beta0*g2bs
  g3bs = G3(s*sqb,beta0)
  # Compute Gauss' Kepler functions:
  f = 1.0 - k*r0inv*g2bs # eqn (25)
  g = r0*g1bs + eta*g2bs # eqn (27)
  if drift_first
    r = norm(f*(x0-h*v0)+g*v0)
  else
    r = norm(f*x0+g*v0)
  end
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv
  if drift_first
    # Drift backwards before Kepler step: (1/22/2018)
    fm1 = -k*r0inv*g2bs
    # This is given in 2/7/2018 notes:
    gmh = k*r0inv*(r0*(g1bs*g2bs-g3bs)+eta*g2bs^2+k*g3bs*g2bs)
  else
    # Drift backwards after Kepler step: (1/24/2018)
    fm1 =  k*rinv*(g2bs-k*r0inv*H1(s*sqb,beta0))
    # This is g-h*dgdt
    gmh = k*rinv*(r0*H2(s*sqb,beta0)+eta*H1(s*sqb,beta0))
  end
  # Compute velocity component functions:
  if drift_first
    dgdtm1 = k*r0inv*rinv*(r0*g0bs*g2bs+eta*g1bs*g2bs+k*g1bs*g3bs)
  else
    dgdtm1 = -k*rinv*g2bs
  end
  for j=1:3
  # Compute difference vectors (finish - start) of step:
    delxv[  j] = fm1*x0[j]+gmh*v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
    delxv[3+j] = dfdt*x0[j]+dgdtm1*v0[j]    # velocity v_ij(t+h)-v_ij(t)
  end
  return delxv
  end

# Use autodiff to compute Jacobian:
delxv_jac = ForwardDiff.jacobian(delx_delv,input)
# Return Jacobian:
return  delxv_jac
end

function kep_drift_ell_hyp!(x0::Array{T,1},v0::Array{T,1},k::T,h::T,
  s0::T,state::Array{T,1},drift_first::Bool) where {T <: Real}
if drift_first
  r0 = norm(x0-h*v0)
else
  r0 = norm(x0)
end
beta0 = 2*k/r0-dot(v0,v0)
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
# Now, solve for s in elliptical Kepler case:
f = zero(T); g=zero(T); dfdt=zero(T); dgdtm1=zero(T); cx=zero(T);sx=zero(T);g1bs=zero(T);g2bs=zero(T);g3bs=zero(T)
s=zero(T); ds = zero(T); r = zero(T);rinv=zero(T); iter=0
if beta0 > zero(T) || beta0 < zero(T)
   s,f,g,dfdt,dgdtm1,cx,sx,g1bs,g2bs,g3bs,r,rinv,ds,iter = solve_kepler_drift!(h,k,x0,v0,beta0,r0,
    s0,state,drift_first)
else
#  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero(T); fill!(state,zero(T)); rinv=zero(T); s=zero(T); ds=zero(T); iter = 0
end
state[8]= r
# These need to be updated. [ ]
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# recompute beta:
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

function kep_drift_ell_hyp!(x0::Array{T,1},v0::Array{T,1},k::T,h::T,
  s0::T,state::Array{T,1},jacobian::Array{T,2},drift_first::Bool,d::Derivatives{T}) where {T <: Real}
if drift_first
  r0 = norm(x0-h*v0)
else
  r0 = norm(x0)
end
beta0 = 2*k/r0-dot(v0,v0)
# Computes the Jacobian as well
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.
r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in elliptical Kepler case:
f = zero(T); g=zero(T); dfdt=zero(T); dgdtm1=zero(T); cx=zero(T);sx=zero(T);g1bs=zero(T);g2bs=zero(T);g3bs=zero(T)
s=zero(T); ds=zero(T); r = zero(T);rinv=zero(T); iter=0
if beta0 > zero(T) || beta0 < zero(T)
   s,f,g,dfdt,dgdtm1,cx,sx,g1bs,g2bs,g3bs,r,rinv,ds,iter = solve_kepler_drift!(h,k,x0,v0,beta0,r0,
    s0,state,drift_first)
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v} & q0 = {x0,v0}.
  delxv_jac = jac_delxv(x0,v0,k,s,beta0,h,drift_first)
#  println("computed jacobian with autodiff")
  # Add in partial derivatives with respect to x0, v0 and k:
  jacobian[1:6,1:7] = delxv_jac[1:6,1:7]
# Add in s and beta0 derivatives:
  compute_jacobian_kep_drift!(h,k,x0,v0,beta0,s,f,g,dfdt,dgdtm1,cx,sx,g1bs,g2bs,r0,r,jacobian,delxv_jac,drift_first,d)
else
#  println("Not elliptic or hyperbolic ",beta0," x0 ",x0)
  r= zero(T); fill!(state,zero(T)); rinv=zero(T); s=zero(T); ds=zero(T); iter = 0
end
# recompute beta:
state[8]= r
# These need to be updated. [ ]
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

function compute_jacobian_kep_drift!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,s::T,
  f::T,g::T,dfdt::T,dgdtm1::T,cx::T,sx::T,g1::T,g2::T,r0::T,r::T,jacobian::Array{T,2},drift_first::Bool,d::Derivatives{T}) where {T <: Real}
# This needs to be updated to incorporate backwards drifts. [ ]
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v,k} & q0 = {x0,v0,k}.
# Now, compute the Jacobian: (9/18/2017 notes)
g0 = one(T)-beta0*g2
# Expand g3 as a series if s is small:
sqb = sqrt(abs(beta0))
g3bs = G3(s*sqb,beta0)
#  g3bs = (s-g1)/beta0
if drift_first
  eta = dot(x0-h*v0,v0)
else
  eta = dot(x0,v0)
end
absv0 = norm(v0)
#dsdbeta = (2h-r0*(s*g0+g1)+k/beta0*(s*g0-g1)-eta*s*g1)/(2beta0*r)
# New expression derived on 8/14/2019:
dsdbeta = (h+eta*g2+2*k*g3bs-s*r)/(2beta0*r)
dsdr0 = -(2k/r0^2*dsdbeta+g1/r)
dsda0 = -g2/r
dsdv0 = -2absv0*dsdbeta
dsdk = 2/r0*dsdbeta-g3bs/r
dbetadr0 = -2k/r0^2
dbetadv0 = -2absv0
dbetadk  = 2/r0
# "p" for partial derivative:
for i=1:3
  d.pxpr0[i] = k/r0^2*g2*x0[i]+g1*v0[i]
  d.pxpa0[i] = g2*v0[i]
  d.pxpk[i]  = -g2/r0*x0[i]
  d.pxps[i]  = -k/r0*g1*x0[i]+(r0*g0+eta*g1)*v0[i]
  d.pxpbeta[i] = -k/(2beta0*r0)*(s*g1-2g2)*x0[i]+1/(2beta0)*(s*r0*g0-r0*g1+s*eta*g1-2*eta*g2)*v0[i]
  d.prvpr0[i] = k*g1/r0^2*x0[i]+g0*v0[i]
  d.prvpa0[i] = g1*v0[i]
  d.prvpk[i] = -g1/r0*x0[i]
  d.prvps[i] = -k*g0/r0*x0[i]+(-beta0*r0*g1+eta*g0)*v0[i]
  d.prvpbeta[i] = -k/(2beta0*r0)*(s*g0-g1)*x0[i]+1/(2beta0)*(-s*r0*beta0*g1+eta*s*g0-eta*g1)*v0[i]
end
prpr0 = g0
prpa0 = g1
prpk  = g2
prps = (k-beta0*r0)*g1+eta*g0
prpbeta = 1/(2beta0)*(s*(k-beta0*r0)*g1+eta*s*g0-eta*g1-2k*g2)
for i=1:3
  d.dxdr0[i] = d.pxps[i]*dsdr0 + d.pxpbeta[i]*dbetadr0 + d.pxpr0[i]
  d.dxda0[i] = d.pxps[i]*dsda0 + d.pxpa0[i]
  d.dxdv0[i] = d.pxps[i]*dsdv0 + d.pxpbeta[i]*dbetadv0
  d.dxdk[i]  = d.pxps[i]*dsdk  + d.pxpbeta[i]*dbetadk + d.pxpk[i]
  d.drvdr0[i] = d.prvps[i]*dsdr0 + d.prvpbeta[i]*dbetadr0 + d.prvpr0[i]
  d.drvda0[i] = d.prvps[i]*dsda0 + d.prvpa0[i]
  d.drvdv0[i] = d.prvps[i]*dsdv0 + d.prvpbeta[i]*dbetadv0
  d.drvdk[i]  = d.prvps[i]*dsdk  + d.prvpbeta[i]*dbetadk + d.prvpk[i]
end
drdr0 = prpr0 + prps*dsdr0 + prpbeta*dbetadr0
drda0 = prpa0 + prps*dsda0
drdv0 = prps*dsdv0 + prpbeta*dbetadv0
drdk  = prpk + prps*dsdk + prpbeta*dbetadk
for i=1:3
  d.vtmp[i] = dfdt*x0[i]+(dgdtm1+1.0)*v0[i]
  d.dvdr0[i] = (d.drvdr0[i]-drdr0*d.vtmp[i])/r
  d.dvda0[i] = (d.drvda0[i]-drda0*d.vtmp[i])/r
  d.dvdv0[i] = (d.drvdv0[i]-drdv0*d.vtmp[i])/r
  d.dvdk[i]  = (d.drvdk[i] -drdk *d.vtmp[i])/r
end
# Now, compute Jacobian:
for i=1:3
  jacobian[  i,  i] = f
  jacobian[  i,3+i] = g
  jacobian[3+i,  i] = dfdt
  jacobian[3+i,3+i] = dgdt+1.0
  jacobian[  i,7] = d.dxdk[i]
  jacobian[3+i,7] = d.dvdk[i]
end
for j=1:3
  for i=1:3
    jacobian[  i,  j] += d.dxdr0[i]*x0[j]/r0 + d.dxda0[i]*v0[j]
    jacobian[  i,3+j] += d.dxdv0[i]*v0[j]/absv0 + d.dxda0[i]*x0[j]
    jacobian[3+i,  j] += d.dvdr0[i]*x0[j]/r0 + d.dvda0[i]*v0[j]
    jacobian[3+i,3+j] += d.dvdv0[i]*v0[j]/absv0 + d.dvda0[i]*x0[j]
  end
end
jacobian[7,7]=one(T)
return
end

function compute_jacobian_kep_drift!(h::T,k::T,x0::Array{T,1},v0::Array{T,1},beta0::T,s::T,
  f::T,g::T,dfdt::T,dgdtm1::T,cx::T,sx::T,g1::T,g2::T,r0::T,r::T,jacobian::Array{T,2},
  delxv_jac::Array{T,2},drift_first::Bool,d::Derivatives{T}) where {T <: Real}
# This needs to be updated to incorporate backwards drifts. [ ]
# Compute the Jacobian.  jacobian[i,j] is derivative of final state variable q[i]
# with respect to initial state variable q0[j], where q = {x,v,k} & q0 = {x0,v0,k}.
# Now, compute the Jacobian: (9/18/2017 notes)
g0 = one(T)-beta0*g2
# Expand g3 as a series if s is small:
sqb = sqrt(abs(beta0))
g3bs = G3(s*sqb,beta0)
#  g3bs = (s-g1)/beta0
if drift_first
  eta = dot(x0-h*v0,v0)
  r0 = norm(x0-h*v0)
else
  eta = dot(x0,v0)
  r0 = norm(x0)
end
absv0 = norm(v0)
r0inv = inv(r0)
#dsdbeta = (2h-r0*(s*g0+g1)+k/beta0*(s*g0-g1)-eta*s*g1)/(2beta0*r)
# New expression derived on 8/14/2019:
dsdbeta = (h+eta*g2+2*k*g3bs-s*r)/(2beta0*r)
dsdr0 = -(2k*r0inv^2*dsdbeta+g1/r)
dsda0 = -g2/r
dsdv0 = -2absv0*dsdbeta
dsdk = 2*r0inv*dsdbeta-g3bs/r
dbetadr0 = -2k*r0inv^2
dbetadv0 = -2absv0
dbetadk  = 2r0inv
# Compute s & beta components of x0 & v0 derivatives:
for i=1:3
  d.dxdr0[i] = delxv_jac[  i,8]*dsdr0 + delxv_jac[  i,9]*dbetadr0
  d.dxda0[i] = delxv_jac[  i,8]*dsda0
  d.dxdv0[i] = delxv_jac[  i,8]*dsdv0 + delxv_jac[  i,9]*dbetadv0
  d.dxdk[i]  = delxv_jac[  i,8]*dsdk  + delxv_jac[  i,9]*dbetadk
  d.dvdr0[i] = delxv_jac[3+i,8]*dsdr0 + delxv_jac[3+i,9]*dbetadr0
  d.dvda0[i] = delxv_jac[3+i,8]*dsda0
  d.dvdv0[i] = delxv_jac[3+i,8]*dsdv0 + delxv_jac[3+i,9]*dbetadv0
  d.dvdk[i]  = delxv_jac[3+i,8]*dsdk  + delxv_jac[3+i,9]*dbetadk
end
# Now, compute Jacobian:
# Add in s & beta components of k derivative:
for i=1:3
  jacobian[  i,7] += d.dxdk[i]
  jacobian[3+i,7] += d.dvdk[i]
end
# Add in s & beta components of x0 & v0 derivatives:
for j=1:3
  for i=1:3
    if drift_first
      jacobian[  i,  j] +=    d.dxdr0[i]*(x0[j]-h*v0[j])*r0inv + d.dxda0[i]*v0[j]
      jacobian[  i,3+j] += -h*d.dxdr0[i]*(x0[j]-h*v0[j])*r0inv + d.dxdv0[i]*v0[j]/absv0 + d.dxda0[i]*(x0[j]-2*h*v0[j])
      jacobian[3+i,  j] +=    d.dvdr0[i]*(x0[j]-h*v0[j])*r0inv + d.dvda0[i]*v0[j]
      jacobian[3+i,3+j] += -h*d.dvdr0[i]*(x0[j]-h*v0[j])*r0inv + d.dvdv0[i]*v0[j]/absv0 + d.dvda0[i]*(x0[j]-2*h*v0[j])
    else
      jacobian[  i,  j] += d.dxdr0[i]*x0[j]*r0inv + d.dxda0[i]*v0[j]
      jacobian[  i,3+j] += d.dxdv0[i]*v0[j]/absv0 + d.dxda0[i]*x0[j]
      jacobian[3+i,  j] += d.dvdr0[i]*x0[j]*r0inv + d.dvda0[i]*v0[j]
      jacobian[3+i,3+j] += d.dvdv0[i]*v0[j]/absv0 + d.dvda0[i]*x0[j]
    end
  end
end
# Mass doesn't change:
jacobian[7,7]=one(T)
return
end


text:
using Oceananigans.Operators: Δxᶠᵃᵃ, Δxᶜᵃᵃ, Δxᶠᶠᵃ, Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δxᶜᶜᵃ
using Oceananigans.Operators: Δyᵃᶠᵃ, Δyᵃᶜᵃ, Δyᶠᶠᵃ, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Δyᶜᶜᵃ

using Oceananigans.Operators: Δzᵃᵃᶜ, Δzᵃᵃᶠ

function test_function_differentiation(T=Float64)
    grid = RectilinearGrid(CPU(), T; size=(3, 3, 3), extent=(3, 3, 3))
    ϕ = rand(T, 3, 3, 3)
    ϕ² = ϕ.^2

    ∂x_ϕ_f = ϕ²[2, 2, 2] - ϕ²[1, 2, 2]
    ∂x_ϕ_c = ϕ²[3, 2, 2] - ϕ²[2, 2, 2]

    ∂y_ϕ_f = ϕ²[2, 2, 2] - ϕ²[2, 1, 2]
    ∂y_ϕ_c = ϕ²[2, 3, 2] - ϕ²[2, 2, 2]

    ∂z_ϕ_f = ϕ²[2, 2, 2] - ϕ²[2, 2, 1]
    ∂z_ϕ_c = ϕ²[2, 2, 3] - ϕ²[2, 2, 2]

    f(i, j, k, grid, ϕ) = ϕ[i, j, k]^2

    assess = true

    for ∂x in (∂xᶜᶜᶜ, ∂xᶜᶜᶠ, ∂xᶜᶠᶜ, ∂xᶜᶠᶠ)
        assess = assess && ∂x(2, 2, 2, grid, f, ϕ) == ∂x_ϕ_c
    end
    for ∂x in (∂xᶠᶜᶜ, ∂xᶠᶜᶠ, ∂xᶠᶠᶜ, ∂xᶠᶠᶠ)
        assess = assess && ∂x(2, 2, 2, grid, f, ϕ) == ∂x_ϕ_f
    end

    for ∂y in (∂yᶜᶜᶜ, ∂yᶜᶜᶠ, ∂yᶠᶜᶜ, ∂yᶠᶜᶠ)
        assess = assess && ∂y(2, 2, 2, grid, f, ϕ) == ∂y_ϕ_c
    end
    for ∂y in (∂yᶜᶠᶜ, ∂yᶠᶠᶜ, ∂yᶜᶠᶠ, ∂yᶠᶠᶠ)
        assess = assess && ∂y(2, 2, 2, grid, f, ϕ) == ∂y_ϕ_f
    end

    for ∂z in (∂zᶜᶜᶜ, ∂zᶜᶠᶜ, ∂zᶠᶜᶜ, ∂zᶠᶠᶜ)
        assess = assess && ∂z(2, 2, 2, grid, f, ϕ) == ∂z_ϕ_c
    end
    for ∂z in (∂zᶜᶜᶠ, ∂zᶜᶠᶠ, ∂zᶠᶜᶠ, ∂zᶠᶠᶠ)
        assess = assess && ∂z(2, 2, 2, grid, f, ϕ) == ∂z_ϕ_f
    end

    stretched_f = [0, 1, 3, 6]
    stretched_c = OffsetArray([-0.5, 0.5, 2, 4.5, 7.5], -1)

    dc(i) = stretched_f[i+1] - stretched_f[i]
    df(i) = stretched_c[i] - stretched_c[i-1]

    grid  = RectilinearGrid(CPU(), T; size=(3, 3, 3), x=stretched_f, y=stretched_f, z=stretched_f, topology = (Bounded, Bounded, Bounded))

    ∂x_f(i, j, k) = (ϕ²[i, j, k]   - ϕ²[i-1, j, k]) / df(i)
    ∂x_c(i, j, k) = (ϕ²[i+1, j, k] - ϕ²[i, j, k])   / dc(i)
    ∂y_f(i, j, k) = (ϕ²[i, j, k]   - ϕ²[i, j-1, k]) / df(j)
    ∂y_c(i, j, k) = (ϕ²[i, j+1, k] - ϕ²[i, j, k])   / dc(j)
    ∂z_f(i, j, k) = (ϕ²[i, j, k]   - ϕ²[i, j, k-1]) / df(k)
    ∂z_c(i, j, k) = (ϕ²[i, j, k+1] - ϕ²[i, j, k])   / dc(k)

    for ∂x in (∂xᶜᶜᶜ, ∂xᶜᶜᶠ, ∂xᶜᶠᶜ, ∂xᶜᶠᶠ)
        assess = assess && ∂x(2, 2, 2, grid, f, ϕ) == ∂x_c(2, 2, 2)
    end
    for ∂x in (∂xᶠᶜᶜ, ∂xᶠᶜᶠ, ∂xᶠᶠᶜ, ∂xᶠᶠᶠ)
        assess = assess && ∂x(2, 2, 2, grid, f, ϕ) == ∂x_f(2, 2, 2)
    end

    for ∂y in (∂yᶜᶜᶜ, ∂yᶜᶜᶠ, ∂yᶠᶜᶜ, ∂yᶠᶜᶠ)
        assess = assess && ∂y(2, 2, 2, grid, f, ϕ) == ∂y_c(2, 2, 2)
    end
    for ∂y in (∂yᶜᶠᶜ, ∂yᶠᶠᶜ, ∂yᶜᶠᶠ, ∂yᶠᶠᶠ)
        assess = assess && ∂y(2, 2, 2, grid, f, ϕ) == ∂y_f(2, 2, 2)
    end

    for ∂z in (∂zᶜᶜᶜ, ∂zᶜᶠᶜ, ∂zᶠᶜᶜ, ∂zᶠᶠᶜ)
        assess = assess && ∂z(2, 2, 2, grid, f, ϕ) == ∂z_c(2, 2, 2)
    end
    for ∂z in (∂zᶜᶜᶠ, ∂zᶜᶠᶠ, ∂zᶠᶜᶠ, ∂zᶠᶠᶠ)
        assess = assess && ∂z(2, 2, 2, grid, f, ϕ) == ∂z_f(2, 2, 2)
    end

    return assess
end

function test_function_interpolation(T=Float64)
    grid = RectilinearGrid(CPU(), T; size=(3, 3, 3), extent=(3, 3, 3))
    ϕ = rand(T, 3, 3, 3)
    ϕ² = ϕ.^2

    ℑx_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[1, 2, 2]) / 2
    ℑx_ϕ_c = (ϕ²[3, 2, 2] + ϕ²[2, 2, 2]) / 2

    ℑy_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[2, 1, 2]) / 2
    ℑy_ϕ_c = (ϕ²[2, 3, 2] + ϕ²[2, 2, 2]) / 2

    ℑz_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[2, 2, 1]) / 2
    ℑz_ϕ_c = (ϕ²[2, 2, 3] + ϕ²[2, 2, 2]) / 2

    f(i, j, k, grid, ϕ) = ϕ[i, j, k]^2

    return (
        ℑxᶜᵃᵃ(2, 2, 2, grid, f, ϕ) == ℑx_ϕ_c &&
        ℑxᶠᵃᵃ(2, 2, 2, grid, f, ϕ) == ℑx_ϕ_f &&

        ℑyᵃᶜᵃ(2, 2, 2, grid, f, ϕ) == ℑy_ϕ_c &&
        ℑyᵃᶠᵃ(2, 2, 2, grid, f, ϕ) == ℑy_ϕ_f &&

        ℑzᵃᵃᶜ(2, 2, 2, grid, f, ϕ) == ℑz_ϕ_c &&
        ℑzᵃᵃᶠ(2, 2, 2, grid, f, ϕ) == ℑz_ϕ_f
    )
end

@testset "Operators" begin
    @info "Testing operators..."

    @testset "Grid lengths, areas, and volume operators" begin
        @info "  Testing grid lengths, areas, and volume operators..."

        x_spacings = ( [eval(Symbol(:Δx, LX, :ᵃ, :ᵃ)) for LX in (:ᶜ, :ᶠ)]...,
                       [eval(Symbol(:Δx, LX, LY, :ᵃ)) for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ)]...,
                       [eval(Symbol(:Δx, LX, LY, LZ)) for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)]...)

        y_spacings = ( [eval(Symbol(:Δy, :ᵃ, LY, :ᵃ)) for LY in (:ᶜ, :ᶠ)]...,
                       [eval(Symbol(:Δy, LX, LY, :ᵃ)) for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ)]...,
                       [eval(Symbol(:Δy, LX, LY, LZ)) for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)]...)

        z_spacings = ( [eval(Symbol(:Δz, :ᵃ, :ᵃ, LZ)) for LZ in (:ᶜ, :ᶠ)]...,
                       [eval(Symbol(:Δz, LX, LY, LZ)) for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)]...)

        FT = Float64
        grid = RectilinearGrid(CPU(), FT, size=(1, 1, 1), extent=(π, 2π, 3π))

        @testset "Easterly lengths" begin
            @info "    Testing easterly lengths..."
            for δ in x_spacings
                @test δ(1, 1, 1, grid) == FT(π)
            end
        end

        @testset "Westerly lengths" begin
            @info "    Testing westerly lengths..."
            for δ in y_spacings
                @test δ(1, 1, 1, grid) == FT(2π)
            end
        end

        @testset "Vertical lengths" begin
            @info "    Testing vertical lengths..."
            for δ in z_spacings
                @test δ(1, 1, 1, grid) == FT(3π)
            end
        end

        @testset "East-normal areas in the yz-plane" begin
            @info "    Testing areas with easterly normal in the yz-plane..."
            for A in (Axᶜᶜᶜ, Axᶠᶜᶜ, Axᶜᶠᶜ, Axᶜᶜᶠ, Axᶠᶠᶠ, Axᶠᶠᶜ, Axᶠᶜᶠ, Axᶜᶠᶠ)
                @test A(1, 1, 1, grid) == FT(6 * π^2)
            end
        end

        @testset "West-normal areas in the xz-plane" begin
            @info "    Testing areas with westerly normal in the xz-plane..."
            for A in (Ayᶜᶜᶜ, Ayᶠᶜᶜ, Ayᶜᶠᶜ, Ayᶜᶜᶠ, Ayᶠᶠᶠ, Ayᶠᶠᶜ, Ayᶠᶜᶠ, Ayᶜᶠᶠ)
                @test A(1, 1, 1, grid) == FT(3 * π^2)
            end
        end

        @testset "Horizontal areas in the xy-plane" begin
            @info "    Testing horizontal areas in the xy-plane..."
            for A in (Azᶜᶜᶜ, Azᶠᶜᶜ, Azᶜᶠᶜ, Azᶜᶜᶠ, Azᶠᶠᶠ, Azᶠᶠᶜ, Azᶠᶜᶠ, Azᶜᶠᶠ)
                @test A(1, 1, 1, grid) == FT(2 * π^2)
            end
        end

        @testset "Volumes" begin
            @info "    Testing volumes..."
            for V in (Vᶜᶜᶜ, Vᶠᶜᶜ, Vᶜᶠᶜ, Vᶜᶜᶠ, Vᶠᶠᶠ, Vᶠᶠᶜ, Vᶠᶜᶠ, Vᶜᶠᶠ)
                @test V(1, 1, 1, grid) == FT(6 * π^3)
            end
        end

    end

    @testset "Function differentiation" begin
        @info "  Testing function differentiation..."
        @test test_function_differentiation()
    end

    @testset "Function interpolation" begin
        @info "  Testing function interpolation..."
        @test test_function_interpolation()
    end

    @testset "2D operators" begin
        @info "  Testing 2D operators..."

        Nx, Ny, Nz = 32, 16, 8
        Lx, Ly, Lz = 100, 100, 100

        arch = CPU()
        grid = RectilinearGrid(CPU(), size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
        bcs = FieldBoundaryConditions(grid, (Center, Center, Center))

        Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz
        Tx, Ty, Tz = Nx+2Hx, Ny+2Hy, Nz+2Hz

        A3 = OffsetArray(zeros(Tx, Ty, Tz), 1-Hx:Nx+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz)
        @. @views A3[1:Nx, 1:Ny, 1:Nz] = rand()
        fill_halo_regions!(A3, bcs, arch, grid)

        # A yz-slice with Nx==1.
        A2yz = OffsetArray(zeros(1+2Hx, Ty, Tz), 1-Hx:1+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz)
        grid_yz = RectilinearGrid(CPU(), size=(1, Ny, Nz), extent=(Lx, Ly, Lz))

        # Manually fill in halos for the slice.
        A2yz[0:2, 0:Ny+1, 1:Nz] .= A3[1:1, 0:Ny+1, 1:Nz]
        A2yz[:, :, 0] .= A2yz[:, :, 1]
        A2yz[:, :, Nz+1] .= A2yz[:, :, Nz]

        # An xz-slice with Ny==1.
        A2xz = OffsetArray(zeros(Tx, 1+2Hy, Tz), 1-Hx:Nx+Hx, 1-Hy:1+Hy, 1-Hz:Nz+Hz)
        grid_xz = RectilinearGrid(CPU(), size=(Nx, 1, Nz), extent=(Lx, Ly, Lz))

        # Manually fill in halos for the slice.
        A2xz[0:Nx+1, 0:2, 1:Nz] .= A3[0:Nx+1, 1:1, 1:Nz]
        A2xz[:, :, 0] .= A2xz[:, :, 1]
        A2xz[:, :, Nz+1] .= A2xz[:, :, Nz]

        test_indices_3d = [(4, 5, 5), (21, 11, 4), (16, 8, 4),  (30, 12, 3), (11, 3, 6), # Interior
                           (2, 10, 4), (31, 5, 6), (10, 2, 4), (17, 15, 5), (17, 10, 2), (23, 5, 7),  # Borderlands
                           (1, 5, 5), (32, 10, 3), (16, 1, 4), (16, 16, 4), (16, 8, 1), (16, 8, 8),  # Edges
                           (1, 1, 1), (32, 16, 8)] # Corners

        test_indices_2d_yz = [(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2),
                              (1, 1, 5), (1, 5, 1), (1, 5, 5), (1, 11, 4),
                              (1, 15, 7), (1, 15, 8), (1, 16, 7), (1, 16, 8)]

        test_indices_2d_xz = [(1, 1, 1), (1, 1, 2), (2, 1, 1), (2, 1, 2),
                              (1, 1, 5), (5, 1, 1), (5, 1, 5), (17, 1, 4),
                              (31, 1, 7), (31, 1, 8), (32, 1, 7), (32, 1, 8)]

        for idx in test_indices_2d_yz
            @test δxᶜᵃᵃ(idx..., grid_yz, A2yz) ≈ 0
            @test δxᶠᵃᵃ(idx..., grid_yz, A2yz) ≈ 0
            @test δyᵃᶜᵃ(idx..., grid_yz, A2yz) ≈ δyᵃᶜᵃ(idx..., grid_yz, A3)
            @test δyᵃᶠᵃ(idx..., grid_yz, A2yz) ≈ δyᵃᶠᵃ(idx..., grid_yz, A3)
            @test δzᵃᵃᶜ(idx..., grid_yz, A2yz) ≈ δzᵃᵃᶜ(idx..., grid_yz, A3)
            @test δzᵃᵃᶠ(idx..., grid_yz, A2yz) ≈ δzᵃᵃᶠ(idx..., grid_yz, A3)
        end

        for idx in test_indices_2d_xz
            @test δxᶜᵃᵃ(idx..., grid_xz, A2xz) ≈ δxᶜᵃᵃ(idx..., grid_xz, A3)
            @test δxᶠᵃᵃ(idx..., grid_xz, A2xz) ≈ δxᶠᵃᵃ(idx..., grid_xz, A3)
            @test δyᵃᶜᵃ(idx..., grid_xz, A2xz) ≈ 0
            @test δyᵃᶠᵃ(idx..., grid_xz, A2xz) ≈ 0
            @test δzᵃᵃᶜ(idx..., grid_xz, A2xz) ≈ δzᵃᵃᶜ(idx..., grid_xz, A3)
            @test δzᵃᵃᶠ(idx..., grid_xz, A2xz) ≈ δzᵃᵃᶠ(idx..., grid_xz, A3)
        end
    end
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for shapes.jl

@testset "Sphere" begin
    @testset "from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Sphere()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(0f0, 0f0, 1f0), NORMAL_Z_false, Vec2D(0f0, 0f0), 1f0, ray, Material())
    end

    @testset "from_behind" begin
        ray = Ray(Point(3f0, 0f0, 0f0), -VEC_X)
        shape = Sphere()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(1f0, 0f0, 0f0), NORMAL_X_false, Vec2D(0f0, 0.5f0), 2f0, ray, Material())
    end

    @testset "from_within" begin
        ray = Ray(Point(0f0, 0f0, 0f0), VEC_X)
        shape = Sphere()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(1f0, 0f0, 0f0), -NORMAL_X_false, Vec2D(0.0, 0.5), 1f0, ray, Material())
    end

    @testset "transposed_from_above" begin
        ray = Ray(Point(10f0, 0f0, 2f0), -VEC_Z)
        shape = Sphere(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(10f0, 0f0, 1f0), NORMAL_Z_false, Vec2D(0f0, 0f0), 1f0, ray, Material())
    end

    @testset "transposed_from_behind" begin
        ray = Ray(Point(13f0, 0f0, 0f0), -VEC_X)
        shape = Sphere(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(11f0, 0f0, 0f0), NORMAL_X_false, Vec2D(0f0, 0.5f0), 2f0, ray, Material())
    end

    @testset "miss_transposed_from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Sphere(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end

    @testset "miss_transposed_from_behind" begin
        ray = Ray(Point(-10f0, 0f0, 0f0), -VEC_X)
        shape = Sphere(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end
end

@testset "Plane" begin
    @testset "from_above" begin
        ray = Ray(Point(0f0, 0f0, 1f0), -VEC_Z)
        shape = Plane()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(0f0, 0f0, 0f0), NORMAL_Z_false, Vec2D(0f0, 0f0), 1f0, ray, Material())
    end

    @testset "laid_on" begin
        ray = Ray(Point(0f0, 0f0, 0f0), VEC_X)
        shape = Plane()
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end

    @testset "from_below" begin
        ray = Ray(Point(0f0, 0f0, -1f0), VEC_Z)
        shape = Plane()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(0f0, 0f0, 0f0), -NORMAL_Z_false, Vec2D(0f0, 0f0), 1f0, ray, Material())
    end

    @testset "transposed_from_above" begin
        ray = Ray(Point(1.5f0, 0.5f0, 3), -VEC_Z)
        shape = Plane(transformation = translation(Vec(0f0, 0f0, 1f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(1.5f0, 0.5f0, 1f0), NORMAL_Z_false, Vec2D(0.5f0, 0.5f0), 2f0, ray, Material())
    end

    @testset "from_diagonal" begin
        ray = Ray(Point(0f0, 0f0, 1f0), -VEC_Z + VEC_X + VEC_Y)
        shape = Plane()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(1f0, 1f0, 0f0), NORMAL_Z_false, Vec2D(0f0, 0f0), 1f0, ray, Material())
    end
end

@testset "Cube" begin
    @testset "from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Cube()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(0f0, 0f0, 5f-1), NORMAL_Z_false, Vec2D(7//8, 1//2), 1.5f0, ray, Material())
    end

    @testset "from_behind" begin
        ray = Ray(Point(3f0, 0f0, 0f0), -VEC_X)
        shape = Cube()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(5f-1, 0f0, 0f0), NORMAL_X_false, Vec2D(5//8, 1//2), 2.5f0, ray, Material())
    end

    @testset "from_within" begin
        ray = Ray(Point(0f0, 0f0, 0f0), VEC_X)
        shape = Cube()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(5f-1, 0f0, 0f0), -NORMAL_X_false, Vec2D(5//8, 1//2), 5f-1, ray, Material())
    end

    @testset "transposed_from_above" begin
        ray = Ray(Point(10f0, 0f0, 2f0), -VEC_Z)
        shape = Cube(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(10f0, 0f0, 5f-1), NORMAL_Z_false, Vec2D(7//8, 1//2), 1.5f0, ray, Material())
    end

    @testset "transposed_from_behind" begin
        ray = Ray(Point(13f0, 0f0, 0f0), -VEC_X)
        shape = Cube(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(10.5f0, 0f0, 0f0), NORMAL_X_false, Vec2D(5//8, 1//2), 2.5f0, ray, Material())
    end

    @testset "miss_transposed_from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Cube(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end

    @testset "miss_transposed_from_behind" begin
        ray = Ray(Point(-10f0, 0f0, 0f0), -VEC_X)
        shape = Cube(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end
end

@testset "Cylinder" begin
    @testset "from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Cylinder()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(0f0, 0f0, 5f-1), NORMAL_Z_false, Vec2D(3//4, 3//4), 1.5f0, ray, Material())
    end

    @testset "from_behind" begin
        ray = Ray(Point(3f0, 0f0, 0f0), -VEC_X)
        shape = Cylinder()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(5f-1, 0f0, 0f0), NORMAL_X_false, Vec2D(1//2, 1//2), 2.5f0, ray, Material())
    end

    @testset "from_within" begin
        ray = Ray(Point(0f0, 0f0, 0f0), VEC_X)
        shape = Cylinder()
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(5f-1, 0f0, 0f0), -NORMAL_X_false, Vec2D(1//2, 1//2), 5f-1, ray, Material())
    end

    @testset "transposed_from_above" begin
        ray = Ray(Point(10f0, 0f0, 2f0), -VEC_Z)
        shape = Cylinder(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(10f0, 0f0, 5f-1), NORMAL_Z_false, Vec2D(3//4, 3//4), 1.5f0, ray, Material())
    end

    @testset "transposed_from_behind" begin
        ray = Ray(Point(13f0, 0f0, 0f0), -VEC_X)
        shape = Cylinder(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect ≈ HitRecord(Point(10.5f0, 0f0, 0f0), NORMAL_X_false, Vec2D(1//2, 1//2), 2.5f0, ray, Material())
    end

    @testset "miss_transposed_from_above" begin
        ray = Ray(Point(0f0, 0f0, 2f0), -VEC_Z)
        shape = Cylinder(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end

    @testset "miss_transposed_from_behind" begin
        ray = Ray(Point(-10f0, 0f0, 0f0), -VEC_X)
        shape = Cylinder(transformation = translation(Vec(10f0, 0f0, 0f0)))
        intersect = ray_intersection(ray, shape)
        @test intersect === nothing
    end
end


text:
# Globals needed to support `entr` and other callbacks

"""
    Revise.revision_event

This `Condition` is used to notify `entr` that one of the watched files has changed.
"""
const revision_event = Condition()

"""
    Revise.user_callbacks_queue

Global variable, `user_callbacks_queue` holds `key` values for which the
file has changed but the user hooks have not yet been called.
"""
const user_callbacks_queue = Set{Any}()

"""
    Revise.user_callbacks_by_file

Global variable, maps files (identified by their absolute path) to the set of
callback keys registered for them.
"""
const user_callbacks_by_file = Dict{String, Set{Any}}()

"""
    Revise.user_callbacks_by_key

Global variable, maps callback keys to user hooks.
"""
const user_callbacks_by_key = Dict{Any, Any}()

"""
    key = Revise.add_callback(f, files, modules=nothing; key=gensym())

Add a user-specified callback, to be executed during the first run of
`revise()` after a file in `files` or a module in `modules` is changed on the
file system. If `all` is set to `true`, also execute the callback whenever any
file already monitored by Revise changes. In an interactive session like the
REPL, Juno or Jupyter, this means the callback executes immediately before
executing a new command / cell.

You can use the return value `key` to remove the callback later
(`Revise.remove_callback`) or to update it using another call
to `Revise.add_callback` with `key=key`.
"""
function add_callback(f, files, modules=nothing; all=false, key=gensym())
    fix_trailing(path) = isdir(path) ? joinpath(path, "") : path   # insert a trailing '/' if missing, see https://github.com/timholy/Revise.jl/issues/470#issuecomment-633298553

    remove_callback(key)

    files = map(fix_trailing, map(abspath, files))
    init_watching(files)

    # in case the `all` kwarg was set:
    # add all files which are already known to Revise
    if all
        for pkgdata in values(pkgdatas)
            append!(files, joinpath.(Ref(basedir(pkgdata)), srcfiles(pkgdata)))
        end
    end

    if modules !== nothing
        for mod in modules
            track(mod)  # Potentially needed for modules like e.g. Base
            id = PkgId(mod)
            pkgdata = pkgdatas[id]
            for file in srcfiles(pkgdata)
                absname = joinpath(basedir(pkgdata), file)
                push!(files, absname)
            end
        end
    end

    # There might be duplicate entries in `files`, but it shouldn't cause any
    # problem with the sort of things we do here
    for file in files
        cb = get!(Set, user_callbacks_by_file, file)
        push!(cb, key)
    end
    user_callbacks_by_key[key] = f

    return key
end

"""
    Revise.remove_callback(key)

Remove a callback previously installed by a call to `Revise.add_callback(...)`.
See its docstring for details.
"""
function remove_callback(key)
    for cbs in values(user_callbacks_by_file)
        delete!(cbs, key)
    end
    delete!(user_callbacks_queue, key)
    delete!(user_callbacks_by_key, key)

    # possible future work: we may stop watching (some of) these files
    # now. But we don't really keep track of what background tasks are running
    # and Julia doesn't have an ergonomic way of task cancellation yet (see
    # e.g.
    #     https://github.com/JuliaLang/Juleps/blob/master/StructuredConcurrency.md
    # so we'll omit this for now. The downside is that in pathological cases,
    # this may exhaust inotify resources.

    nothing
end

function process_user_callbacks!(keys = user_callbacks_queue; throw=false)
    try
        # use (a)sync so any exceptions get nicely collected into CompositeException
        @sync for key in keys
            f = user_callbacks_by_key[key]
            @async Base.invokelatest(f)
        end
    catch err
        if throw
            rethrow(err)
        else
            @warn "[Revise] Ignoring callback errors" err
        end
    finally
        empty!(keys)
    end
end

"""
    entr(f, files; all=false, postpone=false, pause=0.02)
    entr(f, files, modules; all=false, postpone=false, pause=0.02)

Execute `f()` whenever files or directories listed in `files`, or code in `modules`, updates.
If `all` is `true`, also execute `f()` as soon as code updates are detected in
any module tracked by Revise.

`entr` will process updates (and block your command line) until you press Ctrl-C.
Unless `postpone` is `true`, `f()` will be executed also when calling `entr`,
regardless of file changes. The `pause` is the period (in seconds) that `entr`
will wait between being triggered and actually calling `f()`, to handle
clusters of modifications, such as those produced by saving files in certain
text editors.

# Example

```julia
entr(["/tmp/watched.txt"], [Pkg1, Pkg2]) do
    println("update")
end
```
This will print "update" every time `"/tmp/watched.txt"` or any of the code defining
`Pkg1` or `Pkg2` gets updated.
"""
function entr(f::Function, files, modules=nothing; all=false, postpone=false, pause=0.02)
    yield()
    postpone || f()
    key = add_callback(files, modules; all=all) do
        sleep(pause)
        f()
    end
    try
        while true
            wait(revision_event)
            revise(throw=true)
        end
    catch err
        isa(err, InterruptException) || rethrow(err)
    finally
        remove_callback(key)
    end
    nothing
end