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

# --- begin: set_matrix_property inferface --->

@doc """ interface to explicitly specify matrix property in source code"""
function set_matrix_property(
    pmap    :: Dict{Symbol, Int}
)
end

@doc """ specify that matrix A is the lower/upper part or transpose of B"""
function set_matrix_property(
    A       :: Symbol,
    part_of :: Int,
    B       :: Symbol
)
end

@doc """inspect matrix property at run time"""
function inspect_matrix_property(
    M       :: Symbol
)
end

# --- end: set_matrix_property inferface ---

@doc """
A regions's basic information:
"""
type RegionInfo
    region              :: Region
    name                :: AbstractString
    stmts               :: Vector{Statement}
    bblocks             :: Dict{Int, BasicBlock}

    lambda              :: LambdaInfo
    symbol_info         :: Sym2TypeMap
    liveness            :: Liveness
    cfg                 :: CFG

    constants           :: Set
    single_defs         :: Set

    pattern_match_filter :: Dict{Tuple{Expr, ExprPattern}, Int}

    function RegionInfo(
        region          :: Region,
        lambda          :: LambdaInfo,
        symbol_info     :: Sym2TypeMap,
        liveness        :: Liveness,
        cfg             :: CFG,
    )
        info = new()
        info.region = region
        info.lambda = lambda
        info.symbol_info = symbol_info
        info.liveness = liveness
        info.cfg = cfg

        info.pattern_match_filter = Dict{Tuple{Expr, ExprPattern}, Int}()
        info.bblocks = Dict{Int, BasicBlock}()

        if isa(region, LoopRegion)
            first_bb_idx = sort(collect(region.loop.members))[1]
            info.name = "Loop" * string(first_bb_idx)
            bb_idxs = region.loop.members
            for i in bb_idxs
                info.bblocks[i] = cfg.basic_blocks[i]
            end
        else
            info.name = "Func"
            bb_idxs = keys(cfg.basic_blocks)
            info.bblocks = cfg.basic_blocks
        end

        # collect all statements in this region
        info.stmts = []
        for bb_idx in bb_idxs
            bb = cfg.basic_blocks[bb_idx]
            append!(info.stmts, bb.statements)
        end

        info.constants = find_constant_values(region, liveness, cfg)
        info.single_defs = find_single_defs(region, liveness, cfg)

        return info
    end
end

include("symbolic-analysis.jl")
using .SymbolicAnalysis

@doc """
Describe part (lower, upper, diagonal) of a matrix, of which the structure matters.
"""
type MatrixPropertyValues
    constant_valued         :: AbstractSymbol
    constant_sized          :: AbstractSymbol
    constant_structured     :: AbstractSymbol
    maximal_structured      :: AbstractSymbol
    symmetric_valued        :: AbstractSymbol
    symmetric_structured    :: AbstractSymbol
    structure_only          :: AbstractSymbol
    has_dedicated_memory    :: AbstractSymbol # this array (may be a matrix or vector) has a decidated memory space, so no worry of aliases or memory allocation.
    lower_of                :: AbstractSymbol
    upper_of                :: AbstractSymbol
    transpose_of            :: AbstractSymbol
    proxy                   :: Any

    MatrixPropertyValues() = new(TOP_SYMBOL, TOP_SYMBOL,
                            TOP_SYMBOL, TOP_SYMBOL,
                            TOP_SYMBOL, TOP_SYMBOL,
                            TOP_SYMBOL, TOP_SYMBOL,
                            TOP_SYMBOL, TOP_SYMBOL,
                            TOP_SYMBOL, nothing)
end

@doc "sorter for symbol indexed map"
sort_by_str_key = (x) -> sort(collect(keys(x)), by=v->string(v))

function to_string(s :: AbstractSymbol)
    if s == TOP_SYMBOL
        return "-"
    elseif s == BOTTOM_SYMBOL
        return "b"
    else
        return s.value
    end
end

@doc "print out the content of property proxies"
function dprint_property_proxies(
    pmap    :: Dict
)
    val_or_nothing = x -> isa(x, MiddleSymbol)? x.value : "-"
    dprintln(1, 1, "Sym : CV CSize CStruct MStruct SV SS SO DM Lower Upper Trans")
    for k in sort_by_str_key(pmap)
        v = pmap[k]
        dprintln(1, 1, k, " : ",
                    to_string(v.constant_valued), "  ",
                    to_string(v.constant_sized), "  ",
                    to_string(v.constant_structured), "  ",
                    to_string(v.maximal_structured), "  ",
                    to_string(v.symmetric_valued), "  ",
                    to_string(v.symmetric_structured), "  ",
                    to_string(v.structure_only), "  ",
                    to_string(v.has_dedicated_memory), "  ",
                    to_string(v.lower_of), " ",
                    to_string(v.upper_of), " ",
                    to_string(v.transpose_of)
                    )
    end
end

#
const MATRIX_RELATED_TYPES = [SparseMatrixCSC, SparseMatrix.CHOLMOD.Factor, Factorization]

# Symbol types unimportant to analysis
const SKIP_TYPES = [GlobalRef, Int32, Int64, Float64, Bool, QuoteNode, ASCIIString, Complex]

include("structure-analysis-const-size.jl")
include("structure-analysis-maximal-struct.jl")
include("structure-analysis-transpose.jl")
include("structure-analysis-symm-value.jl")
include("structure-analysis-lower.jl")
include("structure-analysis-upper.jl")
include("structure-analysis-struct-only.jl")

const prop_field_const_map = [
           (SA_CONST_VALUED, :constant_valued),
           (SA_CONST_SIZED, :constant_sized),
           (SA_MAXIMAL_STRUCTURED, :maximal_structured),
           (SA_CONST_STRUCTURED, :constant_structured),
           (SA_SYMM_VALUED, :symmetric_valued),
           (SA_SYMM_STRUCTURED, :symmetric_structured),
           (SA_STRUCTURE_ONLY, :structure_only),
           (SA_HAS_DEDICATED_MEMORY, :has_dedicated_memory)]

@doc """
Collect predefined maxtric property according from set_matrix_property statements
in a loop region.
"""
function find_predefined_properties(
    region_info      :: RegionInfo,
)
    property_proxies = Dict{Sym, MatrixPropertyValues}()

    for stmt in region_info.stmts
        expr = stmt.expr
        if typeof(expr) != Expr || expr.head != :call
            continue
        end

        ast = expr
        m, func_name = resolve_call_names(ast.args)
        if func_name == "set_matrix_property"
            if length(ast.args) == 4 # set lower_of/upper_of
                #dump(ast.args)
                assert(isa(ast.args[2], QuoteNode) && isa(ast.args[4], QuoteNode))
                sym = ast.args[2].value
                psym = ast.args[4].value
                if !haskey(property_proxies, sym)
                    property_proxies[sym] = MatrixPropertyValues()
                end
                if ast.args[3].name == :SA_LOWER_OF
                    property_proxies[sym].lower_of = MiddleSymbol(psym)
                    dprintln(1, 1, "Predef: ", sym, " is lower of ", psym)
                elseif ast.args[3].name == :SA_UPPER_OF
                    property_proxies[sym].upper_of = MiddleSymbol(psym)
                    dprintln(1, 1, "Predef: ", sym, " is upper of ", psym)
                else
                    assert(ast.args[3].name == :SA_TRANSPOSE_OF)
                    property_proxies[sym].transpose_of = MiddleSymbol(psym)
                    dprintln(1, 1, "Predef: ", sym, " is transpose of ", psym)
                end
            else # set other properties
                pmap = eval(ast.args[2])
                #dump(pmap)
                assert(isa(pmap, Dict))
                for (sym, p) in pmap
                    if haskey(property_proxies, sym)
                        sp = property_proxies[sym]
                    else
                        sp = MatrixPropertyValues()
                    end

                    for (prop_const_val, prop_field) in prop_field_const_map
                        if (p & prop_const_val) != 0
                            dprintln(1, 1, "Predef: ", sym, ".", prop_field)
                            if prop_const_val == SA_CONST_SIZED
                                sp.(prop_field) = MiddleSymbol((sym, sym))
                            elseif prop_const_val == SA_MAXIMAL_STRUCTURED
                                sp.(prop_field) = MiddleSymbol((sym))
                            else
                                sp.(prop_field) = MiddleSymbol(:true)
                            end
                        end
                    end

                    #add symbol to property map
                    property_proxies[sym] = sp
                end
            end
            #replace the statement with nothing.
            stmt.expr = :(nothing)
        end
    end

    property_proxies
end

