text:
"""
    These functions read particles in a given range over multiple files, based
    on the relevant peano-hilbert keys.
    This code is based on read_particles_in_box.pro by Dr. habil. Klaus Dolag.
"""

using Dates
using Base.Threads

"""
    read_particles_in_box_peano(filename::String, blocks::Vector{String},
                                corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                                parttype::Integer=0, verbose::Bool=true)

Reads all particles within a box defined by a lower left and upper right corner
for a given particle type based on peano hilbert key reading. Returns a dictionary with all requested blocks.
"""
function read_particles_in_box_peano(filename::String, blocks::Vector{String},
                               corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                               parttype::Integer=0, verbose::Bool=true)

    if verbose
        println()
        @info "Running on $(nthreads()) threads"
    end

    filebase = filename
    file_key_index = filebase * ".key.index"

    # if the snapshot does not exist it may be split into multiple files
    if !isfile(filebase)

        if verbose
            @info "File: $filebase not found, looking for sub-files."
        end

        # try reading the first of the distributed snapshots
        filename = select_file(filebase, 0)

        h = head_to_obj(filename)

        if verbose
            @info "$(h.num_files) sub-files found."
        end

        nfiles = h.num_files
    else
        nfiles = 1
        h = head_to_obj(filename)
    end

    blocks, no_mass_block = check_blocks(filename, blocks)

    # read info blocks once here
    snap_info = read_info(filename)
    key_info  = read_info(filename * ".key")

    if verbose
        @info "All requested blocks present!"
        @info "Checking for .key files..."
    end

    # check if key files are present
    file_key = filename * ".key"
    if !block_present(file_key, "KEY")
        error("No .key file present!")
    end

    if verbose
        @info ".key files found!"
        @info "Calculating peano-hilbert keys..."
        t1 = Dates.now()
    end

    # first read the header
    h_key = read_keyheader(file_key)

    # get a list of the required peano-hilbert keys
    keylist = get_keylist(h_key, corner_lowerleft, corner_upperright)

    if verbose
        t2 = Dates.now()
        @info "$(size(keylist,1)) Peano-Hilbert keys found. Took: $(t2 - t1)"
        @info "Looking for relevant files..."
        t1 = Dates.now()
    end

    # find relevant files
    files = find_files_for_keys(filebase, nfiles, keylist)

    N_files = size(files,1)
    if verbose
        t2 = Dates.now()
        @info "$N_files files found. Took: $(t2 - t1)"

        @info "Searching read positions..."
        println()
        t1 = Dates.now()
    end

    # find all the positions where to read data
    file_offset_key, file_part_per_key, file_block_positions = find_read_positions( files, filebase, blocks,
                                                                                parttype, keylist, key_info,
                                                                                verbose)

    N_to_read = 0

    @inbounds for i = 1:N_files
        N_to_read += sum(file_part_per_key[i])
    end

    if verbose
        t2 = Dates.now()
        println()
        @info "Positions read. Took: $(t2 - t1)"
        println()
        @info "Reading $N_to_read particles..."
    end

    # prepare dictionary for particle storage
    d = allocate_data_dict(blocks, N_to_read, snap_info, no_mass_block)

    if verbose
        @info "Reading Blocks..."
        t1 = Dates.now()
    end

    n_read = 0

    for i = 1:N_files

        # select current file
        filename = select_file(filebase, files[i])

        # read header
        h = read_header(filename)

        # read info block
        snap_info = read_info(filename)

        # read blocks in parallel
        @threads for j = 1:size(blocks,1)

            block_info = snap_info[getfield.(snap_info, :block_name) .== blocks[j]][1]

            # add offset of particle types that should not be read
            offset = 0
            for k=1:size(h.npart,1)
                if (block_info.is_present[k] > 0) & (h.npart[k] > 0) & ( k < parttype + 1)
                    offset += h.npart[k]
                end
            end

            # reads data into the dictionary and counts up n_read
            read_block_with_offset!(d[blocks[j]], n_read, filename,
                                    file_block_positions[i][blocks[j]],
                                    block_info, offset, file_offset_key[i],
                                    file_part_per_key[i])

        end # loop over blocks

        n_read += sum(file_part_per_key[i])

        @info "Read $n_read / $N_to_read particles"

    end # for i = 1:size(files,1)

    # finally construct masses of no mass block present
    if no_mass_block
        d["MASS"] = h.massarr[parttype+1] .* ones(Float32, N_to_read)
    end

    if verbose
        t2 = Dates.now()
        @info "Blocks read. Took: $(t2 - t1)"
    end

    return d
end

"""
    read_particles_in_box(filename::String, blocks::String,
                          corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                          parttype::Integer=0, verbose::Bool=true,
                          use_keys::Bool=true)

Reads all particles within a box defined by a lower left and upper right corner
for a given particle type. Returns a dictionary with all requested blocks.
If `use_keys=true` it uses Peano-Hilbert keys to constrain the read-in, otherwise it uses a brute-force read-in with a filter function.
Peano-Hilbert key based read-in is significantly faster.
"""
function read_particles_in_box(filename::String, blocks::Vector{String},
                               corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                               parttype::Integer=0, verbose::Bool=true,
                               use_keys::Bool=true)

    if use_keys
        d = read_particles_in_box_peano(filename, blocks, corner_lowerleft, corner_upperright,
                                        parttype=parttype, verbose=verbose)
    else
        if verbose
            @info "Brute-force read-in."
        end
        filter_function(snap_file) = filter_cube(snap_file, corner_lowerleft, corner_upperright, parttype=parttype)
        d = read_blocks_over_all_files(filename, blocks, filter_function = filter_function, parttype = parttype, verbose = verbose )
    end

    return d
end

"""
    read_particles_in_box(filename::String, blocks::String,
                          corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                          parttype::Integer=0, verbose::Bool=true,
                          use_keys::Bool=true)

Like `read_particles_in_box` but for a single block. Returns the block as an array.

See also: [`read_particles_in_box`](@ref)
"""
function read_particles_in_box(filename::String, blocks::String,
                               corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real};
                               parttype::Integer=0, verbose::Bool=true,
                               use_keys::Bool=true)

    d = read_particles_in_box(filename, [blocks], corner_lowerleft, corner_upperright,
                              parttype=parttype, verbose=verbose, use_keys=use_keys)

    return d[blocks]
end

"""
    read_particles_in_box(filename::String, blocks::Vector{String},
                          center_pos, radius;
                          parttype::Integer=0, verbose::Bool=true,
                          use_keys::Bool=true)

Reads all particles within a box encapsulating a volume defined by center position
and radius for a given particle type. Returns a dictionary with all requested blocks.

See also: [`read_particles_in_box`](@ref)
"""
function read_particles_in_volume(filename::String, blocks::Vector{String},
                                  center_pos::Array{<:Real}, radius::Real;
                                  parttype::Integer=0, verbose::Bool=true,
                                  use_keys::Bool=true)

    # calculate lower left and upper right corner
    x0 = center_pos .- radius
    x1 = center_pos .+ radius

    return read_particles_in_box(filename, blocks, x0, x1, parttype=parttype, verbose=verbose, use_keys=use_keys)
end

"""
    read_particles_in_box(filename::String, blocks::String,
                          center_pos, radius;
                          parttype::Integer=0, verbose::Bool=true)

Like `read_particles_in_volume` but for a single block. Returns the block as an array.

See also: [`read_particles_in_volume`](@ref)
"""
function read_particles_in_volume(filename::String, blocks::String,
                                  center_pos, radius;
                                  parttype::Integer=0, verbose::Bool=true,
                                  use_keys::Bool=true)

    d = read_particles_in_volume(filename, [blocks], center_pos, radius,
                                 parttype=parttype, verbose=verbose,
                                 use_keys=use_keys)

    return d[blocks]
end


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.
=#

@doc """
A function knob is a symbol in compile time. For this symbol, we would like to
record a type for it in the Sym2TypeMap. This abstract type is for this purpose
only.
"""
abstract FunctionKnob

@doc """
The CallSites' extra field for context info discovery.
"""
type ContextExtra
    # Environment before the AST walking
    matrix_properties        :: Symexpr2PropertiesMap
    bb2depth                 :: Dict{BasicBlock, Int}

    # Flag if ast may be under change. We visite the AST twice. It is true
    # the first time (allowing pattern replacement; knobs are generated as well),
    # and false the second time (allowing gathering the knobs of the AST nodes).
    ast_may_change           :: Bool

    # Infomration gathered when AST is walked the first time
    reordering_decider_fknob :: Symbol # Remember fknob instead of AST here as AST may change due to call replacement, etc.
    reordering_decider_power :: Int    # 0 is no power: not to make reordering decision
    reordering_FAR           :: Vector # First Arrays Reordered by the decider
    expr2fknob               :: Dict{Expr, Symbol}
    fknob2expr               :: Dict{Symbol, Expr}
    fknob2mknobs             :: Dict{Symbol, Vector} # Vector[matrix knobs]
    fknob2ctor_dtor          :: Dict{Symbol, Tuple}  # Tuple(fknob_creator, fknob_deletor)
    matrix2mknob             :: Dict{Symexpr, Symbol}
    mknob2matrix             :: Dict{Symbol, Symexpr}
    derivatives              :: Set{Tuple} # Tuple(matrixA, derivative relation, matrixB)

    # The knobs gathered when AST is walked the second time.
    matrix_knobs             :: Set{Symbol}
    function_knobs           :: Set{Symbol}

    # Scratch variables.
    bb                       :: Any              # BasicBlock
    stmt_idx                 :: StatementIndex
    prev_stmt_idx            :: StatementIndex
    live_in_before_prev_expr :: Set{Sym}
    live_in_before_expr      :: Set{Sym}

    # Some patterns (those contain :t1!2, etc.) may hoist some subtrees of the
    # current statement before it. That splits one statement into more than one.
    # Such patterns should be matched at the last, because otherwise, other
    # patterns may not be able to match what they should: they cannot find the
    # subtrees to match, which are no longer in the same statement.
    non_splitting_patterns  :: Vector{Pattern}

    ContextExtra(_matrix_properties, _bb2depth, _ast_may_change) = new(
            _matrix_properties, _bb2depth, _ast_may_change, :nothing, 0, [],
            Dict{Expr, Symbol}(), Dict{Symbol, Expr}(), Dict{Symbol, Vector}(),
            Dict{Symbol, Tuple}(), Dict{Symexpr, Symbol}(),
            Dict{Symbol, Symexpr}(), Set{Tuple}(),
            Set{Symbol}(), Set{Symbol}(), nothing, 0, 0, Set{Sym}(), Set{Sym}(),
            Vector{Pattern}())
end

# Below are the patterns that we would like to replace a function with another,
# which is semantically equivalent, but is context-sensitive
# For short, CS represents "Context Sensitive".

