text:
#  ┌─┬┐
#  ├─┼┤
#  │ ││
#  └─┴┘

function compactstring(x, l = 16)
    str = "$x"
    if length(str) > l
        if l >= 10
            l1 = div(l,2)
            l2 = l - l1 - 1
            return str[1:l1] * "…" * str[end-l2+1:end]
        else
            return str[1:l] * "…"
        end
    else
        return str
    end
end

#compactstring(x::Integer, l) = "$x"

function compactstring(x::Nullable, l = 16)
    if isnull(x)
        #if l == 1
            return "-"
        #else
        #    return "NA"
        #end
    else
        return compactstring(get(x), l)
    end
end

function compactstring(x::Bool, l = 16)
    if l < 5
        x ? (return "T") : (return "F")
        #x ? (return "✓") : (return "✗")
    else
        return "$x"
    end
end

# TODO think hard about this one. Some data is numerical, while some is a label
# Is there a simple way of the user controlling the difference?
function compactstring(x::Integer,l=8)
    if x < 10_000
        return string(x)
    elseif x < 10_000_000
        tmp = div(x + 500, 1000)
        return "$(tmp) k"
    elseif x < 10_000_000_000
        tmp = div(x + 500_000, 1_000_000)
        return "$(tmp) M"
    elseif x < 10_000_000_000_000
        tmp = div(x + 500_000_000, 1_000_000_000)
        return "$(tmp) B"
    else
        tmp = div(x + 500_000_000_000, 1_000_000_000_000)
        return "$(tmp) T"
    end
end

function compactstring(x::Float64,l=16)
    if isnan(x)
        return "NaN"
    elseif isinf(x)
        if x > 0
            return "Inf"
        else
            return "-Inf"
        end
    elseif x == 0.0e0
        return "0.0000"
    else
        if x < 0
            isneg = true
            x = -x
        else
            isneg = false
        end
        if x <= 1e-3 # If it is too small
            exp = 0
            factor = 1e3
            while true
                exp -= 3
                x *= factor
                if x > 1.0
                    break
                end
            end

            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)

            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            #return "$(whole).$(frac_str)e$(exp)"
            return "$("-"^isneg)$(whole).$(frac_str) × 10$(superscript(exp))"

        elseif x >= 1e4 # Or if it is too large
            exp = 0
            factor = 1e-3
            while true
                exp += 3
                x *= factor
                if x < 1e3
                    break
                end
            end

            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)
            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            #return "$(whole).$(frac_str)e$(exp)"
            return "$("-"^isneg)$(whole).$(frac_str) × 10$(superscript(exp))"

        else # Otherwise it is "reasonably" sized
            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)
            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            return "$("-"^isneg)$(whole).$(frac_str)"
        end
    end
end

function compactstring(x::Float32,l=16)
    if isnan(x)
        return "NaN"
    elseif isinf(x)
        if x > 0
            return "Inf"
        else
            return "-Inf"
        end
    elseif x == 0.0f0
        return "0.0000"
    else
        if x < 0
            isneg = true
            x = -x
        else
            isneg = false
        end
        if x <= 1f-3 # If it is too small
            exp = 0
            factor = 1f3
            while true
                exp -= 3
                x *= factor
                if x > 1.0
                    break
                end
            end

            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)

            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            #return "$(whole).$(frac_str)e$(exp)"
            return "$("-"^isneg)$(whole).$(frac_str) × 10$(superscript(exp))"

        elseif x >= 1f4 # Or if it is too large
            exp = 0
            factor = 1f-3
            while true
                exp += 3
                x *= factor
                if x < 1f3
                    break
                end
            end

            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)
            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            #return "$(whole).$(frac_str)e$(exp)"
            return "$("-"^isneg)$(whole).$(frac_str) × 10$(superscript(exp))"

        else # Otherwise it is "reasonably" sized
            tmp = Int32(round(x*10_000))
            (whole,frac) = divrem(tmp, 10_000)
            frac_str = string(frac)
            if frac < 10
                frac_str = "000" * frac_str
            elseif frac < 100
                frac_str = "00" * frac_str
            elseif frac < 1000
                frac_str = "0" * frac_str
            end

            return "$("-"^isneg)$(whole).$(frac_str)"
        end
    end
end

function superscript(x::Integer)
    normalstring = string(x)
    superscriptstring = string()
    for char in normalstring
        superscriptstring = superscriptstring * string(superscript(char))
    end
    return superscriptstring
end

function superscript(char::Char)
    if char == '-'
        return '⁻'
    elseif char == '+'
        return '⁺'
    elseif char == '0'
        return '⁰'
    elseif char == '1'
        return '¹'
    elseif char == '2'
        return '²'
    elseif char == '3'
        return '³'
    elseif char == '4'
        return '⁴'
    elseif char == '5'
        return '⁵'
    elseif char == '6'
        return '⁶'
    elseif char == '7'
        return '⁷'
    elseif char == '8'
        return '⁸'
    elseif char == '9'
        return '⁹'
    end
end

#=
function compactstring(x::Union{Float64,Float32},l=16)
    local str
    if l <= 5
        str = @sprintf("%0.3g",x)
    elseif l == 6
        str = @sprintf("%0.4g",x)
    elseif l == 7
        str = @sprintf("%0.5g",x)
    elseif l == 8
        str = @sprintf("%0.6g",x)
    elseif l == 9
        str = @sprintf("%0.7g",x)
    elseif l == 10
        str = @sprintf("%0.8g",x)
    elseif l == 1
        str = @sprintf("%0.9g",x)
    elseif l == 12
        str = @sprintf("%0.10g",x)
    elseif l == 13
        str = @sprintf("%0.11g",x)
    elseif l == 14
        str = @sprintf("%0.12g",x)
    elseif l == 15
        str = @sprintf("%0.13g",x)
    elseif l == 16
        str = @sprintf("%0.14g",x)
    else
        str = @sprintf("%0.15g",x)
    end
    if str[end] == ' '
        str = str[1:end-1]
    end
    if search(str,'.') == 0 # No decimal point but its a float...
        if search(str,'e') == 0 && search(str,"Inf") == 0:-1 && search(str,"NaN") == 0:-1 # just in case!
            str = str * "."
        end
    end
    return str
end
=#

function compactstring(x::AbstractString,l=16)
    if length(x) > l-2
        return "\"$(x[1:l-3])…\""
    else
        return "\"$(x)\""
    end
end

align_strings(T,x) = x

function align_strings{T <: Integer, S <: AbstractString}(::Union{Type{T},Type{Nullable{T}}}, x::Vector{S})
    isempty(x) && return x
    ls = map(length,x)
    l = maximum(ls)
    for i = 1:length(x)
        x[i] = (" " ^ (l-ls[i])) * x[i]
    end
end

function align_strings{T <: AbstractFloat, S <: AbstractString}(::Union{Type{T},Type{Nullable{T}}}, x::Vector{S})
    isempty(x) && return x

    decimalpoint = [search(x[i],'.') for i = 1:length(x)]
    #if minimum(decimalpoint) == 0
    #    return x # Missing the decimal point somewhere... give up in confusion
    #end
    m = maximum(decimalpoint)

    for i = 1:length(x)
        x[i] = (" " ^ (m-decimalpoint[i])) * x[i]
    end
end

_displaysize(io) = haskey(io, :displaysize) ? io[:displaysize] : _displaysize(io.io)

function printtable(io::IO, header_str::Vector{UTF8String}, data_str::Vector{Vector{UTF8String}}, row_header_str::UTF8String = UTF8String(""), row_str::Vector{UTF8String} = (length(data_str) > 0 ? fill(UTF8String(""), length(data_str[1])) : Vector{UTF8String}(""));
                    dotsbreak = -1, # Insert some dots after this row (set to 0 for at top)
                    max_width = 255,
                    tl = "┌─",
                    t = "─", # length 1
                    tr = "─┐",
                    tsep = "─┬─",
                    hl = "├─",
                    h = "─", # length 1
                    hsep = "─┼─",
                    hr = "─┤",
                    l = "│ ",
                    sep = " │ ",
                    r = " │",
                    bl = "└─",
                    b = "─", # length 1
                    br = "─┘",
                    bsep = "─┴─",
                    pad = " ", # length 1
                    vdots = "⋮",# length 1
                    hdots = "…") # length 1)

    ncols = length(data_str)
    nrows = length(data_str) > 0 ? length(data_str[1]) : 0
    header_str_bak = copy(header_str)

    row_width = max(length(row_header_str), maximum(isempty(row_str) ? 0 : map(length, row_str)))
    widths = [max(length(header_str[i]), maximum(isempty(data_str[i]) ? 0 : map(length, data_str[i]))) for i = 1:ncols]

    n_skipped_cols = 0
    too_wide = row_width + sum(widths) + 2 + 3*length(widths) > max_width
    was_too_wide = too_wide
    while too_wide
        if length(widths) == 1
            break # Show at least one column of data, even if it is ugly
        end
        n_skipped_cols += 1
        pop!(widths)
        pop!(header_str)
        pop!(data_str)
        too_wide = row_width + sum(widths) + 6 + 3*length(widths) > max_width
        ncols = ncols - 1
    end
    if was_too_wide
        push!(widths, 1)
        push!(header_str,hdots)
        push!(data_str,fill(hdots, length(data_str[1])))
        ncols = ncols + 1
    end

    # Top line
    print(io,pad^(row_width+1))
    for c = 1:ncols
        if c == 1
            print(io,tl)
        else
            print(io,tsep)
        end
        print(io, t ^ widths[c])
        if c == ncols
            println(io,tr)
        end
    end

    # Field names
    print(io,row_header_str * (pad ^ (row_width - length(row_header_str) + 1)))
    for c = 1:ncols
        if c == 1
            print(io,l)
        else
            print(io,sep)
        end
        print(io, header_str[c])
        if length(header_str[c]) < widths[c]
            print(io, pad ^ (widths[c] - length(header_str[c])))
        end
        if c == ncols
            println(io,r)
        end
    end

    # Header seperator
    print(io,pad^(row_width+1))
    for c = 1:ncols
        if c == 1
            print(io,hl)
        else
            print(io,hsep)
        end
        print(io, h ^ widths[c])
        if c == ncols
            println(io,hr)
        end
    end

    # Special case for vdots at start of table (e.g. might be great for displaying the tail of a table)
    if dotsbreak == 0
        print(io,pad^(row_width+1))
        for c = 1:ncols
            if c == 1
                print(io,l)
            else
                print(io,sep)
            end
            print(io, pad ^ div(widths[c]-1,2))
            print(io, vdots)
            print(io, pad ^ (widths[c]-1-div(widths[c]-1,2)))
            if c == ncols
                println(io,r)
            end
        end
    end

    # Data
    for i = 1:nrows
        print(io, pad ^ (row_width - length(row_str[i])) * row_str[i] * pad)
        for c = 1:ncols
            if c == 1
                print(io,l)
            else
                print(io,sep)
            end
            print(io, data_str[c][i])
            if length(data_str[c][i]) < widths[c]
                print(io, pad ^ (widths[c] - length(data_str[c][i])))
            end
            if c == ncols
                println(io,r)
            end
        end

        if i == dotsbreak
            print(io,pad^(row_width+1))
            for c = 1:ncols
                if c == 1
                    print(io,l)
                else
                    print(io,sep)
                end
                print(io, pad ^ div(widths[c]-1,2))
                print(io, vdots)
                print(io, pad ^ (widths[c]-1-div(widths[c]-1,2)))
                if c == ncols
                    println(io,r)
                end
            end
        end
    end

    # Bottom line
    print(io,pad^(row_width+1))
    for c = 1:ncols
        if c == 1
            print(io,bl)
        else
            print(io,bsep)
        end
        print(io, b ^ widths[c])
        if c == ncols
            print(io,br)
        end
    end
    if n_skipped_cols > 0
        print(io, "\n" * pad^(row_width) * "+ $n_skipped_cols unshown columns: ")
        for i = length(header_str_bak)+1-n_skipped_cols:length(header_str_bak)
            print(io,"$(header_str_bak[i])")
            if i < length(header_str_bak)
                print(io, ", ")
            end
        end
    end
end

function make_strings(data; maxl::Int = 5, width_suggestion::Int = 20)
    data_str = Vector{UTF8String}()
    nrows = length(data)

    if nrows > 2*maxl
        for i = 1:maxl
            push!(data_str, compactstring(data[i],width_suggestion))
        end
        for i = (nrows-maxl+1):nrows
            push!(data_str, compactstring(data[i],width_suggestion))
        end
    else
        for i = 1:nrows
            push!(data_str, compactstring(data[i],width_suggestion))
        end
    end

    align_strings(eltype(data), data_str)

    return data_str
end

