text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# This file contains all the global constants, variables, routines.

typealias BasicBlock      CompilerTools.CFGs.BasicBlock
typealias Statement       CompilerTools.CFGs.TopLevelStatement
typealias Liveness        CompilerTools.LivenessAnalysis.BlockLiveness
typealias CFG             CompilerTools.CFGs.CFG
typealias DomLoops        CompilerTools.Loops.DomLoops
typealias Loop            CompilerTools.Loops.Loop
typealias GenSymId        Int
typealias BasicBlockIndex Int
typealias StatementIndex  Int
typealias Sym             Union{Symbol, GenSym} # A Symbol or GenSym.
typealias Sym2TypeMap     Dict{Sym, Type}
typealias Symexpr         Union{Symbol, GenSym, Expr} # A Symbol, GenSym or Expr

# Options controlling debugging, performance (library choice, cost model), etc.
@doc """ Enable Sparse Accelerator """
const SA_ENABLE = 1

@doc """ Print verbose dump """
const SA_VERBOSE = 2

@doc """ Use Jula's default sparse matrix functions """
const SA_USE_JULIA = 4

@doc """ Use MKL sparse matrix functions """
const SA_USE_MKL = 8

@doc """
Use Sparse Matrix Pre-processing library (SpMP) functions. SPMP is a
high-performance parallel implementation of BFS/RCM reordering,
Gauss-Seidel smoother, sparse triangular solver, etc.
"""
const SA_USE_SPMP = 16

@doc """"
Pattern match and replace the code that is functionally equivalent to SpMV, dot,
WAXPBY, etc. with calls to the corresponding SPMP library functions.
"""
const SA_REPLACE_CALLS = 32

@doc """
Use splitting patterns during call replacement. Effective only when
SA_REPLACE_CALLS has been specified
"""
const SA_USE_SPLITTING_PATTERNS = 64

@doc """ Enable context-sensitive optimization. """
const SA_CONTEXT = 128

@doc """ Enable context-sensitive optimization for a whole function. """
const SA_CONTEXT_FUNC = 256

@doc """ Enable reordering of arrays. """
const SA_REORDER = 512

@doc """ Disable collective structure prediction. """
const SA_DISABLE_COLLECTIVE_STRUCTURE_PREDICTION = 1024

# The internal booleans corresponding to the above options, and their default values
sparse_acc_enabled                      = false
show_level                              = 256
use_Julia                               = false
use_MKL                                 = false
use_SPMP                                = true
replace_calls_enabled                   = false
use_splitting_patterns                  = false
reorder_enabled                         = false
reorder_when_beneficial                 = true # By default, reordering with benefit-cost analysis
context_sensitive_opt_enabled           = false
context_sensitive_opt_for_func          = false
collective_structure_prediction_enabled = true

# A control used only for trace call replacement internally
trace_call_replacement        = false

@doc """
Set options for Sparso. The arguments can be any one or more
of the following: SA_VERBOSE, SA_USE_JULIA, SA_USE_MKL, SA_USE_SPMP.
They can appear in any order, except that
SA_USE_JULIA, SA_USE_MKL and SA_USE_SPMP are exclusive with each other, and the
last one of them wins.
"""
function set_options(args...)
    for arg in args
        if arg == SA_ENABLE
            if !sparse_acc_enabled
                # Insert sparse accelerator as 1 pass into the optimization framework
                sparse_accelerator_pass = OptFramework.optPass(Sparso.entry, true)
                #OptFramework.setOptPasses([sparse_accelerator_pass])
                CompilerTools.OptFramework.addOptPass(sparse_accelerator_pass)
            end
            global sparse_acc_enabled = true
        elseif arg == SA_VERBOSE
            global show_level = 1
        elseif arg == SA_USE_JULIA
            global use_Julia = true; global use_MKL = false; global use_SPMP = false
        elseif arg == SA_USE_MKL
            global use_Julia = false; global use_MKL = true; global use_SPMP = false
        elseif arg == SA_USE_SPMP
            global use_Julia = false; global use_MKL = false; global use_SPMP = true
        elseif arg == SA_REPLACE_CALLS
            global replace_calls_enabled = true
        elseif arg == SA_USE_SPLITTING_PATTERNS
            global use_splitting_patterns = true
        elseif arg == SA_CONTEXT
            global context_sensitive_opt_enabled = true
        elseif arg == SA_CONTEXT_FUNC
            global context_sensitive_opt_for_func = true
        elseif arg == SA_REORDER
            # Reordering is a kind of context-sensitive optimization: without
            # know a matrix is constant valued, the library might not want to
            # do reordering.
            global context_sensitive_opt_enabled = true
            global reorder_enabled               = true
        elseif arg == SA_DISABLE_COLLECTIVE_STRUCTURE_PREDICTION
            global collective_structure_prediction_enabled = false # control Julia
            set_collective_structure_prediction(false) # control the library
        else
            # TODO: print usage info
        end
    end
end

@doc """
Properties of a variable/argument.
"""
const SA_CONST_VALUED         = 1
const SA_CONST_STRUCTURED     = 2
const SA_SYMM_VALUED          = 4
const SA_SYMM_STRUCTURED      = 8
const SA_STRUCTURE_ONLY       = 16
const SA_LOWER_OF             = 32
const SA_UPPER_OF             = 64
const SA_TRANSPOSE_OF         = 128
const SA_HAS_FREE_MEMORY      = 1024 # Used only in a pattern's skeleton
const SA_HAS_DEDICATED_MEMORY = 2048 # Difference from SA_HAS_FREE_MEMORY: not only the memory is free, but dedicated to an array, not shared/aliased with any other.
const SA_CONST_SIZED          = 4096
const SA_MAXIMAL_STRUCTURED   = 8192

# Create a path to libcsr. This is a CSR(Compressed Sparse Row format)-based
# interface to the SPMP library.
const libcsr = joinpath(dirname(@__FILE__), "..", "lib", "libcsr.so")

@doc """
Build a dictionary for the symbols, including GenSym's, to store the type info
of the typed func_ast. Note that all the symbols' scope is function-wise, even
though in the source code, they might appear to have local or nested scopes:
symbol renaming seems to have been done to make them function-wise.
"""
function build_symbol_dictionary(lambda :: LambdaInfo)
    symbol_info = Sym2TypeMap()

    # Record Symbols' types
    for i in lambda.var_defs
        symbol_info[i[2].name] = i[2].typ
    end

    # Record GenSym's types.
    for id = 1:length(lambda.gen_sym_typs)
        # Note that a GenSym id starts from 0
        symbol_info[GenSym(id - 1)] = lambda.gen_sym_typs[id]
    end
    symbol_info
end

@doc """
Determine the type of an AST node. A Symbol or GenSym gets a type from the
symbol_info. An expression or SymbolNode gets the type stored by Julia type
inference. All the other kinds of AST nodes resort to the default typeof().
"""
function type_of_ast_node(node, symbol_info :: Sym2TypeMap)
    local typ = typeof(node)
    if typ == Symbol
        # Use get() instead [] in case the key (like Symbol "*") does not exist
        # Return Void if no info found
        return get(symbol_info, node, Void)
    elseif typ == GenSym
        return get(symbol_info, node, Void)
    elseif typ == Expr || typ == SymbolNode
        return node.typ
    else
        return typ
    end
end

@doc """ Is the type a scalar (number), or an array? """
function number_or_array(typ :: Type)
    is_number = (typ <: Number)

    # We assume the user program does not use Range for any array
    # computation, although Range is a subtype of AbstractArray
    is_array  = (typ <: AbstractArray && !(typ <: Range))

    is_number, is_array
end

@doc """ Are the types scalars (numbers), or are some of them arrays? """
function numbers_or_arrays(result_type, arg_types :: Tuple)
    all_numbers, some_arrays = number_or_array(result_type)
    for t in arg_types
        is_number, is_array = number_or_array(t)
        all_numbers         = all_numbers && is_number
        some_arrays         = some_arrays || is_array
    end
    all_numbers, some_arrays
end

@doc """ A module (or function)'s name string """
function module_or_function_name(arg :: Any)
    if typeof(arg) == Symbol
        return string(arg)
    elseif typeof(arg) == QuoteNode && typeof(arg.value) == Symbol
        return string(arg.value)
    elseif isdefined(:GlobalRef) && typeof(arg) == GlobalRef
        return string(arg)
    elseif typeof(arg) == Expr && arg.head == :call && length(arg.args) == 3 &&
        arg.args[1] == TopNode(:getfield)
        from_name  = module_or_function_name(arg.args[2])
        field_name = module_or_function_name(arg.args[3])
        return string(from_name, ".", field_name)
    else
        throw(UnhandledModuleOrFunctionName(arg))
    end
end

@doc """
From the given call arguments, figure out the module and function name
"""
function resolve_call_names(call_args :: Vector)
    assert(length(call_args) > 0)
    module_name, function_name = "", ""
    if isdefined(:GlobalRef) && typeof(call_args[1]) == GlobalRef
        return string(call_args[1].mod), string(call_args[1].name)
    elseif typeof(call_args[1]) == Symbol
        # Example: :*
        function_name = string(call_args[1])
    elseif isa(call_args[1], TopNode) && length(call_args) == 3
        # Example: top(getfield), Sparso,:SpMV
        # special case: top(getfield), GenSym, ...
        if call_args[1] == TopNode(:getfield) && !isa(call_args[2], GenSym)
            module_name   = module_or_function_name(call_args[2])
            function_name = module_or_function_name(call_args[3])
        else
            function_name = module_or_function_name(call_args[1].name)
        end
    elseif  isa(call_args[1], TopNode) && call_args[1] == TopNode(:setfield!)
            function_name = "setfield!"
    elseif  isa(call_args[1], TopNode) && call_args[1] == TopNode(:typeof)
            function_name = "typeof"
    elseif  isa(call_args[1], TopNode) && call_args[1] == TopNode(:fieldtype)
            function_name = "fieldtype"
    elseif  isa(call_args[1], TopNode) && call_args[1] == TopNode(:tuple)
            function_name = "tuple"
    elseif  isa(call_args[1], TopNode) && call_args[1] == TopNode(:apply_type)
            function_name = "apply_type"
    elseif isa(call_args[1], Expr) && call_args[1].head == :call
        # The function may be represented in an Expr. Example 1:
        #           Expr call [SparseMatrixCSC{Float64,Int32}]
        #                Expr call [Type{SparseMatrixCSC{Float64,Int32}}]
        #                    top(apply_type) [TopNode]
        #                    Main.SparseMatrixCSC [GlobalRef]
        #                    Main.Cdouble [GlobalRef]
        #                    Main.Cint [GlobalRef]
        # In this apply_type case, the function can be treated as an assignment.
        # Example 2: (:call, top(getfield), Sparso,:SpMV)
        return resolve_call_names(call_args[1].args)
    end
    module_name, function_name
end

@doc """
Search in a database for information of the function specified by the module
name, function name, and argument types.
"""
function look_for_function(
    database       :: Vector,
    module_name    :: AbstractString,
    function_name  :: AbstractString,
    argument_types :: Tuple
)
    for item in database
        if module_name    == item.module_name   &&
           function_name  == item.function_name &&
           length(argument_types) == length(item.argument_types)
            found = true
            for i in 1:length(argument_types)
               if !(argument_types[i] <: item.argument_types[i])
                   found = false
                   break
                end
            end
            if found
               return item
            end
        end
    end
    return nothing
end

@doc """
For a SymbolNode, return its symbol (So that, for example, in building a dictionary,
a variable can be indexed by only a Symbol, instead of by a Symbol sometimes,
and SymbolNode sometimes). Otherwise, return the original input.
"""
get_symexpr(s) = (typeof(s) == SymbolNode) ? s.name : s

@doc """ Abstract type for a pattern to match and/or replace AST nodes. """
abstract Pattern

@doc """ Abstract type for an action to take in transforming CFG. """
abstract Action

@doc """
Abstract type for a region. A region can be an arbitrary part of user code.
"""
abstract Region

@doc """ Dummy function that does nothing. Used for empty call back. """
function do_nothing()
end

@doc """
A call site of interesting functions (like SpMV, triangular solver, etc.).
"""
type CallSite
    ast :: Expr
end

@doc """"
Properties of a matrix in a region (The properties may be different in another
region).

constant_valued       : The matrix is a constant in value(and thus of course
                        constant in structure).
constant_structured   : The matrix has always the same structure, even though its
                        value may change.
is_symmetric          : The matrix is symmetric in value (and of course symmetric
                        in structure
is_structure_symmetric: The matrix is symmetric in structure.
is_structure_only     : Only the structure of matrix is to be used.
is_single_def         : The matrix is statically defined only once.
has_dedicated_memory  : The matrix (actually for vector as well) has a decidated
                        memory space, so no worry of aliases or memory allocation.
"""
type MatrixProperties
    constant_valued        :: Bool
    constant_structured    :: Bool
    is_symmetric           :: Bool
    is_structure_symmetric :: Bool
    is_structure_only      :: Bool
    is_single_def          :: Bool
    has_dedicated_memory   :: Bool
    lower_of               :: Any
    upper_of               :: Any
    transpose_of           :: Any

    MatrixProperties() = new(false, false, false, false, false, false, false, nothing, nothing, nothing)