@doc """
Pre-processing function of CS_fwd/bwdTriSolve!_pattern. It checks if the L/U
matrix has a proxy structure available or not. That proxy matrix has to be a
symmetric matrix (not just lower or upper part). This is due to the special
requirement of the library to fast build a dependence graph. If the matrix is just
lower or upper part, the library building the graph would be twice slower.
"""
function CS_fwdBwdTriSolve!_check(
    ast           :: Expr,
    call_sites    :: CallSites,
    fknob_creator :: AbstractString,
    fknob_deletor :: AbstractString
)
    arg2 = ast.args[2]
    if typeof(arg2) == SymbolNode
        L_or_U = arg2.name
    else
        L_or_U = arg2
    end
    assert(typeof(L_or_U) <: Symexpr)

    if !haskey(call_sites.extra.matrix_properties, L_or_U) ||
        (call_sites.extra.matrix_properties[L_or_U].lower_of == nothing &&
         call_sites.extra.matrix_properties[L_or_U].upper_of == nothing)
        return false
    end

    # It is not possible to be a lower of some matrix, and also an upper of some
    # (other) matrix.
    assert(call_sites.extra.matrix_properties[L_or_U].lower_of == nothing ||
           call_sites.extra.matrix_properties[L_or_U].upper_of == nothing)

    # Check it is lower or upper part of a symmetric matrix.
    if call_sites.extra.matrix_properties[L_or_U].lower_of != nothing
        M = get_symexpr(call_sites.extra.matrix_properties[L_or_U].lower_of)
    else
        M = get_symexpr(call_sites.extra.matrix_properties[L_or_U].upper_of)
    end

    if !haskey(call_sites.extra.matrix_properties, M) ||
       !call_sites.extra.matrix_properties[M].constant_structured
       !call_sites.extra.matrix_properties[M].is_structure_symmetric
        return false
    end

    if typeof(L_or_U) == Symbol && typeof(M) == Symbol
        push!(call_sites.extra.derivatives, (L_or_U, DERIVATIVE_TYPE_SYMMETRIC, M))
    end

    return true
end

@doc """
Post-processing function of a pattern: Gather context sensitive info of the call
site at the given AST.
"""
function gather_context_sensitive_info(
    ast               :: Expr,
    call_sites        :: CallSites,
    fknob_creator     :: AbstractString,
    fknob_deletor     :: AbstractString,
    matrices_to_track :: Tuple,
    reordering_power  :: Int,
    reordering_FAR    :: Tuple
)
    assert(call_sites.extra.ast_may_change)

    # Create a function-specific knob.
    fknob                                   = new_symbol("fknob")
    call_sites.extra.expr2fknob[ast]        = fknob
    call_sites.extra.fknob2expr[fknob]      = ast
    call_sites.extra.fknob2ctor_dtor[fknob] = (fknob_creator, fknob_deletor)
    call_sites.extra.fknob2mknobs[fknob]    = []

    # Record the new symbol with a FunctionKnob type so that we may distinguish
    # it from a normal argument.
    symbol_info        = call_sites.symbol_info
    symbol_info[fknob] = FunctionKnob

    # Get the previous and the current expression that would be used for
    # argument replacement below.
    bb            = call_sites.extra.bb
    stmt_idx      = call_sites.extra.stmt_idx
    prev_stmt_idx = call_sites.extra.prev_stmt_idx
    cur_expr      = bb.statements[stmt_idx].expr
    prev_expr     = prev_stmt_idx == 0 ? nothing : bb.statements[prev_stmt_idx].expr

    # Create matrix-specific knobs for the matrices inputs.
    args = ast.args
    for arg in matrices_to_track
        if arg == :r
            M = ast
        else
            new_arg = replacement_arg(arg, args, ast, call_sites, prev_expr, cur_expr)
            assert(typeof(new_arg) == Symbol || typeof(new_arg) == SymbolNode ||
                   typeof(new_arg) == GenSym)
            M = ((typeof(new_arg) == SymbolNode) ? new_arg.name : new_arg)
        end
        if !haskey(call_sites.extra.matrix2mknob, M)
            mknob = new_symbol(string("mknob", typeof(M) == Expr ? "Expr" : string(M)))
            call_sites.extra.matrix2mknob[M]     = mknob
            call_sites.extra.mknob2matrix[mknob] = M
        end
        push!(call_sites.extra.fknob2mknobs[fknob], call_sites.extra.matrix2mknob[M])
    end

    bb     = call_sites.extra.bb
    depth  = call_sites.extra.bb2depth[bb]
    weight = 2 << depth
    if reorder_enabled && depth > 0 && call_sites.extra.reordering_decider_power < (reordering_power * weight)
        # The decider should be in a basic block that is an immediate member of
        # a loop region.
        # The decider should not be in the function region: we want one
        # reordering of data to impact many loop iterations, not sequential code.
        call_sites.extra.reordering_decider_fknob = fknob
        call_sites.extra.reordering_decider_power = reordering_power * weight
        call_sites.extra.reordering_FAR           = []
        for arg in reordering_FAR
            new_arg = replacement_arg(arg, args, ast, call_sites, prev_expr, cur_expr)
            assert(typeof(new_arg) == SymbolNode || typeof(new_arg) <: Symexpr)
            M = ((typeof(new_arg) == SymbolNode) ? new_arg.name : new_arg)
            push!(call_sites.extra.reordering_FAR, M)
        end
    end

    return true
end

@doc """
Pre-processing function of CS_ADAT_pattern: check it is A*D*A'; D is a diagonal
array; and A has constant value.
"""
function CS_ADAT_check(
    ast           :: Expr,
    call_sites    :: CallSites,
    fknob_creator :: AbstractString,
    fknob_deletor :: AbstractString
)
    A = ast.args[2]
    assert(typeof(A) == SymbolNode)
    assert(typeof(ast.args[4].args[2]) == SymbolNode)
    if A.name != ast.args[4].args[2].name
        return false
    end

    if !haskey(call_sites.extra.matrix_properties, A.name) ||
        !call_sites.extra.matrix_properties[A.name].constant_valued
        return false
    end

    # TODO: check D is diagonal after structure analysis is done.
    return true
end

@doc """
Post-processing function of CS_ADAT_pattern.
This will create an action to transpose A before the loop head. Usually we do
not create any action during this phase (the first walk of the AST). This is
an exception.
"""
function CS_ADAT_post_replacement(
    ast               :: Expr,
    call_sites        :: CallSites,
    fknob_creator     :: AbstractString,
    fknob_deletor     :: AbstractString,
    matrices_to_track :: Tuple,
    reordering_power  :: Int,
    reordering_FAR    :: Tuple
)
    # We need to replace A * D * A' into ADB(A, D, A', fknob). At this
    # moment, it is in the middle form of ADB(A, D, A). That means:
    # (1) Make a symbol AT. Insert AT = A' before the loop (This is to hoist A'
    #     out of loop)
    # (2) Replace ADB(A, D, A) as ADB(A, D, AT).
    # fknob is not added for now, which will be added automatically later.
    action = (typeof(call_sites.region) == FunctionRegion) ?
        InsertBeforeOrAfterStatement(Vector{Statement}(), true, call_sites.region.entry, 1) :
        InsertBeforeLoopHead(Vector{Statement}(), call_sites.region.loop, true)
    push!(call_sites.actions, action)

    A = ast.args[4]
    assert(typeof(A) == SymbolNode)
    properties = call_sites.extra.matrix_properties[A.name]
    AT = Symbol(string("__", string(A.name), "T__"))
    stmt = Statement(0, Expr(:(=), AT, Expr(:call, GlobalRef(Main, :ctranspose), A)))
    push!(action.new_stmts, stmt)

    ast.args[4] = SymbolNode(AT, A.typ)

    # AT has the same properties as A.
    call_sites.extra.matrix_properties[AT]               = properties
    call_sites.extra.matrix_properties[AT].is_single_def = false # AT definition is out of the region, and not counted as def.

    push!(call_sites.extra.derivatives, (A.name, DERIVATIVE_TYPE_TRANSPOSE, AT))

    # The result is a symmetric matrix.
    # ISSUE: in analysis like constant value/structure analysis, should we
    # consider each AST node instead of only symbols? If that is the case, we
    # do not have to make properties here.
    call_sites.extra.matrix_properties[ast]                        = MatrixProperties()
    call_sites.extra.matrix_properties[ast].constant_valued        = false #TODO: should determined by if A and D are constant valued
    call_sites.extra.matrix_properties[ast].constant_structured    = true
    call_sites.extra.matrix_properties[ast].is_symmetric           = false #true
    call_sites.extra.matrix_properties[ast].is_structure_symmetric = false #true
    call_sites.extra.matrix_properties[ast].is_structure_only      = false
    call_sites.extra.matrix_properties[ast].is_single_def          = true

    return gather_context_sensitive_info(ast, call_sites, fknob_creator,
                                         fknob_deletor, matrices_to_track,
                                         reordering_power, reordering_FAR)
end

@doc """
Pre-processing function of SpSquareWithEps_pattern1: check it is X*X',
and X is symmetric valued.
"""
function SpSquareWithEps_check(
    ast           :: Expr,
    call_sites    :: CallSites,
    fknob_creator :: AbstractString,
    fknob_deletor :: AbstractString
)
    X = ast.args[2]
    assert(typeof(X) == SymbolNode)
    assert(typeof(ast.args[3].args[2]) == SymbolNode)
    if X.name != ast.args[3].args[2].name
        return false
    end

    if !haskey(call_sites.extra.matrix_properties, X.name) ||
        !call_sites.extra.matrix_properties[X.name].is_symmetric
        return false
    end

    return true
end

@doc """
The following patterns are to match the following source code
    for
        ...
        B = A*D*A' // A and D::SparseMatrixCSC{Float64,Int32}, D is diagnoal matrix
        ...
        R = cholfact_int32(B)
        ...
        dy = R\t2

Then for each function call site, create a fknob. For each of them, for each
SparseMatrix input and output, add a mknob to represent them. Even though
a matrix might have not been allocated yet, we can still represent it
with its mknob. The purpose of the fknobs are to let us find the mkobs shared
accross functions.

Then for each function call site, we check some properties of the input mknobs;
if they hold, we check th existence of some data of the output mknob; if they do
not exist, create them; otherwise, directly use them.

This procedure is general, and not limited to this specific case.

The above code is transformed into
    # dy is used in cholmod_factor_inverse_divide(), and the pattern knows that
    # we need to create dy ahead of time, if dy is NOT live befor the loop.
    dy = Array(Cdouble, m)

    create mknob for A, B, R
    create fknob for ADB, cholfact_int32, and cholmod_factor_inverse_divide
    add mknob_A, B into fknob_ADB
    add mknob_B, R into fknob_cholfact_int32
    add mknob_R into fknob_cholmod_factor_inverse_divide
    for
        ...
        # Sparso.ADB(A', D, A, fknob_ADB)
        #   if fknob_ADB == NULL
        #       return A * D * A'
        #   get the input mknob_A from fknob_ADB
        #   if !mknob_A->constant_structured (This check is actually already done by pattern matching)
        #       return A * D * A'
        #   get the output mknob_B from fknob_ADB
        #   if mknob_B->structure_proxy == NULL
        #       ADAT = adb_inspect(A', A),
        #       mknob_B->structure_proxy = ADAT
        #       mknob_B->matrix = ADAT
        #   d = diag(D)
        #   adb_execute!(mknob_B->matrix, A', A, d)
        B = Sparso.ADB(A', D, A, fknob_ADB)
        ...
        # Sparso.cholfact_int32(B, fknobcholfact_int32)
        #   if fknobcholfact_int32 == NULL
        #       return cholfact_int32(B)
        #   get the input mknob_B from fknob_cholfact_int32
        #   if mknob_B->structure_proxy == NULL || mknob_B->matrix == NULL
        #       return cholfact_int32(B)
        #   get the output mknob_R from fknob_cholfact_int32
        #   if mknob_R->dss_handle == NULL
        #       dss_handle = dss_analyze(mknob_B->structure_proxy)
        #       mknob_R->dss_handle = dss_handle
        #   dss_factor(mknob_R->dss_handle, mknob_B->matrix)
        #   BUT HOW TO RETURN R?

        R = Sparso.cholfact_int32(B, fknob_cholfact_int32)
        ...

        # cholmod_factor_inverse_divide(dy, R, t2, fknob_cholmod_factor_inverse_divide)
        #   if fknob_cholmod_factor_inverse_divide == NULL
        #       return R \ t2
        #   get the input mknob_R from fknob_cholmod_factor_inverse_divide
        #   if mknob_R->dss_handle == NULL
        #       return R \ t2
        #   else
        #       opt = MKL_DSS_DEFAULTS
        #       dss_solve!(mknob_R->dss_handle, t2, dy)

        dy = cholmod_factor_inverse_divide(R, t2, fknob_cholmod_factor_inverse_divide)
"""