function Base.show(io::IO,table::Table)
    if ncol(table) == 0
        print(io,"Empty Table")
        return
    else
        ncols = ncol(table)
        nrows = length(table)

        println(io, "$nrows-row × $ncols-column Table:")
    end

    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5
    maxl = max(5,div(s[1],5)) # Maximum number of lines to show (head, then tail)

    # First we format all of our output and determine its size
    # Rows (on side)
    row_header_str = UTF8String("Row")
    if nrows > 2*maxl
        row_str = UTF8String[string(i) for i = vcat(1:maxl, nrows-maxl+1:nrows)]
        dotsbreak = maxl
    else
        row_str = UTF8String[string(i) for i = 1:nrows]
        dotsbreak = -1
    end

    # header....
    col_names = colnames(table)
    header_str = [UTF8String(string(col_names[i])) for i = 1:ncols]

    data_str = [Vector{UTF8String}() for i = 1:ncols]
    for c = 1:ncols
        width_suggestion = 20
        if eltype(table.(c)) <: Union{Bool,Nullable{Bool},Float64,Float32} && length(header_str[c]) < width_suggestion
            width_suggestion = length(header_str[c])
        end

        data_str[c] = make_strings(table.(c); maxl = maxl, width_suggestion=width_suggestion)
    end

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str, row_header_str, row_str; max_width = s[2],
                                                                  dotsbreak = dotsbreak,
                                                                  tl = "╔═",
                                                                  t = "═",
                                                                  tr = "═╗",
                                                                  tsep = "═╤═",
                                                                  l = "║ ",
                                                                  r = " ║",
                                                                  hl = "╟─",
                                                                  hr = "─╢",
                                                                  bl = "╚═",
                                                                  b = "═",
                                                                  br = "═╝",
                                                                  bsep = "═╧═")
end

function Base.showall(io::IO,table::Table)
    if ncol(table) == 0
        print(io,"Empty Table")
        return
    else
        ncols = ncol(table)
        nrows = length(table)

        println(io, "$nrows-row × $ncols-column Table:")
    end

    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5
    maxl = max(5,div(s[1],5)) # Maximum number of lines to show (head, then tail)

    # First we format all of our output and determine its size
    # Rows (on side)
    row_header_str = UTF8String("Row")
    row_str = UTF8String[string(i) for i = 1:nrows]

    # header....
    col_names = colnames(table)
    header_str = [UTF8String(string(col_names[i])) for i = 1:ncols]

    data_str = [Vector{UTF8String}() for i = 1:ncols]
    for c = 1:ncols
        width_suggestion = 20
        if eltype(table.(c)) <: Union{Bool,Nullable{Bool},Float64,Float32} && length(header_str[c]) < width_suggestion
            width_suggestion = length(header_str[c])
        end

        data_str[c] = make_strings(table.(c); maxl = nrows, width_suggestion=width_suggestion)
    end

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str, row_header_str, row_str; max_width = s[2],
                                                                  tl = "╔═",
                                                                  t = "═",
                                                                  tr = "═╗",
                                                                  tsep = "═╤═",
                                                                  l = "║ ",
                                                                  r = " ║",
                                                                  hl = "╟─",
                                                                  hr = "─╢",
                                                                  bl = "╚═",
                                                                  b = "═",
                                                                  br = "═╝",
                                                                  bsep = "═╧═")
end

function Base.show(io::IO, col::Column)
    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5
    maxl = max(5,div(s[1],5)) # Maximum number of lines to show (head, then tail)
    nrows = length(col)

    println(io, "$nrows-row Column:")

    # First we format all of our output and determine its size
    # Rows (on side)
    row_header_str = UTF8String("Row")
    if nrows > 2*maxl
        row_str = UTF8String[string(i) for i = vcat(1:maxl, nrows-maxl+1:nrows)]
        dotsbreak = maxl
    else
        row_str = UTF8String[string(i) for i = 1:nrows]
        dotsbreak = -1
    end

    # header....
    header_str = [UTF8String(string(colname(col)))]

    data_str = Vector{UTF8String}[make_strings(col.(1); maxl = maxl, width_suggestion=60),]

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str, row_header_str, row_str; max_width = s[2],
                                                                  dotsbreak = dotsbreak,
                                                                  tl = "╒═",
                                                                  t = "═",
                                                                  tr = "═╕",
                                                                  tsep = "═╤═",
                                                                  bl = "╘═",
                                                                  b = "═",
                                                                  br = "═╛",
                                                                  bsep = "═╧═")
end

function Base.showall(io::IO, col::Column)
    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5
    nrows = length(col)

    println(io, "$nrows-element Column:")

    # First we format all of our output and determine its size
    # Rows (on side)
    row_header_str = UTF8String("Row")
    row_str = UTF8String[string(i) for i = 1:nrows]

    # header....
    header_str = [UTF8String(string(colname(col)))]

    data_str = Vector{UTF8String}[make_strings(col.(1); maxl = nrows, width_suggestion=60),]

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str, row_header_str, row_str; max_width = s[2],
                                                                  tl = "╒═",
                                                                  t = "═",
                                                                  tr = "═╕",
                                                                  tsep = "═╤═",
                                                                  bl = "╘═",
                                                                  b = "═",
                                                                  br = "═╛",
                                                                  bsep = "═╧═")
end

function Base.show(io::IO,row::Row)
    if ncol(row) == 0
        print(io,"Empty Row")
        return
    end

    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5
    ncols = ncol(row)

    println(io, "$ncols-element Row:")

    # First we format all of our output and determine its size

    # header....
    col_names = colnames(row)
    header_str = [UTF8String(string(col_names[i])) for i = 1:ncols]

    data_str = [Vector{UTF8String}() for i = 1:ncols]
    for c = 1:ncols
        width_suggestion = 20
        if eltype(row.(c)) <: Union{Bool,Nullable{Bool},Float64,Float32} && length(header_str[c]) < width_suggestion
            width_suggestion = length(header_str[c])
        end

        data_str[c] = make_strings([row.(c)]; width_suggestion=width_suggestion)
    end

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str; max_width = s[2],
                                         tl = "╓─",
                                         tr = "─╖",
                                         l = "║ ",
                                         r = " ║",
                                         hl = "╟─",
                                         hr = "─╢",
                                         bl = "╙─",
                                         br = "─╜")
end

showall(io::IO, row::Row) = show(io, row)

function Base.show(io::IO, cell::Cell)
    println(io, "Cell:")

    s = Base.tty_size() # [height, width] in characters TODO fix for Julia 0.5

    # First we format all of our output and determine its size

    # header....
    header_str = [UTF8String(string(colname(cell)))]

    data_str = Vector{UTF8String}[make_strings([cell.(1)]; width_suggestion=60),]

    # Now we show the table using computed widths for decorations
    printtable(io, header_str, data_str; max_width = s[2])
end

showall(io::IO, cell::Cell) = show(io, cell)


text:
const TU = typeof(1.0f0s)
const TI = typeof(1.0f0s .. 2.0f0s)
const default_fs = 1000.0f0Hz

struct TrialData
    trialrange::UnitRange{Int}
    soundinterval::TI
    cuetime::TU
    odorinterval::TI
    odordirection::Int
    lickonsets::Vector{TU}
end
TrialData(tr::AbstractUnitRange, si::AbstractInterval, ct::Quantity, oi::AbstractInterval, od) = TrialData(tr, si, ct, oi, od, TU[])

Base.:(==)(td1::TrialData, td2::TrialData) = td1.trialrange == td2.trialrange &&
                                             td1.soundinterval == td2.soundinterval &&
                                             td1.cuetime == td2.cuetime &&
                                             td1.odorinterval == td2.odorinterval &&
                                             td1.odordirection == td2.odordirection &&
                                             td1.lickonsets == td2.lickonsets

function Base.show(io::IO, td::TrialData)
    println(io, "Trial over indices ", td.trialrange, ':')
    println(io, "  sound: ", td.soundinterval)
    println(io, "  cue: ", td.cuetime)
    println(io, "  odor delivery: ", td.odorinterval)
    println(io, "  odor direction: ", td.odordirection)
    print(io, "  lick times: ")
    Base.show_delim_array(io, td.lickonsets, '[', ",", ']', false)  # skips printing of unit
end

struct DataFrameRate
    df::DataFrame
    fs::typeof(default_fs)
end

function Base.show(io::IO, dfr::DataFrameRate)
    show(io, dfr.df)
    print(io, "\nRecorded at ", dfr.fs)
end

function Base.getproperty(dfr::DataFrameRate, name::Symbol)
    name === :fs && return getfield(dfr, :fs)
    df = getfield(dfr, :df)
    name === :df && return df
    return getproperty(df, name)
end

ext: jl


text:
module RODClf

using Images

export img_load,
    integral_sum,
    imgobj_extract,
    ImageFilter,
    ImageObject,
    FilterShape,
    FilterGeometry

type ImageObject
  properties::Dict
  image::Array
  ndims::Tuple
end

type FilterGeometry
  dx_resolution::Signed #total number of 'dx' units for image
  dy_resolution::Signed #total number of 'dy' units for image
  dx_units::Signed #total 'dx' units needed
  dy_units::Signed #total 'dy' units needed
  y_units::Signed #number of panels on 'y' axis
  x_units::Signed #number of panels on 'x'axis
  y_offset::Signed #offset-coefficient for y-axis
  x_offset::Signed #offset-coefficient for x-axis
end

function img_load(imgpath::String)
  f = separate(imread(imgpath)::Image)
  imgobj = ImageObject(f.properties, f.data, size(f.data))::ImageObject
  return imgobj
end

function integral_sum(image_array)
  integral_image = 0x00::Uint8
  for ui in image_array
    integral_image += ui
    continue
  end
  return integral_image
end

function ImageFilter(filter_shape::FilterGeometry, imgobj::ImageObject)
    dx_resolution = filter_shape.dx_resolution
    dy_resolution = filter_shape.dy_resolution
    dx_units = filter_shape.dx_units
    dy_units = filter_shape.dy_units
    y_units = filter_shape.x_units
    x_units = filter_shape.y_units
    dx = signed(imgobj.ndims[1]/dx_resolution-0.5)
    dy = signed(imgobj.ndims[2]/dy_resolution-0.5);
    y_offset = filter_shape.y_offset
    x_offset = filter_shape.x_offset

    x = Array[]
    for i in [1:x_units]
        for j in [1:y_units]
            push!(x, imgobj.image[((i+y_offset)*dx)-dx+1:(i+y_offset)*dx,((j+x_offset)*dy)-dy+1:(j+x_offset)*dy,:])
        end
    end
    return x
end

function imgobj_extract(imgobj::ImageObject)
  shape = imgobj.ndims::Tuple
  properties = imgobj.properties::Dict
  dx, dy = signed(shape[2]/8)::Signed, signed(shape[3]/8)::Signed
  testA = imgobj.image[:,dx:dx*2, dy:dy*8]::Array
  testB = imgobj.image[:,dx*2:dx*3, dy:dy*8]::Array
  return testA, testB
end

end

max_stars_repo_name: torrange/RODClf

max_issues_repo_name: torrange/RODClf

max_forks_repo_name: torrange/RODClf


text:
export visualize!