end
typealias Symexpr2PropertiesMap Dict{Symexpr, MatrixProperties}

@doc """
Call sites of interesting functions (like SpMV, triangular solver, etc.). The
function's result and argument types are figured out with the help of symbol_info.
Some patterns may need matrix properties in order to match.
The patterns to match are specified by the specific analysis. In addition to the
direct change of the call site ASTs due to replacement, there might be additional
actions resulted (like hoisting some computation out of a loop, etc.)
Additional information specific to each analysis can be attached to the "extra"
field.
"""
type CallSites
    sites       :: Set{CallSite}
    region      :: Region
    lambda      :: LambdaInfo
    symbol_info :: Sym2TypeMap
    liveness    :: Liveness
    patterns    :: Vector{Pattern}
    actions     :: Vector{Action}
    extra       :: Any
end

@doc """ Insert new statements to a basic block before or after a statement. """
type InsertBeforeOrAfterStatement <: Action
    new_stmts   :: Vector{Statement}
    before      :: Bool
    bb          :: BasicBlock
    stmt_idx    :: StatementIndex
end

@doc """ Insert new statements before a loop's head block """
type InsertBeforeLoopHead <: Action
    new_stmts    :: Vector{Statement}
    loop         :: Loop # The loop.
    outside_loop :: Bool # Outside_loop == true/false would make the statements
                         # inserted outside/inside the loop.
end

@doc """ Insert new statements on an edge """
type InsertOnEdge <: Action
    new_stmts    :: Vector{Statement}
    from_bb      :: BasicBlock
    to_bb        :: BasicBlock
end

@doc """ The total number of new symbols generated """
new_symbol_counter = 0

@doc """ Generate a new symbol with the given string name + new_symbol_counter """
function new_symbol(s::ASCIIString)
    x = symbol(string("__", s, string(new_symbol_counter), "__"))
    global new_symbol_counter += 1
    x
end

@doc """
Entry of Sparso.
"""
function entry(func_ast :: Expr, func_arg_types :: Tuple, func_args)
    old_ast = copy(func_ast)
    new_ast = nothing
    try
        dprintln(1, 0, "******************************* Sparso ******************************")
        dprintln(1, 0, "Signature:")
        for i = 1 : length(func_args)
            if i == 1
                dprint(1, 1, "(", func_args[i], "::", func_arg_types[i])
            else
                dprint(1, 0, ", ", func_args[i], "::", func_arg_types[i])
            end
        end
        dprintln(1, 0, ")\n\nAST:")
        dprintln(1, 1, func_ast)

        assert(func_ast.head == :lambda)

        # Build the common facilities: symbols' type dictionary, liveness
        # info, and control flow graph.
        LivenessAnalysis.set_use_inplace_naming_convention()
        lambda      = lambdaExprToLambdaInfo(func_ast)
        symbol_info = build_symbol_dictionary(lambda)
        liveness    = LivenessAnalysis.from_expr(func_ast, no_mod = LivenessAnalysis.create_unmodified_args_dict())
        cfg         = liveness.cfg
        loop_info   = Loops.compute_dom_loops(cfg)

        dprintln(1, 0, "\nFunction body showing structures:")
        dsprintln(1, 1, symbol_info, liveness, func_ast)

        func_region = FunctionRegion(func_ast)
        regions, bb2depth = region_formation(func_region, cfg, loop_info)

        if context_sensitive_opt_enabled
            # Reordering and context-sensitive optimization: Do all analyses, and
            # put their intended transformation actions into a list. Then
            # transform the code with the list on the CFG.
            # This phase may change the AST. It does not change the CFG, but
            # records the actions for changing the CFG.
            # REQUIREMENT: even though the AST and CFG can be changed, liveness
            # information regarding each statement that originally exists
            # should still remain valid. Even though the expression contained in
            # an original statement can be changed, the statement (a container)
            # should not be touched. Therefore, we should still consult liveness
            # info based on the original statement, no matter the expression in
            # it has been changed or not. This is important: otherwise, we would
            # have to rebuild liveness in order for subsequent optimizations (
            # like call replacement) to continue.
            actions = Vector{Action}()
            actions = AST_context_sensitive_transformation(actions, func_region, regions, bb2depth, lambda, symbol_info, liveness, cfg)
            CFG_transformation(actions, cfg)
        end

        # Call replacement happens at the last, as we would like context-
        # optimizations, which are more important, to happen first.
        # Call replacement relies only on type info and liveness, and it affects
        # only one or two statements. It should not change liveness of statements
        # either.
        if replace_calls_enabled
            actions = replace_calls(func_region, regions, lambda, symbol_info, liveness, cfg)
            CFG_transformation(actions, cfg)
        end

        # Now create a new function based on the CFG
        body_reconstructed   = CFGs.createFunctionBody(cfg)
         func_ast.args[3].args = body_reconstructed
         new_ast = func_ast
        #new_ast = lambdaInfoToLambdaExpr(lambda, Expr(:body, body_reconstructed...))

        dprintln(1, 0, "\nNew AST:")
        dprintln(1, 1, new_ast)
        dprintln(1, 0, "********************************************************************************")
        Libc.flush_cstdio()
        flush(STDOUT)
    catch ex
        # In case any exception happen in the printing, try
        try
            Libc.flush_cstdio()
            flush(STDOUT)
            dprintln(256, 0, "Exception! Sparse Accelerator skips optimizing the call.")
            dprintln(256, 1, ex)
            dprintln(256, 0, "********************************************************************************")
            Libc.flush_cstdio()
            flush(STDOUT)
        catch
            # Do nothing
        end

        # Return the original AST without any change
        new_ast = old_ast
    finally
        return new_ast
    end
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
function LagSeries(j::Int64,z::Float64,x::AbstractVector{T}) where T # need to investigate stability
  c = Lag(j,z) # c[2] gives L_0^(1), c[i] gives L_{i-2}^(1)
  v0 = x
  out = zeros(Complex{Float64},length(x),j)
  out[1:end,1] = v0*c[2]
  @inbounds for i = 2:j
    out[1:end,i] = (1.0 .+ x).*out[1:end,i-1] .+ x*(c[i+1]-c[i])
  end

  return out
end

# old version, still here for testing
function LagSeries_old(z::Float64,x::AbstractVector{T},cfs::AbstractVector{S}) where {T,S} # need to investigate stability
  j = length(cfs)
  c = Lag(j,z) # c[2] gives L_0^(1), c[i] gives L_{i-2}^(1)
  v0 = x
  #out = zeros(Complex{Float64},length(x))
  out = v0*c[2]
  pm = -1
  sum = (pm*cfs[1])*out
  @inbounds for i = 2:j
    pm = pm*(-1)
    out = (1.0 .+ x).*out .+ x*(c[i+1]-c[i])
    sum = sum + (pm*cfs[i])*out
  end
  return sum
end

function LagSeries(z::Complex{Float64},x::Vector{Complex{Float64}},cfs::Vector{Complex{Float64}}) # need to investigate stability
  j = length(cfs)
  #c = Lag(j,z) # c[2] gives L_0^(1), c[i] gives L_{i-2}^(1)
  #v0 = x
  #out = zeros(Complex{Float64},length(x))
  c0 = 0.
  c1 = exp(-z/2) # c[2]
  out = x*c1
  k = 0
  pm = -1
  sum = (pm*cfs[1])*out
  @inbounds for i = 2:j
    pm = pm*(-1)
    c2 = (2*k+2-z)*c1/(k+1) - c0 # c2 = c[i + 1]
    #out = (1.0 .+ x).*out .+ x*(c2-c1)
    out .*= (1.0 .+ x)
    out .+= x*(c2-c1)
    #@inbounds sum = sum + (pm*cfs[i])*out
    axpy!(pm*cfs[i],out,sum)
    k += 1
    c0 = c1
    c1 = c2
  end
  return sum
end

function LagSeries(z::BigFloat,x::Vector{Complex{BigFloat}},cfs::Vector{Complex{BigFloat}}) # need to investigate stability
  j = length(cfs)
  #c = Lag(j,z) # c[2] gives L_0^(1), c[i] gives L_{i-2}^(1)
  #v0 = x
  #out = zeros(Complex{Float64},length(x))
  c0 = 0.
  c1 = exp(-z/2) # c[2]
  out = x*c1
  k = 0
  pm = -1
  sum = (pm*cfs[1])*out
  @inbounds for i = 2:j
    pm = pm*(-1)
    c2 = (2*k+2-z)*c1/(k+1) - c0 # c2 = c[i + 1]
    #out = (1.0 .+ x).*out .+ x*(c2-c1)
    out .*= (1.0 .+ x)
    out .+= x*(c2-c1)
    #@inbounds sum = sum + (pm*cfs[i])*out
    axpy!(pm*cfs[i],out,sum)
    k += 1
    c0 = c1
    c1 = c2
  end
  return sum
end

function Res(j::Int64,α::Float64,β::Float64,z::Vector{Complex{Float64}})
  x = -2im*sign(j)*β./(z.+1im*sign(j)*β)
  y = -2*sign(j)*α*β
  return -LagSeries(abs(j),y,x)
end

function Res(j::Integer,α::Complex{Float64},β::Float64,z::Vector{Float64},cfs::Vector{Complex{Float64}})
  x = -2im*sign(j)*β./(z.+1im*sign(j)*β)
  y = -2*sign(j)*α*β
  return -LagSeries(y,x,cfs)
end

function Res(j::Integer,α::BigFloat,β::BigFloat,z::Vector{Complex{BigFloat}},cfs::Vector{Complex{BigFloat}})
  x = -2im*sign(j)*β./(z.+1im*sign(j)*β)
  y = -2*sign(j)*α*β
  return -LagSeries(y,x,cfs)
end

function CauchyPNO(f::Fun{OscLaurent{DD,RR}}) where {DD,RR <: BigFloat}# returns the non-oscillatory portion, assumes decay
    α = space(f).exp
    sp = OscLaurent(domain(space(f)),0.)
    if α ≈ 0.
      b = copy(f.coefficients)
      b[[1;2:2:end]] = zeros(typeof(f.coefficients[1]),length(b[[1;2:2:end]])) # zero positive
      b[1] = -fsum(b)
      return Fun(sp,b)
    elseif real(α) < 0. && length(f.coefficients) >= 3
      #m = length(f.coefficients[3:2:end])
      #c = Res(m,α,domain(space(f)).L,points(f))
      return Fun(sp,ApproxFun.transform(sp,-Res(1,α,BigFloat(domain(space(f)).L),points(f),f.coefficients[3:2:end])))
      #return Fun(sp,ApproxFun.transform(sp,-c*([(-1)^i for i = 1:m].*f.coefficients[3:2:end])))
    elseif real(α) <0.
      return Fun(sp,[0.0im,0.0im])
    elseif length(f.coefficients) >= 2
      #m = length(f.coefficients[2:2:end])
      #c = Res(-m,α,domain(space(f)).L,points(f))
      return Fun(sp,ApproxFun.transform(sp,Res(-1,α,BigFloat(domain(space(f)).L),points(f),f.coefficients[2:2:end])))
      #return Fun(sp,ApproxFun.transform(sp,c*([(-1)^i for i = 1:m].*f.coefficients[2:2:end])))
    else
      return Fun(sp,[0.0im])
    end
end

function CauchyPNO(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}# returns the non-oscillatory portion, assumes decay
    α = space(f).exp
    sp = OscLaurent(domain(space(f)),0.)
    if α ≈ 0.
      b = copy(f.coefficients)
      b[[1;2:2:end]] = zeros(typeof(f.coefficients[1]),length(b[[1;2:2:end]])) # zero positive
      b[1] = -fsum(b)
      return Fun(sp,b)
    elseif real(α) < 0. && length(f.coefficients) >= 3
      #m = length(f.coefficients[3:2:end])
      #c = Res(m,α,domain(space(f)).L,points(f))
      return Fun(sp,ApproxFun.transform(sp,-Res(1,α,domain(space(f)).L,points(f),f.coefficients[3:2:end])))
      #return Fun(sp,ApproxFun.transform(sp,-c*([(-1)^i for i = 1:m].*f.coefficients[3:2:end])))
    elseif real(α) <0.
      return Fun(sp,[0.0im,0.0im])
    elseif length(f.coefficients) >= 2
      #m = length(f.coefficients[2:2:end])
      #c = Res(-m,α,domain(space(f)).L,points(f))
      return Fun(sp,ApproxFun.transform(sp,Res(-1,α,domain(space(f)).L,points(f),f.coefficients[2:2:end])))
      #return Fun(sp,ApproxFun.transform(sp,c*([(-1)^i for i = 1:m].*f.coefficients[2:2:end])))
    else
      return Fun(sp,[0.0im])
    end
end