@doc """ Find the properties of all the matrices in the region.
A region is currently defined as a loop region.
"""
function find_properties_of_matrices(
    region          :: Region,
    lambda          :: LambdaInfo,
    symbol_info     :: Sym2TypeMap,
    liveness        :: Liveness,
    cfg             :: CFG,
)
    region_info = RegionInfo(region, lambda, symbol_info, liveness, cfg)
    region_name = region_info.name

    dprintln(1, 0, "\n---- Matrix Structure Property Analysis for " * region_info.name * " -----\n")

    # symbol does not have a constant structure?
    property_proxies = find_predefined_properties(region_info)
    predefined_symbols = keys(property_proxies)

    region.property_proxies = property_proxies

    prop_fields = [:constant_valued, :constant_structured,
                   :constant_sized, :maximal_structured,
                   :symmetric_valued, :symmetric_structured,
                   :structure_only, :has_dedicated_memory,
                   :lower_of, :upper_of, :transpose_of ]

    #inherite positive properties from the parent
    if isa(region, LoopRegion)
        for (sym, p) in region.parent.property_proxies
            if !haskey(property_proxies, sym)
                property_proxies[sym] = MatrixPropertyValues()
            end
            sp = property_proxies[sym]

            for prop in prop_fields
                p_val = getfield(p, prop)
                if isa(p_val, MiddleSymbol)
                    sp.(prop) = MiddleSymbol(copy(p_val.value))
                end
            end
        end
    end

    # do constant value analysis here
    for c in region_info.constants
        if !haskey(property_proxies, c)
            property_proxies[c] = MatrixPropertyValues()
        end
        property_proxies[c].constant_valued = MiddleSymbol(:true)
    end

    dprintln(1, 0, "\nInitial structure proxies:")
    dprint_property_proxies(property_proxies)

    structure_property_passes = (
         StructureAnalysisConstSize.pass_info,
         StructureAnalysisMaximalStruct.pass_info,
         StructureAnalysisSymmValue.pass_info,
         StructureAnalysisTranspose.pass_info,
         StructureAnalysisLower.pass_info,
         StructureAnalysisUpper.pass_info,
         StructureAnalysisStructOnly.pass_info,
    )

    for pass in structure_property_passes
        assert(isa(pass[1], AbstractString)) # name
        assert(isa(pass[2], Tuple)) # rules
        assert(isa(pass[3], Union{Function, Void})) # preprocess
        assert(isa(pass[4], Union{Function, Void})) # postprocess
        assert(isa(pass[5], Function)) # symbolizer

        dprintln(1, 0, "\nBegin ", pass[1], " pass:")

        rules = check_and_gen_transfer_rules(pass[2])
        analyzer = SymbolicAnalyzer(pass[1], rules, pass[5])

        # preprocess
        if pass[3] == nothing
            predefined = []
        else
            predefined = pass[3](region.property_proxies, symbol_info)
        end

        #
        res = run_analyzer(analyzer, region_info, predefined, 1)
        # creat entry for k if k is missing
        for k in keys(res)
            if !haskey(property_proxies, k)
                property_proxies[k] = MatrixPropertyValues()
            end
        end

        # postprocess
        if pass[4] != nothing
            pass[4](res, region.property_proxies, symbol_info)
        end
    end

    dprintln(1, 0, "\nFinal structure proxies:")
    dprint_property_proxies(property_proxies)

    is_matrix_type = (x) ->  any(t -> (haskey(symbol_info, x) && symbol_info[x]<:t), MATRIX_RELATED_TYPES)
    mfilter = (x) -> filter(is_matrix_type, x)
    has_pos_val = (v) -> isa(v, MiddleSymbol)

    dprintln(1, 0, "\n" * region_name * " Matrix structures discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(property_proxies)))

    dprintln(1, 0, "\n" * region_name * " Constant value discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.constant_valued), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Constant size discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.constant_sized), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Constant structures discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.maximal_structured), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Maximal structures discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.maximal_structured), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Value symmetry discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.symmetric_valued), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Structure symmetry discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.symmetric_structured), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Structure only discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.structure_only), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Has_dedicated_memory discovered:")
    dprintln(1, 1, mfilter(sort_by_str_key(filter((k, v) -> has_pos_val(v.has_dedicated_memory), property_proxies))))

    dprintln(1, 0, "\n" * region_name * " Upper/Lower matrix discovered:")
    for k in sort_by_str_key(property_proxies)
        v = property_proxies[k]
        if has_pos_val(v.lower_of)
            dprintln(1, 1, k, " is lower of ", v.lower_of.value)
        end
        if has_pos_val(v.upper_of)
            dprintln(1, 1, k, " is upper of ", v.upper_of.value)
        end
    end

    dprintln(1, 0, "\n" * region_name * " Transpose matrix discovered:")
    for k in sort_by_str_key(property_proxies)
        v = property_proxies[k]
        if has_pos_val(v.transpose_of)
            dprintln(1, 1, k, " is transpose of ", v.transpose_of.value)
        end
    end

    # Merge with structure info
    symbols = union(region_info.single_defs, collect(keys(property_proxies)))

    matrix_properties = Symexpr2PropertiesMap()

    is_middle_symbol = (v) -> isa(v, MiddleSymbol)
    get_sym_or_nothing = (v) -> isa(v, MiddleSymbol) ? v.value : nothing

    for s in symbols
        ## These symbols are not necessary related with matrices. But
        ## some symbols may be of subtypes of Array, some other may be
        ## not (like CHOLMOD.Factor, which is not a matrix, but is realted
        ## with matrix). So we'd better not to filter out any symbol.
        matrix_properties[s] = MatrixProperties()
        if in(s, region_info.single_defs)
            matrix_properties[s].is_single_def = true
        end

        if haskey(property_proxies, s)
            p = property_proxies[s]
            matrix_properties[s].constant_valued = is_middle_symbol(p.constant_valued)
            matrix_properties[s].constant_structured = is_middle_symbol(p.maximal_structured)
            #matrix_properties[s].constant_sized = is_middle_symbol(p.constant_sized)
            matrix_properties[s].is_symmetric = is_middle_symbol(p.symmetric_valued)
            matrix_properties[s].is_structure_symmetric = is_middle_symbol(p.symmetric_structured)
            matrix_properties[s].is_structure_only = is_middle_symbol(p.structure_only)
            matrix_properties[s].has_dedicated_memory = is_middle_symbol(p.has_dedicated_memory)
            matrix_properties[s].lower_of = get_sym_or_nothing(p.lower_of)
            matrix_properties[s].upper_of = get_sym_or_nothing(p.upper_of)
            matrix_properties[s].transpose_of = get_sym_or_nothing(p.transpose_of)
        end
    end

    region.symbol_property = matrix_properties

    matrix_properties
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
### ODEProblem from AbstractReactionNetwork ###
function DiffEqBase.ODEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union{AbstractArray, Number}, args...; kwargs...)
    isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing ODEProblems")
    ODEProblem(rn.odefun, u0::Union{AbstractArray, Number}, args...; kwargs...)
end

### SDEProblem from AbstractReactionNetwork ###
function DiffEqBase.SDEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union{AbstractArray, Number}, args...; kwargs...)
    isa(rn.sdefun, SDEFunction) || error("Call addsdes! before constructing SDEProblems")
    SDEProblem(rn.sdefun, rn.g::Function, u0, args...;noise_rate_prototype=rn.p_matrix, kwargs...)
end

### DiscreteProblem, passing through syms
function DiffEqBase.DiscreteProblem(rn::DiffEqBase.AbstractReactionNetwork, u0, tspan::Tuple, p=nothing; kwargs...)
    f = (du,u,p,t) -> du.=u    # identity function to make syms works
    df = DiscreteFunction(f,syms=rn.syms)
    DiscreteProblem(df, u0, tspan, p; kwargs...)
end

### JumpProblem ###
function build_jump_problem(prob, aggregator, rn, kwargs...)
    if typeof(prob)<:DiscreteProblem && any(x->typeof(x) <: VariableRateJump, rn.jumps)
        error("When using time dependant reaction rates a DiscreteProblem should not be used (try an ODEProblem). Also, use a continious solver.")
    end

    # get a JumpSet of the possible jumps
    jset = network_to_jumpset(rn, prob.p)

    # construct map from species index to indices of reactions that depend on it
    if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator)
        rxidxs_to_jidxs   = rxidxs_to_jidxs_map(rn, get_num_majumps(jset))
        spec_to_jumps_set = spec_to_dep_jumps_map(rn, rxidxs_to_jidxs)
        spec_to_jumps_vec = [sort!(collect(specset)) for specset in spec_to_jumps_set]
        jump_to_specs_vec = jump_to_dep_specs_map(rn, rxidxs_to_jidxs)
    else
        rxidxs_to_jidxs   = nothing
        spec_to_jumps_set = nothing
        spec_to_jumps_vec = nothing
        jump_to_specs_vec = nothing
    end

    # construct reaction dependency graph
    if needs_depgraph(aggregator)
        dep_graph = depgraph_from_network(rn, jset, rxidxs_to_jidxs, spec_to_jumps_set)
    else
        dep_graph = nothing
    end

    JumpProblem(prob, aggregator, jset; dep_graph=dep_graph,
                                        vartojumps_map=spec_to_jumps_vec,
                                        jumptovars_map=jump_to_specs_vec,
                                        kwargs...)
end

### JumpProblem from AbstractReactionNetwork ###
function DiffEqJump.JumpProblem(prob, aggregator, rn::DiffEqBase.AbstractReactionNetwork; kwargs...)
    (rn.jumps !== nothing) || error("Call addjumps! before constructing JumpProblems")
    build_jump_problem(prob, aggregator, rn, kwargs...)
end

### SteadyStateProblems from AbstractReactionNetwork ###
function DiffEqBase.SteadyStateProblem(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...)
    isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing SteadyStateProblems")
    SteadyStateProblem(rn.odefun, args...; kwargs...)
end

function DiffEqBase.SteadyStateProblem{isinplace}(rn::DiffEqBase.AbstractReactionNetwork, args...; kwargs...) where isinplace
    isa(rn.odefun, ODEFunction) || error("Call addodes! before constructing SteadyStateProblems")
    SteadyStateProblem{isinplace}(rn.odefun, args...; kwargs...)
end

ext: jl


text:
import ..Common.test_comparison

function run_test_comparison()
    test_comparison(Atom(atom_name="H1"), Atom(atom_name="H1"), true, "Compare :atom_name")
    test_comparison(Atom(atom_name="H1"), Atom(atom_name="H2"), false, "Compare :atom_name")
    test_comparison(Atom(coordinate=[1., 0., 0.]), Atom(coordinate=[1., 0., 0.]), true, "Compare :coordinate")
    test_comparison(Atom(coordinate=[1., 0., 0.001]), Atom(coordinate=[1., 0., 0.]), false, "Compare :coordinate")
end


text:
module test_bukdu_assoc

using Test
using Bukdu: Assoc

assoc = Assoc("A" => "B")
@test !isempty(assoc)

@test pairs(["A" => "B"]) == pairs(assoc)

params = Assoc("x" => "2")
@test get(params, :x, 0) == 2
@test get(params, :y, 3) == 3
@test !haskey(params, :a)
@test params.a == ""

params[:x] = "3"
@test params.x == "3"
params[:y] = "5"
@test params.y == "5"
@test length(params) == 2
@test keys(params) == ["x", "y"]
@test values(params) == ["3", "5"]

empty!(params)
@test !haskey(params, "x")

