text:
defaultTableHelpText = """
PgUp/PgDn,
Arrow keys : standard navigation
<spc>,<rtn>: toggle leaf expansion
Home       : jump to the start
End        : jump to the end
ctrl_left/right arrow: paginate to left/right
[, ]       : make current column narrower/wider
ctrl_up    : move up to the start of the current branch or previous branch
ctrl_down  : move down to the next branch
+, -       : Expand or collapse 1 level everywhere
p          : Change pivot
c          : Change columns/order
v          : Switch preset views
/          : search text (on string columns) on exposed rows
?          : DEEP search text (on string columns), may expand nodes
ctrl_n     : search next occurence on exposed rows
ctrl_p     : search previous occurence on exposed rows
N          : search next occurence on any row. may expand more nodes
F6         : popup window for value
F7         : popup window for leaf & root node stats
"""
defaultTableBottomText = "F1:help p:Pivot c:ColOrd v:Views"
if VERSION < v"0.4.0-dev+1930"
    @lintpragma( "DataFrame is a container type" )
    @lintpragma( "AbstractDataFrame is a container type" )
else
    @lintpragma( "DataFrames.DataFrame is a container type" )
    @lintpragma( "DataFrames.AbstractDataFrame is a container type" )
end

type TwTableColInfo
    name::Symbol
    displayname::UTF8String
    format::FormatHints
    aggr::Any
end

type TwDfTableNode
    parent::WeakRef
    context::WeakRef
    pivotcols::Array{ Symbol, 1 }
    pivotvals::Tuple
    isOpen::Bool
    children::Array{ Any, 1} # only root node can have empty children
    subdataframe::Any
    subdataframesorted::Any
    colvalcache::Dict{ Symbol, Any } # order indep aggregated values
    function TwDfTableNode()
        new( WeakRef(), WeakRef(), Symbol[], (), false, Any[], nothing, nothing, Dict{Symbol,Any}() )
    end
    function TwDfTableNode( p::TwDfTableNode )
        x = TwDfTableNode()
        x.parent = WeakRef( p )
        x.context = p.context
        x
    end
end

function getindex( n::TwDfTableNode, c::Symbol )
    if haskey( n.colvalcache, c )
        return n.colvalcache[ c ]
    else
        context = n.context.value
        aggr = context.allcolInfo[ c ].aggr
        f = liftAggrSpecToFunc( c, aggr )
        ret = f( n.subdataframe )
        if typeof( ret ) <: AbstractDataFrame
            ret = ret[1][1] # first col, first row
        end

        n.colvalcache[c] = ret
        return ret
    end
end

type TwTableView
    name::UTF8String
    pivots::Array{Symbol,1}
    sortorder::Array{ (@compat Tuple{Symbol,Symbol}), 1 } # [ (:col1, :asc ), (:col2, :desc), ... ]
    columns::Array{ Symbol, 1 }
    initdepth::Int
end

# convenient functions to construct views
function TwTableView( df::AbstractDataFrame, name::UTF8String;
    pivots = Symbol[], initdepth=1,
    colorder = Any[ "*" ],
    hidecols = Any[], sortorder = Any[] )

    # construct visible columns in the right order
    function move_columns( targetarray::Array{ Symbol,1 }, sourcearray::Array, remaincols::Array{Symbol,1} )
        local i = 1
        local n = length( sourcearray )
        while !isempty( remaincols ) && i <= n
            local s = sourcearray[i]
            if s == "*"
                error( "illegal extra *")
            elseif typeof( s ) == Regex
                j = 1
                while j <= length( remaincols )
                    if match( s, string( remaincols[j] ) ) != nothing
                        push!( targetarray, remaincols[j] )
                        deleteat!( remaincols, j )
                        continue
                    end
                    j += 1
                end
            elseif typeof( s ) == Symbol
                if !in( s, remaincols )
                    error( "Column " * string( s ) * " not in " * string( remaincols )  * " src:" * string( sourcearray ) )
                end
                push!( targetarray, s )
                local idx = findfirst( remaincols, s )
                deleteat!( remaincols, idx )
            else
                error( "illegal " * string( s ) * " found" )
            end

            i += 1
        end
    end

    remaincols = names( df )
    if in( "*", colorder )
        colspreast = Symbol[]
        colspostast = Symbol[]
        local idx = findfirst( colorder, "*" )
        move_columns( colspreast, colorder[ 1:(idx-1) ], remaincols )
        finalcolorder = colspreast
        move_columns( colspostast, colorder[ (idx+1):length( colorder ) ], remaincols )
        append!( finalcolorder, remaincols )
        append!( finalcolorder, colspostast )
    else
        finalcolorder = Symbol[]
        move_columns( finalcolorder, colorder, remaincols )
    end

    # remove anything in hidecols
    removed = Symbol[]
    move_columns( removed, hidecols, finalcolorder )

    if eltype( sortorder ) == Symbol
        actualsortorder = (@compat Tuple{Symbol,Symbol})[ (s,:asc) for s in sortorder ]
    elseif eltype( sortorder ) == @compat Tuple{Symbol,Symbol}
        actualsortorder = sortorder
    else
        error( "sortorder eltype expects Symbol, or Tuple{Symbol,Symbol}: " * string( eltype( sortorder ) ) )
    end

    TwTableView( utf8(name), pivots, actualsortorder, finalcolorder, initdepth )
end

# this is the widget data. all subnodes hold a weakref back to this to
# facilitate aggregation, ordering and output
type TwDfTableData
    rootnode::TwDfTableNode
    pivots::Array{ Symbol, 1 }
    sortorder::Array{ (@compat Tuple{Symbol,Symbol}), 1 } # [ (:col1, :asc ), (:col2, :desc), ... ]
    datalist::Array{Any, 1} # (tuple{symbol}, 0) -> node, (tuple{symbol},#) -> row within that sub-df
    datalistlen::Int
    datatreewidth::Int
    headerlines::Int # how many lines does the header occupy, usually just 1
    currentTop::Int
    currentLine::Int
    currentCol::Int
    currentLeft::Int # left most on-screen column
    currentRight::Int # right most on-screen column
    colInfo::Array{ TwTableColInfo, 1 } # only the visible ones, maybe off-screen
    allcolInfo::Dict{ Symbol, TwTableColInfo } # including invisible ones
    bottomText::UTF8String
    helpText::UTF8String
    initdepth::Int
    views::Array{ TwTableView, 1 }
    calcpivots::Dict{ Symbol, CalcPivot }
    searchText::UTF8String
    # calculated dimension
    TwDfTableData() = new( TwDfTableNode(),
        Symbol[], Tuple{Symbol,Symbol}[], Any[], 0, 10, 1, 1, 1, 1, 1, 1, TwTableColInfo[],
        Dict{Symbol,TwTableColInfo}(), "", defaultTableHelpText, 1, TwTableView[], Dict{Symbol,CalcPivot}(),utf8("") )
end

#TODO: allow Regex in formatHints and aggrHints
function newTwDfTable( scr::TwObj, df::DataFrame;
        height::Real=1.0, width::Real=1.0, posy::Any=:center, posx::Any=:center,
        pivots = Symbol[],
        initdepth = 1,
        colorder = Any[ "*" ], # mix of symbol, regex, and "*" (the rest), "*" can be in the middle
        hidecols = Any[], # anything here trumps colorder, Symbol, or Regex
        sortorder = Tuple{Symbol,Symbol}[],
        title = "DataFrame",
        formatHints = Dict{Any,FormatHints}(), # Symbol/Type -> FormatHints
        aggrHints = Dict{Any,Any}(), # Symbol/Type -> string/symbol/expr/function
        widthHints = Dict{Symbol,Int}(),
        headerHints = Dict{Symbol,UTF8String}(),
        bottomText = defaultTableBottomText,
        views = Dict{Symbol,Any}[],
        calcpivots = Dict{Symbol,CalcPivot}() )
    obj = TwObj( TwDfTableData(), Val{:DfTable} )
    obj.value = df
    obj.title = title
    obj.box = true
    obj.borderSizeV= 1
    obj.borderSizeH= 2
    obj.data.rootnode.subdataframe = df
    obj.data.rootnode.context = WeakRef( obj.data )

    mainV = TwTableView( df, utf8( "#Main"), pivots = pivots,
        initdepth=initdepth, sortorder=sortorder, colorder=colorder, hidecols=hidecols )

    obj.data.pivots = mainV.pivots
    obj.data.initdepth = mainV.initdepth
    obj.data.sortorder = mainV.sortorder
    obj.data.calcpivots = calcpivots
    finalcolorder = mainV.columns

    push!( obj.data.views, mainV )
    for (i,d) in enumerate( views )
        if isempty( d )
            error( "nothing in view #" * string( i ) )
        end
        vname = get( d, :name, string( "v#" * string( i ) ) )
        vpivots = get( d, :pivots, pivots )
        vinitdepth = get( d, :initdepth, initdepth )
        vcolorder = get( d, :colorder, colorder )
        vhidecols = get( d, :hidecols, hidecols )
        vsortorder = get( d, :sortorder, sortorder )
        v = TwTableView( df, utf8( vname ), pivots = vpivots, initdepth = vinitdepth,
            sortorder=vsortorder, colorder = vcolorder, hidecols = vhidecols )
        push!( obj.data.views, v )
    end

    # construct colInfo for each col in finalcolorder
    allcols = names( df )
    for c in allcols

        if haskey( calcpivots, c )
            error( "calcpivots interfere with an existing column " * string( c ) )
        end

        t = eltype( df[ c ] )

        hdr = get( headerHints, c, string( c ) )
        fmt = get( formatHints, c,
                get( formatHints, t, deepcopy( FormatHints( t ) ) ) )
        if haskey( widthHints, c )
            fmt.width = widthHints[c]
        end
        agr = get( aggrHints, c,
                get( aggrHints, t, defaultAggr( t ) ) )
        ci = TwTableColInfo( c, hdr, fmt, agr )
        obj.data.allcolInfo[c] = ci
    end

    for c in finalcolorder
        ci = obj.data.allcolInfo[c]
        push!( obj.data.colInfo, ci )
    end

    expandnode( obj.data.rootnode, initdepth )
    ordernode( obj.data.rootnode )
    builddatalist( obj.data )

    updateTableDimensions( obj )
    obj.data.bottomText = bottomText
    link_parent_child( scr, obj, height, width, posy, posx )
    obj
end

function expandnode( n::TwDfTableNode, depth::Int=1 )
    if n.isOpen # nothing to do
        if depth > 1
            for r in n.children
                expandnode( r, depth-1 )
            end
        end
        return
    end
    pivots = n.context.value.pivots
    npivots = n.pivotcols # this is the node's pivot. length <= length( pivots )
    if length( npivots ) < length( pivots ) # populate children nodes
        if isempty( n.children )
            nextpivots = deepcopy( npivots )
            nextpivot = pivots[ length( npivots ) + 1 ]
            push!( nextpivots, nextpivot )

            if haskey( n.context.value.calcpivots, nextpivot )
                calcpvt = n.context.value.calcpivots[ nextpivot ]
                pvtspec = calcpvt.spec
                # if we have already pivoted :a, then including it
                # again in *by* would not do anything.
                pvtby   = setdiff( calcpvt.by, npivots )
                f = liftCalcPivotToFunc( pvtspec, pvtby )
                if isempty( pvtby )
                    colvalues = f( n.subdataframe )
                    # Note that setindex! doesn't work for subdataframe
                    # And we most certainly don't want to mutate the original
                    # dataframe (if the node n here is the rootnode)
                    # some sort of composit data frame is needed to
                    # avoid inefficient copying (or is it that bad?)
                    ns = names( n.subdataframe )
                    localdf = DataFrame( [ n.subdataframe[i] for i in 1:length(ns)], ns )
                    localdf[ nextpivot ] = colvalues
                    gd = DataFrames.groupby( localdf, nextpivots )
                else
                    # figure out the aggregation dependency
                    # the lift function just now ensures we have this cache.
                    aggrs = CalcPivotAggrDepCache[ (pvtspec, pvtby) ]
                    kwargs = Any[]
                    for a in aggrs
                        push!( kwargs, ( a, n.context.value.allcolInfo[ a ].aggr ) )
                    end
                    # the lifted function expects us to provide
                    # the aggregation spec on all needed columns,
                    # as keyword arguments
                    df = f( n.subdataframe, nextpivot; kwargs... )
                    gd = DataFrames.groupby( join( n.subdataframe, df, on=pvtby, kind=:left ), nextpivots )
                end
            else
                gd = DataFrames.groupby( n.subdataframe, nextpivots )
            end
            for g in gd
                dfr = DataFrameRow( g[ gd.cols ], 1 )
                valtuple = tuple( [ x[2] for x in dfr ]... )
                r = TwDfTableNode( n )
                r.pivotcols = nextpivots
                r.pivotvals = valtuple
                r.subdataframe = g
                push!( n.children, r )
            end
        end
        if depth > 1
            for r in n.children
                expandnode( r, depth-1 )
            end
        end
    end
    n.isOpen = true