function CauchyP(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
  α = space(f).exp
  if real(α) <= 0.
    return CauchyPNO(f)
  else
    return CauchyPNO(f) + f
  end
end

function CauchyP_SumFun(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
  α = space(f).exp
  if real(α) <= 0.
    return SumFun(CauchyPNO(f))
  else
    return SumFun([CauchyPNO(f),f])
  end
end

function CauchyM(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
  α = space(f).exp
  if real(α) >= 0.
    return CauchyMNO(f)
  else
    return CauchyMNO(f) - f
  end
end

function CauchyM_SumFun(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
  α = space(f).exp
  if real(α) >= 0.
    return SumFun(CauchyMNO(f))
  else
    return SumFun([CauchyMNO(f),-f])
  end
end

function CauchyMNO(f::Fun{OscLaurent{DD,RR}}) where {DD,RR <: BigFloat}# returns the non-oscillatory portion
    α = space(f).exp
    sp = OscLaurent(domain(space(f)),0.)
    if α ≈ 0.
      b = -copy(f.coefficients)
      b[[1;3:2:end]] = zeros(typeof(f.coefficients[1]),length(b[[1;3:2:end]])) # zero positive
      b[1] = -fsum(b)
      return Fun(sp,b)
    elseif real(α) < 0. && length(f.coefficients) >= 3
      #m = length(f.coefficients[3:2:end])
      #c = Res(m,α,domain(space(f)).L,points(f))
      #return Fun(sp,ApproxFun.transform(sp,-c*([(-1)^i for i = 1:m].*f.coefficients[3:2:end])))
      return Fun(sp,ApproxFun.transform(sp,-Res(1,α,BigFloat(domain(space(f)).L),points(f),f.coefficients[3:2:end])))
    elseif real(α) < 0.
      return Fun(sp,[0.0im,0.0im])
    elseif length(f.coefficients) >= 2
      #m = length(f.coefficients[2:2:end])
      #c = Res(-m,α,domain(space(f)).L,points(f))
      #return Fun(sp,ApproxFun.transform(sp,c*([(-1)^i for i = 1:m].*f.coefficients[2:2:end])))
      return Fun(sp,ApproxFun.transform(sp,Res(-1,α,BigFloat(domain(space(f)).L),points(f),f.coefficients[2:2:end])))
    else
      return Fun(sp,[0.0im])
    end
end

function CauchyMNO(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}# returns the non-oscillatory portion
    α = space(f).exp
    sp = OscLaurent(domain(space(f)),0.)
    if α ≈ 0.
      b = -copy(f.coefficients)
      b[[1;3:2:end]] = zeros(typeof(f.coefficients[1]),length(b[[1;3:2:end]])) # zero positive
      b[1] = -fsum(b)
      return Fun(sp,b)
    elseif real(α) < 0. && length(f.coefficients) >= 3
      #m = length(f.coefficients[3:2:end])
      #c = Res(m,α,domain(space(f)).L,points(f))
      #return Fun(sp,ApproxFun.transform(sp,-c*([(-1)^i for i = 1:m].*f.coefficients[3:2:end])))
      return Fun(sp,ApproxFun.transform(sp,-Res(1,α,domain(space(f)).L,points(f),f.coefficients[3:2:end])))
    elseif real(α) < 0.
      return Fun(sp,[0.0im,0.0im])
    elseif length(f.coefficients) >= 2
      #m = length(f.coefficients[2:2:end])
      #c = Res(-m,α,domain(space(f)).L,points(f))
      #return Fun(sp,ApproxFun.transform(sp,c*([(-1)^i for i = 1:m].*f.coefficients[2:2:end])))
      return Fun(sp,ApproxFun.transform(sp,Res(-1,α,domain(space(f)).L,points(f),f.coefficients[2:2:end])))
    else
      return Fun(sp,[0.0im])
    end
end

for cauchy in (:CauchyP,:CauchyM)
  @eval begin
    function $cauchy(f::Fun{T}) where {T<:SumSpace}
      sum(map(x -> $cauchy(x),components(f)))
    end

    function $cauchy(f::Fun{T}) where {T<:PiecewiseSpace}
      sum(map(x -> $cauchy(x),components(f)))
    end

    function $cauchy(f::Fun{T}) where {T<:ApproxFun.ArraySpace{J,1}} where {J}
      Fun(map($cauchy,Array(f)))
    end
  end
end

abstract type SIOperator{S,OT,T} <: Operator{T} end

macro SI_operator(Op)
    ConcOp = Meta.parse("Concrete"*string(Op))
    OpWrap = Meta.parse(string(Op)*"Wrapper")
    return esc(quote
        abstract type $Op{SSS,OT,TTT} <: SIOperator{SSS,OT,TTT} end

        struct $ConcOp{S<:Space,OT,T} <: $Op{S,OT,T}
            space::S
            sign::OT
        end

        struct $OpWrap{BT<:Operator,S<:Space,OT,T} <: $Op{S,OT,T}
            op::BT
            sign::OT
        end

        @wrapper $OpWrap

        $ConcOp(sp::Space,k) = $ConcOp{typeof(sp),typeof(k),prectype(sp)}(sp,k)
        $ConcOp(sp::Space) = $ConcOp(sp,1)

        $Op(sp::UnsetSpace,k) = $ConcOp(sp,k)
        $Op(sp::UnsetSpace,k::Integer) = $ConcOp(sp,k)

        # not needed yet
        function Base.convert(::Type{Operator{T}},D::$ConcOp) where T
            if T==eltype(D)
            D
            else
                $ConcOp{typeof(D.space),typeof(D.sign),T}(D.space,D.sign)
            end
        end

        function Base.convert(::Type{Operator{T}},D::$OpWrap) where T
            if T==eltype(D)
            D
            else
                # work around typeinfernece bug
                op=convert(Operator{T},D.op)
                $OpWrap{typeof(op),typeof(domainspace(op)),T}(op)
            end
        end

        ApproxFunBase.domain(D::$ConcOp) = domain(D.space)
        ApproxFunBase.domainspace(D::$ConcOp) = D.space
        Base.getindex(::$ConcOp{UnsetSpace,T},k::Integer,j::Integer) where {OT,T} =
            error("Spaces cannot be inferred for operator")
        ApproxFunBase.rangespace(D::$ConcOp{UnsetSpace,OT,T}) where {OT,T} = UnsetSpace()

        ApproxFunBase.promotedomainspace(D::$Op,sp::UnsetSpace) = D
        function ApproxFunBase.promotedomainspace(D::$Op,sp::Space)
            if isambiguous(domain(sp))
                $Op(typeof(sp)(domain(D)),D.sign)
            else
                $Op(sp,D.sign)
            end
        end
    end)
end

#Cauchy(σ,f) = σ == 1 ? CauchyP(f) : CauchyM(f)
@SI_operator(CauchyOperator)

function rangespace(D::ConcreteCauchyOperator{S,OT,T}) where {S<:OscLaurent,OT,T}
    α = D.space.exp
    L = D.space.domain.L
    return SumSpace(OscLaurent(α,L),OscLaurent(0.0,L))
end
bandwidths(D::ConcreteCauchyOperator{OscLaurent{DD,RR},OT,T}) where {DD,RR,OT,T} = ∞,∞
getindex(D::ConcreteCauchyOperator{OscLaurent{DD,RR},OT,T},k::Integer,j::Integer) where {DD,RR,OT,T} = one(T)

##TODO: Finish implementation of CauchyOperator

abstract type AbstractCauchyOperator end

struct Cauchy <: AbstractCauchyOperator
  pm::Integer
end

struct EvaluateCauchy <: AbstractCauchyOperator
  pm::Integer
  x::Complex{Float64}
end

Cauchy(pm,x::Complex{Float64}) = EvaluateCauchy(pm,x)
Cauchy(pm,x::Float64) = Cauchy(pm,convert(Complex{Float64},x))

function *(C::EvaluateCauchy,f::Fun{OscLaurent{DD,RR}}) where {DD,RR}
    if f.space.exp != 0.0
      @warn "non-zero exponent, cauchy eval is invalid"
    end
    S = f.space
    zz = mappoint(domain(S),Circle(),C.x)
    if abs(zz) < 1.0 || (abs(zz) ≈ 1 && C.pm == 1)
      return horner(f.coefficients,1:2:length(f.coefficients),zz)
    elseif abs(zz) > 1.0 || (abs(zz) ≈ 1 && C.pm == -1)
      invz = 1/zz
      return horner(f.coefficients,2:2:length(f.coefficients),invz)*invz + f.coefficients[1]
    elseif isnan(zz)
      return  f.coefficients[1]
    end
end

*(C::Cauchy,F::Fun) = (C.pm == 1) ? CauchyP(F) : CauchyM(F)
*(C::Cauchy,F::SumFun) = (C.pm == 1) ? +(map(CauchyP_SumFun,F.funs)...) : +(map(CauchyM_SumFun,F.funs)...)
*(C::Cauchy,F::Array{T,1}) where T = map(x -> C*x, F)

ext: jl


text:
_defaultNamesMixtures(N::Int) = ((Symbol[Symbol("c$i") for i in 1:N])...,)

"""
$(TYPEDEF)

A `Mixture` object for use with either a `<: AbstractPrior` or `<: AbstractRelative`.

Notes
- The internal data representation is a `::NamedTuple`, which allows total type-stability for all component types.
- Various construction helpers can accept a variety of inputs, including `<: AbstractArray` and `Tuple`.
- `N` is the number of components used to make the mixture, so two bumps from two Normal components means `N=2`.

DevNotes
- FIXME swap API order so Mixture of distibutions works like a distribtion, see Caesar.jl #808
  - Should not have field mechanics.
- TODO on sampling see #1099 and #1094 and #1069

Example
```juila
# prior factor
msp = Mixture(Prior,
              [Normal(0,0.1), Uniform(-pi/1,pi/2)],
              [0.5;0.5])

addFactor!(fg, [:head], msp, tags=[:MAGNETOMETER;])

# Or relative
mlr = Mixture(LinearRelative,
              (correlator=AliasingScalarSampler(...), naive=Normal(0.5,5), lucky=Uniform(0,10)),
              [0.5;0.4;0.1])

addFactor!(fg, [:x0;:x1], mlr)
```
"""
struct Mixture{N, F<:AbstractFactor, S, T<:Tuple} <: AbstractFactor
  """ factor mechanics """
  mechanics::F
  components::NamedTuple{S,T}
  diversity::Distributions.Categorical
  """ dimension of factor, so range measurement would be dims=1 """
  dims::Int
  labels::Vector{Int}
end

Mixture(f::Type{F},
        z::NamedTuple{S,T},
        c::Distributions.DiscreteNonParametric ) where {F<:AbstractFactor, S, T} = Mixture{length(z),F,S,T}(f(LinearAlgebra.I), z, c, size( rand(z[1],1), 1), zeros(Int, 0))
Mixture(f::F,
        z::NamedTuple{S,T},
        c::Distributions.DiscreteNonParametric ) where {F<:AbstractFactor, S, T} = Mixture{length(z),F,S,T}(f, z, c, size( rand(z[1],1), 1), zeros(Int, 0))
Mixture(f::Union{F,Type{F}},z::NamedTuple{S,T},
        c::AbstractVector{<:Real}) where {F<:AbstractFactor,S,T} = Mixture(f, z, Categorical([c...]) )
Mixture(f::Union{F,Type{F}},
        z::NamedTuple{S,T},
        c::NTuple{N,<:Real}) where {N,F<:AbstractFactor,S,T} = Mixture(f, z, [c...] )
Mixture(f::Union{F,Type{F}},
        z::Tuple,
        c::Union{<:Distributions.DiscreteNonParametric, <:AbstractVector{<:Real}, <:NTuple{N,<:Real}} ) where {F<:AbstractFactor, N} = Mixture(f,NamedTuple{_defaultNamesMixtures(length(z))}(z), c )
Mixture(f::Union{F,Type{F}},
        z::AbstractVector{<:SamplableBelief},
        c::Union{<:Distributions.DiscreteNonParametric, <:AbstractVector{<:Real}, <:NTuple{N,<:Real}} ) where {F <: AbstractFactor, N} = Mixture(f,(z...,), c )

function Base.resize!(mp::Mixture, s::Int)
  resize!(mp.labels, s)
end

_lengthOrNothing(val) = length(val)
_lengthOrNothing(val::Nothing) = 0

getManifold(m::Mixture) = getManifold(m.mechanics)

# TODO make in-place memory version
function sampleFactor( cf::CalcFactor{<:Mixture}, N::Int=1)
  #
  # TODO consolidate #927, case if mechanics has a special sampler
  # TODO slight bit of waste in computation, but easiest way to ensure special tricks in s.mechanics::F are included
  ## example case is old FluxModelsPose2Pose2 requiring velocity
  # FIXME better consolidation of when to pass down .mechanics, also see #1099 and #1094 and #1069

  cf_ = CalcFactor( cf.factor.mechanics, cf.metadata, 0, _lengthOrNothing(cf._legacyMeas), cf._legacyMeas, cf._legacyParams, cf._allowThreads, cf.cache)
  smpls = [getSample(cf_) for _=1:N]
    # smpls = Array{Float64,2}(undef,s.dims,N)
  #out memory should be right size first
  length(cf.factor.labels) != N ? resize!(cf.factor.labels, N) : nothing
  cf.factor.labels .= rand(cf.factor.diversity, N)
  for i in 1:N
    mixComponent = cf.factor.components[cf.factor.labels[i]]
    # measurements relate to the factor's manifold (either tangent vector or manifold point)
    setPointsMani!(smpls[i], rand(mixComponent,1))
  end

  # TODO only does first element of meas::Tuple at this stage, see #1099
  smpls
end

function DistributedFactorGraphs.isPrior(::Mixture{N, F, S, T}) where {N,F,S,T}
  return F <: AbstractPrior
end

"""
$(TYPEDEF)

Serialization type for `Mixture`.
"""
Base.@kwdef mutable struct PackedMixture <: AbstractPackedFactor
  N::Int
  # store the packed type for later unpacking
  F_::String
  S::Vector{String}
  components::Vector{PackedSamplableBelief}
  diversity::PackedSamplableBelief
end

function convert(::Type{<:PackedMixture}, obj::Mixture{N,F,S,T}) where {N,F,S,T}
  allcomp = PackedSamplableBelief[]
  for val in obj.components
    dtr_ = convert(PackedSamplableBelief, val)
    # FIXME ON FIRE, likely to be difficult for non-standard "Samplable" types -- e.g. Flux models in RoME
    push!(allcomp, dtr_)
  end
  pm = DFG.convertPackedType(obj.mechanics)
  pm_ = convert(pm, obj.mechanics)
  sT = string(typeof(pm_))
  dvst = convert(PackedSamplableBelief, obj.diversity)
  PackedMixture( N, sT, string.(collect(S)), allcomp, dvst )
end

function convert(::Type{<:Mixture}, obj::PackedMixture)
  N = obj.N
  F1 = getfield(Main, Symbol(obj.F_))
  S = (Symbol.(obj.S)...,)
  F2 = DFG.convertStructType(F1)

  components = map(c->convert(SamplableBelief,c), obj.components)
  diversity = convert(SamplableBelief, obj.diversity)
  # tupcomp = (components...,)
  ntup = NamedTuple{S}(components) # ,typeof(tupcomp)
  Mixture(F2, ntup, diversity)
end

  #


text:
using AlphaShapes
using Test

@testset "Test Volume" begin
   p = [0. 1.;1. 0.;1. 1.]
   @test AlphaShapes.SimplexVolume(p) == 0.5
end


text:
#####
##### update nutrient fields
#####

##### apply tendency to nutrient field
@kernel function apply_tendency_kernel!(nut, Gc, consume, ΔT, g::AbstractGrid)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds nut[ii, jj, kk] += Gc[ii, jj, kk] + consume[ii, jj, kk] /volume(ii, jj, kk, g)
end
function apply_tendency!(nuts, Gcs, consume, ΔT, g::AbstractGrid, arch::Architecture)
    kernel! = apply_tendency_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(nuts[name].data, Gcs[name].data, consume[name].data, ΔT, g, dependencies=barrier)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

function nut_update!(nutrients, Gcs, nut_temp, arch::Architecture, g::AbstractGrid, params, vel, consume, ΔT, iter)
    ##### compute advection tendency
    nut_advection!(nutrients, nut_temp, Gcs, vel, g, ΔT, arch)

    ##### compute nutrient diffusion,for each time step
    nut_diffusion!(Gcs, arch, g, nutrients, params["κh"], ΔT)

    ##### compute biogeochemical forcings of nutrients,for each time step
    zero_fields!(nut_temp)
    nut_forcing!(Gcs, nut_temp, nutrients, params, ΔT)

    ##### apply boundary conditions
    apply_bcs!(Gcs, nutrients, g, iter, ΔT, arch)

    ##### apply diffusion and forcing tendency
    apply_tendency!(nutrients, Gcs, consume, ΔT, g, arch)

    fill_halo_nut!(nutrients, g)
end

function nut_update!(nutrients, Gcs, nut_temp, arch::Architecture, g::AbstractGrid, params, consume, ΔT, iter)
    # compute biogeochemical forcings of nutrients,for each time step
    nut_forcing!(Gcs, nut_temp, nutrients, params, ΔT)

    ##### apply boundary conditions
    apply_bcs!(Gcs, nutrients, g, iter, ΔT, arch)

    # apply forcing tendency
    apply_tendency!(nutrients, Gcs, consume, ΔT, g, arch)

    fill_halo_nut!(nutrients, g)
end


text:
using ModelingToolkit
using ModelingToolkit: value
using Test

@parameters t
@variables x(t) y(t) z(t)

null_op = 0*t
@test isequal(simplify(null_op), 0)

one_op = 1*t
@test isequal(simplify(one_op), t)

identity_op = Num(Term(identity,[x.val]))
@test isequal(simplify(identity_op), x)

minus_op = -x
@test isequal(simplify(minus_op), -1x)
simplify(minus_op)

@variables x

@test toexpr(expand_derivatives(Differential(x)((x-2)^2))) == :($(*)(2, $(+)(-2, x)))
@test toexpr(expand_derivatives(Differential(x)((x-2)^3))) == :($(*)(3, $(^)($(+)(-2, x), 2)))
@test toexpr(simplify(x+2+3)) == :($(+)(5, x))

d1 = Differential(x)((-2 + x)^2)
d2 = Differential(x)(d1)
d3 = Differential(x)(d2)

@test toexpr(expand_derivatives(d3)) == :(0)
@test toexpr(simplify(x^0)) == :(1)

@test ModelingToolkit.substitute(value(2x + y == 1), Dict(x => 0.0, y => 0.0)) === false
@test ModelingToolkit.substitute(value(2x + y == 1), Dict(x => 0.0, y => 1.0)) === true

# 699
using SymbolicUtils: substitute
@parameters t a(t) b(t)

# back and forth substitution does not work for parameters with dependencies
term = value(a)
term2 = substitute(term, a=>b)
@test term2 isa Term{ModelingToolkit.Parameter{Real}}
@test isequal(term2, b)
term3 = substitute(term2, b=>a)
@test term3 isa Term{ModelingToolkit.Parameter{Real}}
@test isequal(term3, a)


text:
# Work in progress here

# not going to work proper yet
function propagate(b1Xl::BallTreeDensity, bDxb1::Array{Float64,1}, s::Array{Float64,1})
  pts1 = getPoints(b1Xl)
  pts = zeros(size(pts1))
  for i in 1:size(pts,2)
    ent = [s[1]*randn();s[2]*randn();s[3]*randn()]
    b1Tb = inv(SE2(bDxb1+ent))
    pt = vec(pts1[:,i])
    bTl = SE2([pt[1];pt[2];0.0])
    pts[:,i] = vec(se2vee(b1Tb*bTl))[1:2]
  end
  return kde!( pts, "lcv");
end

# not going to work proper yet
function update(bhatXl::BallTreeDensity, z::Array{Float64,1}, s::Array{Float64,1}; N::Int=75)

  bXl = p2cPtsKDE(z,s, N=N)
  # take the product between predicted and measured position
  dummy = kde!(rand(2,N),[1.0])
  pGM, = prodAppxMSGibbsS(dummy, [bhatXl, bXl], Union{}, Union{}, 5)

  # error("update -- Not ready")
  return kde!(pGM, "lcv")
end
function updatelin(bhatXl::BallTreeDensity, z::Array{Float64,1}, s::Array{Float64,1}; N::Int=75)

  bXl = resample(kde!((z')',s),N)
  # take the product between predicted and measured position
  dummy = kde!(rand(length(z),N),[1.0])
  pGM, = prodAppxMSGibbsS(dummy, [bhatXl, bXl], Union{}, Union{}, 5)

  # error("update -- Not ready")
  return kde!(pGM, "lcv")
end
function updatelin(bhatXl::BallTreeDensity, bXl::BallTreeDensity; N=75)
  dummy = kde!(rand(bhatXl.bt.dims,N),[1.0])
  pGM, = prodAppxMSGibbsS(dummy, [bhatXl, bXl], Union{}, Union{}, 5)
  # error("update -- Not ready")
  return kde!(pGM, "lcv")
end


text:
# Source/receiver geometry structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

export Geometry, compareGeometry, GeometryIC, GeometryOOC, get_nsrc, n_samples

abstract type Geometry{T} end

function rotateXY!(x::Array{<:Union{Float32,Float64},1}, y::Array{<:Union{Float32,Float64},1}, ϕ::Number)
    for i in 1:length(x)
        ynew = sin(ϕ)*x[i] + cos(ϕ)*y[i]
        x[i] = cos(ϕ)*x[i] - sin(ϕ)*y[i]
        y[i] = ynew
    end
end

const CoordT = Union{Array{T, 1}, Array{Array{T, 1}, 1}} where T<:Number
(::Type{CoordT})(x::Vector{Any}) = rebuild_maybe_jld(x)

# In-core geometry structure for seismic header information
mutable struct GeometryIC{T} <: Geometry{T}
    xloc::CoordT  # Array of receiver positions (fixed for all experiments)
    yloc::CoordT
    zloc::CoordT
    dt::Array{T,1}
    nt::Array{Integer,1}
    t::Array{T,1}
end

# Out-of-core geometry structure, contains look-up table instead of coordinates
mutable struct GeometryOOC{T} <: Geometry{T}
    container::Array{SegyIO.SeisCon,1}
    dt::Array{T,1}
    nt::Array{Integer,1}
    t::Array{T,1}
    nsamples::Array{Integer,1}
    key::String
    segy_depth_key::String

    # h5geo members
    h5geo::PyCall.PyObject
    crs::String
    lengthUnits::String
    temporalUnits::String
    model_origin_x::Float64
    model_origin_y::Float64
    model_orientation::Float64
end

######################## shapes easy access ################################
get_nsrc(g::GeometryIC) = length(g.xloc)
get_nsrc(g::GeometryOOC) = length(g.container)

n_samples(g::GeometryOOC, ::Info) = sum(g.nsamples)
n_samples(g::GeometryIC, info::Info) = sum([length(g.xloc[j])*g.nt[j] for j=1:info.nsrc])
n_samples(g::GeometryOOC, ::Integer) = sum(g.nsamples)
n_samples(g::GeometryIC, nsrc::Integer) = sum([length(g.xloc[j])*g.nt[j] for j=1:nsrc])

################################################ Constructors ####################################################################

"""
    Geometry
        xloc::Array{Array{T, 1},1}
        yloc::Array{Array{T, 1},1}
        zloc::Array{Array{T, 1},1}
        dt::Array{T,1}
        nt::Array{Integer,1}
        t::Array{T,1}

Geometry structure for seismic sources or receivers. Each field is a cell array, where individual cell entries\\
contain values or arrays with coordinates and sampling information for the corresponding shot position. The \\
first three entries are in meters and the last three entries in milliseconds.

Constructors
============

Only pass `dt` and `n` and automatically set `t`:

    Geometry(xloc, yloc, zloc; dt=[], nt=[])

Pass single array as coordinates/parameters for all `nsrc` experiments:

    Geometry(xloc, yloc, zloc, dt=[], nt=[], nsrc=1)

Create geometry structure for either source or receivers from a SegyIO.SeisBlock object.\\
`segy_depth_key` is the SegyIO keyword that contains the depth coordinate and `key` is \\
set to either `source` for source geometry or `receiver` for receiver geometry:

    Geometry(SeisBlock; key="source", segy_depth_key="")

Create geometry structure for from a SegyIO.SeisCon object (seismic data container):

    Geometry(SeisCon; key="source", segy_depth_key="")

Examples
========

(1) Set up receiver geometry for 2D experiment with 4 source locations and 80 fixed receivers:

    xrec = range(100,stop=900,length=80)
    yrec = range(0,stop=0,length=80)
    zrec = range(50,stop=50,length=80)
    dt = 4f0
    t = 1000f0

    rec_geometry = Geometry(xrec, yrec, zrec; dt=dt, t=t, nsrc=4)

(2) Set up corresponding source geometry (coordinates can be of type `linspace` or regular arrays):

    xsrc = [200,400,600,800]
    ysrc = [0,0,0,0]
    zsrc = [50,50,50,50]

    src_geometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=t, nsrc=4)

(3) Read source and receiver geometries from SEG-Y file:

    using SegyIO
    seis_block = segy_read("test_file.segy")
    rec_geometry = Geometry(seis_block; key="receiver", segy_depth_key="RecGroupElevation")
    src_geometry = Geometry(seis_block; key="source", segy_depth_key="SourceDepth")

Check the seis_block's header entries to findall out which keywords contain the depth coordinates.\\
The source depth keyword is either `SourceDepth` or `SourceSurfaceElevation`. The receiver depth \\
keyword is typically `RecGroupElevation`.

(4) Read source and receiver geometries from out-of-core SEG-Y files (for large data sets). Returns an out-of-core \\
geometry object `GeometryOOC` without the source/receiver coordinates, but a lookup table instead:

    using SegyIO
    seis_container = segy_scan("/path/to/data/directory","filenames",["GroupX","GroupY","RecGroupElevation","SourceDepth","dt"])
    rec_geometry = Geometry(seis_container; key="receiver", segy_depth_key="RecGroupElevation")
    src_geometry = Geometry(seis_container; key="source", segy_depth_key="SourceDepth")

"""
Geometry(xloc::CoordT, yloc::CoordT, zloc::CoordT, dt::Array{T,1}, nt::Array{Integer,1}, t::Array{T,1}) where T = GeometryIC{T}(xloc,yloc,zloc,dt,nt,t)

# Fallback constructors for non standard input types
function Geometry(xloc, yloc, zloc; dt=[], t=[], nsrc=nothing)
    if any(typeof(x) <: StepRangeLen for x=[xloc, yloc, zloc])
        args = [typeof(x) <: StepRangeLen ? collect(x) : x for x=[xloc, yloc, zloc]]
        isnothing(nsrc) && (return Geometry(args...; dt=dt, t=t))
        return Geometry(args...; dt=dt, t=t, nsrc=nsrc)
    end
    isnothing(nsrc) && (return Geometry(tof32(xloc), tof32(yloc), tof32(zloc); dt=dt, t=t))
    return Geometry(tof32(xloc), tof32(yloc), tof32(zloc); dt=dt, t=t, nsrc=nsrc)
end

# Constructor if nt is not passed
function Geometry(xloc::Array{Array{T, 1},1}, yloc::CoordT, zloc::Array{Array{T, 1},1};dt=[],t=[]) where T
    nsrc = length(xloc)
    # Check if single dt was passed
    dtCell = typeof(t) <: Real ? [T(dt) for j=1:nsrc] : T.(dt)
    # Check if single t was passed
    tCell = typeof(t) <: Real ? [T(t) for j=1:nsrc] : T.(t)

    # Calculate number of time steps
    ntCell = typeof(t) <: Real ? [floor(Int, t / dt) + 1 for j=1:nsrc] : floor.(Int, tCell ./ dtCell) .+ 1
    return GeometryIC{T}(xloc, yloc, zloc, dtCell, ntCell, tCell)
end

# Constructor if coordinates are not passed as a cell arrays
function Geometry(xloc::Array{T, 1}, yloc::CoordT, zloc::Array{T, 1}; dt=[], t=[], nsrc::Integer=1) where T
    xlocCell = [xloc for j=1:nsrc]
    ylocCell = [yloc for j=1:nsrc]
    zlocCell = [zloc for j=1:nsrc]
    dtCell = [T(dt) for j=1:nsrc]
    ntCell = [floor(Int, t/dt)+1 for j=1:nsrc]
    tCell = [T(t) for j=1:nsrc]
    return GeometryIC{T}(xlocCell, ylocCell, zlocCell, dtCell, ntCell, tCell)
end

################################################ Constructors from SEGY data  ####################################################

# Set up source geometry object from in-core data container
function Geometry(data::SegyIO.SeisBlock; key="source", segy_depth_key="")
    src = get_header(data,"FieldRecord")
    nsrc = length(unique(src))
    if key=="source"
        isempty(segy_depth_key) && (segy_depth_key="SourceSurfaceElevation")
        params = ["SourceX","SourceY",segy_depth_key]
        gt = Float32
    elseif key=="receiver"
        isempty(segy_depth_key) && (segy_depth_key="RecGroupElevation")
        params = ["GroupX","GroupY",segy_depth_key]
        gt = Array{Float32, 1}
    else
        throw("Specified keyword not supported")
    end
    xloc = Array{gt, 1}(undef, nsrc)
    yloc = Array{gt, 1}(undef, nsrc)
    zloc = Array{gt, 1}(undef, nsrc)
    dt = Array{Float32}(undef, nsrc)
    nt = Array{Integer}(undef, nsrc)
    t = Array{Float32}(undef, nsrc)

    xloc_full = get_header(data, params[1])
    yloc_full = get_header(data, params[2])
    zloc_full = get_header(data, params[3])
    dt_full = get_header(data, "dt")[1]
    nt_full = get_header(data, "ns")[1]

    for j=1:nsrc
        traces = findall(src .== unique(src)[j])
        if key=="source"    # assume same source location for all traces within one shot record
            xloc[j] = convert(gt, xloc_full[traces][1])
            yloc[j] = convert(gt,yloc_full[traces][1])
            zloc[j] = abs.(convert(gt,zloc_full[traces][1]))
        else
            xloc[j] = convert(gt, xloc_full[traces])
            yloc[j] = convert(gt, yloc_full[traces])
            zloc[j] = abs.(convert(gt, zloc_full[traces]))
        end
        dt[j] = Float32(dt_full/1f3)
        nt[j] = convert(Integer,nt_full)
        t[j] =  Float32((nt[j]-1)*dt[j])
    end

    if key == "source"
        xloc = convertToCell(xloc)
        yloc = convertToCell(yloc)
        zloc = convertToCell(zloc)
    end
    return GeometryIC{Float32}(xloc,yloc,zloc,dt,nt,t)
end

# Set up geometry summary from out-of-core data container
function Geometry(data::SegyIO.SeisCon; key="source", segy_depth_key="")

    if key=="source"
        isempty(segy_depth_key) && (segy_depth_key="SourceSurfaceElevation")
    elseif key=="receiver"
        isempty(segy_depth_key) && (segy_depth_key="RecGroupElevation")
    else
        throw("Specified keyword not supported")
    end

    # read either source or receiver geometry
    nsrc = length(data)
    container = Array{SegyIO.SeisCon}(undef, nsrc)
    dt = Array{Float32}(undef, nsrc)
    nt = Array{Integer}(undef, nsrc)
    t = Array{Float32}(undef, nsrc)
    nsamples = Array{Integer}(undef, nsrc)
    for j=1:nsrc
        container[j] = split(data,j)
        dt[j] = data.blocks[j].summary["dt"][1]/1f3
        nt[j] = data.ns
        t[j] = (nt[j]-1)*dt[j]
        key=="source" ? nsamples[j] = data.ns : nsamples[j] = Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4)*data.ns)
    end
    return  GeometryOOC{Float32}(container,dt,nt,t,nsamples,key,segy_depth_key,nothing,"","","",0.0,0.0,0.0)
end

# Set up geometry summary from out-of-core data container passed as cell array
function Geometry(data::Array{SegyIO.SeisCon,1}; key="source", segy_depth_key="")

    if key=="source"
        isempty(segy_depth_key) && (segy_depth_key="SourceSurfaceElevation")
    elseif key=="receiver"
        isempty(segy_depth_key) && (segy_depth_key="RecGroupElevation")
    else
        throw("Specified keyword not supported")
    end

    nsrc = length(data)
    container = Array{SegyIO.SeisCon}(undef, nsrc)
    dt = Array{Float32}(undef, nsrc); nt = Array{Integer}(undef, nsrc); t = Array{Float32}(undef, nsrc)
    nsamples = Array{Integer}(undef, nsrc)
    for j=1:nsrc
        container[j] = data[j]
        dt[j] = data[j].blocks[1].summary["dt"][1]/1f3
        nt[j] = data[j].ns
        t[j] = (nt[j]-1)*dt[j]
        key=="source" ? nsamples[j] = data[j].ns : nsamples[j] = Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4)*data[j].ns)
    end
    return  GeometryOOC{Float32}(container,dt,nt,t,nsamples,key,segy_depth_key,nothing,"","","",0.0,0.0,0.0)
end

# Load geometry from out-of-core Geometry container
function Geometry(geometry::GeometryOOC)
    nsrc = length(geometry.container)

    # read either source or receiver geometry
    if geometry.key=="source"
        params = ["SourceX","SourceY",geometry.segy_depth_key,"dt","ns"]
        gt = Float32
    elseif geometry.key=="receiver"
        params = ["GroupX","GroupY",geometry.segy_depth_key,"dt","ns"]
        gt = Array{Float32, 1}
    else
        throw("Specified keyword not supported")
    end
    xloc = Array{gt, 1}(undef, nsrc)
    yloc = Array{gt, 1}(undef, nsrc)
    zloc = Array{gt, 1}(undef, nsrc)
    dt = Array{Float32}(undef, nsrc); nt = Array{Integer}(undef, nsrc); t = Array{Float32}(undef, nsrc)

    for j=1:nsrc
        header = read_con_headers(geometry.container[j], params, 1)
        if geometry.key=="source"
            xloc[j] = convert(gt, get_header(header, params[1])[1])
            yloc[j] = convert(gt, get_header(header, params[2])[1])
            zloc[j] = abs.(convert(gt,get_header(header, params[3])[1]))
        else
            xloc[j] = convert(gt, get_header(header, params[1]))
            yloc[j] = convert(gt, get_header(header, params[2]))
            zloc[j] = abs.(convert(gt, get_header(header, params[3])))
        end
        if isnothing(geometry.h5geo)
            dt[j] = get_header(header, params[4])[1]/1f3
        else
            dt[j] = get_header(header, params[4])[1]
        end
        nt[j] = convert(Integer, get_header(header, params[5])[1])
        t[j] =  (nt[j]-1)*dt[j]
    end
    if geometry.key == "source"
        xloc = convertToCell(xloc)
        yloc = convertToCell(yloc)
        zloc = convertToCell(zloc)
    end

    if isnothing(geometry.h5geo)
        return GeometryIC(xloc,yloc,zloc,dt,nt,t)
    end

    dt = geometry.h5geo.sr.convertUnits(dt, geometry.temporalUnits, "ms")
    t = geometry.h5geo.sr.convertUnits(t, geometry.temporalUnits, "ms")
    if isnothing(dt) || length(dt) < 1 || isnothing(t) || length(t) < 1
        @error "Unable to convert temporal units"
    end

    for i in 1:length(xloc)
        xloc[i], yloc[i] = geometry.h5geo.sr.transformCoordFrom(
                            xloc[i], yloc[i], geometry.lengthUnits, geometry.crs)
        xloc[i] -= geometry.model_origin_x
        yloc[i] -= geometry.model_origin_y
        rotateXY!(xloc[i], yloc[i], geometry.model_orientation)
        xloc[i] += geometry.model_origin_x
        yloc[i] += geometry.model_origin_y
        xloc[i] = geometry.h5geo.sr.convertUnits(xloc[i], geometry.lengthUnits, "m")
        yloc[i] = geometry.h5geo.sr.convertUnits(yloc[i], geometry.lengthUnits, "m")
    end

    return GeometryIC(xloc,yloc,zloc,dt,nt,t)
end

Geometry(geometry::GeometryIC) = geometry
Geometry(::Nothing) = nothing

###########################################################################################################################################

# Subsample in-core geometry structure
function subsample(geometry::GeometryIC,srcnum)
    if length(srcnum)==1
        srcnum = srcnum[1]
        geometry = Geometry(geometry.xloc[srcnum], geometry.yloc[srcnum], geometry.zloc[srcnum];
                            dt=geometry.dt[srcnum],t=geometry.t[srcnum],nsrc=1)
    else
        geometry = Geometry(geometry.xloc[srcnum], geometry.yloc[srcnum], geometry.zloc[srcnum],
                            geometry.dt[srcnum], geometry.nt[srcnum], geometry.t[srcnum])
    end
    return geometry
end

# Subsample out-of-core geometry structure
subsample(geometry::GeometryOOC, srcnum) = Geometry(geometry.container[srcnum]; key=geometry.key, segy_depth_key=geometry.segy_depth_key)

# Compare if geometries match
function compareGeometry(geometry_A::Geometry, geometry_B::Geometry)
    if isequal(geometry_A.xloc, geometry_B.xloc) && isequal(geometry_A.yloc, geometry_B.yloc) && isequal(geometry_A.zloc, geometry_B.zloc) &&
    isequal(geometry_A.dt, geometry_B.dt) && isequal(geometry_A.nt, geometry_B.nt)
        return true
    else
        return false
    end
end

isequal(geometry_A::Geometry, geometry_B::Geometry) = compareGeometry(geometry_A, geometry_B)

function compareGeometry(geometry_A::GeometryOOC, geometry_B::GeometryOOC)
    check = true
    for j=1:length(geometry_A.container)
        if ~isequal(geometry_A.container[j].blocks[1].summary["GroupX"], geometry_B.container[j].blocks[1].summary["GroupX"]) ||
        ~isequal(geometry_A.container[j].blocks[1].summary["GroupY"], geometry_B.container[j].blocks[1].summary["GroupY"]) ||
        ~isequal(geometry_A.container[j].blocks[1].summary["SourceX"], geometry_B.container[j].blocks[1].summary["SourceX"]) ||
        ~isequal(geometry_A.container[j].blocks[1].summary["SourceY"], geometry_B.container[j].blocks[1].summary["SourceY"]) ||
        ~isequal(geometry_A.container[j].blocks[1].summary["dt"], geometry_B.container[j].blocks[1].summary["dt"])
            check = false
        end
    end
    return check
end

isequal(geometry_A::GeometryOOC, geometry_B::GeometryOOC) = compareGeometry(geometry_A, geometry_B)

compareGeometry(geometry_A::GeometryOOC, geometry_B::Geometry) = true
compareGeometry(geometry_A::Geometry, geometry_B::GeometryOOC) = true

for G in [GeometryOOC, GeometryIC]
    @eval function push!(G1::$G, G2::$G)
        for k in fieldnames($G)
            pushfield!(getfield(G1, k), getfield(G2, k))
        end
    end
end

pushfield!(a::Array, b::Array) = append!(a, b)
pushfield!(a, b) = nothing


text:
@testset "BLP inversion: $it" for it in (SimpleFixedPointIteration(), SQUAREMIteration())
    # This random seed gives enough of an outside option so the contraction can swiftly converge
    Random.seed!(123)

    J = 50
    I = 1000

    sigma = [2.0 0.0; 4.0 1.0]
    K2 = size(sigma, 1)

    x2 = exp.(randn((J, K2)))
    nodes = randn((I, K2))
    weights = fill(1/I, I)

    delta = exp.(randn(J))
    mu = x2 * sigma * nodes'

    shares = choice_probabilities(delta .+ mu) * weights

    market = Market(x2, shares, weights, nodes)

    result = compute_delta(market, it, sigma)

    @test result.status == ITERATION_CONVERGED
    @test result.delta ≈ delta
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function create_cache(::Type{IndicatorHennemannGassner}, equations::AbstractEquations{1}, basis::LobattoLegendreBasis)

  alpha = Vector{real(basis)}()
  alpha_tmp = similar(alpha)

  A = Array{real(basis), ndims(equations)}
  indicator_threaded  = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]
  modal_threaded      = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]

  return (; alpha, alpha_tmp, indicator_threaded, modal_threaded)