const CS_transpose_pattern = ExprPattern(
    "CS_transpose_pattern",
    (:call, GlobalRef(Main, :ctranspose), SparseMatrixCSC),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:NO_CHANGE, ),
    do_nothing,
    "",
    "",
    (),
    0,
    ()
)

const CS_ADAT_pattern = ExprPattern(
    "CS_ADAT_pattern",
    (:call, GlobalRef(Main, :*), SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC),
    (nothing, nothing, nothing, nothing, CS_transpose_pattern),
    CS_ADAT_check,
    (:call, GlobalRef(Sparso, :ADB), :a2, :a3, :a2),
    CS_ADAT_post_replacement,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:r, :a2, :a3, :a4),
    0,
    ()
)

const CS_cholfact_int32_pattern = ExprPattern(
    "CS_cholfact_int32_pattern",
    (:call, GlobalRef(Main, :cholfact_int32), SparseMatrixCSC{Float64, Int32}),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :cholfact_int32), :a2),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a2,),
    0,
    ()
)

const CS_cholsolve_pattern = ExprPattern(
    "CS_cholsolve_pattern",
    (:call, GlobalRef(Main, :\), Factorization{Float64}, Any),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :cholfact_inverse_divide), :a2, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a2,),
    0,
    ()
)

@doc """
Pattern for foward triangular solver.
"""
const CS_fwdTriSolve!_pattern = ExprPattern(
    "CS_fwdTriSolve!_pattern",
    (:call, TypedExprNode(Function, :call, TopNode(:getfield), GlobalRef(Base, :SparseMatrix), QuoteNode(:fwdTriSolve!)),
      SparseMatrixCSC, Vector),
    (:NO_SUB_PATTERNS,),
    CS_fwdBwdTriSolve!_check,
    (:call, GlobalRef(Sparso, :fwdTriSolve!), :a2, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a2,),
    100,
    (:a2, :a3)
)

@doc """
Pattern for forward triangular solver (A lower triangular matrix \ a vector). The
matrix's structure is the lower part of a constant symmetric matrix's structure.
"""
const CS_fwdTriSolve!_backslash_pattern = ExprPattern(
    "CS_fwdTriSolve!_backslash_pattern",
    (:(=), Vector,
      Expr(:call, GlobalRef(Main, :\),
            AD(SparseMatrixCSC, (SA_LOWER_OF, SA_CONST_STRUCTURED | SA_SYMM_STRUCTURED)),
            Vector)),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :fwdTriSolve!), :a1, :a2_2, :a2_3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    100,
    (:a3, :a2)
)

@doc """
Pattern for backward triangular solver.
"""
const CS_bwdTriSolve!_pattern = ExprPattern(
    "CS_bwdTriSolve!_pattern",
    (:call, TypedExprNode(Function, :call, TopNode(:getfield), GlobalRef(Base, :SparseMatrix), QuoteNode(:bwdTriSolve!)),
      Any, Vector), #SparseMatrixCSC, Vector), # Somehow, type inference does not figure out U's type is SparseMatrixCSC
    (:NO_SUB_PATTERNS,),
    CS_fwdBwdTriSolve!_check,
    (:call, GlobalRef(Sparso, :bwdTriSolve!), :a2, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a2,),
    90,
    (:a2, :a3)
)

@doc """
Pattern for backward triangular solver (A upper triangular matrix \ a vector). The
matrix's structure is the upper part of a constant symmetric matrix's structure.
"""
const CS_bwdTriSolve!_backslash_pattern = ExprPattern(
    "CS_bwdTriSolve!_backslash_pattern",
    (:(=), Vector,
      Expr(:call, GlobalRef(Main, :\),
            AD(SparseMatrixCSC, (SA_UPPER_OF, SA_CONST_STRUCTURED | SA_SYMM_STRUCTURED)),
            Vector)),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :bwdTriSolve!), :a1, :a2_2, :a2_3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    100,
    (:a3, :a2)
)

@doc """
The following patterns are for SpMV. The difference from the SpMV patterns in
call-replacement.jl is that here we will add mknobs and fknob to SpMV.
"""

@doc """ A * x => SpMV(1, A, x) """
const CS_SpMV_pattern1 = ExprPattern(
    "CS_SpMV_pattern1",
    (:call, GlobalRef(Main, :*), SparseMatrixCSC, Vector),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV), 1, :a2, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    40,
    (:a3, :a4, :r)
)

@doc """ a * A * x + g => SpMV(a, A, x, 0, x, g) """
const CS_SpMV_pattern2 = ExprPattern(
    "CS_SpMV_pattern2",
    (:call, GlobalRef(Main, :+), Vector, Number),
    (nothing, nothing, number_times_matrix_vector_pattern, nothing),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV), :a2_2, :a2_3, :a2_4, 0, :a2_4, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    40,
    (:a3, :a4, :r)
)

@doc """ a * A * x => SpMV(a, A, x) """
const CS_SpMV_pattern3 = ExprPattern(
    "CS_SpMV_pattern3",
    (:call, GlobalRef(Main, :*), Number, SparseMatrixCSC, Vector),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV), :a2, :a3, :a4),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    40,
    (:a3, :a4, :r)
)

@doc """
SpMV(a, A, x) + g => SpMV(a, A, x, 0, x, g)
Note: its reordering power is bigger than CS_SpMV_pattern3, because
CS_SpMV_pattern3 generates a SpMV_3_parameters_pattern, which is a sub-pattern
of this pattern. Making reordering power of this pattern bigger ensures that
when the sub-pattern is replaced, the reordering decider becomes the AST
generated by this pattern, and no more being that by the sub-pattern.
"""
const CS_SpMV_pattern4 = ExprPattern(
    "CS_SpMV_pattern4",
    (:call, GlobalRef(Main, :+), Vector, Number),
    (nothing, nothing, SpMV_3_parameters_pattern, nothing),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV),
     :a2_2, :a2_3, :a2_4, 0, :a2_4, :a3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    50,
    (:a3, :a4, :r)
)

@doc """
SpMV(a, A, x) + b => SpMV(a, A, x, 1, b, 0)
Note: its reordering power is bigger than CS_SpMV_pattern3, because
CS_SpMV_pattern3 generates a SpMV_3_parameters_pattern, which is a sub-pattern
of this pattern. Making reordering power of this pattern bigger ensures that
when the sub-pattern is replaced, the reordering decider becomes the AST
generated by this pattern, and no more being that by the sub-pattern.
"""
const CS_SpMV_pattern5 = ExprPattern(
    "CS_SpMV_pattern5",
    (:call, GlobalRef(Main, :+), Vector, Vector),
    (nothing, nothing, SpMV_3_parameters_pattern, nothing),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV),
     :a2_2, :a2_3, :a2_4, 1, :a3, 0),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    50,
    (:a3, :a4, :a6, :r)
)

@doc """
SpMV(a, A, x) - b => SpMV(a, A, x, -1, b, 0)
Note: its reordering power is bigger than CS_SpMV_pattern3, because
CS_SpMV_pattern3 generates a SpMV_3_parameters_pattern, which is a sub-pattern
of this pattern. Making reordering power of this pattern bigger ensures that
when the sub-pattern is replaced, the reordering decider becomes the AST
generated by this pattern, and no more being that by the sub-pattern.
"""
const CS_SpMV_pattern6 = ExprPattern(
    "CS_SpMV_pattern6",
    (:call, GlobalRef(Main, :-), Vector, Vector),
    (nothing, nothing, SpMV_3_parameters_pattern, nothing),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV),
     :a2_2, :a2_3, :a2_4, -1, :a3, 0),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    50,
    (:a3, :a4, :a6, :r)
)

@doc """
    (-1/a)*A*x + b*y => SpMV(-1/a, A, x, b, y)
"""
const CS_SpMV_pattern7 = ExprPattern(
    "CS_SpMV_pattern7",
    (:call, GlobalRef(Main, :+),
      Expr(:call, GlobalRef(Main, :/),
        Expr(:call, GlobalRef(Main, :-),
            Expr(:call, GlobalRef(Main, :*), SparseMatrixCSC, Vector)),
        Number),
      Expr(:call, GlobalRef(Main, :*), Number, Vector)),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV),
      TypedExprNode(Float64, :call, GlobalRef(Main, :/), -1, :a2_3), :a2_2_2_2, :a2_2_2_3, :a3_2, :a3_3),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    70,
    (:a3, :a4, :a6, :r)
)

@doc """
The following patterns are potentially related with SpMV!. In SpMP library, SpMV!
is implemented in such a way that the result and the inputs are written and read
in parallel. Thus the result and the inputs cannot share space to avoid any race
condition.

There are 3 cases:

(Case 1): y = A *x ==> y = SpMV(A, x)

This is equivalent to
        result = allocate memory of size(A, 1)
        SpMV!(result, A, x)
        y = result
The first two statements have been realized by SpMV. Thus we can simply generate
code like this
        y = SpMV(A, x)

One might think it is great to hoist the allocation out of the current region (
a loop). However, that can cause subtle issues. Consider a simplified pagerank:
        y = A * x
        if |y - x| < err exit
        x = y

If we hoist the allocation out, the code becomes:
    Before the loop:
        result = allocate memory of size(A, 1)
    In the loop
        SpMV!(result, A, x)
        y = result
        if |y - x| < err exit
        x = y
After the first iteration, x, y, and result share the same memory. Thus |y-x|
is always 0 afterwards.

(Case 2) x = A * x  ==> x = SpMV(A, x) or SpMV!(x, A, x) with library extension.

Here we are sure that x has been allocated memory. There are three ways we can
generate code:

CodeGen 21: generate
        x = SpMV(A, x)
The overhead is: memory allocation and copy every time SpMV! is called.

CodeGen 22: generate
        SpMV!(x, A, x)
This requires that the library is extended to allow that the result and
the input are aliased. To do that, library identifies such an alias, and
creates a temporary memory space t as result of A * x, and run SpMV!(t, A, x),
where it can write and read inputs in parallel, and finally copy t's content to
x. The overhead is the same: memory allocation and copy every time SpMV! is
called.

CodeGen 23: generate
    Before the loop:
        result = allocate memory of size(A, 1)
        assert(result and x have different spaces)
    In the loop:
        SpMV!(result, A, x)
        tmp = x
        x = result
        result = tmp
        assert(result and x have different spaces)
This has 1 allocation, and no copy back.

So far, we do not have Case 2 in any benchmark yet. So let us leave it for future.
Simply follow CodeGen 21 or 22 for now.

(Case 3) A_mul_B!(y, A, x) ==> SpMV!(y, A, x) with library extension.

This is equivalent to y = A * x, but y is sure to have been allocated memory
before.

If y is lexically x, it is a special case A_mul_B!(x, A, x). Same as Case 2.

Otherwise, we need to know if y and x are aliased (share the same memory).
If they are aliased, that is Case 2.
If they are not, simply call SpMV!(y, A, x) as usual.

So far, we do not rely on static alias analysis. So probably we should extend
libary to identify and allow the result and the input being aliased. See notes
for CodeGen 22.
"""