end

# order the children, or if the terminal node, order the subdataframe
# also order opened children (recursive)
function ordernode( n::TwDfTableNode )
    pivots = n.context.value.pivots
    npivots = n.pivotcols # this is the node's pivot. length <= length( pivots )
    sortorder = n.context.value.sortorder
    if length( npivots ) < length( pivots ) # populate children nodes
        if length( sortorder ) > 0
            sort!( n.children, lt = (x,y) -> begin
                for sc in sortorder
                    if isna( x[sc[1]] )
                        if !isna( y[sc[1]] )
                            return false
                        else
                            continue
                        end
                    elseif isna( y[sc[1]] )
                        return true
                    end
                    if x[sc[1]] == y[sc[1]]
                        continue
                    end

                    if sc[2] == :desc
                        return y[sc[1]] < x[sc[1]]
                    else
                        return x[sc[1]] < y[sc[1]]
                    end
                end
                return false
            end )
        else
            sort!( n.children, lt = (x,y) -> x.pivotvals[end] < y.pivotvals[end] )
        end
        for c in n.children
            if c.isOpen
                ordernode( c )
            end
        end
    else
        if length( sortorder ) == 0
            n.subdataframesorted = n.subdataframe
        else
            n.subdataframesorted = sort(n.subdataframe,
                cols=Symbol[x[1] for x in sortorder],
                rev=Bool[x[2]==:desc for x in sortorder ] )
        end
    end
end

function builddatalist( o::TwDfTableData )
    o.datalist = Any[]
    pivots = o.pivots

    function presublist( subn::TwDfTableNode, substack::Array{Int,1}, skiplines::Array{Int,1}, islast::Bool )
        if islast
            newskip = copy(skiplines)
            push!( newskip, length(substack))
            sublist( subn, substack, newskip )
        else
            sublist( subn, substack, skiplines )
        end
    end

    function sublist( n::TwDfTableNode, stack::Array{Int,1}, skiplines::Array{Int,1} )
        if isempty( n.pivotvals )
            name = "Root"
        else
            name = string( n.pivotvals[end] )
        end

        push!( o.datalist, ( name, stack, n.isOpen ? :open : :close, skiplines, n ) )
        if n.isOpen
            if length( n.pivotcols ) < length( pivots )
                len = length( n.children )
                for (i,c) in enumerate( n.children )
                    substack = copy( stack )
                    push!( substack, i )
                    presublist( c, substack, skiplines, i==len )
                end
            else
                for i in 1:size( n.subdataframe,1 )
                    substack = copy( stack )
                    push!( substack, i )
                    push!( o.datalist, ( string(i), substack, :single, skiplines, n ) )
                end
            end
        end
    end

    sublist( o.rootnode, Int[], Int[] )
    o.currentLine = min( o.currentLine, length( o.datalist ) )
end

function updateTableDimensions( o::TwObj )
    o.data.datalistlen = length( o.data.datalist )
    o.data.headerlines = maximum( map( x->length( split( x.displayname, "\n" ) ), o.data.colInfo ) )
    # reminder: (name, stack, exphints, skiplines, node )
    o.data.datatreewidth = maximum( map( d -> 2*(length(d[2])+1) + length(d[1])+1,
        o.data.datalist ) )
end

function draw( o::TwObj{TwDfTableData} )
    updateTableDimensions( o )
    viewContentHeight = o.height - 2 * o.borderSizeV - o.data.headerlines
    viewContentWidth  = o.width - 2 * o.borderSizeH

    box( o.window, 0,0 )
    if !isempty( o.title )
        titlestr = o.title
        mvwprintw( o.window, 0, (@compat round(Int, ( o.width - length(titlestr) )/2 )), "%s", titlestr )
    end
    if o.data.datalistlen <= viewContentHeight
        msg = "ALL"
    else
        msg = @sprintf( "%d/%d %5.1f%%", o.data.currentLine, o.data.datalistlen,
            o.data.currentLine / o.data.datalistlen * 100 )
    end
    mvwprintw( o.window, 0, o.width - length(msg)-3, "%s", msg )
    updateTableDimensions( o )

    # header row(s)
    wattron( o.window, COLOR_PAIR(3) )
    startx = 1+o.data.datatreewidth
    lastcol = 1
    lastwidth = 8
    for col = o.data.currentLeft:length( o.data.colInfo )
        if startx > viewContentWidth
            break
        end
        ci = o.data.colInfo[col]
        lines = split( ci.displayname, "\n" )
        nlines= length(lines)
        width = ci.format.width

        islastcol = ( startx+width+6 > viewContentWidth ) || col == length( o.data.colInfo )
        if islastcol
            width = min( width, viewContentWidth-startx+2 )
        end
        o.data.currentRight = lastcol = col
        lastwidth = width

        if o.data.currentCol == col
            wattron( o.window, A_REVERSE )
        end
        for (i,line) in enumerate( lines )
            s = ensure_length( line, width )
            if i == nlines
                wattron( o.window, A_UNDERLINE )
            end
            mvwprintw( o.window, i+o.data.headerlines-nlines, startx, "%s", s )
            if i == nlines
                wattroff( o.window, A_UNDERLINE )
            end
        end
        if o.data.currentCol == col
            wattroff( o.window, A_REVERSE )
        end
        if !islastcol
            for i = 1:o.data.headerlines
                mvwaddch( o.window, i, startx+width, get_acs_val( 'x' ) )
            end
        end
        startx += width + 1
    end
    wattroff( o.window, COLOR_PAIR(3) )
    # reminder: (name, stack, exphints, skiplines, node )
    for r in o.data.currentTop:min( o.data.currentTop + viewContentHeight - 1, o.data.datalistlen )
        stacklen = length( o.data.datalist[r][2])

        # treecolume is always shown
        s = ensure_length( repeat( " ", 2*stacklen + 1) * o.data.datalist[r][1], o.data.datatreewidth-2 )

        if r == o.data.currentLine
            wattron( o.window, A_BOLD | COLOR_PAIR(o.hasFocus ? 15 : 30 ) )
        end
        mvwprintw( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2, "%s", s )
        for i in 1:stacklen - 1
            if !in( i, o.data.datalist[r][4] ) # skiplines
                mvwaddch( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2*i, get_acs_val( 'x' ) ) # vertical line
            end
        end
        if stacklen != 0
            contchar = get_acs_val('t') # tee pointing right
            if r == o.data.datalistlen ||  # end of the whole thing
                length(o.data.datalist[r+1][2]) < stacklen || # next one is going back in level
                ( length(o.data.datalist[r+1][2]) > stacklen && in( stacklen, o.data.datalist[r+1][4] ) ) # going deeper in level
                contchar = get_acs_val( 'm' ) # LL corner
            end
            mvwaddch( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2*stacklen, contchar )
            mvwaddch( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2*stacklen+1, get_acs_val('q') ) # horizontal line
        end
        if o.data.datalist[r][3] == :close
            mvwprintw( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2*stacklen+2, "%s", string( @compat Char( 0x25b8 ) ) ) # right-pointing small triangle
        elseif o.data.datalist[r][3] == :open
            mvwprintw( o.window, o.data.headerlines + 1+r-o.data.currentTop, 2*stacklen+2, "%s", string( @compat Char( 0x25be ) ) ) # down-pointing small triangle
        end

        if r == o.data.currentLine
            wattroff( o.window, A_BOLD | COLOR_PAIR(o.hasFocus ? 15 : 30 ) )
        end
        mvwaddch( o.window, o.data.headerlines+1+r-o.data.currentTop, o.data.datatreewidth, get_acs_val( 'x' ) )

        # other columns
        # get the node or DataFrameRow first
        node = o.data.datalist[r][5]
        isnode = (o.data.datalist[r][3] != :single)
        startx = 1+o.data.datatreewidth
        underline =  r < o.data.datalistlen && length( o.data.datalist[r+1][2] ) < stacklen
        for col = o.data.currentLeft:lastcol
            cn = o.data.colInfo[ col ].name
            if isnode
                v = node[ cn ]
            else
                v = node.subdataframesorted[ cn ][ o.data.datalist[r][2][end] ]
            end
            width = ( col == lastcol ? lastwidth : o.data.colInfo[ col ].format.width )
            isred = false
            if typeof( v ) == NAtype
                str = ensure_length( "", width )
            elseif typeof( v ) <: Real
                str = applyformat( v, o.data.colInfo[col].format )
                if length(str) > width # all ####
                    str = repeat( "#", width )
                elseif length(str) < width
                    str = repeat( " ", width - length(str) ) * str
                end
                isred =  v < 0
            else
                str = applyformat( v, o.data.colInfo[col].format )
                str = ensure_length( str, width )
            end
            flags = underline ? A_UNDERLINE : 0
            if col == o.data.currentCol && r == o.data.currentLine
                flags |= A_BOLD
                if isred
                    flags |= COLOR_PAIR(o.hasFocus ? 9 : 31)
                else
                    flags |= COLOR_PAIR(o.hasFocus ? 15 : 30 )
                end
            elseif isnode
                flags |= A_BOLD
                if  mod( length( o.data.datalist[r][2] ), 2 ) == 0
                    if isred
                        flags |= COLOR_PAIR(1)
                    else
                        flags |= COLOR_PAIR(7)
                    end
                else
                    if isred
                        flags |= COLOR_PAIR(29)
                    else
                        flags |= COLOR_PAIR(13)
                    end
                end
            else
                if isred
                    flags |= COLOR_PAIR(1)
                else
                    flags |= COLOR_PAIR(7)
                end
            end
            wattron( o.window, flags )
            mvwprintw( o.window, o.data.headerlines + 1+r-o.data.currentTop, startx, "%s", str )
            wattroff( o.window, flags )
            if col != lastcol
                mvwaddch( o.window, o.data.headerlines + 1+ r-o.data.currentTop, startx+width, get_acs_val( 'x' ) )
            end
            startx += width + 1
        end
    end
    bottomtext = o.data.bottomText
    pivottext = ""
    if !isempty( o.data.pivots )
        pivottext = "▾"*join( o.data.pivots, "▾" )
    end
    if length( bottomtext ) + length( pivottext ) > o.width-4
        if length( pivottext ) < o.width - 12 # just shorten helptext
            bottomtext = ensure_length( bottomtext, o.width-4-length(pivottext ) )
        else
            bottomtext = ensure_length( bottomtext, 8 )
            pivottext = ensure_length( pivottext, o.width-12 )
        end
    end
    if length( bottomtext ) != 0
        mvwprintw( o.window, o.height-1, 3, "%s", bottomtext )
    end
    if length( pivottext ) != 0
        mvwprintw( o.window, o.height-1, o.width - length( pivottext ) - 1, "%s", pivottext )
    end
end