buf = IOBuffer()
show(buf, MIME"text/plain"(), params)
@test String(take!(buf)) == "Bukdu.Assoc()\n"

end # module test_bukdu_assoc

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
include("tools.jl")
import .DATools: FixedSizeBinaryMaxHeap, length, pop!, push!, top, getindex

function create_resp_ranks(num_props::Int, num_resps::Int, resp_prefs::Vector{Vector{Int}})
    resp_ranks = zeros(UInt16, num_props, num_resps)
    for j in 1:num_resps
        for (r, i) in enumerate(resp_prefs[j])
            @inbounds resp_ranks[i, j] = r
        end
    end
    return resp_ranks
end

function convert_into_matched(num_props::Int,
                              resp_prefs::Vector{Vector{Int}},
                              resp_matched_ranks::Vector{FixedSizeBinaryMaxHeap},
                              caps::Vector{Int})
    prop_matched = zeros(Int, num_props)
    resp_matched = Vector{Int}(sum(caps))
    ctr = 1
    for j in 1:length(resp_matched_ranks)
        for k in 1:caps[j]
            if k <= length(resp_matched_ranks[j])
                p = resp_prefs[j][resp_matched_ranks[j][k]]
                prop_matched[p] = j
                resp_matched[ctr] = p
            else
                resp_matched[ctr] = 0
            end
            ctr += 1
        end
    end
    return prop_matched, resp_matched
end

function deferred_acceptance(prop_prefs::Vector{Vector{Int}}, resp_prefs::Vector{Vector{Int}})
    prop_matched, resp_matched, _ = deferred_acceptance(prop_prefs, resp_prefs, ones(Int, length(resp_prefs)))
    return prop_matched, resp_matched
end

function deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                             resp_prefs::Vector{Vector{Int}},
                             caps::Vector{Int})
    # Set Up
    num_props = length(prop_prefs)
    num_resps = length(resp_prefs)

    resp_ranks = create_resp_ranks(num_props, num_resps, resp_prefs)

    prop_ptrs = ones(UInt16, num_props)
    resp_matched_ranks = FixedSizeBinaryMaxHeap[FixedSizeBinaryMaxHeap(caps[j]) for j in 1:num_resps]
    prop_unmatched = Vector{Int}(num_props)
    for i in 1:num_props
        prop_unmatched[i] = i
    end
    remaining::Int = num_props
    p::Int = 0
    p_rank::UInt16 = 0
    least_p::Int = 0
    least_p_rank::UInt16 = 0
    r::Int = 0

    # Main Loop
    while remaining > 0
        p = prop_unmatched[remaining]
        if prop_ptrs[p] > length(prop_prefs[p])#p no longer wants to find a partner
            remaining -= 1
            continue
        end
        r = prop_prefs[p][prop_ptrs[p]]
        resp_matched_rank = resp_matched_ranks[r]
        least_p_rank = top(resp_matched_rank)
        p_rank = resp_ranks[p, r]
        if p_rank == 0
            prop_ptrs[p] += 1
        elseif length(resp_matched_rank) < caps[r]
            push!(resp_matched_rank, p_rank)
            remaining -= 1
        elseif least_p_rank > p_rank
            pop!(resp_matched_rank)
            push!(resp_matched_rank, p_rank)
            least_p = resp_prefs[r][least_p_rank]
            prop_unmatched[remaining] = least_p
            prop_ptrs[least_p] += 1
        else
            prop_ptrs[p] += 1
        end
    end

    prop_matched, resp_matched = convert_into_matched(num_props, resp_prefs, resp_matched_ranks, caps)
    indptr = Vector{Int}(num_resps+1)
    indptr[1] = 1
    for j in 1:num_resps
        indptr[j+1] = indptr[j] + caps[j]
    end
    return prop_matched, resp_matched, indptr
end

function deferred_acceptance(prop_prefs::Vector{Vector{Int}},
                             resp_prefs::Vector{Vector{Int}},
                             resp_caps::Vector{Int},
                             prop_caps::Vector{Int})
    v_num_props = sum(prop_caps)# number of virtual props
    num_props = length(prop_prefs)
    refs = Vector{Vector{Int}}(v_num_props)
    refs_rev = Vector{Int}(v_num_props)
    prop_indptr = Vector{Int}(num_props+1)

    prop_indptr[1] = 1
    for i in 1:num_props
        prop_indptr[i+1] = prop_indptr[i] + prop_caps[i]
    end

    for (i, prop_cap) in enumerate(prop_caps)
        refs[i] = collect(prop_indptr[i]:prop_indptr[i+1]-1)
        refs_rev[prop_indptr[i]:prop_indptr[i+1]-1] = i
    end

    v_prop_prefs = Vector{Vector{Int}}(v_num_props)
    for (i, prop_pref) in enumerate(prop_prefs)
        for v_i in refs[i]
            v_prop_prefs[v_i] = prop_pref
        end
    end

    v_resp_prefs = deepcopy(resp_prefs)
    for j in 1:length(resp_prefs)
        for ind in length(resp_prefs[j]):-1:1
            i = resp_prefs[j][ind]
            v_is = refs[i]
            splice!(v_resp_prefs[j], ind, v_is)
        end
    end

    prop_matched, resp_matched, resp_indptr = deferred_acceptance(v_prop_prefs, v_resp_prefs, resp_caps)

    for k in 1:length(resp_matched)
        if resp_matched[k] != 0
            resp_matched[k] = refs_rev[resp_matched[k]]
        end
    end

    return prop_matched, resp_matched, resp_indptr, prop_indptr
end

function generate_random_prefs{T <: Integer}(num_props::T, num_resps::T; max_prop_pref::Int=num_resps, max_resp_pref::Int=num_props)
    prop_prefs = Vector{Vector{Int}}(num_props)
    resp_prefs = Vector{Vector{Int}}(num_resps)
    for i in 1:num_props
        prop_prefs[i] = collect(Iterators.take(shuffle(collect(1:num_resps)), max_prop_pref))
    end
    for j in 1:num_resps
        resp_prefs[j] = collect(Iterators.take(shuffle(collect(1:num_props)), max_resp_pref))
    end
    return prop_prefs, resp_prefs
end


text:
_fill_tuples_expr(inner::Function, SZ::Tuple{Int}, inds...) =
    :(tuple($(ntuple(i->inner(i, inds...), SZ[1])...)))
_fill_tuples_expr{N}(inner::Function, SZ::NTuple{N, Int}, inds...) =
    :(tuple($(ntuple(i->_fill_tuples_expr(inner, SZ[1:end-1], i, inds...),SZ[end])...)))
fill_tuples_expr(inner::Function, SZ::Tuple) = _fill_tuples_expr(inner, SZ)

_fill_tuples(inner, originalSZ, SZ::Tuple{Int}, inds::Int...) =
    ntuple(i->inner(SZ, i, inds...), Val{SZ[1]})
_fill_tuples{N}(inner, originalSZ, SZ::NTuple{N, Int}, inds::Int...) =
    ntuple(i->_fill_tuples(inner, originalSZ, SZ[1:end-1], i, inds...), Val{SZ[end]})
fill_tuples{N}(inner, SZ::NTuple{N, Int}) = _fill_tuples(inner, SZ, SZ)

#=
Constructors for homogenous non tuple arguments
Is unrolled for the first 4, since a::T... leads to slow performance
=#
call{FSA <: FixedArray, T}(::Type{FSA}, a::T) = FSA(NTuple{1,T}((a,)))
call{FSA <: FixedArray, T}(::Type{FSA}, a::T, b::T) = FSA(NTuple{2,T}((a,b)))
call{FSA <: FixedArray, T}(::Type{FSA}, a::T, b::T, c::T) = FSA(NTuple{3,T}((a,b,c)))
call{FSA <: FixedArray, T}(::Type{FSA}, a::T, b::T, c::T, d::T) = FSA(NTuple{4,T}((a,b,c,d)))
call{FSA <: FixedArray, T}(::Type{FSA}, a::T...) = FSA(a)

immutable ParseFunctor{T, S <: AbstractString} <: Func{1}
    t::Type{T}
    a::Vector{S}
end
call{T}(pf::ParseFunctor{T}, i::Int) = parse(T, pf.a[i])
call(pf::ParseFunctor{Void}, i::Int) = parse(pf.a[i])

"""
Constructs a fixedsize array from a Base.Array
"""
@generated function call{FSA <: FixedArray, T <: Array}(::Type{FSA}, a::T)
    if eltype(a) <: AbstractString
        ElType = eltype_or(FSA, Void)
        return :(map(ParseFunctor($ElType, a), FSA)) # can't be defined in another method as it leads to lots of ambigouity with the default constructor
    end
    SZ     = size_or(FSA, :(size(a)))
    ElType = eltype_or(FSA, eltype(a))
    if isa(SZ, Expr)
        expr = :($FSA(fill_tuples((sz, i...)->$ElType(a[i...]), $SZ)))
    else
        tupexpr = fill_tuples_expr((i,inds...) -> :($ElType(a[$i, $(inds...)])), SZ)
        expr = :($FSA($tupexpr))
    end
    quote
        $SZ != size(a) && throw(DimensionMismatch("size of $FSA is not fitting array $(typeof(a)), with size: $(size(a))"))
        $expr
    end
end

"""
Constructor for singular arguments.
Can be a tuple, is not declared as that, because it will generate ambigouities
and overwrites the default constructor.
"""
@generated function call{FSA <: FixedArray, X}(::Type{FSA}, a::X)
    ND = ndims(FSA)
    if X <: Tuple
        types_svec = a.parameters
        if all(x-> x <: Tuple, types_svec) && ND == 1
            return :(throw(
                DimensionMismatch("tried to construct $FSA from $a. I can't allow that!")
            ))
        end
        orlen = length(types_svec)
        ortyp = promote_type(types_svec...)
    else
        orlen = 1
        ortyp = a
    end
    SZ      = size_or(FSA, (orlen, ntuple(x->1, ND-1)...))
    T       = eltype_or(FSA, ortyp)
    FSAT    = similar(FSA, T, SZ)
    if X <: Tuple
        expr = fill_tuples_expr((inds...)->:($T(a[$(inds[1])])), SZ)
    else
        expr = fill_tuples_expr((inds...)->:($T(a)), SZ)
    end
    return :($FSAT($expr))