@doc """
A_mul_B!(y, A, x) = SpMV!(y, A, x)
"""
const CS_SpMV!_pattern1 = ExprPattern(
    "CS_SpMV!_pattern1",
    (:call, GlobalRef(Main, :A_mul_B!), Vector, SparseMatrixCSC, Vector),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :a2, :a3, :a4),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a3,),
    50,
    (:a3, :a4, :a2) # Put seed (A) as the first
)

@doc """ z = SpMV(a, A, x, b, y, g), z is x or y => SpMV!(z, a, A, x, b, y, g) """
const CS_SpMV!_pattern2 = ExprPattern(
    "CS_SpMV!_pattern2",
    (:(=), Vector, Vector),
    (nothing, nothing, SpMV_6_parameters_pattern),
    # NOTE: this is a HACK! We should not do this without knowing that z has memory
    # allocated. This is only to enable using SpMV! for performance.
    # TODO: uncomment this condition below after Linxiang's analysis works
    # LHS_in_RHS,
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :a1, :a2_2, :a2_3, :a2_4, :a2_5, :a2_6, :a2_7),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a4,),
    80,
    (:a4, :a2, :a5, :a7) # Put seed (A) at the first
)

@doc """ w = a * A * x => SpMV!(w, a, A, x, 0, x, 0) """
const CS_SpMV!_pattern3 = ExprPattern(
    "CS_SpMV!_pattern3",
    (:(=), AD(Vector, SA_HAS_FREE_MEMORY),
        Expr(:call, GlobalRef(Sparso, :SpMV),
              Number, SparseMatrixCSC, Vector)),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :a1, :a2_2, :a2_3, :a2_4, 0.0, :a1, 0.0),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a4,),
    60,
    (:a4, :a2, :a5) # Put seed (A) at the first
)

@doc """ A_mul_B!(a, A, x, b, y) => SpMV!(y, a, A, x, b, y , 0) """
const CS_SpMV!_pattern4 = ExprPattern(
    "CS_SpMV!_pattern4",
    (:call, GlobalRef(Main, :A_mul_B!), Number, SparseMatrixCSC, Vector, Number, Vector),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :a6, :a2, :a3, :a4, :a5, :a6, 0.0),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a4,),
    70,
    (:a4, :a2, :a5, :a7) # Put seed (A) at the first
)

@doc """ z = SpMV(a, A, x, b, y), z has dedicated memory => SpMV!(z, a, A, x, b, y) """
const CS_SpMV!_pattern5 = ExprPattern(
    "CS_SpMV!_pattern5",
    (:(=), AD(Vector, SA_HAS_FREE_MEMORY),
        Expr(:call, GlobalRef(Sparso, :SpMV),
             Number, SparseMatrixCSC, Vector, Number, Vector)),
    (:NO_SUB_PATTERNS,),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :a1, :a2_2, :a2_3, :a2_4, :a2_5, :a2_6, 0),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a4,),
    80,
    (:a4, :a2, :a5, :a7) # Put seed (A) at the first
)

@doc """ spmatmul_witheps(X, X', eps), X is symmetric => SpSquareWithEps(X, eps) """
const SpSquareWithEps_pattern1 = ExprPattern(
    "SpSquareWithEps_pattern1",
    (:call,   GlobalRef(Main, :spmatmul_witheps), SparseMatrixCSC, SparseMatrixCSC,      Number),
    (nothing, nothing,                            nothing,         CS_transpose_pattern, nothing),
    SpSquareWithEps_check,
    (:call, GlobalRef(Sparso, :SpSquareWithEps),
     :a2, :a4),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:a2,),
    0,
    ()
)

@doc """" Patterns that will actually transform the code. """
CS_transformation_patterns = [
    CS_ADAT_pattern,
    CS_cholfact_int32_pattern,
    CS_cholsolve_pattern,
    CS_fwdTriSolve!_pattern,
    CS_fwdTriSolve!_backslash_pattern,
    CS_bwdTriSolve!_pattern,
    CS_bwdTriSolve!_backslash_pattern,
    CS_SpMV_pattern7,
    CS_SpMV_pattern1,
    CS_SpMV_pattern2,
    CS_SpMV_pattern3,
    CS_SpMV_pattern4,
    CS_SpMV_pattern5,
    CS_SpMV_pattern6,
    CS_SpMV!_pattern1,
    CS_SpMV!_pattern2,
    CS_SpMV!_pattern3,
    # CS_SpMV!_pattern4] # Do not handle y = A * x + y case for now
    CS_SpMV!_pattern5,
    SpSquareWithEps_pattern1
]

@doc """
    Sparso.SpMV!(z, a, A, x, b, y, r)
    z = z .* u
=>
    Sparso.SpMV!(z, a, A, x, b, y, r, u)
    z = z # Useless
"""
const CS_SpMV!_two_statements_pattern1 = TwoStatementsPattern(
    "CS_SpMV!_two_statements_pattern1",
    (:call, GlobalRef(Sparso, :SpMV!),
      Vector, Number, SparseMatrixCSC, Vector, Number, Vector, Number),
    (:(=), :f2, Expr(:call, GlobalRef(Main, :.*), :f2, Vector)),
    do_nothing,
    (:call, GlobalRef(Sparso, :SpMV!),
     :f2, :f3, :f4, :f5, :f6, :f7, :f8, :s2_3),
    (:(=), :f2, :f2),
    gather_context_sensitive_info,
    "NewFunctionKnob",
    "DeleteFunctionKnob",
    (:f4,),
    100,
    (:f4, :f2, :f5, :f7, :f9) # Put seed (A) at the first
)

CS_two_statements_patterns = [
    CS_SpMV!_two_statements_pattern1,
]