"""
    visualize!(simulation; save_as, fps, frame_nums, frame_size, xlim, ylim, particle_colors, multithreaded)

# Arguments
- `simulation::Simulation`: simulation to visualize
- `save_as::Tuple{Symbol, String}` = (:gif, "temp.gif"):
- `fps::Int64 = 10`: frames per second if creating an animation
- `frame_nums::Union{Array{Int64, 1}, Symbol} = :all`: frame numbers to use.
- `frame_size::Tuple{Int64, Int64} = (600, 600)`: size of frame in pixels
- `xlim::Union{Array{Float64, 1}, Symbol} = :all`: frame limits along x-axis
- `ylim::Union{Array{Float64, 1}, Symbol} = :all`: frame limits along y-axis
- `particle_colors::Union{DefaultDict{Symbol, String, String}, Dict{Symbol, String}} = DefaultDict{Symbol, String}("black")`: particle colors
- `multithreaded::Bool = false`: use threads to update frames
"""
function visualize!(simulation::Simulation; save_as::Tuple{Symbol, String} = (:gif, "temp.gif"),
                                            fps::Int64 = 10,
                                            frame_nums::Union{Array{Int64, 1}, Symbol} = :all,
                                            frame_size::NTuple{2, Int64} = (600, 600),
                                            xlim::Union{Array{Float64, 1}, Symbol} = :all, ylim::Union{Array{Float64, 1}, Symbol} = :all,
                                            particle_colors::Union{DefaultDict{Symbol, String, String}, Dict{Symbol, String}} = DefaultDict{Symbol, String}("black"),
                                            multithreaded::Bool = false)
    save_type, save_str = save_as
    if !isdir(dirname(save_str))
        mkpath(dirname(save_str))
    end

    if frame_nums == :all
        frame_nums = [f for f = 1 : length(simulation.history)]
    elseif frame_nums == :start
        frame_nums = [1]
    elseif frame_nums == :end
        frame_nums = [length(simulation.history)]
    elseif frame_nums == :startend
        frame_nums = [1, length(simulation.history)]
    end

    if xlim == :all
        xlim = [0.0, simulation.L_x]
    end
    if ylim == :all
        ylim = [0.0, simulation.L_y]
    end

    if save_type == :gif
        print_message("GENERATING GIF")
        animation = Animation()
    elseif save_type == :images
        print_message("GENERATING IMAGES")
    end

    N_particles = length(simulation.things_to_save.particles)
    N_bonds = length(simulation.things_to_save.bonds)
    scene = plot(size = frame_size, legend = false, axis = false, grid = false, aspect_ratio = 1, xlim = xlim, ylim = ylim)
    for (f, frame_num) = enumerate(frame_nums)
        if f == 1
            θ = LinRange(0.0, 2 * pi, 20)
            unit_circle = (xs = cos.(θ), ys = sin.(θ))
            for particle in simulation.history[frame_num].particles
                cxs, cys = particle.R .* unit_circle.xs .+ particle.x, particle.R .* unit_circle.ys .+ particle.y
                color = particle_colors[particle.ptype]
                plot!(cxs, cys, seriestype = [:shape,], color = color, linecolor = color, fillalpha = 0.3)
            end

            for (particle_1, particle_2) in simulation.history[frame_num].bonds
                x_1, y_1 = particle_1.x, particle_1.y
                x_2, y_2 = particle_2.x, particle_2.y
                Δx = wrap_displacement(x_1 - x_2; period = simulation.periodic_in_x ? simulation.L_x : -1.0)
                Δy = wrap_displacement(y_1 - y_2; period = simulation.periodic_in_y ? simulation.L_y : -1.0)
                plot!([x_2, x_2 + Δx], [y_2, y_2 + Δy], color = "black", linealpha = 0.3)
                plot!([x_1, x_1 - Δx], [y_1, y_1 - Δy], color = "black", linealpha = 0.3)
            end
        else
            prev_frame_num = frame_nums[f - 1]
            for n = 1 : N_particles
                prev_particle, particle = simulation.history[prev_frame_num].particles[n], simulation.history[frame_num].particles[n]
                Δx, Δy = particle.x - prev_particle.x, particle.y - prev_particle.y
                scene.series_list[n][:x] .+= Δx
                scene.series_list[n][:y] .+= Δy
            end
            for n = N_particles + 1 : 2 : N_particles + 2 * N_bonds
                particle_1, particle_2 = simulation.history[frame_num].bonds[div(n - N_particles - 1, 2) + 1]
                x_1, y_1 = particle_1.x, particle_1.y
                x_2, y_2 = particle_2.x, particle_2.y
                Δx = wrap_displacement(x_1 - x_2; period = simulation.periodic_in_x ? simulation.L_x : -1.0)
                Δy = wrap_displacement(y_1 - y_2; period = simulation.periodic_in_y ? simulation.L_y : -1.0)
                scene.series_list[n][:x], scene.series_list[n][:y] = [x_2, x_2 + Δx], [y_2, y_2 + Δy]
                scene.series_list[n + 1][:x], scene.series_list[n + 1][:y] = [x_1, x_1 - Δx], [y_1, y_1 - Δy]
            end
        end

        if save_type == :gif
            frame(animation)
        elseif save_type == :images
            savefig(string(save_str, "Frame ", frame_num, ".png"))
        end
        println("Frame ", frame_num, " done")
    end

    if save_type == :gif
        gif(animation, save_str, fps = fps)
        print_message("GIF GENERATED")
    elseif save_type == :images
        print_message("IMAGES GENERATED")
    end
end


text:
module Cobweb

using DefaultApplication: DefaultApplication
using Scratch: @get_scratch!

#-----------------------------------------------------------------------------# init
struct CobwebDisplay <: AbstractDisplay end

function __init__()
    global DIR = @get_scratch!("CobWeb")
    pushdisplay(CobwebDisplay())
end

#-----------------------------------------------------------------------------# Node
struct Node
    tag::String
    attrs::Dict{String,String}
    children::Vector
    function Node(tag, attrs, children)
        new(tag, attrs, children)
    end
end

function Base.:(==)(a::Node, b::Node)
    a.tag == b.tag &&
        a.attrs == b.attrs &&
        length(a.children) == length(b.children) &&
        all(ac == bc for (ac,bc) in zip(a.children, b.children))
end

function Base.getproperty(node::Node, class::String)
    node.attrs["class"] = class
    node
end

get_attrs(kw) = Dict(string(k) => string(v) for (k,v) in kw)

(node::Node)(children...; kw...) = Node(node.tag, merge(node.attrs, get_attrs(kw)), vcat(node.children, children...))

#-----------------------------------------------------------------------------# h
h(tag, children...; kw...) = Node(tag, get_attrs(kw), collect(children))

h(tag, attrs::Dict, children...) = Node(tag, attrs, collect(children))

function Base.getproperty(::typeof(h), tag::Symbol)
    f(children...; kw...) = h(String(tag), children...; kw...)
    return f
end

# TODO: Something smarter than flat-out replacing symbols
macro h(ex)
    esc(_h(ex))
end

function _h(ex::Expr)
    for i in 1:length(ex.args)
        x = ex.args[i]
        if x isa Expr
            ex.args[i] = _h(x)
        elseif x isa Symbol
            ex.args[i] = Expr(:., :(Cobweb.h), QuoteNode(ex.args[1]))
        end
    end
    ex
end