end
"""
Constructors for heterogenous multiple arguments.
E.g. 1, 2f0, 4.0
"""
@generated function call{FSA <: FixedArray}(::Type{FSA}, a...)
    SZ     = size_or(FSA, (length(a),))
    ElType = eltype_or(FSA, promote_type(a...))
    all(x-> x <: Tuple, a) && return :( $FSA(a) ) # TODO be smarter about this
    any(x-> x != ElType, a) && return :($FSA(map($ElType, a)))
    return :($FSA(a))
end

"""
Construction from other FixedSizeArrays + X arguments
E.g. Vec4f0(Vec3f0(1), 0)
"""
@generated function call{FSA <: FixedArray, T1 <: FixedArray}(::Type{FSA}, a::T1, b...)
    if isempty(b) # this is the conversion constructor for 2 FSA's
        #easiest way is to just construct from the tuple
        expr = :(FSA(get_tuple(a)))
        if size_or(FSA, nothing) == nothing # no concrete size
            return expr
        else #has a size
            len1 = size(FSA, 1)
            len2 = size(T1, 1)
            if len1 < len2 # we need to shrink
                return :(FSA(get_tuple(a)[1:$len1]))
            elseif len1==len2
                return expr
            else
                return :(throw(DimensionMismatch(
                    "tried to create $FSA from $T1. The latter has too many elements"
                )))
            end
        end
    end
    SZ      = size_or(FSA, nothing)
    ElType  = eltype_or(FSA, eltype(T1))
    a_expr  = :( a )
    if SZ != nothing
        if (prod(SZ) > (length(T1) + length(b)))
            throw(DimensionMismatch("$FSA is too small, can not be constructed from array $a and $b arguments"))
        elseif prod(SZ) < length(T1) && isempty(b) #shrinking constructor, e.g. Vec3(Vec4(...))
            a_expr = :( a[1:$(prod(SZ))] )
        end
    end
    return :( $FSA(a..., b...) )
end

@inline zero{FSA <: FixedArray}(::Type{FSA}) = map(ConstFunctor(zero(eltype(FSA))), FSA)
zero(fsa::FixedArray) = zero(typeof(fsa))
@inline one{FSA <: FixedArray}(::Type{FSA})  = map(ConstFunctor(one(eltype(FSA))), FSA)
@inline eye{FSA <: FixedArray}(::Type{FSA})  = map(EyeFunc{eltype(FSA)}, FSA)
@inline unit{FSA <: FixedVector}(::Type{FSA}, i::Integer) = map(UnitFunctor(i, eltype(FSA)), FSA)

@inline rand{FSA <: FixedArray}(m::MersenneTwister, x::Type{FSA}) = map(MersenneFunctor{eltype(FSA)}(m), FSA)
@inline function rand{FSA <: FixedArray}(x::Type{FSA}, range::Range)
    T = eltype(FSA)
    map(RandFunctor(T(first(range)):T(step(range)):T(last(range))), FSA) # there's no easy way to convert eltypes of ranges (I think)
end
"""
Macro `fsa` helps to create fixed size arrays like Julia arrays.
E.g.
```
@fsa([1 2 3;4 5 6])
@fsa([a,2,3]) # you can also use variables
@fsa([a 2 3])
```
"""
macro fsa(expr)
    if expr.head == :vect
        result = Expr(:call, :Vec, expr.args...)
    elseif expr.head == :hcat
        result = Expr(:call, :Mat, [Expr(:tuple, a) for a in expr.args]...)
    elseif expr.head == :vcat
        if isa(expr.args[1], Expr) && expr.args[1].head == :row
            cols = [Any[] for _=1:length(expr.args[1].args)]
            for row in expr.args
                for (j, a) in enumerate(row.args)
                    push!(cols[j], a)
                end
            end
            result = Expr(:call, :Mat, Expr(:tuple, [Expr(:tuple, col...) for col in cols]...))
        else
            result = Expr(:call, :Vec, expr.args...)
        end
    end
    esc(result)
end


text:
## namedarraytypes.jl.
## (c) 2013--2018 David A. van Leeuwen

## Julia type that implements a drop-in replacement of Array with named dimensions.

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## DT is a tuple of Dicts, characterized by the types of the keys.
## This way NamedArray is dependent on the dictionary type of each dimensions.
## The inner constructor checks for consistency, the values must all be 1:d

using DataStructures: OrderedDict

if ! @isdefined NamedArray

struct Name{T}
    name::T
end

Base.show(io::IO, name::Name) = print(io, name.name)

function checkdict(dict::AbstractDict)
    pairs = Pair[]
    union = Union{}
    n = length(dict)
    covered = falses(n)
    for (key, value) in dict
        if isa(key, Integer)
            key = Name(key)
        end
        union = Union{union, typeof(key)}
        push!(pairs, key => value)
        if isa(value, Integer) && 1 ≤ value ≤ n
            covered[value] = true
        end
    end
    all(covered) || error("Not all target indices are covered")
    return OrderedDict{union, Int}(pairs)
end

mutable struct NamedArray{T,N,AT,DT} <: AbstractArray{T,N}
    array::AT
    dicts::DT
    dimnames::NTuple{N, Any}
    function (::Type{S})(array::AbstractArray{T, N},
                         dicts::NTuple{N, OrderedDict},
                         dimnames::NTuple{N, Any}) where {S<:NamedArray, T, N}
        size(array) == map(length, dicts) || error("Inconsistent dictionary sizes")
        ## dicts = map(dict -> checkdict(dict), dicts)
        new{T, N, typeof(array), typeof(dicts)}(array, dicts, dimnames)
    end
end

const NamedVector{T} = NamedArray{T,1}
const NamedMatrix{T} = NamedArray{T,2}
const NamedVecOrMat{T} = Union{NamedVector{T},NamedMatrix{T}}
const ArrayOrNamed{T,N} = Union{Array{T,N}, NamedArray{T,N,Array}}

if @isdefined RowVector
    NamedRowVector{T,RVT<:AbstractVector} = NamedArray{T,2,RowVector{T,RVT}}
end

end


text:
"""
    solve_multigrad(problem, gradients)

Solve `problem` for multiple magnetic field gradients.
"""
function solve_multigrad(problem, gradients)
    itertimes = zeros(0)
    results = Vector{output_type(problem)}()
    for grad ∈ gradients
        @info "Solving" typeof(problem) grad
        starttime = time()
        push!(results, solve(problem, grad))
        push!(itertimes, time() - starttime)
    end
    results, itertimes
end


text:
using DataStructures

include("tree_helpers.jl")

mutable struct Node
    neighbors::Vector{Node}
    ID::Int64
    Node(ID::Int64) = new(Vector{Node}(), ID)
end

function DFS_recursive(n::Node)
    # Here we are doing something...
    println(n.ID)

    for neighbor in n.neighbors
        DFS_recursive(neighbor)
    end
end

function DFS_stack(n::Node)
    s = Stack{Node}()
    push!(s, n)

    while(length(s) > 0)
        println(top(s).ID)
        temp = pop!(s)
        for neighbor in temp.neighbors
            push!(s, neighbor)
        end
    end
end

function BFS_queue(n::Node)
    q = Queue{Node}()
    enqueue!(q, n)

    while(length(q) > 0)
        println(front(q).ID)
        temp = dequeue!(q)
        for neighbor in temp.neighbors
            enqueue!(q, neighbor)
        end
    end

end


text:
using Test, MirroredArrayViews, Random

const A1 = [1 2
            3 4]

const A2 = [1, 2, 3, 4, 5]

@testset "mirror_vector" begin
    @test MirroredArrayView(A2, 1) == [5, 4, 3, 2, 1]
    @test_throws BoundsError MirroredArrayView(A2, 2)
    @test MirroredArrayView(MirroredArrayView(A2, 1), 1) == A2
end

@testset "mirror_matrix" begin
    @test MirroredArrayView(A1, 1) == [3 4
                                       1 2]

    @test MirroredArrayView(A1, 2) == [2 1
                                       4 3]

    @test MirroredArrayView(A1, 1, 2) == [4 3
                                          2 1]
    @test MirroredArrayView(A1, 2, 1) == MirroredArrayView(A1, 1, 2)

    @test_throws BoundsError MirroredArrayView(A1, 3)
    @test_throws BoundsError MirroredArrayView(A1, 1, 2, 3)

    @test begin
        B = MirroredArrayView(A1, 1)
        C = similar(B)
        C .= B
        MirroredArrayView(C, 1) == A1
    end

    @test begin
        B = MirroredArrayView(A1, 2)
        C = similar(B)
        C .= B
        MirroredArrayView(C, 2) == A1
    end

    @test begin
        B = MirroredArrayView(A1, 1, 2)
        C = similar(B)
        C .= B
        MirroredArrayView(C, 2, 1) == A1
    end
end

@testset "general array mirrors" begin
    # Test that mirroring an array of arbitrary dimension (<= 10 for
    # these tests to keep it reasonable) in some number of dimensions and
    # then mirroring back in each dimension preserves the parent array.

    # Number of base arrays.
    outer_iters = 10

    # Number of combinations of dimensions to test for each array
    inner_iters = 5

    # Limit to array sizes
    maxdims = 10
    maxsize = 5

    for iter = 1:outer_iters
        dim = rand(3:maxdims)
        A = rand(Int, ( (rand(1:maxsize) for i = 1:dim)..., ))
        println("A has size $(size(A))")

        for iiter = 1:inner_iters
            dims = ( shuffle!([i for i = 1:rand(1:dim)])..., )
            println("    Mirroring A along $(dims)")
            B = MirroredArrayView(A, dims...)

            for i = length(dims):-1:1
                unmirror = dims[i:length(dims)]
                mirror = dims[1:i-1]
                println("      Testing mirror(A, $mirror) == mirror(mirror(A, $dims), $unmirror)")
                @test MirroredArrayView(A, mirror...) ==
                      MirroredArrayView(B, unmirror...)
            end
        end
    end
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using GPUifyLoops, Cthulhu, CuArrays, CUDAnative