end

# this method is used when the indicator is constructed as for AMR
function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, equations::AbstractEquations{1}, dg::DGSEM, cache)
  create_cache(typ, equations, dg.basis)
end

function (indicator_hg::IndicatorHennemannGassner)(u::AbstractArray{<:Any,3},
                                                   equations, dg::DGSEM, cache;
                                                   kwargs...)
  @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg
  @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache
  # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR?
  #       Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)`
  #       or just `resize!` whenever we call the relevant methods as we do now?
  resize!(alpha, nelements(dg, cache))
  if alpha_smooth
    resize!(alpha_tmp, nelements(dg, cache))
  end

  # magic parameters
  threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25)
  parameter_s = log((1 - 0.0001)/0.0001)

  @threaded for element in eachelement(dg, cache)
    indicator  = indicator_threaded[Threads.threadid()]
    modal      = modal_threaded[Threads.threadid()]

    # Calculate indicator variables at Gauss-Lobatto nodes
    for i in eachnode(dg)
      u_local = get_node_vars(u, equations, dg, i, element)
      indicator[i] = indicator_hg.variable(u_local, equations)
    end

    # Convert to modal representation
    multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, indicator)

    # Calculate total energies for all modes, without highest, without two highest
    total_energy = zero(eltype(modal))
    for i in 1:nnodes(dg)
      total_energy += modal[i]^2
    end
    total_energy_clip1 = zero(eltype(modal))
    for i in 1:(nnodes(dg)-1)
      total_energy_clip1 += modal[i]^2
    end
    total_energy_clip2 = zero(eltype(modal))
    for i in 1:(nnodes(dg)-2)
      total_energy_clip2 += modal[i]^2
    end

    # Calculate energy in lower modes
    energy = max((total_energy - total_energy_clip1) / total_energy,
                 (total_energy_clip1 - total_energy_clip2) / total_energy_clip1)

    alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold)))

    # Take care of the case close to pure DG
    if alpha_element < alpha_min
      alpha_element = zero(alpha_element)
    end

    # Take care of the case close to pure FV
    if alpha_element > 1 - alpha_min
      alpha_element = one(alpha_element)
    end

    # Clip the maximum amount of FV allowed
    alpha[element] = min(alpha_max, alpha_element)
  end

  if (alpha_smooth)
    # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha
    # Copy alpha values such that smoothing is indpedenent of the element access order
    alpha_tmp .= alpha

    # Loop over interfaces
    for interface in eachinterface(dg, cache)
      # Get neighboring element ids
      left  = cache.interfaces.neighbor_ids[1, interface]
      right = cache.interfaces.neighbor_ids[2, interface]

      # Apply smoothing
      alpha[left]  = max(alpha_tmp[left],  0.5 * alpha_tmp[right], alpha[left])
      alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left],  alpha[right])
    end
  end

  return alpha
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function create_cache(::Type{IndicatorLöhner}, equations::AbstractEquations{1}, basis::LobattoLegendreBasis)

  alpha = Vector{real(basis)}()

  A = Array{real(basis), ndims(equations)}
  indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]

  return (; alpha, indicator_threaded)
end

# this method is used when the indicator is constructed as for AMR
function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEquations{1}, dg::DGSEM, cache)
  create_cache(typ, equations, dg.basis)
end

function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any,3},
                                   equations, dg::DGSEM, cache;
                                   kwargs...)
  @assert nnodes(dg) >= 3 "IndicatorLöhner only works for nnodes >= 3 (polydeg > 1)"
  @unpack alpha, indicator_threaded = löhner.cache
  resize!(alpha, nelements(dg, cache))

  @threaded for element in eachelement(dg, cache)
    indicator = indicator_threaded[Threads.threadid()]

    # Calculate indicator variables at Gauss-Lobatto nodes
    for i in eachnode(dg)
      u_local = get_node_vars(u, equations, dg, i, element)
      indicator[i] = löhner.variable(u_local, equations)
    end

    estimate = zero(real(dg))
    for i in 2:nnodes(dg)-1
      # x direction
      u0 = indicator[i  ]
      up = indicator[i+1]
      um = indicator[i-1]
      estimate = max(estimate, local_löhner_estimate(um, u0, up, löhner))
    end

    # use the maximum as DG element indicator
    alpha[element] = estimate
  end

  return alpha
end

# this method is used when the indicator is constructed as for shock-capturing volume integrals
function create_cache(::Type{IndicatorMax}, equations::AbstractEquations{1}, basis::LobattoLegendreBasis)

  alpha = Vector{real(basis)}()

  A = Array{real(basis), ndims(equations)}
  indicator_threaded = [A(undef, nnodes(basis)) for _ in 1:Threads.nthreads()]

  return (; alpha, indicator_threaded)