#-----------------------------------------------------------------------------# escapeHTML
# Taken from HTTPCommon.jl (ref: http://stackoverflow.com/a/7382028/3822752)
function escape_html(x::String)
    s = replace(x,  '&' => "&")
    s = replace(s, '"' => """)
    s = replace(s, ''' => "'")
    s = replace(s, '<' => "<")
    replace(s, '>' => ">")
end

#-----------------------------------------------------------------------------# show (html)
pretty(x) = (io = IOBuffer(); pretty(io, x); String(take!(io)))

pretty(io::IO, node::Node) = show(IOContext(io, :pretty=>true), node)

function Base.show(io::IO, node::Node)
    color = get(io, :tagcolor, 1)
    pretty = get(io, :pretty, false)
    level = pretty ? get(io, :level, 0) : 0
    indent = pretty ? ' ' ^ get(io, :indent, 2) : ""
    p(args...) = printstyled(io, args...; color)
    p(indent ^ level, '<', node.tag)
    for (k,v) in node.attrs
        if v == "true"
            p(' ', k)
        elseif v == "false"
        else
            p(' ', k, '=', '"', v, '"')
        end
    end
    p('>')
    has_node_children = any(x -> x isa Node, node.children)
    pretty && has_node_children && p('\n')
    n_nodes = 0
    for (i, child) in enumerate(node.children)
        if child isa String
            has_node_children ?
                p(indent ^ (level + 1), escape_html(child), '\n') :
                p(escape_html(child))
        elseif child isa Node
            show(IOContext(io, :tagcolor => color + i + n_nodes, :level => level+1), MIME("text/html"), child)
            n_nodes = sum(x -> x isa Node, child.children, init=0)
        else
            show(io, MIME("text/html"), child)
        end
    end
    if has_node_children
        p(indent ^ level, "</", node.tag, '>')
    else
        p("</", node.tag, '>')
    end
    pretty && println(io)
end
Base.show(io::IO, ::MIME"text/html", node::Node) = show(io, node)
Base.show(io::IO, ::MIME"text/xml", node::Node) = show(io, node)
Base.show(io::IO, ::MIME"application/xml", node::Node) = show(io, node)
Base.show(io::IO, ::MIME"text/xml", node::Node) = show(io, node)

#-----------------------------------------------------------------------------# show (javascript)
struct Javascript
    x::String
end
Base.show(io::IO, ::MIME"text/javascript", j::Javascript) = print(io, j.x)
function Base.show(io::IO, ::MIME"text/html", j::Javascript)
    print(io, "<script>")
    print(io, j.x)
    print(io, "</script>")
end

#-----------------------------------------------------------------------------# CSS
struct CSS
    content::Dict{String, Dict{String,String}}
    function CSS(o::AbstractDict)
        new(Dict(string(k) => Dict(string(k2) => string(v2) for (k2,v2) in pairs(v)) for (k,v) in pairs(o)))
    end
end
function Base.show(io::IO, o::CSS)
    for (k,v) in o.content
        println(io, k, " {")
        for (k2, v2) in v
            println(io, "  ", k2, ':', ' ', v2, ';')
        end
        println(io, '}')
    end
end
Base.show(io::IO, ::MIME"text/css", o::CSS) = print(io, o)
function Base.show(io::IO, ::MIME"text/html", o::CSS)
    println(io, "<style>")
    print(io, o)
    println(io, "</style>")
end
save(file::String, o::CSS) = save(o, file)
save(o::CSS, file::String) = open(io -> show(io, x), touch(file), "w")

#-----------------------------------------------------------------------------# Page
struct Page
    content
end
Page(pg::Page) = pg

save(file::String, page::Page) = save(page, file)

function save(page::Page, file=joinpath(DIR, "index.html"))
    Base.open(touch(file), "w") do io
        println(io, "<!doctype html>")
        show(io, MIME("text/html"), page.content)
    end
    file
end

Base.display(::CobwebDisplay, page::Page) = DefaultApplication.open(save(page))

end #module


text:
"""
    non_relativistic_orbital(n, ℓ[, Z = 1])

Returns a callable object to evaluate the _reduced_ radial orbital ``P_{nℓ}(r)`` of a
hydrogen-like atom in non-relativistic theory, given by

```math
P_{nℓ}(r) =
\\sqrt{
    \\left(\\frac{2Z}{n}\\right)^3
    \\frac{(n-\\ell-1)!}{2n(n+\\ell)!}
}
\\left(\\frac{2Z}{n} \\cdot r\\right)^\\ell
L^{(2\\ell+1)}_{n-\\ell-1}\\left(\\frac{2Z}{n} \\cdot r\\right)
\\exp\\left(-\\frac{Z}{n} \\cdot r\\right)
```

``P_{nℓ}(r)`` is defined through

```math
\\psi_{nlm}(r,\\theta,\\varphi) = \\frac{1}{r} P_{nℓ}(r) Y_{ℓm}(\\theta,\\varphi)
```

where ``\\psi_{nlm}`` are the normalized solutions to the non-relativistic single-particle
Schrödinger equation for a point-like nucleus with charge `Z`, and where ``Y_{ℓm}`` are the
spherical harmonics.
"""
function non_relativistic_orbital(n::Integer, ℓ::Integer, Z = 1)
    n > 0 || throw(ArgumentError("Invalid n = $n, must be n ≥ 1"))
    0 ≤ ℓ < n || throw(ArgumentError("Invalid ℓ = $ℓ for n = $n, must be 0 ≤ ℓ < n"))
    # Normalization factor.
    # Note: n! = Γ(n+1), so (n-ℓ-1)!/(2n(n+ℓ)!) = Γ(n+1)/(2nΓ(n+ℓ+1))
    N = sqrt((2Z/n)^3*gamma(n-ℓ)/(2n*gamma(n+ℓ+1)))
    # Generalized Laguerre polynomial L^{(2ℓ+1)}_{n-ℓ-1}(q), where q = 2Zr/n:
    L = convert(Polynomial, Laguerre{2ℓ+1}([zeros(n-ℓ-1); 1], :q))
    # Single Polynomial object for the polynomial part, multiplied by r to get the
    # reduced polynomial:
    C = let r = Polynomial([0, 1], :r)
        N * (2*Z*r/n)^ℓ * L(2*Z*r/n) * r
    end
    # Multiply with the exponential exp(Zr/n) and return a callable object:
    q = -Z/n # pre-compute the coefficient in the exponential
    r -> C(r) * exp(q*r)
end

export non_relativistic_orbital


text:
# Does this return a function? I hope so

#== Lark architecture (Python)

For each grammar rule, Lark creates a chain of
transformations associated with that rule. These are kept in a Callback
object, with each rule indexed by an internal name and
returning the combined function object. In order to
make it possible to configure, but not execute, the functions,
they are created as callable objects which are initialised
during callback creation.

==#

ExpandSingleChild(node_builder::Function) = function(c)
    if length(c) == 1
        return c[1]
    else
        return node_builder(c)
    end
end

PropagatePositions(node_builder::Function) = function(children)
    res = node_builder(children)
    #println("Res is $res")
    if res isa Tree
        meta(res).empty = true
        for c in children
            #println("Child: $c")
            if c isa Tree && c.children != [] && !(meta(c).empty)
                #println("Setting values from $c: $(meta(c))")
                meta(res).line = meta(c).line
                meta(res).column = meta(c).column
                meta(res).start_pos = meta(c).start_pos
                meta(res).empty = false
                break
            elseif c isa Token
                #println("Setting token values from $c")
                meta(res).line = c.line
                meta(res).column = c.column
                meta(res).start_pos = c.pos_in_stream
                meta(res).empty = false
                break
            end
        end

        for c in reverse(children)
            if c isa Tree && c.children != [] && !(meta(c).empty)
                #println("Setting end values from $c to $(meta(c))")
                meta(res).end_line = meta(c).end_line
                meta(res).end_column = meta(c).end_column
                meta(res).end_pos = meta(c).end_pos
                meta(res).empty = false
                break
            elseif c isa Token
                #println("Setting token end values from $c")
                meta(res).end_line = c.end_line
                meta(res).end_column = c.end_column
                meta(res).end_pos = c.pos_in_stream + length(c.value)
            break
            end
        end
    end
    return res
end

# Note: we explicitly give the type of `filtered`, otherwise
# it will adopt the type of c[i].children, which may be
# a subarray and therefore not allow later push! operations.

ChildFilter(to_include,append_none,node_builder) = function(c)
    filtered = Union{Token,Tree,Nothing}[]
    for (i,to_expand,add_none) in to_include
        if add_none > 0
            append!(filtered,fill(nothing,add_none))
        end
        if to_expand
            append!(filtered,c[i].children)
        else
            push!(filtered,c[i])
        end
    end

    if append_none > 0
        append!(filtered,fill(nothing,append_none))
    end
    node_builder(filtered)
end

ChildFilterLALR(to_include,node_builder) = function(c)
    #println("Filtering children: $to_include")
    filtered = Union{Token,Tree}[]
    for (i, to_expand) in to_include
            if to_expand
                append!(filtered, c[i].children)
            else
                push!(filtered,c[i])
            end
    end
    return node_builder(filtered)
end

ChildFilterLALR_NoPlaceholders(to_include,node_builder) = function(c)
    filtered = Union{Token,Tree}[]
    for (i, to_expand) in to_include
        if to_expand
            if length(filtered) > 0
                append!(filtered,c[i].children)
            else
                filtered = c[i].children
            end
        else
            push!(filtered,c[i])
        end
    end
    return node_builder(filtered)
end

_should_expand(sym) = !is_terminal(sym) && first(sym.name) == '_'

# TODO: possible bug from 0-indexing in python going to 1-indexing in Julia
maybe_create_child_filter(expansion,keep_all_tokens,ambiguous,_empty_indices) = begin
    if _empty_indices !== nothing
        @assert count(x->!x,_empty_indices) == length(expansion)
        s = join([x-> x ? "1" : "0" , _empty_indices],"")
        empty_indices = [length(ones) for ones in split(s,"0")]
        @assert length(empty_indices) == length(expansion) + 1
    else
        empty_indices = fill(0,length(expansion)+1)
    end
    to_include = []
    nones_to_add = 0
    for (i,sym) in enumerate(expansion)
        nones_to_add += empty_indices[i]
        if keep_all_tokens || !(is_terminal(sym) && sym.filter_out)
            push!(to_include, (i,_should_expand(sym), nones_to_add))
            nones_to_add = 0
        end
    end
    nones_to_add += empty_indices[length(expansion)+1]

    if _empty_indices !== nothing || length(to_include) < length(expansion) || any([to_expand for (i,to_expand) in to_include])
        if _empty_indices !== nothing || ambiguous
            return partial(ambiguous ? ChildFilter : ChildFilterLALR, to_include, nones_to_add)
        else
            return partial(ChildFilterLALR_NoPlaceholders, [(i,x) for (i,x,_) in to_include])
        end
    end
    return false
end

AmbiguousExpander(to_expand,node_builder) = function(children)
    _is_ambig_tree(t) = hasproperty(t,:data) && t.data == "_ambig"
    ambiguous = [i for i in to_expand if _is_ambig_tree(children[i])]
    if !isempty(ambiguous)
        expand = map(enumerate(children)) do (i,child)
            if i in ambiguous child.children else repeated(child) end
        end
        return Tree("_ambig",[node_builder(collect(f[1])) for f in product(zip(expand...))])
    end
    return node_builder(children)
end

maybe_create_ambiguous_expander(expansion,keep_all_tokens) = begin
    to_expand = [i for (i,sym) in enumerate(expansion)
                 if keep_all_tokens || ((!(is_terminal(sym) && sym.filter_out)) && _should_expand(sym))]
    if length(to_expand)>0
        return partial(AmbiguousExpander, to_expand)
    else
        false
    end
end

"""
    Propagate ambiguous intermediate nodes and their derivations up to the
    current rule.

    In general, converts

    rule
      _iambig
        _inter
          someChildren1
          ...
        _inter
          someChildren2
          ...
      someChildren3
      ...

    to

    _ambig
      rule
        someChildren1
        ...
        someChildren3
        ...
      rule
        someChildren2
        ...
        someChildren3
        ...
      rule
        childrenFromNestedIambigs
        ...
        someChildren3
        ...
      ...

    propagating up any nested '_iambig' nodes along the way.
    """
AmbiguousIntermediateExpander(node_builder,children) = begin
    _is_iambig_tree(child) = begin
        return hasproperty(child, :data) && child.data == "_iambig"
    end

    _collapse_iambig(children) = begin
        #==
            Recursively flatten the derivations of the parent of an '_iambig'
            node. Returns a list of '_inter' nodes guaranteed not
            to contain any nested '_iambig' nodes, or None if children does
            not contain an '_iambig' node.
        ==#

        # Due to the structure of the SPPF,
        # an '_iambig' node can only appear as the first child
        if !isempty(children) && _is_iambig_tree(children[1])
            iambig_node = children[1]
            result = []
            for grandchild in iambig_node.children
                collapsed = _collapse_iambig(grandchild.children)
                if !isempty(collapsed)
                    for child in collapsed
                        append!(child.children, children[2:end])
                        append!(result, collapsed)
                    end
                else
                    new_tree = Tree("_inter", append!(copy(grandchild.children), children[2:end]))
                    push!(result,new_tree)
                end
            end
            return result
        end
    end

    collapsed = _collapse_iambig(children)
    if !isempty(collapsed)
        processed_nodes = [node_builder(c.children) for c in collapsed]
        return Tree("_ambig", processed_nodes)
    end
    return node_builder(children)
end

abstract type Callback end

# No need for the @wrap decorator from Python as that is just
# for function naming and doc strings
ptb_inline_args(func) = begin
    return children -> func(children...)
    end

inplace_transformer(func) = begin
    return children -> func(Tree(Symbol(func),children))
end

apply_visit_wrapper(func,name,wrapper) = begin
    return children -> wrapper(func,name,children,nothing)
end

struct ParseTreeBuilder
    rule_builders::Array{Tuple{Rule,Array{Function}},1}
end

ParseTreeBuilder(rules::Array{Rule,1};propagate_positions=false,ambiguous=false,maybe_placeholders=false)= begin
    rule_builders = _init_builders(rules,ambiguous,propagate_positions,maybe_placeholders)
    ParseTreeBuilder(rule_builders)
end

# The wrapper chain provides a set of operations to perform based on options, very compact

## Note the Python behaviour of 'and' and 'or'; they return the most recent argument
# that has been evaluated

_init_builders(rules,ambiguous,propagate_positions,maybe_placeholders) = begin
    build_list = []
    for rule in rules
        options = rule.options
        keep_all_tokens = options !== nothing ? options.keep_all_tokens : false
        expand_single_child = options.expand1
        wrapper_chain =
            filter(x-> typeof(x) != Bool, [
                if (expand_single_child && rule.alias == nothing) ExpandSingleChild else false end,
                maybe_create_child_filter(rule.expansion, keep_all_tokens, ambiguous, maybe_placeholders ? options.empty_indices : nothing),
                if propagate_positions PropagatePositions else false end,
                if ambiguous
                maybe_create_ambiguous_expander(rule.expansion,keep_all_tokens)
                else false
                end,
                if ambiguous AmbiguousIntermediateExpander else false end
            ])
        push!(build_list, (rule,wrapper_chain))
    end
    return build_list
end

# These callbacks are executed during parsing, as each rule is
# matched (i.e. the reduction step). If a transformer has been
# specified, the appropriate method is called; by default a
# parse tree is constructed using Tree, which can later be
# visited by transformers or visitors.

# Note the use of partial to store partial information about
# the appropriate transformer to call, for later application.
# We have to supply a dummy "nothing" argument for meta.

# TODO: integrate transformer callbacks into the callback chain
# while retaining type stability. Probably requires transformers
# to be user-typed

create_callback(ptb::ParseTreeBuilder;transformer=nothing) = begin
    callbacks = Dict()
    i = 0
    for (rule,wrapper_chain) in ptb.rule_builders
        #println("Rule: $(rule.origin)")
        # Now find the actual transformer function ...
        user_callback_name =
            if rule.alias !== nothing rule.alias
            elseif rule.options.template_source !== nothing rule.options.template_source
            else rule.origin.name
            end
        f = partial(transformer_func,transformer,Val{Symbol(user_callback_name)}(),Meta())
        if typeof(transformer) <: Transformer_InPlace
            f = inplace_transformer(f)
        end

        for w in wrapper_chain
            if w == nothing
                println("Rule $rule has a nothing in the wrapper!!")
                @assert w!=nothing
            end
            f = w(f)  #ExpandSingleChild, maybe_create_child_filter,propagate positions
        end
        if rule in keys(callbacks)
            throw(GrammarError("Rule '$rule' already exists"))
        end
        callbacks[rule]=f
    end
    return callbacks
end


text:
function ssdDistance(referenceImage::regImage,templateImage::regImage,
                     transformedGrid::scaledArray;
                     doDerivative::Bool=false,doHessian::Bool=false,options::regOptions=regOptions(),centeredGrid::Array{Float64,1}=zeros(1))

    return ssdDistance(referenceImage,templateImage,
        transformedGrid.data;
        doDerivative=doDerivative,doHessian=doHessian,options=options,centeredGrid=centeredGrid)

end

function ssdDistance(referenceImage::regImage,templateImage::regImage,
                     transformedGrid::Array{Float64,1};
                     doDerivative::Bool=false,doHessian::Bool=false,options::regOptions=regOptions(),centeredGrid::Array{Float64,1}=zeros(1), interpolationScheme=BSpline(Cubic(Line())))

  # get relevant options
  parametricOnly = options.parametricOnly
  matrixFree = options.matrixFree

  # check centered grid
  if(size(centeredGrid)[1]==1)
    centeredGrid = getCellCenteredGrid(referenceImage).data
  end

  # check grid
  if(checkForOddNumberOfGridPoints(referenceImage.data,transformedGrid))
    @printf("Dimension mismatch: 2* image data: %s (%s) <-> transformed grid: %s", size(referenceImage.data),size(referenceImage.data[:]), size(transformedGrid))
    error("StaggeredGrids are not supported. Please transform to cell centered first.")
  end

  # interpolation of the template image at transformed grid points
  transformedImage, dX_transformedImage, dY_transformedImage =
      interpolateImage(templateImage,transformedGrid,doDerivative=true, interpolationScheme=interpolationScheme)

  # measure the ssd distance
  N = prod(size(referenceImage.data)) # product of sizes

  residual = Array(transformedImage .- referenceImage.data)[:]
  prodH = prod(referenceImage.voxelsize)
  functionValue = 0.5 * prodH * residual' * residual
  # calculate ssd derivatives matrix free?
  if(doDerivative)
      if(matrixFree)
          dFunctionValue,d2FunctionValue = ssdDerivativesMatrixFree(dX_transformedImage,dY_transformedImage,residual,centeredGrid,prodH,N,parametricOnly,doDerivative,doHessian)
      else
          dFunctionValue,d2FunctionValue = ssdDerivativesMatrixBased(dX_transformedImage,dY_transformedImage,residual,centeredGrid,prodH,N,parametricOnly,doDerivative,doHessian)
      end
  else
      dFunctionValue,d2FunctionValue = 0,0
  end

  return [functionValue,dFunctionValue,d2FunctionValue,(dX_transformedImage,dY_transformedImage)]

end

function ssdDerivativesMatrixBased(dX_transformedImage::Array{Float64,1},
                                   dY_transformedImage::Array{Float64,1},
                                   residual::Array{Float64,1},
                                   centeredGrid::Array{Float64,1},
                                   pixelSize::Number,N::Number,
                                   parametricOnly::Bool,doDerivative::Bool,doHessian::Bool)

  if(doDerivative)
    dTransformedImage = spdiagm((dX_transformedImage,dY_transformedImage),[0,N])
    if(parametricOnly)
      Q = sparse([centeredGrid[1:N] centeredGrid[N+1:end] ones(N)])
      Q = [Q spzeros(size(Q)[1],size(Q)[2]); spzeros(size(Q)[1],size(Q)[2]) Q]
      dTransformedImage = dTransformedImage * Q
    end
    dFunctionValue = pixelSize .* dTransformedImage' * residual
  else
    dFunctionValue = 0.0
    dTransformedImage = 0.0
  end

  if(doHessian)
    d2FunctionValue = pixelSize .* dTransformedImage' * dTransformedImage
  else
    d2FunctionValue = 0.0
  end

  return dFunctionValue, d2FunctionValue

end

function ssdDerivativesMatrixFree(dX_transformedImage::Array{Float64,1},
                                  dY_transformedImage::Array{Float64,1},
                                  residual::Array{Float64,1},
                                  centeredGrid::Array{Float64,1},
                                  pixelSize::Number,N::Number,
                                  parametricOnly::Bool,doDerivative::Bool,doHessian::Bool)

  if(doDerivative)
    dFunctionValue = BLAS.scal(2*N,pixelSize,dTtransposedMultiplication(N,dX_transformedImage,dY_transformedImage,residual),1)
    if(parametricOnly)
      dFunctionValue = QtransposedMultiplication(N,dFunctionValue,centeredGrid)
    end
  else
    dFunctionValue = 0.0
    dTransformedImage = 0.0
  end

  if(doHessian)
    if(parametricOnly)
      d2FunctionValue = computeParametricHessian(N,pixelSize,dX_transformedImage,dY_transformedImage,centeredGrid)
    else
      d2FunctionValue(x) = hessianFunction(N,pixelSize,dX_transformedImage,dY_transformedImage,x)
    end
  else
    d2FunctionValue = 0.0
  end

  return dFunctionValue,d2FunctionValue

end

function sampleArrayToGridWithSmoothing(a::scaledArray, newGrid::scaledArray)
    # @assert (a.dimensions[1] > newGrid.dimensions[1] &&
    #     a.dimensions[2] > newGrid.dimensions[2]) "new Grid has to be coarser then current grid"
    width = newGrid.voxelsize[1]/a.voxelsize[1]
    #print(width)
    adata = reshape(a.data, a.dimensions)
    # adata = smoothArray(adata,3,width)
    adata = smoothArray(adata,3, 3.0, kernel=1/3*ones(3))
    #a_int = interpolateArray(0.5*adata[1:1:end-1, 1:1:end-1] + 0.5*adata[2:1:end, 2:1:end], a.voxelsize, a.shift, newGrid.data, newGrid.dimensions)[1]
    a_int = interpolateArray(adata, a.voxelsize, a.shift, newGrid.data, newGrid.dimensions)[1]

    return a_int
end

function ssdDistanceArbitraryGrid(referenceImage::regImage,templateImage::regImage,
                     transformedGrid::scaledArray;
                     doDerivative::Bool=false, doHessian::Bool=false,
                     options::regOptions=regOptions(),
                     interpolationScheme=BSpline(Cubic(Line())), mask::Array{Float64,2}=ones(size(referenceImage.data)))

  parametricOnly = options.parametricOnly

  targetGrid = getCellCenteredGrid(referenceImage)

  # interpolation of the template image at transformed grid points
  transformedImage, dX_transformedImage, dY_transformedImage =
      interpolateImage(templateImage, transformedGrid, targetGrid,
        doDerivative=doDerivative, interpolationScheme=interpolationScheme)

  # measure the ssd distance
  residual = mask[:].*(transformedImage .- referenceImage.data)[:]
  prodh = prod(targetGrid.voxelsize)
  functionValue = 0.5 * prodh * sum(residual.^2)

  coarseCenteredGrid = getCellCenteredGrid(transformedGrid)

  residual  = scaledArray(residual,  targetGrid.dimensions, targetGrid.voxelsize, targetGrid.shift)
  residual = sampleArrayToGridWithSmoothing(residual , coarseCenteredGrid)[:]
  prodH = prod(coarseCenteredGrid.voxelsize)

  if doDerivative
      dX_transformedImage  = scaledArray(dX_transformedImage, targetGrid.dimensions, targetGrid.voxelsize, targetGrid.shift)
      dX_transformedImage  = sampleArrayToGridWithSmoothing(dX_transformedImage , coarseCenteredGrid)[:]
      dY_transformedImage  = scaledArray(dY_transformedImage, targetGrid.dimensions, targetGrid.voxelsize, targetGrid.shift)
      dY_transformedImage  = sampleArrayToGridWithSmoothing(dY_transformedImage , coarseCenteredGrid)[:]
  end

  if(doDerivative)
      N = prod(coarseCenteredGrid.dimensions) # product of sizes
      dFunctionValue,d2FunctionValue = DeformableRegistration.Distance.ssdDerivativesMatrixFree(dX_transformedImage, dY_transformedImage ,residual ,[1.0],prodH,N,false,doDerivative,doHessian)
      return [functionValue, dFunctionValue, d2FunctionValue, (dX_transformedImage, dY_transformedImage)]
  end

  return [functionValue, 0, 0, (0, 0)]

end


text:
"""
    wlmetric_ball(S::AbstractVector{<:GroupElem}
        [, center=one(first(S)); radius=2, op=*])
Compute metric ball as a list of elements of non-decreasing length, given the
word-length metric on the group generated by `S`. The ball is centered at `center`
(by default: the identity element). `radius` and `op` keywords specify the
radius and multiplication operation to be used.
"""
function wlmetric_ball_serial(S::AbstractVector{T}, center::T=one(first(S)); radius = 2, op = *) where {T}
    @assert radius >= 1
    old = unique!([center, [center*s for s in S]...])
    return _wlmetric_ball(S, old, radius, op, collect, unique!)
end

function wlmetric_ball_thr(S::AbstractVector{T}, center::T=one(first(S)); radius = 2, op = *) where {T}
    @assert radius >= 1
    old = unique!([center, [center*s for s in S]...])
    return _wlmetric_ball(S, old, radius, op, ThreadsX.collect, ThreadsX.unique)
end

function _wlmetric_ball(S, old, radius, op, collect, unique)
    sizes = [1, length(old)]
    for r in 2:radius
        old = let old = old, S = S,
            new = collect(
                (g = op(o, s); hash(g); g)
                for o in @view(old[sizes[end-1]:end]) for s in S
            )
            append!(old, new)
            unique(old)
        end
        push!(sizes, length(old))
    end
    return old, sizes[2:end]
end

function wlmetric_ball(
    S::AbstractVector{T},
    center::T = one(first(S));
    radius = 2,
    op = *,
    threading = true,
) where {T}
    threading && return wlmetric_ball_thr(S, center, radius = radius, op = op)
    return wlmetric_ball_serial(S, center, radius = radius, op = op)
end


text:
module Deepgroebner

include("model.jl")
# include("train.jl")
# export term, gt, LT, |, div, /, *, S, update!, select, minus, mdiv, buchberger

end


text:
function neural_ode(model,x,tspan,
                    args...;kwargs...)
  Tracker.istracked(x) && error("u0 is not currently differentiable.")
  p = destructure(model)
  dudt_(du,u::TrackedArray,p,t) = du .= restructure(model,p)(u)
  dudt_(du,u::AbstractArray,p,t) = du .= Flux.data(restructure(model,p)(u))
  prob = ODEProblem(dudt_,x,tspan,p)
  return diffeq_adjoint(p,prob,args...;kwargs...)
end

neural_msde(x,model,mp,tspan,args...;kwargs...) = neural_msde(x,model,mp,tspan,
                                                         diffeq_fd,
                                                         args...;kwargs...)
function neural_msde(model,x,mp,tspan,
                    ad_func::Function,
                    args...;kwargs...)
  p = Flux.data(destructure(model))
  dudt_(du,u::TrackedArray,p,t) = du .= restructure(model,p)(u)
  dudt_(du,u::AbstractArray,p,t) = du .= Flux.data(restructure(model,p)(u))
  g(du,u,p,t) = du .= mp.*u
  prob = SDEProblem(dudt_,g,x,tspan,p)

  if ad_func === diffeq_adjoint
    return ad_func(p,prob,args...;kwargs...)
  elseif ad_func === diffeq_fd
    return ad_func(p,Array,length(p),prob,args...;kwargs...)
  else
    return ad_func(p,Array,prob,args...;kwargs...)
  end
end


text:
"""
Polymer flooding model
injecting polymer/water mixture into a reservoir
The adsorption of polymer and change of the injected water viscosity due to the
change of concentration is included in the model.
Note: currently, the adsorption term is not included; it will be added later
"""
function polymer_implicit_upwind(core_props, fluids_pw, fluids_fw, rel_perms_fw,
  rel_perms_pw, core_flood; pv_inj = 0.3)
  c_insitu = 1.0
  c_inj = 0
  # pv_inj = 0.3

  C_high_sal=1.0
  C_low_sal=0.0
  SF_high_sal=0.0 # similar to oil wet
  SF_low_sal=1.0 # similar to water wet

  Nx = 500 # number of cells in x direction
  H = core_props.length # [m] length of the domain in y direction
  D_core = core_props.diameter
  A_core=π*D_core^2/4
  V_core=H*A_core
  m = createMesh1D(Nx, H)
  # m = createMesh2D(Nx, Ny, H, H)
  # m = createMeshCylindrical2D(Nx, Ny, W, H) # creates a 2D mesh
  ## define the physical parametrs
  # all the Corey-type relperm parametrs are defined for the oil-wet and water-wet
  # cases
  k0 = core_props.permeability # [m^2] average reservoir permeability
  phi0 = core_props.porosity # average porosity
  mu_oil_pw = fluids_pw.oil_viscosity # [Pa.s] oil viscosity
  mu_water_pw = fluids_pw.water_viscosity # [Pa.s] water viscosity
  mu_oil_fw = fluids_fw.oil_viscosity # [Pa.s] oil viscosity
  mu_water_fw = fluids_fw.water_viscosity # [Pa.s] water viscosity
  krw0_ww = rel_perms_pw.krw0
  krw0_ow = rel_perms_fw.krw0
  kro0_ww = rel_perms_pw.kro0
  kro0_ow = rel_perms_fw.kro0
  nw_ww = rel_perms_pw.nw
  nw_ow= rel_perms_fw.nw
  no_ww = rel_perms_pw.no
  no_ow= rel_perms_fw.no
  sor_ww=rel_perms_pw.sor
  sor_ow=rel_perms_fw.sor
  swc_ww=rel_perms_pw.swc
  swc_ow=rel_perms_fw.swc

  # SF0=(c_insitu-C_high_sal)/(C_low_sal-C_high_sal)
  # SF=createFaceVariable(m, SF0) # 1 is water wet, 0 is oil wet
  # krw0=krw0_ww*SF+krw0_ow*(1-SF)
  # kro0=kro0_ww*SF+kro0_ow*(1-SF)
  # sor=sor_ww*SF+sor_ow*(1-SF)
  # swc=swc_ww*SF+swc_ow*(1-SF)
  # no= no_ww*SF+no_ow*(1-SF)
  # nw= nw_ww*SF+nw_ow*(1-SF)

  # sor0=a*cos(teta0)^2+b*cos(teta0)+c

  p0 = core_flood.back_pressure # [bar] pressure
  u_in= core_flood.injection_velocity # [m/s] equal to 1 m/day
  sw0 = core_flood.initial_water_saturation
  # sw0(10:end-10, 10:end-10)=swc+0.2
  # sw0 = swc+0.1 # initial water saturation
  sw_in = core_flood.injected_water_saturation

  eps1=1e-5
  perm_val= k0 #field2d(Nx,Ny,k0,V_dp,clx,cly)

  k=createCellVariable(m, perm_val)
  phi=createCellVariable(m, phi0)

  mu_water_pw_cell= createCellVariable(m, mu_water_pw)
  mu_water_fw_cell= createCellVariable(m, mu_water_fw)
  mu_water_pw_face= arithmeticMean(mu_water_pw_cell)
  mu_water_fw_face= arithmeticMean(mu_water_fw_cell)
  k_face = geometricMean(k)
  lw = k_face/mu_water_fw_face
  lo = k_face/mu_oil_fw

  ## Define the boundaries: all fixed Sw=1, fixed pressure everywhere(?)
  BCp = createBC(m) # Neumann BC for pressure
  BCs = createBC(m) # Neumann BC for saturation
  BCc = createBC(m) # Neumann BC for salinity
  # left boundary pressure gradient
  # BCp.left.a[:]=(krw(sw_in)*lw.xvalue(1,:)+kro(sw_in)*lo.xvalue(1,:)) BCp.left.b[:]=0 BCp.left.c[:]=-u_in
  # change the right boandary to constant pressure (Dirichlet)
  BCp.left.a[:].=-k0/mu_water_pw
  BCp.left.b[:].=0
  BCp.left.c[:].=u_in
  BCp.right.a[:].=0.0
  BCp.right.b[:].=1.0
  BCp.right.c[:].=p0
  # BCp.top.a[:]=0.0
  # BCp.top.b[:]=1.0
  # BCp.top.c[:]=p0
  # BCp.bottom.a[:]=0.0
  # BCp.bottom.b[:]=1.0
  # BCp.bottom.c[:]=p0
  # change the left boundary to constant saturation (Dirichlet)
  # BCs.left.a[:]=0 BCs.left.b[:]=1 BCs.left.c[:]=1.0-sor
  # calculate the water saturation at the boundaries
  # sw_b=fzero(f1, 0.5)
  # println(sw_b)
  # for test, no capillary end effect
  # BCs.right.a[:]=0.0
  # BCs.right.b[:]=1.0
  # BCs.right.c[:]=sw_b
  BCs.left.a[:].=0.0
  BCs.left.b[:].=1.0
  BCs.left.c[:].=sw_in
  # BCs.top.a[:]=0.0
  # BCs.top.b[:]=1.0
  # BCs.top.c[:]=sw_b
  # BCs.bottom.a[:]=0.0
  # BCs.bottom.b[:]=1.0
  # BCs.bottom.c[:]=sw_b
  # Boundary conditions for salinity
  BCc.left.a[:].=0.0
  BCc.left.b[:].=1.0
  BCc.left.c[:].=c_inj

  ## define the time step and solver properties
  # dt = 1000 # [s] time step
  # dt=(W/Nx)/u_in/20 # [s]
  # pv_inj=5
  t_end = pv_inj*V_core*phi0/(u_in*A_core) # [s] final time
  dt=H/Nx/u_in/20
  SW0 = sw0
  eps_p = 1e-3 # pressure accuracy
  eps_sw = 1e-6 # saturation accuracy
  ## define the variables
  sw_old = createCellVariable(m, SW0, BCs)
  sw_new=createCellVariable(m, SW0, BCs)
  p_old = createCellVariable(m, p0, BCp)
  p_new= createCellVariable(m, p0, BCp)
  c_old = createCellVariable(m, c_insitu, BCc)
  c_face=arithmeticMean(c_old) # just to initialize
  SF=(1/(C_low_sal-C_high_sal))*(c_face-C_high_sal)
  # calculate the new relperm curves
  krw0=krw0_ww*SF+krw0_ow*(1-SF)
  kro0=kro0_ww*SF+kro0_ow*(1-SF)
  sor=sor_ww*SF+sor_ow*(1-SF)
  swc=swc_ww*SF+swc_ow*(1-SF)
  no= no_ww*SF+no_ow*(1-SF)
  nw= nw_ww*SF+nw_ow*(1-SF)

  # rel_perm_av = oil_water_rel_perms(krw0 = krw0, kro0=kro0, swc=swc, sor=sor, nw=nw, no=no)

  sw = copyCell(sw_old)
  oil_init=domainInt(1-sw_old)
  p = copyCell(p_old)
  uw = -gradientTerm(p_old) # an estimation of the water velocity
  ## start the main loop
  # generate intial pressure profile (necessary to initialize the fully
  # implicit solver)
  rec_fact=zeros(1)
  c_out_sal=zeros(1)
  water_cut=zeros(1)
  c_out_sal[1]=c_face.xvalue[end]
  t_sec=zeros(1)
  t = 0.0
  dt0=dt
  dsw_alwd= 0.005
  dp_alwd= 100.0 # Pa
  dc_alwd= 0.005
  labdaw=createFaceVariable(m,1.0) # initialize
  labdao=createFaceVariable(m,1.0) # initialize
  FL=fluxLimiter("SUPERBEE")

  prog_1=Progress(100, 1)

  while (t<t_end)
  # print("progress is $((t/t_end*100)) %", "\u1b[1G")
  ProgressMeter.update!(prog_1, floor(Int, t/t_end*100))
  # Implicit loop
  error_p = 1e5
  error_sw = 1e5
  # implicit solver
  loop_count=0
  while ((error_sw>eps_sw))
    loop_count+=1
    if loop_count>10
      break
    end
    # calculate parameters
    pgrad = gradientTerm(p)
    uw = -labdaw*pgrad
  #         pcgrad=gradientTerm(pc(sw))
    sw_face = upwindMean(sw, uw) # average value of water saturation
    # sw_ave=arithmeticMean(sw)
    # solve for pressure at known Sw
    labdao = lo*faceEval(CF.kro, sw_face, kro0, sor, swc, no)
    labdaw = lw*faceEval(CF.krw, sw_face, krw0, sor, swc, nw)
    dlabdaodsw = lo*faceEval(CF.dkrodsw, sw_face, kro0, sor, swc, no)
    dlabdawdsw = lw*faceEval(CF.dkrwdsw, sw_face, krw0, sor, swc, nw)
    labda = labdao+labdaw
    dlabdadsw = dlabdaodsw+dlabdawdsw
    # compute [Jacobian] matrices
    Mdiffp1 = diffusionTerm(-labda)
    Mdiffp2 = diffusionTerm(-labdaw)
    # Mconvsw1 = convectionUpwindTerm(-dlabdadsw.*pgrad, uw)
    # Mconvsw2 = convectionUpwindTerm(-dlabdawdsw.*pgrad, uw)
    Mconvsw1 = convectionUpwindTerm(-dlabdadsw*pgrad, uw)
    Mconvsw2 = convectionUpwindTerm(-dlabdawdsw*pgrad, uw)
    Mtranssw2, RHStrans2 = transientTerm(sw_old, dt, phi)
    # Compute RHS values
    RHS1 = Mconvsw1*sw.value # divergenceTerm(-dlabdadsw.*sw_face.*pgrad)
    RHS2 = Mconvsw2*sw.value # divergenceTerm(-dlabdawdsw.*sw_face.*pgrad)
    # include boundary conditions
    Mbcp, RHSbcp = boundaryConditionTerm(BCp)
    Mbcsw, RHSbcsw = boundaryConditionTerm(BCs)
    # Couple the equations; BC goes into the block on the main diagonal
    M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw]
    RHS = [RHS1+RHSbcp; RHS2+RHStrans2+RHSbcsw]
    x = M\RHS
    p_new.value[:] .= (x[1:(Nx+2)])
    sw_new.value[:] .= (x[(Nx+2)+1:end])

    error_p = maximum(abs.((internalCells(p_new).-internalCells(p))./internalCells(p_new)))
    error_sw = maximum(abs.(internalCells(sw_new).-internalCells(sw)))
    # dt_new=dt*min(dp_alwd/error_p, dsw_alwd/error_sw)
    # print("sw_error = $error_sw \n")
    # assign new values of p and sw
    p.value[:]=p_new.value[:]
    sw.value[:]=sw_new.value[:]
  end # end of implicit loop
  if loop_count>10
    p=copyCell(p_old)
    sw=copyCell(sw_old)
    dt=dt/5
    continue
  end
  # solve the advection equation
  # ϕ d/dt(c sw)=[sw dc/dt + c dsw/dt]
  dswdt=(sw_new-sw_old)/dt
  Msc=linearSourceTerm(phi*dswdt)
  sw_face = upwindMean(sw, -labdaw*pgrad) # average value of water saturation
  pgrad = gradientTerm(p_new)
  labdaw = lw*faceEval(CF.krw, sw_face, krw0, sor, swc, nw)
  uw=-labdaw*pgrad
  Mtc, RHStc=transientTerm(c_old, dt, phi*sw)
  Mbcc, RHSbcc=boundaryConditionTerm(BCc)
  c_temp=copyCell(c_old)
  for j in 1:4
    Mconvc, RHSconv=convectionTvdTerm(uw, c_temp, FL)
    c_temp=solvePDE(m, Msc+Mtc+Mconvc+Mbcc, RHStc+RHSbcc+RHSconv)
  end
  c_new=copyCell(c_temp)
  c_face=tvdMean(c_new, uw, FL)

  # calculate the wettability modifier
  SF=(1/(C_low_sal-C_high_sal))*(c_face-C_high_sal)
  # SF.xvalue[SF.xvalue.>1.0]=1.0
  # SF.xvalue[SF.xvalue.<0.0]=0.0

  # calculate the new relperm curves
  krw0=krw0_ww*SF+krw0_ow*(1-SF)
  kro0=kro0_ww*SF+kro0_ow*(1-SF)
  sor=sor_ww*SF+sor_ow*(1-SF)
  swc=swc_ww*SF+swc_ow*(1-SF)
  no= no_ww*SF+no_ow*(1-SF)
  nw= nw_ww*SF+nw_ow*(1-SF)

  # Calculate the viscosity
  mu_w_face = SF*mu_water_pw_face+(1-SF)*mu_water_fw_face
  lw = k_face/mu_w_face

  # adaptive time step
  dp = maximum(abs.((internalCells(p_new)-internalCells(p_old))./internalCells(p_new)))
  dsw = maximum(abs.(internalCells(sw_new)-internalCells(sw_old)))
  dc = maximum(abs.(internalCells(c_new)-internalCells(c_old)))
  # dt=min(dt*(dsw_alwd/dsw), dt*(dc_alwd/dc), 15*dt, t_end-t-eps())

  # update the variables
  t+=dt
  c_old=copyCell(c_new)
  p_old = copyCell(p)
  sw_old = copyCell(sw)

  # calculate the recovery factor
  push!(rec_fact, (oil_init-domainInt(1-sw))/oil_init)
  push!(t_sec, t)
  push!(c_out_sal, c_face.xvalue[end])
  push!(water_cut, labdaw.xvalue[end]/(labdaw.xvalue[end]+labdao.xvalue[end]))

  # print(t)
  #GR.plot(sw.value[2:end-1])
  # GR.plot(c_new.value[2:end-1])
  #GR.plot(SF.xvalue)
  # GR.imshow(sw.value[2:end-1,2:end-1])
  # visualizeCells(1-sw)
  # GR.plot(t_sec/3600/24, rec_fact)
  # xlabel('time [day]')
  # ylabel('recovery factor')
  # title([num2str(t/3600/24) ' day']) drawnow
  end
  pv = t_sec*u_in/(H*phi0)
  sw_face = linearMean(sw)
  c_face  = linearMean(c_old)
  x = m.facecenters.x
  t_sec, pv, rec_fact, x, sw_face.xvalue, c_face.xvalue, c_out_sal, water_cut
end # imb_impes

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
### For simplicity we assume the same group to be ussed everywhere.
G = CryptoGroups.MODP160Group()
signdata(data,signer) = DSASignature(hash("$data"),signer)
verifydata(data,signature) = verify(signature,G) && signature.hash==hash("$data")

serverkey = Signer(G)
serverid = hash(serverkey.pubkey)

maintainerkey = Signer(G)
maintainerid = hash(maintainerkey.pubkey)

user1key = Signer(G)
user2key = Signer(G)
user3key = Signer(G)

config = ServerConfig(1999,maintainerid,G)
server = Server(config,data->signdata(data,serverkey),verifydata)

sleep(1.0)

maintainer = Maintainer(1999,serverid,data->signdata(data,maintainerkey),verifydata)

ballotboxcfg = BallotBoxConfig(G)
start(maintainer,ballotboxcfg,2001)
push!(maintainer,serverid,2001)

broute = BallotBoxRoute(2001,serverid)
gkconfig = GateKeeperConfig(3,broute,G)
start(maintainer,gkconfig,2000)

push!(maintainer,hash(user1key.pubkey),2000)
push!(maintainer,hash(user2key.pubkey),2000)
push!(maintainer,hash(user3key.pubkey),2000)

### Users do:
sleep(1.0)

gk = GateKeeperRoute(2000,serverid,serverid)

@async vote("msg1",gk,data->signdata(data,user1key),verifydata)
@async vote("msg2",gk,data->signdata(data,user2key),verifydata)
@async vote("msg3",gk,data->signdata(data,user3key),verifydata)

### After that maintainer gets a ballot

sleep(5.0)

@show takeballot!(maintainer,2000)

### Now let's test closing

stop(maintainer,2001)
stop(maintainer,2000)

stop(server)


text:
function approximation_term(
    ::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::IteratedLagrangianRefinementAlgorithm,
) where {T, U}
    return NaN
end
function approximation_ratio(
    ::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::IteratedLagrangianRefinementAlgorithm,
) where {T, U}
    return 0.5
end

function solve(
    i::MinimumBudget{SpanningTreeInstance{T, Maximise}, U},
    ::IteratedLagrangianRefinementAlgorithm,
) where {T, U}
    # Approximately solve the following problem:
    #     \max_{x spanning tree} rewards x  s.t.  weights x >= budget
    # This algorithm provides a multiplicative approximation to this problem. If x* is the optimum solution and x~ the one
    # returned by this algorithm,
    #     weights x* >= budget   and   weights x~ >= budget                 (the returned solution is feasible)
    #     rewards x~ >= rewards x* / 2                                      (multiplicative approximation)

    # Based on the same ideas as https://doi.org/10.1007/s10107-009-0307-4

    # If there are too few edges, not much to do.
    if ne(i.instance.graph) <= 2
        return solve(i, LagrangianRefinementAlgorithm())
    end

    # For each pair of edges, force these two edges to be part of the solution and discard all edges with a higher value.
    best_sol = nothing
    best_sol_value = -Inf
    for e1 in edges(i.instance.graph)
        for e2 in edges(i.instance.graph)
            # (e1, e2) must be a pair of distinct edges.
            if (src(e1) == src(e2) && dst(e1) == dst(e2)) ||
               (src(e1) == dst(e2) && src(e1) == dst(e2))
                continue
            end

            # Filter out the edges that have a higher value than any of these two edges. Give a very large reward to them both.
            cutoff = min(i.instance.rewards[e1], i.instance.rewards[e2])
            rewards = filter(kv -> kv[2] < cutoff, i.instance.rewards)
            rewards[e1] = rewards[e2] = prevfloat(Inf)

            graph = SimpleGraph(nv(i.instance.graph))
            for e in keys(rewards)
                add_edge!(graph, e)
            end

            weights = Dict(e => i.weights[e] for e in keys(rewards))

            # No other simplification can be made, unfortunately: all other edges might participate in the optimum solution.
            # (Only exception: e1 and e2 are incident to the same vertex; any edge linking the other extremities of these edges
            # can be removed, but that's only one edge at most.)
            # As e1 and e2 have the best reward (as they are really bumped), they must be in any optimum solution.

            # Solve this subproblem.
            bsti = MinimumBudget(
                SpanningTreeInstance(graph, rewards),
                weights,
                i.min_budget,
            )
            sol = solve(bsti, LagrangianRefinementAlgorithm())

            # This subproblem is infeasible. Maybe it's because the overall problem is infeasible or just because too many
            # edges were removed.
            if length(sol.variables) == 0
                continue
            end

            # Impossible to have a feasible solution with these two edges, probably because of the budget constraint.
            # Due to the direction of the budget constraint (>= budget), it is not possible to check for feasibility
            # before solving an instance.
            if !(e1 in sol.variables) || !(e2 in sol.variables)
                continue
            end

            # Only keep the best solution.
            sol_value = _budgeted_spanning_tree_compute_value(i, sol.variables)
            if sol_value > best_sol_value
                # sol's instance is the one used internally for the subproblems.
                best_sol = SimpleBudgetedSpanningTreeSolution(i, sol.variables)
                best_sol_value = sol_value
            end
        end
    end

    if best_sol !== nothing
        return best_sol
    else
        return SimpleBudgetedSpanningTreeSolution(i)
    end
end


text:
function _check_constraint_bounds(bounds::PSI.ConstraintBounds, valid_bounds::NamedTuple)
    @test bounds.coefficient.min == valid_bounds.coefficient.min
    @test bounds.coefficient.max == valid_bounds.coefficient.max
    @test bounds.rhs.min == valid_bounds.rhs.min
    @test bounds.rhs.max == valid_bounds.rhs.max
end

function _check_variable_bounds(bounds::PSI.VariableBounds, valid_bounds::NamedTuple)
    @test bounds.bounds.min == valid_bounds.min
    @test bounds.bounds.max == valid_bounds.max
end

@testset "Test Numerical Stability of Constraints" begin
    template = get_thermal_dispatch_template_network()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    valid_bounds =
        (coefficient = (min = 1.0, max = 1.0), rhs = (min = 0.4, max = 9.930296584))
    model = OperationsProblem(template, c_sys5; optimizer = GLPK_optimizer)
    @test build!(model; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT

    bounds = PSI.get_constraint_numerical_bounds(model; verbose = false)
    _check_constraint_bounds(bounds, valid_bounds)

    model_bounds = PSI.get_constraint_numerical_bounds(model; verbose = true)
    valid_model_bounds = Dict(
        :CopperPlateBalance => (
            coefficient = (min = 1.0, max = 1.0),
            rhs = (min = 6.434489705000001, max = 9.930296584),
        ),
        :P_lb__ThermalStandard__RangeConstraint =>
            (coefficient = (min = 1.0, max = 1.0), rhs = (min = Inf, max = -Inf)),
        :P_ub__ThermalStandard__RangeConstraint =>
            (coefficient = (min = 1.0, max = 1.0), rhs = (min = 0.4, max = 6.0)),
    )
    for (constriant_key, constriant_bounds) in model_bounds
        _check_constraint_bounds(constriant_bounds, valid_model_bounds[constriant_key])
    end
end

@testset "Test Numerical Stability of Variables" begin
    template = get_template_basic_uc_simulation()
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5_uc")
    valid_bounds = (min = 0.0, max = 6.0)
    model = OperationsProblem(template, c_sys5; optimizer = GLPK_optimizer)
    @test build!(model; output_dir = mktempdir(cleanup = true)) == PSI.BuildStatus.BUILT

    bounds = PSI.get_variable_numerical_bounds(model; verbose = false)
    _check_variable_bounds(bounds, valid_bounds)

    model_bounds = PSI.get_variable_numerical_bounds(model; verbose = true)
    valid_model_bounds = Dict(
        :stop__ThermalStandard => (min = 0.0, max = 1.0),
        :start__ThermalStandard => (min = 0.0, max = 1.0),
        :P__ThermalStandard => (min = 0.4, max = 6.0),
        :On__ThermalStandard => (min = 0.0, max = 1.0),
    )
    for (variable_key, variable_bounds) in model_bounds
        _check_variable_bounds(variable_bounds, valid_model_bounds[variable_key])
    end
end


text:
using DashBootstrapComponents

modal = html_div([
    dbc_button("Open", id = "open-centered"),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Header"), close_button = true),
            dbc_modalbody("This modal is vertically centered"),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id = "close-centered",
                    className = "ms-auto",
                    n_clicks = 0,
                ),
            ),
        ],
        id = "modal-centered",
        centered = true,
        is_open = false,
    ),
]);

callback!(
    app,
    Output("modal-centered", "is_open"),
    Input("open-centered", "n_clicks"),
    Input("close-centered", "n_clicks"),
    State("modal-centered", "is_open"),
) do n1, n2, is_open
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# XML Document
# ============

"""
An XML/HTML document type.
"""
struct Document
    node::Node

    function Document(ptr::Ptr{_Node})
        @assert ptr != C_NULL
        ntype = unsafe_load(ptr).typ
        @assert ntype ∈ (DOCUMENT_NODE, HTML_DOCUMENT_NODE)
        return new(Node(ptr))
    end
end

# Document properties
# -------------------

@inline function Base.getproperty(doc::Document, name::Symbol)
    name == :version  ? (hasversion(doc)  ? version(doc)  : nothing) :
    name == :encoding ? (hasencoding(doc) ? encoding(doc) : nothing) :
    name == :root     ? (hasroot(doc)     ? root(doc)     : nothing) :
    name == :dtd      ? (hasdtd(doc)      ? dtd(doc)      : nothing) :
    Core.getfield(doc, name)
end

# Document constructors
# ---------------------

"""
    XMLDocument(version="1.0")

Create an XML document.
"""
function XMLDocument(version::AbstractString="1.0")
    node = XMLDocumentNode(version)
    return Document(node.ptr)
end

"""
    HTMLDocument(uri=nothing, externalID=nothing)

Create an HTML document.
"""
function HTMLDocument(uri=nothing, externalID=nothing)
    node = HTMLDocumentNode(uri, externalID)
    return Document(node.ptr)
end

function Base.print(io::IO, doc::Document)
    print(io, doc.node)
end

"""
    prettyprint([io], doc::Document)

Print `doc` with formatting.
"""
function prettyprint(doc::Document)
    prettyprint(stdout, doc)
end

function prettyprint(io::IO, doc::Document)
    prettyprint(io, doc.node)
end

"""
    parsexml(xmlstring)

Parse `xmlstring` and create an XML document.
"""
function parsexml(xmlstring::AbstractString)
    if isempty(xmlstring)
        throw(ArgumentError("empty XML string"))
    end
    doc_ptr = @check ccall(
        (:xmlParseMemory, libxml2),
        Ptr{_Node},
        (Cstring, Cint),
        xmlstring, sizeof(xmlstring)) != C_NULL
    return Document(doc_ptr)
end

function parsexml(xmldata::Vector{UInt8})
    return parsexml(String(xmldata))
end

function parsexml(xmldata::Base.CodeUnits{UInt8,String})
    return parsexml(String(xmldata))
end

"""
    parsehtml(htmlstring)

Parse `htmlstring` and create an HTML document.
"""
function parsehtml(htmlstring::AbstractString)
    if isempty(htmlstring)
        throw(ArgumentError("empty HTML string"))
    end
    url = C_NULL
    encoding = C_NULL
    options = 1
    doc_ptr = @check ccall(
        (:htmlReadMemory, libxml2),
        Ptr{_Node},
        (Cstring, Cint, Cstring, Cstring, Cint),
        htmlstring, sizeof(htmlstring), url, encoding, options) != C_NULL
    show_warnings()
    return Document(doc_ptr)
end

function parsehtml(htmldata::Vector{UInt8})
    return parsehtml(String(htmldata))
end

function parsehtml(htmldata::Base.CodeUnits{UInt8,String})
    return parsehtml(String(htmldata))
end

"""
    readxml(filename)

Read `filename` and create an XML document.
"""
function readxml(filename::AbstractString)
    encoding = C_NULL
    # Do not reuse the context dictionary.
    options = 1 << 12
    doc_ptr = @check ccall(
        (:xmlReadFile, libxml2),
        Ptr{_Node},
        (Cstring, Ptr{UInt8}, Cint),
        filename, encoding, options) != C_NULL
    return Document(doc_ptr)
end

"""
    readxml(input::IO)

Read `input` and create an XML document.
"""
function readxml(input::IO)
    readcb = make_read_callback()
    closecb = C_NULL
    context = pointer_from_objref(input)
    uri = C_NULL
    encoding = C_NULL
    options = 0
    doc_ptr = @check ccall(
        (:xmlReadIO, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring, Cint),
        readcb, closecb, context, uri, encoding, options) != C_NULL
    return Document(doc_ptr)
end

"""
    readhtml(filename)

Read `filename` and create an HTML document.
"""
function readhtml(filename::AbstractString)
    encoding = C_NULL
    options = 0
    doc_ptr = @check ccall(
        (:htmlReadFile, libxml2),
        Ptr{_Node},
        (Cstring, Cstring, Cint),
        filename, encoding, options) != C_NULL
    show_warnings()
    return Document(doc_ptr)
end

"""
    readhtml(input::IO)

Read `input` and create an HTML document.
"""
function readhtml(input::IO)
    readcb = make_read_callback()
    closecb = C_NULL
    context = pointer_from_objref(input)
    uri = C_NULL
    encoding = C_NULL
    options = 0
    doc_ptr = @check ccall(
        (:htmlReadIO, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring, Cint),
        readcb, closecb, context, uri, encoding, options) != C_NULL
    show_warnings()
    return Document(doc_ptr)
end

function Base.write(filename::AbstractString, doc::Document)
    format = 0
    encoding = "UTF-8"
    ret = @check ccall(
        (:xmlSaveFormatFileEnc, libxml2),
        Cint,
        (Cstring, Ptr{Cvoid}, Cstring, Cint),
        filename, doc.node.ptr, encoding, format) != -1
    return Int(ret)
end

function make_read_callback()
    # Passing an input stream as an argument is impossible to create a callback
    # because Julia does not support C-callable closures yet.
    return @cfunction(Cint, (Ref{IO}, Ptr{UInt8}, Cint)) do context, buffer, len
        input = context
        avail = min(bytesavailable(input), len)
        if avail > 0
            unsafe_read(input, buffer, avail)
            read = avail
        elseif len > 0 && !eof(input)
            # An input stream may return bytesavailable = 0 before reading data.
            # So, read a byte to kick it ready.
            unsafe_store!(buffer, Base.read(input, UInt8))
            read = 1
        else
            read = 0
        end
        # debug
        # @show unsafe_string(buffer, read)
        return Cint(read)
    end
end

# Properties
# ----------

"""
    hasversion(doc::Document)

Return if `doc` has a version string.
"""
function hasversion(doc::Document)
    doc_ptr = Ptr{_Document}(doc.node.ptr)
    @assert doc_ptr != C_NULL
    return unsafe_load(doc_ptr).version != C_NULL
end

"""
    version(doc::Document)

Return the version string of `doc`.
"""
function version(doc::Document)
    if !hasversion(doc)
        throw(ArgumentError("no version string"))
    end
    doc_ptr = Ptr{_Document}(doc.node.ptr)
    @assert doc_ptr != C_NULL
    return unsafe_string(unsafe_load(doc_ptr).version)
end

"""
    hasencoding(doc::Document)

Return if `doc` has an encoding string.
"""
function hasencoding(doc::Document)
    doc_ptr = Ptr{_Document}(doc.node.ptr)
    @assert doc_ptr != C_NULL
    return unsafe_load(doc_ptr).encoding != C_NULL
end

"""
    encoding(doc::Document)

Return the encoding string of `doc`.
"""
function encoding(doc::Document)
    if !hasencoding(doc)
        throw(ArgumentError("no encoding string"))
    end
    doc_ptr = Ptr{_Document}(doc.node.ptr)
    @assert doc_ptr != C_NULL
    return unsafe_string(unsafe_load(doc_ptr).encoding)
end

# Root
# ----

"""
    hasroot(doc::Document)

Return if `doc` has a root element.
"""
function hasroot(doc::Document)
    ptr = ccall(
        (:xmlDocGetRootElement, libxml2),
        Ptr{Cvoid},
        (Ptr{Cvoid},),
        doc.node.ptr)
    return ptr != C_NULL
end

"""
    root(doc::Document)

Return the root element of `doc`.
"""
function root(doc::Document)
    if !hasroot(doc)
        throw(ArgumentError("no root element"))
    end
    root_ptr = @check ccall(
        (:xmlDocGetRootElement, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid},),
        doc.node.ptr) != C_NULL
    return Node(root_ptr)
end

"""
    setroot!(doc::Document, node::Node)

Set the root element of `doc` to `node` and return the root element.
"""
function setroot!(doc::Document, root::Node)
    if !iselement(root)
        throw(ArgumentError("not an element node"))
    end
    old_root_ptr = ccall(
        (:xmlDocSetRootElement, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid}, Ptr{Cvoid}),
        doc.node.ptr, root.ptr)
    update_owners!(root, doc.node)
    if old_root_ptr != C_NULL
        old_root = Node(old_root_ptr)
        update_owners!(old_root, old_root)
    end
    return root
end

# DTD and validation
# ------------------

"""
    hasdtd(doc::Document)

Return if `doc` has a DTD node.
"""
function hasdtd(doc::Document)
    dtd_ptr = ccall(
        (:xmlGetIntSubset, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid},),
        doc.node.ptr)
    return dtd_ptr != C_NULL
end

"""
    dtd(doc::Document)

Return the DTD node of `doc`.
"""
function dtd(doc::Document)
    if !hasdtd(doc)
        throw(ArgumentError("no DTD"))
    end
    dtd_ptr = @check ccall(
        (:xmlGetIntSubset, libxml2),
        Ptr{_Node},
        (Ptr{Cvoid},),
        doc.node.ptr) != C_NULL
    return Node(dtd_ptr)
end

"""
    setdtd!(doc::Document, node::Node)

Set the DTD node of `doc` to `node` and return the DTD node.
"""
function setdtd!(doc::Document, node::Node)
    if !isdtd(node)
        throw(ArgumentError("not a DTD node"))
    elseif hasdtd(doc)
        unlink!(dtd(doc))
    end
    # Insert `node` as the first child of `doc.node`.
    if hasnode(doc.node)
        linkprev!(firstnode(doc.node), node)
    else
        link!(doc.node, node)
    end
    return node
end

struct _ValidCtxt
    # type tag
end

"""
    readdtd(filename::AbstractString)

Read `filename` and create a DTD node.
"""
function readdtd(filename::AbstractString)
    dtd_ptr = @check ccall(
        (:xmlParseDTD, libxml2),
        Ptr{_Node},
        (Cstring, Cstring),
        C_NULL, filename) != C_NULL
    return Node(dtd_ptr)
end

"""
    validate(doc::Document, [dtd::Node])

Validate `doc` against `dtd` and return the validation log.

The validation log is empty if and only if `doc` is valid. The DTD node in `doc`
will be used if `dtd` is not passed.
"""
function validate(doc::Document)
    ctxt_ptr = new_valid_context()
    @assert isempty(XML_GLOBAL_ERROR_STACK)
    valid = ccall(
        (:xmlValidateDocument, libxml2),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}),
        ctxt_ptr, doc.node.ptr)
    free(ctxt_ptr)
    @assert (valid == 1) == isempty(XML_GLOBAL_ERROR_STACK)
    log = copy(XML_GLOBAL_ERROR_STACK)
    empty!(XML_GLOBAL_ERROR_STACK)
    return log
end

function validate(doc::Document, dtd::Node)
    ctxt_ptr = new_valid_context()
    @assert isempty(XML_GLOBAL_ERROR_STACK)
    valid = ccall(
        (:xmlValidateDtd, libxml2),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        ctxt_ptr, doc.node.ptr, dtd.ptr)
    free(ctxt_ptr)
    @assert (valid == 1) == isempty(XML_GLOBAL_ERROR_STACK)
    log = copy(XML_GLOBAL_ERROR_STACK)
    empty!(XML_GLOBAL_ERROR_STACK)
    return log
end

function new_valid_context()
    ctxt_ptr = @check ccall(
        (:xmlNewValidCtxt, libxml2),
        Ptr{_ValidCtxt},
        ()) != C_NULL
    return ctxt_ptr
end

function free(ptr::Ptr{_ValidCtxt})
    ccall(
        (:xmlFreeValidCtxt, libxml2),
        Cvoid,
        (Ptr{Cvoid},),
        ptr)
end


text:
using Test
using SimpleTensorNetworks
using CUDA
using SimpleTensorNetworks: c2l, l2c
using Random

CUDA.allowscalar(false)

@testset "c2l" begin
    for i=1:100
        shape = (4,rand(1:5),rand(1:7),5,19)
        target = ([rand(1:s) for s in shape]...,)
        @test c2l(shape, target) == LinearIndices(shape)[target...]
    end
    for i=1:100
        shape = (4,rand(1:5),rand(1:12),15,19)
        ci = CartesianIndices(shape)
        i = rand(1:prod(shape))
        @test l2c(shape, i) == ci[i].I
    end
    @test l2c((), 1) == ()
    @test c2l((), ()) == 1
end

@testset "permutedims" begin
    a = randn(rand(1:3, 20)...)
    A = CuArray(a)
    p = randperm(20)
    @test Array(permutedims(A, p)) ≈ permutedims(a, p)
end

@testset "tensor contract - GPU" begin
    A = zeros(Float64, 10, 32, 21);
    B = zeros(Float64, 32, 11, 5, 2, 41, 10);
    tA = LabeledTensor(A, [1,2,3])
    tB = LabeledTensor(B, [2,4,5,6,7,1])
    tOut = tA * tB
    tnet = TensorNetwork([tA, tB]) |> togpu

    tOut2, contracted_labels = contract_label!(tnet, 1)
    @test Array(tnet.tensors[].array) ≈ tOut.array
    @test Array(tnet.tensors[].array) ≈ Array(tOut2.array)
    @test contracted_labels == [1, 2]
end


text:
std_rnorm = Normal()

abstract type ConditionalDistribution end

struct GaussianConditional <: ConditionalDistribution
    mu_intercept::Real
    mu_slope::Real
    sigma::Real
end

struct DiscreteConditional <: ConditionalDistribution
    p_0::Real
end

function sample(conditional::GaussianConditional, condition)
    x = rand(std_rnorm, 1)
    return x[1] * conditional.sigma + conditional.mu_intercept +
        conditional.mu_slope * condition
    return x
end

function sample(conditional::DiscreteConditional, condition)

end

function gibbs_sampler(px_y, py_x; y0=0, num_samples=10000, burn_in=100, thin=10)
    x_vec = zeros(0)
    y_vec = zeros(0)
    y = y0
    for i in 1:num_samples
        x = sample(px_y, y)
        y = sample(py_x, x)
        if((i > burn_in) && (i % thin == 0))
            push!(x_vec, x)
            push!(y_vec, y)
        end
    end
    return x_vec, y_vec
end

function check_discrete_compatibility(px_y, py_x)

end

function check_gaussian_compatibility(px_y, py_x)

end

function inconsistent_gaussian_gibbs()
    px_y = GaussianConditional(0, 1.5, 2)
    py_x = GaussianConditional(0, 0.5, 1)
    x,y = gibbs_sampler(px_y, py_x)
    p = plot(x, y, seriestype=:scatter, xlabel="X", ylabel="Y", legend=false)
    display(p)
    return x,y
end


text:
# bc_solvers.jl for advection

"""
flux1

Calculates the boundary flux for the advection equation. It works at the nodal
level.

**Inputs**

*  `u_sbp_`: The entry from u_sbp for this node
*  `dxidx` : The jacobian for this node
*  `nrm`   : nrm is the normal vector
*  `net_flux`:
*  `params`: the equation ParamType

**Outputs**

*  None
"""
function flux1(params::ParamType2, u_sbp_, dxidx, nrm, net_flux )
  # This function works at the nodal level
  # u_sbp_ is the entry from u_sbp for this node
  # dxi_dx is the jacobian for this node
  # nrm is the normal vector

  u_xi = dxidx[1,1]*alpha_x + dxidx[1,2]*alpha_y
  u_eta = dxidx[2,1]*alpha_x + dxidx[2,2]*alpha_y
  mag_x = u_xi*nrm[1]  # x flow magnitude
  mag_y = u_eta*nrm[2]  # y flow magnitude
  mag = mag_x + mag_y

  u_bc = sin(-1)  # boundary condition (for all sides)
  # println("u_bc = ", u_bc)

  # because mag is scalar,not vector, can combine these if statements
  if mag < 0  # inflow condition, added factor of nrm[2]
    flx_xi = u_xi*u_bc
  else
    flx_xi = u_xi*u_sbp_
  end

  if mag < 0  # inflow condition, added factor of nrm[1]
    flx_eta = u_eta*u_bc
  else
    flx_eta = u_eta*u_sbp_
  end

  net_flux = flx_xi*nrm[1] + flx_eta*nrm[2]

  return nothing
end # end function flux1

#TODO: consolidate this with LFFlux?
"""
### AdvectionEquationMod.RoeSolver

Roe solver for the advection equations. It determines the boundary flux on
each boundary. It is called at the nodal level

**Inputs**

*  `u`    : Solution of advection equation at a particular node
*  `u_bc` : Prescribed solution value at the boundary
*  `params`: the equation ParamType object
*  `nrm`  : scaled face normal vector in x-y space

**Outputs**

*  `bndryflux` : Boundary flux at the particular node

"""
function RoeSolver(params::ParamType2, u::Tsol, u_bc, nrm) where Tsol
    alpha_x = params.alpha_x
    alpha_y = params.alpha_y
    #=
    alpha_xi = dxidx[1,1]*alpha_x + dxidx[1,2]*alpha_y
    alpha_eta = dxidx[2,1]*alpha_x + dxidx[2,2]*alpha_y
    alpha_n  = alpha_xi*nrm[1] + alpha_eta*nrm[2]
    =#
    alpha_n = alpha_x*nrm[1] + alpha_y*nrm[2]
    bndryflux = 0.5*alpha_n.*(u_bc + u) - 0.5*absvalue(alpha_n).*(u_bc - u)

  return bndryflux
end # end function RoeSolver

function RoeSolver(params::ParamType3, u::Tsol, u_bc, nrm) where Tsol
    alpha_x = params.alpha_x
    alpha_y = params.alpha_y
    alpha_z = params.alpha_z
    #=
    alpha_xi = dxidx[1,1]*alpha_x + dxidx[1,2]*alpha_y + dxidx[1,3]*alpha_z
    alpha_eta = dxidx[2,1]*alpha_x + dxidx[2,2]*alpha_y + dxidx[2,3]*alpha_z
    alpha_psi = dxidx[3,1]*alpha_x + dxidx[3,2]*alpha_y + dxidx[3,3]*alpha_z

    alpha_n  = alpha_xi*nrm[1] + alpha_eta*nrm[2] + alpha_psi*nrm[3]
    =#
    alpha_n = alpha_x*nrm[1] + alpha_y*nrm[2] + alpha_z*nrm[3]
    bndryflux = 0.5*alpha_n.*(u_bc + u) - 0.5*absvalue(alpha_n).*(u_bc - u)

  return bndryflux
end # end function RoeSolver


text:
"""

$(TYPEDEF)

Unit conversions.

$(TYPEDFIELDS)

"""
@with_kw struct Units
    mole = 6.022140857e23
    Angs3tocm3 = 1e24
    Angs3toL = 1e27
    Angs3tocm3permol = mole / Angs3tocm3
    Angs3toLpermol = mole / Angs3toL
    SitesperAngs3tomolperL = Angs3toL / mole
end
units = Units()