function kernel!(::Val{nbcol}, A, B) where {nbcol}
    @inbounds @loop for i in (1:size(A,1);
                              (blockIdx().x-1)*blockDim().x + threadIdx().x)
        a = zero(eltype(B))
        for j = 1:nbcol
            a += B[i, j]
        end
        A[i] = a
    end
    nothing
end

@eval function kernel!(A::CuArray, B::CuArray)
  threads = 1024
  blocks = ceil(Int, size(A,1)/threads)
  @launch(CUDA(), threads=threads, blocks=blocks, kernel!(Val(size(B,2)), A, B))
end

b = rand(Float32, 10^8, 10)
a = similar(b, size(b,1))
cb = CuArray(a)
ca = similar(cb, size(b,1))

@static if isdefined(Base, :active_repl)
    # Example of how to decend using Cthulhu
    args = (Val(size(cb,2)), ca, cb)
    cuargs = map(cudaconvert, args)
    tt = map(typeof, cuargs)
    descend(kernel!, Tuple{tt...})
end

# Run the kernel
kernel!(ca, cb)
kernel!(ca, cb)
kernel!(ca, cb)
kernel!(ca, cb)

# Dump intermediate files into the tmp directory
@device_code dir="tmp" kernel!(ca, cb)


text:
using NashGames
using Test

@testset "NashGames.jl" begin
  # Write your tests here.
end


text:
function freqTwo(a::Array{Int, 1}, x, y, left=1, right=length(a))
	if right == left
		return 0
	end
	if right-left == 1
		if a[left] == x && a[right] == y
			return 1
		end
	end
	mid = (left+right) ÷ 2
	increment = 0
	if a[mid] == x && a[mid+1] == y
		increment = 1
	end
	return increment +
		   freqTwo(a, x, y, left, mid) +
		   freqTwo(a, x, y, mid+1, right)
end