function inject( o::TwObj{TwDfTableData}, token )
    dorefresh = false
    retcode = :got_it # default behavior is that we know what to do with it
    viewContentHeight = o.height - 2 * o.borderSizeV - o.data.headerlines
    viewContentWidth = o.width - 2* o.borderSizeH - o.data.datatreewidth
    widths = map( x->x.format.width, o.data.colInfo )

    update_tree_data = ()->begin
        builddatalist( o.data )
    end

    checkTop = () -> begin
        o.data.datalistlen = length( o.data.datalist )
        if o.data.currentLine > o.data.datalistlen
            o.data.currentLine = o.data.datalistlen
        end
        if o.data.currentTop < 1
            o.data.currentTop = 1
        end
        if o.data.currentTop > o.data.datalistlen - viewContentHeight + 1
            o.data.currentTop = max(1,o.data.datalistlen - viewContentHeight + 1)
        end
        if o.data.currentTop > o.data.currentLine
            o.data.currentTop = o.data.currentLine
        end
        if o.data.currentLine - o.data.currentTop > viewContentHeight-1
            o.data.currentTop = o.data.currentLine - viewContentHeight+1
        end
    end
    checkLeft = () -> begin
        if o.data.currentLeft < 1
            o.data.currentLeft = 1
        else # check if we have enough width to show from currentLeft to currentCol
            revcumwidths = cumsum( map( x->x+1, reverse( widths ) ) ) # with boundary
            widthrng = searchsorted( revcumwidths, viewContentWidth )
            if o.data.currentLeft > length( o.data.colInfo ) - widthrng.stop +1
                o.data.currentLeft = length( o.data.colInfo ) - widthrng.stop +1
            end
        end
        if o.data.currentLeft > o.data.currentCol
            o.data.currentLeft = o.data.currentCol
        else
            revcumwidths = cumsum( map( x->x+1, reverse( widths[o.data.currentLeft:o.data.currentCol] ) ) ) # with boundary
            widthrng = searchsorted( revcumwidths, viewContentWidth )
            if o.data.currentLeft < o.data.currentCol - widthrng.stop +1
                o.data.currentLeft = o.data.currentCol - widthrng.stop +1
            end
        end
    end
    movevertical = n -> begin
        oldline = o.data.currentLine
        o.data.currentLine = max(1, min( o.data.datalistlen, o.data.currentLine + n) )
        if oldline != o.data.currentLine
            checkTop()
            return true
        else
            beep()
            return false
        end
    end
    movehorizontal = n -> begin
        oldcol = o.data.currentCol
        o.data.currentCol = max(1, min( length( o.data.colInfo ), o.data.currentCol + n) )
        if oldcol != o.data.currentCol
            checkLeft()
            return true
        else
            beep()
            return false
        end
    end

    function searchNext( step::Int, trivialstop::Bool )
        local st = o.data.currentLine
        o.data.searchText = lowercase(o.data.searchText)
        i = trivialstop ? st : ( mod( st-1+step, o.data.datalistlen ) + 1 )
        while true
            node = o.data.datalist[i][5]
            isnode = (o.data.datalist[i][3] != :single)
            # only visible columns and exposed rows are search
            for col in 1:length( o.data.colInfo )
                cn = o.data.colInfo[ col ].name
                if isnode
                    v = node[ cn ]
                else
                    v = node.subdataframesorted[ cn ][ o.data.datalist[i][2][end] ]
                end
                if !( typeof( v ) <: AbstractString )
                    continue
                end

                if contains( lowercase( v ), o.data.searchText )
                    o.data.currentLine = i
                    o.data.currentCol = col
                    checkTop()
                    checkLeft()
                    return i
                end
            end
            i = mod( i-1+step, o.data.datalistlen ) + 1
            if i == st
                beep()
                return 0
            end
        end
    end

    # deep search always go forward
    function searchNextDeep( trivialstop::Bool )
        local st = o.data.currentLine
        local stp = 1

        o.data.searchText = lowercase(o.data.searchText)
        i = trivialstop ? st : ( mod( st-1+stp, o.data.datalistlen ) + 1 )
        ncols = length( o.data.colInfo )
        function checknode( nd::TwDfTableNode, substack::Array{Int,1} )
            function searchdfstring( df::AbstractDataFrame )
                for j in 1:nrow( df )
                    for col in 1:ncols
                        cn = o.data.colInfo[ col ].name
                        v = df[cn][j]
                        if !(typeof(v)<:AbstractString)
                            continue
                        end
                        if contains( lowercase( v ), o.data.searchText )
                            o.data.currentCol = col
                            return j
                        end
                    end
                end
                return 0
            end
            if searchdfstring( nd.subdataframe ) != 0
                if !nd.isOpen
                    expandnode( nd )
                    ordernode( nd )
                end
                if !isempty( nd.children )
                    substacktmp = copy( substack )
                    push!( substacktmp, 0 )
                    for (k,c) in enumerate( nd.children )
                        substacktmp[end] = k
                        ret = checknode( c, substacktmp )
                        if ret != nothing
                            return ret
                        end
                    end
                else
                    therow = searchdfstring( nd.subdataframesorted )
                    substacktmp = copy( substack )
                    push!( substacktmp, therow )
                    #log( "found on " * string( substacktmp ) )
                    return substacktmp
                end
            end
            return nothing
        end
        while true
            node = o.data.datalist[i][5]
            isnode = (o.data.datalist[i][3] != :single)
            if isnode && !node.isOpen
                foundstack = checknode( node, o.data.datalist[i][2] )
                if foundstack != nothing
                    builddatalist( o.data )
                    o.data.currentLine = searchsortedfirst( o.data.datalist, Any[ nothing, foundstack ], by=y->y[2], lt=(s1,s2)->begin
                        #log( "comparing " * string(s1) * " with " * string(s2) )
                        for j in 1:min( length(s1), length(s2) )
                            if s1[j] == s2[j]
                                continue
                            end
                            #log( "Pivot level " * string( j ) * " index is different" )
                            return s1[j] < s2[j]
                        end
                        #log( "All common indices are equal. Pivot lengths are compared" )
                        return length(s1) < length(s2)
                    end )
                    checkTop()
                    checkLeft()
                    return o.data.currentLine
                end
            end
            if !isnode
                for col in 1:ncols
                    cn = o.data.colInfo[ col ].name
                    v = node.subdataframesorted[ cn ][ o.data.datalist[i][2][end] ]
                    if !( typeof( v ) <: AbstractString )
                        continue
                    end

                    if contains( lowercase( v ), o.data.searchText )
                        o.data.currentLine = i
                        o.data.currentCol = col
                        checkTop()
                        checkLeft()
                        return i
                    end
                end
            end
            i = mod( i-1+stp, o.data.datalistlen ) + 1
            if i == st
                beep()
                return 0
            end
        end
    end

    # reminder: (name, stack, exphints, skiplines, node )
    if token == :esc
        retcode = :exit_nothing
    elseif ( token == " " || token == Symbol( "return" ) || token == :enter ) && o.data.datalist[ o.data.currentLine ][3] != :single
        expandhint = o.data.datalist[ o.data.currentLine ][3]
        node = o.data.datalist[ o.data.currentLine][5]
        if node.isOpen
            node.isOpen = false
        else
            expandnode( node )
            ordernode( node )
        end
        update_tree_data()
        checkTop()
        dorefresh = true
    elseif token == "+"
        somethingchanged=false
        for r in 1:length( o.data.datalist )
            if o.data.datalist[r][3] == :close
                node = o.data.datalist[r][5]
                expandnode( node )
                ordernode( node )
                somethingchanged=true
            end
        end
        if somethingchanged
            update_tree_data()
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif token == "-"
        somethingchanged=false
        for r in 1:length( o.data.datalist )
            if o.data.datalist[r][3] == :open
                node = o.data.datalist[r][5]
                if !isempty( node.children )
                    if all( x->!x.isOpen, node.children )
                        node.isOpen = false
                        somethingchanged=true
                    end
                else
                    node.isOpen = false
                    somethingchanged=true
                end
            end
        end
        if somethingchanged
            update_tree_data()
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif token == :up
        dorefresh = movevertical(-1)
    elseif token == :down
        dorefresh = movevertical(1)
    elseif token == :left
        dorefresh = movehorizontal(-1)
    elseif token == :ctrl_left
        if o.data.currentCol != o.data.currentLeft
            o.data.currentCol = o.data.currentLeft
            dorefresh = true
        elseif o.data.currentLeft == 1
            beep()
        else
            # page left
            revcumwidths = cumsum( map( x->x+1, reverse( widths[1:o.data.currentLeft] ) ) ) # with boundary
            widthrng = searchsorted( revcumwidths, viewContentWidth )
            o.data.currentLeft = o.data.currentCol = max( 1, o.data.currentLeft - widthrng.start + 1 )
            checkLeft()
            dorefresh=true
        end
    elseif token == :right
        dorefresh = movehorizontal(1)
    elseif token == :ctrl_right
        if o.data.currentCol != o.data.currentRight
            o.data.currentCol = o.data.currentRight
            dorefresh = true
        elseif o.data.currentRight == length( o.data.colInfo )
            beep()
        else
            cumwidths = cumsum( map( x->x+1, reverse( widths[o.data.currentRight:end] ) ) ) # with boundary
            widthrng = searchsorted( cumwidths, viewContentWidth )
            o.data.currentRight = o.data.currentCol = min( o.data.currentRight + widthrng.stop, length( o.data.colInfo ) )
            checkLeft()
            dorefresh=true
        end
    elseif token == :pageup
        dorefresh = movevertical( -viewContentHeight + o.data.headerlines)
    elseif token == :pagedown
        dorefresh = movevertical( viewContentHeight - o.data.headerlines )
    elseif  token == :home
        if o.data.currentTop != 1 || o.data.currentLeft != 1 || o.data.currentLine != 1 || o.data.currentCol != 1
            o.data.currentTop = 1
            o.data.currentLeft = 1
            o.data.currentLine = 1
            o.data.currentCol = 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :KEY_MOUSE
        (mstate,x,y, bs ) = getmouse()
        if mstate == :scroll_up
            dorefresh = movevertical( -(@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :scroll_down
            dorefresh = movevertical( (@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :button1_pressed
            rely, relx = screen_to_relative( o.window, y, x )
            if 1<=relx<o.width-1 && o.data.headerlines<rely<o.height-1
                o.data.currentLine = min( o.data.datalistlen, o.data.currentTop + rely - o.borderSizeH + 1 - o.data.headerlines )
                checkTop()
                dorefresh = true
            end
            if o.data.datatreewidth+1<relx<o.width-1 && o.data.headerlines<=rely<o.height-1
                cumwidths = cumsum( map( _->_+1, widths[o.data.currentLeft:end] ) ) # with boundary
                widthrng = searchsorted( cumwidths, relx - o.data.datatreewidth - 1)
                o.data.currentCol = min( length( o.data.colInfo ), o.data.currentLeft + widthrng.start - 1 )
                checkLeft()
                dorefresh = true
            end
            if !dorefresh
                retcode = :pass
            end
        end
    elseif in( token, Any[ Symbol("end") ] )
        if o.data.currentTop + viewContentHeight -1 < o.data.datalistlen
            o.data.currentTop = o.data.datalistlen - viewContentHeight + 1
            o.data.currentLine = o.data.datalistlen
            dorefresh = true
        else
            beep()
        end
    elseif token == "["
        width = o.data.colInfo[ o.data.currentCol ].format.width
        if width > 4
            width -=1
            o.data.colInfo[ o.data.currentCol ].format.width = width
            dorefresh = true
        else
            beep()
        end
    elseif token == "]"
        width = o.data.colInfo[ o.data.currentCol ].format.width
        if width < viewContentWidth-1
            width +=1
            o.data.colInfo[ o.data.currentCol ].format.width = width
            dorefresh = true
        else
            beep()
        end
    elseif token == :ctrl_down
        curr = o.data.currentLine
        stck = o.data.datalist[curr][2]
        if isempty( stck )
            beep()
        else
            if o.data.datalist[curr][3] == :single
                tmpstack = copy( stck )
                pop!( tmpstack )
            else
                tmpstack = stck
            end
            for r in curr+1:length(o.data.datalist )
                rstack = o.data.datalist[r][2]
                if o.data.datalist[r][3] != :single && ( length(rstack) <= length( tmpstack) || rstack[ 1:length( tmpstack ) ] == tmpstack )
                    if r != curr
                        o.data.currentLine = r
                        checkTop()
                        dorefresh = true
                    end
                    break
                end
            end
        end
    elseif token == :ctrl_up
        curr = o.data.currentLine
        stck = o.data.datalist[curr][2]
        if isempty( stck )
            beep()
        else
            tmpstack = copy( stck )
            pop!( tmpstack )
            for r in curr-1:-1:1
                rstack = o.data.datalist[r][2]
                if o.data.datalist[r][3] != :single && length(rstack) <= length( tmpstack)
                    if r != curr
                        o.data.currentLine = r
                        checkTop()
                        dorefresh = true
                    end
                    break
                end
            end
        end
    elseif token == "p"
        allcols = map(_->utf8(string(_)), names( o.data.rootnode.subdataframe ) )
        append!( allcols, map( _->utf8(string(_)), collect( keys( o.data.calcpivots ) ) ) )
        pvts = map( _->utf8(string(_)), o.data.pivots )
        helper = newTwMultiSelect( o.screen.value, allcols, selected = pvts, title="Pivot order", orderable=true, substrsearch=true )
        newpivots = activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        if newpivots != nothing && newpivots != pvts
            o.data.pivots = Symbol[ Symbol( _ ) for _ in newpivots ]
            o.data.rootnode.children = Any[]
            o.data.rootnode.isOpen = false
            expandnode( o.data.rootnode, o.data.initdepth )
            ordernode( o.data.rootnode )
            update_tree_data()
            o.data.currentLine = 1
            checkTop()
        end
        dorefresh = true
    elseif token == "c"
        allcols = map(_->utf8(string(_)), names( o.data.rootnode.subdataframe ) )
        visiblecols = map( _->utf8(string(_.name)), o.data.colInfo )
        helper = newTwMultiSelect( o.screen.value, allcols, selected = visiblecols, title="Visible columns & their order", orderable=true, substrsearch=true )
        newcols = activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        if newcols != nothing && newcols != visiblecols
            o.data.colInfo = TwTableColInfo[]
            for c in newcols
                push!( o.data.colInfo, o.data.allcolInfo[ Symbol( c ) ] )
            end
        end
        dorefresh = true
    elseif token == "v"
        allviews = map( _->_.name, o.data.views )
        helper = newTwPopup( o.screen.value, allviews, substrsearch=true, title = "Views" )
        vname = activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        if vname != nothing
            idx = findfirst( _->_.name == vname, o.data.views )
            v = o.data.views[idx]
            o.data.colInfo = TwTableColInfo[]
            o.data.pivots = v.pivots
            o.data.sortorder = v.sortorder
            o.data.initdepth = v.initdepth
            for c in v.columns
                push!( o.data.colInfo, o.data.allcolInfo[ Symbol( c ) ] )
            end
            o.data.rootnode.children = Any[]
            o.data.rootnode.isOpen = false
            expandnode( o.data.rootnode, o.data.initdepth )
            ordernode( o.data.rootnode )
            update_tree_data()
            o.data.currentLine = 1
            checkTop()
            dorefresh=true
        end
    elseif token == "/"
        helper = newTwEntry( o.screen.value, UTF8String; width=30, posy=:center, posx=:center, title = "Search: " )
        helper.data.inputText = o.data.searchText
        s = activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        if s != nothing
            if s != "" && o.data.searchText != s
                o.data.searchText = s
                searchNext( 1, true )
            end
        end
        dorefresh = true
    elseif token == "?"
        helper = newTwEntry( o.screen.value, UTF8String; width=30, posy=:center, posx=:center, title = "Search: " )
        helper.data.inputText = o.data.searchText
        s = activateTwObj( helper )
        unregisterTwObj( o.screen.value, helper )
        if s != nothing
            o.data.searchText = s
            searchNextDeep( true )
        end
        dorefresh = true
    elseif token == :ctrl_n
        searchNext( 1, false )
        dorefresh = true
    elseif token == :ctrl_p
        searchNext( -1, false )
        dorefresh = true
    elseif token == "N"
        searchNextDeep( false )
        dorefresh = true
    elseif token == :F6
        colsym = o.data.colInfo[ o.data.currentCol ].name
        node = o.data.datalist[o.data.currentLine][5]
        isnode = (o.data.datalist[o.data.currentLine][3] != :single)
        if isnode
            v = node[ colsym ]
        else
            v = node.subdataframesorted[ colsym ][ o.data.datalist[o.data.currentLine][2][end] ]
        end
        if typeof( v ) != NAtype && !in( v, [ nothing, Void, Any ] )
            tshow( v; title = string( colsym ), posx=:center, posy=:center )
            dorefresh = true
        end
    elseif token == :F7
        colsym = o.data.colInfo[ o.data.currentCol ].name
        node = o.data.datalist[o.data.currentLine][5]
        out = IOBuffer()
        describe( out, node.subdataframe[ colsym ] )

        if node != o.data.rootnode
            println( out, "\nRoot table stats" )
            describe( out, o.data.rootnode.subdataframe[ colsym ] )
        end
        tshow( takebuf_string( out ); title = string( colsym )  * " stats", posx=:center,posy=:center)
        dorefresh = true
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        refresh(o)
    end

    return retcode
end

helptext( o::TwObj{TwDfTableData} ) = o.data.helpText


text:
using Tensors: ∇, ∇∇, Δ
import ForwardDiff

function Ψ(C, μ, Kb)
    detC = det(C)
    J = sqrt(detC)
    Ĉ = detC^(-1/3)*C
    return 1/2*(μ * (tr(Ĉ)- 3) + Kb*(J-1)^2)
end

function S(C, μ, Kb)
    I = one(C)
    J = sqrt(det(C))
    invC = inv(C)
    return μ * det(C)^(-1/3)*(I - 1/3*tr(C)*invC) + Kb*(J-1)*J*invC
end

const μ = 1e10;
const Kb = 1.66e11;
Ψ(C) = Ψ(C, μ, Kb)
S(C) = S(C, μ, Kb)

@testsection "AD" begin
    for dim in 1:3
        @testsection "dim $dim" begin
        F = one(Tensor{2,dim}) + rand(Tensor{2,dim});
        C = tdot(F);
        C2 = F' ⋅ F;

        @test 2(@inferred ∇(Ψ, C))::typeof(C) ≈ S(C)
        @test 2(@inferred ∇(Ψ, C2))::typeof(C2) ≈ S(C2)

        b = rand(SymmetricTensor{2, dim})
        @test 2 * (@inferred ∇∇(Ψ, C))::SymmetricTensor{4, dim} ⊡ b ≈ ∇(S, C) ⊡ b
        @test 2 * (@inferred ∇∇(Ψ, C2))::Tensor{4, dim} ⊡ b ≈ ∇(S, C2) ⊡ b

        @test ∇(Ψ, C) ≈ ∇(Ψ, C2)
        @test ∇(S, C) ⊡ b ≈ ∇(S, C2) ⊡ b
        @test ∇∇(Ψ, C) ⊡ b ≈ ∇∇(Ψ, C2) ⊡ b

        for T in (Float32, Float64)
            Random.seed!(1234) # needed for getting "good" tensors for calculating det and friends
            A = rand(Tensor{2, dim, T})
            B = rand(Tensor{2, dim, T})
            A_sym = rand(SymmetricTensor{2, dim, T})
            B_sym = rand(SymmetricTensor{2, dim, T})
            v = rand(Vec{dim, T})
            II = one(Tensor{4, dim, T})
            II_sym = one(SymmetricTensor{4, dim, T})

            # Gradient of scalars
            @testsection "scalar grad" begin
                @test (@inferred ∇(norm, v))::typeof(v) ≈ ((@inferred ∇(norm, v, :all))[1])::typeof(v) ≈ v / norm(v)
                @test (∇(norm, v, :all)[2])::T ≈ norm(v)
                @test (@inferred ∇(norm, A))::typeof(A) ≈ ((@inferred ∇(norm, A, :all))[1])::typeof(A) ≈ A / norm(A)
                @test (∇(norm, A, :all)[2])::T ≈ norm(A)
                @test (@inferred ∇(norm, A_sym))::typeof(A_sym) ≈ ((@inferred ∇(norm, A_sym, :all))[1])::typeof(A_sym) ≈ A_sym / norm(A_sym)
                @test (∇(norm, A_sym, :all)[2])::T ≈ norm(A_sym)

                @test ∇(v -> 3*v, v) ≈ ∇(v -> 3*v, v, :all)[1] ≈ diagm(Tensor{2, dim}, 3.0)
                @test ∇(v -> 3*v, v, :all)[2] ≈ 3*v
                # function does not return dual
                @test ∇(A -> T(1), A)::typeof(A) ≈ ∇(A -> T(1), A, :all)[1] ≈ 0*A
                @test ∇(A -> T(1), A, :all)[2] == 1
                @test ∇(A -> T(1), A_sym)::typeof(A_sym) ≈ ∇(A -> T(1), A, :all)[1] ≈ 0*A_sym
                @test ∇(A -> T(1), A_sym, :all)[2] == 1
            end

            @testsection "2nd tensor grad" begin
            # Gradient of second order tensors
            # https://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivatives_of_the_invariants_of_a_second-order_tensor
                I1, DI1 = A -> tr(A), A -> one(A)
                I2, DI2 = A -> (tr(A)^2 - tr(A⋅A)) / 2, A -> I1(A) * one(A) - A'
                I3, DI3 = A -> det(A), A -> det(A) * inv(A)'

                @test (@inferred ∇(I1, A))::typeof(A) ≈ ((@inferred ∇(I1, A, :all))[1])::typeof(A) ≈ DI1(A)
                @test (∇(I1, A, :all)[2])::T ≈ I1(A)
                @test (@inferred ∇(I2, A))::typeof(A) ≈ ((@inferred ∇(I2, A, :all))[1])::typeof(A) ≈ DI2(A)
                @test (∇(I2, A, :all)[2])::T ≈ I2(A)
                @test (@inferred ∇(I3, A))::typeof(A) ≈ ((@inferred ∇(I3, A, :all))[1])::typeof(A) ≈ DI3(A)
                @test (∇(I3, A, :all)[2])::T ≈ I3(A)
                @test (@inferred ∇(I1, A_sym))::typeof(A_sym) ≈ ((@inferred ∇(I1, A_sym, :all))[1])::typeof(A_sym) ≈ DI1(A_sym)
                @test (∇(I1, A_sym, :all)[2])::T ≈ I1(A_sym)
                @test (@inferred ∇(I2, A_sym))::typeof(A_sym) ≈ ((@inferred ∇(I2, A_sym, :all))[1])::typeof(A_sym) ≈ DI2(A_sym)
                @test (∇(I2, A_sym, :all)[2])::T ≈ I2(A_sym)
                @test (@inferred ∇(I3, A_sym))::typeof(A_sym) ≈ ((@inferred ∇(I3, A_sym, :all))[1])::typeof(A_sym) ≈ DI3(A_sym)
                @test (∇(I3, A_sym, :all)[2])::T ≈ I3(A_sym)

                @test (@inferred ∇(identity, A))::Tensor{4, dim, T} ≈ ((@inferred ∇(identity, A, :all))[1])::Tensor{4, dim, T} ≈ II
                @test (∇(identity, A, :all)[2])::typeof(A) ≈ A
                @test (@inferred ∇(identity, A_sym))::SymmetricTensor{4, dim, T} ≈ ((@inferred ∇(identity, A_sym, :all))[1])::SymmetricTensor{4, dim, T} ≈ II_sym
                @test (∇(identity, A_sym, :all)[2])::typeof(A_sym) ≈ A_sym
                @test (@inferred ∇(transpose, A))::Tensor{4, dim, T} ⊡ B ≈ ((@inferred ∇(transpose, A, :all))[1])::Tensor{4, dim, T} ⊡ B ≈ B'
                @test (∇(transpose, A, :all)[2])::typeof(A) ≈ A'
                @test (@inferred ∇(transpose, A_sym))::SymmetricTensor{4, dim, T} ⊡ B_sym ≈ ((@inferred ∇(transpose, A_sym, :all))[1])::SymmetricTensor{4, dim, T} ⊡ B_sym ≈ B_sym'
                @test (∇(transpose, A_sym, :all)[2])::typeof(A_sym) ≈ A_sym'
                @test (@inferred ∇(inv, A))::Tensor{4, dim, T} ⊡ B ≈ ((@inferred ∇(inv, A, :all))[1])::Tensor{4, dim, T} ⊡ B ≈ - inv(A) ⋅ B ⋅ inv(A)
                @test (∇(inv, A, :all)[2])::typeof(A) ≈ inv(A)
                @test (@inferred ∇(inv, A_sym))::SymmetricTensor{4, dim, T} ⊡ B_sym ≈ ((@inferred ∇(inv, A_sym, :all))[1])::SymmetricTensor{4, dim, T} ⊡ B_sym ≈ - inv(A_sym) ⋅ B_sym ⋅ inv(A_sym)
                @test (∇(inv, A_sym, :all)[2])::typeof(A_sym) ≈ inv(A_sym)
                # function does not return dual
                @test ∇(A -> B, A) ≈ ∇(A -> B, A, :all)[1] ≈ zero(Tensor{4, dim, T})
                @test ∇(A -> B, A, :all)[2] ≈ B
                @test isa(∇(A -> B, A), Tensor{4, dim, T})
            end

            # Hessians of scalars
            @testsection "hessian" begin
                @test (@inferred ∇∇(norm, A))::Tensor{4, dim, T} ≈ ((@inferred ∇∇(norm, A, :all))[1])::Tensor{4, dim, T} ≈ reshape(ForwardDiff.hessian(x -> sqrt(sum(abs2, x)), A), (dim, dim, dim, dim))
                @test (∇∇(norm, A, :all)[2])::typeof(A) ≈ reshape(ForwardDiff.gradient(x -> sqrt(sum(abs2, x)), A), (dim, dim))
                @test (∇∇(norm, A, :all)[3])::T ≈ norm(A)
                @test (@inferred ∇∇(norm, A_sym))::SymmetricTensor{4, dim, T} ≈ ((@inferred ∇∇(norm, A_sym, :all))[1])::SymmetricTensor{4, dim, T}
                @test (∇∇(norm, A_sym, :all)[2])::typeof(A_sym) ≈ ∇(norm, A_sym)
                @test (∇∇(norm, A_sym, :all)[3])::T ≈ norm(A_sym)
                # function does not return dual
                @test ∇∇(A -> T(1), A)::Tensor{4, dim, T} ≈ ∇∇(A -> T(1), A, :all)[1] ≈ 0*II
                @test ∇∇(A -> T(1), A, :all)[2] ≈ 0*A
                @test ∇∇(A -> T(1), A, :all)[3] == T(1)
            end
            end # loop T
        end # testsection
    end # loop dim

    @testsection "vector calculus identities" begin
        Random.seed!(1234)
        φ(x) = norm(x)^4
        ϕ(x) = sum(x)
        A(x) = Vec{3}((x[1]*x[2]^3*x[3], x[1]*x[2]*x[3]^3, x[1]^3*x[2]*x[3]))
        B(x) = Vec{3}((x[1]*x[1], x[1]*x[2], x[1]*x[3]))
        for T in (Float32, Float64)
            x = rand(Vec{3, T})
            # gradient
            @test gradient(x -> φ(x) + ϕ(x), x) ≈ gradient(φ, x) + gradient(ϕ, x)
            @test gradient(x -> φ(x) * ϕ(x), x) ≈ φ(x) * gradient(ϕ, x) + gradient(φ, x) * ϕ(x)
            @test gradient(x -> A(x) ⋅ B(x), x) ≈ gradient(B, x)⋅A(x) + gradient(A, x)⋅B(x) + A(x)×curl(B, x) + B(x)×curl(A, x)
            # divergence
            @test divergence(x -> A(x) + B(x), x) ≈ divergence(A, x) + divergence(B, x)
            @test divergence(x -> φ(x) * A(x), x) ≈ φ(x)*divergence(A, x) + gradient(φ, x)⋅A(x)
            @test divergence(x -> A(x) × B(x), x) ≈ B(x)⋅curl(A, x) - A(x)⋅curl(B, x)
            # curl
            @test curl(x -> A(x) + B(x), x) ≈ curl(A, x) + curl(B, x)
            @test curl(x -> φ(x) * A(x), x) ≈ φ(x)*curl(A, x) + gradient(φ, x)×A(x)
            @test curl(x -> A(x) × B(x), x) ≈ A(x)*divergence(B, x) - B(x)*divergence(A, x) + gradient(A, x)⋅B(x) - gradient(B, x)⋅A(x)
            # second derivatives
            @test divergence(x -> curl(A, x), x) ≈ 0 atol = eps(T)
            @test curl(x -> gradient(φ, x), x) ≈ zero(Vec{3}) atol = 10eps(T)
            @test divergence(x -> gradient(φ, x), x) ≈ laplace(φ, x)
            @test gradient(x -> divergence(A, x), x) ≈ curl(x -> curl(A, x), x) + laplace.(A, x)
            @test divergence(x -> ϕ(x)*gradient(φ, x), x) ≈ ϕ(x)*laplace(φ, x) + gradient(ϕ, x)⋅gradient(φ, x)
            @test divergence(x -> φ(x)*gradient(ϕ, x) - gradient(φ, x)*ϕ(x), x) ≈ φ(x)*laplace(ϕ, x) - laplace(φ, x)*ϕ(x)
            @test laplace(x -> ϕ(x)*φ(x), x) ≈ laplace(ϕ, x)*φ(x) + 2*gradient(ϕ, x)⋅gradient(φ, x) + ϕ(x)*laplace(φ, x)
            @test laplace.(x -> φ(x)*A(x), x) ≈ A(x)*laplace(φ, x) + 2*(gradient(A, x)⋅gradient(φ, x)) + φ(x)*laplace.(A, x)
            @test laplace(x -> A(x)⋅B(x), x) ≈ A(x)⋅laplace.(B, x) - B(x)⋅laplace.(A, x) + 2*divergence(x -> gradient(A, x)⋅B(x) + B(x)×curl(A, x), x)
            # third derivatives
            @test laplace.(x -> gradient(φ, x), x) ≈ gradient(x -> divergence(x -> gradient(φ, x), x), x) ≈ gradient(x -> laplace(φ, x), x)
            @test laplace(x -> divergence(A, x), x) ≈ divergence(x -> gradient(x->divergence(A, x), x), x) ≈ divergence(x -> laplace.(A, x), x)
            @test laplace.(x -> curl(A, x), x) ≈ -curl(x -> curl(x -> curl(A, x), x), x) ≈ curl(x -> laplace.(A, x), x)
        end # loop T
    end # testsection

    @testsection "f: scalar -> scalar" begin
        f(x) = 2 * x^3
        x = 3.0
        @test gradient(f, x) ≈ gradient(f, x, :all)[1] ≈ 6 * x^2
        @test gradient(f, x, :all)[2] ≈ f(x) ≈ 2 * x^3
        @test hessian(f, x) ≈ hessian(f, x, :all)[1] ≈ 12 * x
        @test hessian(f, x, :all)[2] ≈ gradient(f, x) ≈ 6 * x^2
        @test hessian(f, x, :all)[3] ≈ gradient(f, x, :all)[2] ≈ 2 * x^3
    end
    @testsection "f: scalar -> AbstractTensor" begin
        for dim in 1:3, S in (rand(Vec{dim}), rand(Tensor{2,dim}), rand(SymmetricTensor{2,dim}))
            f(x) = x^3 * S
            x = 3.0
            @test gradient(f, x) ≈ gradient(f, x, :all)[1] ≈ 3 * x^2 * S
            @test gradient(f, x, :all)[2] ≈ f(x) ≈ x^3 * S
            @test hessian(f, x) ≈ hessian(f, x, :all)[1] ≈ 6 * x * S
            @test hessian(f, x, :all)[2] ≈ gradient(f, x) ≈ 3 * x^2 * S
            @test hessian(f, x, :all)[3] ≈ f(x) ≈ x^3 * S
        end
    end
    @testsection "mixed scalar/vec" begin
        f(v::Vec, s::Number) = s * v[1] * v[2]
        v = Vec((2.0, 3.0))
        s = 4.0
        @test gradient(x -> f(v, x), s)::Float64 ≈ v[1] * v[2]
        @test gradient(x -> f(x, s), v)::Vec{2} ≈ Vec((s * v[2], s * v[1]))
        @test gradient(y -> gradient(x -> f(y, x), s), v)::Vec{2} ≈ Vec((v[2], v[1]))
        @test gradient(y -> gradient(x -> f(x, y), v), s)::Vec{2} ≈ Vec((v[2], v[1]))
    end
end # testsection

ext: jl


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    cubes, permutations = Dict{String, Int}(), Dict{String, Int}()
    for i ∈ 1:10 ^ n
        permutation = join(sort(collect(string(i ^ 3))))
        if haskey(cubes, permutation)
            permutations[permutation] += 1
            if permutations[permutation] == n
                return cubes[permutation]
            end
        else
            cubes[permutation] = i ^ 3
            permutations[permutation] = 1
        end
    end
end

compute(3)

compute(5)

@benchmark compute(5)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using WebAssemblyText
using Test

@testset "WebAssemblyText.jl" begin
    a = 1
    b = 1
    @test a == b
end


text:
#! /usr/bin/julia

# comments are allowed before YAML frontmatter

# ---
# title: Write your demo in julia
# id: juliademocard_example
# cover: assets/literate.png
# date: 2020-09-13
# author: "[Johnny Chen](https://github.com/johnnychen94)"
# julia: 1.0.1
# description: This demo shows you how to write your demo in julia
# ---

# Different from markdown format demo, source demo files are preprocessed so that it generates:
#
# 1. assets such as cover image
# 2. julia source file
# 3. mardkown file
# 4. jupyter notebook file
#
# Links to nbviewer and source files are provided as well.

# !!! note
#     This entails every source file being executed three times: 1) asset generation, 2) notebook
#     generation, and 3) Documenter HTML generation. If the generation time is too long for your
#     application, you could then choose to write your demo in markdown format. Or you can disable
#     the notebook generation via the `notebook` keyword in your demos, this helps reduce the
#     runtime to 2x.

#
# The conversions from source demo files to `.jl`, `.md` and `.ipynb` files are mainly handled by
# [`Literate.jl`](https://github.com/fredrikekre/Literate.jl). The syntax to control/filter the
# outputs can be found [here](https://fredrikekre.github.io/Literate.jl/stable/fileformat/)

x = 1//3
y = 2//5
x + y

#-

# Images can be loaded and displayed
using TestImages, ImageShow
img = testimage("lena")

# The frontmatter is also written in YAML with only a leading `#`, for example:

# ```
# # ---
# # title: <title>
# # cover: <cover>
# # id: <id>
# # date: 2020-09-13
# # author: Johnny Chen
# # julia: 1.3
# # description: <description>
# # ---
# ```

# In addition to the keywords supported by markdown files, Julia format demos accept one extra
# frontmatter keyword: `julia`. It allows you to specify the compat version of your demo. With this,
# `DemoCards` would:
#
#   - throw a warning if your demos are generated using a lower version of Julia
#   - insert a compat version badge.
#
# The warning is something like "The running Julia version `1.0.5` is older than the declared
# compatible version `1.3.0`."

# Another extra keyword is `notebook`, it allows you to control whether you needs to generate
# jupyter notebook `.ipynb` for the demo. The valid values are `true` and `false`. See also
# [Override default values with properties](@id page_sec_properties) on how to disable all notebook
# generation via the properties entry.

# !!! warning
#     You should be careful about the leading whitespaces after the first `#`. Frontmatter as weird as the
#     following is not guaranteed to work and it is very likely to hit a YAML parsing error.

#     ```yaml
#     #---
#     # title: <title>
#     #  cover: <cover>
#     #  id: <id>
#     # description: <description>
#     #---
#     ```

# !!! tip
#     Comments are allowed before frontmatter, but it would only be appeared in the julia source
#     codes. Normally, you may only want to add magic comments and license information before the
#     YAML frontmatter.


text:
using Base.Intrinsics

import Base: convert, ==, !=, .==, zero, show, length,
             getindex, setindex!, rand, string

# Residues
# ========

"""
Most of the **MIToS** design is created around the `Residue` bitstype. It represents the 20 natural amino acids and a GAP value to represent insertion, deletion but also missing data: ambiguous residues and non natural amino acids.
Each residue is encoded as an integer number, this allows fast indexing operation using Residues of probability or frequency matrices.

**Residue creation and conversion**

Creation and `convert`ion of `Residue`s should be treated carefully.
`Residue` is encoded as an 8 bits type similar to `Int8`, to get faster indexing using `Int(x::Residue)`.
In this way, `Int`, `Int8` and other signed integers returns the integer value encoded by the residue.
Conversions to and from `Char`s and `Uint8` are different, to use the `Char`acter representation in IO operations.

```julia

julia> alanine = Residue('A')
A

julia> Int(alanine)
1

julia> Char(alanine)
'A'

julia> UInt8(alanine) # 0x41 == 65 == 'A'
0x41

julia> for residue in res"ARNDCQEGHILKMFPSTWYV-"
           println(residue, " ", Int(residue))
       end
A 1
R 2
N 3
D 4
C 5
Q 6
E 7
G 8
H 9
I 10
L 11
K 12
M 13
F 14
P 15
S 16
T 17
W 18
Y 19
V 20
- 21

```
"""
bitstype 8 Residue

# Conversion from/to integers
# ---------------------------
#
# This conversions are used for indexing, comparison and operations

convert(::Type{Residue}, x::Int8) = box(Residue,unbox(Int8,x))
convert(::Type{Int8}, x::Residue) = box(Int8,unbox(Residue,x))

for typ in (Int16, Int32, Int64)
    @eval convert(::Type{$typ}, x::Residue) = convert($typ,Int8(x))
    @eval convert(::Type{Residue},  x::$typ)= convert(Residue,Int8(x))
end

Residue(x) = convert(Residue, x)

# Gaps
# ----

"""
`GAP` is the character/number representation on **MIToS** for gaps (also for non standard residues).
Lowercase characters and dots are also encoded as `GAP` in conversion from `String`s and `Char`s.
This `Residue` constant is encoded as `Residue(21)`.
"""
const GAP = Residue(21)

# Conversion from/to Char/Uint8
# -----------------------------
#
# Conversions for input and output.
# Inserts (lowercase and dots) and invalid characters will be converted into gaps ('-').

const _to_char  = Char['A','R','N','D','C','Q','E','G','H','I','L',
                       'K','M','F','P','S','T','W','Y','V','-']

convert(::Type{Char}, x::Residue) = _to_char[ Int( x ) ]

const _to_uint8 = UInt8[ UInt8(ch) for ch in  _to_char ]

convert(::Type{UInt8}, x::Residue) = _to_uint8[ Int( x ) ]

const _to_residue = fill(GAP,256)

for i in 1:length( _to_char )
    _to_residue[ Int(_to_char[i]) ] = Residue(i)
end

convert(::Type{Residue}, x::UInt8) = _to_residue[ Int(x) ];
convert(::Type{Residue}, x::Char)  = _to_residue[ Int(x) ];

# Conversion from/to ASCIIString
# ------------------------------
#
# They are useful for IO and creation of Vector{Residue}

convert(::Type{Residue}, str::ASCIIString) = Residue[ Residue(char) for char in str.data ]
convert(::Type{Residue}, str::Vector{UInt8}) = convert(Vector{Residue}, str)
convert(::Type{ASCIIString}, seq::AbstractVector{Residue}) = ASCIIString(
    UInt8[ UInt8(res) for res in seq ] )

string(seq::AbstractVector{Residue}) = ascii(seq) # "AR..." instead of the standar "[A,R,..."

# For convert a MSA stored as Vector{ASCIIString} to Matrix{Residue}
function convert(::Type{Matrix{Residue}}, sequences::Array{ASCIIString,1})
    nseq = length(sequences)
    nres = length(sequences[1])
    aln = Array(Residue,nseq,nres)
    @inbounds for i in 1:nseq
        seq = sequences[i].data
        if length(seq) == nres
            aln[i,:] = seq
        else
            throw(ErrorException(string(
            "There is an aligned sequence with different number of columns [ ",
            length(seq), " != ", nres," ]: ", ascii(seq) )))
        end
    end
    aln
end

"""
The MIToS macro `@res_str` takes a string and returns a `Vector` of `Residues` (sequence).

```julia

julia> res"MIToS"
5-element Array{MIToS.MSA.Residue,1}:
 M
 I
 T
 -
 S

```
"""
macro res_str(str)
    Residue(str)
end

# Show
# ----

show(io::IO, x::Residue) = write(io, convert(Char, x))

# Comparisons
# -----------

for fun in [:(==), :(!=), :(.==)] # , :(<), :(<=)
    @eval $(fun)(x::Residue,y::Residue) = $(fun)(Int8(x), Int8(y))
end

length(x::Residue) = length(UInt8(x))

# Random
# ------

"""
It chooses from the 20 natural residues (it doesn't generate gaps).

```julia

julia> rand(Residue)
T

julia> rand(Residue, 4, 4)
4x4 Array{MIToS.MSA.Residue,2}:
 L  E  F  L
 R  Y  L  K
 K  V  S  G
 Q  V  M  T

```
"""
rand(r, ::Type{Residue}) = Residue( rand(r, Int8(1):Int8(20)) )

# Three letters (for PDB)
# =======================

const _res2three = ASCIIString[ "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL" ]

# Thanks to Elin- for this list
const _three2res = Dict{ASCIIString, Residue}(
			                 "GLY"=>'G',
                             "ALA"=>'A',
                             "LEU"=>'L',
                             "MET"=>'M',
                             "PHE"=>'F',
                             "TRP"=>'W',
                             "LYS"=>'K',
                             "GLN"=>'Q',
                             "GLU"=>'E',
                             "SER"=>'S',
                             "PRO"=>'P',
                             "VAL"=>'V',
                             "ILE"=>'I',
                             "CYS"=>'C',
                             "TYR"=>'Y',
                             "HIS"=>'H',
                             "ARG"=>'R',
                             "ASN"=>'N',
                             "ASP"=>'D',
                             "THR"=>'T',
                             "MSE"=>'M',
                             "CSS"=>'C',
                             "2AS"=>'D',
                             "3AH"=>'H',
                             "5HP"=>'E',
                             "ACL"=>'R',
                             "AIB"=>'A',
                             "ALM"=>'A',
                             "ALO"=>'T',
                             "ALY"=>'K',
                             "ARM"=>'R',
                             "ASA"=>'D',
                             "ASB"=>'D',
                             "ASK"=>'D',
                             "ASL"=>'D',
                             "ASQ"=>'D',
                             "AYA"=>'A',
                             "BCS"=>'C',
                             "BHD"=>'D',
                             "BMT"=>'T',
                             "BNN"=>'A',
                             "BUC"=>'C',
                             "BUG"=>'L',
                             "C5C"=>'C',
                             "C6C"=>'C',
                             "CCS"=>'C',
                             "CEA"=>'C',
                             "CHG"=>'A',
                             "CLE"=>'L',
                             "CME"=>'C',
                             "CSD"=>'A',
                             "CSO"=>'C',
                             "CSP"=>'C',
                             "CSS"=>'C',
                             "CSW"=>'C',
                             "CXM"=>'M',
                             "CY1"=>'C',
                             "CY3"=>'C',
                             "CYG"=>'C',
                             "CYM"=>'C',
                             "CYQ"=>'C',
                             "DAH"=>'F',
                             "DAL"=>'A',
                             "DAR"=>'R',
                             "DAS"=>'D',
                             "DCY"=>'C',
                             "DGL"=>'E',
                             "DGN"=>'Q',
                             "DHA"=>'A',
                             "DHI"=>'H',
                             "DIL"=>'I',
                             "DIV"=>'V',
                             "DLE"=>'L',
                             "DLY"=>'K',
                             "DNP"=>'A',
                             "DPN"=>'F',
                             "DPR"=>'P',
                             "DSN"=>'S',
                             "DSP"=>'D',
                             "DTH"=>'T',
                             "DTR"=>'W',
                             "DTY"=>'Y',
                             "DVA"=>'V',
                             "EFC"=>'C',
                             "FLA"=>'A',
                             "FME"=>'M',
                             "GGL"=>'E',
                             "GLZ"=>'G',
                             "GMA"=>'E',
                             "GSC"=>'G',
                             "HAC"=>'A',
                             "HAR"=>'R',
                             "HIC"=>'H',
                             "HIP"=>'H',
                             "HMR"=>'R',
                             "HPQ"=>'F',
                             "HTR"=>'W',
                             "HYP"=>'P',
                             "IIL"=>'I',
                             "IYR"=>'Y',
                             "KCX"=>'K',
                             "LLP"=>'K',
                             "LLY"=>'K',
                             "LTR"=>'W',
                             "LYM"=>'K',
                             "LYZ"=>'K',
                             "MAA"=>'A',
                             "MEN"=>'N',
                             "MHS"=>'H',
                             "MIS"=>'S',
                             "MLE"=>'L',
                             "MPQ"=>'G',
                             "MSA"=>'G',
                             "MSE"=>'M',
                             "MVA"=>'V',
                             "NEM"=>'H',
                             "NEP"=>'H',
                             "NLE"=>'L',
                             "NLN"=>'L',
                             "NLP"=>'L',
                             "NMC"=>'G',
                             "OAS"=>'S',
                             "OCS"=>'C',
                             "OMT"=>'M',
                             "PAQ"=>'Y',
                             "PCA"=>'E',
                             "PEC"=>'C',
                             "PHI"=>'F',
                             "PHL"=>'F',
                             "PR3"=>'C',
                             "PRR"=>'A',
                             "PTR"=>'Y',
                             "SAC"=>'S',
                             "SAR"=>'G',
                             "SCH"=>'C',
                             "SCS"=>'C',
                             "SCY"=>'C',
                             "SEL"=>'S',
                             "SEP"=>'S',
                             "SET"=>'S',
                             "SHC"=>'C',
                             "SHR"=>'K',
                             "SOC"=>'C',
                             "STY"=>'Y',
                             "SVA"=>'S',
                             "TIH"=>'A',
                             "TPL"=>'W',
                             "TPO"=>'T',
                             "TPQ"=>'A',
                             "TRG"=>'K',
                             "TRO"=>'W',
                             "TYB"=>'Y',
                             "TYQ"=>'Y',
                             "TYS"=>'Y',
                             "TYY"=>'Y',
                             "AGM"=>'R',
                             "GL3"=>'G',
                             "SMC"=>'C',
                             "ASX"=>'B',
                             "CGU"=>'E',
                             "CSX"=>'C',
                             "GLX"=>'Z',
                             "UNK"=>'X'
                            )

"""
This function returns the three letter name of the `Residue`.

```julia
julia> residue2three(Residue('G'))
"GLY"

```
"""
residue2three(res::Residue) = _res2three[Int(res)]

"""
It takes a three letter residue name and returns the `Residue`.

```julia
julia> three2residue("ALA")
A

```
"""
function three2residue(res::ASCIIString)
    if length(res) == 3
        get(_three2res, uppercase(res), GAP)
    else
        throw(ErrorException("The residue name should have 3 letters."))
    end
end


text:
function igyp02a()
#
#  generalization of Igarash and Ypma
#
    igyp02(8)
end


text:
export AnEnKS

struct AnEnKS

    np::Int64

end

"""
    data_assimilation( yo, da)

Apply stochastic and sequential data assimilation technics using
model forecasting or analog forecasting.
"""
function forecast(da::DataAssimilation, yo::TimeSeries, mc::AnEnKS; progress = true)

    # dimensions
    nt = yo.nt        # number of observations
    np = mc.np        # number of particles
    nv = yo.nv        # number of variables (dimensions of problem)

    # initialization
    x̂ = TimeSeries(nt, nv)

    m_xa_part = [zeros(Float64, (nv, np)) for i = 1:nt]
    xf_part = [zeros(Float64, (nv, np)) for i = 1:nt]
    part = [zeros(Float64, (nv, np)) for i = 1:nt]
    pf = [zeros(Float64, (nv, nv)) for i = 1:nt]
    xf = zeros(Float64, (nv, np))
    m_xa_part_tmp = similar(xf)
    xf_mean = similar(xf)
    ef = similar(xf)
    Ks = zeros(Float64, (nv, nv))

    rng = MersenneTwister(1234)

    if progress
        p = Progress(nt)
    end

    sampler1 = MvNormal(da.xb, da.B)

    for k = 1:nt
        if progress
            next!(p)
        end

        # update step (compute forecasts)
        if k == 1
            rand!(rng, sampler1, xf)
        else
            xf, xf_mean = da.m(part[k-1])
            m_xa_part_tmp .= xf_mean
            m_xa_part[k] .= xf_mean
        end

        xf_part[k] .= xf

        ef .= xf * (Matrix(I, np, np) .- 1 / np)
        pf[k] .= (ef * ef') ./ (np - 1)

        # analysis step (correct forecasts with observations)
        ivar_obs = findall(.!isnan.(yo.u[k]))
        n = length(ivar_obs)

        if n > 0

            μ = zeros(Float64, n)
            σ = da.R[ivar_obs, ivar_obs]
            sampler2 = MvNormal(μ, σ)
            eps = rand(rng, sampler2, np)
            yf = da.H[ivar_obs, :] * xf

            Σ = (da.H[ivar_obs, :] * pf[k]) * da.H[ivar_obs, :]'
            Σ .+= da.R[ivar_obs, ivar_obs]

            invΣ = inv(Σ)
            K = (pf[k] * da.H[ivar_obs, :]') * invΣ
            d = yo.u[k][ivar_obs] .- yf .+ eps
            part[k] .= xf .+ K * d
            # compute likelihood
            innov_ll = vec(mean(yo.u[k][ivar_obs] .- yf, dims = 2))
            loglik = -0.5 * innov_ll'invΣ * innov_ll .- 0.5 * (n * log(2π) + log(det(Σ)))

        else

            part[k] .= xf

        end

        x̂.u[k] .= vec(sum(part[k] ./ np, dims = 2))

    end

    if progress
        p = Progress(nt)
    end

    for k = nt:-1:1
        if progress
            next!(p)
        end

        if k == nt
            part[k] .= part[nt]
        else
            m_xa_part_tmp = m_xa_part[k+1]
            tej, m_xa_tmp = da.m(mean(part[k], dims = 2))
            tmp1 = part[k] .- mean(part[k], dims = 2)
            tmp2 = m_xa_part_tmp .- m_xa_tmp
            Ks .= ((tmp1 * tmp2') * pinv(pf[k+1], rtol = 1e-6)) ./ (np - 1)
            part[k] .+= Ks * (part[k+1] .- xf_part[k+1])
        end
        x̂.u[k] .= vec(sum(part[k] ./ np, dims = 2))

    end

    x̂

end


text:
function axis_aligned_bounding_box(	particle_type::String,
									R::Array{Float64, 1},
									x::Float64,
									y::Float64,
									z::Float64,
									q0::Float64,
									q1::Float64,
									q2::Float64,
									q3::Float64)
	lbx::Float64 = 0.0
	ubx::Float64 = 0.0
	lby::Float64 = 0.0
	uby::Float64 = 0.0
	lbz::Float64 = 0.0
	ubz::Float64 = 0.0

	a11::Float64 = 0.0
	a12::Float64 = 0.0
	a13::Float64 = 0.0
	a21::Float64 = 0.0
	a22::Float64 = 0.0
	a23::Float64 = 0.0
	a31::Float64 = 0.0
	a32::Float64 = 0.0
	a33::Float64 = 0.0

	if particle_type == "sphere"
		lbx = x - R[1]
		ubx = x + R[1]
		lby = y - R[1]
		uby = y + R[1]
		lbz = z - R[1]
		ubz = z + R[1]
	elseif particle_type == "ellipse"
		# Rotation matrix.
		(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(q0, q1, q2, q3)

		# Unrotated corner coordinates of the ellipse (effectively 2-D) AABB. We don't subtract the 1 from the z coordinates because we only use four corners now.
		C = 2.0 * [[0.0, 0.0, 0.0] [0.0, 1.0, 0.0] [1.0, 0.0, 0.0] [1.0, 1.0, 0.0]]
		C[1, :] += - 1.0
		C[2, :] += - 1.0

		# Unrotated but correctly scaled corner coordinates of the ellipse AABB.
		for i = 1:2
			C[i, :] *= R[i]
		end

		# Rotated corner coordinates, yielding the ellipse OBB.
		for i = 1:8
			C[:, i] = [a11 a12 a13 ; a21 a22 a23 ; a31 a32 a33] * C[:, i]
		end

		# Lower and upper bounds of OBB.
		lbx = x + minimum(C[1, :])
		ubx = x + maximum(C[1, :])
		lby = y + minimum(C[2, :])
		uby = y + maximum(C[2, :])
		lbz = z + minimum(C[3, :])
		ubz = z + maximum(C[3, :])
	elseif particle_type == "ellipsoid"
		# Rotation matrix.
		(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(q0, q1, q2, q3)

		# Unrotated corner coordinates of the ellipsoids AABB.
		C = 2.0 * [[0.0, 0.0, 0.0] [0.0, 0.0, 1.0] [0.0, 1.0, 0.0] [0.0, 1.0, 1.0] [1.0, 0.0, 0.0] [1.0, 0.0, 1.0] [1.0, 1.0, 0.0] [1.0, 1.0, 1.0]] - 1.0

		# Unrotated but correctly scaled corner coordinates of the ellipsoids AABB.
		for i = 1:3
			C[i, :] *= R[i]
		end

		# Rotated corner coordinates, yielding the ellipsoids OBB.
		for i = 1:8
			C[:, i] = [a11 a12 a13 ; a21 a22 a23 ; a31 a32 a33] * C[:, i]
		end

		# Lower and upper bounds of OBB.
		lbx = x + minimum(C[1, :])
		ubx = x + maximum(C[1, :])
		lby = y + minimum(C[2, :])
		uby = y + maximum(C[2, :])
		lbz = z + minimum(C[3, :])
		ubz = z + maximum(C[3, :])
	elseif particle_type == "cuboid"
		# Rotation matrix.
		(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(q0, q1, q2, q3)

		# Unrotated corner coordinates.
		C = 2.0 * [[0.0, 0.0, 0.0] [0.0, 0.0, 1.0] [0.0, 1.0, 0.0] [0.0, 1.0, 1.0] [1.0, 0.0, 0.0] [1.0, 0.0, 1.0] [1.0, 1.0, 0.0] [1.0, 1.0, 1.0]] - 1.0

		# Unrotated but correctly scaled corner coordinates.
		for i = 1:3
			C[i, :] *= R[i]
		end

		# Rotated corner coordinates.
		for i = 1:8
			C[:, i] = [a11 a12 a13 ; a21 a22 a23 ; a31 a32 a33] * C[:, i]
		end

		# Lower and upper bounds of AABB.
		lbx = x + minimum(C[1, :])
		ubx = x + maximum(C[1, :])
		lby = y + minimum(C[2, :])
		uby = y + maximum(C[2, :])
		lbz = z + minimum(C[3, :])
		ubz = z + maximum(C[3, :])
	end

	return (lbx, ubx, lby, uby, lbz, ubz)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    aitchison(c₁, c₂)

Return the Aitchison distance between compositions `c₁` and `c₂`.
"""
aitchison(c₁::Composition, c₂::Composition) = norm(c₁ - c₂)

"""
    Aitchison()

Aitchison distance following the Distances.jl API.
"""
struct Aitchison <: Metric end

(::Aitchison)(c₁::Composition, c₂::Composition) = aitchison(c₁, c₂)

(d::Aitchison)(w₁, w₂) = d(Composition(w₁), Composition(w₂))

result_type(::Aitchison, w₁, w₂) = Float64


text:
export BilinearLens

@doc doc"""

    BilinearLens(ϕ)

`BilinearLens` is a lensing operator that computes lensing with
bilinear interpolation. The action of the operator, as well as its
adjoint, inverse, inverse-adjoint, and gradient w.r.t. `ϕ` can all be
computed. The log-determinant of the operation is non-zero and can't
be computed.

Internally, `BilinearLens` forms a sparse matrix with the
interpolation weights, which can be applied and adjoint-ed extremely
fast (e.g. at least an order of magnitude faster than
[`LenseFlow`](@ref)). Inverse and inverse-adjoint lensing is somewhat
slower since it requires an iterative solve, here performed with the
[preconditioned generalized minimal
residual](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method)
algorithm.

"""
mutable struct BilinearLens{T,Φ<:Field{<:Any,T},S} <: ImplicitOp{T}
    ϕ :: Φ
    sparse_repr :: S
    anti_lensing_sparse_repr :: Union{S, Nothing}
end

function BilinearLens(ϕ::FlatField)

    # if ϕ == 0 then just return identity operator
    if norm(ϕ) == 0
        return BilinearLens(ϕ,I,I)
    end

    @unpack Nbatch,Nx,Ny,Δx = ϕ
    T = real(ϕ.T)
    Nbatch > 1 && error("BilinearLens with batched ϕ not implemented yet.")

    # the (i,j)-th pixel is deflected to (ĩs[i],j̃s[j])
    j̃s,ĩs = getindex.((∇*ϕ)./Δx, :Ix)
    ĩs .=  ĩs  .+ (1:Ny)
    j̃s .= (j̃s' .+ (1:Nx))'

    # sub2ind converts a 2D index to 1D index, including wrapping at edges
    indexwrap(i,N) = mod(i - 1, N) + 1
    sub2ind(i,j) = Base._sub2ind((Ny,Nx),indexwrap(i,Ny),indexwrap(j,Nx))

    # compute the 4 non-zero entries in L[I,:] (ie the Ith row of the sparse
    # lensing representation, L) and add these to the sparse constructor
    # matrices, M, and V, accordingly. this function is split off so it can be
    # called directly or used as a CUDA kernel
    function compute_row!(I, ĩ, j̃, M, V)

        # (i,j) indices of the 4 nearest neighbors
        left,right = floor(Int,ĩ) .+ (0, 1)
        top,bottom = floor(Int,j̃) .+ (0, 1)

        # 1-D indices of the 4 nearest neighbors
        M[4I-3:4I] .= @SVector[sub2ind(left,top), sub2ind(right,top), sub2ind(left,bottom), sub2ind(right,bottom)]

        # weights of these neighbors in the bilinear interpolation
        Δx⁻, Δx⁺ = ((left,right) .- ĩ)
        Δy⁻, Δy⁺ = ((top,bottom) .- j̃)
        A = @SMatrix[
            1 Δx⁻ Δy⁻ Δx⁻*Δy⁻;
            1 Δx⁺ Δy⁻ Δx⁺*Δy⁻;
            1 Δx⁻ Δy⁺ Δx⁻*Δy⁺;
            1 Δx⁺ Δy⁺ Δx⁺*Δy⁺
        ]
        V[4I-3:4I] .= inv(A)[1,:]

    end

    # a surprisingly large fraction of the computation for large Nside, so memoize it:
    @memoize getK(Nx,Ny) = Int32.((4:4*Nx*Ny+3) .÷ 4)

    # CPU
    function compute_sparse_repr(is_gpu_backed::Val{false})
        K = Vector{Int32}(getK(Nx,Ny))
        M = similar(K)
        V = similar(K,T)
        for I in 1:length(ĩs)
            compute_row!(I, ĩs[I], j̃s[I], M, V)
        end
        sparse(K,M,V,Nx*Ny,Nx*Ny)
    end

    # GPU
    function compute_sparse_repr(is_gpu_backed::Val{true})
        K = CuVector{Cint}(getK(Nx,Ny))
        M = similar(K)
        V = similar(K,T)
        cuda(ĩs, j̃s, M, V; threads=256) do ĩs, j̃s, M, V
            index = threadIdx().x
            stride = blockDim().x
            for I in index:stride:length(ĩs)
                compute_row!(I, ĩs[I], j̃s[I], M, V)
            end
        end
        CuSparseMatrixCSR(CuSparseMatrixCOO{T}(K,M,V,(Nx*Ny,Nx*Ny)))
    end

    BilinearLens(ϕ, compute_sparse_repr(Val(is_gpu_backed(ϕ))), nothing)

end

# lazily computing the sparse representation for anti-lensing

function get_anti_lensing_sparse_repr!(Lϕ::BilinearLens)
    if Lϕ.anti_lensing_sparse_repr == nothing
        Lϕ.anti_lensing_sparse_repr = BilinearLens(-Lϕ.ϕ).sparse_repr
    end
    Lϕ.anti_lensing_sparse_repr
end

getϕ(Lϕ::BilinearLens) = Lϕ.ϕ
(Lϕ::BilinearLens)(ϕ::FlatField) = BilinearLens(ϕ)
hash(L::BilinearLens, h::UInt64) = foldr(hash, (typeof(L), getϕ(L)), init=h)

# applying various forms of the operator

function *(Lϕ::BilinearLens, f::FlatField)
    Lϕ.sparse_repr===I && return f
    Łf = Ł(f)
    f̃ = similar(Łf)
    for batch in 1:size(f.arr,4), pol in 1:size(f.arr,3)
        mul!(@views(f̃.arr[:,:,pol,batch][:]), Lϕ.sparse_repr, @views(Łf.arr[:,:,pol,batch][:]))
    end
    f̃
end

function *(Lϕ::Adjoint{<:Any,<:BilinearLens}, f::FlatField)
    parent(Lϕ).sparse_repr===I && return f
    Łf = Ł(f)
    f̃ = similar(Łf)
    for batch in 1:size(f.arr,4), pol in 1:size(f.arr,3)
        mul!(@views(f̃.arr[:,:,pol,batch][:]), parent(Lϕ).sparse_repr', @views(Łf.arr[:,:,pol,batch][:]))
    end
    f̃
end

function \(Lϕ::BilinearLens, f̃::FlatField)
    Lϕ.sparse_repr===I && return f̃
    Łf̃ = Ł(f̃)
    f = similar(Łf̃)
    for batch in 1:size(f.arr,4), pol in 1:size(f.arr,3)
        @views(f.arr[:,:,pol,batch][:]) .= gmres(
            Lϕ.sparse_repr, @views(Łf̃.arr[:,:,pol,batch][:]),
            Pl = get_anti_lensing_sparse_repr!(Lϕ), maxiter = 5
        )
    end
    f
end

function \(Lϕ::Adjoint{<:Any,<:BilinearLens}, f̃::FlatField)
    parent(Lϕ).sparse_repr===I && return f̃
    Łf̃ = Ł(f̃)
    f = similar(Łf̃)
    for batch in 1:size(f.arr,4), pol in 1:size(f.arr,3)
        @views(f.arr[:,:,pol,batch][:]) .= gmres(
            parent(Lϕ).sparse_repr', @views(Łf̃.arr[:,:,pol,batch][:]),
            Pl = get_anti_lensing_sparse_repr!(parent(Lϕ))', maxiter = 5
        )
    end
    f
end

for op in (:*, :\)
    @eval function ($op)(Lϕ::Union{BilinearLens, Adjoint{<:Any,<:BilinearLens}}, f::FieldTuple)
        FieldTuple(map(f->($op)(Lϕ,f), f.fs))
    end
end

# gradients

@adjoint BilinearLens(ϕ) = BilinearLens(ϕ), Δ -> (Δ,)

@adjoint function *(Lϕ::BilinearLens, f::Field{B}) where {B}
    f̃ = Lϕ * f
    function back(Δ)
        (∇' * (Ref(spin_adjoint(Ł(Δ))) .* Ł(∇*f̃))), B(Lϕ*Δ)
    end
    f̃, back
end

# gpu

adapt_structure(storage, Lϕ::BilinearLens) = BilinearLens(adapt(storage, fieldvalues(Lϕ))...)


text:
export uniform,
       uniform_left,
       uniform_right,
       Roberts_left,
       Roberts_right,
       Roberts_both,
       cluster

function uniform_direction(hstart::FT, Δh::FT, N::IT, dir::IT) where {FT, IT}
  @assert (N>0)
  hn = [hstart+FT(dir)*FT(i-1)*Δh for i in 1:N+1]
  hn = dir>0 ? hn : reverse(hn)
  return hn
end

function assert_monotonically_increasing!(h::VFT, hmin::FT, hmax::FT, β::FT, caller::AbstractString) where {FT,VFT<:AbstractVector{FT}}
  s = length(h)
  for i in 1:s-2
    R = (h[i+2]-h[i+1])/(h[i+1]-h[i])
    if R < FT(1)
      println("Error: coordinates not monotonically increasing in ", caller)
      println("hmin, hmax, β = ", hmin, hmax, β)
      println("R = ", R)
      println("h = ", h)
      throw(ArgumentError("coordinates not monotonically increasing"))
    end
  end
end

function insist_monotonically_decreasing!(h::VFT, hmin::FT, hmax::FT, β::FT, caller::AbstractString) where {FT,VFT<:AbstractVector{FT}}
  s = length(h)
  for i in 1:s-2
    R = (h[i+2]-h[i+1])/(h[i+1]-h[i])
    if R > FT(1)
     println("Error: coordinates not monotonically decreasing in ", caller)
     println("in stretch_funcs.jl")
     println("hmin, hmax, β = ", hmin, hmax, β)
     println("R = ", R)
     println("h = ", h)
     throw(ArgumentError("coordinates not monotonically increasing"))
    end
  end
end

"""
    transformation1(hmin::FT, hmax::FT, N::FT, β::FT) where FT

The coordinates and differences of a Robert's stretching
function as described in section 5.6 (page 333) of
Computational Fluid Mechanics and Heat Transfer,
2nd edition, J. Tannehill et al. (Transformation 1)

# Input
```
    hmin     = minimum value
    hmax     = maximum value
    N        = N segments of Δh
    β     = 1.0 <= β <= infinity = stretching factor
             = larger -> less stretching
  y=0                         y=h
                               |
   |-|--|---|-------|----------|
   |------> y
```

!!! note
    I have abused notation a bit to provide consistent notation
    with the reference as well as generalize the returned grid
    so that it need not start at y=0.
"""
function transformation1(hmin::FT, hmax::FT, N::IT, β::FT) where {FT, IT}
  a = FT(1); b = β
  g = (b+FT(1))/(b-FT(1))
  Δh = (hmax - hmin)/FT(N)
  hnbar = FT[hmin+FT(i-1)*Δh for i in 1:N+1]
  hnbar = (hnbar .- hmin)./(hmax-hmin)
  hn = [((b+a)-(b-a)*g^(a-hnbar[i]))/(g^(a-hnbar[i])+a) for i in 1:N+1]
  assert_monotonically_increasing!(hn, hmin, hmax, β, "transformation1")
  hn = hmin .+ (hmax - hmin).*hn
  return hn
end

"""
    transformation2

Coordinates and differences of a Robert's
stretching function as described in section 5.6 (page 333) of
Computational Fluid Mechanics and Heat Transfer,
2nd edition, J. Tannehill et al.

# Input
```
hmin     = minimum value
hmax     = maximum value
N        = N segments of Δh
α        = 0      stretching at y=h only
α        = 0.5    stretching at y=0 and y=hmax
β        = 1.0 <= β <= infinity = stretching factor
         = larger -> less stretching
```

Here is a picture illustration for α = 0:
```
                             y=h
                              |
  |----------|-------|---|--|-|
  |------> y
```

Note that this must be used in reverse for the lid driven
cavity geometry for the 'front' and 'back' walls.

!!! note
    I have abused notation a bit to provide consistent notation
    with the reference as well as generalize the returned grid
    so that it need not start at y=0.
"""
function transformation2(hmin::FT, hmax::FT, N::IT, α::FT, β::FT) where {FT, IT}
  a = α; b = β
  g = (b+FT(1))/(b-FT(1))
  Δh = (hmax - hmin)/FT(N)
  hnbar = FT[hmin+FT(i-1)*Δh for i in 1:N+1]
  hnbar = (hnbar .- hmin)./(hmax-hmin)
  hn = [ ((b+FT(2)*a)*g^((hnbar[i]-a)/(FT(1)-a))-
          b+FT(2)*a)/((FT(2)*a+FT(1))*(FT(1)+
          g^((hnbar[i]-a)/(FT(1)-a)))) for i in 1:N+1]
  if abs(α)<eps(FT)
   insist_monotonically_decreasing!(hn, hmin, hmax, β, "transformation2")
  end
  hn .= hmin .+ (hmax - hmin).*hn
  return hn
end

"""
    transformation3

This function returns the coordinates and differences of a Robert's
stretching function as described in section 5.6 (page 333) of
Computational Fluid Mechanics and Heat Transfer,
2nd edition, J. Tannehill et al. (Transformation 3)

# Input
```
hmin     = minimum value
hmax     = maximum value
N        = N segments of Δh
τ        = 0 <= τ <= infinity = stretching factor

τ        = 0            --> no stretching
τ        = large values --> strong stretching
```

Here is a picture illustration for α = 0:

```
                                 y=yc                        y=h
                                  |
      |----------|-------|---|--|-|-|--|---|-------|----------|
      |------> y
```

Note that this must be used in reverse for the lid driven
cavity geometry for the 'front' and 'back' walls.

!!! note
    Notation has been abused a bit to provide consistent notation
    with the reference as well as generalize the returned grid
    so that it need not start at y=0.
"""
function transformation3(hmin::FT, hmax::FT, N::IT, yc::FT, τ::FT) where {FT, IT}
  a = FT(1); c = FT(2)
  Δh = (hmax - hmin)/FT(N)
  hnbar = [hmin+FT(i-1)*Δh for i in 1:N+1]
  hnbar = (hnbar .- hmin) ./ (hmax-hmin)
  d = yc/(hmax-hmin); e = FT(exp(τ))
  B = a/(c*τ)*log((a+(e-a)*d)/(a+(a/e-a)*d))
  hn = [d*(a+FT(sinh(τ*(hnbar[i]-B)))/FT(sinh(τ*B))) for i in 1:N+1]
  hn .= hmin .+ (hmax - hmin).*hn
  return hn
end

uniform_left(hstart::FT, Δh::FT, N::IT) where {FT, IT} = uniform_direction(hstart, Δh, N,  IT(-1))
uniform_right(hstart::FT, Δh::FT, N::IT) where {FT, IT} = uniform_direction(hstart, Δh, N, IT( 1))

uniform(hmin::FT, hmax::FT, N::IT) where {FT, IT} = [hmin+FT(i-1)*(hmax - hmin)/FT(N) for i in 1:N+1]
Roberts_left(hmin::FT, hmax::FT, N::IT, β::FT) where {FT, IT} = transformation1(hmin, hmax, N, β)
Roberts_right(hmin::FT, hmax::FT, N::IT, β::FT) where {FT, IT} = transformation2(hmin, hmax, N, FT(0), β)
Roberts_both(hmin::FT, hmax::FT, N::IT, β::FT) where {FT, IT} = transformation2(hmin, hmax, N, FT(1/2), β)
cluster(hmin::FT, hmax::FT, N::IT, yc::FT, tau::FT) where {FT, IT} = transformation3(hmin, hmax, N, yc, tau)