end

# this method is used when the indicator is constructed as for AMR
function create_cache(typ::Type{IndicatorMax}, mesh, equations::AbstractEquations{1}, dg::DGSEM, cache)
  cache = create_cache(typ, equations, dg.basis)
end

function (indicator_max::IndicatorMax)(u::AbstractArray{<:Any,3},
                                       equations, dg::DGSEM, cache;
                                       kwargs...)
  @unpack alpha, indicator_threaded = indicator_max.cache
  resize!(alpha, nelements(dg, cache))

  @threaded for element in eachelement(dg, cache)
    indicator = indicator_threaded[Threads.threadid()]

    # Calculate indicator variables at Gauss-Lobatto nodes
    for i in eachnode(dg)
      u_local = get_node_vars(u, equations, dg, i, element)
      indicator[i] = indicator_max.variable(u_local, equations)
    end

    alpha[element] = maximum(indicator)
  end

  return alpha
end

end # @muladd


text:
using Documenter, MLMetrics

# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

makedocs(
    modules = [MLMetrics],
    clean = false,
    format = :html,
    assets = [
        joinpath("assets", "favicon.ico"),
        joinpath("assets", "style.css"),
    ],
    sitename = "MLMetrics.jl",
    authors = "Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        "gettingstarted.md",
        hide("classification.md", Any[
            joinpath.("counts", readdir(joinpath(@__DIR__, "src", "counts")))...,
            joinpath.("fractions", readdir(joinpath(@__DIR__, "src", "fractions")))...,
        ]),
        hide("Indices" => "indices.md"),
        "LICENSE.md",
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/JuliaML/MLMetrics.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
"""
    Network{T<:Integer}

Network of HTM [`Regions`](@ref) connected with feedforward (proximal) and/or contextual (distal) connections.
"""
struct Network{T<:Integer}
  regions::Vector{Region}
  "`region_α₋` is the cell activity of each region at the previous time step"
  region_α₋::Vector{<:HierarchicalTemporalMemory.CellActivity}
  """
  `circuit_proximal` is the proximal connection graph between regions.
  There are `length(regions) + 2` vertices, because at the end of the list an extra vertex is added for the input and the output.
  """
  circuit_proximal::SimpleDiGraph{T}
  circuit_distal::SimpleDiGraph{T}
  in_vertex::Int
  out_vertex::Int
end

"""
    Network(regions; connect_forward, connect_context)

Create a network of HTM [`Regions`](@ref).
The given regions are connected with feedforward (proximal) and/or contextual (distal) connections according to the given adjacency matrices.
One region serves as the input point and another one as the output, designated by their index on the `regions` vector.

## Inputs

- `regions::Vector{Regions}`: HTM [`Regions`](@ref) that will form the network
- `connect_forward: [N+2 × N+2]`: adjacency matrix for feedforward connections between regions and input/output nodes
- `connect_context: [N+2 × N+2]`: adjacency matrix for contextual connections between regions and input/output nodes

TODO: remove I/O from connect_* and move them to separate inputs:
TODO: MIMO, extend I/O to vectors of regions

- `in_forward`: index of the region that will receive external feedforward input
- `in_context`: index of the region that will receive external contextual input
- `out`: index of the region that will emit network output

## Adjacency matrix

The adjacency matrix has size [N+2 × N+2] because the 2 last nodes represent the input and the output connections.
Since the connections are directional, the adjacency matrix has an order: pre/post.

- The first N rows/columns refer to the connections between regions.
- Row N+1 is the input connection
- Column N+2 is the output connection

Example:

```
# vertices:
# 1 2 3 I O
c= [
  0 0 0 0 0  # 1
  1 0 0 0 0  # 2
  0 1 0 0 1  # 3
  1 0 0 0 0  # In
  0 0 0 0 0  # Out
]
# Connections:
c[2,1]  # 2 -> 1
c[3,2]  # 3 -> 2
c[4,1]  # in -> 1
c[3,5]  # 3 -> out
```
"""
Network(regions; connect_forward, connect_context=zeros(length(regions).*(1,1)))= begin
  circuit_proximal= SimpleDiGraph(length(regions)+2)
  circuit_distal= SimpleDiGraph(length(regions)+2)
  add_edges!(circuit_proximal, connect_forward)
  add_edges!(circuit_distal, connect_context)
  # Initialize the state of each region to `falses`
  Network(regions, map(r-> falses(length(r)), regions),
      circuit_proximal, circuit_distal,
      length(regions)+1, length(regions)+2)
end

reset!(n::Network)= begin
  foreach(reset!, n.regions)
  n.region_α₋ .= map(r-> falses(length(r)), n.regions)
  n
end

"""
    networkPlot(network, plotFunction, extraArgs...)

Uses `plotFunction` to plot a network graph. Assumes that `plotFunction == GraphMakie.graphplot`, but might work with others.
`extraArgs` are passed as-is to the plotting function.
"""
networkPlot(network, plotFunction, extraArgs...)= begin
  allconnections= network.circuit_proximal ∪ network.circuit_distal
  # TODO: find the edges that come from the distal synapses!
  edgecolors= [
    e ∈ network.circuit_proximal|>edges ? :black : :cyan for e in allconnections|>edges
  ]
  plotFunction(allconnections,
      nlabels= [string.(1:length(network.regions));"in";"out"],
      node_color= [map(_->:black,1:length(network.regions)); :grey; :grey],  # Can't seem to get this to work
      node_marker= [map(_->:circle,1:length(network.regions)); :diamond; :diamond],  # Can't seem to get this to work
      edge_color= edgecolors, arrow_size= 12, node_size= 10,
      extraArgs...
  )
end

"`add_edges!(g,adjacency)` inserts the connections from an adjacency matrix into the graph"
add_edges!(g,adjacency)= map(findall(adjacency .> 0)) do edge
		add_edge!(g, edge.I...)
	end

"Evaluate a \"vertex map\", applying the Network's state to each region"
(n::Network)(proximal, distal=falses(0))= activateNetwork(n,activate,proximal,distal)
step!(n::Network, proximal, distal=falses(0))= activateNetwork(n,step!,proximal,distal)

getindexOverflow(a,overflow,i)= getindexOverflow.(Ref(a),Ref(overflow),i)
getindexOverflow(a,overflow,i::Integer)= i .<= length(a) ? a[i] : overflow

"`activateNetwork(n,f,proximal,distal)` calculates each region's input and applies function `f`. It's the common denominator of network activation functions."
activateNetwork(n,f,proximal,distal)= begin
  # For each region, calculate its activation based on the previous activation of every neighbor
  n.region_α₋.= map(n.regions|> enumerate) do (i,r)
    proximal_input= getindexOverflow(n.region_α₋, proximal, inneighbors(n.circuit_proximal,i))
    distal_input= getindexOverflow(n.region_α₋, distal, inneighbors(n.circuit_distal,i))
    f(r,proximal_input,distal_input).active
  end
  # Emit output activation
  getindexOverflow(n.region_α₋, proximal, inneighbors(n.circuit_proximal,n.out_vertex))|> gateCombine
end

"""
    propagation_delay(n::Network)

The shortest time after repeated stimulation with `a` that `n(a)` will be related to `a` (the minimum number of repeated stimulations).
The network will always produce some activity, but until that point it will be unrelated to the input `a`.
The signal traverses 1 Region inside the Network with each stimulation, therefore the number of regions between the input and the output
of the network is the minimum time necessary before the output has a relationship with the input.

The network might not be stable at this point and continued stimulation with `a` might allow longer paths to the output to contribute.
"""
propagation_delay(n::Network)= length(a_star(n.circuit_proximal, n.in_vertex, n.out_vertex)) - 1


text:
@testset "Testing Corner Cases" begin
    Local, Remote = CachedArrays.Local, CachedArrays.Remote
    manager = CachedArrays.CacheManager(
        CachedArrays.AlignedAllocator(),
        CachedArrays.AlignedAllocator();
        localsize = 100_000_000,
        remotesize = 500_000_000,
        minallocation = 12, # 4096 Bytes
    )

    # Case 1: Requesting movement of a Region that is in the free buffer.
    let
        _a = CachedArrays.CachedArray{Float32}(undef, manager, 10_000)
        ### Unsafe operation.
        # Run finalizers on the `Region` backing `_a`.
        # This will put it on the manager's freebuffer.
        # Don't use `_a` after this operation.
        finalize(_a.object)

        # When we manually invoke an eviction, the manager will discover during allocation
        # that this object is queued to be freed.
        #
        # As such, it will not even move the array and simple return.
        CachedArrays.@spinlock CachedArrays.alloc_lock(manager) begin
            CachedArrays.evict!(_a, manager.policy, manager)
        end

        @test length(CachedArrays.visible_ids(manager, Local)) == 0
        @test length(CachedArrays.visible_ids(manager, Remote)) == 0

        @test CachedArrays.check(manager.remote_heap)
        @test CachedArrays.check(manager.local_heap)
        @test CachedArrays.check(manager)
    end

    # Case 2: Try the same movement corner case, but going the other direction.
    let
        _a = CachedArrays.CachedArray{Float32}(undef, manager, 10_000)

        # Move this block into the remote heap.
        # Make sure that it actually resides in the remote heap and that it doesn't
        # have any siblings in the local heap for some reason.
        CachedArrays.evict!(_a)
        @test CachedArrays.getpool(CachedArrays.metadata(_a)) == CachedArrays.Remote
        @test CachedArrays.getsibling(CachedArrays.metadata(_a)) === nothing

        ### Unsafe operation.
        # Run finalizers on the `Region` backing `_a`.
        # This will put it on the manager's freebuffer.
        # Don't use `_a` after this operation.
        finalize(_a.object)

        # When we manually invoke an eviction, the manager will discover during allocation
        # that this object is queued to be freed.
        #
        # As such, it will not even move the array and simple return.
        CachedArrays.prefetch!(_a)

        @test length(CachedArrays.visible_ids(manager, Local)) == 0
        @test length(CachedArrays.visible_ids(manager, Remote)) == 0

        @test CachedArrays.check(manager.remote_heap)
        @test CachedArrays.check(manager.local_heap)
        @test CachedArrays.check(manager)
    end
end


text:
@with_kw struct NoisyIFParameter <: AbstractIFParameter
    σ::SNNFloat = 0
end

@with_kw mutable struct NoisyIF <: AbstractIF
    param::NoisyIFParameter = NoisyIFParameter()
    randncache::Vector{SNNFloat} = randn(N)
end

function integrate!(p::NoisyIF, param::NoisyIFParameter, dt::SNNFloat)
    randn!(randncache)
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i] + σ / √dt * randncache[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
module BSplines

import LinearAlgebra

using Base: @propagate_inbounds, has_offset_axes
using LinearAlgebra: I, diagind, ldiv!, lmul!, rdiv!, rmul!
using OffsetArrays: OffsetArray

export
    # Types
    AllDerivatives,
    BSpline,
    BSplineBasis,
    Derivative,
    Spline,

    # Functions
    approximate,
    averagebasis,
    basis,
    basismatrix,
    basismatrix!,
    breakpoints,
    bsplines,
    bsplines!,
    coeffs,
    interpolate,
    intervalindex,
    intervalindices,
    knotaverages,
    knotaverages!,
    knots,
    order,
    splinevalue,
    support

include("vectortypes.jl")
include("derivatives.jl")
include("bsplinebasis.jl")
include("spline.jl")
include("plotting.jl")

bspline_returntype(basis::BSplineBasis, x::Number) = bspline_returntype(basis, typeof(x))
bspline_returntype(basis::BSplineBasis, types::Type...) =
    bspline_returntype(eltype(breakpoints(basis)), types...)

bspline_returntype(spline::Spline, x::Number) = bspline_returntype(spline, typeof(x))
bspline_returntype(spline::Spline, xtype::Type) =
    bspline_returntype(basis(spline), eltype(coeffs(spline)), xtype)

function bspline_returntype(knottype::Type, types::Type...)
    T = promote_type(knottype, types...)
    typeof(oneunit(T)/oneunit(knottype))
end

"""
    approximate(f, basis::BSplineBasis; indices=eachindex(basis)) -> Spline

Approximate the function `f` in the B-spline basis `basis`. If `indices` is supplied, only
the basis functions at the given indices of `basis` are used.

The approximation is calculated by interpolating samples of `f` at the
[`knotaverages`](@ref) of the basis.

See also: [`interpolate`](@ref)

# Examples

```jldoctest
julia> basis = BSplineBasis(6, 0:0.25:1);

julia> spl = approximate(sin, basis, indices=2:length(basis))
$(Spline{BSplineBasis{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}, Vector{Float64}}):
 basis: 9-element $(BSplineBasis{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}):
  order: 6
  breakpoints: 0.0:0.25:1.0
 coeffs: [0.0, 0.05, 0.15, 0.298438, 0.486979, 0.651299, 0.755166, 0.814456, 0.841471]

julia> spl(π/4) ≈ sqrt(1/2)
true
```
"""
function approximate(f, basis::BSplineBasis; indices::Union{AbstractUnitRange,Colon}=Colon())
    xvalues = knotaverages(basis, indices=indices)
    interpolate(basis, xvalues, f.(xvalues); indices=indices)
end

"""
    interpolate(basis::BSplineBasis, xvalues, yvalues; indices=eachindex(basis)) -> Spline

Interpolate the data given by `xvalues` and `yvalues` in the B-spline basis `basis`. If
`indices` is supplied, only the basis functions at the given indices of `basis` are used.

The spline interpolation is calculated by creating the matrix
`B = [basis[i](x) for x=xvalues, i=indices]` and then calculating `B\\yvalues`.

See also: [`approximate`](@ref)

# Examples

```jldoctest
julia> basis = BSplineBasis(5, 1:10);

julia> xs = range(1, stop=10, length=length(basis)); ys = log.(xs);

julia> spl = interpolate(basis, xs, ys)
$(Spline{BSplineBasis{UnitRange{Int64}}, Vector{Float64}}):
 basis: 13-element BSplineBasis{UnitRange{Int64}}:
  order: 5
  breakpoints: 1:10
 coeffs: [0.0, 0.248019, 0.596872, 0.946671, 1.26894, 1.51405, 1.71149, 1.87666, 2.01856, 2.14292, 2.22592, 2.27758, 2.30259]

julia> spl(float(ℯ))
0.9999766059171411
```
"""
function interpolate(basis::BSplineBasis, xvalues::AbstractVector, yvalues::AbstractVector;
                     indices::Union{AbstractUnitRange,Colon}=Colon())
    if !checkbounds(Bool, basis, indices)
        throw(ArgumentError("indices must be valid indices for basis."))
    end
    if length(xvalues) != length(yvalues)
        throw(DimensionMismatch("lengths of x and y values do not match."))
    end
    bmatrix = basismatrix(basis, xvalues, indices=indices)
    _interpolate(basis, bmatrix, yvalues, indices)
end

function _interpolate(basis, basismatrix, yvalues, indices::AbstractUnitRange)
    indcoeffs = basismatrix \ yvalues
    coeffs = zeros(eltype(indcoeffs), length(basis))
    coeffs[indices] = indcoeffs
    Spline(basis, coeffs)
end

_interpolate(basis, basismatrix, yvalues, ::Colon) = Spline(basis, basismatrix\yvalues)

"""
    knotaverages(basis::BSplineBasis; indices=eachindex(basis))

Return the knot averages `τ[i] = mean(knots[i+1:i+order-1])` for `i ∈ indices` and the
knots and order of `basis`. The knot averages are recommended in [^deBoor1978] (p. 214) as
data points for interpolation.

See also: [`knotaverages!`](@ref)

[^deBoor1978]:
    Carl de Boor, *A Practical Guide to Splines*, New York, N.Y.: Springer, 1978.

# Examples

```jldoctest
julia> knotaverages(BSplineBasis(3, 0:5))
7-element $(Vector{Float64}):
 0.0
 0.5
 1.5
 2.5
 3.5
 4.5
 5.0

julia> knotaverages(BSplineBasis(4, [1, 3//2, 5//2, 4]), indices=2:6)
5-element $(Vector{Rational{Int64}}):
 7//6
 5//3
 8//3
 7//2
 4//1
```
"""
function knotaverages(basis::BSplineBasis; indices=Colon())
    if !checkbounds(Bool, basis, indices)
        throw(ArgumentError("invalid indices for basis: $indices"))
    end
    T = bspline_returntype(basis, Bool)
    dest = Vector{T}(undef, length_indices(basis, indices))
    _knotaverages!(dest, basis, range_indices(basis, indices))
end

"""
    knotaverages!(dest, basis::BSplineBasis; indices=eachindex(basis))

Calculate the knot averages `τ[i] = mean(knots[i+1:i+order-1])` for `i ∈ indices` and the
knots and order of `basis` and store the result in `dest`. The knot averages are recommended
in [^deBoor1978] (p. 214) as data points for interpolation.

See also: [`knotaverages`](@ref)

[^deBoor1978]:
    Carl de Boor, *A Practical Guide to Splines*, New York, N.Y.: Springer, 1978.

# Examples

```jldoctest
julia> dest = Vector{Float64}(undef, 7);

julia> knotaverages!(dest, BSplineBasis(3, 0:5))
7-element $(Vector{Float64}):
 0.0
 0.5
 1.5
 2.5
 3.5
 4.5
 5.0

julia> dest = Vector{Rational{Int}}(undef, 5);

julia> knotaverages!(dest, BSplineBasis(3, 0:5), indices=2:6)
5-element $(Vector{Rational{Int64}}):
 1//2
 3//2
 5//2
 7//2
 9//2
```
"""
function knotaverages!(dest, basis::BSplineBasis; indices=Colon())
    has_offset_axes(dest) && throw(ArgumentError("destination vector may not have offset axes."))
    if !checkbounds(Bool, basis, indices)
        throw(ArgumentError("invalid indices for basis: $indices"))
    end
    len = length_indices(basis, indices)
    if length(dest) != len
        throw(DimensionMismatch("length of destination, $(length(dest)), " *
                                "does not match number of indices, $len."))
    end
    _knotaverages!(dest, basis, range_indices(basis, indices))
end

@propagate_inbounds function _knotaverages!(dest, basis, indices::AbstractUnitRange)
    km1 = order(basis) - 1
    moving_average!(dest, view(knots(basis), first(indices)+1:last(indices)+km1), km1)
end

@propagate_inbounds function moving_average!(dest, xs, n)
    n⁻¹ = inv(convert(eltype(dest), n))
    moving_sum = sum(xs[j] for j=1:n)
    dest[1] = n⁻¹ * moving_sum
    for i = firstindex(dest)+1:lastindex(dest)
        moving_sum -= xs[i-1]
        moving_sum += xs[i-1+n]
        dest[i] = n⁻¹ * moving_sum
    end
    dest
end

"""
    averagebasis(order, datapoints) -> BSplineBasis

Returns a B-spline basis with the specified `order` that is well-suited for interpolation on
the given `datapoints`. The `datapoints` vector is assumed to be sorted.

The calculated breakpoints are described in [^deBoor1978], p. 219, as a “reasonable
alternative” to the optimal breakpoint sequence since they are “often very close to the
optimum” and are computationally inexpensive.

[^deBoor1978]:
    Carl de Boor, *A Practical Guide to Splines*, New York, N.Y.: Springer, 1978.

# Examples

```jldoctest
julia> averagebasis(5, 0:10)
11-element $(BSplineBasis{Vector{Float64}}):
 order: 5
 breakpoints: [0.0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 10.0]
```
"""
function averagebasis(order::Integer, datapoints::AbstractVector{<:Real})
    if order > length(datapoints)
        throw(ArgumentError("order cannot not be greater than number of datapoints."))
    end
    BSplineBasis(order, averagebasis_breakpoints(order, datapoints))
end

function averagebasis_breakpoints(k, datapoints)
    T = bspline_returntype(Int, eltype(datapoints))
    bps = Vector{T}(undef, length(datapoints)-k+2)
    bps[1] = first(datapoints)
    bps[end] = last(datapoints)
    if length(bps) > 2
        moving_average!(@view(bps[2:end-1]), @view(datapoints[2:end-1]), k-1)
    end
    bps
end

end # module BSplines


text:
using Colors, FixedPointNumbers, Test, InteractiveUtils

@testset "Utilities" begin

    @test hex(RGB(1,0.5,0)) == "FF8000"
    @test hex(RGBA(1,0.5,0,0.25)) == "40FF8000"

    # test utility function weighted_color_mean
    parametric2 = [GrayA,AGray32,AGray]
    parametric3 = ColorTypes.parametric3
    parametric4A = setdiff(subtypes(ColorAlpha),[GrayA])
    parametricA4 = setdiff(subtypes(AlphaColor),[AGray32,AGray])
    colortypes = vcat(parametric3,parametric4A,parametricA4)
    colorElementTypes = [Float16,Float32,Float64,BigFloat,N0f8,N6f10,N4f12,N2f14,N0f16]

    for T in colorElementTypes
        c1 = Gray{T}(1)
        c2 = Gray{T}(0)
        @test weighted_color_mean(0.5,c1,c2) == Gray{T}(0.5)
    end

    for C in parametric2
        for T in colorElementTypes
            C == AGray32 && (T=N0f8)
        c1 = C(T(0),T(1))
            c2 = C(T(1),T(0))
        @test weighted_color_mean(0.5,c1,c2) == C(T(1)-T(0.5),T(0.5))
        end
    end

    for C in colortypes
        for T in (Float16,Float32,Float64,BigFloat)
        if C<:Color
                c1 = C(T(1),T(1),T(0))
                c2 = C(T(0),T(1),T(1))
            @test weighted_color_mean(0.5,c1,c2) == C(T(0.5),T(0.5)+T(1)-T(0.5),T(1)-T(0.5))
        else
                C == ARGB32 && (T=N0f8)
                c1 = C(T(1),T(1),T(0),T(1))
                c2 = C(T(0),T(1),T(1),T(0))
            @test weighted_color_mean(0.5,c1,c2) == C(T(0.5),T(0.5)+T(1)-T(0.5),T(1)-T(0.5),T(0.5))
        end
        end
    end

    # test utility function range
    # range uses weighted_color_mean which is extensively tested.
    # Therefore it suffices to test the function using gray colors.
    for T in colorElementTypes
        c1 = Gray(T(1))
        c2 = Gray(T(0))
        linc1c2 = range(c1,stop=c2,length=43)
        @test length(linc1c2) == 43
        @test linc1c2[1] == c1
        @test linc1c2[end] == c2
        @test linc1c2[22] == Gray(T(0.5))
        @test typeof(linc1c2) == Array{Gray{T},1}
        if VERSION >= v"1.1"
            @test range(c1,c2,length=43) == range(c1,stop=c2,length=43)
        end
    end
end


text:
Random.seed!(0)

# Dynamics
dt = 0.001
J = [1; 5; 5]
function dynamics(x,u)
    xp = x + dt*u./J
    xp[1] += dt*(J[2]-J[3])/J[1]*x[2]*x[3]
    xp[2] += dt*(J[3]-J[1])/J[2]*x[3]*x[1]
    xp[3] += dt*(J[1]-J[2])/J[3]*x[1]*x[2]
    return xp
end

W = 0.001*Matrix{Float64}(I,3,3)

dmodel = NonlinearDynamicsModel(dynamics,W);

# Measurement
c=10
function observe(x,u)
    y = x
    y = min.(y,c)
    y = max.(y,-c)
    return y
end

V = 0.3*Matrix{Float64}(I,3,3)

omodel = NonlinearObservationModel(observe,V)

# Build EKF
ekf = ExtendedKalmanFilter(dmodel,omodel)

# Simulation
times = 0:dt:0.1
action_sequence = [[0.0,0.0,0.0] for t in times]

b0 = GaussianBelief([10.0,0.0,0.0], Matrix{Float64}(I,3,3))

sim_states, sim_measurements = simulation(ekf,b0,action_sequence)

# Filter
filtered_beliefs = run_filter(ekf, b0, action_sequence, sim_measurements)
μ, Σ = unpack(filtered_beliefs)

# Tests
@test ekf.d isa NonlinearDynamicsModel
@test ekf.o isa NonlinearObservationModel
@test length(sim_measurements) == length(action_sequence)

first_update = update(ekf, b0, action_sequence[1], sim_measurements[1])
first_predict = predict(ekf, b0, action_sequence[1])
first_predict_measure = measure(ekf, first_predict, sim_measurements[1], u = action_sequence[1])
array_isapprox(first_update.μ, first_predict_measure.μ)
array_isapprox(first_update.Σ, first_predict_measure.Σ)

@test length(filtered_beliefs[end].μ) == length(b0.μ)
@test size(filtered_beliefs[end].Σ) == size(b0.Σ)
@test size(μ) == (length(filtered_beliefs), length(filtered_beliefs[1].μ))
@test size(Σ) == (length(filtered_beliefs), size(filtered_beliefs[1].Σ)...)

# change with algorithm / rng changes

array_isapprox(sim_states[end], [10.7664, 0.213008, -0.463593];atol=0.001)
array_isapprox(sim_measurements[end], [10.9978, 0.491584, -0.33985];atol=0.001)
array_isapprox(filtered_beliefs[end].μ, [10.051, -0.0285158, -0.575077];
                atol=0.001)
array_isapprox(filtered_beliefs[end].Σ, [0.329097 -0.00281178 0.000218104;
                                        -0.00281178 0.0168616 -1.87166e-6;
                                        0.000218104 -1.87166e-6 0.0168371]; atol=0.001)


text:
import Base: delete!

# `Attributes` is a collection of named attributes

"Return all attribute names"
function listAtt(ncid,varid)
    natts = nc_inq_varnatts(ncid,varid)
    names = Vector{String}(undef,natts)

    for attnum = 0:natts-1
        names[attnum+1] = nc_inq_attname(ncid,varid,attnum)
    end

    return names
end

function Base.get(a::BaseAttributes, name::SymbolOrString,default)
    if haskey(a,name)
        return a[name]
    else
        return default
    end
end

"""
    getindex(a::Attributes,name::SymbolOrString)

Return the value of the attribute called `name` from the
attribute list `a`. Generally the attributes are loaded by
indexing, for example:

```julia
ds = NCDataset("file.nc")
title = ds.attrib["title"]
```
"""
function Base.getindex(a::Attributes,name::SymbolOrString)
    return nc_get_att(a.ds.ncid,a.varid,name)
end

"""
    Base.setindex!(a::Attributes,data,name::SymbolOrString)

Set the attribute called `name` to the value `data` in the
attribute list `a`. `data` can be a vector or a scalar. A scalar
is handeld as a vector with one element in the NetCDF data model.

Generally the attributes are defined by indexing, for example:

```julia
ds = NCDataset("file.nc","c")
ds.attrib["title"] = "my title"
close(ds)
```

If `data` is a string, then attribute is saved as a list of
NetCDF characters (`NC_CHAR`) with the appropriate length. To save the attribute
as a string (`NC_STRING`) you can use the following:

```julia
ds = NCDataset("file.nc","c")
ds.attrib["title"] = ["my title"]
close(ds)
```

"""
function Base.setindex!(a::Attributes,data,name::SymbolOrString)
    defmode(a.ds) # make sure that the file is in define mode
    return nc_put_att(a.ds.ncid,a.varid,name,data)
end

"""
    Base.keys(a::Attributes)

Return a list of the names of all attributes.
"""
Base.keys(a::Attributes) = listAtt(a.ds.ncid,a.varid)

"""
    Base.delete!(a::Attributes, name)

Delete the attribute `name` from the attribute list `a`.
"""
function Base.delete!(a::Attributes,name::SymbolOrString)
    defmode(a.ds)
    nc_del_att(a.ds.ncid,a.varid,name)
    return nothing
end

function Base.show(io::IO, a::BaseAttributes; indent = "  ")
    try
        # use the same order of attributes than in the NetCDF file
        for (attname,attval) in a
            print(io,indent,@sprintf("%-20s = ",attname))
            printstyled(io, @sprintf("%s",attval),color=attribute_color())
            print(io,"\n")
        end
    catch err
        if isa(err,NetCDFError)
            if err.code == NC_EBADID
                print(io,"NetCDF attributes (file closed)")
                return
            end
        end
        rethrow()
    end
end


text:
# bcbootlp.m
# Luca Brugnolini 2014
# Based on Kilian code
# This function gives bias-corrected blocks-of-blocks bootstrap 95%inverval
# for IRF by LP.
# Bias is estimated by bootstrap from the original data.
# For structural-form IRFs, bootstrap replication of SIGMA* is drawn from
# the asymptotic distribution of SIGMA hat estimated by VAR

function bcbootlp(z, rfirfl, SIGMA, paic, block, H)

# Input : data, reduced form IRF estimated by LP
#         SIGMA hat estimated from VAR, lag-order (Hx1 vector)
#         size of blocks for consecutive draws
# Output: bias-corrected blocks-of-blocks bootstrap 95% confidence interval

(T,K) = size(z);

#  Bootstrap for bias estimation
bias = zeros(K^2, H+1);

    for h = 1:H

        p = paic[h];                      # lag-order for horizon h
        p = itrunc(p);

        # Creating blocks of data from which bootstrap sample will be drawn
        # To preserve serial correlation of original data
        ys = z[p+h:T,:];                  # LHS variable: y(t+1)~y(t+H)
        yt = z[p:T-h,:];                  # RHS variable of interest: y(t)
        x = ones(T-p-h+1,1);                # constant term

        for i = 2:p
            x = [x z[p+1-i:T-i-h+1,:]];   # other regressor:y(t-1)~y(t-p+1)

        end

        bias[:,h+1] = boot_biaslp(ys, yt, x, rfirfl[:,h+1], block, 0);

    end

#  Bootstrap simulation for CI
nrep = 100;                             # number of bootstrap replications 2000
IRFrmat = zeros(nrep,K^2*(H+1));       # IRF for each bootstrap replication

     for j = 1:nrep                         # Starting bootstrap cycle
         IRFr = zeros(K^2,H+1);

         # Drawing SIGMA* from asymptotic distribution of SIGMA hat
         D = duplication(K);                                                    # duplication matrix
         Dplus = (D'*D)\D';                                                     # ***
         SIGSIG = 2*Dplus*kron(SIGMA,SIGMA)*Dplus';                             # SIGMA(vech(SIGMA))
         m = sum(1:K);                                                          # number of draws
         vechSIGMAr = chol(SIGSIG)'*randn(m,1)/sqrt(T)+vech(SIGMA);             # SIGMA(vech(SIGMA*))
         vecSIGMAr = duplication(K)*vechSIGMAr;                                 # vec(SIGMA*)
         SIGMAr = devec(vecSIGMAr,K);                                           # recover SIGMA*
         Pr = chol(SIGMAr)';                                                    # A0inv*
         IRFr[:,1] = vec(Pr');                                                  # IRF for horizon 0

         for h = 1:H                           # Horizons cycle
             p = paic[h];                      # lag-order for horizon h
             p = itrunc(p);

             # Creating blocks of data from which bootstrap sample will be drawn
             # To preserve serial correlation of original data
             ys = z[p+h:T,:];                  # LHS variable: y(t+1)~y(t+H)
             yt = z[p:T-h,:];                  # RHS variable of interest: y(t)
             x = ones(T-p-h+1,1);              # constant term

             for i = 2:p
                 x = [x z[p+1-i:T-i-h+1,:]];   # other regressor:y(t-1)~y(t-p+1)

             end

         # Drawing blocks of bootstrap data from the blocks of original data
         # ysr = []; ytr = []; xr = [];
         pos = ceil(rand()*(size(ys,1)-block))+1;
         ysr = ys[pos:pos+block-1,:];
         ytr = yt[pos:pos+block-1,:];
         xr  = x[pos:pos+block-1,:];

         for ii = 1:ceil(size(ys,1)/block)                               # Here was length(ys)
             pos = ceil(rand()*(size(ys,1)-block))+1;
             ysr = [ysr; ys[pos:pos+block-1,:]];
             ytr = [ytr; yt[pos:pos+block-1,:]];
             xr  = [xr; x[pos:pos+block-1,:]];

         end

        IRFr[:,h+1] = vec(Pr'*(irflpest(ysr, ytr, xr)-reshape(bias[:,h+1], K, K)));     # bias correction

     end                                                 # End horizons cycle

     IRFrmat[j,:] = vec(IRFr')';

end                                                   # End bootstrap cycle

# Calculate 95 percent interval endpoints
CILl = zeros(1,size(IRFrmat,2));
CIHl = zeros(1,size(IRFrmat,2));

for i = 1:size(IRFrmat,2)
    CILl[:,i] = quantile(vec(IRFrmat[:,i]),0.025);
    CIHl[:,i] = quantile(vec(IRFrmat[:,i]),0.975);

end

CI = [CILl; CIHl];

        return CI

end