text:
# Julia wrapper for header: /usr/include/scip/cons_bivariate.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrBivariate(scip)
    ccall((:SCIPincludeConshdlrBivariate, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsBivariate(scip, cons, name, f, convextype, z, zcoef, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsBivariate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_EXPRTREE}, SCIP_BIVAR_CONVEXITY, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, f, convextype, z, zcoef, lhs, rhs, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicBivariate(scip, cons, name, f, convextype, z, zcoef, lhs, rhs)
    ccall((:SCIPcreateConsBasicBivariate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_EXPRTREE}, SCIP_BIVAR_CONVEXITY, Ptr{SCIP_VAR}, Cdouble, Cdouble, Cdouble), scip, cons, name, f, convextype, z, zcoef, lhs, rhs)
end

function SCIPgetLinearVarBivariate(scip, cons)
    ccall((:SCIPgetLinearVarBivariate, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLinearCoefBivariate(scip, cons)
    ccall((:SCIPgetLinearCoefBivariate, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetExprtreeBivariate(scip, cons)
    ccall((:SCIPgetExprtreeBivariate, libscip), Ptr{SCIP_EXPRTREE}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetLhsBivariate(scip, cons)
    ccall((:SCIPgetLhsBivariate, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRhsBivariate(scip, cons)
    ccall((:SCIPgetRhsBivariate, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x = -2π:0.001:2π
linCol = ["#000000", "#E69F00", "#56B4E9", "#009E73","#F0E442"]
@gp x cos.(x) "w l t 'dt 1' lw 2 dt 1 lc '$(linCol[1])'"
@gp :- "set grid" "set auto fix"
@gp :- "set offsets graph .05, graph .05, graph .05, graph .05"
@gp :- ylabel="Y label" xlab="X label"
@gp :- "set key b r font ',10' title 'dashtype' opaque"
for (idx,ϕ) in enumerate(0.3:0.3:1.3)
    @gp :- x cos.(x .+ 2ϕ) "w l t 'dt $(idx+1)' lw 2 dt $(idx+1) lc '$(linCol[idx+1])'"
end
save(term="pngcairo font 'Consolas, 12' size 600,400", output="./code/plt2_ex6.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# Contiguous on-device arrays (host side representation)

import Base: length, size, copy!, unsafe_convert, pointer, Array

export
    CuArray, free

type CuArray{T,N} <: AbstractArray{T,N}
    ptr::DevicePtr{T}
    shape::NTuple{N,Int}
    len::Int

    function CuArray(shape::NTuple{N,Int})
        if !isbits(T)
            # non-isbits types results in an array with references to CPU objects
            throw(ArgumentError("CuArray with non-bit element type not supported"))
        elseif (sizeof(T) == 0)
            throw(ArgumentError("CuArray with zero-sized element types does not make sense"))
        end
        len = prod(shape)
        ptr = cualloc(T, len)
        new(ptr, shape, len)
    end

    function CuArray(shape::NTuple{N,Int}, ptr::DevicePtr{T})
        len = prod(shape)
        new(ptr, shape, len)
    end
end

(::Type{CuArray{T}}){T,N}(shape::NTuple{N,Int}) = CuArray{T,N}(shape)
(::Type{CuArray{T}}){T}(len::Int)               = CuArray{T,1}((len,))

(::Type{CuArray{T}}){T,N}(shape::NTuple{N,Int}, p::DevicePtr{T}) = CuArray{T,N}(shape, p)
(::Type{CuArray{T}}){T}(len::Int, p::DevicePtr{T})               = CuArray{T,1}((len,), p)

# deprecated
CuArray{T,N}(::Type{T}, shape::NTuple{N,Int}) = CuArray{T,N}(shape)
CuArray{T}(::Type{T}, len::Int)               = CuArray{T,1}((len,))

unsafe_convert{T,N}(::Type{DevicePtr{T}}, a::CuArray{T,N}) = a.ptr
pointer{T}(x::CuArray{T}) = unsafe_convert(DevicePtr{T}, x)

length(g::CuArray) = g.len
size(g::CuArray) = g.shape

"Free GPU memory allocated to the pointer"
function free(g::CuArray)
    if !isnull(g.ptr)
        free(g.ptr)
        g.ptr = DevicePtr{eltype(g.ptr)}()
    end
end

"Copy an array from device to host in place"
function copy!{T}(dst::Array{T}, src::CuArray{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @apicall(:cuMemcpyDtoH, (Ptr{Void}, Ptr{Void}, Csize_t),
                            pointer(dst), src.ptr.inner, nbytes)
    return dst
end

"Copy an array from host to device in place"
function copy!{T}(dst::CuArray{T}, src::Array{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @apicall(:cuMemcpyHtoD, (Ptr{Void}, Ptr{Void}, Csize_t),
                            dst.ptr.inner, pointer(src), nbytes)
    return dst
end

"Transfer an array from host to device, returning a pointer on the device"
CuArray{T,N}(a::Array{T,N}) = copy!(CuArray(T, size(a)), a)

"Transfer an array on the device to host"
Array{T}(g::CuArray{T}) = copy!(Array(T, size(g)), g)


text:
@testset "371.sum-of-two-integers.jl" begin
    @test get_sum371(1, -1) == 0
    @test get_sum371(20, 23) == 43
end


text:
# Transcoding Error
# =================

"""
Container of transcoding error.

An object of this type is used to notify the caller of an exception that
happened inside a transcoding method.  The `error` field is undefined at first
but will be filled when data processing failed. The error should be set by
calling the `setindex!` method (e.g. `error[] = ErrorException("error!")`).
"""
mutable struct Error
    error::Exception

    function Error()
        return new()
    end
end

# Test if an exception is set.
function haserror(error::Error)
    return isdefined(error, :error)
end

function Base.setindex!(error::Error, ex::Exception)
    @assert !haserror(error) "an error is already set"
    error.error = ex
    return error
end

function Base.getindex(error::Error)
    @assert haserror(error) "no error is set"
    return error.error
end


text:
line_defaults() = Dict{Symbol, Any}([
    :labelV => nothing
    ]);

"""
	$(SIGNATURES)

Line graph. Simple wrapper around `lines`.
"""
function line_plot(xV, yV :: AbstractVector{F};
    fig = blank_plot(), pos = (1,1), kwargs ...) where F

    args = merge(line_defaults(), kwargs);
    ax = make_axis(fig, pos; args...);
    lines!(ax, xV, yV; args...);
    return fig, ax
end

function line_plot(yV :: AbstractVector{F};
    fig = blank_plot(), pos = (1,1), kwargs ...) where F
    return line_plot(1 : length(yV), yV; fig = fig, pos = pos, kwargs...);
end

"""
	$(SIGNATURES)

Multi-line plot from data matrix. Columns are series.

This is really the same as `Makie.series`.

# Arguments

- legPos: `:below` places a legend below the graph (in a new frame).
- labelV: labels for legend (optional).

Fix: Legend colors are wrong +++

# Example

```julia
fig, ax = line_plot(1:4, rand(4,3); labelV = ["Lbl \$j" for j = 1 : 3]);
axislegend()
```
"""
function line_plot(xV, yM :: AbstractMatrix{F};
    fig = blank_plot(), pos = (1,1),
    legPos = :none,
    kwargs...) where F

    args = merge(line_defaults(), kwargs);
    ax = make_axis(fig, pos; args...);
    nr, nc = size(yM);
    @assert nr == length(xV);
    for j = 1 : nc
        add_line!(ax, xV, yM[:, j];
            label = get_idx(args[:labelV], j),
            color = get_colors(j, nc),  # (fill(j, nr), nc),
            # color = fill(j, nr), colorrange = (1, nc),
            args...);
    end
    if !isnothing(args[:labelV])  &&  (legPos == :below)
        legPos = (pos[1] + 1, pos[2]);
        fig[legPos...] = Legend(fig, ax;
            orientation = :horizontal,
            tellwidth = false, tellheight = true);
    end
    return fig, ax
end

"""
	$(SIGNATURES)

Add line to a plot.
"""
function add_line!(ax :: Axis, x, y; kwargs...)
    lines!(ax, x, y; kwargs...);
end

add_line!(p :: Makie.FigureAxisPlot, x, y; kwargs...) =
    add_line!(p.axis, x, y; kwargs...);

"""
	$(SIGNATURES)

Add error band.

How to get these colored is not clear. Adding
`color = (get_colors(j, n), 0.2)`
achieves transparency, but no color (all grey).
"""
add_error_band!(ax :: Axis, x, y, errorV; kwargs...) =
    band!(ax, x, y .- errorV, y .+ errorV; kwargs...)

add_error_band!(ax :: Axis, x, y, errorV :: Nothing; kwargs...) = nothing;

# -----------


text:
c1,c2,c3,c4=c
A11,A12,A13 = A[1,:]
A22,A23= A[2,2:3]
A33=A[3,3]

#A11
C1111=A11*c1*c1
C1112=A11*c1*c2
C1113=A11*c1*c3
C1114=A11*c1*c4
C1122=A11*c2*c2
C1123=A11*c2*c3
C1124=A11*c2*c4
C1133=A11*c3*c3
C1134=A11*c3*c4
C1144=A11*c4*c4

#A12
C1211=A12*c1*c1
C1212=A12*c1*c2
C1213=A12*c1*c3
C1214=A12*c1*c4
C1222=A12*c2*c2
C1223=A12*c2*c3
C1224=A12*c2*c4
C1233=A12*c3*c3
C1234=A12*c3*c4
C1244=A12*c4*c4

#A13
C1311=A13*c1*c1
C1312=A13*c1*c2
C1313=A13*c1*c3
C1314=A13*c1*c4
C1322=A13*c2*c2
C1323=A13*c2*c3
C1324=A13*c2*c4
C1333=A13*c3*c3
C1334=A13*c3*c4
C1344=A13*c4*c4

#A22
C2211=A22*c1*c1
C2212=A22*c1*c2
C2213=A22*c1*c3
C2214=A22*c1*c4
C2222=A22*c2*c2
C2223=A22*c2*c3
C2224=A22*c2*c4
C2233=A22*c3*c3
C2234=A22*c3*c4
C2244=A22*c4*c4

#A23
C2311=A23*c1*c1
C2312=A23*c1*c2
C2313=A23*c1*c3
C2314=A23*c1*c4
C2322=A23*c2*c2
C2323=A23*c2*c3
C2324=A23*c2*c4
C2333=A23*c3*c3
C2334=A23*c3*c4
C2344=A23*c4*c4

#A33
C3311=A33*c1*c1
C3312=A33*c1*c2
C3313=A33*c1*c3
C3314=A33*c1*c4
C3322=A33*c2*c2
C3323=A33*c2*c3
C3324=A33*c2*c4
C3333=A33*c3*c3
C3334=A33*c3*c4
C3344=A33*c4*c4


text:
# mcmc.jl

export mcmc!, adapt!; #, est_imp;

"""
mcmc!(model, n_keep[, monitor, report_filename="out_progress.txt",
thin=1, report_freq=10000])
"""
function mcmc!(model::Model_BNP_WMReg_Joint, n_keep::Int,
    updatevars::Updatevars_BNP_WMReg_Joint=Updatevars_BNP_WMReg_Joint(true, true, false, true, true, true),
    monitor::Monitor_BNP_WMReg_Joint=Monitor_BNP_WMReg_Joint(true, false, false, false),
    report_filename::Union{String, Nothing}="out_progress.txt", thin::Int=1,
    report_freq::Int=100)

    writeout = typeof(report_filename) == String

    ## output files
    if writeout
        report_file = open(report_filename, "a+")
        write(report_file, "Commencing MCMC at $(Dates.now()) from iteration $(model.state.iter) for $(n_keep * thin) iterations.\n\n")
        close(report_file)
    end

    sims, symb_monitor = postSimsInit_BNP_WMReg_Joint(monitor, n_keep, model.state)
    start_accpt = deepcopy(model.state.accpt)
    start_iter = deepcopy(model.state.iter)
    prev_accpt = deepcopy(model.state.accpt) # set every report_freq iterations

    yX = hcat(model.y, model.X) # useful for allocation update

    lfc_γon = 0.0 .* model.state.γ

    ## Initialize lNX and lωNX_vec
    if model.state.iter == 0
        if model.γ_type in (:fixed, :global)
            model.state.lNX, model.state.lωNX_vec = lNXmat_lωNXvec(model, model.state.γ)
        elseif model.γ_type == :local
            model.state.lNX, model.state.lωNX_vec = lNXmat_lωNXvec(model, trues(model.K))
            lNXmat_lωNXvec!(model, model.state.γ) # update each mixture component h in 1:H
        end
    end

    ## sampling
    for i in 1:n_keep
        for j in 1:thin

            if updatevars.γ && (model.γ_type != :fixed)

                if model.γ_type == :local
                    update_ξ!(model)
                    update_π_γ!(model)
                end

                lfc_γon = update_γ!(model)
            end

            if updatevars.η
                Λβ0star_ηy = model.state.Λ0star_ηy * model.state.β0star_ηy
                βΛβ0star_ηy = PDMats.quad(model.state.Λ0star_ηy, model.state.β0star_ηy)

                # update_ηy = (j == thin) # only update η_y on collection scan
                update_ηy = true

                if model.K > 1
                    for h = 1:model.H
                        update_η_h_Met!(model, h, Λβ0star_ηy, βΛβ0star_ηy, update_ηy=update_ηy)
                    end
                else
                    for h = 1:model.H
                        update_η_h_Met_K1!(model, h, Λβ0star_ηy, βΛβ0star_ηy, update_ηy=update_ηy)
                    end
                end
            end

            if updatevars.lω
                update_vlω_mvSlice!(model)
            end

            if updatevars.α
                model.state.α = rand(BayesInference.post_alphaDP(model.H, model.state.lω[model.H],
                    model.prior.α_sh, model.prior.α_rate))
            end

            if updatevars.G0
                update_G0!(model)
            end

            ## do last, followed by llik calculation (no longer necessary because marginalized allocation update does not calculate llik numerator)

            if updatevars.S
                update_alloc!(model)
            end

            if j == thin # llik doesn't get used in any updates, so you can calculate it for inference samples only
                llik_num_mat = llik_numerator(model)
                model.state.llik = llik_BNP_WMReg_Joint(llik_num_mat, model.state.lωNX_vec)
            end

            model.state.iter += 1
            if model.state.iter % report_freq == 0
                if writeout
                    report_file = open(report_filename, "a+")
                    write(report_file, "Iter $(model.state.iter) at $(Dates.now())\n")
                    write(report_file, "Log-likelihood (last inference sample) $(model.state.llik)\n")
                    write(report_file, "Current Metropolis acceptance rates: $(float((model.state.accpt - prev_accpt) / report_freq))\n")
                    write(report_file, "Current allocation: $(counts(model.state.S, 1:model.H))\n")
                    write(report_file, "Current final weight: $(exp(model.state.lω[model.H]))\n")
                    if model.γ_type == :global
                        write(report_file, "Current variable selection: $(1*model.state.γ)\n\n")
                    elseif model.γ_type == :local
                        γ_occupied = counts(model.state.S, 1:model.H)'model.state.γ ./ float(model.n)
                        write(report_file, "Current variable selection of occupied clusters: $(γ_occupied)\n\n")
                    end
                    close(report_file)
                end
                prev_accpt = deepcopy(model.state.accpt)
            end
        end

        sims[i] = BayesInference.deepcopyFields(model.state, symb_monitor)
        sims[i][:Scounts] = counts(model.state.S, 1:model.H)

        if updatevars.γ
            # sims[i][:lwimp] = model.state.lwimp
            if model.γ_type == :global
                sims[i][:lfc_on] = deepcopy(lfc_γon)
            elseif model.γ_type == :local
                sims[i][:γ_occupied] = vec(sims[i][:Scounts]'model.state.γ ./ float(model.n))
            end
        end

    end

    return (sims, float((model.state.accpt - start_accpt)/(model.state.iter - start_iter)) )
end

function adjust_from_accptr(accptr::T, target::T, adjust_bnds::Array{T,1}) where T <: Real
    if accptr < target
        out = (1.0 - adjust_bnds[1]) * accptr / target + adjust_bnds[1]
    else
        out = (adjust_bnds[2] - 1.0) * (accptr - target) / (1.0 - target) + 1.0
    end

    return out
end

function adapt!(model::Model_BNP_WMReg_Joint;
    n_iter_collectSS::Int=1000, n_iter_scale::Int=500,
    adapt_thin::Int=1,
    updatevars::Updatevars_BNP_WMReg_Joint,
    report_filename::String="out_progress.txt",
    maxtries::Int=50,
    accptr_bnds::Vector{T}=[0.23, 0.44], adjust_bnds::Vector{T}=[0.01, 10.0],
    tune_only::Bool=false) where T <: Real

    target = StatsBase.mean(accptr_bnds)
    d = size(model.state.runningsum_ηlδx, 2)
    collect_scale = 2.38^2 / float(d)

    ## initial runs
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 1 of 4 (initial scaling) at $(Dates.now())\n\n")
    close(report_file)

    model.state.adapt = false
    reset_adapt!(model)
    tries = 0

    fails = trues(model.H)

    while any(fails)
        tries += 1
        if tries > maxtries
            println("Exceeded maximum adaptation attempts, Phase 1\n")
            report_file = open(report_filename, "a+")
            write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 1\n\n")
            close(report_file)
            break
        end

        sims, accptr = mcmc!(model, n_iter_scale,
            updatevars,
            Monitor_BNP_WMReg_Joint(false, false, false, false),
            tries == 1 ? report_filename : nothing,
            1, n_iter_scale)

        for h = 1:model.H
            fails[h] = (accptr[h] < accptr_bnds[1])
            if fails[h]
                model.state.cSig_ηlδx[h] *= adjust_from_accptr(accptr[h], target, adjust_bnds)
            end
        end

    end

    ## local scaling
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 2 of 4 (local scaling) at $(Dates.now())\n\n")
    close(report_file)

    model.state.adapt = false
    reset_adapt!(model)

    localtarget = target * 1.0

    for ii = 1:3

        for group in keys(model.indx_ηx)

            ig = model.indx_ηx[group]

            tries = 0
            fails = trues(model.H)
            while any(fails)
                tries += 1
                if tries > maxtries
                    println("Exceeded maximum adaptation attempts, Phase 2\n")
                    report_file = open(report_filename, "a+")
                    write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 2\n\n")
                    close(report_file)
                    break
                end

                sims, accptr = mcmc!(model, n_iter_scale, updatevars,
                    Monitor_BNP_WMReg_Joint(false, false, false, false),
                    nothing,
                    1, n_iter_scale)

                for h = 1:model.H
                    too_low = accptr[h] < (accptr_bnds[1] * 0.5)
                    too_high = accptr[h] > (accptr_bnds[2])

                    if too_low || too_high

                        fails[h] = true

                        tmp = Matrix(model.state.cSig_ηlδx[h])
                        σ = sqrt.(LinearAlgebra.diag(tmp))
                        ρ = StatsBase.cov2cor(tmp, σ)

                        σ[ig] *= adjust_from_accptr(accptr[h], localtarget, adjust_bnds)
                        tmp = StatsBase.cor2cov(ρ, σ)
                        tmp += Diagonal(fill(0.1*minimum(σ), size(tmp,1)))
                        model.state.cSig_ηlδx[h] = PDMat_adj(tmp)

                    else
                        fails[h] = false
                    end

                end

            end

        end
    end

    if !tune_only
        ## cΣ collection
        report_file = open(report_filename, "a+")
        write(report_file, "\n\nBeginning Adaptation Phase 3 of 4 (covariance collection) at $(Dates.now())\n\n")
        close(report_file)

        reset_adapt!(model)
        model.state.adapt = true
        model.state.adapt_thin = deepcopy(adapt_thin)

        sims, accptr = mcmc!(model, n_iter_collectSS, updatevars,
            Monitor_BNP_WMReg_Joint(false, false, false, false),
            report_filename, 1, n_iter_collectSS)

        for h = 1:model.H
            Sighat = model.state.runningSS_ηlδx[h,:,:] / float(model.state.adapt_iter)
            minSighat = minimum(abs.(Sighat))
            SighatPD = Sighat + Matrix(Diagonal(fill(0.1*minSighat, d)))
            model.state.cSig_ηlδx[h] = PDMat_adj(collect_scale * SighatPD)
        end

        ## final scaling
        report_file = open(report_filename, "a+")
        write(report_file, "\n\nBeginning Adaptation Phase 4 of 4 (final scaling) at $(Dates.now())\n\n")
        close(report_file)

        model.state.adapt = false
        reset_adapt!(model)
        tries = 0

        fails = trues(model.H)

        while any(fails)
            tries += 1
            if tries > maxtries
                println("Exceeded maximum adaptation attempts, Phase 4\n")
                report_file = open(report_filename, "a+")
                write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 4\n\n")
                close(report_file)
                break
            end

            sims, accptr = mcmc!(model, n_iter_scale, updatevars,
                Monitor_BNP_WMReg_Joint(false, false, false, false),
                tries == 1 ? report_filename : nothing,
                1, n_iter_scale)

            for h = 1:model.H
                too_low = accptr[h] < accptr_bnds[1]
                too_high = accptr[h] > accptr_bnds[2]

                if too_low || too_high
                    fails[h] = true
                    model.state.cSig_ηlδx[h] *= adjust_from_accptr(accptr[h], target, adjust_bnds)
                else
                    fails[h] = false
                end

            end

        end
    end

    report_file = open(report_filename, "a+")
    write(report_file, "\n\nAdaptation concluded at $(Dates.now())\n\n")
    close(report_file)

    ## note that mcmc! also closes the report file
    reset_adapt!(model)
    model.state.adapt = false

    return nothing
end

## estimate quantities with importance weights
# function est_imp(x::Vector{Float64}, weights::Vector{Float64})
#     sumw = sum(weights)
#     return sum(x .* weights) / sumw
# end


text:
using Pkg, Documenter, Literate, SmoothedParticles

#
# Replace SOURCE_URL marker with github url of source
#
function replace_source_url(input,source_url)
    lines_in = collect(eachline(IOBuffer(input)))
    lines_out=IOBuffer()
    for line in lines_in
        println(lines_out,replace(line,"SOURCE_URL" => source_url))
    end
    return String(take!(lines_out))
end

function make_all()
    #
    # Generate Markdown pages from examples
    #
    examples = [
        "static_container.jl",
        "collapse_dry.jl",
        "collapse_dry_implicit.jl",
        "cavity_flow.jl",
        "cylinder.jl",
        "drop.jl",
        "collapse3d.jl",
        "collapse_symplectic.jl",
        "Kepler_vortex.jl"
    ]
    example_md_dir  = joinpath(@__DIR__,"src","examples")
    generated_examples = []

    for example in examples
        base,ext=splitext(example)
        if ext==".jl"
            Literate.markdown(joinpath(@__DIR__,"..","examples",example),
                              example_md_dir,
                              documenter=false,
                              info=false
                             )
            push!(generated_examples, joinpath("examples", base*".md"))
        end
    end

    #generated_examples=joinpath.("examples",readdir(example_md_dir, sort = false))
    println(generated_examples)

    makedocs(
        sitename="SmoothedParticles.jl",
        clean = true,
        doctest = true,
        authors = "O. Kincl",
        modules = [SmoothedParticles],
        format = Documenter.HTML(prettyurls = false),
        pages=[
            "Home"=>"index.md",
            "Examples" => generated_examples,
            "API Documentation" => [
                "core.md",
                "kernels.md",
                "structs.md",
                "geometry.md",
                "IO.md",
                "grids.md"
            ],
        ]
    )
end

make_all()

deploydocs(
    repo = "github.com/OndrejKincl/SmoothedParticles.jl.git",
)


text:
using HackerNewsAnalysis
using Test

@testset "HackerNewsAnalysis.jl" begin
end


text:
using FactCheck

include("models.jl")

exitstatus()


text:
# The serialization header can change between Julia versions. When the header does change
# then we can expect the "legacy deserialization" testset to fail and we'll need to
# conditionally set the value of the header here.
const SERIALIZED_HEADER = "7JL\n\x04\0\0"  # Header for Julia 1.5+

# Declare a new `isnan` and `isinf` functions to avoid type piracy. These new function work
# with `Char` and `Period` types as well as other types.
# Note: A generic `isfinite` is declared in Intervals.
isnan(x) = (x != x)::Bool
isnan(x::Real) = Base.isnan(x)

isinf(x) = !isnan(x) && !isfinite(x)
isinf(x::Real) = Base.isinf(x)


text:
module selection
importall common, species
using core, fitness, individual
export Selection, select, select_ids, SelectionDefinition

"""
The base type used by all selection operator definitions.
"""
abstract SelectionDefinition <: OperatorDefinition

"""
The base type used by all selection operators.
"""
abstract Selection <: Operator

"""
Selects a given number of individuals from a collection of candidate individuals
according to a given selection method, returning the selected individuals as a
new individual collection.
"""
select(s::Selection, buffer::IndividualCollection, sp::Species, ic::IndividualCollection, n::Integer) =
  write_to_individual_collection(buffer, ic, select_ids(s, sp.fitness, indexed_fitnesses(ic), n))

"""
Selects a specified number of individuals from a list of candidates, provided
as a vector of tuples, specifying the id and fitness of each candidate.

Returns the internal IDs of the selected individuals.
"""
select_ids{F}(s::Selection, ::Species, ::Vector{Tuple{Int, F}}, ::Integer) =
  error("Unimplemented `select_ids' method for this selection operator: $(typeof(s)).")

"""
Prepares a given selection operator for the breeding process.
"""
prepare(s::Selection, ic::IndividualCollection) = ic

# Load each of the selection operators.
include("selection/random.jl")
include("selection/tournament.jl")
end


text:
# Autogenerated wrapper script for primecount_jll for i686-w64-mingw32-cxx03
export libprimecount, libprimesieve, primecount

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("primecount")
JLLWrappers.@declare_library_product(libprimecount, "libprimecount.dll")
JLLWrappers.@declare_library_product(libprimesieve, "libprimesieve.dll")
JLLWrappers.@declare_executable_product(primecount)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libprimecount,
        "bin\\libprimecount.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libprimesieve,
        "bin\\libprimesieve.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        primecount,
        "bin\\primecount.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestDecisionTree
using Test
using StrayCopses
using StrayCopses: feature_best_split, best_split, Node, DecisionTreeContainer, process_node, create_containers
using StableRNGs

@testset "Basic split" begin
    X = reshape([1., 2., 3., 4.], :, 1)
    y = [1, 1, 2, 2]

    containers = create_containers(2, y)
    @test feature_best_split(containers, X, y, 2, 1) == (val = 3.000, impurity = 0.5)

    rng = StableRNG(2020)
    N = 5
    n_classes = 3
    X = reshape(rand(rng, 2*N), :, 2)
    y = rand(rng, 1:n_classes, N)
    containers = StrayCopses.create_containers(n_classes, y)

    res = feature_best_split(containers, X, y, n_classes, 1)
    @test res.val ≈ 0.257183060127947
    @test res.impurity ≈ 0.2599999999999999

    res = feature_best_split(containers, X, y, n_classes, 2)
    @test res.val ≈ 0.9442017794271014
    @test res.impurity ≈ 0.36
end

@testset "Best split" begin
    X = collect(reshape([1., 2., 3., 4., 1., 2., 2., 4.], :, 2))
    y = [1, 1, 2, 2]

    containers = create_containers(2, y)
    @test best_split(X, y, 2, 1:2) == (feature = 1, val = 3.0000)
    @test best_split(X, y, 2, 1:1) == (feature = 1, val = 3.0000)
    @test best_split(X, y, 2, [2]) == (feature = 2, val = 2.0000)

    X = [6.5 2.8 4.6 1.5; 6.5 3.0 5.8 2.2]
    y = [2, 3]
    @test best_split(X, y, 3, 1:4) == (feature = 2, val = 3.0)
end

@testset "random dataset" begin
    rng = StableRNG(2020)
    X = rand(rng, 10, 10)
    X = floor.(X .* 100)
    y = rand(rng, 1:3, 10)
    T = eltype(X)

    root = Node{T}()
    dtc = DecisionTreeContainer{T}(root, 10, 3, 6, 1)
    process_node(dtc, root, X, y, rng, 1:10)
    io = IOBuffer()
    print(io, root)
    tree = String(take!(io))
    @test tree == "[X4 < 66.0]\n↦[X10 < 17.0]\n⎸↦[X3 < 78.0]\n⎸⎸↳[1]\n⎸⎸↳[2]\n⎸↳[1]\n↦[X4 < 83.0]\n ↳[2]\n ↳[3]\n"
end

end # module


text:
function printd(msg::String, level::PrintL, debug_level::PrintL)
    if level <= debug_level && debug_level > mute
        println(msg)
    end
end

dxy(x,y) = abs(x-y)

"""
image_fom (imge figure of merit) computes the sum of the
differences between the pixels of two images, defined as 2-d Float64 arrays
providing a (naive) metric of the difference between both.
"""
function image_fom(img1::Array{Float64,2}, img2::Array{Float64,2})
    dff = sum(broadcast(dxy, img1, img2))
    return dff /sum(img1)
end

"""
Find the (i,j) indexes defining the center of an array
"""
function array_center_indexes(a::AbstractArray)
    is = size(a)./2     # index center
    ic = map(ceil, is)  # take the ceiling
    return map(x -> convert(Int, x), ic) # and cast to Int
end

"""
Return a view of img. The view is taken from the center
with a width hw
"""
function center_view(img::Array{Float64,2}, hw::Int)
    imc = array_center_indexes(img)
    il  = imc[1] - hw
    ir  = imc[2] + hw
    return view(img, il:ir, il:ir)
end

"""
Flips an array: Equivalent to numpy.flip
"""
function flip2d(x::AbstractArray)
    y = copy(x)
    for d in 1:2
        y = reverse(y, dims=d)
    end
    return y
end

"""
In order to avoid zero divsion, replace any element smaller than eps()
---a function that gives the small number accesible for the type with
the value of eps(). Also, one may want to clip the image (between 0 and 1)
Both tasks can be performed with the function below
"""
function limit_boundaries(a::AbstractArray, xmin::Float64 = eps(), xmax::Float64 = 1e+9)
    map!(x -> x>xmin ? x : xmin, a, a)
    map!(x -> x < xmax ? x : xmax, a, a)
end

"""
A simple 1d gaussian kernel
(yields the same result than the image package kernel)
"""
gxy((x,y)) = exp(-(x^2/2 + y^2/2))
function gaussianxy(x)
    g = gxy.(x)
    s = sum(g)
    g/s
end

max_stars_repo_name: jjgomezcadenas/JImages

max_issues_repo_name: jjgomezcadenas/JImages

max_forks_repo_name: jjgomezcadenas/JImages


text:
# module for simulate and draw graph based slam

module AnimeGraphBasedSlam
  # packages
  using Combinatorics
  using LinearAlgebra
  using Statistics
  using Plots
  pyplot()

  # external modules
  include(joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/obsrv_edge.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/motion_edge.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/map_edge.jl"))

  # read logged trajectory and observation data
  function read_data()
    pose_list = Dict() # trajectory
    obsrv_list = Dict() # observation
    delta = 0.0
    input_list = Dict()

    log_path = "src/slam/graph_based_slam/traj_obsrv_edge_input_log.txt"
    open(joinpath(split(@__FILE__, "src")[1], log_path), "r") do fp
      for line in eachline(fp)
        tmp = split(line)

        type = tmp[1]
        step = tmp[2]
        if type == "x" # [x, y, theta]
          pose_list[step] = [parse(Float64, tmp[3]), parse(Float64, tmp[4]), parse(Float64, tmp[5])]
        elseif type == "z" # (id, [distance, direction])
          if !haskey(obsrv_list, step)
            obsrv_list[step] = []
          end
          push!(obsrv_list[step], (parse(Int64, tmp[3]), [parse(Float64, tmp[4]), parse(Float64, tmp[5])]))
        elseif type == "delta" # time interval
          delta = parse(Float64, tmp[2])
        elseif type == "u" # input
          input_list[step] = [parse(Float64, tmp[3]), parse(Float64, tmp[4])]
        end
      end
    end

    return pose_list, obsrv_list, input_list, delta
  end

  function make_axis()
    plot([], [], aspect_ratio=true, xlabel="X", ylabel="Y",
         xlims=(-5.0, 5.0), ylims=(-5.0, 5.0), legend=false)
  end

  function draw_observation(pose_list, obsrv_list)
    for step in 0:length(pose_list)-1
      step_str = string(step)
      if haskey(obsrv_list, step_str)
        for obsrv in obsrv_list[step_str]
          x, y, theta = pose_list[step_str][1], pose_list[step_str][2], pose_list[step_str][3]
          dist, dir = obsrv[2][1], obsrv[2][2]
          mx = x + dist * cos(theta + dir)
          my = y + dist * sin(theta + dir)
          plot!([x, mx], [y, my], color="pink", legend=false,
                aspect_ratio=true, alpha=0.5)
        end
      end
    end
  end

  function draw_trajectory(pose_list)
    x_list = [pose_list[string(step)][1] for step in 0:length(pose_list)-1]
    y_list = [pose_list[string(step)][2] for step in 0:length(pose_list)-1]
    scatter!(x_list, y_list, markersize=2, markercolor="black")
    plot!(x_list, y_list, linewidth=0.5, color="black")
  end

  function make_edges(pose_list, obsrv_list)
    # record time and observation
    landmark_keys_list = Dict()
    for step in 0:length(pose_list)-1
      step_str = string(step)
      if haskey(obsrv_list, step_str)
        for obsrv in obsrv_list[step_str]
          id_str = string(obsrv[1])
          if !haskey(landmark_keys_list, id_str)
            landmark_keys_list[id_str] = []
          end
          push!(landmark_keys_list[id_str], (step, obsrv))
        end
      end
    end

    edges = []
    for key_id in keys(landmark_keys_list)
      step_pairs =[pair for pair in combinations(landmark_keys_list[key_id], 2)]
      for pair in step_pairs
        push!(edges, ObsrvEdge(pair[1][1], pair[2][1], pair[1][2], pair[2][2], pose_list))
      end
    end

    return edges, landmark_keys_list
  end

  function add_edge(edge, omega, xi)
    # index of precision matrix, xi vector
    f1, f2 = edge.t1*3+1, edge.t2*3+1
    t1, t2 = f1+2, f2+2
    # set matrix
    omega[f1:t1, f1:t1] += edge.omega_upper_left
    omega[f1:t1, f2:t2] += edge.omega_upper_right
    omega[f2:t2, f1:t1] += edge.omega_lower_left
    omega[f2:t2, f2:t2] += edge.omega_lower_right
    # set xi
    xi[f1:t1] += edge.xi_upper
    xi[f2:t2] += edge.xi_lower
  end

  function draw_edges(edges)
    for e in edges
      plot!([e.x1[1], e.x2[1]], [e.x1[2], e.x2[2]], color="red", alpha=0.5)
    end
  end

  function draw_landmarks(landmark_list)
    scatter!([m[1] for m in landmark_list],
             [m[2] for m in landmark_list],
             markershape=:star,
             markercolor=:blue,
             markersize=15)
  end

  function draw(pose_list, obsrv_list, edges, landmark_list, n)
    make_axis()
    draw_observation(pose_list, obsrv_list)
    draw_trajectory(pose_list)
    draw_landmarks(landmark_list)
    annotate!(-3.5, 4.5, "iterate: $(n)")
  end

  function main(is_test=false)
    pose_list, obsrv_list, input_list, delta = read_data()

    # dimension of trajectory vector
    dim = length(pose_list) * 3

    # iterate for optimization
    anime = @animate for n in 1:10000
      edges, landmark_keys_list = make_edges(pose_list, obsrv_list)

      # add motion edges
      for i in 0:length(pose_list)-2
        for j in 1:100
          push!(edges, MotionEdge(i, i+1, pose_list, input_list, delta))
        end
      end

      # add landmark objects
      landmark_list = []
      for lm_pair in landmark_keys_list
        map_edges = []
        head_z = landmark_keys_list[lm_pair[1]][1]
        for z in landmark_keys_list[lm_pair[1]]
          push!(map_edges, MapEdge(z[1], z[2][2], head_z[1], head_z[2][2], pose_list))
        end

        omega_map = zeros(2, 2)
        xi_map = zeros(2)
        for me in map_edges
          omega_map += me.omega
          xi_map += me.xi
        end
        push!(landmark_list, inv(omega_map)*xi_map)
      end

      # initialize precision matrix for graph
      omega_graph = zeros(dim, dim)
      omega_graph[1:3, 1:3] += Matrix{Float64}(I,3,3).*1000000

      # initialize coefficient vector for graph
      xi_graph = zeros(dim)

      # calculate difference of trajectory
      for e in edges
        add_edge(e, omega_graph, xi_graph)
      end
      delta_traj = inv(omega_graph) * xi_graph

      # update trajectory
      for i in 0:length(pose_list)-1
        pose_list[string(i)] += delta_traj[(i*3+1):((i+1)*3)]
      end

      # convergence
      diff = norm(delta_traj)
      draw(pose_list, obsrv_list, edges, landmark_list, n)
      if diff < 0.01
        break
      end
    end

    if is_test == false
      save_path = joinpath(split(@__FILE__, "src")[1], "src/slam/graph_based_slam/est_poses_map_2vars.gif")
      gif(anime, fps=1, save_path)
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy