text:
module BrowseTables

export write_html_table, open_html_table, TableOptions

using ArgCheck: @argcheck
import DefaultApplication
using DocStringExtensions: SIGNATURES
import Tables

# options and customizations

const HTMLHEADSTART = """
<!DOCTYPE html>
  <meta content="text/html; charset=utf-8" http-equiv="Content-Type">
  <meta content="width=device-width, initial-scale=1, shrink-to-fit=no" name="viewport">
"""

const DEFAULTCSSPATH = abspath(@__DIR__, "..", "assets", "BrowseTables.css")

Base.@kwdef struct TableOptions
    css_path::AbstractString = DEFAULTCSSPATH
    css_inline::Bool = true
end

# low-level HTML construction

function writeescaped(io, str::AbstractString)
    for char in str
        if char == '&'
            write(io, "&")
        elseif char == '<'
            write(io, "<")
        elseif char == '>'
            write(io, ">")
        elseif char == '"'
            write(io, """)
        else
            write(io, char)
        end
    end
    nothing
end

escapestring(str::AbstractString) = (io = IOBuffer(); writeescaped(io, str); String(take!(io)))

function writetags(f, io::IO, tag::AbstractString; bropen = false, brclose = bropen,
                   attributes = NamedTuple())
    write(io, "<", tag)
    for (k, v) in pairs(attributes)
        write(io, " ", string(k), "=\"")
        writeescaped(io, v)
        write(io, "\"")
    end
    write(io, ">")
    bropen && println(io)
    f(io)
    write(io, "</", tag, ">")
    brclose && println(io)
    nothing
end

writecell(io, content::AbstractString; kind = "td", attributes = NamedTuple()) =
    writetags(io -> writeescaped(io, content), io, kind; attributes = attributes)

# cell formatting

"""
$(SIGNATURES)

HTML representation of a cell, with attributes. Content should be escaped properly, as it is
written to HTML directly.
"""
cellwithattributes(::TableOptions, x) = escapestring(string(x)), NamedTuple()

cellwithattributes(::TableOptions, x::Real) =
    escapestring(string(x)), isfinite(x) ? NamedTuple() : (class = "nonfinite", )

cellwithattributes(::TableOptions, str::AbstractString) =
    escapestring(str), (class = "alignleft", )

cellwithattributes(::TableOptions, x::Union{Missing,Nothing}) =
    escapestring(repr(x)), (class = "likemissing", )

struct RowId
    id::Int
end

cellwithattributes(::TableOptions, rowid::RowId) =
    escapestring(string(rowid.id)), (class = "rowid", )

# table structure

function writerow(io, options::TableOptions, nt::NamedTuple)
    writetags(io, "tr"; brclose = true) do io
        for x in values(nt)
            htmlcontent, attributes = cellwithattributes(options, x)
            writecell(io, htmlcontent; attributes = attributes)
        end
    end
end

writecaption(io, ::Nothing) = nothing

writecaption(io, str::AbstractString) = writetags(io -> writeescaped(io, str), io, "caption")

writeschema(io, ::Nothing) = nothing

function writeschema(io, sch::Tables.Schema)
    writetags(io, "thead"; brclose = true) do io
        writetags(io, "tr") do io
            writecell(io, "#"; kind = "th") # row id
            foreach(x -> writecell(io, escapestring(string(x)); kind = "th"), sch.names)
        end
    end
end

function writestyle(io::IO, path::AbstractString, inline::Bool)
    if inline
        writetags(io -> write(io, read(path, String)), io, "style")
    else
        println(io, raw"<link href=\"", path, raw"\" rel=\"stylesheet\" type=\"text/css\">")
    end
    nothing
end

# high-level API

function write_html_table(filename::AbstractString, table;
                          title = "Table", caption = nothing,
                          options::TableOptions = TableOptions())
    @argcheck Tables.istable(table) "The table should support the interface of Tables.jl."
    open(filename, "w") do io
        write(io, HTMLHEADSTART)
        writetags(io -> writeescaped(io, title), io, "title")
        writestyle(io, options.css_path, options.css_inline)
        println(io, "")    # close manually, opened in HTMLHEADSTART
        writetags(io, "body"; bropen = true) do io
            writetags(io, "table"; bropen = true) do io
                writecaption(io, caption)
                rows = Tables.rows(table)
                writeschema(io, Tables.schema(rows))
                writetags(io, "tbody"; bropen = true) do io
                    for (id, row) in enumerate(rows)
                        writerow(io, options, merge((rowid = RowId(id),), row))
                    end
                end
            end
        end
        println(io, "")
    end
    nothing
end

function open_html_table(table; filename = tempname() * ".html", kwargs...)
    write_html_table(filename, table; kwargs...)
    DefaultApplication.open(filename)
end

end # module


text:
#ubound-multiplication.jl

################################################################################
## multiplication

@universal function mul(a::Unum, b::Ubound)
  mul(b, a)
end

@universal function mul(a::Ubound, b::Unum)

  neg_b = is_negative(b)

  lower_result = resolve_lower((neg_b ? a.upper : a.lower) * b)
  upper_result = resolve_upper((neg_b ? a.lower : a.upper) * b)

  @ejectsolutions

  if is_ulp(lower_result) && is_ulp(upper_result)
    resolve_as_utype!(lower_result, upper_result)
  else
    B(lower_result, upper_result)
  end
end

@universal function mul(a::Ubound, b::Ubound)

  signcode::UInt16 = 0

  is_negative(a.lower) && (signcode += 1)
  is_negative(a.upper) && (signcode += 2)
  is_negative(b.lower) && (signcode += 4)
  is_negative(b.upper) && (signcode += 8)

  if (signcode == 0) #everything is positive

    lower_result = resolve_lower(a.lower * b.lower)
    upper_result = resolve_upper(a.upper * b.upper)

    @ejectsolutions

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 1) #only a.lowbound is negative
    lower_result = resolve_lower(a.lower * b.upper)
    upper_result = resolve_upper(a.upper * b.upper)

    @ejectsolutions

    B(lower_result, upper_result)
  #signcode 2 is not possible
  elseif (signcode == 3) #a is negative and b is positive
    lower_result = resolve_lower(a.lower * b.upper)
    upper_result = resolve_upper(a.upper * b.lower)

    @ejectsolutions

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 4) #only b.lowbound is negative
    lower_result = resolve_lower(b.lower * a.upper)
    upper_result = resolve_upper(b.upper * a.lower)

    @ejectsolutions

    B(lower_result, upper_result)
  elseif (signcode == 5) #a.lowbound and b.lowbound are negative
    minchoice1 = resolve_lower(b.lower * a.upper)
    minchoice2 = resolve_lower(b.upper * a.lower)
    maxchoice1 = resolve_upper(b.lower * a.lower)
    maxchoice2 = resolve_upper(b.upper * a.upper)

    (isnan(minchoice1) || isnan(minchoice2) || isnan(maxchoice1) || isnan(maxchoice2)) && return nan(U)

    B(min(minchoice1, minchoice2), max(maxchoice1, maxchoice2))
  #signcode 6 is not possible
  elseif (signcode == 7) #only b.highbound is positive
    lower_result = resolve_lower(b.upper * a.lower)
    upper_result = resolve_upper(b.lower * a.lower)

    @ejectsolutions

    B(lower_result, upper_result)
  #signcode 8, 9, 10, 11 are not possible
  elseif (signcode == 12) #b is negative, a is positive
    lower_result = resolve_lower(b.lower * a.upper)
    upper_result = resolve_upper(b.upper * a.lower)

    @ejectsolutions

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 13) #b is negative, a straddles
    lower_result = resolve_lower(b.lower * a.upper)
    upper_result = resolve_upper(b.lower * a.lower)

    @ejectsolutions

    B(lower_result, upper_result)
  #signcode 14 is not possible
  elseif (signcode == 15) #everything is negative
    lower_result = resolve_lower(a.upper * b.upper)
    upper_result = resolve_upper(a.lower * b.lower)

    @ejectsolutions

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  else
    throw(ArgumentError("error multiplying ubounds $a and $b, throws invalid signcode $signcode."))
  end
end

ext: jl

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module test_julia_types

using Test

@test Base.uniontypes(Int) == [Int]
@test Base.uniontypes(Union{Int,String}) == [Int,String]

@test typeof(Type{String}) === DataType
@test typeof(Type{<:AbstractString}) === UnionAll
@test typeof(String) === DataType

@test typeof(UnionAll) === DataType

@test typeof(Union) === DataType
@test typeof(Union{}) === Core.TypeofBottom
@test Union{} === Base.Bottom
@test Base.Bottom isa Core.TypeofBottom
@test typeof(==) === typeof(==)
@test Union{Int,Int} === Union{Int} === Int
@test typeintersect(Int, Bool) === Union{}

@test typeof(Vector{}) === typeof(Vector) === UnionAll
@test Float64 isa Type{Float64}

# https://docs.julialang.org/en/v1/devdocs/types/
TV, NV = TypeVar(:T), TypeVar(:N)
if VERSION >= v"1.6.0-DEV.849"
    @test string(Array{TV,NV}) == "Array{T, N}"
end

@test Any === (T where T) === (Int where Int) === (where where where)
@test (+)::typeof(+) === +

end # module test_julia_types

module test_julia_datatype

using Test
using InteractiveUtils # subtypes

@test typeintersect(DataType, Type) === DataType

@test supertype(DataType) === Type.body
@test subtypes(DataType) == []

@test Type <: Type
@test Type{Int} <: Type
@test Type{<:Integer} <: Type
@test Type{Int} <: Type{<:Integer}
@test Type{Signed} <: Type{<:Integer}
@test !(Int <: Type{<:Integer})

@test supertype(Type) === Any
@test subtypes(Type) == [Core.TypeofBottom, DataType, Union, UnionAll]

end # module test_julia_datatype

module test_julia_types_eltype

using Test

abstract type AbstractT{V} end

struct T{V, A} <: AbstractT{V}
    data::Vector{V}
    T(data::A) where {A} = new{eltype(A), A}(data)
end

@test T([1, 2, 3]) isa T{Int,Vector{Int}}
@test T(['a']) isa T{Char,Vector{Char}}

end # module test_julia_types_eltype

module test_julia_types_isimmutable

using Test

@test isimmutable(0)
@test isimmutable(pi)
@test isimmutable((2,))
@test isimmutable((a=2,))
@test isimmutable(:a=>2)
@test isimmutable(nothing)

# mutable types
@test !isimmutable(:abc)
@test !isimmutable("abc")
@test !isimmutable([])
@test !isimmutable(Dict())

mutable struct A
end
@test !isimmutable(A())

end # module test_julia_types_isimmutable

module test_julia_types_tuple

using Test

@test Tuple{Int,Int} <: Tuple
@test Tuple{Int,Int} <: Tuple{T,T} where T
@test !(Tuple{Real,Real} <: Tuple{T,T} where T) # diagonal variables
@test Tuple{Real} <: Tuple{T} where T # not diagonal
@test Tuple{Int,Int} <: (Tuple{T,T} where T) <: Tuple
@test NTuple{2,Int} <: Tuple{T,T} where T
@test NTuple{2,Int} <: Tuple{T,T} where T
@test !(Tuple{Int,String} <: Tuple{T,T} where T)

@test (Tuple{T} where T) <: Tuple{Vararg{T} where T}
@test (Tuple{T,T} where T) <: Tuple{Vararg{T}} where T
@test Tuple{Int,String} <: Tuple{Vararg{T} where T}

f(x::T, y::T) where {T} = 0
@test methods(f, Tuple{T,T} where T).mt.name === :f

end # module test_julia_types_tuple

module test_julia_types_union

using Test

f(x::Int, y::Int) = :int
f(x::Union{Int,String}, y::Union{Int,String}) = :union

@test f(1, 2) === :int
@test f(1, "") === :union

@test ((+)::typeof(+))(1, 2) == 3

@test Base.uniontypes(Union{Int, String}) == [Int, String]

end # module test_julia_types_union

module test_julia_types_param

using Test

abstract type Result end

struct OK <: Result
    result
end

struct Fail <: Result
    result
end

struct A{S}
end

f(::A{:a}) = OK(:a)
f(::A{:b}) = OK(:b)

@test f(A{:a}()) == OK(:a)
@test f(A{:b}()) == OK(:b)

f(::A{X}) where X = Fail(X)

@test f(A{:c}()) == Fail(:c)

end # module test_julia_types_param

module test_julia_types_struct

using Test
using ..test_julia_types_param

struct A
    a::test_julia_types_param.OK
end

struct B
    f::typeof(push!)
end

@test B(push!).f === push!

end # module test_julia_types_struct

module test_julia_types_isstructtype

using Test

@test !isstructtype(Int)
@test isstructtype(Vector)

end # module test_julia_types_isstructtype

module test_julia_types_subtyping

using Test

abstract type AbstractT end

const AT = AbstractT

struct ConcretT <: AT
end

@test supertype(ConcretT) === AbstractT
@test supertype(typeof(ConcretT())) === AbstractT

end # module test_julia_types_subtyping


text:
export ModelSymbol

const doc_macro = MacroTools.unblock(quote
    "hello"
    world
end).args[1]

struct ModelSymbol
    doc::String
    name::Symbol
    type::Symbol
    ModelSymbol(d::String, n::Symbol, t::Symbol) =
                t ∉ (:steady, :lin, :log, :shock) ?
                    throw(ArgumentError("Invalid symbol type $t.")) :
                    new(d, n, t)
end

ModelSymbol(s::Symbol) = ModelSymbol("", s, :lin)
ModelSymbol(d::String, s::Symbol) = ModelSymbol(d, s, :lin)
ModelSymbol(s::Symbol, t::Symbol) = ModelSymbol("", s, t)

function ModelSymbol(s::Expr)
    s = MacroTools.unblock(s)
    if MacroTools.isexpr(s, :macrocall) && s.args[1] == doc_macro
        return ModelSymbol(s.args[3], s.args[4])
    else
        return ModelSymbol("", s)
    end
end

function ModelSymbol(doc::String, s::Expr)
    s = MacroTools.unblock(s)
    if MacroTools.isexpr(s, :macrocall)
        t = Symbol(String(s.args[1])[2:end])
        return ModelSymbol(doc, s.args[3], t)
    else
        throw(ArgumentError("Invalid variable or shock expression $s."))
    end
end

for sym ∈ (:shock, :log, :lin, :steady)
    to_sym = Symbol("to_", sym)
    issym = Symbol("is", sym)
    eval(quote
        @inline $(to_sym)(s::ModelSymbol) = ModelSymbol(s.doc, s.name, $(QuoteNode(sym)))
        @inline $(to_sym)(any) = $(to_sym)(convert(ModelSymbol, any))
        export $(to_sym)
        @inline $(issym)(any) = false
        @inline $(issym)(s::ModelSymbol) = s.type == $(QuoteNode(sym))
        export $(issym)
    end)
end

Symbol(s::ModelSymbol) = s.name
Base.convert(::Type{Symbol}, v::ModelSymbol) = v.name
Base.convert(::Type{ModelSymbol}, v::Symbol) = ModelSymbol(v)
Base.convert(::Type{ModelSymbol}, v::Expr) = ModelSymbol(v)
Base.:(==)(a::ModelSymbol, b::ModelSymbol) = a.name == b.name
Base.:(==)(a::ModelSymbol, b::Symbol) = a.name == b
Base.:(==)(a::Symbol, b::ModelSymbol) = a == b.name

# The hash must be the same as the hash of the symbol, so that we can use
# ModelSymbol as index in a Dict with Symbol keys
Base.hash(v::ModelSymbol, h::UInt) = hash(v.name, h)

function Base.show(io::IO, v::ModelSymbol)
    if get(io, :compact, false)
        print(io, v.name)
    else
        doc = isempty(v.doc) ? "" : "\"$(v.doc)\" "
        type = v.type ∈ (:lin, :shock) ? "" : "@$(v.type) "
        print(io, doc, type, v.name)
    end
end


text:
crawl(h::Val{:function}, signature, body) = :(
    JSAST(
        :function,
        $(_crawl_function_signature(signature)),
        $(crawl(body)),
    )
)
function deparse(::Val{:function}, signature, body)
    body = jsstring(
        "function ", deparse(signature), " ",
        deparse(body),
    )
    # We need to wrap anonymous functions in parens (otherwise, it's a syntax
    # error in JavaScript).
    if signature.args[1] == JSTerminal(Symbol(""))
        return jsstring("(", body, ")")
    end
    return body
end

"""
Crawl the signature of a function.

In JS,
```js
function (...args) { /* ... */ }
```
is valid syntax, but in Julia, since there's no function name, the function
signature gets parsed as a `:tuple` instead of a `:call`. Since we convert
tuples to arrays in JS, we need to make sure the function signature actually is
a `:call` expression before crawling it.
"""
function _crawl_function_signature(signature::Expr)
    if signature.head == :call
        return crawl(signature)
    elseif signature.head == :tuple
        return crawl(Expr(:call, Symbol(""), signature.args...))
    end
    error("Invalid function signature (expected :call or :tuple Expr).")
end
function _crawl_function_signature(signature)
    error("Invalid JS function signature: $(signature).")
end

# NOTE: We need to treat arrow functions differently than normal functions
# because of the way `this` binding works in JS using arrow functions (i.e.,
# its invalid to transform () -> ... into JavaScript that uses the function
# keyword instead of the => JavaScript syntax).
function crawl(::Val{:(->)}, arguments, body)
    return :(JSAST(
        :arrowfunction,
        $(_crawl_arrow_signature(arguments)),
        $(crawl(body)),
    ))
end
function deparse(::Val{:arrowfunction}, signature, body)
    return jsstring(
        deparse(signature), " => ",
        _deparse_returnify_block(body),
    )
end

function crawl(h::Val{:return}, ex)
    return :(JSAST(
        :return,
        $(crawl(ex)),
    ))
end
function deparse(h::Val{:return}, ex)
    return jsstring("return ", deparse(ex))
end

function crawl(::Val{:block}, body...)
    return :(JSAST(
        :block,
        $(crawl.(body)...),
    ))
end
function deparse(::Val{:block}, body...)
    if length(body) == 0
        return jsstring("{}")
    end
    return jsstring("{ ", join(deparse.(body), "; "), "; }")
end

function returnify_js_statement(statement::JSNode)
    if !isa(statement, JSAST)
        return JSAST(:return, statement)
    end

    head = statement.head
    if head == :return
        return statement
    elseif head == :if
        statement.args[2] = returnify_js_statement(statement.args[2])
        if length(statement.args) > 2
            statement.args[3] = returnify_js_statement(statement.args[3])
        end
        return statement
    elseif head == :block
        statement.args[end] = returnify_js_statement(statement.args[end])
        return statement
    elseif in(head, (:while, :for))
        # We give up trying to add implicit returns if it's "too hard".
        return statement
    end
    return JSAST(:return, statement)
end

"""
Deparse a block, transforming its final statement into a return statement.

This ensures that Julia semantics wherein the last expression is implicitly the
return value is translated to the generated JavaScript code (as is necessary
for arrow functions to work as expected).
"""
function _deparse_returnify_block(block::JSNode)
    if !(isa(block, JSAST) && block.head == :block)
        @show block
        error("Argument to returnify_block must be a JSAST(:block, ...).")
    end
    if isempty(block.args)
        return deparse(block)
    end
    block.args[end] = returnify_js_statement(block.args[end])
    return deparse(block)
end

"""
    _crawl_arrow_signature(expr)

Convert arrow function arguments to a `JSAST` with head `call`.
This works by turning a tuple of arguments (or a symbol in the case of a single
argument) into a `call` whose function name is an empty string (which matches
the javascript syntax where the function signature looks like a call but with
an optional function name).
"""
_crawl_arrow_signature(ex::Expr) = _crawl_function_signature(ex)
_crawl_arrow_signature(s::Symbol) = _crawl_function_signature(:($s, ))


text:
module Conductivity

# This implements different methods to specify conductivity of rocks
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam, AbstractMaterialParamsStruct
import Base.show

abstract type AbstractConductivity <: AbstractMaterialParam end

export  ComputeConductivity,                # calculation routines
        ComputeConductivity!,
        ConstantConductivity,               # constant
        T_Conductivity_Whittacker,          # T-dependent heat capacity
        TP_Conductivity,                    # TP dependent conductivity
        Set_TP_Conductivity                 # Routine to set pre-defined parameters

# Constant Conductivity -------------------------------------------------------
"""
    ConstantConductivity(k=3.0W/m/K)

Set a constant conductivity
```math
    k  = cst
```
where ``k`` is the thermal conductivity [``W/m/K``].
"""
@with_kw_noshow mutable struct ConstantConductivity <: AbstractConductivity
    equation::LaTeXString   =   L"k = cst"
    k::GeoUnit              =   3.0Watt/m/K
end

# Calculation routine
function ComputeConductivity(P, T, s::ConstantConductivity)
    @unpack k   = s

    if length(T)>1
        return Value(k).*ones(size(T))
    else
        return Value(k)
    end
end

"""
    ComputeConductivity(k_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantConductivity) where N

In-place routine to compute constant conductivity
"""
function ComputeConductivity!(k_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::ConstantConductivity) where N
    @unpack k   = s

    k_array .= ustrip(Value(k))

    return nothing
end

# Print info
function show(io::IO, g::ConstantConductivity)
    print(io, "Constant conductivity: k=$(g.k.val)")
end
#-------------------------------------------------------------------------

# Temperature dependent conductivity -------------------------------
"""
    T_Conductivity_Whittacker()

Sets a temperature-dependent conductivity following the parameterization of *Whittacker et al. (2009), Nature.*
Their parameterization is originally given for the thermal diffusivity, together with a parameterization for thermal conductivity, which allows us
```math
    Cp = a + b T - c/T^2
```
```math
    \\kappa = d/T - e, if T<=846K
```
```math
    \\kappa = f - g*T, if T>846K
```
```math
    \\rho = 2700 kg/m3
```
```math
    k = \\kappa Cp \\rho
```

where ``Cp`` is the heat capacity [``J/mol/K``], and ``a,b,c`` are parameters that dependent on the temperature `T`:
- a = 199.50 J/mol/K    if T<= 846 K
- a = 199.50 J/mol/K    if T> 846 K
- b = 0.0857J/mol/K^2   if T<= 846 K
- b = 0.0323J/mol/K^2   if T> 846 K
- c = 5e6J/mol*K        if T<= 846 K
- c = 47.9e-6J/mol*K    if T> 846 K
- d = 576.3m^2/s*K
- e = 0.062m^2/s
- f = 0.732m^2/s
- g = 0.000135m^2/s/K
"""
@with_kw_noshow mutable struct T_Conductivity_Whittacker <: AbstractConductivity
    # Note: the resulting curve of k was visually compared with Fig. 2 of the paper
    equation::LaTeXString   =   L"k = f(T) "
    a0::GeoUnit             =   199.5J/mol/K                # prefactor for low T       (T<= 846 K)
    a1::GeoUnit             =   229.32J/mol/K               # prefactor for high T      (T>  846 K)
    b0::GeoUnit             =   0.0857J/mol/K^2             # linear term for low T     (T<= 846 K)
    b1::GeoUnit             =   0.0323J/mol/K^2             # linear term for high T    (T>  846 K)
    c0::GeoUnit             =   5e6J/mol*K                  # quadratic term for low T  (T<= 846 K)
    c1::GeoUnit             =   47.9e-6J/mol*K              # quadratic term for high T (T>  846 K)
    molmass::GeoUnit        =   0.22178kg/mol               # average molar mass
    Tcutoff::GeoUnit        =   846K                        # cutoff temperature
    rho::GeoUnit            =   2700kg/m^3                  # Density they use for an average crust
    d::GeoUnit              =   576.3*1e-6m^2/s*K           # diffusivity parameterization
    e::GeoUnit              =   0.062*1e-6m^2/s             # diffusivity parameterization
    f::GeoUnit              =   0.732*1e-6m^2/s             # diffusivity parameterization
    g::GeoUnit              =   0.000135*1e-6m^2/s/K        # diffusivity parameterization
end

# Calculation routine
function ComputeConductivity(P,T, s::T_Conductivity_Whittacker)
    @unpack a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g   = s

    ρ  = Value(rho)
    k  = zeros(size(T))*   Value(a0)/Value(molmass)*ρ*Value(e)  # the last multiplication ensures the correct units even for non-dimensional cases

    for i in eachindex(T)
        if T[i] <= Value(Tcutoff)
            a,b,c = Value(a0),Value(b0),Value(c0)
            κ     = Value(d)/T[i] - Value(e)
        else
            a,b,c = Value(a1),Value(b1),Value(c1)
            κ     = Value(f) - Value(g)*T[i]
        end

        cp = (a + b*T[i] - c/T[i]^2)/molmass # conductivity

        k[i] = κ*ρ*cp       # compute conductivity from diffusivity

    end

    return k
end

"""
    ComputeConductivity(k_array::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::T_Conductivity_Whittacker) where N

In-place routine to compute temperature-dependent conductivity
"""
function ComputeConductivity!(k::AbstractArray{<:AbstractFloat,N},P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, s::T_Conductivity_Whittacker) where N
    @unpack a0,a1,b0,b1,c0,c1,molmass,Tcutoff,rho,d,e,f,g   = s
    a0,b0,c0    =   NumValue(a0),NumValue(b0),NumValue(c0)
    a1,b1,c1    =   NumValue(a1),NumValue(b1),NumValue(c1)
    ρ           =   NumValue(rho)
    d,e,f,g     =   NumValue(d),NumValue(e),NumValue(f), NumValue(g)
    Tcutoff     =   NumValue(Tcutoff)
    molmass     =   NumValue(molmass)

    ind         =   (T .<= Tcutoff)
    T_local     =   view(T, ind )
    k[ind]      =   (a0 .+ b0*T_local - c0./T_local.^2)/molmass .* (d./T_local .- e) .* ρ

    ind         =   (T .> Tcutoff)
    T_local     =   view(T, ind )
    k[ind]      =   (a1 .+ b1*T_local - c1./T_local.^2)/molmass  .* (f .- g.*T_local ) .* ρ

    return nothing
end

# Print info
function show(io::IO, g::T_Conductivity_Whittacker)
    print(io, "T-dependent conductivity following Whittacker et al. (2009) for average crust). \n");
end
#-------------------------------------------------------------------------

# Temperature (& Pressure) dependent conductivity -------------------------------
"""
    TP_Conductivity()

Sets a temperature (and pressure)-dependent conductivity parameterization as described in Gerya, Numerical Geodynamics (2nd edition, Table 21.2).
The general for

```math
    k = \\left( a_k +  {b_k \\over {T + c_k}} \\right) (1 + d_k P)
```

where ``k`` is the conductivity [``W/K/m``], and ``a_k,b_k,c_k,d_k`` are parameters that dependent on the temperature `T` and pressure `P`:
- ``a_k`` = 1.18Watt/K/m
- ``b_k`` = 474Watt/m
- ``c_k`` = 77K
- ``d_k`` = 0/MPa
"""
@with_kw_noshow mutable struct TP_Conductivity <: AbstractConductivity
    equation::LaTeXString   =   L"k = \left(a_k + {b_k/{T + c_k}} \right)*(1 + d_k*P) "
    a::GeoUnit              =   1.18Watt/K/m        # empirical fitting term
    b::GeoUnit              =   474Watt/m           # empirical fitting term
    c::GeoUnit              =   77K                 # empirical fitting term
    d::GeoUnit              =   0/MPa               # empirical fitting term
    Comment::String         =   ""                  # Some remarks you want to add about this creep law implementation
    BibTex_Reference        =   ""                  # BibTeX reference
end

"""
    Set_TP_Conductivity["Name of temperature(-pressure) dependent conductivity"]

This is a dictionary with pre-defined laws:
- "UpperCrust"
- "LowerCrust"
- "OceanicCrust"
- "Mantle"

# Example
```julia
julia> k=Set_TP_Conductivity["Mantle"]
T/P dependent conductivity: k = (0.73 W K⁻¹ m⁻¹ + 1293 W m⁻¹/(T + 77 K))*(1 + 4.0e-5 MPa⁻¹*P)
```

"""
Set_TP_Conductivity = Dict([
    ("UpperCrust",
        TP_Conductivity( a=0.64Watt/K/m, b=807Watt/m, c=77K, d=0/MPa,
            Comment="Sediment/upper crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!")
    )

    ("LowerCrust",
        TP_Conductivity( a=1.18Watt/K/m, b=474Watt/m, c=77K, d=0/MPa,
            Comment="Lower crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!")
    )

    ("OceanicCrust",
        TP_Conductivity( a=1.18Watt/K/m, b=474Watt/m, c=77K, d=0/MPa,
            Comment="Oceanic crust T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!")
    )

    ("Mantle",
        TP_Conductivity( a=0.73Watt/K/m, b=1293Watt/m, c=77K, d=0.00004/MPa,
            Comment="Mantle T-dependent conductivity, as listed in table 21.2 of Gerya et al. | Reference still to be verified!")
    )

])

# Calculation routine
function ComputeConductivity(P,T, s::TP_Conductivity)
    @unpack a,b,c,d   = s

    a_k, b_k = Value(a), Value(b)
    c_k, d_k = Value(c), Value(d)

    k  = zeros(size(T))*   a_k  # the last multiplication ensures the correct units even for non-dimensional cases

    if ustrip(d_k)==0
        for i in eachindex(T)
            k[i] = a_k + b_k/(T[i] + c_k)
        end
    else
        if size(T) != size(P)
            error("Size of P and T arrays should be the same")
        end

        for i in eachindex(T)
            k[i] = (a_k + b_k/(T[i] + c_k))*(1 + d_k*P[i])
        end
    end

    return k
end

# Calculation routine
function ComputeConductivity!(K::AbstractArray{T, N}, P::AbstractArray{T, N},Temp::AbstractArray{T, N}, s::TP_Conductivity) where{T<:AbstractFloat, N}
    @unpack a,b,c,d   = s

    a_k, b_k = NumValue(a), NumValue(b)
    c_k, d_k = NumValue(c), NumValue(d)

    if d_k==0
        K[:] = a_k .+ b_k./(Temp .+ c_k)
    else
        K[:] = (a_k .+ b_k./(Temp .+ c_k)).*(1.0 .+ d_k.*P)
    end

    return nothing
end

# Print info
function show(io::IO, g::TP_Conductivity)
    if ustrip(Value(g.d))==0
        print(io, "T/P dependent conductivity: k = $(g.a.val) + $(g.b.val)/(T + $(g.c.val))  \n");
    else
        print(io, "T/P dependent conductivity: k = ($(g.a.val) + $(g.b.val)/(T + $(g.c.val)))*(1 + $(g.d.val)*P)  \n");
    end
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
    k = ComputeConductivity(P, T, s:<AbstractConductivity)

Returns the thermal conductivity `k` at any temperature `T` and pressure `P` using any of the parameterizations implemented.

Currently available:
- ConstantConductivity
- T\\_Conductivity_Whittacker
- TP\\_Conductivity

# Example
Using dimensional units
```julia
julia> T  = (250:100:1250)*K;
julia> cp = T_HeatCapacity_Whittacker()
julia> Cp = ComputeHeatCapacity(0,T,cp)
```

"""
ComputeConductivity()

"""
    k = ComputeConductivity(T::Any, s::AbstractConductivity)

Computes conductivity if only temperature (and not pressure) is specified
"""
ComputeConductivity(T::Any, s::AbstractConductivity) =  ComputeConductivity(0,T, s)

"""
    k = ComputeConductivity(s::ConstantConductivity)

Returns conductivity if we are sure that we will only employ constant values throughout the simulation
"""
ComputeConductivity(s::ConstantConductivity) =  ComputeConductivity(0,0, s)

"""
    ComputeConductivity!(K::AbstractArray{<:AbstractFloat}, Phases::AbstractArray{<:Integer}, P::AbstractArray{<:AbstractFloat},Temp::AbstractArray{<:AbstractFloat}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of conductivity `K` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `Temp` arrays.
This assumes that the `Phase` of every point is specified as an Integer in the `Phases` array.

"""
function ComputeConductivity!(K::AbstractArray{T, N}, Phases::AbstractArray{<:Integer, N}, P::AbstractArray{T, N},Temp::AbstractArray{T, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {T<:AbstractFloat,N}

    for i = 1:length(MatParam)

        if !isnothing(MatParam[i].Conductivity)
            # Create views into arrays (so we don't have to allocate)
            ind = Phases .== MatParam[i].Phase;
            K_local     =   view(K   , ind )
            P_local     =   view(P   , ind )
            T_local     =   view(Temp, ind )

            ComputeConductivity!(K_local, P_local, T_local, MatParam[i].Conductivity[1] )
        end

    end

end

"""
    ComputeConductivity!(k::AbstractArray{T,N}, PhaseRatios::AbstractArray{T, M}, P::AbstractArray{<:AbstractFloat,N},T::AbstractArray{<:AbstractFloat,N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct})

In-place computation of density `rho` for the whole domain and all phases, in case a vector with phase properties `MatParam` is provided, along with `P` and `T` arrays.
This assumes that the `PhaseRatio` of every point is specified as an Integer in the `PhaseRatios` array, which has one dimension more than the data arrays (and has a phase fraction between 0-1)

"""
function ComputeConductivity!(k::AbstractArray{T, N}, PhaseRatios::AbstractArray{T, M}, P::AbstractArray{T, N},Temp::AbstractArray{T, N}, MatParam::AbstractArray{<:AbstractMaterialParamsStruct, 1}) where {T<:AbstractFloat, N,M}

    if M!=(N+1)
        error("The PhaseRatios array should have one dimension more than the other arrays")
    end

    k .= 0.0;
    k_local     = zeros(size(k))
    for i = 1:length(MatParam)
        k_local .= 0.0
        Fraction    = selectdim(PhaseRatios,M,i);
        if (maximum(Fraction)>0.0) & (!isnothing(MatParam[i].Conductivity))

            ComputeConductivity!(k_local, P, Temp, MatParam[i].Conductivity[1] )

            k .= k .+ k_local.*Fraction
        end

    end

end

end


text:
# Definition of the abstract Domain type and its interface

# The type Domain{T} is defined in IntervalSets.jl

eltype(::Type{<:Domain{T}}) where {T} = T
prectype(::Type{<:Domain{T}}) where {T} = prectype(T)
numtype(::Type{<:Domain{T}}) where {T} = numtype(T)

convert_numtype(d::Domain{T}, ::Type{U}) where {T,U} = convert(Domain{to_numtype(T,U)}, d)
convert_prectype(d::Domain{T}, ::Type{U}) where {T,U} = convert(Domain{to_prectype(T,U)}, d)

Domain(d) = convert(Domain, d)

# Concrete types can implement similardomain(d, ::Type{T}) where {T}
# to support convert(Domain{T}, d) functionality.
convert(::Type{Domain{T}}, d::Domain{T}) where {T} = d
convert(::Type{Domain{T}}, d::Domain{S}) where {S,T} = similardomain(d, T)

"Can the domains be promoted without throwing an error?"
promotable_domains(d1, d2) = promotable_eltypes(eltype(d1),eltype(d2))
promotable_eltypes(::Type{S}, ::Type{T}) where {S,T} =
    isconcretetype(promote_type(S, T))
promotable_eltypes(::Type{S}, ::Type{T}) where {S<:AbstractVector,T<:AbstractVector} =
    promotable_eltypes(eltype(S), eltype(T))

"Promote the given domains to have a common element type."
promote_domains() = ()
promote_domains(domains...) = promote_domains(domains)
promote_domains(domains) = convert_eltype.(mapreduce(eltype, promote_type, domains), domains)

promote_domains(domains::AbstractSet{<:Domain{T}}) where {T} = domains
promote_domains(domains::AbstractSet{<:Domain}) = Set(promote_domains(collect(domains)))

convert_eltype(::Type{T}, d::Domain) where {T} = convert(Domain{T}, d)
convert_eltype(::Type{T}, d) where {T} = _convert_eltype(T, d, eltype(d))
_convert_eltype(::Type{T}, d, ::Type{T}) where {T} = d
_convert_eltype(::Type{T}, d, ::Type{S}) where {S,T} =
    error("Don't know how to convert the `eltype` of $(d).")
# Some standard cases
convert_eltype(::Type{T}, d::AbstractArray) where {T} = convert(AbstractArray{T}, d)
convert_eltype(::Type{T}, d::AbstractRange) where {T} = map(T, d)
convert_eltype(::Type{T}, d::Set) where {T} = convert(Set{T}, d)

promote(d1::Domain, d2::Domain) = promote_domains((d1,d2))
promote(d1::Domain, d2) = promote_domains((d1,d2))
promote(d1, d2::Domain) = promote_domains((d1,d2))

"A `EuclideanDomain` is any domain whose eltype is `<:StaticVector{N,T}`."
const EuclideanDomain{N,T} = Domain{<:StaticVector{N,T}}

"A `VectorDomain` is any domain whose eltype is `Vector{T}`."
const VectorDomain{T} = Domain{Vector{T}}

const AbstractVectorDomain{T} = Domain{<:AbstractVector{T}}

CompositeTypes.Display.displaysymbol(d::Domain) = 'D'

"What is the Euclidean dimension of the domain?"
dimension(::Domain{T}) where {T} = euclideandimension(T)

"Is the given combination of point and domain compatible?"
iscompatiblepair(x, d) = _iscompatiblepair(x, d, typeof(x), eltype(d))
_iscompatiblepair(x, d, ::Type{S}, ::Type{T}) where {S,T} =
    _iscompatiblepair(x, d, S, T, promote_type(S,T))
_iscompatiblepair(x, d, ::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U} = true
_iscompatiblepair(x, d, ::Type{S}, ::Type{T}, ::Type{Any}) where {S,T} = false
_iscompatiblepair(x, d, ::Type{S}, ::Type{Any}, ::Type{Any}) where {S} = true

# Some generic cases where we can be sure:
iscompatiblepair(x::SVector{N}, ::EuclideanDomain{N}) where {N} = true
iscompatiblepair(x::SVector{N}, ::EuclideanDomain{M}) where {N,M} = false
iscompatiblepair(x::AbstractVector, ::EuclideanDomain{N}) where {N} = length(x)==N

# Note: there are cases where this warning reveals a bug, and cases where it is
# annoying. In cases where it is annoying, the domain may want to specialize `in`.
compatible_or_false(x, domain) =
    iscompatiblepair(x, domain) ? true : (@warn "`in`: incompatible combination of point: $(typeof(x)) and domain eltype: $(eltype(domain)). Returning false."; false)

compatible_or_false(x::AbstractVector, domain::AbstractVectorDomain) =
    iscompatiblepair(x, domain) ? true : (@warn "`in`: incompatible combination of vector with length $(length(x)) and domain '$(domain)' with dimension $(dimension(domain)). Returning false."; false)

"Promote point and domain to compatible types."
promote_pair(x, d) = _promote_pair(x, d, promote_type(typeof(x),eltype(d)))
_promote_pair(x, d, ::Type{T}) where {T} = convert(T, x), convert(Domain{T}, d)
_promote_pair(x, d, ::Type{Any}) = x, d
# Some exceptions:
# - matching types: avoid promotion just in case it is expensive
promote_pair(x::T, d::Domain{T}) where {T} = x, d
# - `Any` domain
promote_pair(x, d::Domain{Any}) = x, d
# - tuples: these are typically composite domains and the elements may be promoted later on
promote_pair(x::Tuple, d::Domain{<:Tuple}) = x, d
# - abstract vectors: promotion may be expensive
promote_pair(x::AbstractVector, d::AbstractVectorDomain) = x, d
# - SVector: promotion is likely cheap
promote_pair(x::AbstractVector{S}, d::EuclideanDomain{N,T}) where {N,S,T} =
    _promote_pair(x, d, SVector{N,promote_type(S,T)})

# At the level of Domain we attempt to promote the arguments to compatible
# types, then we invoke indomain. Concrete subtypes should implement
# indomain. They may also implement `in` in order to accept more types.
#
# Note that if the type of x and the element type of d don't match, then
# both x and the domain may be promoted (using convert(Domain{T}, d)) syntax).
in(x, d::Domain) = compatible_or_false(x, d) && indomain(promote_pair(x, d)...)

"""
Return a suitable tolerance to use for verifying whether a point is close to
a domain. Typically, the tolerance is close to the precision limit of the numeric
type associated with the domain.
"""
default_tolerance(d::Domain) = default_tolerance(prectype(d))
default_tolerance(::Type{T}) where {T <: AbstractFloat} = 100eps(T)

"""
`approx_in(x, domain::Domain [, tolerance])`

Verify whether a point lies in the given domain with a certain tolerance.

The tolerance has to be positive. The meaning of the tolerance, in relation
to the possible distance of the point to the domain, is domain-dependent.
Usually, if the outcome is true, it means that the distance of the point to
the domain is smaller than a constant times the tolerance. That constant may
depend on the domain.

Up to inexact computations due to floating point numbers, it should also be
the case that `approx_in(x, d, 0) == in(x,d)`. This implies that `approx_in`
reflects whether a domain is open or closed.
"""
approx_in(x, d::Domain) = approx_in(x, d, default_tolerance(d))

function compatible_or_false(x, d, tol)
    tol >= 0 || error("Tolerance has to be positive in `approx_in`.")
    compatible_or_false(x, d)
end

approx_in(x, d::Domain, tol) =
    compatible_or_false(x, d, tol) && approx_indomain(promote_pair(x, d)..., tol)

# Fallback to `in`
approx_indomain(x, d::Domain, tol) = in(x, d)

isapprox(d1::Domain, d2::Domain; kwds...) = d1 == d2

isreal(d::Domain) = isreal(eltype(d))

infimum(d::Domain) = minimum(d)  # if the minimum exists, then it is also the infimum
supremum(d::Domain) = maximum(d)  # if the maximum exists, then it is also the supremum

"""
Return a bounding box of the given domain.

A bounding box is an interval, a hyperrectangle or the full space. It is such that
each point in the domain also lies in the bounding box.
"""
boundingbox(d) = FullSpace{eltype(d)}()

"Return the boundary of the given domain as a domain."
function boundary end
const ∂ = boundary

"""
Return the normal of the domain at the point `x`.

It is assumed that `x` is a point on the boundary of the domain.
"""
function normal end

"""
Return the tangents of the domain at the point `x`. The tangents form a
basis for the tangent plane, perpendicular to the normal direction at `x`.
"""
function tangents end

# "Lazy representation of the boundary of a domain."
# struct Boundary{T,D} <: Domain{T}
#     domain  :: D
# end
#
# domain(d::Boundary) = d.domain
#
# Boundary(domain) = Boundary{eltype(domain)}(domain)
# Boundary{T}(domain) where {T} = Boundary{T,typeof(domain)}(domain)
# Boundary{T}(domain::Domain{T}) where {T} = Boundary{T,typeof(domain)}(domain)
# Boundary{T}(domain::Domain{S}) where {S,T} = Boundary{T}(convert(Domain{T}, domain))
#
# indomain(x, d::Boundary) = in(x, boundary(domain(d)))


text:
using Tensorial: fromvoigt, ⊡, ⋅

import Tensorial: contraction, double_contraction

export contraction, double_contraction

Base.inv(c::StiffnessTensor) = ComplianceTensor(inv(c.data))
Base.inv(s::ComplianceTensor) = StiffnessTensor(inv(s.data))
Base.inv(c::StiffnessMatrix) = ComplianceMatrix(inv(c.data))
Base.inv(s::ComplianceMatrix) = StiffnessMatrix(inv(s.data))

Base.convert(::Type{TensorStress{T}}, s::EngineeringStress{T}) where {T} =
    TensorStress(s[1], s[6], s[5], s[2], s[4], s[3])
Base.convert(::Type{EngineeringStress{T}}, σ::TensorStress{T}) where {T} =
    EngineeringStress(σ[1, 1], σ[2, 2], σ[3, 3], σ[2, 3], σ[1, 3], σ[1, 2])
Base.convert(::Type{TensorStrain{T}}, ϵ::EngineeringStrain{T}) where {T} =
    TensorStrain(ϵ[1], ϵ[6] / 2, ϵ[5] / 2, ϵ[2], ϵ[4] / 2, ϵ[3])
Base.convert(::Type{EngineeringStrain{T}}, ε::TensorStrain{T}) where {T} =
    EngineeringStrain(ε[1, 1], ε[2, 2], ε[3, 3], 2ε[2, 3], 2ε[1, 3], 2ε[1, 2])
Base.convert(::Type{StiffnessMatrix{T}}, c::StiffnessTensor{T}) where {T} = StiffnessMatrix(
    c[1, 1, 1, 1],
    c[1, 1, 2, 2],
    c[1, 1, 3, 3],
    c[1, 1, 2, 3],
    c[1, 1, 1, 3],
    c[1, 1, 1, 2],
    c[2, 2, 2, 2],
    c[2, 2, 3, 3],
    c[2, 2, 2, 3],
    c[2, 2, 1, 3],
    c[2, 2, 1, 2],
    c[3, 3, 3, 3],
    c[3, 3, 2, 3],
    c[3, 3, 1, 3],
    c[3, 3, 1, 2],
    c[2, 3, 2, 3],
    c[2, 3, 1, 3],
    c[2, 3, 1, 2],
    c[1, 3, 1, 3],
    c[1, 3, 1, 2],
    c[1, 2, 1, 2],
)
Base.convert(::Type{ComplianceMatrix{T}}, s::ComplianceTensor{T}) where {T} =
    ComplianceMatrix(
        s[1, 1, 1, 1],
        s[1, 1, 2, 2],
        s[1, 1, 3, 3],
        2s[1, 1, 2, 3],
        2s[1, 1, 1, 3],
        2s[1, 1, 1, 2],
        s[2, 2, 2, 2],
        s[2, 2, 3, 3],
        2s[2, 2, 2, 3],
        2s[2, 2, 1, 3],
        2s[2, 2, 1, 2],
        s[3, 3, 3, 3],
        2s[3, 3, 2, 3],
        2s[3, 3, 1, 3],
        2s[3, 3, 1, 2],
        4s[2, 3, 2, 3],
        4s[2, 3, 1, 3],
        4s[2, 3, 1, 2],
        4s[1, 3, 1, 3],
        4s[1, 3, 1, 2],
        4s[1, 2, 1, 2],
    )
Base.convert(::Type{StiffnessTensor{T}}, c::StiffnessMatrix{T}) where {T} =
    StiffnessTensor(fromvoigt(SymmetricFourthOrderTensor{3,T}, c.data))
function Base.convert(::Type{ComplianceTensor{T}}, s::ComplianceMatrix{T}) where {T}
    ComplianceTensor(
        SymmetricFourthOrderTensor{3,T}(function (i, j, k, l)
            if i == j && k == l
                return s[i, k]
            elseif i != j && k != l  # 4 = 9 - (2+3), 5 = 9 - (1+3), 6 = 9 - (1+2)
                return s[9-(i+j), 9-(k+l)] / 4
            elseif i == j && k != l
                return s[i, 9-(k+l)] / 2
            else  # i != j && k == l
                return s[9-(i+j), k] / 2
            end
        end),
    )
end

Base.:*(c::StiffnessMatrix, ϵ::EngineeringStrain) = EngineeringStress(c.data ⋅ ϵ.data)
Base.:*(s::ComplianceMatrix, σ::EngineeringStress) = EngineeringStrain(s.data ⋅ σ.data)

contraction(c::StiffnessTensor, ε::TensorStrain, ::Val{2}) = TensorStress(c.data ⊡ ε.data)
contraction(s::ComplianceTensor, σ::TensorStress, ::Val{2}) = TensorStrain(s.data ⊡ σ.data)

@inline double_contraction(c::StiffnessTensor, ε::TensorStrain) = contraction(c, ε, Val(2))
@inline double_contraction(s::ComplianceTensor, σ::TensorStress) = contraction(s, σ, Val(2))


text:
###########################
# FunctionCell Interface  #
###########################
"""
	scaler!([f,] opts)

Returns a cell that scales data piped to it according to the scaling options specified
in `opts`. If the dataset is labeled, an additional function `f` can be specified to
obtain the labels calling `LearnBase.targets(f,data)`. The argument `opts` can be a
string specifying the scaling method or a `Dict(idx=>method)` if different methods
are to be used for different variables. In this case,

* `idx` can be an `Int`, `Vector{Int}` or `UnitRange` and specifies the variable indices
* `method` is a string that specifies the method. Available methods: `"mean"`, `"variance"`,
`"domain"`,`"c-mean"`, `"c-variance"` and `"2-sigma"`. Unrecognized options will be ignored.

# Examples
```
julia> a=[1.0 0 0; 0 1.0 1.0];

julia> w = a |> scaler!("mean")
Data scaler! (mean), 2 -> 2, trained

julia> a |> w
2×3 Array{Float64,2}:
 0.666667  -0.333333  -0.333333
 -0.666667   0.333333   0.333333

julia> a=datacell([1.0 -1 0 0 0; 5 0 1.0 1.0 1.0; 1 2 3 4 5], [0, 0, 1, 1, 1]);

julia> +a
3×5 Array{Float64,2}:
 1.0  -1.0  0.0  0.0  0.0
 5.0   0.0  1.0  1.0  1.0
 1.0   2.0  3.0  4.0  5.0

julia> w=a |> scaler!(Dict(1=>"mean", 2=>"2-sigma"))
Data scaler! (mixed), 3 -> 3, trained

julia> a|>w; +a
3×5 Array{Float64,2}:
 1.0       -1.0       0.0       0.0       0.0
 0.880132   0.321115  0.432918  0.432918  0.432918
 1.0        2.0       3.0       4.0       5.0
```
"""
scaler!(f::Function, opts::T where T<:AbstractString) = FunctionCell(scaler!, (f,opts), ModelProperties(), "Data scaler! ("*opts*")")
scaler!(opts::T where T<:AbstractString) = scaler!(identity, opts)

scaler!(f::Function, opts::T where T<:Dict) = FunctionCell(scaler!, (f, opts), ModelProperties(), "Data scaler! (mixed)")
scaler!(opts::T where T<:Dict) = scaler!(identity, opts)

############################
# DataCell/Array Interface #
############################
"""
	scaler!(data, [f,] opts)

Scales `data` according to the scaling options specified in `opts`.
If the dataset is labeled, an additional function `f` can be specified to
obtain the labels calling `LearnBase.targets(data,f)`. """
# Training
scaler!(x::T where T<:CellData, opts) = scaler!(strip(x), identity, opts)
scaler!(x::T where T<:CellData, f::Function, opts) = scaler!(strip(x), f, opts)

scaler!(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, opts) = scaler!(x[1], x[2], identity, opts)
scaler!(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, f::Function, opts) = scaler!(x[1], x[2], f, opts)

scaler!(x::T where T<:AbstractArray, opts) = scaler!(x, Void[], identity, opts)
scaler!(x::T where T<:AbstractArray, f::Function, opts) = scaler!(x, Void[], f, opts)

scaler!(x::T where T<:AbstractArray, y::S where S<: AbstractArray, f::Function, opts) = begin

	# Get dictionary or construct a proper one from original;
	# Inputs: scaling option and total number of variables
	_sdict_(x::T, m) where T<:String = Dict(i=>x for i in 1:m)
	_sdict_(x::T, m) where T<:Dict = begin
		out = Dict{Int, String}()
		for (k,v) in x
			if(k isa Int)
				@assert (k>=1)&&(k<= m) "[scaler!] Index $(k) for a variable out of bounds."
				push!(out, k=>v)
			elseif (k isa UnitRange{Int} || k isa Vector{Int})
				@assert (minimum(k)>=1)&&(maximum(k)<= m) "[scaler!] Index $(k) for a variable out of bounds."
				push!(out, (ki=>v for ki in k)...)
			else
				error("[scaler!] Unsupported variable index syntax.")
			end
		end
		return out
	end
	dopts = _sdict_(opts, nvars(x))

	# Determine scaler name (for information only)
	_scalername_(x) = begin
		if isempty(x) return "Data scaler! ()"
		elseif (length(unique(x)) == 1) return "Data scaler! ("*collect(x)[1]*")"
		else return "Data scaler! (mixed)"
		end
	end
	scalername = _scalername_(values(dopts))

	# Get targets (nothing for unlabeled data cells and arrays)
	_targets_(f,y) = targets(f,y)
	_targets_(f,::AbstractArray{Void}) = nothing
	labels = _targets_(f,y)
	@assert !any(isna.(labels))
	@assert !any(isnan.(labels))
	@assert (labels isa Vector)||(labels isa Void) "[scaler!] Labels need to be a Vector or ::Void."

	# Get priors
	priors = countmapn(labels)

	# Create model properties
	modelprops = ModelProperties(nvars(x), nvars(x), nothing,
			Dict("priors" => priors,						# Priors
			"variables_to_process" => collect(keys(dopts))				# Which columns to process (can be any sort of vector)
			)
	)

	# Define processing functions (x - data vector, t - targets array, p - class priors), returns Tuple(mean, variance, clip)
	# - Scalers not based on class information have also an overloaded version for labeled datasets which reverts to one where
	#   no labels are present
	# - Class based scalers expect targets to be a vector (not matrix) or default to non-class versions if possible
	_mean_(x::T, ::Void) where T<:AbstractVector = (1.0, -mean(x), false)
	_mean_(x::T, t::S) where T<:AbstractVector where S<:AbstractArray = _mean_(x,nothing)

	_variance_(x::T, ::Void) where T<:AbstractVector = (1/std(x), -mean(x)/std(x), false)
	_variance_(x::T, t::S) where T<:AbstractVector where S<:AbstractArray = _variance_(x,nothing)

	_domain_(x::T, ::Void) where T<:AbstractVector = begin
		maxv = maximum(x)
		minv = minimum(x)
		return (1/(maxv-minv+eps()), -minv/(maxv-minv+eps()), false)
	end
	_domain_(x::T, t::S) where T<:AbstractVector where S<:AbstractArray = _domain_(x,nothing)

	_cmean_(x::T, ::Void, p) where T<:AbstractVector = _mean_(x,nothing)
	_cmean_(x::T, t::S, p) where T<:AbstractVector where S<:AbstractVector = begin
		offset = 0.0
		@inbounds for c in keys(p)
			d = x[isequal.(t,c)]
			offset += p[c] * mean(d)
		end
		return (1.0, -offset, false)
	end

	_cvariance_(x::T, ::Void, p) where T<:AbstractVector = _variance_(x,nothing)
	_cvariance_(x::T, t::S, p) where T<:AbstractVector where S<:AbstractArray = begin
		offset = 0.0
		scale = 0.0
		@inbounds for c in keys(p)
			d = x[isequal.(t,c)]
			offset += p[c] * mean(d)
			scale += p[c] * var(d)
		end
		return (1.0/scale, -offset/scale, false)
	end

	_2sigma_(x::T, ::Void, p) where T<:AbstractVector = error("[scaler!] 2-sigma scaling needs labels")
	_2sigma_(x::T, t::S, p) where T<:AbstractVector where S<:AbstractVector = begin
		offset = 0.0
		scale = 0.0
		@inbounds for c in keys(p)
			d = x[isequal.(t,c)]
			offset += p[c] * mean(d)
			scale += p[c] * var(d)
		end
		scale = 4*sqrt(scale)
		offset = offset - 0.5*scale
		return (1.0/scale, -offset/scale, true)
	end

	# Create dictionary for each processed variable Dict(column=>(scale, offset, labels))
	modeldata = Dict{Int, Tuple{eltype(1.0),eltype(1.0), Bool}}()
	for (idx, v) in dopts
		variable = _variable_(x,idx)
		if v=="mean"
			s = _mean_(variable, labels)
		elseif v=="variance"
			s = _variance_(variable, labels)
		elseif v=="domain"
			s = _domain_(variable, labels)
		elseif v=="c-mean"
			s = _cmean_(variable, labels, priors)
		elseif v=="c-variance"
			s = _cvariance_(variable, labels, priors)
		elseif v=="2-sigma"
			s = _2sigma_(variable, labels, priors)
		else # do not scale if option not recognized
			warn("[scaler!] Option $(v) not recognized, will not process variable $(idx)")
			s = (1.0, 0.0, false)
		end
		push!(modeldata, idx => s )
	end

	# Return trained cell
	return FunctionCell(scaler!, Model(modeldata, modelprops), scalername)
end

# Execution
scaler!(x::T where T<:CellData, model::Model) = datacell(scaler!(strip(x), model))

scaler!(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractArray, model::Model) = (scaler!(x[1], model),x[2])

scaler!(x::T where T<:AbstractArray, model::Model) = begin

	@inbounds for (idx, (scale, offset, clip)) in model.data
		variable = _variable_(x,idx) 	# Assign temp variable

		# Process variable vector
		variable .*= scale
		variable .+= offset
		ifelse(clip, clamp!(variable,0.0,1.0), nothing)
	end
	return x
end


text:
using Yota
using Distributions
using GradDescent
using MLDataUtils
using MLDatasets
using StatsBase
using ImageView

include("common.jl")

# variational autoencoder with Gaussian observed and latent variables
mutable struct VAE{T}
    # encoder / recognizer
    We1::AbstractMatrix{T}  # encoder: layer 1 weight
    be1::AbstractVector{T}  # encoder: layer 1 bias
    We2::AbstractMatrix{T}  # encoder: layer 2 weight
    be2::AbstractVector{T}  # encoder: layer 2 bias
    We3::AbstractMatrix{T}  # encoder: layer 3 mu weight
    be3::AbstractVector{T}  # encoder: layer 3 mu bias
    We4::AbstractMatrix{T}  # encoder: layer 3 log(sigma^2) weight
    be4::AbstractVector{T}  # encoder: layer 3 log(sigma^2) bias
    # decoder / generator
    Wd1::AbstractMatrix{T}  # decoder: layer 1 weight
    bd1::AbstractVector{T}  # decoder: layer 1 bias
    Wd2::AbstractMatrix{T}  # decoder: layer 2 weight
    bd2::AbstractVector{T}  # decoder: layer 2 bias
    Wd3::AbstractMatrix{T}  # decoder: layer 3 weight
    bd3::AbstractVector{T}  # decoder: layer 3 bias
end

function Base.show(io::IO, m::VAE{T}) where T
    print(io, "VAE{$T}($(size(m.We1,2)), $(size(m.We1,1)), $(size(m.We2,1)), " *
          "$(size(m.We3,1)), $(size(m.Wd1,1)), $(size(m.Wd2,1)), $(size(m.Wd3,1)))")
end

VAE{T}(n_inp, n_he1, n_he2, n_z, n_hd1, n_hd2, n_out) where T =
    VAE{T}(
        # encoder
        xavier_init(n_he1, n_inp),
        zeros(n_he1),
        xavier_init(n_he2, n_he1),
        zeros(n_he2),
        xavier_init(n_z, n_he2),
        zeros(n_z),
        xavier_init(n_z, n_he2),
        zeros(n_z),
        # decoder
        xavier_init(n_hd1, n_z),
        zeros(n_hd1),
        xavier_init(n_hd2, n_hd1),
        zeros(n_hd2),
        xavier_init(n_out, n_hd2),
        zeros(n_out)
    )

function encode(m::VAE, x)
    he1 = softplus.(m.We1 * x .+ m.be1)
    he2 = softplus.(m.We2 * he1 .+ m.be2)
    mu = m.We3 * he2 .+ m.be3
    log_sigma2 = m.We4 * he2 .+ m.be4
    return mu, log_sigma2
end

function decode(m::VAE, z)
    hd1 = softplus.(m.Wd1 * z .+ m.bd1)
    hd2 = softplus.(m.Wd2 * hd1 .+ m.bd2)
    x_rec = logistic.(m.Wd3 * hd2 .+ m.bd3)
    return x_rec
end

function vae_cost(m::VAE, eps, x)
    mu, log_sigma2 = encode(m, x)
    z = mu .+ sqrt.(exp.(log_sigma2)) .* eps
    x_rec = decode(m, z)
    # loss
    rec_loss = -sum(x .* log.(1e-10 .+ x_rec) .+ (1 .- x) .* log.(1e-10 + 1.0 .- x_rec); dims=1)
    KLD = -0.5 .* sum(1 .+ log_sigma2 .- mu .^ 2.0f0 - exp.(log_sigma2); dims=1)
    cost = mean(rec_loss .+ KLD)
end

function StatsBase.fit!(m::VAE{T}, X::AbstractMatrix{T};
              n_epochs=50, batch_size=100) where T
    for epoch in 1:n_epochs
        empty!(Yota.GRAD_CACHE)
        print("Epoch $epoch: ")
        epoch_cost = 0
        t = @elapsed for (i, x) in enumerate(eachbatch(X, size=batch_size))
            eps = typeof(x)(rand(Normal(0, 1), size(m.We3, 1), batch_size))
            cost, g = Yota._grad(vae_cost, m, eps, x)
            update!(m, g[1], (x, gx) -> x .- 0.01gx)
            epoch_cost += cost
            gW = g[1][(:We1,)]
            @info "epoch = $epoch; loss = $cost; norm(∇We1) = $(sum(abs2, gW)); norm(We1) = $(sum(abs2, m.We1))"
        end
        println("avg_cost=$(epoch_cost / (size(X,2) / batch_size)), elapsed=$t")
    end
    return m
end

function reconstruct(m::VAE, x::AbstractVector)
    x = reshape(x, length(x), 1)
    mu, _ = encode(m, x)
    z = mu
    x_rec = decode(m, z)
    return x_rec
end

function show_pic(x)
    reshape(x, 28, 28)' |> imshow
end

function show_recon(m, x)
    x_ = reconstruct(m, x)
    show_pic(x)
    show_pic(x_)
end

function main()
    m = VAE{Float64}(784, 500, 500, 20, 500, 500, 784) #  |> to_cuda

    X, _ = MNIST.traindata()
    X = convert(Matrix{Float64}, reshape(X, 784, 60000))
    @time m = fit!(m, X)

    # check reconstructed image
    for i=1:2:10
        show_recon(m, X[:, i])
    end
end


text:
# Importing functions from, Adding functions to, and renaming functions in LinearAlgebra

module LinearAlgebraUtilities
    import LinearAlgebra.norm,
           LinearAlgebra.checksquare

    export norm,
           isSquare,
           isGeneric,
           isSymplectic,
           Omega, Ω,
           dsum, ⊕

    Omega(n::Int)::AbstractMatrix = cat(fill([0 1; -1 0], n)...; dims=(1,2))
    Omega(matrix::AbstractMatrix)::AbstractMatrix = Omega(size(matrix)[1] ÷ 2)
    Ω = Omega

    function isSquare(matrix::AbstractMatrix)::Bool
        try
            n = checksquare(matrix)
            return true
        catch err
            return false
        end
    end
    isGeneric(matrix::AbstractMatrix, tolerance::Real = 1e-6)::Bool = all(x -> abs(x) > tolerance, matrix )
    function isSymplectic(matrix::AbstractMatrix, tolerance = 1e-6::Real)::Bool
        if !(isSquare(matrix)) return false end
        if (checksquare(matrix) % 2 != 0) return false end
        if (norm(matrix*Omega(checksquare(matrix) ÷ 2)*matrix' - Omega(checksquare(matrix) ÷ 2)) >= tolerance) return false end
        return true
    end

    dsum(M1::AbstractMatrix, M2::AbstractMatrix)::AbstractMatrix = cat(M1, M2; dims=(1,2))
    dsum(M1::AbstractMatrix, M2::AbstractMatrix, Ms::Vararg{AbstractMatrix})::AbstractMatrix = dsum(dsum(M1, M2),  mapreduce(identity, dsum, Ms ))
    ⊕ = dsum
end


text:
# adapted from http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_hw5a.py

import GDALUtils; const GU = GDALUtils

GU.registerdrivers() do
    GU.read("ospy/data5/aster.img") do inDS
        # get image size
        rows = GU.height(inDS)
        cols = GU.width(inDS)
        bands = GU.nraster(inDS)

        # get the band and block sizes
        inband2 = GU.fetchband(inDS, 2)
        inband3 = GU.fetchband(inDS, 3)
        (xbsize, ybsize) = GU.getblocksize(inband2)

        buffer2 = Array(Float32, ybsize, xbsize)
        buffer3 = Array(Float32, ybsize, xbsize)
        ndvi    = Array(Float32, ybsize, xbsize)
        # create the output image
        GU.create("tmp/ndvi.img", GU.getdriver(inDS),
                  width=cols, height=rows, nbands=1, dtype=Float32) do outDS
            for ((i,j),(nrows,ncols)) in GU.BlockIterator(inband2)
                GU.rasterio!(inband2, buffer2, ncols, nrows, j, i)
                GU.rasterio!(inband3, buffer3, ncols, nrows, j, i)
                data2 = buffer2[1:nrows,1:ncols]
                data3 = buffer3[1:nrows,1:ncols]
                for row in 1:nrows, col in 1:ncols
                    denominator = data2[row, col] + data3[row, col]
                    if denominator > 0
                        numerator = data3[row, col] - data2[row, col]
                        ndvi[row, col] = numerator / denominator
                    else
                        ndvi[row, col] = -99
                    end
                end
                # write the data
                GU.update!(outDS, ndvi, 1, ncols, nrows, j, i)
            end
            # flush data to disk, set the NoData value and calculate stats
            outband = GU.fetchband(outDS, 1)
            GU.flushcache(outband)
            GU.setnodatavalue(outband, -99)
            # stats = outBand.GetStatistics(0, 1)

            # georeference the image and set the projection
            GU.setgeotransform!(outDS, GU.getgeotransform(inDS))
            GU.setproj(outDS, GU.getproj(inDS))

            # build pyramids
            # gdal.SetConfigOption('HFA_USE_RRD', 'YES')
            GU.buildoverviews(outDS,
                              Cint[2,4,8,16,32,64,128], # overview list
                                                        # bandlist (omit to include all bands)
                              resampling="NEAREST")     # resampling method (default: nearest)
        end
    end
end


text:
include("framework.jl")
using CustomTest
using Base.Test

my_tests = ["test_statespace",
            "test_transferfunction",
            "test_generalizedtf",
            "test_zpk",
            "test_analysis",
            "test_connections",
            "test_discrete",
            "test_linalg",
            "test_simplification",
            "test_freqresp",
            "test_synthesis",
            "test_matrix_comps",
            "test_timeresp",
            "test_conversion"]

try
    Pkg.installed("ControlExamplePlots")
    push!(my_tests, "test_plots")
catch
    warn("The unregistered package ControlExamplePlots is currently needed to test plots, install using:
            Pkg.clone(\"https://github.com/JuliaControl/ControlExamplePlots.jl.git\")")
end

run_tests(my_tests)


text:
"""
    post_plot_rotor_sweep(aircraft, parameters, step_range, step_symbol)

# Arguments:

* `aircraft::Aircraft`: aircraft system struct
* `parameters <: Parameters`: inherits from the `Parameters` type; object containing data required by `<: Action` functions
* `step_range::AbstractArray`: range of times defining simulation
* `step_symbol::String`: defines the step, e.g. `alpha` or `time`

# Modifies:

* `aircraft::Aircraft`
* `parameters <: Parameters`

# Returns:

* `flag::Bool`: true if an action experiences errors

`parameters <: Parameters` requires:

* `Js::AbstractArray`: array of advance ratios of the ith rotor at the jth step
* `CTs::Array{Float64,2}`: array of thrust coefficients of the ith rotor at the jth step
* `CQs::Array{Float64,2}`: array of torque coefficients of the ith rotor at the jth step
* `ηs::Array{Float64,2}`: array of propulsive efficiencies of the ith rotor at the jth step
* `rotor_names::Vector{String}`: vector of rotor names for use in plot legend
* `plot_directory::String`: directory where plots are saved
* `plot_base_name::String`: first portion of the saved figure file name
* `plot_extension::String`: extension of saved figure files

"""
function post_plot_rotor_sweep(aircraft, parameters, step_range, step_symbol)

    # extract info
    CTs = parameters.CTs
    CQs = parameters.CQs
    ηs = parameters.ηs
    data = [CTs, CQs, ηs]
    rotor_names = parameters.rotor_names
    plot_directory = parameters.plot_directory
    plot_base_name = parameters.plot_base_name
    plot_extension = parameters.plot_extension
    nrotors = length(aircraft.rotor_system.index)

    # check data
    @assert length(rotor_names) == nrotors "length of rotor_names is inconsistent: got $(length(rotor_names)); expected $nrotors"
    @assert isdir(plot_directory) "plot_directory does not exist"

    # prepare figure
    fig = plt.figure(plot_base_name * "_rotor_sweep")
    fig.clear()
    fig.add_subplot(311, ylabel = L"C_T")
    fig.add_subplot(312, ylabel = L"C_Q")
    fig.add_subplot(313, ylabel = L"\eta", xlabel = L"J")
    axs = fig.get_axes()

    # plot
    for (i,ax) in enumerate(axs)
        for jrotor in 1:nrotors
            cratio = jrotor / nrotors
            Js = parameters.Js[jrotor,:]
            ax.plot(Js, data[i][jrotor,:], color = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label = rotor_names[jrotor])
        end
    end

    axs[1].legend(loc="upper left", bbox_to_anchor=(1.01,1))
    fig.set_size_inches(7, 4, forward=true)
    fig.tight_layout()

    # save
    savepath = joinpath(plot_directory, plot_base_name * "_rotor_sweep" * plot_extension)
    fig.savefig(savepath, bbox_inches = "tight")

    # also plot power as a function of J
    fig_power = plt.figure(plot_base_name * "_rotor_power")
    fig_power.clear()
    fig_power.add_subplot(121, ylabel = L"P", xlabel = L"J")
    fig_power.add_subplot(122, xlabel = L"n")
    axs_power = fig_power.get_axes()

    # plot
    for jrotor in 1:nrotors
        cratio = jrotor / nrotors
        Js = parameters.Js[jrotor,:]
        ns = parameters.omegas[jrotor,:] ./ (2*pi)
        Ps = parameters.Ps[jrotor,:]
        axs_power[1].plot(Js, Ps, color = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label = rotor_names[jrotor])
        axs_power[2].plot(ns, Ps, color = (0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label = rotor_names[jrotor])
    end

    axs_power[2].legend(loc="upper left", bbox_to_anchor=(1.01,1))
    fig_power.set_size_inches(9, 3, forward=true)
    fig_power.tight_layout()

    # save
    savepath = joinpath(plot_directory, plot_base_name * "_rotor_power" * plot_extension)
    fig.savefig(savepath, bbox_inches = "tight")

    return false
end

"""
    post_plot_rotor_sweep(system, step_range) <: PostAction

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft`: system to be simulated
* `step_range::AbstractArray`: defines each step of the simulation

# Returns:

* `Js::AbstractArray`: array of advance ratios of the ith rotor at the jth step
* `CTs::Array{Float64,2}`: array of thrust coefficients of the ith rotor at the jth step
* `CQs::Array{Float64,2}`: array of torque coefficients of the ith rotor at the jth step
* `ηs::Array{Float64,2}`: array of propulsive efficiencies of the ith rotor at the jth step
* `rotor_names::Vector{String}`: vector of rotor names for use in plot legend
* `plot_directory::String`: directory where plots are saved
* `plot_base_name::String`: first portion of the saved figure file name
* `plot_extension::String`: extension of saved figure files

"""
function post_plot_rotor_sweep(aircraft, step_range)

    _, Js, _, _, CTs, CQs, ηs, _, _ = solve_rotor_nondimensional(aircraft, step_range)
    nrotors = length(aircraft.rotor_system.index)
    rotor_names = ["rotor$i" for i in 1:nrotors]
    plot_directory = ""
    plot_base_name = "default"
    plot_extension = ".pdf"

    return Js, CTs, CQs, ηs, rotor_names, plot_directory, plot_base_name, plot_extension
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
module IntegrationTests

using STLCutters
using Gridap
import Gridap: ∇
using GridapEmbedded
using Test

using STLCutters: compute_stl_model
using STLCutters: read_stl, merge_nodes, get_bounding_box

# Manufactured solution
u0(x) = x[1] + x[2] - x[3]

X,T,N = read_stl(joinpath(@__DIR__,"data/cube.stl"))
stl = compute_stl_model(T,X)
stl = merge_nodes(stl)
# writevtk(stl,"geo")
n = 10
δ = 0.2
pmin,pmax = get_bounding_box(stl)
diagonal = pmax-pmin
pmin = pmin - diagonal*δ
pmax = pmax + diagonal*δ
partition = (n,n,n)

geo = STLGeometry(stl)
bgmodel = CartesianDiscreteModel(pmin,pmax,partition)

# Cut the background model

cutgeo, = cut(bgmodel,geo)

# Setup integration meshes
Ω = Triangulation(cutgeo,PHYSICAL)
Γd = EmbeddedBoundary(cutgeo)
Γg = GhostSkeleton(cutgeo)

# Setup normal vectors
n_Γd = get_normal_vector(Γd)
n_Γg = get_normal_vector(Γg)

#writevtk(Ω,"trian_O")
#writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Γg,"trian_Gg",cellfields=["normal"=>n_Γg])
#writevtk(Triangulation(bgmodel),"bgtrian")

# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
dΓd = Measure(Γd,degree)
dΓg = Measure(Γg,degree)

vol = sum( ∫(1)*dΩ  )
surf = sum( ∫(1)*dΓd )

@test vol ≈ 1
@test surf ≈ 6

# Setup FESpace
Ω_act = Triangulation(cutgeo,ACTIVE,geo)
V = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
v = FEFunction(V,rand(num_free_dofs(V)))
u = interpolate(u0,V)

# Check divergence theorem
a = sum( ∫( ∇(v)⋅∇(u) ) * dΩ )
b = sum( ∫( v⋅n_Γd⋅∇(u) ) * dΓd )
@test abs( a-b ) < 1e-9

end


text:
using Test
using ReferenceFrameRotations
using LinearAlgebra
using Random
using StaticArrays

import Base: isapprox

################################################################################
#                             Auxiliary functions
################################################################################

# Function to uniformly sample an angle in [-π, π].
_rand_ang(T = Float64) = -T(π) + T(2π) * rand(T)

# Function to uniformly sample an angle in [-π / 2, +π / 2].
_rand_ang2(T = Float64) = -T(π) / 2 + T(π) * rand(T)

# Function to uniformly sample an angle in [0.1, 1.5].
_rand_ang3(T = Float64) = T(0.1) + T(1.4) * rand(T)

# Normalize angle between [-π, +π].
function _norm_ang(α::T) where T
    Tf = float(T)
    αr = mod(α, Tf(2π))  # Make sure α ∈ [0, 2π].
    if αr ≤ π
        return αr
    else
        return αr - Tf(2π)
    end
end

# Available rotations.
const valid_rot_seqs_2angles = [
    :XY,
    :XZ,
    :YX,
    :YZ,
    :ZX,
    :ZY,
]

const valid_rot_seqs = [
    :XYX,
    :XYZ,
    :XZX,
    :XZY,
    :YXY,
    :YXZ,
    :YZX,
    :YZY,
    :ZXY,
    :ZXZ,
    :ZYX,
    :ZYZ
]

# Define the function `isapprox` for `EulerAngleAxis` to make the comparison
# easier.
function isapprox(x::EulerAngleAxis, y::EulerAngleAxis; keys...)
    a = isapprox(x.a, y.a; keys...)
    v = isapprox(x.v, y.v; keys...)

    return a && v
end

# Define the function `isapprox` for `EulerAngles` to make the comparison
# easier.
function isapprox(x::EulerAngles, y::EulerAngles; keys...)
    a1 = isapprox(x.a1, y.a1; keys...)
    a2 = isapprox(x.a2, y.a2; keys...)
    a3 = isapprox(x.a3, y.a3; keys...)
    r  = x.rot_seq == y.rot_seq

    return a1 && a2 && a3 && r
end

@time @testset "Direction cosine matrices" verbose = true begin
    include("./dcm/misc.jl")
    include("./dcm/kinematics.jl")
    include("./dcm/orthonormalize.jl")
end
println("")

@time @testset "Euler angle and axis" verbose = true begin
    include("./angleaxis/functions.jl")
    include("./angleaxis/operations.jl")
end
println("")

@time @testset "Euler angles" verbose = true begin
    include("./angle/functions.jl")
    include("./angle/operations.jl")
end
println("")

@time @testset "Quaternions" verbose = true begin
    include("./quaternion/constructors.jl")
    include("./quaternion/functions.jl")
    include("./quaternion/julia_api.jl")
    include("./quaternion/kinematics.jl")
    include("./quaternion/operations.jl")
end
println("")

@time @testset "Conversions" verbose = true begin
    include("./conversions/angle_to_angle.jl")
    include("./conversions/angle_to_angleaxis.jl")
    include("./conversions/angle_to_dcm.jl")
    include("./conversions/angle_to_quat.jl")
    include("./conversions/angle_to_rot.jl")
    include("./conversions/angleaxis_to_angle.jl")
    include("./conversions/angleaxis_to_dcm.jl")
    include("./conversions/angleaxis_to_quat.jl")
    include("./conversions/api.jl")
    include("./conversions/dcm_to_angleaxis.jl")
    include("./conversions/dcm_to_angle.jl")
    include("./conversions/dcm_to_quaternion.jl")
    include("./conversions/quaternion_to_angle.jl")
    include("./conversions/quaternion_to_angleaxis.jl")
    include("./conversions/quaternion_to_dcm.jl")
end
println("")

@time @testset "Compose rotations" verbose = true begin
    include("./compose_rotations.jl")
end
println("")

@time @testset "Invert rotations" verbose = true begin
    include("./inv_rotations.jl")
end
println("")

@time @testset "Random rotations" verbose = true begin
    include("./random.jl")
end
println("")


text:
using IndexedTables
using FileIO
using SurrealNumbers
out_dir = "Data/"

# subtraction table
x = convert.(SurrealFinite, [0, 1/2, 1, 2])
A = zeros(SurrealFinite, length(x), length(x))
B = zeros(Int, length(x), length(x))
for i=1:length(x)
    for j=1:length(x)
        A[i,j] = x[i] - x[j]
        B[i,j] = generation(x[i]) + generation(x[j])
     end
end
float.(A)
A
generation.(A)
# A .== canonicalise.(A)
iscanonical.(A)
B
C = dag_stats.(A)
n = [C[i,j].nodes for i=1:size(C)[1], j=1:size(C)[2] ]

# subtraction table 2
ns1 = [4; 4; 4; 4; 4; 4; 1; 2; 3]
ns2 = [0; 1; 2; 3; 4; 5; 1; 2; 3]
n = length(ns1)
m = [10; 10; 10; 10; 10; 10; 10; 10; 10]
x = Array{SurrealFinite,1}(n)
t = Array{Float64,1}(n)
bytes = Array{Float64,1}(n)
gctime = Array{Float64,1}(n)
c_plus = Array{Integer,1}(n)
c_times = Array{Integer,1}(n)
c_created = Array{Integer,1}(n)
c_equals = Array{Integer,1}(n)
c_leq = Array{Integer,1}(n)
i = 0
for i=1:n
    println(ns1[i], " - ", ns2[i])
    t[i] = 0
    b = 0
    g = 0
    mem = 0
    val = 0
    for j=1:m[i]
        clearcache()
        val, time, b, g, mem = @timed convert(SurrealFinite, ns1[i]) - convert(SurrealFinite, ns2[i])
        t[i] += time
    end
    x[i] = val
    bytes[i] = b
    gctime[i] = g

    c_plus[i] = Count['+']
    c_times[i] = Count['*']
    c_created[i] = Count['c']
    c_equals[i] = Count['=']
    c_leq[i] = Count['≦']
end
s = dag_stats.( x )
no = [s[i].nodes for i=1:length(s)]
ed = [s[i].edges for i=1:length(s)]
p = [s[i].paths for i=1:length(s)]
d = (e+1) ./ no
g = [s[i].generation for i=1:length(s)]
t = t./m[1:n]
bytes = bytes

t1 = table(@NT(i=0:n-1,
               n1=ns1,
               n2=ns2,
               nodes=no,
               edges=ed,
               paths=p,
               density=d,
               generation=g,
               time=t,
               bytes=bytes,
               adds=c_plus,
               prods=c_times,
               created=c_created,
               equals=c_equals,
               leq=c_leq);
           pkey = [:i])

save("$(out_dir)subtraction_tables.csv", t1)

# examples with diagrams
a21 = convert(SurrealFinite, 1) - convert(SurrealFinite, 1)
file = "$(out_dir)subtraction_ex_1.dot"
FID = open(file, "w")
surreal2dag(FID, a21)
close(FID)
run(`dot -Tpdf -O $file`)

a22 = convert(SurrealFinite, 2) - convert(SurrealFinite, 2)
file = "$(out_dir)subtraction_ex_2.dot"
FID = open(file, "w")
surreal2dag(FID, a22)
close(FID)
run(`dot -Tpdf -O $file`)

a22 = convert(SurrealFinite, 2) - convert(SurrealFinite, 2)
file = "$(out_dir)subtraction_ex_2.dot"
FID = open(file, "w")
surreal2dag(FID, a22)
close(FID)
run(`dot -Tpdf -O $file`)

a23 = convert(SurrealFinite, 1/2) - convert(SurrealFinite, 1/2)
file = "$(out_dir)subtraction_ex_3.dot"
FID = open(file, "w")
surreal2dag(FID, a23)
close(FID)
run(`dot -Tpdf -O $file`)

a41 = convert(SurrealFinite, 4) - convert(SurrealFinite, 1)
file = "$(out_dir)subtraction_ex_41.dot"
FID = open(file, "w")
surreal2dag(FID, a41)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

a42 = convert(SurrealFinite, 4) - convert(SurrealFinite, 2)
file = "$(out_dir)subtraction_ex_42.dot"
FID = open(file, "w")
surreal2dag(FID, a42)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

a43 = convert(SurrealFinite, 4) - convert(SurrealFinite, 3)
file = "$(out_dir)subtraction_ex_43.dot"
FID = open(file, "w")
surreal2dag(FID, a43)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)

a5 = convert(SurrealFinite, 3) - (convert(SurrealFinite, 1)-convert(SurrealFinite, 1))
file = "$(out_dir)subtraction_ex_5.dot"
FID = open(file, "w")
surreal2dag(FID, a5)
close(FID)
run(`dot -Tpdf -O $file`)
run(`dot -Tsvg -O $file`)


text:
#
# Structures required for the computation of the 3D density of the gmd function
#

struct GMD3DGrid
  x :: Vector{Float64}
  atom :: Int64
  dmin :: Float64
  rho :: Float64
end

mutable struct MutableGMD3DGrid
  x :: Vector{Float64}
  atom :: Int64
  dmin :: Float64
  rho :: Float64
end

MutableGMD3DGrid() = MutableGMD3DGrid(zeros(3),0,0.,0.)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

const MortarSegmentation = Dict{Int, Vector{Tuple{Int, Vector{Float64}}}}

"""
    calculate_segments(slave_element_ids::Vector{Int},
                       master_element_ids::Vector{Int},
                       elements::Dict{Int, Vector{Int}},
                       element_types::Dict{Int, Symbol},
                       coords::Dict{Int, Vector{Float64}},
                       normals::Dict{Int, Vector{Float64}})

Given slave surface elements, master surface elements, nodal coordinates and
normal direction on nodes of slave surface elements, calculate contact
segments.

Return type is a dictionary, where key is slave element id and values is a
list of master elements giving contribution to that slave elements and
xi-coordinates of slave side element.

# Example
```jldoctest
elements = Dict(1 => [1, 2], 2 => [3, 4])
element_types = Dict(1 => :Seg2, 2 => :Seg2)
coords = Dict(
    1 => [1.0, 2.0],
    2 => [3.0, 2.0],
    3 => [0.0, 2.0],
    4 => [2.0, 2.0])
normals = Dict(
    1 => [0.0, -1.0],
    2 => [0.0, -1.0])
slave_ids = [1]
master_ids = [2]
segments = calculate_segments(slave_ids, master_ids, elements,
                              element_types, coords, normals)

# output

Dict{Int64,Array{Tuple{Int64,Array{Float64,1}},1}} with 1 entry:
  1 => Tuple{Int64,Array{Float64,1}}[(2, [-1.0, -0.0])]

```

Here, output result means that slave element #1 has segment with master
element(s) #2 with dimensionless slave element coordinate xi = [-1, 0].
That is, the start and end point of projection in physical coordinate
system is:

    x_start = 1/2*(1-xi[1])*xs1 + 1/2*(1+xi[1])*xs2
    x_stop = 1/2*(1-xi[2])*xs1 + 1/2*(1+xi[2])*xs2

"""
function calculate_segments(slave_element_ids, master_element_ids, elements,
                            element_types, coords, normals)
    S = MortarSegmentation()
    for sid in slave_element_ids
        @assert element_types[sid] == :Seg2
        haskey(S, sid) || (S[sid] = [])
        scon = elements[sid]
        xs1 = coords[scon[1]]
        xs2 = coords[scon[2]]
        ns1 = normals[scon[1]]
        ns2 = normals[scon[2]]
        for mid in master_element_ids
            @assert element_types[mid] == :Seg2
            mcon = elements[mid]
            xm1 = coords[mcon[1]]
            xm2 = coords[mcon[2]]
            # first project from slave to master, to find out
            # are we completely outside of domain
            xi2a = project_from_slave_to_master(Val{:Seg2}, xs1, ns1, xm1, xm2)
            xi2b = project_from_slave_to_master(Val{:Seg2}, xs2, ns2, xm1, xm2)
            xi2a >  1.0 && xi2b >  1.0 && continue
            xi2a < -1.0 && xi2b < -1.0 && continue
            xi1a = project_from_master_to_slave(Val{:Seg2}, xm1, xs1, xs2, ns1, ns2)
            xi1b = project_from_master_to_slave(Val{:Seg2}, xm2, xs1, xs2, ns1, ns2)
            xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
            isapprox(abs(xi1[2]-xi1[1]), 0.0) && continue
            push!(S[sid], (mid, xi1))
        end
    end
    return S
end


text:
"""
    function get_Pg_thermal_bindingPg(args... ;kwargs...)

Identify generators which Active Power is binding (Skips Slack) and returns a matrix with
their Name, Power output and Bus number (sorted by Bus).

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

"""
function get_Pg_thermal_bindingPg(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_gen_tol::Float64
    )

    # Define elements
    gen_loc = 0
    ngen_bindPg = 0 #number of gen binding (start of counter)
    gens_thermal_bindPg_Pg = Array{Float64}(undef,0)
    gens_thermal_bindPg_names = Array{String}(undef,0)
    gens_thermal_bindPg_busnumbers = Array{Int64}(undef,0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised active power

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if  bindub || bindlb #if binding dual is !≈  0.0 then
                ngen_bindPg = ngen_bindPg +1
                gen_thermal_bindPg_name = gen_thermal.name
                resize!(gens_thermal_bindPg_Pg,ngen_bindPg)
                resize!(gens_thermal_bindPg_names,ngen_bindPg)
                resize!(gens_thermal_bindPg_busnumbers,ngen_bindPg)
                gens_thermal_bindPg_Pg[ngen_bindPg] =
                    all_PGenThermal[1, gen_thermal_bindPg_name]
                gens_thermal_bindPg_names[ngen_bindPg] = gen_thermal_bindPg_name
                gens_thermal_bindPg_busnumbers[ngen_bindPg] = gen_thermal.bus.number
            end
        end
    end

    # Matrix of generators with their information sorted
    gens_thermal_bindPg_Matrix = Array{Any}(undef,(ngen_bindPg, 3));
    gens_thermal_bindPg_Matrix[:, 1] = gens_thermal_bindPg_names;
    gens_thermal_bindPg_Matrix[:, 2] = gens_thermal_bindPg_Pg;
    gens_thermal_bindPg_Matrix[:, 3] = gens_thermal_bindPg_busnumbers;
    # Sort by bus to match with PTDF structure
    gens_thermal_bindPg_Matrix =
     gens_thermal_bindPg_Matrix[sortperm(gens_thermal_bindPg_Matrix[:, 3]), :]

    return gens_thermal_bindPg_Matrix
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# Function for loading sample data

"""
    get_data(dataset, poollength, datalength, horizon)

Function for importing one of the sample datasets in the repository. `dataset` can be one of
`:solar`, `:traffic`, `:exchange_rate` or `:electricity`. `poollength` gives the number of
timesteps to pool for the model, `datalength` determines the batch size of the output data,
and `horizon` determines the number of time steps that should be forecasted by the model.
"""
function get_data(dataset, poollength, datalength, horizon)
    admissible = [:solar, :traffic, :exchange_rate, :electricity]
    dataset in admissible || error("Sample data not implemented")

    (dataset == :solar) && (BSON.@load "../data/solar_AL.bson" inp_raw)
    (dataset == :traffic) && (BSON.@load "../data/traffic.bson" inp_raw)
    (dataset == :exchange_rate) && (BSON.@load "../data/exchange_rate.bson" inp_raw)
    (dataset == :electricity) && (BSON.@load "../data/electricity.bson" inp_raw)

    lostIdx = horizon + poollength
    datalength = min(datalength,size(inp_raw,1) - lostIdx)
    out_ft = similar(inp_raw, size(inp_raw,2), poollength, 1, datalength)
    for i=1:poollength
        for j=1:datalength
            out_ft[:,i,1,j] = inp_raw[poollength + j - i,:]
        end
    end
    out_lb = inp_raw[(1 + lostIdx :datalength + lostIdx ),1]
    return out_ft, out_lb
end

max_stars_repo_name: AndyTGardner/FluxArchitectures

max_issues_repo_name: AndyTGardner/FluxArchitectures

max_forks_repo_name: AndyTGardner/FluxArchitectures


text:
using Test
using Dash
using Dash.TableFormat
import JSON3
@testset "named values" begin
   @test TableFormat.Align.left.value == "<"
   a = TableFormat.Align.default
   @test a == TableFormat.Align.default
   @test TableFormat.possible_values(TableFormat.Align) == "Align.default, Align.left, Align.right, Align.center, Align.right_sign"
   @test a != TableFormat.Align.left
   @test a != TableFormat.Scheme.default
end

@testset "specifier" begin
    f = Format()

    @test f.specifier[:align] == TableFormat.Align.default.value
    TableFormat.align!(f, TableFormat.Align.left)
    @test f.specifier[:align] == TableFormat.Align.left.value
    @test_throws ArgumentError TableFormat.align!(f, TableFormat.Group.no)

    TableFormat.fill!(f, '-')
    @test f.specifier[:fill] == "-"
    TableFormat.fill!(f, "+")
    @test f.specifier[:fill] == "+"
    @test_throws ArgumentError TableFormat.fill!(f, "++")

    TableFormat.group!(f, TableFormat.Group.yes)
    @test f.specifier[:group] == TableFormat.Group.yes.value

    TableFormat.group!(f, true)
    @test f.specifier[:group] == TableFormat.Group.yes.value
    TableFormat.group!(f, false)
    @test f.specifier[:group] == TableFormat.Group.no.value
    @test_throws ArgumentError TableFormat.group!(f, "ffff")

    TableFormat.padding!(f, TableFormat.Padding.yes)
    @test f.specifier[:padding] == TableFormat.Padding.yes.value
    TableFormat.padding!(f, true)
    @test f.specifier[:padding] == TableFormat.Padding.yes.value
    TableFormat.padding!(f, false)
    @test f.specifier[:padding] == TableFormat.Padding.no.value
    @test_throws ArgumentError TableFormat.padding!(f, "ffff")

    TableFormat.padding_width!(f, 100)
    @test f.specifier[:width] == 100
    TableFormat.padding_width!(f, nothing)
    @test f.specifier[:width] == ""
    @test_throws ArgumentError TableFormat.padding_width!(f, "ffff")
    @test_throws ArgumentError TableFormat.padding_width!(f, -10)

    TableFormat.precision!(f, 3)
    @test f.specifier[:precision] == ".3"
    TableFormat.precision!(f, nothing)
    @test f.specifier[:precision] == ""
    @test_throws ArgumentError TableFormat.precision!(f, "ffff")
    @test_throws ArgumentError TableFormat.precision!(f, -10)

    @test f.specifier[:type] == TableFormat.Scheme.default.value
    TableFormat.scheme!(f, TableFormat.Scheme.decimal)
    @test f.specifier[:type] == TableFormat.Scheme.decimal.value
    @test_throws ArgumentError TableFormat.scheme!(f, "ccc")

    @test f.specifier[:sign] == TableFormat.Sign.default.value
    TableFormat.sign!(f, TableFormat.Sign.negative)
    @test f.specifier[:sign] == TableFormat.Sign.negative.value
    @test_throws ArgumentError TableFormat.sign!(f, "ccc")

    @test f.specifier[:symbol] == TableFormat.DSymbol.no.value
    TableFormat.symbol!(f, TableFormat.DSymbol.yes)
    @test f.specifier[:symbol] == TableFormat.DSymbol.yes.value
    @test_throws ArgumentError TableFormat.symbol!(f, "ccc")

    @test f.specifier[:trim] == TableFormat.Trim.no.value
    TableFormat.trim!(f, TableFormat.Trim.yes)
    @test f.specifier[:trim] == TableFormat.Trim.yes.value
    TableFormat.trim!(f, false)
    @test f.specifier[:trim] == TableFormat.Trim.no.value
    @test_throws ArgumentError TableFormat.trim!(f, "ccc")

end

@testset "locale" begin
    f = Format()
    @test isempty(f.locale)
    TableFormat.symbol_prefix!(f, "lll")
    @test haskey(f.locale, :symbol)
    @test f.locale[:symbol][1] == "lll"
    @test f.locale[:symbol][2] == ""
    TableFormat.symbol_prefix!(f, "rrr")
    @test f.locale[:symbol][1] == "rrr"
    @test f.locale[:symbol][2] == ""

    f = Format()
    @test isempty(f.locale)
    TableFormat.symbol_suffix!(f, "lll")
    @test haskey(f.locale, :symbol)
    @test f.locale[:symbol][1] == ""
    @test f.locale[:symbol][2] == "lll"
    TableFormat.symbol_suffix!(f, "rrr")
    @test f.locale[:symbol][1] == ""
    @test f.locale[:symbol][2] == "rrr"

    TableFormat.symbol_prefix!(f, "kkk")
    @test f.locale[:symbol][1] == "kkk"
    @test f.locale[:symbol][2] == "rrr"

    f = Format()
    TableFormat.decimal_delimiter!(f, '|')
    @test f.locale[:decimal] == "|"
    TableFormat.group_delimiter!(f, ';')
    @test f.locale[:group] == ";"

    TableFormat.groups!(f, [3,3,3])
    @test f.locale[:grouping] == [3,3,3]

    TableFormat.groups!(f, [5])
    @test f.locale[:grouping] == [5]

    @test_throws ArgumentError TableFormat.groups!(f, Int[])
    @test_throws ArgumentError TableFormat.groups!(f, Int[-1,2])
end

@testset "nully & prefix" begin
    f = Format()
    TableFormat.nully!(f, "gggg")
    @test f.nully == "gggg"
    TableFormat.si_prefix!(f, TableFormat.Prefix.femto)
    @test f.prefix == TableFormat.Prefix.femto.value
end

@testset "Format creation" begin
    f = Format(
        align = Align.left,
        fill = "+",
        group = true,
        padding = Padding.yes,
        padding_width = 10,
        precision = 2,
        scheme = Scheme.decimal,
        sign = Sign.negative,
        symbol = DSymbol.yes,
        trim = Trim.yes,

        symbol_prefix = "tt",
        symbol_suffix = "kk",
        decimal_delimiter = ";",
        group_delimiter = ",",
        groups = [2, 2],
        nully = "f",
        si_prefix = Prefix.femto
    )

    @test f.specifier == Dict{Symbol, Any}(
        :align => Align.left.value,
        :fill => "+",
        :group => Group.yes.value,
        :padding => Padding.yes.value,
        :width => 10,
        :precision => ".2",
        :type => Scheme.decimal.value,
        :sign => Sign.negative.value,
        :symbol => DSymbol.yes.value,
        :trim => Trim.yes.value
    )

    @test f.nully == "f"
    @test f.prefix == Prefix.femto.value

    @test f.locale[:symbol] == ["tt", "kk"]
    @test f.locale[:decimal] == ";"
    @test f.locale[:group] == ","
    @test f.locale[:grouping] == [2,2]

end


text:
using IEEE802_15_4
using ZMQ

sink  = PacketSink( BPSK, 5, diff_enc = true, verbosity = 0 )
ctx   = Context( 1 )
zsock = Socket( ctx, PULL )

connect( zsock, "tcp://localhost:5555" )

while true
    msg = recv( zsock )
    out = convert( IOStream, msg )
    bytes = takebuf_array( out )
    exec( sink, bytes )
end


text:
#
# Database and bindings.
#

type Database <: AbstractDatabase
    loccache::Dict{Type, BindingTable}
    globcache::BindingTable

    Database() = new(Dict{Type, BindingTable}(), NO_BINDING_TABLE)
end

const DB = Database()

globals(db::Database) = db.globcache
locals(db::Database, dom::Domain) =
    locals(db, dom.desc)

locals(db::Database, desc::Type) =
    get(db.loccache, desc, NO_BINDING_TABLE)

locals(db::Database, desc::Symbol) =
    locals(db, Entity{desc})

function locals(db::Database, desc::Vector{Output})
    tbl = NO_BINDING_TABLE
    for pos in eachindex(desc)
        fld = FieldQuery(Domain(desc), pos)
        tbl = assoc(tbl, (Symbol(pos), 0), QueryBinding(fld))
        tag = decoration(desc[pos], :tag, Symbol, Symbol(""))
        if tag != Symbol("")
            tbl = assoc(tbl, (tag, 0), QueryBinding(fld))
        end
    end
    return tbl
end

function attach!(db::Database, name::Symbol, prim::AbstractPrimitive)
    T = datatype(domain(input(prim)))
    tbl = get(db.loccache, T, NO_BINDING_TABLE)
    tbl = assoc(tbl, (name, 0), PrimitiveBinding(prim))
    db.loccache[T] = tbl
end

attach!(name::Symbol, prim::AbstractPrimitive) =
    attach!(DB, name, prim)

function attach!(db::Database, name::Symbol, q::Query)
    T = datatype(domain(input(q)))
    tbl = get(db.loccache, T, NO_BINDING_TABLE)
    tbl = assoc(tbl, (name, 0), QueryBinding(q))
    db.loccache[T] = tbl
end

attach!(name::Symbol, q::Query) =
    attach!(DB, name, q)

function attachglobal!(db::Database, name::Symbol, arity::Int, binding::AbstractBinding)
    db.globcache = assoc(db.globcache, (name, arity), binding)
end

attachglobal!(db::Database, name::Symbol, arity::Int=-1) =
    attachglobal!(db, name, arity, StaticBinding{name}())

attachglobal!(name::Symbol, arity::Int=-1) =
    attachglobal!(DB, name, arity)

immutable LiteralBinding <: AbstractBinding
    val
end

immutable StaticBinding{N} <: AbstractBinding
end

compile{N}(binding::StaticBinding{N}, scope::Scope, args::Vector{Syntax}) =
    compile(binding, scope, (compile(scope, arg) for arg in args)...)

immutable PrimitiveBinding <: AbstractBinding
    prim::AbstractPrimitive
end

compile(binding::PrimitiveBinding, scope::Scope, args::Vector{Syntax}) =
    begin
        @assert isempty(args)
        compile(binding, scope)
    end

compile(binding::PrimitiveBinding, scope::Scope) =
    Query(binding.prim) |> setnamespace(scope.db)

immutable QueryBinding <: AbstractBinding
    query::Query
end

compile(binding::QueryBinding, scope::Scope, args::Vector{Syntax}) =
    begin
        @assert isempty(args)
        q = binding.query
        if q.ns === NO_NAMESPACE
            ns = Namespace(scope.db, locals(scope.db, domain(q.ity)), locals(scope.db, domain(q.oty)), scope.globs)
            q = q |> setnamespace(ns)
        end
        q
    end

compile(syntax::Syntax) =
    compile(DB, syntax)

compile(db::AbstractDatabase, syntax::Syntax) =
    compile(Scope(db, Void), syntax)

function compile(scope::Scope, syntax::Syntax)
    binding =
        if !haslabel(syntax)
            QueryBinding(ConstQuery(syntax.val))
        else
            maybe_binding = lookup(scope, syntax.label, length(syntax.args))
            if isnull(maybe_binding)
                error("invalid identifier: $(syntax.label)")
            end
            get(maybe_binding)
        end
    return compile(binding, scope, syntax.args)
end

query(db::AbstractDatabase, syntax::Syntax, paramtypes::Vector{Pair}=Pair[]) =
    let scope = Scope(db, Void)
        globs = scope.globs
        for (tag, oty) in paramtypes
            tag  = Symbol(tag)
            oty = convert(Output, oty)
            globs = assoc(globs, (tag, 0), QueryBinding(SlotQuery(tag, oty) |> setnamespace(db)))
        end
        scope = Scope(scope.db, scope.dom, globs, scope.locs)
        cse(format(compile(scope, syntax)))
    end

query(db::AbstractDatabase, expr, paramtypes=Pair[]) =
    query(db, expr2syntax(expr), paramtypes)

query(syntax, paramtypes=Pair{Symbol}[]) =
    query(DB, syntax, paramtypes)

macro query(args...)
    L = endof(args)
    while L > 0 && isa(args[L], Expr) && args[L].head == :(::)
        L = L-1
    end
    @assert 1 <= L <= 2
    if L == 2
        db, expr = args[1], args[2]
    else
        db, expr = DB, args[1]
    end
    expr = expr2syntax(expr)
    params = [:( $(string(arg.args[1])) => $(arg.args[2]) ) for arg in args[L+1:end]]
    return quote
        query($(esc(db)), $expr, Pair[$(params...)])
    end
end

include("bind/static.jl")

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module SearchLightSessionAdapter

using Sessions, Genie, Logger, Genie.Configuration, App, Migration, JSON, SearchLight

type StorageSession <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  ### fields
  id::Nullable{SearchLight.DbId}
  name::String
  val::String

  ### constructor
  StorageSession(;
    id = Nullable{SearchLight.DbId}(),
    name = "",
    val = ""
  ) = new("storage_sessions", "id", id, name, val)
end

"""
    write(session::Sessions.Session) :: Sessions.Session

Persists the `Session` object to the DB, using the configured SearchLight DB and returns it.
"""
function write(session::Sessions.Session) :: Sessions.Session
  try
    SearchLight.update_by_or_create!!(StorageSession(name = session.id, val = JSON.json(session.data)), :name)
  catch ex
    Logger.log("Error when serializing session in $(@__FILE__):$(@__LINE__)", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    rethrow(ex)
  end

  session
end

"""
    read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
    read(session::Sessions.Session) :: Nullable{Sessions.Session}

Attempts to read from DB the session object.
"""
function read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
  try
    session_info = SearchLight.find_one_by!!(StorageSession, :name, session_id)
    session = Sessions.Session(session_info.name, JSON.parse(session_info.val))

    return isnull(session) ? Nullable{Sessions.Session}() : Nullable{Sessions.Session}(session)
  catch ex
    Logger.log("Can't read session", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    return Nullable{Sessions.Session}(write(Sessions.Session(session_id)))
  end
end
function read(session::Sessions.Session) :: Nullable{Sessions.Session}
  read(session.id)
end

### Migrations module

module Migrations

module CreateTableStorageSessions

import Migration: create_table, column, column_id, add_index, drop_table

using Configuration, App

"""

"""
function up()
  create_table(App.config.session_table) do
    [
      column_id()
      column(:name, :string, "UNIQUE", limit = 64)
      column(:val, :text)
    ]
  end

  add_index(App.config.session_table, :name)
end

function down()
  drop_table(App.config.session_table)
end

end # module CreateTableStorageSessions

end # module Migrations

end


text:
module Day21

using AdventOfCode2019

function day21(input::String = readInput(joinpath(@__DIR__, "..", "data", "day21.txt")))
    program = parse.(Int, split(input, ','))

    # logic: (¬A || ¬B || ¬C) && D
    instructionsPart1 = [
        "NOT A T",
        "NOT B J",
        "OR T J",
        "NOT C T",
        "OR T J",
        "AND D J",
        "WALK\n",
    ]

    # logic: (¬A || ¬B || (¬C && H)) && D
    instructionsPart2 = [
        "NOT A T",
        "NOT B J",
        "OR T J",
        "NOT C T",
        "AND H T",
        "OR T J",
        "AND D J",
        "RUN\n",
    ]

    p1 = run(program, instructionsPart1)
    p2 = run(program, instructionsPart2)

    return [p1, p2]
end

function run(program::Array{Int,1}, instructions)
    ascii = [Int(x) for x in join(instructions, "\n")]
    inp = Channel{Int}(100)
    out = Channel{Int}(1)

    t = @async AdventOfCode2019.IntCode.run_program!(copy(program), inp, out, nothing)
    for a in ascii
        put!(inp, a)
    end
    val = 0
    while !istaskdone(t)
        val = take!(out)
    end
    close(inp)
    close(out)
    return val
end

end # module


text:
using ShapeDistributions
using GeometryBasics
using LinearAlgebra
using Statistics
using Test
include("test_distances.jl")
include("test_sampling.jl")
include("test_shape_function.jl")
include("test_utils.jl")

@testset "ShapeDistributions.jl" begin
    test_point_sampling()
    test_shape_function()
end


text:
export AbstractGridWorldAction, MoveForward, TurnLeft, TurnRight
export MOVE_FORWARD, TURN_LEFT, TURN_RIGHT

#####
# Actions
#####

abstract type AbstractGridWorldAction end

struct MoveForward <: AbstractGridWorldAction end
const MOVE_FORWARD = MoveForward()

struct TurnRight <: AbstractGridWorldAction end
const TURN_RIGHT = TurnRight()

struct TurnLeft <: AbstractGridWorldAction end
const TURN_LEFT = TurnLeft()

(x::TurnRight)(::Left) = UP
(x::TurnRight)(::Up) = RIGHT
(x::TurnRight)(::Right) = DOWN
(x::TurnRight)(::Down) = LEFT
(x::TurnLeft)(::Left) = DOWN
(x::TurnLeft)(::Up) = LEFT
(x::TurnLeft)(::Right) = UP
(x::TurnLeft)(::Down) = RIGHT

struct Pickup <: AbstractGridWorldAction end
const PICK_UP = Pickup()

struct Drop <: AbstractGridWorldAction end
const DROP = Drop()