@doc """
Create statements that will create a matrix knob for matrix M.
"""
function create_new_matrix_knob(
    new_stmts  :: Vector{Statement},
    mknob      :: Symbol,
    M          :: Symexpr,
    call_sites :: CallSites
)
    properties = call_sites.extra.matrix_properties[M]
    new_stmt = Expr(:(=), mknob,
                    Expr(:call, GlobalRef(Sparso, :new_matrix_knob),
                         typeof(M) == Expr ? QuoteNode(symbol("Expr")) : QuoteNode(M),
                         properties.constant_valued,
                         properties.constant_structured,
                         properties.is_symmetric,
                         properties.is_structure_symmetric,
                         properties.is_structure_only,
                         properties.is_single_def))
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Create a statement to propagate matrix info from one matrix knob to another.
"""
function create_propagate_matrix_info(
    to_mknob   :: Symbol,
    from_mknob :: Symbol,
    bb         :: BasicBlock,
    stmt_idx   :: StatementIndex,
    call_sites :: CallSites
)
    # Insert the propagation after the (assignment) statement.
    action = InsertBeforeOrAfterStatement(Vector{Statement}(), false, bb, stmt_idx)
    stmt = Statement(0,
                     Expr(:call, GlobalRef(Sparso, :propagate_matrix_info),
                           to_mknob, from_mknob))
    push!(action.new_stmts, stmt)
    push!(call_sites.actions, action)
end

@doc """
Create statements that will create a function knob for the call site, and add
the function knob to the call as a parameter.
"""
function create_new_function_knob(
    new_stmts     :: Vector{Statement},
    fknob         :: Symbol,
    fknob_creator :: AbstractString,
)
    assert(fknob_creator != "")

    if fknob_creator == "NewFunctionKnob"
        # Create a function knob. It has no private info specific to that function.
        # Call the parameterless version of new_function_knob for cleaner code.
        new_stmt = Expr(:(=), fknob,
                        Expr(:call, GlobalRef(Sparso, :new_function_knob))
                   )
    else
        # So far, we do not create any special knob. So this case is never reached.
        assert(false)
        new_stmt = Expr(:(=), fknob,
                        Expr(:call, GlobalRef(Sparso, :new_function_knob),
                              fknob_creator
                        )
                   )
    end
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Create statements that add the matrix knob to the function knob.
"""
function create_add_mknob_to_fknob(
    new_stmts :: Vector{Statement},
    mknob     :: Symbol,
    fknob     :: Symbol
)
    new_stmt = Expr(:call, GlobalRef(Sparso, :add_mknob_to_fknob), mknob, fknob)
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Create statements that will delete the function knob.
"""
function create_delete_function_knob(
    new_stmts     :: Vector{Statement},
    fknob         :: Symbol,
    fknob_deletor :: AbstractString
)
    assert(fknob_deletor != "")

    if fknob_deletor == "DeleteFunctionKnob"
        # Delete a function knob, which has no private info specific to that function.
        # Call the 1-parameter version of delete_function_knob for cleaner code.
        new_stmt = Expr(:call, GlobalRef(Sparso, :delete_function_knob),
                         fknob)
    else
        # So far, we have not created any special knob. So this case is never reached.
        assert(false)
        new_stmt = Expr(:call, GlobalRef(Sparso, :delete_function_knob),
                         fknob_deletor, fknob)
    end
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Create statements that will delete the matrix knob.
"""
function create_delete_matrix_knob(
    new_stmts :: Vector{Statement},
    mknob     :: Symbol
)
    new_stmt = Expr(:call, GlobalRef(Sparso, :delete_matrix_knob), mknob)
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Create statements that set mknob2 as a derivative of mknob1.
"""
function create_set_derivative(
    new_stmts :: Vector{Statement},
    mknob1    :: Symbol,
    relation  :: Int,
    mknob2    :: Symbol
)
    new_stmt = Expr(:call, GlobalRef(Sparso, :set_derivative),
                     mknob1, int2derivative_map[relation], mknob2)
    push!(new_stmts, Statement(0, new_stmt))
end

@doc """
Visit each expression the loop region. If recursive is true, visit each AST node
of the expression and handle each with the handler. Otherwise, handle the
expression with the handler.
"""
function visit_expressions(
    region                 :: Region,
    cfg                    :: CFG,
    call_sites             :: CallSites,
    recursive              :: Bool,
    handler                :: Function,
    two_statements_handler :: Function = do_nothing
)
    members     = (typeof(region) == FunctionRegion) ? region.members : region.loop.members
    blocks      = cfg.basic_blocks
    symbol_info = call_sites.symbol_info
    liveness    = call_sites.liveness
    for bb_idx in members
        bb                             = blocks[bb_idx]
        statements                     = bb.statements
        call_sites.extra.bb            = bb
        prev_expr                      = nothing
        call_sites.extra.prev_stmt_idx = 0
        for stmt_idx in 1 : length(statements)
            call_sites.extra.stmt_idx = stmt_idx
            stmt                      = statements[stmt_idx]
            expr                      = stmt.expr
            if typeof(expr) != Expr
                continue
            end

            call_sites.extra.live_in_before_expr = LivenessAnalysis.live_in(stmt, liveness)

            if recursive
                CompilerTools.AstWalker.AstWalk(expr, handler, call_sites)
            else
                handler(expr, call_sites)
            end

            if prev_expr != nothing && two_statements_handler != do_nothing
                # Try to merge this and the previous expression
                two_statements_handler(prev_expr, expr, call_sites, CS_two_statements_patterns)
            end
            prev_expr                                 = expr
            call_sites.extra.prev_stmt_idx            = stmt_idx
           call_sites.extra.live_in_before_prev_expr = call_sites.extra.live_in_before_expr
        end
    end
end

@doc """"
Add the function and matrix knobs of the AST node to
call_sites.extra.matrix_knobs and function_knobs.
"""
function gather_knobs(
    ast        :: Any,
    call_sites :: CallSites
)
    assert(!call_sites.extra.ast_may_change)

    if typeof(ast) != Expr
        return
    end

    if haskey(call_sites.extra.expr2fknob, ast)
        fknob = call_sites.extra.expr2fknob[ast]
        push!(call_sites.extra.function_knobs, fknob)

        mknobs = call_sites.extra.fknob2mknobs[fknob]
        union!(call_sites.extra.matrix_knobs, mknobs)

        site = CallSite(ast)
        push!(call_sites.sites, site)
    end
    return nothing
end

@doc """ A handler to visit_expressions(). """
gather_knobs(ast, call_sites :: CallSites, top_level_number, is_top_level, read) =
    gather_knobs(ast, call_sites)

@doc """ Create statements that will delete matrix and function knobs. """
function delete_knobs(
    action         :: Action,
    matrix_knobs   :: Set{Symbol},
    function_knobs :: Set{Symbol},
    call_sites     :: CallSites
)
    for mknob in matrix_knobs
        create_delete_matrix_knob(action.new_stmts, mknob)
    end

    for fknob in function_knobs
        fknob_creator, fknob_deletor = call_sites.extra.fknob2ctor_dtor[fknob]
        create_delete_function_knob(action.new_stmts, fknob, fknob_deletor)
    end
end

@doc """
Create statements that will create and intialize a knob for a matrix/function
before the loop, and statements that will delete all the matrix/function
knobs at each exit of the loop.
"""
function generate_and_delete_knobs(
    region     :: Region,
    call_sites :: CallSites
)
    region         = call_sites.region
    matrix_knobs   = call_sites.extra.matrix_knobs
    function_knobs = call_sites.extra.function_knobs
    derivatives    = call_sites.extra.derivatives
    symbol_info    = call_sites.symbol_info

    action_before_region = (typeof(region) == FunctionRegion) ?
        InsertBeforeOrAfterStatement(Vector{Statement}(), true, region.entry, 1) :
        InsertBeforeLoopHead(Vector{Statement}(), region.loop, true)
    push!(call_sites.actions, action_before_region)

    # Add statements that generates mknobs and fknobs
    for mknob in matrix_knobs
        M = call_sites.extra.mknob2matrix[mknob]
        create_new_matrix_knob(action_before_region.new_stmts, mknob, M, call_sites)
    end

    for fknob in function_knobs
        fknob_creator, fknob_deletor = call_sites.extra.fknob2ctor_dtor[fknob]
        create_new_function_knob(action_before_region.new_stmts, fknob, fknob_creator)
    end

    # Derivative information is already included in matrix_knob, so we copy it into derivatives
    for mknob in matrix_knobs
        M = call_sites.extra.mknob2matrix[mknob]
        prop = call_sites.extra.matrix_properties[M]
        if prop.lower_of != nothing && haskey(call_sites.extra.matrix2mknob, prop.lower_of)
            assert(prop.upper_of == nothing)
            push!(derivatives, (M, DERIVATIVE_TYPE_LOWER_TRIANGULAR, prop.lower_of))

            # ISSUE: is the following condition right for setting symmetric derivative? Do we
            # need both prop1.is_symmetric && prop1.is_structure_symmetric? Is this derivative
            # symmetric in structure or value?
            # ISSUE: without setting the symmetric derivative, pcg-illustration causes spmp
            # segmentation fault
            prop1 = call_sites.extra.matrix_properties[prop.lower_of]
            if prop1.is_symmetric && prop1.is_structure_symmetric
                push!(derivatives, (M, DERIVATIVE_TYPE_SYMMETRIC, prop.lower_of))
            end
        elseif prop.upper_of != nothing && haskey(call_sites.extra.matrix2mknob, prop.upper_of)
            push!(derivatives, (M, DERIVATIVE_TYPE_UPPER_TRIANGULAR, prop.upper_of))

            # ISSUE: is the following condition right for setting symmetric derivative? Do we
            # need both prop1.is_symmetric && prop1.is_structure_symmetric? Is this derivative
            # symmetric in structure or value?
            # ISSUE: without setting the symmetric derivative, pcg-illustration causes spmp
            # segmentation fault
            prop1 = call_sites.extra.matrix_properties[prop.upper_of]
            if prop1.is_symmetric && prop1.is_structure_symmetric
                push!(derivatives, (M, DERIVATIVE_TYPE_SYMMETRIC, prop.upper_of))
            end
        end
        if prop.transpose_of != nothing && haskey(call_sites.extra.matrix2mknob, prop.transpose_of)
            push!(derivatives, (M, DERIVATIVE_TYPE_TRANSPOSE, prop.transpose_of))
        end
    end

    # Add statements that define the relationship between matrix knobs
    for (M1, relation, M2) in derivatives
        mknob1 = call_sites.extra.matrix2mknob[M1]
        mknob2 = call_sites.extra.matrix2mknob[M2]
        create_set_derivative(action_before_region.new_stmts, mknob1, relation, mknob2)
    end

    if (typeof(region) == FunctionRegion)
        for pred_bb in region.exit.preds
            len = length(pred_bb.statements)
            assert(len > 0)
            last_stmt = pred_bb.statements[len]
            last_expr = last_stmt.expr
            assert(typeof(last_expr) <: Expr && last_expr.head == :return)

            action  = InsertBeforeOrAfterStatement(Vector{Statement}(), true, pred_bb, len)
            push!(call_sites.actions, action)
            delete_knobs(action, matrix_knobs, function_knobs, call_sites)
        end
    else
        # Create statements that will delete the knobs at region exits
        for exit in region.exits
            action  = InsertOnEdge(Vector{Statement}(), exit.from_bb, exit.to_bb)
            push!(call_sites.actions, action)
            delete_knobs(action, matrix_knobs, function_knobs, call_sites)
        end
    end
end

@doc """"
Create statements before the region that add matrix knobs into the function knob.
"""
function add_matrix_knobs_to_function_knobs(
    call_sites :: CallSites
)
    assert(!call_sites.extra.ast_may_change)

    region               = call_sites.region
    action_before_region = (typeof(region) == FunctionRegion) ?
        InsertBeforeOrAfterStatement(Vector{Statement}(), true, region.entry, 1) :
        InsertBeforeLoopHead(Vector{Statement}(), region.loop, true)
    push!(call_sites.actions, action_before_region)

    for fknob in call_sites.extra.function_knobs
        mknobs = call_sites.extra.fknob2mknobs[fknob]
        for mknob in mknobs
            create_add_mknob_to_fknob(action_before_region.new_stmts,
                                      mknob, fknob)
        end
    end
end

@doc """"
Add function knob to the expressions (library function calls) in the AST.
"""
function add_function_knobs_to_calls(
    call_sites :: CallSites
)
    assert(!call_sites.extra.ast_may_change)

    for fknob in call_sites.extra.function_knobs
        expr = call_sites.extra.fknob2expr[fknob]
        push!(expr.args, fknob)

        # It seems that Julia looks up a dictionary's key by value, not by
        # address (It does not have the concept of pointer anyway), and thus
        # expr is considered changed with the new fknob added. So update the
        # mappings between the expr and the fknob:
        call_sites.extra.expr2fknob[expr]  = fknob
        call_sites.extra.fknob2expr[fknob] = expr
    end
end

@doc """
Create statements to propagate matrix info from one matrix knob to another if
the expression is an assignment statement.
"""
function propagate_matrix_info(
    expr        :: Any,
    call_sites :: CallSites
)
    assert(typeof(expr) == Expr)

    bb           = call_sites.extra.bb
    stmt_idx     = call_sites.extra.stmt_idx
    matrix_knobs = call_sites.extra.matrix_knobs
    matrix2mknob = call_sites.extra.matrix2mknob
    if expr.head == :(=) && length(expr.args) == 2
        # TODO: enable propagation for other kinds of assignment like >>=
        if haskey(matrix2mknob, expr.args[1]) && haskey(matrix2mknob, expr.args[2])
            mknob1 = matrix2mknob[expr.args[1]]
            mknob2 = matrix2mknob[expr.args[2]]
            create_propagate_matrix_info(mknob1, mknob2, bb, stmt_idx, call_sites)
        end
    end
    return nothing
end

@doc """ A handler to visit_expressions(). """
propagate_matrix_info(ast, call_sites :: CallSites, top_level_number, is_top_level, read) =
    propagate_matrix_info(ast, call_sites)

@doc """
Discover context-sensitive function calls in the loop region. Insert matrix and
function-specific context info (mknobs and fknobs).
"""
function context_sensitive_transformation(
    actions         :: Vector{Action},
    region          :: Region,
    bb2depth        :: Dict{BasicBlock, Int},
    lambda          :: LambdaInfo,
    symbol_info     :: Sym2TypeMap,
    liveness        :: Liveness,
    cfg             :: CFG
)
    # We need to walk the AST recursively twice: first, match/replace and gather
    # context info. We cannot add context info into the AST yet, because
    # replacement happens recursively: a subtree might be replaced, and then
    # its parent tree is replaced subsequently, and during which the subtree is
    # removed. If we add context info into the subtree, some of its effect (the
    # actions to create them) remains, even though the subtree is gone.
    # Another reason that we should not add context info in this step is that if
    # we do it, a subtree contains fknob, and then matching its parent tree needs
    # a pattern containing a fknob. That is rather problematic.
    # Therefore, we'd better walk the AST twice. After the AST is stablized after
    # the first walk, the second walk can add the context info.
    dprintln(1, 0, "\nContext-sensitive transformation:")
    recursive      = true
    ast_may_change = true
    call_sites = CallSites(Set{CallSite}(), region, lambda, symbol_info, liveness,
                           CS_transformation_patterns,
                           actions, ContextExtra(region.symbol_property, bb2depth, ast_may_change))

#    global trace_call_replacement = true

    # All patterns are non-splittable.
    call_sites.extra.non_splitting_patterns = copy(CS_transformation_patterns)
    visit_expressions(region, cfg, call_sites, recursive, match_replace_an_expr_pattern, match_replace_an_two_statements_pattern)

#    global trace_call_replacement = false

    # The second walk: gather the knobs of the nodes in the AST. The knobs
    # were generated in the first walk. Not all those generated are gathered,
    # since some ast nodes have been replaced, and thus won't be visited this
    # time.
    call_sites.extra.ast_may_change = false
    visit_expressions(region, cfg, call_sites, recursive, gather_knobs)

    # Generate/delete the knobs at the beginning/exit of the loop
    generate_and_delete_knobs(region, call_sites)

    # Now at the beginning of the loop, associate function knobs with their
    # matrix knobs.
    add_matrix_knobs_to_function_knobs(call_sites)

    # Propagate matrix info from one matrix knob to another for assignment statements
    recursive  = false
    visit_expressions(region, cfg, call_sites, recursive, propagate_matrix_info)

    # Add fknobs to function calls in the AST. This changes AST, and thus should
    # be put at the last, because otherwise, the maps in the extra including
    # matrix2mknobs and expr2fknob whose keys are AST nodes, might become
    # out of date.
    add_function_knobs_to_calls(call_sites)

    # Reordering is now part of context-sensitive optimization.
    if reorder_enabled
        reordering(actions, region, symbol_info, liveness, cfg, call_sites)
    end

    #global trace_call_replacement = false
end

@doc """
Discover context-sensitive function calls in the loop regions. Insert matrix and
function-specific context info (mknobs and fknobs) into actions.
"""
function AST_context_sensitive_transformation(
    actions     :: Vector{Action},
    func_region :: FunctionRegion,
    regions     :: Vector{LoopRegion},
    bb2depth    :: Dict{BasicBlock, Int},
    lambda      :: LambdaInfo,
    symbol_info :: Sym2TypeMap,
    liveness    :: Liveness,
    cfg         :: CFG
)
    # Find properties for the function and loop regions.
    func_region.symbol_property = find_properties_of_matrices(func_region, lambda, symbol_info, liveness, cfg)
    for region in regions
        region.symbol_property = find_properties_of_matrices(region, lambda, symbol_info, liveness, cfg)
    end

    if context_sensitive_opt_for_func
        # Context opt for the function
        context_sensitive_transformation(actions, func_region, bb2depth, lambda, symbol_info, liveness, cfg)
    else
        # Context opt for each loop region
       for region in regions
            context_sensitive_transformation(actions, region, bb2depth, lambda, symbol_info, liveness, cfg)
       end
    end

    dprintln(1, 0, "\nContext-sensitive actions to take:")
    dprintln(1, 1, "", actions)

    actions
end

ext: jl

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using REPL
using REPL.LineEdit
# FIXME: Should refactor all REPL related functions into a struct that keeps track
#        of global state (terminal size, current prompt, current module etc).
# FIXME: Find a way to reprint what's currently entered in the REPL after changing
#        the module (or delete it in the buffer).

using Logging: with_logger
using .Progress: JunoProgressLogger

function get_main_mode()
  mode = Base.active_repl.interface.modes[1]
  mode isa LineEdit.Prompt || error("no julia repl mode found")
  mode
end

function isREPL(; before_run_repl = false)
  isdefined(Base, :active_repl) &&
  isdefined(Base.active_repl, :interface) &&
  isdefined(Base.active_repl.interface, :modes) &&
  if before_run_repl
    true
  else
    isdefined(Base.active_repl, :mistate) &&
    isa(Base.active_repl.mistate, REPL.LineEdit.MIState)
  end
end

handle("changeprompt") do prompt
  isREPL() || return

  if !(isempty(prompt))
    changeREPLprompt(prompt)
  end
  nothing
end

handle("changemodule") do data
  isREPL() || return

  @destruct [mod || ""] = data
  if !isempty(mod) && !isdebugging()
    parts = split(mod, '.')
    if length(parts) > 1 && parts[1] == "Main"
      popfirst!(parts)
    end
    changeREPLmodule(mod)
  end
  nothing
end

handle("fullpath") do uri
  return fullREPLpath(uri)
end

handle("validatepath") do uri
  return isfile′(fullREPLpath(uri)[1])
end

const juliaprompt = "julia> "

handle("resetprompt") do linebreak
  isREPL() || return
  linebreak && println()
  changeREPLprompt(juliaprompt)
  nothing
end

const current_prompt = Ref{String}(juliaprompt)

function hideprompt(f)
  isREPL() || return f()

  repl = Base.active_repl
  mistate = repl.mistate
  mode = mistate.current_mode

  buf = String(take!(copy(LineEdit.buffer(mistate))))

  # clear input buffer
  truncate(LineEdit.buffer(mistate), 0)
  LineEdit.refresh_multi_line(mistate)

  print(stdout, "\e[1K\r")
  r = f()

  flush(stdout)
  flush(stderr)
  sleep(0.05)

  pos = @rpc cursorpos()
  pos[1] != 0 && println()

  # restore prompt
  if applicable(LineEdit.write_prompt, stdout, mode)
    LineEdit.write_prompt(stdout, mode)
  elseif mode isa LineEdit.PrefixHistoryPrompt || :parent_prompt in fieldnames(typeof(mode))
    LineEdit.write_prompt(stdout, mode.parent_prompt)
  else
    printstyled(stdout, current_prompt[], color=:green)
  end

  truncate(LineEdit.buffer(mistate), 0)

  # restore input buffer
  LineEdit.edit_insert(LineEdit.buffer(mistate), buf)
  LineEdit.refresh_multi_line(mistate)
  r
end

function changeREPLprompt(prompt; color = :green, write = true)
  if strip(prompt) == strip(current_prompt[])
    return nothing
  end

  main_mode = get_main_mode()
  main_mode.prompt = prompt
  main_mode.prompt_prefix = color isa Symbol ?
    sprint() do io
      printstyled(IOContext(io, :color => true), " ", color=color)
    end |> split |> first |> string :
    color
  if Base.active_repl.mistate isa REPL.LineEdit.MIState &&
     Base.active_repl.mistate.current_mode == main_mode &&
     write
    print(stdout, "\e[1K\r")
    REPL.LineEdit.write_prompt(stdout, main_mode)
  end
  current_prompt[] = prompt
  nothing
end

# basically the same as Base's `display_error`, just with different frames removed
function display_error(io, err, st)
  ind = findfirst(frame -> frame.file == Symbol(@__FILE__) && frame.func == :repleval, st)
  st = st[1:(ind == nothing ? end : ind - 2)]
  printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
  showerror(IOContext(io, :limit => true), err, st)
  println(io)
end

function display_parse_error(io, err::Meta.ParseError)
  printstyled(io, "ERROR: "; bold=true, color=Base.error_color())
  printstyled(io, "syntax: ", err.msg; color=Base.error_color())
  println(io)
end

# don't inline this so we can find it in the stacktrace
@noinline repleval(mod, line) = Core.eval(mod, line)

const inREPL = Ref{Bool}(false)

function evalrepl(mod, line)
  global ans
  try
    lock(evallock)
    msg("working")
    inREPL[] = true
    fixjunodisplays()
    # this is slow:
    errored = false
    with_logger(JunoProgressLogger()) do
      try
         line = Meta.parse(line)
      catch err
        errored = true
        display_parse_error(stderr, err)
      end
      errored && return nothing
      try
        ans = repleval(mod, line)
      catch err
        # #FIXME: This is a bit weird (there shouldn't be any printing done here), but
        # seems to work just fine.

        errored = true
        display_error(stderr, err, stacktrace(catch_backtrace()))
      end
    end
    errored ? nothing : ans
  catch err
    # This is for internal errors only.
    display_error(stderr, err, stacktrace(catch_backtrace()))
  finally
    inREPL[] = false
    unlock(evallock)
    @async begin
      msg("doneWorking")
      msg("updateWorkspace")
    end
    nothing
  end
end

function changeREPLmodule(mod)
  islocked(evallock) && return nothing

  mod = getmodule(mod)

  # change the repl module in the context of both evaluation and completions
  # NOTE: help_mode doesn't allow custom provider yet:
  # https://github.com/JuliaLang/julia/pull/33102 may help with them all
  main_mode = get_main_mode()
  main_mode.on_done = REPL.respond(Base.active_repl, main_mode; pass_empty = false) do line
    if !isempty(line)
      quote
        $evalrepl($mod, $line)
      end
    end
  end
  main_mode.complete = JunoREPLCompletionProvider(mod)
end

function reset_repl_history()
  mode = get_main_mode()
  hp = mode.hist
  replc = mode.complete

  if Base.active_repl.history_file
    try
      hist_path = REPL.find_hist_file()
      mkpath(dirname(hist_path))
      f = open(hist_path, read=true, write=true, create=true)
      atexit(() -> close(f))
      REPL.hist_from_file(hp, f, hist_path)
      REPL.history_reset_state(hp)
    catch e
    end
  end
end

# make sure DisplayHook() is higher than REPLDisplay() in the display stack
function fixdisplayorder()
  if isREPL()
    Media.unsetdisplay(Editor(), Any)
    fixjunodisplays()
  end
end

function fixjunodisplays()
  for d in reverse(Base.Multimedia.displays)
    if d isa JunoDisplay
      popdisplay(JunoDisplay())
    end
  end
  pushdisplay(JunoDisplay())
end

# completions

struct JunoREPLCompletionProvider <: REPL.CompletionProvider
  mod::Module
end

function LineEdit.complete_line(c::JunoREPLCompletionProvider, s)
  partial = REPL.beforecursor(s.input_buffer)
  full = LineEdit.input_string(s)

  # module-aware repl backend completions
  comps, range, should_complete = REPLCompletions.completions(full, lastindex(partial), c.mod)
  ret = map(REPLCompletions.completion_text, comps) |> unique!

  return ret, partial[range], should_complete
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, HistogramThresholding
makedocs(sitename="Documentation",
            Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"))
deploydocs(repo = "github.com/zygmuntszpak/HistogramThresholding.jl.git")


text:
################
# Secrets Type #
################

@ghdef mutable struct Secret
    name::Union{String, Nothing}
    created_at::Union{Dates.DateTime, Nothing}
    updated_at::Union{Dates.DateTime, Nothing}
end

Secret(name::AbstractString) = Secret(Dict("name" => name))

namefield(secret::Secret) = secret.name

##################
# PublicKey Type #
##################

mutable struct PublicKey <: GitHubType
    key_id::Union{String, Nothing}
    key::Union{String, Nothing}
end

PublicKey(data::Dict) = json2github(PublicKey, data)

namefield(key::PublicKey) = key.key

function publickey(api::GitHubAPI, repo; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/actions/secrets/public-key"; options...)
    return PublicKey(result)
end

###############
# API Methods #
###############

@api_default function secret(api::GitHubAPI, repo, secret_obj; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/actions/secrets/$(name(secret_obj))"; options...)
    return Secret(result)
end

@api_default function secrets(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/actions/secrets"; options...)
    return map(Secret, results["secrets"]), page_data
end

@api_default function create_secret(api::GitHubAPI, repo, secretname; value, options...)
    key = publickey(api, repo; options...)
    k = SodiumSeal.KeyPair(key.key)
    encrypted = SodiumSeal.seal(base64encode(value), k)
    result = gh_put_json(api, "/repos/$(name(repo))/actions/secrets/$secretname"; params=Dict("encrypted_value"=>encrypted, "key_id"=>key.key_id), options...)
    return nothing
end

@api_default function delete_secret(api::GitHubAPI, repo, secretname; options...)
    gh_delete(api, "/repos/$(name(repo))/actions/secrets/$(name(secretname))"; options...)
    return nothing
end

"""
    secrets(repo; auth::Authorization) -> list::Vector{Secret}, page_data

List the names of secrets for `repo`. Requires that you [`authenticate`](@ref).

Secret values cannot be queried, only the names and the times of creation and most recent updates.
"""
secrets

"""
    create_secret(repo, name; value::String, auth::Authorization)

Create GitHub Action secret `name` with value `value`. Requires that you [`authenticate`](@ref).

`value` should be supplied as plain-text and will be encrypted via [libsodium](https://doc.libsodium.org/) for transmission.
"""
create_secret


text:
module UniversalConstants
export getUniversalConstant,getUniversalConstantList

UniversalConstantDict = Dict{String,Real}("standardAccelerationOfGravity" => 9.80665,
"avogadroConstant" => 6.02214199e+23,
"boltzmannConstant" => 1.3806503e-23,
"molarGasConstant" => 8.314472,
"plankConstant" =>6.626070150e-34,
"speedOfLightInVacuum"=> 299792458,
"standardAccelerationOfGravity"=> 9.80665,
"idealGasStateReferencePressure"=> 101325)

getUniversalConstant(str::String) = UniversalConstantDict[str]
getUniversalConstantList()= collect(keys(UniversalConstantDict))
end

__precompile__()


text:
# Define the dynamics model of the model.
struct UnicycleGame{N,M,P,SVu,SVx,SVz} <: AbstractGameModel
	# s = [x,y,θ,v]
    n::Int  # Number of states
    m::Int  # Number of controls
	p::Int  # Number of players
	ni::Vector{Int}  # Number of states for each player
	mi::Vector{Int}  # Number of controls for each player
	pu::SVu # Indices of the each player's controls
	px::SVx # Indices of the each player's x and y positions
	pz::SVz # Indices of the each player's states
end

function UnicycleGame(;p::Int=2)
	# p = number of players
	n = 4p
	m = 2p
	pu = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	px = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	pz = [SVector{4,Int}([i + (j-1)*p for j=1:4]) for i=1:p]
	TYPE = typeof.((pu,px,pz))
	ni = 4*ones(Int,p)
	mi = 2*ones(Int,p)
	return UnicycleGame{n,m,p,TYPE...}(n,m,p,ni,mi,pu,px,pz)
end

@generated function RobotDynamics.dynamics(model::UnicycleGame{N,M,P}, x, u) where {N,M,P}
	xd  = [:(cos(x[$i])*x[$i+P]) for i=M+1:M+P]
	yd  = [:(sin(x[$i])*x[$i+P]) for i=M+1:M+P]
	qdd = [:(u[$i]) for i=1:M]
	return :(SVector{$N}($(xd...), $(yd...), $(qdd...)))
end

dim(model::UnicycleGame) = 2


text:
# AbstractCodeGenAction
function takeModule(x::T) where {T<:AbstractCodeGenAction}
    @check_ptrs x
    m = clang_CodeGenAction_takeModule(x)
    m == C_NULL && error("failed to generate IR.")
    return LLVM.Module(m)
end

dispose(x::AbstractCodeGenAction) = clang_CodeGenAction_dispose(x)

# LLVMOnlyAction
LLVMOnlyAction(ctx::Context) = LLVMOnlyAction(create_emit_llvm_only_action(ctx))

function create_emit_llvm_only_action(ctx::Context)
    status = Ref{CXInit_Error}(CXInit_NoError)
    act = clang_EmitLLVMOnlyAction_create(status, ctx.ref)
    @assert status[] == CXInit_NoError
    return act
end


text:
#!/usr/bin/env julia

using VibronicToolkit

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "--conf"
        metavar = "FILE"
        help = "path to config file"
        required = true
    "--beta"
        metavar = "T"
        help = "reciprocal temperature"
        arg_type = Float64
        required = true
    "--uncoupled"
        help = "remove inter-surface and quadratic coupling from system"
        action = :store_true
end
c = parse_args(ARGS, s, as_symbols=true)

sys = read(c[:conf], DenseSystem)
beta = c[:beta]

if c[:uncoupled]
    sys = simplify(diag(sys))
else
    try
        global sys = DiagonalSystem(sys)
    catch ex
        ex isa SurfaceCouplingException || rethrow()
        error("Analytical solution only applies to uncoupled systems")
    end
end

analytical = Analytical(sys, beta)

println("Z: $(analytical.Z)")
println("E: $(analytical.E)")
println("Cv: $(analytical.Cv)")


text:
mutable struct PreFilteringWrapper <: Sampler
    sampler::Sampler
    threshold_strat::Union{ThresholdStrategy, Nothing}
    threshold::Union{Real, Nothing}
    kde_cache::Union{KDECache, Nothing}
    gamma::Union{Symbol, Real}
    function PreFilteringWrapper(sampler::Sampler, kde_cache::Union{KDECache, Nothing}, gamma::Union{Symbol, Real}, threshold_strat::Union{ThresholdStrategy, Nothing}=nothing)
        return new(sampler, threshold_strat, nothing, kde_cache, gamma)
    end
end

function PreFilteringWrapper(sampler::Sampler, kde_cache::KDECache, threshold_strat::Union{ThresholdStrategy, Nothing}=nothing)
    return PreFilteringWrapper(sampler, kde_cache, :precomputed, threshold_strat)
end

function PreFilteringWrapper(sampler::Sampler, gamma::Union{Symbol, Real}, threshold_strat::Union{ThresholdStrategy, Nothing}=nothing)
    return PreFilteringWrapper(sampler, nothing, gamma, threshold_strat)
end

function sample(sampler::PreFilteringWrapper, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    if sampler.threshold_strat === nothing
        return sample(sampler.sampler, data, labels)
    end
    if sampler.kde_cache === nothing || sampler.kde_cache.data != data
        sampler.kde_cache = KDECache(data, sampler.gamma)
    end
    debug(LOGGER, "Prepruning data set with density threshold $(sampler.threshold_strat).")
    sampler.threshold, inlier_mask, outlier_mask = split_masks(sampler.threshold_strat,
                                                               sampler.kde_cache,
                                                               labels)
    sample_mask = inlier_mask
    sub_mask = sample(sampler.sampler, data[:, sample_mask], labels[sample_mask])
    sample_mask[sample_mask] = sub_mask
    return sample_mask
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

mutable struct HParamInfosList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HParamInfosList(; kwargs...)
        obj = new(meta(HParamInfosList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct HParamInfosList
const __meta_HParamInfosList = Ref{ProtoMeta}()
function meta(::Type{HParamInfosList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HParamInfosList)
            __meta_HParamInfosList[] = target = ProtoMeta(HParamInfosList)
            allflds = Pair{Symbol,Union{Type,String}}[:hparam_infos => Base.Vector{HParamInfo}]
            meta(target, HParamInfosList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_HParamInfosList[]
    end
end
function Base.getproperty(obj::HParamInfosList, name::Symbol)
    if name === :hparam_infos
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{HParamInfo}
    else
        getfield(obj, name)
    end
end

mutable struct MetricInfosList <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function MetricInfosList(; kwargs...)
        obj = new(meta(MetricInfosList), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct MetricInfosList
const __meta_MetricInfosList = Ref{ProtoMeta}()
function meta(::Type{MetricInfosList})
    ProtoBuf.metalock() do
        if !isassigned(__meta_MetricInfosList)
            __meta_MetricInfosList[] = target = ProtoMeta(MetricInfosList)
            allflds = Pair{Symbol,Union{Type,String}}[:metric_infos => Base.Vector{MetricInfo}]
            meta(target, MetricInfosList, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_MetricInfosList[]
    end
end
function Base.getproperty(obj::MetricInfosList, name::Symbol)
    if name === :metric_infos
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{MetricInfo}
    else
        getfield(obj, name)
    end
end

mutable struct HParams_HparamsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HParams_HparamsEntry(; kwargs...)
        obj = new(meta(HParams_HparamsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct HParams_HparamsEntry (mapentry)
const __meta_HParams_HparamsEntry = Ref{ProtoMeta}()
function meta(::Type{HParams_HparamsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HParams_HparamsEntry)
            __meta_HParams_HparamsEntry[] = target = ProtoMeta(HParams_HparamsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => ProtoBuf.google.protobuf.Value]
            meta(target, HParams_HparamsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_HParams_HparamsEntry[]
    end
end
function Base.getproperty(obj::HParams_HparamsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::ProtoBuf.google.protobuf.Value
    else
        getfield(obj, name)
    end
end

mutable struct HParams <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HParams(; kwargs...)
        obj = new(meta(HParams), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct HParams
const __meta_HParams = Ref{ProtoMeta}()
function meta(::Type{HParams})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HParams)
            __meta_HParams[] = target = ProtoMeta(HParams)
            allflds = Pair{Symbol,Union{Type,String}}[:hparams => Base.Dict{AbstractString,ProtoBuf.google.protobuf.Value}]
            meta(target, HParams, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_HParams[]
    end
end
function Base.getproperty(obj::HParams, name::Symbol)
    if name === :hparams
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,ProtoBuf.google.protobuf.Value}
    else
        getfield(obj, name)
    end
end

export HParamInfosList, MetricInfosList, HParams_HparamsEntry, HParams
# mapentries: "HParams_HparamsEntry" => ("AbstractString", "ProtoBuf.google.protobuf.Value")


text:
# We wrap the discrete Fourier and cosine transforms in an operator.
# For Float64 and Complex{Float64} we can use FFTW. The plans that FFTW computes
# support multiplication, and hence we can store them in a MultiplicationOperator.
#
# For other types (like BigFloat) we have to resort to a different implementation.
# We make specific operator types that call fft and dct from FastTransforms.jl
# (dct moved there from this code).
#
# This situation will improve once the pure-julia implementation of FFT lands (#6193).

# The types T for which we can use FFTW's fft routines (with element type Complex{T})
FFTW_TYPES = Union{Float32,Float64}

# The types for which we can use FFTW's dct routines
DCT_FFTW_TYPES = Union{Float32,Float64,Complex{Float32},Complex{Float64}}

#############################
# The Fast Fourier transform
#############################

inverse_fourier_operator(src, dest; options...) =
    inverse_fourier_operator(src, dest, operatoreltype(src,dest); options...)
forward_fourier_operator(src, dest; options...) =
    forward_fourier_operator(src, dest, operatoreltype(src,dest); options...)

# For any type but the FFTW types
inverse_fourier_operator(src, dest, ::Type{Complex{T}}; options...) where {T} =
    length(dest) * FunctionOperator{T}(src, dest, ifft)
forward_fourier_operator(src, dest, ::Type{Complex{T}}; options...) where {T} =
    inv(inverse_fourier_operator(dest, src, Complex{T}; options...))

inverse_fourier_operator(src, dest, ::Type{Complex{T}};
            fftwflags = FFTW.MEASURE, options...) where {T <: FFTW_TYPES} =
    bfftw_operator(src, dest, Complex{T}; options...)

function forward_fourier_operator(src, dest, ::Type{Complex{T}};
            fftwflags = FFTW.MEASURE, options...) where {T <: FFTW_TYPES}
    t_op = fftw_operator(src, dest, Complex{T}; options...)
    1/convert(T, length(src)) * t_op
end

function bfftw_operator(src, dest, ::Type{T}; fftwflags = FFTW.MEASURE, options...) where {T}
    dims = 1:dimension(src)
    plan = plan_bfft!(zeros(T, dest), dims; flags = fftwflags)
    MultiplicationOperator{T}(src, dest, plan; inplace = true)
end

function fftw_operator(src, dest, ::Type{T}; fftwflags = FFTW.MEASURE, options...) where {T}
    dims = 1:dimension(src)
    plan = plan_fft!(zeros(T, dest), dims; flags = fftwflags)
    MultiplicationOperator{T}(src, dest, plan; inplace = true)
end

# We have to know the precise type of the FFT plans in order to intercept calls to
# adjoint and inv.
const FFTPLAN{T,N} = FFTW.cFFTWPlan{T,-1,true,N,UnitRange{Int64}} where {T,N}
const IFFTPLAN{T,N} = FFTW.cFFTWPlan{T,1,true,N,UnitRange{Int64}} where {T,N}

# We define adjoints ourselves, since FFTW doesn't
adjoint_multiplication(op::MultiplicationOperator, object::FFTPLAN) =
    bfftw_operator(dest(op), src(op), eltype(object); fftwflags=object.flags)

adjoint_multiplication(op::MultiplicationOperator, object::IFFTPLAN) =
    fftw_operator(dest(op), src(op), eltype(object); fftwflags=object.flags)

# Explicitly return an inverse in order to avoid creating a scaled FFT plan.
inv_multiplication(op::MultiplicationOperator, object::FFTPLAN) =
    1/convert(eltype(object), length(src(op))) * bfftw_operator(dest(op), src(op), eltype(object); fftwflags=object.flags)

inv_multiplication(op::MultiplicationOperator, object::IFFTPLAN) =
    1/convert(eltype(object), length(src(op))) * fftw_operator(dest(op), src(op), eltype(object); fftwflags=object.flags)

adjoint_function(op::FunctionOperator, fun::typeof(fft)) =
    FunctionOperator(dest(op), src(op), ifft) * ScalingOperator(dest(op), length(dest(op)))

adjoint_function(op::FunctionOperator, fun::typeof(ifft)) =
    FunctionOperator(dest(op), src(op), fft) * ScalingOperator(dest(op), 1/convert(eltype(op), length(dest(op))))

# The two functions below are used when computing the inv of a
# FunctionOperator that uses fft/ifft:
inv(fun::typeof(fft)) = ifft
inv(fun::typeof(ifft)) = fft

##################################
# The Discrete Cosine transform
##################################

# We invoke the standard type II DCT and less standard type I DCT here.
# The type II DCT is used for the discrete Chebyshev transform from Chebyshev roots
# The type I DCT is used for the discrete Chebyshev transform from Chebyshev extrema
# See Strang's paper for the different versions:
# http://www-math.mit.edu/~gs/papers/dct.pdf

inverse_chebyshev_operator(src, dest; options...) =
    inverse_chebyshev_operator(src, dest, operatoreltype(src,dest); options...)
forward_chebyshev_operator(src, dest; options...) =
    forward_chebyshev_operator(src, dest, operatoreltype(src,dest); options...)

inverse_chebyshevI_operator(src, dest; options...) =
    inverse_chebyshevI_operator(src, dest, operatoreltype(src,dest); options...)
forward_chebyshevI_operator(src, dest; options...) =
    forward_chebyshevI_operator(src, dest, operatoreltype(src,dest); options...)

# The generic routines for the DCTII. They rely on dct and idct being available for the
# types of coefficients.
inverse_chebyshev_operator(src, dest, ::Type{T}; options...) where {T} =
    FunctionOperator{T}(src, dest, idct)
forward_chebyshev_operator(src, dest, ::Type{T}; options...) where {T} =
    FunctionOperator{T}(src, dest, dct)

# The generic routines for the DCTI. They they are not yet available
inverse_chebyshevI_operator(src, dest, ::Type{T}; options...) where {T} =
    error("The inverse DCTI-type transform is not available for type ", eltype(src))
forward_chebyshevI_operator(src, dest, ::Type{T}; options...) where {T} =
    error("The DCTI-type transform is not available for type ", eltype(src))

# For the FFTW routines, we use a MultiplicationOperator with the DCT plans.
inverse_chebyshev_operator(src, dest, ::Type{T};
            fftwflags = FFTW.MEASURE, options...) where {T <: DCT_FFTW_TYPES} =
    idctw_operator(src, dest, T; fftwflags = fftwflags, options...)

forward_chebyshev_operator(src, dest, ::Type{T};
            fftwflags = FFTW.MEASURE, options...) where {T <: DCT_FFTW_TYPES} =
    dctw_operator(src, dest, T; fftwflags = fftwflags, options...)

function idctw_operator(src, dest, ::Type{T}; fftwflags = FFTW.MEASURE, options...) where {T}
    dims = 1:dimension(src)
    plan = plan_idct!(zeros(T, dest), dims; flags = fftwflags)
    MultiplicationOperator{T}(src, dest, plan; inplace = true)
end

function dctw_operator(src, dest, ::Type{T}; fftwflags = FFTW.MEASURE, options...) where {T}
    dims = 1:dimension(src)
    plan = plan_dct!(zeros(T, dest), dims; flags = fftwflags)
    MultiplicationOperator{T}(src, dest, plan; inplace = true)
end

inverse_chebyshevI_operator(src, dest, ::Type{T};
            fftwflags = FFTW.MEASURE, options...) where {T <: DCT_FFTW_TYPES} =
    dctIw_operator(src, dest, T; fftwflags = fftwflags, options...)

forward_chebyshevI_operator(src, dest, ::Type{T};
            fftwflags = FFTW.MEASURE, options...) where {T <: DCT_FFTW_TYPES} =
    dctIw_operator(dest, src, T; fftwflags = fftwflags, options...)

function dctIw_operator(src, dest, ::Type{T}; fftwflags = FFTW.MEASURE, options...) where {T}
    dims = 1:dimension(src)
    plan = FFTW.plan_r2r!(zeros(T, src), FFTW.REDFT00, dims; flags = fftwflags)
    MultiplicationOperator{T}(src, dest, plan; inplace = true)
end

# We have to know the precise type of the DCT plans in order to intercept calls to
# adjoint and inv.
const DCTII_PLAN{T} = FFTW.DCTPlan{T,5,true} where {T}
const INV_DCTII_PLAN{T} = FFTW.DCTPlan{T,4,true} where {T}
const DCTI_PLAN{T} = FFTW.r2rFFTWPlan{T,(3,),true,1} where {T}

const DCTPLANS = Union{DCTII_PLAN,INV_DCTII_PLAN,DCTI_PLAN}

inv_multiplication(op::MultiplicationOperator, object::DCTII_PLAN) =
    idctw_operator(dest(op), src(op), eltype(object))

inv_multiplication(op::MultiplicationOperator, object::INV_DCTII_PLAN) =
    dctw_operator(dest(op), src(op), eltype(object))

inv_multiplication(op::MultiplicationOperator, object::DCTI_PLAN) =
    1/(2*convert(eltype(object), size(op,1))-2) * dctIw_operator(dest(op), src(op), eltype(object))

adjoint_multiplication(op::MultiplicationOperator, object::INV_DCTII_PLAN) =
    inv_multiplication(op, object)

adjoint_multiplication(op::MultiplicationOperator, object::DCTII_PLAN) =
    inv_multiplication(op, object)

function adjoint_multiplication(op::MultiplicationOperator, object::DCTI_PLAN)
    T = eltype(object)
    diag = ones(T, size(op,1))
    diag[1] = 2
    diag[end] = 2
    D1 = DiagonalOperator(dest(op), dest(op), diag)
    D2 = DiagonalOperator(src(op), src(op), inv.(diag))
    D2 * op * D1
end

# The four functions below are used when computing the adjoint or inv of a
# FunctionOperator that uses dct/idct:
inv(fun::typeof(dct)) = idct
inv(fun::typeof(idct)) = dct

adjoint(fun::typeof(dct)) = idct
adjoint(fun::typeof(idct)) = dct


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
println("Finished loading libraries."); flush(stdout)

println("ARGS: ", ARGS)

marker = Symbol(ARGS[1])
println("marker: ", marker); flush(stdout)

K = parse(Int, ARGS[2])
println("K: ", K)

stm = parse(Bool, ARGS[3])
println("skewtmix: ", stm); flush(stdout)

sims = dict_list(Dict(:K => K, :marker => marker, :skewtmix => stm))

# Run simulations in parallel.
println("Running analyses ..."); flush(stdout)
@time run(sims[1])

# Post process
println("Post process ..."); flush(stdout)
@time postprocess(sims[1])

# Print number of small clusters
@time print_number_of_small_clusters(sims[1], p=0.02)


text:
export bissect_secA
function bissect_secA(h :: AbstractNLPModel,
                      nlpstop :: AbstractStopping;
                      t1 = h.meta.x0, t2 = 100.0,
                      verbose :: Bool = false)

    (t0, t1) = trouve_intervalle(h, t1, 2.5)
    γ = 0.8; t = t1; tp = t0; tqnp = t0
    hk = [0.0]; hkm1 = [0.0]; gkm1 = [0.0]; hplus = [0.0];
    iter = 0

    hk   = obj(h, t);    gk   = grad(h, t)
    hkm1 = obj(h, tqnp); gkm1 = grad(h, tqnp)

    OK = update_and_start!(nlpstop, x = t, gx = gk, g0 = copy(gk))

    verbose &&
      @printf(" iter        tp        t         dN         gk          ")
    verbose &&
      @printf(" gplus        \n")
    verbose &&
      @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                iter, tp[1], t[1], 0.0, gk[1], 0.0)

    while !OK
        s = t .- tqnp; y = gk .- gkm1

        Γ = [3.0] .* (gk .+ gkm1) .* s .- [6.0] .* (hk .- hkm1)
        if y .* s .+ Γ < eps.(Float64) .* (s .^ 2)
          yt = y .+ Γ ./ s
        else
          yt = y #.+ Γ ./ s
        end
        dN = vec(-gk .* s ./ yt)

        if  true in ((tp - t) .* dN .> 0.0) &&  true in (dN ./ (tp - t) .< γ)
          tplus = t .+ dN
          hplus = obj(h, tplus)
          gplus = grad(h, tplus)
        else
          tplus = (t .+ tp) ./ 2
          hplus = obj(h, tplus)
          gplus = grad(h, tplus)
        end

        if (true in (t .> tp))
          if (true in (gplus .< 0.0))
            tp = t; tqnp = t; t = tplus
          else
            tqnp = t; t = tplus
          end
        else
          if (true in (gplus .> 0.0))
            tp = t; tqnp = t; t = tplus
          else
            tqnp = t; t = tplus
          end
        end

        # update values
        hkm1 = hk; gkm1 = gk; hk = hplus; gk = gplus
        OK = update_and_stop!(nlpstop, x = tp, gx = gk)
        iter = iter + 1

        verbose &&
          @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                  iter, tp[1], t[1], dN[1], gk[1], gplus[1])
      end

      optimal = nlpstop.meta.optimal
      return optimal, nlpstop
end


text:
using CalibrationErrors

@testset "Two-dimensional example" begin
    skce = BiasedSKCE(UniformScalingKernel(4, SquaredExponentialKernel(2)))

    # only two predictions, i.e., one term in the estimator
    @test calibrationerror(skce, ([1 0; 0 1], [1, 2])) ≈ 0
    @test calibrationerror(skce, ([1 0; 0 1], [1, 1])) ≈ 2
    @test calibrationerror(skce, ([1 0; 0 1], [2, 1])) ≈ 4 - 4 * exp(-4)
    @test calibrationerror(skce, ([1 0; 0 1], [2, 2])) ≈ 2
end