text:
#== lexer.py

from .utils import Str, classify  #Str not needed

from .exceptions import UnexpectedCharacters, LexError
==#
# Patterns.  We assume all have a value and flags field

abstract type Pattern end

get_flags(p1::Pattern) = p1.flags
get_value(p1::Pattern) = p1.value

Base.:(==)(p1::Pattern, p2::Pattern) = begin
    if typeof(p1) == typeof(p2) && p1.value == p2.value && p1.flags == p2.flags
        return true
    end
    return false
end

Base.hash(p1::Pattern) = begin
    return hash((typeof(p1),p1.value,p1.flags))
end

_get_flags(p1::Pattern,value) = begin
    add_flags(get_flags(p1),value)
end

"""
Add the characters in `flags` as regexp flags for
regular expression `value`
"""
add_flags(flags,value) = begin
    for f in flags
        value = ("(?$f)"*value)
    end
    return value
end

struct PatternStr <: Pattern
    value::String
    flags::Set{Char}
end

struct PatternRE <: Pattern
    value::String
    flags::Set{Char}
    _width
end

PatternRE(value;flags=Set()) = PatternRE(value,flags)
PatternRE(value,flags::String) = PatternRE(value,Set(flags))
PatternStr(value,flags::String) = PatternStr(value,Set(flags))
PatternStr(value,flags::Char) = PatternStr(value,[flags])
PatternRE(value,flags::Char) = PatternRE(value,[flags])
PatternRE(value,flags::Set) = PatternRE(value,flags,get_regexp_width(add_flags(flags,value)))

to_regexp(pre::PatternRE) = add_flags(pre.flags,pre.value)
to_regexp(pstr::PatternStr) = _get_flags(pstr,escape_re_string(pstr.value))
min_width(pre::PatternRE) = pre._width[1]
max_width(pre::PatternRE) = pre._width[2]

min_width(pstr::PatternStr) = length(pstr.value)
max_width(pstr::PatternStr) = min_width(pstr)

struct TerminalDef
    name::String
    pattern::Pattern
    priority::Int64
end

TerminalDef(name,pattern; priority=1) = TerminalDef(name,pattern,priority)

"""
    A string with meta-information, that is produced by the lexer.

    When parsing text, the resulting chunks of the input that haven't been discarded,
    will end up in the tree as Token instances. Token is a subclass of AbstractString
    so normal string comparisons and operations will work as expected.

    Attributes:
        type: Name of the token (as specified in grammar)
        value: Value of the token (redundant, as ``token.value == token`` will always be true)
        pos_in_stream: The index of the token in the text
        line: The line of the token in the text (starting with 1)
        column: The column of the token in the text (starting with 1)
        end_line: The line where the token ends
        end_column: The next column after the end of the token. For example,
            if the token is a single character with a column value of 4,
            end_column will be 5.
        end_pos: the index where the token ends (basically ``pos_in_stream + length(token)``)
"""
mutable struct Token <: AbstractString
    type_::String  #type is a Julia keyword
    pos_in_stream::Union{Int,Nothing}
    value::String
    line::Union{Int,Nothing}
    column::Union{Int,Nothing}
    end_line::Union{Int,Nothing}
    end_column::Union{Int,Nothing}
    end_pos::Union{Int,Nothing}
end

Token(type_,value;pos_in_stream=nothing,line=nothing,column=nothing) =
    Token(type_,pos_in_stream,String(value),line,column,nothing,nothing,nothing)

# Reimplement some string methods for later use...
Base.ncodeunits(t::Token) = ncodeunits(t.value)
Base.codeunit(t::Token) = codeunit(t.value)
Base.codeunit(t::Token,i::Integer) = codeunit(t.value,i)
Base.isvalid(t::Token,i::Int64) = isvalid(t.value,i)
Base.iterate(t::Token) = Base.iterate(t.value)
Base.iterate(t::Token,i::Int64) = Base.iterate(t.value,i)

# Lexer.py implements a generic equality, where any comparison other
# than with another token defaults to string comparison
Base.:(==)(t1::Token,t2::Token) = begin
    if t1.type_ != t2.type_
        return false
    end
    return t1.value == t2.value
end

Base.:(==)(t1::Token,t2::String) = begin
    return isequal(t1.value,t2)
end

Base.:(==)(t1::String,t2::Token) = begin
    return isequal(t1,t2.value)
end

Base.hash(t1::Token) = hash(t1.value)

Base.show(io::IO,t::Token) = print(io,"Token($(t.type_), $(t.value))")

update(t::Token;type_=nothing,value=nothing) = begin
    return new_borrow_pos(type_!=nothing ? type_ : t.type_,
                          value!=nothing ? value : t.value,
                          t)
end

new_borrow_pos(type_,value, borrow_t) = Token(type_,borrow_t.pos_in_stream,value,
                                              borrow_t.line, borrow_t.column,borrow_t.end_line,
                                              borrow_t.end_column,borrow_t.end_pos)

## __reduce__ not reproduced; need to check if we need it.

mutable struct LineCounter
    newline_char::Char
    char_pos::Int #1 is first pos unlike Python
    line::Int
    column::Int
    line_start_pos::Int
end

LineCounter() = LineCounter('\n',1,1,1,1)

LineCounter(newline_char) = LineCounter(newline_char,1,1,1,1)
"""
feed

Consume a token and calculate the new line & column.
As an optional optimization, set test_newline=false if token doesn't contain a newline.
"""
feed!(lc::LineCounter,token;test_newline=true) = begin
    #println("Handling $token")
    if test_newline
        newlines = count(x -> x == lc.newline_char, token)
        if newlines > 0
            lc.line = lc.line + newlines
            lc.line_start_pos = lc.char_pos + first(findlast("$(lc.newline_char)",token))
            #println("Line $(lc.line) starts at $(lc.line_start_pos)")
        end
    end
    lc.char_pos += ncodeunits(token)  #Unicode
    lc.column = lc.char_pos - lc.line_start_pos + 1
    #println("At end column is $(lc.column) for char pos $(lc.char_pos) token length $(length(token))")
end

#==
mutable struct _Lex
    lexer
    state
end

_Lex(lexer;state=nothing) = _Lex(lexer,state)
==#

abstract type Lexer end

make_lexer_state(text) = begin
    line_ctr = LineCounter('\n')
    return LexerState(text,line_ctr)
end

# UnlessCallback appears to be a closure?
#struct UnlessCallback
#    mres
#end

# This and `_create_unless` check for strings that will also match
# REs. If a string pattern is matched by a re pattern, and the
# string has a greater priority, then even if the re pattern
# matches, the actual match will be the string pattern.
#
unless_callback(mres) = function (t)
    # Note that mres is a single tuple as we are not splitting
    # REs into 100-pattern chunks unlike Python
    mre, type_from_index = mres
    #println("Unless...$mre for $(t.value)")
    m = Base.match(mre,t.value,1)
    if m != nothing
        t.type_ = type_from_index[lastmatch(m)]
    end
    return t
end

_create_unless(terminals, g_regex_flags) = begin
    tokens_by_type = classify(terminals,key = t->typeof(t.pattern))
    embedded_strs = Set()
    callback = Dict()
    for retok in get(tokens_by_type,PatternRE,[])
        unless = []
        for strtok in get(tokens_by_type,PatternStr,[])
            if strtok.priority > retok.priority
                continue
            end
            s = strtok.pattern.value
            m = Base.match(Regex(to_regexp(retok.pattern),g_regex_flags),s)
            if m!=nothing && m.match == s
                push!(unless,strtok)
                if get_flags(strtok.pattern) <= get_flags(retok.pattern)
                    push!(embedded_strs,strtok)
                end
            end
        end
        if !isempty(unless)
            callback[retok.name] = unless_callback(build_mres(unless, g_regex_flags, match_whole = true))
        end
    end
    terminals = [t for t in terminals if !(t in embedded_strs)]
    return terminals, callback
end

# Match whole means that the regex
# must match the whole expression e.g. when checking tokens for matches with
# other things. "\A" forces the match to match at the beginning.
build_mres(terminals,g_regex_flags;match_whole=false) = begin
    postfix = ")"
    prefix = raw"\A("
    if match_whole
        postfix = "\$"
        prefix = "^"
    end
    prep_string = join(["(?P<$(t.name)>$(prefix*to_regexp(t.pattern)*postfix))" for t in terminals],"|")
    #println("String for regex: $prep_string")
    mres = Regex(prep_string,g_regex_flags)
    names_by_idx = Base.PCRE.capture_names(mres.regex)
    #println("All regexes now $mres")
    return mres,names_by_idx
end

_regexp_has_newline(r) = begin
    return occursin("\n", r) || occursin("\\n", r) || occursin("[^", r) ||
         occursin("\\s",r) || (occursin("(?s", r) && occursin('.', r))
end

struct TraditionalLexer <: Lexer
    newline_types::Array{String,1}
    ignore_types::Set{String}
    callback::Dict{String,Function}
    terminals::Array{TerminalDef}
    mres::Tuple{Regex,Dict{Int,String}}
end

TraditionalLexer(conf) = begin
    terminals = collect(conf.tokens)
    @assert all(x-> x isa TerminalDef,terminals)
    if !conf.skip_validation
        for t in terminals
            try
                Regex(to_regexp(t.pattern),conf.g_regex_flags)
            catch
                throw(LexError("Cannot compile token $(t.name): $(t.pattern)"))
            end
            if min_width(t.pattern) == 0
                throw(LexError("Lexer does not allow zero-width terminals. ($(t.name): $(t.pattern))"))
            end
        end
        @assert Set(conf.ignore) <= Set([t.name for t in terminals])
    end
    newline_types = [t.name for t in terminals if _regexp_has_newline(to_regexp(t.pattern))]
    ignore_types = Set(conf.ignore)
    sort!(terminals, by= x -> (-x.priority, -max_width(x.pattern),-length(x.pattern.value),x.name))
    # _build called instead in Lark 0.11.1
    temp_terminals, callback = _create_unless(terminals,conf.g_regex_flags)
    for (type_,f) in conf.callbacks
        @assert !(type_ in keys(callback)) "TODO: add call chain" #Lark has CallChain here
        callback[type_] = f
    end
    TraditionalLexer(newline_types,ignore_types,callback,terminals,build_mres(temp_terminals,conf.g_regex_flags))
end

# Note that an empty string will match an empty regex, instead of returning nothing
match(tl::TraditionalLexer,text,pos) = begin
    if tl.mres[1] == r"" return nothing end
    m = Base.match(tl.mres[1],SubString(text,pos))
    if !(m===nothing)
        match_num = findfirst(v -> v!=nothing,m.captures)
        return m.match, tl.mres[2][match_num]
    end
end

next_token(tl::TraditionalLexer,lex_state) = begin
    line_ctr = lex_state.line_ctr
    #println("Pos: $(line_ctr.char_pos) out of $(lex_state.text_length), remaining text\n"*lex_state.text[line_ctr.char_pos:end]*"EOF")
    while line_ctr.char_pos <= lex_state.text_length
        res = match(tl,lex_state.text,line_ctr.char_pos)
        if res === nothing
            allowed = setdiff(Set(values(tl.mres[2])), tl.ignore_types)
            if length(allowed) == 0
                allowed = Set(["<END-OF-FILE>"])
            end
            throw(UnexpectedCharacters(lex_state.text,line_ctr.char_pos,line_ctr.line,line_ctr.column,allowed=allowed,token_history=lex_state.last_token !== nothing ? [lex_state.last_token] : nothing))
        end
        value,type_ = res
        if !(type_ in tl.ignore_types)
            t = Token(type_,value,pos_in_stream=line_ctr.char_pos,line=line_ctr.line,column=line_ctr.column)
            feed!(line_ctr,value, test_newline=type_ in tl.newline_types)
            t.end_line = line_ctr.line
            t.end_column = line_ctr.column
            t.end_pos = line_ctr.char_pos
            if type_ in keys(tl.callback)
                t = tl.callback[type_](t)
                if !(t isa Token) throw(error("Callbacks must return a token (returned $t)")) end
            end
            lex_state.last_token = t
            return t
        else
            if type_ in keys(tl.callback)
                t2 = Token(type_, value, pos_in_stream=line_ctr.char_pos, line=line_ctr.line, column=line_ctr.column)
                tl.callback[type_](t2)
            end
            feed!(line_ctr,value,test_newline=type_ in tl.newline_types)
        end
    end
    return nothing
end

next_token(tl::TraditionalLexer,ls,dummy) = next_token(tl,ls)

include("common.jl")   #definitions for LexerConf and ParserConf

mutable struct LexerState
    text::String
    line_ctr::LineCounter
    last_token::Union{Token,Nothing}
    text_length::Int64   #optimisation for speed
end

LexerState(text,line_ctr;last_token=nothing) = begin
    LexerState(text,line_ctr,last_token,ncodeunits(text))
end

struct ContextualLexer <: Lexer
    lexers::Dict{Any,TraditionalLexer}
    root_lexer::TraditionalLexer
end

ContextualLexer(conf::LexerConf,states;always_accept=()) = begin
    terminals = collect(conf.tokens)
    tokens_by_name = Dict()
    for t in terminals
        @assert !(t.name in keys(tokens_by_name))
        tokens_by_name[t.name] = t
    end

    trad_conf = reconfigure(conf,terminals)
    lexer_by_tokens = Dict()
    lexers = Dict()
    lexer = missing  #define in outer scope
    for (state,accepts) in states
        #println("State $state accepts $accepts")
        key = Set(accepts)
        try
            lexer = lexer_by_tokens[key]
        catch KeyError
            accepts = union( Set(accepts), conf.ignore, always_accept)
            state_tokens = [tokens_by_name[n] for n in accepts if (n != nothing && n in keys(tokens_by_name))]
            lexer_conf = reconfigure(trad_conf,state_tokens)
            lexer = TraditionalLexer(lexer_conf)
            lexer_by_tokens[key] = lexer
            #println("Added lexer that accepts $key")
        end
        lexers[state] = lexer
    end
    @assert trad_conf.tokens == terminals
    root_lexer = TraditionalLexer(trad_conf)
    ContextualLexer(lexers,root_lexer)
end

# Was lex
next_token(cl::ContextualLexer,lexer_state,parser_state) = begin
    try
        lexer = cl.lexers[position(parser_state)]
        return next_token(lexer,lexer_state)
    catch e
        if e isa EOFError
            return nothing
        elseif e isa UnexpectedCharacters
            token = next_token(cl.root_lexer,lexer_state)
            throw(UnexpectedToken(token,e.allowed,state=position(parser_state)))
        else
            rethrow(e)
        end
    end
end

struct LexerThread
    lexer::Lexer
    state::LexerState
end

LexerThread(lexer,text::String) = LexerThread(lexer,make_lexer_state(text))

lex(l::LexerThread,parser_state) = l.lexer.lex(l.state,parser_state)

### Not in Lark: an iterator-based approach
#==
The parser state is passed in to the constructor. It would have been
created by the parser, and the current state is altered by the
parser externally.
==#
struct LexerParser
    lexer::Lexer
    lexer_state::LexerState
    parser_state
end

LexerParser(lexer::Lexer,text::AbstractString,ps) = begin
    LexerParser(lexer,make_lexer_state(text),ps)
end

LexerParser(lt::LexerThread,ps) = begin
    LexerParser(lt.lexer,lt.state,ps)
end

Base.iterate(l::LexerParser) = begin
    n = next_token(l.lexer,l.lexer_state,l.parser_state)
    if n === nothing return n end
    return (n,0)
end

Base.iterate(l::LexerParser,dummy) = begin
    n = next_token(l.lexer,l.lexer_state,l.parser_state)
    if n === nothing return n end
    return (n,0)
end


text:
module Sparsity

import MultivariatePolynomials
const MP = MultivariatePolynomials
import MultivariateBases
const MB = MultivariateBases
using SemialgebraicSets

include("ChordalExtensionGraph.jl")
using .ChordalExtensionGraph: ChordalCompletion, ClusterCompletion
export ChordalCompletion, ClusterCompletion

import SumOfSquares

export SignSymmetry
abstract type Pattern end
struct NoPattern <: Pattern end

include("xor_space.jl")
include("sign.jl")
include("variable.jl")
include("monomial.jl")

include("preorder.jl")
include("ideal.jl")

end

ext: jl

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# generate and open an example table for tweaking CSS

using BrowseTables

table = collect(NamedTuple{(:a,:b,:c), Tuple{Union{Int,Missing},Float64,Union{String,Nothing}}},
                iszero(i % 10) ? (a = missing, b = NaN, c = nothing) :
                (a = i, b = Float64(i), c = ('a'+(i%26))^((i%7+1))) for i in 1:100)

open_html_table(table; options = TableOptions(; css_inline = false))


text:
using DataStructures: PriorityQueue, enqueue!, dequeue_pair!
using DataFrames: DataFrame, columnindex
using CategoricalArrays: CategoricalArray

struct CausalSearchResult
    S                               ::Union{Vector{Int64}, Vector{Symbol}}
    confint                         ::Matrix{Float64}
    trace_confint_min               ::Matrix{Float64}
    trace_confint_max               ::Matrix{Float64}
    trace_p_values                  ::Vector{Float64}
    α                               ::Float64
    p                               ::Int64
    variables_considered            ::Union{Vector{Int64}, Vector{Symbol}}
    selection_only                  ::Bool
    model_reject                    ::Bool
    function CausalSearchResult(S, confint, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only=false)
        new(collect(S), confint, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered, selection_only, false)
    end
    function CausalSearchResult(S, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only=true)
        confint = Matrix{Float64}(undef, length(S), 2)
        confint[:, 1] .= -Inf
        confint[:, 2] .= Inf
        new(collect(S), confint, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered, selection_only, false)
    end
    function CausalSearchResult(trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only=true)
        new(Vector{Int64}(), Matrix{Float64}(undef, 0, 2), trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered, selection_only, true)
    end
end

function Base.show(io::IO, result::CausalSearchResult)
    if result.model_reject
        printstyled("\n * The whole model is rejected! i.e., Y | (all variables) is not invariant.\n", color=:light_magenta, bold=true)
    elseif isempty(result.S)
        # no causal variable found
        printstyled("\n * Found no causal variable (empty intersection).\n", color=:light_magenta, bold=true)
    else
        printstyled("\n * Causal variables include: $(result.S)\n", color=:green, bold=true)
        if result.selection_only
            printstyled(" * No confidence intervals produced in selection_only mode\n", color=:light_blue)
        else
            println(@sprintf("\n   %-10s \t %-3s %% \t\t %-3s %%", "variable", result.α * 100, (1 - result.α) * 100))
            for i in result.S
                if isa(i, Int64)
                    j = i
                else
                    j = find(result.variables_considered.==i)[1]
                end
                println(@sprintf("   %-10s \t% -04.4f \t% -04.4f", i, result.confint[j, 1], result.confint[j, 2]))
            end
        end
    end
    if length(result.variables_considered) == result.p
        println("\n ⋅ Variables considered include 1:$(result.p)")
    else
        println("\n ⋅ Variables considered include $(result.variables_considered)")
    end
end

"""
    causalSearch(X, y, env, [ S=1:size(X,2) ];
                 α=0.01, method="chow", p_max=8, screen="auto", verbose=true,
                 selection_only=false, iterate_all=false, n_max_for_exact=5000, max_num_true_causes=Inf)

    causalSearch(df::DataFrame, target::(Int64 or Symbol), env; ...)

Searching over subsets in `X[,S]` for direct causes of `y`

# Arguments
* `X`:                 either an n x p matrix or a `DataFrames.DataFrame`
* `y`:                 vector of n (`X` and `y` can be alternatively specified by `df` and its column `target`)
* `env`:               environment indicators (rows of X): 1, 2, ...
* `S`:                 set of variables (col indices of X) to search, can be a Vector or a Set
* `α`:                 significance level (e.g. 0.01)
* `method`:
     + `"chow"` for Gaussian linear regression, combined two-sample chow test
     + `"logistic-LR"` for logistic regression (`y` consists of 0 and 1), combined likelihood-ratio test
     + `"logistic-SF"`  for logistic regression (`y` consists of 0 and 1), combined Sukhatme-Fisher test
* `p_max`:             maximum number of variables to consider.
                       will method in `screen` to screen out `p_max` number of variables if `p_max < |S|`.
                       (set to `Inf` if want no screening)
* `screen`:
     + `"lasso"`: with lasso (from glmnet) solution path (see `screen_lasso`)
     + `"HOLP"`: "High dimensional ordinary least squares projection" method of Wang & Leng, only when p ≧ n (see `screen_HOLP`)
     + `"auto"`: use `"HOLP"` when p > n, and `"lasso"` otherwise
* `verbose`:           if true, will print each subset tested
* `selection_only`:    if true, will prune supersets of an invariant set;
                       but not able to produce valid confidence intervals
* `iterate_all`:       if true, will iterate over all subsets to ensure validity of confidence intervals
                       (if model is not rejected)
* `n_max_for_exact`:   maximum number of observations of an environment for exact testing;
                       otherwise a subsample of n_max rows will be used
* `max_num_true_causes`: maximum number of true causal variables; if specified to smaller than `|S|`,
                        it will skip testing subsets with bigger size than `max_num_true_causes`.
"""
function causalSearch(X::Union{Matrix{Float64}, DataFrame}, y::Vector{Float64}, env::Vector{Int64}, S=1:size(X,2);
                      α=0.01, method="chow", p_max=8, screen="auto", verbose=true,
                      selection_only=false, iterate_all=false,
                      n_max_for_exact=5000, max_num_true_causes=Inf)
    @assert size(X, 1) == length(y) == length(env)
    if method=="chow"
        model = "linear"
        if isa(X, DataFrame)
            X = Matrix{Float64}(X)  # note: current linear fitting has to work with Matrix{Float64}
        end
    elseif method=="logistic-LR" || method=="logistic-SF"
        model = "logistic"
        # combine into a DataFrame (note: GLM.jl has to work with DataFrame)
        @assert all((y.==1) .| (y.==0))
        df = DataFrame(isa(X, DataFrame) ? hcat(X, y, makeunique=true) : hcat(X, y))
        for _col in propertynames(df)
            if isa(df[!, _col], CategoricalArray)
                @assert length(unique(df[!, _col])) == 2 "categorical variable $_col should be recoded to binary"
            end
        end
        target = propertynames(df)[end]  # target is the last column
    else
        error("method must be one of: `chow`, `logistic-LR`, `logistic-SF`")
    end
    S = collect(S)
    S = unique(S)
    p = size(X, 2)
    if p_max < length(S)
        @assert model=="linear" "screening unsupported for GLM"
        q = length(S)
        if screen == "auto"
            if q <= size(X, 1)
                screen = "lasso"
            else
                screen = "HOLP"
            end
        end
        if screen == "lasso"
            S = S[screen_lasso(X[:, S], y, p_max)]
        elseif screen == "HOLP"
            S = S[screen_HOLP(X[:, S], y, p_max)]
        else
            error("screen must be one of: `auto`, `lasso`, `HOLP`")
        end
        printstyled("$(length(S)) variables are screened out from $q variables with $screen: $S\n", color=:blue)
    end
    variables_considered = S[:]
    if max_num_true_causes < length(S)
        printstyled("the size of |S| is restricted to ≦ $max_num_true_causes \n", color=:blue)
    else
        max_num_true_causes = length(S)
    end
    if iterate_all
        selection_only = false
    end
    accepted_sets = Dict{Union{Vector{Symbol},Vector{Int64}}, Float64}()
    n_tested_sets = 0
    running_intersection = S
    running_confintervals = zeros(p, 2)
    running_confintervals[:, 1] .= Inf
    running_confintervals[:, 2] .= -Inf
    _trace_confint_max = Vector{Vector{Float64}}()
    _trace_confint_min = Vector{Vector{Float64}}()
    trace_p_values = Vector{Float64}()
    n_env = maximum(env)
    if method == "chow"
        for i in 1:n_env
            ni = sum(env.==i)
            if ni > n_max_for_exact
                printstyled(@sprintf("environment %d has %d obs, subsample of %d is used\n", i, ni, n_max_for_exact), color=:blue)
            end
        end
    end
    println(@sprintf "Causal invariance search across %d environments with at α=%s (|S| = %d, method = %s, model = %s)\n" n_env α length(S) method model)
    max_num_true_causes < length(S) && printstyled("|S| is restricted to subsets with size ≦ $max_num_true_causes.\n", color=:blue)
    # priority queue: S -> -p.value (so sets with higher p-values are tested sooner)
    candidate_sets = PriorityQueue{Union{Set{Int64}, Set{Symbol}}, Float64}()
    enqueue!(candidate_sets, Set{typeof(S[1])}()=>0.)
    size_current_set = 0
    while size_current_set <= max_num_true_causes
        base_sets = Dict{Union{Set{Int64}, Set{Symbol}}, Float64}()  # set -> (-p_value)
        while !isempty(candidate_sets)
            # dequeue the set with highest p-value
            _S, neg_p_value = dequeue_pair!(candidate_sets)
            _S_vec = collect(_S)  # convert from Set{Int64} to Vector{Int64}
            # skip supersets under `selection_only`
            selection_only && length(running_intersection)!=length(S) && issubset(running_intersection, _S) && continue
            # test conditional invariance
            n_tested_sets += 1
            if method == "chow"
                rej, p_value, conf_intervals = conditional_inv_test_chow(X[:,_S_vec], y, env, n_env, α=α, n_max_for_exact=n_max_for_exact)
            elseif method == "logistic-LR"
                # target is the last column of df
                rej, p_value, conf_intervals = conditional_inv_test_logistic(df, target, _S_vec, env, n_env, α=α,
                                                                             add_intercept=true, method="logistic-LR")
            elseif method == "logistic-SF"
                rej, p_value, conf_intervals = conditional_inv_test_logistic(df, target, _S_vec, env, n_env, α=α,
                                                                             add_intercept=true, method="logistic-SF")
            end
            base_sets[_S] = -p_value
            push!(trace_p_values, p_value)
            if !rej
                # _S is an invariant set
                accepted_sets[_S_vec] = p_value
                # running ∩
                running_intersection = intersect(running_intersection, _S_vec)
                # running ∪
                conf_intervals_padded = zeros(size(running_confintervals))  # unincluded variables have [0,0] as confint
                if isa(_S_vec, Vector{Int64})
                    conf_intervals_padded[_S_vec, :] = conf_intervals
                else
                    _idx_vec = [columnindex(df, z) for z in _S_vec]
                    conf_intervals_padded[_idx_vec, :] = conf_intervals
                end
                running_confintervals[:, 1] = min.(running_confintervals[:, 1], conf_intervals_padded[:, 1])
                running_confintervals[:, 2] = max.(running_confintervals[:, 2], conf_intervals_padded[:, 2])
                # keep trace
                push!(_trace_confint_min, conf_intervals_padded[:, 1])
                push!(_trace_confint_max, conf_intervals_padded[:, 2])
                if verbose
                    println(@sprintf "S = %-40s: p-value = %-1.4f [%1s] ⋂ = %s" (isempty(_S_vec) ? "[]" : _S_vec) p_value (rej ? " " : "*") running_intersection)
                end
                if isempty(running_intersection) && (!iterate_all)
                    break   # empty, terminate
                end
            else
                if verbose
                    println(@sprintf "S = %-40s: p-value = %-1.4f [%1s] ⋂ = %s" (isempty(_S_vec) ? "[]" : _S_vec) p_value (rej ? " " : "*") running_intersection)
                end
            end
        end
        isempty(running_intersection) && (!iterate_all) && break  # have to break twice
        # generate sets with size ++
        for (_base_S, _neg_p_value) in base_sets
            for i in setdiff(S, _base_S)
                S_new = union(_base_S, Set([i]))
                # skipping supersets of running intersection in `selection_only` mode
                selection_only && length(running_intersection)!=length(S) && issubset(running_intersection, S_new) && continue
                if haskey(candidate_sets, S_new)
                    candidate_sets[S_new] = min(candidate_sets[S_new], _neg_p_value)
                else
                    candidate_sets[S_new] = _neg_p_value
                end
            end
        end
        size_current_set += 1
    end
    println("\nTested $n_tested_sets sets: $(length(accepted_sets)) sets are accepted.")
    trace_confint_min = hcat(_trace_confint_min...)'
    trace_confint_max = hcat(_trace_confint_max...)'
    if isempty(accepted_sets)
        # model rejected
        return CausalSearchResult(trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only = selection_only)
    elseif selection_only
        return CausalSearchResult(running_intersection, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only = true)
    else
        return CausalSearchResult(running_intersection, running_confintervals, trace_confint_min, trace_confint_max, trace_p_values, α, p, variables_considered; selection_only = false)
    end
end

function causalSearch(df::DataFrame, target::Union{Int64, Symbol}, env::Vector{Int64};
                      α=0.01, method="chow", screen="auto", p_max=8, verbose=true,
                      selection_only=false, iterate_all=false,
                      n_max_for_exact=5000, max_num_true_causes=Inf)
    if isa(target, Int64)
        target = propertynames(df)[target]
    end
    S = setdiff(propertynames(df), [target])
    X = df[:, S]
    y = df[!, target]
    causalSearch(X, y, env, S, α=α, method=method, screen=screen, p_max=p_max, verbose=verbose,
                 selection_only=selection_only, iterate_all=iterate_all,
                 n_max_for_exact=n_max_for_exact, max_num_true_causes=max_num_true_causes)
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
# f(x,y)=x'+y

FXtransposePY(x::Float64,y::Float64)=Fxpy(x,y)
FXtransposePY(x::Float64,y)=Fxpy(x,y)
FXtransposePY(x,y::Float64)=Fxpy(x',y)
FXtransposePY(x,y)=Fxpy(x',y)
FXtransposePY_inplace(value,auxvalue,x,y)=Fxpy_inplace(value,auxvalue,x',y)

function DXtransposePY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c',grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    end
end

if PROC=="GPU"
    function FXtransposePY_inplace(value::CudaArray,auxvalue,x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            gfill!(value,x); axpy!(1.0,y,value)
        elseif size(y)==(1,1)
            gfill!(value,y)
            CUBLAS.geam!('T','N',1.0,x,1.0,value,value)
        else
            CUBLAS.geam!('T','N',1.0,x,1.0,y,value)
        end
    end

    function DXtransposePY(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            if size(x)==(1,1)
                axpy!(1.0,sum(grad_c),grad_n)
            else
                tmp=CudaArray(Float64,size(grad_c))
                CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
                axpy!(1.0,tmp,grad_n); free(tmp)
            end
        elseif derivativeIDX==2
            if size(y)==(1,1)
                axpy!(1.0,sum(grad_c),grad_n)
            else
                axpy!(1.0,grad_c,grad_n)
            end
        end
    end
end

Derivative[FXtransposePY]=DXtransposePY
Inplace[FXtransposePY]=FXtransposePY_inplace

export FXtransposePY

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Printf

include("ylmr2.jl")

function test_main()

  Ng = 5
  lmax = 3

  Gvec = zeros(3,Ng)

  for ig = 1:Ng
    Gvec[:,ig] = [-2.1, 2.0, 1.0]
  end

  ylm = ylmr2( lmax, Gvec )

  for lm = 1:(lmax+1)^2
    @printf("\n")
    @printf("lm = %5d\n", lm)
    for ig = 1:Ng
      @printf("%5d %18.10f\n", ig, ylm[ig,lm])
    end
  end

end

test_main()


text:
module LBFGSB

using Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("LBFGSB not installed properly, run Pkg.build(\"LBFGSB\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

include("subroutine.jl")
export setulb, timer

include("wrapper.jl")
export L_BFGS_B

end # module


text:
module EMSOModelLibrary

# package code goes here
  export controllers,costs,electrical,heat_exchangers,mixers_splitters,pressure_changers,reactors
  export stage_separators,streams,types,UserModels,power_plant

  module types
    export EMLtypes
    include("../JuliaEMSOModels/types.jl")
  end

  module streams
    export EMLstreams
    include("../JuliaEMSOModels/streams.jl")
  end

  module controllers
    export EMLComparator,EMLHiLoSelect,EMLiae,EMLise
    export EMLlag_1,EMLlead_lag,EMLmultiply
    export EMLPIDIncr,EMLPIDs,EMLratio,EMLsum
    include("../JuliaEMSOModels/controllers/Comparator.jl");
    include("../JuliaEMSOModels/controllers/flow_measurement.jl");
    include("../JuliaEMSOModels/controllers/heat_flow.jl");
    include("../JuliaEMSOModels/controllers/HiLoSelect.jl");
    include("../JuliaEMSOModels/controllers/iae.jl");
    include("../JuliaEMSOModels/controllers/ise.jl");
    include("../JuliaEMSOModels/controllers/lag_1.jl");
    include("../JuliaEMSOModels/controllers/lead_lag.jl");
    include("../JuliaEMSOModels/controllers/multiply.jl");
    include("../JuliaEMSOModels/controllers/pi.jl");
    include("../JuliaEMSOModels/controllers/PIDIncr.jl");
    include("../JuliaEMSOModels/controllers/PIDs.jl");
    include("../JuliaEMSOModels/controllers/ratio.jl");
    include("../JuliaEMSOModels/controllers/sum.jl");
  end

  module mixers_splitters
    export EMLmixer,EMLsepComp,EMLsplitter
    include("../JuliaEMSOModels/mixers_splitters/mixer.jl")
    include("../JuliaEMSOModels/mixers_splitters/sepComp.jl")
    include("../JuliaEMSOModels/mixers_splitters/splitter.jl")
  end

  module pressure_changers
    export EMLcompressor,EMLexpander,EMLpump,EMLturbine,EMLvalve
    include("../JuliaEMSOModels/pressure_changers/compressor.jl")
    include("../JuliaEMSOModels/pressure_changers/expander.jl")
    include("../JuliaEMSOModels/pressure_changers/pump.jl")
    include("../JuliaEMSOModels/pressure_changers/turbine.jl")
    include("../JuliaEMSOModels/pressure_changers/valve.jl")
  end

  module stage_separators
    export EMLbatch_dist,EMLcolumn,EMLcolumn_Eff,EMLcondenser
    export EMLflash,EMLreboiler,EMLtank,EMLtray,EMLtray_Eff
    include("../JuliaEMSOModels/stage_separators/batch_dist.jl")
    include("../JuliaEMSOModels/stage_separators/tray.jl")
    include("../JuliaEMSOModels/stage_separators/tank.jl")
    include("../JuliaEMSOModels/stage_separators/reboiler.jl")
    include("../JuliaEMSOModels/stage_separators/condenser.jl")
    include("../JuliaEMSOModels/stage_separators/column.jl")
    include("../JuliaEMSOModels/stage_separators/flash.jl")
  end

  module heat_exchangers
    export EMLDoublePipe,EMLDoublePipeIncr,EMLHairpin,EMLHairpinIncr
    export EMLheater,EMLHeatex,EMLHeatExchangerDetailed,EMLHEX_Engine,EMLPHE
    include("../JuliaEMSOModels/heat_exchangers/HEX_Engine.jl")
    include("../JuliaEMSOModels/heat_exchangers/DoublePipe.jl")
    include("../JuliaEMSOModels/heat_exchangers/DoublePipeIncr.jl")
    include("../JuliaEMSOModels/heat_exchangers/Hairpin.jl")
    include("../JuliaEMSOModels/heat_exchangers/HairpinIncr.jl")
    include("../JuliaEMSOModels/heat_exchangers/heater.jl")
    include("../JuliaEMSOModels/heat_exchangers/Heatex.jl")
    include("../JuliaEMSOModels/heat_exchangers/HeatExchangerDetailed.jl")
    include("../JuliaEMSOModels/heat_exchangers/PHE.jl")
  end

  module costs
    export EMLcolumn_cost,EMLflash_cost,EMLHeatExchangerDetailed_cost
    export EMLHeatExchangerSimplified_cost,EMLreboiler_cost,EMLtank_cost
    include("../JuliaEMSOModels/costs/column_cost.jl")
    include("../JuliaEMSOModels/costs/flash_cost.jl")
    include("../JuliaEMSOModels/costs/HeatExchangerDetailed_cost.jl")
    include("../JuliaEMSOModels/costs/HeatExchangerSimplified_cost.jl")
    include("../JuliaEMSOModels/costs/reboiler_cost.jl")
    include("../JuliaEMSOModels/costs/tank_cost.jl")
  end

  module electrical
    export EMLelectrical
    include("../JuliaEMSOModels/electrical/electrical.jl")
  end

  module reactors
    export EMLbatch,EMLcstr,EMLequil,EMLgibbs,EMLpfr
    export EMLstoic,EMLtank_basic,EMLvol_tank,EMLyield
    include("../JuliaEMSOModels/reactors/batch.jl")
    include("../JuliaEMSOModels/reactors/cstr.jl")
    include("../JuliaEMSOModels/reactors/vol_tank.jl")
    include("../JuliaEMSOModels/reactors/tank_basic.jl")
    include("../JuliaEMSOModels/reactors/equil.jl")
    include("../JuliaEMSOModels/reactors/gibbs.jl")
    include("../JuliaEMSOModels/reactors/pfr.jl")
    include("../JuliaEMSOModels/reactors/stoic.jl")
    include("../JuliaEMSOModels/reactors/yield.jl")
  end

  module power_plant
    export EMLpower_plant
    include("../JuliaEMSOModels/water_steam/power_plant.jl")
  end

  module UserModels
    export EMLUserTray,EMLUserColumn
    include("../JuliaEMSOModels/UserModels/UserTray.jl")
    #include("../JuliaEMSOModels/UserModels/UserColumn.jl")
  end
end # module


text:
function BasicConstantThrustPropulsionModel(arg0::jdouble, arg1::jdouble, arg2::Vector3D, arg3::JString)
    return BasicConstantThrustPropulsionModel((jdouble, jdouble, Vector3D, JString), arg0, arg1, arg2, arg3)
end

function get_flow_rate(obj::BasicConstantThrustPropulsionModel)
    return jcall(obj, "getFlowRate", jdouble, ())
end

function get_flow_rate(obj::BasicConstantThrustPropulsionModel, arg0::Vector{RealFieldElement})
    return jcall(obj, "getFlowRate", RealFieldElement, (Vector{RealFieldElement},), arg0)
end

function get_flow_rate(obj::BasicConstantThrustPropulsionModel, arg0::Vector{jdouble})
    return jcall(obj, "getFlowRate", jdouble, (Vector{jdouble},), arg0)
end

function get_parameters_drivers(obj::BasicConstantThrustPropulsionModel)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_thrust_vector(obj::BasicConstantThrustPropulsionModel)
    return jcall(obj, "getThrustVector", Vector3D, ())
end

function get_thrust_vector(obj::BasicConstantThrustPropulsionModel, arg0::Vector{RealFieldElement})
    return jcall(obj, "getThrustVector", FieldVector3D, (Vector{RealFieldElement},), arg0)
end

function get_thrust_vector(obj::BasicConstantThrustPropulsionModel, arg0::Vector{jdouble})
    return jcall(obj, "getThrustVector", Vector3D, (Vector{jdouble},), arg0)
end


text:
using Parameters

"""
    TestFunction

Data structure to store test functions.
"""
@with_kw struct TestFunction
    fun::Function
    sr::Array{Tuple{Float64,Float64},1}
    min_val::Float64
    min_loc::Array{Float64,2}
    max_val::Float64
    max_loc::Array{Float64,2}
end

function rosenbrock_2D(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function rotatedHyperElipsoid_2D(x)
    out = 0.0
    xy = [x[1],x[2]]
    for i = 1:2
        inner = 0.0
        for j = 1:i
            inner = inner + xy[j]^2
        end
        out += inner
    end
    return out
end

function styblinskiTang_2D(x)
    out = 0.0
    xy = [x[1],x[2]]
    for i = 1:2
        out += xy[i]^4 - 16*xy[i]^2 + 5*xy[i]
    end
    out *= 0.5
    return out
end

function threeHumpCamel_2D(x)
    out = 2*x[1]^2 - 1.05*x[1]^4 + (x[1]^6)/6 + x[1]*x[2] + x[2]^2
end

function bohachevsky_2D(x)
    out = x[1]^2 + 2*x[2]^2 - 0.3*cos(3*π*x[1]) - 0.4*cos(4*π*x[2]) + 0.7
end

function permfunc_2D(x)
    out = 0.0
    xy = x
    for i = 1:2
        inner = 0.0
        for j = 1:2
            xj = xy[j]
            inner = inner + (j)*(xj^i-(1/j)^i)
        end
        out = out + inner^2
    end
    return out
end

function sphereFunction_2D(x)
    out = 0.0
    xy = x
    for i = 1:2
        out += xy[i]^2
    end
    return out
end

function sumOfDifferentPowers_2D(x)
    out = 0.0
    xy = x
    for i = 1:2
        out += abs(xy[i])^(i+1)
    end
    return out
end

function sumOfSquares_2D(x)
    out = 0.0
    xy = x
    for i = 1:2
        out += i*xy[i]^2
    end
    return out
end

function trid_2D(x)
    xy = x
    out = (xy[1]-1)^2
    inner = 0.0
    for i = 2:2
        out += (xy[i]-1)^2
        inner += xy[i]*xy[i-1]
    end
    return out - inner
end

function booth_2D(x)
    out = (x[1]+2*x[2]-1)^2 + (2*x[1]+x[2]-5)^2
end

function matyas_2D(x)
    out = 0.26*(x[1]^2+x[2]^2) - 0.48*x[1]*x[2]
end

function mccormick_2D(x)
    out = sin(x[1]+x[2]) + (x[1]-x[2])^2 - 1.5*x[1] + 2.5*x[2] + 1
end

function powersum_2D(x)
    out = 0.0
    b = [8.0,18.0]
    xy = x
    for i = 1:2
        inner = 0.0
        for j = 1:2
            inner += xy[j]^i
        end
        out += (inner-b[i])^2
    end
    return out
end

function zakharov_2D(x)
    out = 0.0
    tmp = 0.0
    xy = x
    for i = 1:2
        out += xy[i]^2
        tmp += 0.5*i*xy[i]
    end
    return out + tmp^2 + tmp^4
end

function sixHumpCamel_2D(x)
    out = (4-2.1*x[1]^2 + (x[1]^4)/3)*x[1]^2 + x[1]*x[2] + (-4 + 4*x[2]^2)*x[2]^2
end

function dixonPrice_2D(x)
    out = 0.0
    xy = x
    out += (x[1]-1)^2
    for i = 2:2
        out += i*(2*xy[i]^2-xy[i]-1)^2
    end
    return out
end

function beale_2D(x)
    out = (1.5-x[1]+x[1]*x[2])^2 + (2.25-x[1]+x[1]*x[2]^2)^2 + (2.625-x[1]+x[1]*x[2]^3)^2
end

function branin_2D(x)
    out = 1*(x[2] - 5.1/(4*π^2)*x[1]^2 + 5/π*x[1] - 6)^2 + 10*(1-1/8/π)*cos(x[1]) + 10
end

function goldsteinPrice_2D(x)
    out = (1 + ((x[1]+x[2]+1)^2)*(19-14*x[1]+3*x[1]^2-14*x[2]+6*x[1]*x[2]+3*x[2]^2)) *
                      (30 + ((2*x[1]-3*x[2])^2)*(18-32*x[1]+12*x[1]^2+48*x[2]-36*x[1]*x[2]+27x[2]^2))
end

function permdbeta_2D(x)
    out = 0.0
    xy = x
    for i = 1:2
        inner = 0.0
        for j = 1:2
            inner += (j^i+0.5)*((xy[j]/j)^i-1)
        end
        out += inner^2
    end
    return out
end

function rastrigin_2D(x)
    sum = 0
    for ii = 1:2
        xi = x[ii]
        sum = sum + (xi^2 - 10*cos(2*pi*xi))
    end
    out = 10*2 + sum
end

function hart_4D(x)

    alpha = [1.0 1.2 3.0 3.2]
    A = [10 3 17 3.5 1.7 8;
        0.05 10 17 0.1 8 14;
        3 3.5 1.7 10 17 8;
        17 8 0.05 10 0.1 14]
    P = 10^(-4) * [1312 1696 5569 124 8283 5886;
                2329 4135 8307 3736 1004 9991;
                2348 1451 3522 2883 3047 6650;
                4047 8828 8732 5743 1091 381]

    outer = 0
    for ii = 1:4
        inner = 0
        for jj = 1:4
            xj = x[jj]
            Aij = A[ii, jj]
            Pij = P[ii, jj]
            inner += Aij*(xj-Pij)^2
        end
        outer += alpha[ii] * exp(-inner)
    end

    out = (1.1 - outer) / 0.839

end

function hart_6D(x)
    alpha = [1.0 1.2 3.0 3.2]
    A = [10 3 17 3.5 1.7 8;
         0.05 10 17 0.1 8 14;
         3 3.5 1.7 10 17 8;
         17 8 0.05 10 0.1 14];
    P = 10^(-4) * [1312 1696 5569 124 8283 5886;
                   2329 4135 8307 3736 1004 9991;
                   2348 1451 3522 2883 3047 6650;
                   4047 8828 8732 5743 1091 381]

    outer = 0
    for ii = 1:4
        inner = 0
        for jj = 1:6
            xj = x[jj]
            Aij = A[ii, jj]
            Pij = P[ii, jj]
            inner = inner + Aij*(xj-Pij)^2
        end
        new = alpha[ii] * exp(-inner)
        outer = outer + new
    end

    y = -outer
end

function rosenbrock_ND(x)
    d = length(x)
    sum = 0.0
    for ii = 1:(d-1)
        xi = x[ii]
        xnext = x[ii+1]
        sum += 100*(xnext-xi^2)^2 + (xi-1)^2
    end
    out = sum
end

##### Note, the maximums are approximate. The functions may have several optima, only one is stored.
test_funs = Dict(
    :rosenbrock_2D => TestFunction( fun = rosenbrock_2D,
                                    sr = [(-5.0,5.0),(-5.0,5.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.0; 1.0]'),
                                    max_val = 90036.0,
                                    max_loc = permutedims([-5.0; -5.0]')),

    :rotatedHyperElipsoid_2D => TestFunction( fun = rotatedHyperElipsoid_2D,
                                    sr = [(-65.536, 65.536), (-65.536, 65.536)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 12884.901888,
                                    max_loc = permutedims([65.536; -65.536]')),

    :styblinskiTang_2D => TestFunction( fun = styblinskiTang_2D,
                                    sr = [(-5.0, 5.0), (-5.0, 5.0)],
                                    min_val = -78.33233140754285,
                                    min_loc = permutedims([-2.90353; -2.90353]'),
                                    max_val = 250.0,
                                    max_loc = permutedims([5.0; 5.0]')),

    :threeHumpCamel_2D => TestFunction( fun = threeHumpCamel_2D,
                                    sr = [(-5.0, 5.0), (-5.0, 5.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 2047.9166666666665,
                                    max_loc = permutedims([5.0; 5.0]')),

    :bohachevsky_2D => TestFunction( fun = bohachevsky_2D,
                                    sr = [(-100.0, 100.0), (-100.0, 100.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 30000.0,
                                    max_loc = permutedims([-100.0; -100.0]')),

    :permfunc_2D => TestFunction( fun = permfunc_2D,
                                    sr = [(-2.0, 2.0), (-2.0, 2.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.0; 0.5]'),
                                    max_val = 174.25,
                                    max_loc = permutedims([-2.0; -2.0]')),

    :sphereFunction_2D => TestFunction( fun = sphereFunction_2D,
                                    sr = [(-5.12, 5.12), (-5.12, 5.12)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 52.4288,
                                    max_loc = permutedims([5.12; 5.12]')),

    :sumOfDifferentPowers_2D => TestFunction( fun = sumOfDifferentPowers_2D,
                                    sr = [(-1.0, 1.0), (-1.0, 1.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 2.0,
                                    max_loc = permutedims([1.0; -1.0]')),

    :sumOfSquares_2D => TestFunction( fun = sumOfSquares_2D,
                                    sr = [(-10.0, 10.0), (-10.0, 10.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0;0.0]'),
                                    max_val = 300.0,
                                    max_loc = permutedims([-10.0; 10.0]')),

    :trid_2D => TestFunction( fun = trid_2D,
                                    sr = [(-4.0, 4.0), (-4.0, 4.0)],
                                    min_val = -2.0,
                                    min_loc = permutedims([2.0; 2.0]'),
                                    max_val = 50.0,
                                    max_loc = permutedims([-4.0; 4.0]')),

    :booth_2D => TestFunction( fun = booth_2D,
                                    sr = [(-10.0, 10.0), (-10.0, 10.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([3.0; -1.0]'),
                                    max_val = 2186.0,
                                    max_loc = permutedims([-10.0; -10.0]')),

    :matyas_2D => TestFunction( fun = matyas_2D,
                                    sr = [(-10.0, 10.0), (-10.0, 10.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 100.0,
                                    max_loc = permutedims([10.0; -10.0]')),

    :mccormick_2D => TestFunction( fun = mccormick_2D,
                                    sr = [(-1.5, 4.0), (-3.0, 4.0)],
                                    min_val = -1.9132229549810367,
                                    min_loc = permutedims([-0.547198; -1.5472]'),
                                    max_val = 44.09847214410396,
                                    max_loc = permutedims([-1.5; 4.0]')),

    :powersum_2D => TestFunction( fun = powersum_2D,
                                    sr = [(0.0, 2.0), (-0.0, 2.0)],
                                    min_val = 116.0,
                                    min_loc = permutedims([2.0; 2.0]'),
                                    max_val = 388.0,
                                    max_loc = permutedims([0.0;0.0]')),

    :zakharov_2D => TestFunction( fun = zakharov_2D,
                                    sr = [(-5.0, 10.0), (-5.0, 10.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([0.0; 0.0]'),
                                    max_val = 51050.0,
                                    max_loc = permutedims([10.0; 10.0]')),

    :sixHumpCamel_2D => TestFunction( fun = sixHumpCamel_2D,
                                    sr = [(-3.0, 3.0), (-2.0, 2.0)],
                                    min_val = -1.0316284534898774,
                                    min_loc = permutedims([-0.089842; 0.712656]'),
                                    max_val = 162.89999999999998,
                                    max_loc = permutedims([3.0; 2.0]')),

    :dixonPrice_2D => TestFunction( fun = dixonPrice_2D,
                                    sr = [(-10.0, 10.0), (-10.0, 10.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.0; -0.5]'),
                                    max_val = 87483.0,
                                    max_loc = permutedims([-10.0; -10.0]')),

    :beale_2D => TestFunction( fun = beale_2D,
                                    sr = [(-4.5, 4.5), (-4.5, 4.5)],
                                    min_val = 0.0,
                                    min_loc = permutedims([3.0; -0.5]'),
                                    max_val = 181853.61328125,
                                    max_loc = permutedims([-4.5; -4.5]')),

    :branin_2D => TestFunction( fun = branin_2D,
                                    sr = [(-5.0, 10.0), (0.0, 15.0)],
                                    min_val = 0.39788735772973816,
                                    min_loc = permutedims([3.14159; 2.275]'),
                                    max_val = 308.12909601160663,
                                    max_loc = permutedims([-5.0; 1.12981e-15]')),

    :goldsteinPrice_2D => TestFunction( fun = goldsteinPrice_2D,
                                    sr = [(-2.0, 2.0), (-2.0, 2.0)],
                                    min_val = 3,
                                    min_loc = permutedims([0; -1.0]'),
                                    max_val = 1.0156902717980599e6,
                                    max_loc = permutedims([-1.73737; 2.0]')),

    :permdbeta_2D => TestFunction( fun = permdbeta_2D,
                                    sr = [(-2.0, 2.0), (-2.0, 2.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.0; 2.0]'),
                                    max_val = 110.5,
                                    max_loc = permutedims([-2.0; -2.0]')),

    :rastrigin_2D => TestFunction( fun = rastrigin_2D,
                                    sr = [(-5.12, 5.12), (-5.12, 5.12)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.05917e-9; 1.72311e-9]'),
                                    max_val = 80.70658038767792,
                                    max_loc = permutedims([-4.52299; 4.52299]')),

    :hart_4D => TestFunction( fun = hart_4D,
                                    sr = [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],
                                    min_val = -3.1344941412224,
                                    min_loc = permutedims([0.187395; 0.194152; 0.557918; 0.26478]'),
                                    max_val = 1.3095406214559708,
                                    max_loc = permutedims([1.0; 1.0; 1.42362e-13; 1.0]')),

    :hart_6D => TestFunction( fun = hart_6D,
                                    sr = [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],
                                    min_val = -3.322368011415515,
                                    min_loc = permutedims([0.20169; 0.150011; 0.476874; 0.275332; 0.311652; 0.657301]'),
                                    max_val = -2.8124505439686604e-8,
                                    max_loc = permutedims([1.0; 1.0; 3.13887e-17; 1.0; 1.0; 1.0]')),

    :rosenbrock_9D => TestFunction( fun = rosenbrock_ND,
                                    sr = [(-5.0, 5.0) for i in 1:9],
                                    min_val = 0.0,
                                    min_loc = [1.0 for i in 1:9, j in 1:1],
                                    max_val = 720288.0,
                                    max_loc = [-5.0 for i in 1:9, j in 1:1]),

    :rosenbrock_12D => TestFunction( fun = rosenbrock_ND,
                                    sr = [(-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0)],
                                    min_val = 0.0,
                                    min_loc = permutedims([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0]'),
                                    max_val = 990396.0,
                                    max_loc = permutedims([-5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0; -5.0]')),
                                    );


text:
using SparseArrays
using MAT
using MatrixNetworks
using LinearAlgebra

# Controlled growth from a seed set R in a hypergraph. Look at the one hop
# neighborhood, and order all of those nodes by how many hyperedges they are in
# that have nodes from R. Order that last, and take the top k
function TopNeighbors(H::SparseMatrixCSC{Float64,Int64},R::Vector{Int64},R1hop::Vector{Int64},k::Int64)

    if length(R1hop) > k
        # Get all edges touching R
        HR = H[:,R]
        rp = HR.rowval
        edges = unique(rp)

        # Consider how many touch the 1-hop neighborhood
        HL = H[edges,R1hop]

        # For each node in R1hop, compute the number of edges it has that touch R
        d2R = vec(sum(HL,dims=1))

        # order = sortperm(d2R, rev=true)
        b = partialsortperm(d2R, 1:k, rev=true)
        Rmore = R1hop[b]
    else
        Rmore = R1hop
    end

    return union(R, Rmore)
end

# Controlled growth from a seed set R in a hypergraph. Look at the one hop
# neighborhood, and order all of those nodes by what percent of their
# edges touch R. Order that, and take the top k
function BestNeighbors(H::SparseMatrixCSC{Float64,Int64},d::Vector{Float64},R::Vector{Int64},R1hop::Vector{Int64},k::Int64)

    if length(R1hop) > k
        # Get all edges touching R
        HR = H[:,R]
        rp = HR.rowval
        edges = unique(rp)

        # Consider how many touch the 1-hop neighborhood
        HL = H[edges,R1hop]

        # For each node in R1hop, compute the number of edges it has that touch R
        d1 = d[R1hop]
        d2 = vec(sum(HL,dims=1))

        # order = sortperm(d2R, rev=true)
        b = partialsortperm(d2./d1, 1:k, rev=true)
        Rmore = R1hop[b]
    else
        Rmore = R1hop
    end
    return union(R, Rmore)
end

# Simple function for returning ALL the indices where we find a maximum
function findallmax(v)

    l = length(v)
    m = minimum(v)
    M = maximum(v)
    if M == m
        return collect(1:l)
    else
        Inds = Vector{Int64}()
        a,b = findmax(v)
        while a == M
            push!(Inds,b)
            v[b] = m
            a,b = findmax(v)
        end
        return Inds
    end

end

## Delta-Linear (tl = thresholded linear) conductance computation.
# e.g. tl_cond(H,S,d,delta,volA,order)
function tl_cond(H::SparseMatrixCSC,S::Vector{Int64},d::Vector{Float64},delta::Float64,volA::Float64,order::Vector{Int64})

    if volA == 0.0
        volA = sum(d)
    end
    n = length(d)
    volS = sum(d[S])
    cut = tl_cut(H,S,delta,order)

    cond = cut/min(volS, volA-volS)

    return cond, volS, cut

end

## Delta-Linear (thresholded linear) normalized Cut computation.
# e.g. tl_ncut(H,S,d,delta,volA,order)
function tl_ncut(H::SparseMatrixCSC,S::Vector{Int64},d::Vector{Float64},delta::Float64,volA::Float64,order::Vector{Int64})

    if volA == 0.0
        volA = sum(d)
    end
    n = length(d)
    volS = sum(d[S])
    cut = tl_cut(H,S,delta,order)

    cond = cut/min(volS, volA-volS)
    ncut = cut/(volS) + cut/(volA-volS)

    # rncut = round(Int64,ncut)
    # rcut = round(Int64,cut)
    # rcond = round(cond,digits = 4)
    # rvol = round(Int64,volS)

    return cond, ncut, volS, cut

end

# Thresholded linear cut value for a set
# calling e.g. tl_cut(H,S,delta,order)
function tl_cut(H::SparseMatrixCSC{Float64,Int64}, S::Vector{Int64}, delta::Float64,order::Vector{Int64})

    # Check the cut
    HS = H[:,S]
    sumHS = sum(HS,dims = 2)  # Count number of S nodes in each hyperedge
    inds = findall(x->x>0,sumHS)    # Extract hyperedges with > 0 nodes from S
    ES = sumHS[inds]
    verts = order[inds]               # Get the size of these hyperedges

    # Find number of nodes on small side of cut
    SmallSide = round.(Int64,min.(ES, verts-ES))
    # Compute the cardinality-based cut score
    cutval = 0.0
    for j = 1:length(SmallSide)
        sm = SmallSide[j]
        if sm > 0
            if sm < delta
                cutval += sm
            else
                cutval += delta
            end
        end
    end

    return cutval
end

# For a set S in a hypergraph, return the hypergraph local conductance
# score with thresholded linear splitting penalty
# e.g. hlc_tl(H,order,R,S,d,volA,epsilon,delta)
function hlc_tl(H::SparseMatrixCSC{Float64,Int64},order::Vector{Int64},R::Vector{Int64},
    S::Vector{Int64},d::Vector{Float64},volA::Float64,epsilon::Float64,
    delta::Float64,)

    volS = sum(d[S])
    RnS = intersect(R,S)
    volRnS = sum(d[RnS])
    cut = tl_cut(H,S,delta,order)

    lcond = cut/((1+epsilon)*volRnS - epsilon*volS)

    return lcond

end

# Expand a hypergraph using the thresholded linear splitting function.
#
#   Hyperedges = Hyperedge list
#   delta = TL splitting function parameter
function tl_expansion(Hyperedges::Vector{Vector{Int64}}, order::Vector{Int64}, delta::Float64,n::Int64)

        BigEdges = length(findall(x->x>3,order))
        N = n + 2*BigEdges

        ## Build the adjacency matrix
        ap = n+1   # "auxiliary node pointer", points to next "free" aux node

        # Build the sparse matrix
        U = Vector{Int64}()
        V = Vector{Int64}()
        vals = Vector{Float64}()

        for edge = Hyperedges
            nv = length(edge)
            if nv == 2
                i = edge[1]; j = edge[2]
                #A[i,j] += 1; A[j,i] += 1
                push!(U,i); push!(V,j); push!(vals,1)
                push!(U,j); push!(V,i); push!(vals,1)
            elseif nv == 3
                i = edge[1]; j = edge[2]; k = edge[3]
                # A[i,j] += 1/2; A[j,i] += 1/2
                # A[k,j] += 1/2; A[j,k] += 1/2
                # A[k,i] += 1/2; A[i,k] += 1/2
                push!(U,i); push!(V,j); push!(vals,1/2)
                push!(U,j); push!(V,i); push!(vals,1/2)
                push!(U,i); push!(V,k); push!(vals,1/2)
                push!(U,k); push!(V,i); push!(vals,1/2)
                push!(U,j); push!(V,k); push!(vals,1/2)
                push!(U,k); push!(V,j); push!(vals,1/2)
            else
                # We need to add auxiliary vertices
                for i = edge
                    # A[i,auxpointer] = 1
                    # A[auxpointer+1,i] = 1
                    # A[auxpointer,auxpointer+1] = w2
                    push!(U,i); push!(V,ap); push!(vals,1)
                    push!(U,ap+1); push!(V,i); push!(vals,1)
                end
                push!(U,ap); push!(V,ap+1); push!(vals,delta)
                ap += 2
            end

        end
        @show maximum(U), maximum(V), N
        A = sparse(U,V,vals,N,N)
        return A
end

# Given an incidence matrix for a hypergraph and its transpose (having both
# handy makes different parts of the code faster), and a set of nodes R,
# return the immediate neighbors of R that don't include R itself
function get_immediate_neighbors(H::SparseMatrixCSC{Float64,Int64},
    Ht::SparseMatrixCSC{Float64,Int64},R::Vector{Int64})

    Hr = H[:,R]
    rp_r = Hr.rowval
    R_edges = unique(rp_r)

    He = Ht[:,R_edges]
    rp_e = He.rowval
    Rneighbs = unique(rp_e)
    Rn = setdiff(Rneighbs,R)

    return Rn

end

function neighborhood(H::SparseMatrixCSC{Float64,Int64},
    Ht::SparseMatrixCSC{Float64,Int64},R::Vector{Int64})
    Hr = H[:,R]
    rp_r = Hr.rowval
    R_edges = unique(rp_r)

    He = Ht[:,R_edges]
    rp_e = He.rowval
    Rn = unique(rp_e)

    return Rn
end

function neighborlist(H::SparseMatrixCSC{Float64,Int64},
    Ht::SparseMatrixCSC{Float64,Int64})

    Neighbs = Dict()
    n = size(H,2)
    t1 = 0
    t2 = 0
    t3 = 0
    for i = 1:n
        # s = time()
        ivec = H[:,i]
        #n_edges = findnz(ivec)[1]
        n_edges = ivec.nzind    # get neighboring edges
        # t1 += time()-s

        # s = time()
        He = Ht[:,n_edges]      # nodes touching those edges
        rp_e = He.rowval
        neighbs_i = unique(rp_e)
        # t2 += time()-s

        # s = time()
        push!(Neighbs,neighbs_i)
        # t3 += time()-s
    end
    # @show t1, t2, t3
    return Neighbs
end

# Expand a hypergraph using the thresholded linear splitting function.
#
#   H = |E| x |V| binary incidence matrix for the hypergraph
#   delta = TL splitting function parameter
function tl_expansion_inc(H::SparseMatrixCSC{Float64,Int64}, order::Vector{Int64}, delta::Float64)

        n = size(H,2)
        BigEdges = length(findall(x->x>3,order))
        N = n + 2*BigEdges

        Hyperedges = incidence2elist(H)

        ## Build the adjacency matrix
        ap = n+1   # "auxiliary node pointer", points to next "free" aux node

        # Build the sparse matrix
        U = Vector{Int64}()
        V = Vector{Int64}()
        vals = Vector{Float64}()

        for ee = 1:length(Hyperedges)
            edge = Hyperedges[ee]
            nv = length(edge)
            # if order[ee] != nv
            #     @show ee, nv, order[ee], edge
            # end
            if nv == 1
                # ignore
                # println("This")
            elseif nv == 2
                i = edge[1]; j = edge[2]
                #A[i,j] += 1; A[j,i] += 1
                push!(U,i); push!(V,j); push!(vals,1)
                push!(U,j); push!(V,i); push!(vals,1)
            elseif nv == 3
                i = edge[1]; j = edge[2]; k = edge[3]
                # A[i,j] += 1/2; A[j,i] += 1/2
                # A[k,j] += 1/2; A[j,k] += 1/2
                # A[k,i] += 1/2; A[i,k] += 1/2
                push!(U,i); push!(V,j); push!(vals,1/2)
                push!(U,j); push!(V,i); push!(vals,1/2)
                push!(U,i); push!(V,k); push!(vals,1/2)
                push!(U,k); push!(V,i); push!(vals,1/2)
                push!(U,j); push!(V,k); push!(vals,1/2)
                push!(U,k); push!(V,j); push!(vals,1/2)
            else
                # We need to add auxiliary vertices
                for i = edge
                    # A[i,auxpointer] = 1
                    # A[auxpointer+1,i] = 1
                    # A[auxpointer,auxpointer+1] = delta
                    push!(U,i); push!(V,ap); push!(vals,1)
                    push!(U,ap+1); push!(V,i); push!(vals,1)
                end
                push!(U,ap); push!(V,ap+1); push!(vals,delta)
                ap += 2
            end

        end
        # @show maximum(U), maximum(V), length(U), length(V), N, ap
        A = sparse(U,V,vals,N,N)
        return A
end

# Convert a hyperedge list to a hypergraph binary incidence matrix
function incidence2elist(Hin::SparseMatrixCSC{Float64,Int64})

    H = sparse(Hin')
    rp = H.rowval
    ci = H.colptr
    Hyperedges = Vector{Vector{Int64}}()
    n,m = size(H)

    for i = 1:m
        startedge = ci[i]
        endedge = ci[i+1]-1
        edge = rp[startedge:endedge]
        push!(Hyperedges,edge)
    end
    return Hyperedges
end

# Take a list of hyperedges and turn it into a hyperedge incidence matrix
# H. N is the number of nodes in the hypergraph.
#  H(e,u) = 1  iff node u is in hyperedge e
function elist2incidence(Hyperedges::Vector{Vector{Int64}}, N::Int64)
    U = Vector{Int64}()
    E = Vector{Int64}()
    M = length(Hyperedges)
    for enum = 1:length(Hyperedges)
        e = Hyperedges[enum]
        for node in e
            push!(U,node)
            push!(E,enum)
        end
    end

    H = sparse(E,U,ones(length(U)),M,N)
    return H
end

# This computes the precision, recall, and F1 score for a set Returned
# compared against a Target set
function PRF(Target,Returned)

    if length(Returned) == 0
        pr = 0; re = 0; F1 = 0
    else
        TruePos = intersect(Returned,Target)
        pr = length(TruePos)/length(Returned)
        re = length(TruePos)/length(Target)
        F1 = 2*(pr*re)/(pr+re)

        if length(TruePos) == 0
            F1 = 0
        end
    end

    return pr, re, F1

end

## Given a binary incidence matrix H for a hypergraph, find the one-hop
# neighborhood of a set of nodes S
function hyper_neighborhood(H::SparseMatrixCSC{Float64,Int64},S::Vector{Int64})

    A = H'*H
    n = size(A,1)
    for i = 1:n
        A[i,i] = 0
    end
    dropzeros!(A)
    return neighborhood(A,S,1)

end

## Given a binary incidence matrix H for a hypergraph, find the one-hop
# neighborhood of a set of nodes S, when considering only hyperedges with
# a maximum number of M nodes
function hyper_neighborhood(H::SparseMatrixCSC{Float64,Int64},S::Vector{Int64},order::Vector{Int64},M::Int64)

    good = findall(x->x<=M,order)
    H = H[good,:]

    ##
    A = H'*H
    n = size(A,1)
    for i = 1:n
        A[i,i] = 0
    end
    dropzeros!(A)
    return neighborhood(A,S,1)

end

## Simple Clique Expansion
#   A[i,j] = number of hyperedges nodes i and j share
function SimpleCliqueExp(H::SparseMatrixCSC{Float64,Int64})

    A = H'*H
    for i = 1:size(A,1)
        A[i,i] = 0.0
    end
    dropzeros!(A)
    return A
end

## Weighted Clique Expansion
#  When performing the clique expansion, for each hyperedge expanded into a
#   clique, multiply each edge in the expansion by 1/order(e)
function WeightedCliqueExpansion(H::SparseMatrixCSC{Float64,Int64}, order::Vector{Int64})

    m,n = size(H)
    I = Vector{Int64}()
    J = Vector{Int64}()
    vals = Vector{Float64}()
    Hyperedges = incidence2elist(H)
    for e = 1:m
        Edge = Hyperedges[e]
        Ord = order[e]
        for ii = 1:length(Edge)
            for jj = ii+1:length(Edge)
                i = Edge[ii]
                j = Edge[jj]
                push!(I,i); push!(J,j); push!(vals,1/Ord)
            end
        end
        if mod(e,10000)==0
            println("$e")
        end
    end

    A = sparse(I,J,vals,n,n)
    A = sparse(A+A')
    return A
end

max_stars_repo_name: nveldt/SparseCardDSFM

max_issues_repo_name: nveldt/SparseCardDSFM

max_forks_repo_name: nveldt/SparseCardDSFM


text:
struct HailstoneSeq{T<:Integer}
	start::T
end

Base.eltype(::HailstoneSeq{T}) where T = T

Base.start(hs::HailstoneSeq) = (-1, hs.start)
Base.done(::HailstoneSeq, state) = state == (1, 4)
function Base.next(::HailstoneSeq, state)
	_, s2 = state
	s1 = s2
	if iseven(s2)
		s2 = s2 ÷ 2
	else
		s2 = 3s2 + 1
	end
	return s1, (s1, s2)
end

function Base.length(hs::HailstoneSeq)
	r = 0
	for _ in hs
		r += 1
	end
	return r
end

function Base.show(io::IO, hs::HailstoneSeq)
	f5 = collect(Iterators.take(hs, 5))
	print(io, "HailstoneSeq(", join(f5, ", "), "...)")
end

hs = HailstoneSeq(27)
println("Collection of the Hailstone sequence from 27: $hs")
cl = collect(hs)
println("First 5 elements: ", join(cl[1:5], ", "))
println("Last 5 elements: ", join(cl[end-4:end], ", "))

Base.isless(h::HailstoneSeq, s::HailstoneSeq) = length(h) < length(s)
println("The number with the longest sequence under 100,000 is: ", maximum(HailstoneSeq.(1:100_000)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Base.Test
# using Lint
using Mocking
Mocking.enable(force=true)

using Playground
using FilePaths
using Memento

Memento.config("debug"; fmt="[ {level} ] {msg}")

TEST_DIR = cwd()

#TEST_TMP_DIR = mktempdir("../test/")
TEST_TMP_DIR = join(TEST_DIR, p"tmp")
TEST_PLAYGROUND_DIR = join(TEST_TMP_DIR, p"playground")
TEST_CONFIG = Config(DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)

mkdir(TEST_TMP_DIR; recursive=true, exist_ok=true)
mkdir(TEST_PLAYGROUND_DIR; recursive=true, exist_ok=true)

Playground.init(TEST_CONFIG)

# Order matters.
tests = [
    # "lint",
    "utils",
    "list",
    "parsing",
    "install",
    "create",
    "activate",
    "shell",
    "execute",
    "clean"
]

for t in tests
    tfile = string(t, ".jl")
    include(join(TEST_DIR, Path(tfile)))
end

@testset "main" begin
    test_main(args) = main(args, DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)

    test_main(["--debug", "install", "download", "0.5", "--labels", "julia-0.5"])
    test_main(["--debug", "install", "link", join(string(TEST_CONFIG.bin), "julia-0.5"), "--labels", "julia-stable-dir"])
    test_main(["--debug", "create"])
    test_main(["--debug", "exec", "ls -al", join(string(TEST_DIR), ".playground")])
    test_main(["--debug", "list"])
    test_main(["--debug", "clean"])
    test_main(["--debug", "rm", "--dir", join(string(TEST_DIR), ".playground")])
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(domain::StructuredGrid{T,N}, data::AbstractVector) where {N,T}
  X = coordinates(domain)

  if N == 1
    seriestype --> :path
    X[1,:], data
  elseif N == 2
    seriestype --> :scatter
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:]
  elseif N == 3
    seriestype --> :scatter
    aspect_ratio --> :equal
    marker_z --> data
    colorbar --> true
    X[1,:], X[2,:], X[3,:]
  else
    @error "cannot plot in more than 3 dimensions"
  end
end

@recipe function f(domain::StructuredGrid{T,N}) where {N,T}
  X  = coordinates(domain)
  sz = size(domain)

  seriescolor --> :black
  legend --> false

  if N == 1
    @series begin
      seriestype --> :scatterpath
      X[1,:], fill(zero(T), sz[1])
    end
  elseif N == 2
    aspect_ratio --> :equal
    linear = LinearIndices(sz)
    @series begin
      seriestype --> :scatter
      X[1,:], X[2,:]
    end
    for i in 1:sz[1]
      @series begin
        seriestype --> :path
        primary --> false
        inds = [linear[i,j] for j in 1:sz[2]]
        X[1,inds], X[2,inds]
      end
    end
    for j in 1:sz[2]
      @series begin
        seriestype --> :path
        primary --> false
        inds = [linear[i,j] for i in 1:sz[1]]
        X[1,inds], X[2,inds]
      end
    end
  elseif N == 3
    aspect_ratio --> :equal
    linear = LinearIndices(sz)
    @series begin
      seriestype --> :scatter
      X[1,:], X[2,:], X[3,:]
    end
    for i in 1:sz[1], j in 1:sz[2]
      @series begin
        seriestype --> :path
        primary --> false
        inds = [linear[i,j,k] for k in 1:sz[3]]
        X[1,inds], X[2,inds], X[3,inds]
      end
    end
    for i in 1:sz[1], k in 1:sz[3]
      @series begin
        seriestype --> :path
        primary --> false
        inds = [linear[i,j,k] for j in 1:sz[2]]
        X[1,inds], X[2,inds], X[3,inds]
      end
    end
    for j in 1:sz[2], k in 1:sz[3]
      @series begin
        seriestype --> :path
        primary --> false
        inds = [linear[i,j,k] for i in 1:sz[1]]
        X[1,inds], X[2,inds], X[3,inds]
      end
    end
  else
    @error "cannot plot in more than 3 dimensions"
  end
end


text:
module TestStructuredCospans
using Test

using Catlab, Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra
using Catlab.Graphs.BasicGraphs: AbstractGraph, TheoryGraph

# Structured cospans of C-sets
##############################

const OpenGraphOb, OpenGraph = OpenCSetTypes(Graph, :V)

# Directed fork as open graph with one input and two outputs.
g0 = Graph(4)
add_edges!(g0, [1,2,2], [2,3,4])
g = OpenGraph(g0, FinFunction([1],4), FinFunction([3,4],4))
@test apex(g) == g0
@test dom.(legs(g)) == [Graph(1), Graph(2)]
@test feet(g) == [FinSet(1), FinSet(2)]
@test dom(g) == OpenGraphOb(FinSet(1))
@test codom(g) == OpenGraphOb(FinSet(2))

# Opposite of previous graph.
h0 = Graph(4)
add_edges!(h0, [1,2,3], [3,3,4])
h = OpenGraph(h0, FinFunction([1,2],4), FinFunction([4],4))

# Structured multicospans.
g′ = OpenGraph(g0, FinFunction([1],4), FinFunction([3],4), FinFunction([4],4))
@test feet(g′) == fill(FinSet(1), 3)
@test bundle_legs(g′, [1, (2,3)]) == g

# Category
#---------

# Composition.
k = compose(g, h)
@test dom(k) == dom(g)
@test codom(k) == codom(h)
k0 = apex(k)
@test (nv(k0), ne(k0)) == (6, 6)
@test (src(k0), tgt(k0)) == ([1,2,2,3,4,5], [2,3,4,5,5,6])

# Identities.
a = OpenGraphOb(FinSet(3))
@test dom(id(a)) == a
@test codom(id(a)) == a
@test compose(g, id(codom(g))) == g
@test compose(id(dom(g)), g) == g

# Symmetric monoidal category
#----------------------------

a, b = OpenGraphOb(FinSet(3)), OpenGraphOb(FinSet(2))
@test otimes(a, b) == OpenGraphOb(FinSet(5))
@test munit(OpenGraphOb) == OpenGraphOb(FinSet(0))
@test dom(braid(a, b)) == a⊗b
@test codom(braid(b, a)) == b⊗a
@test force(braid(a,b)⋅braid(b,a)) == force(id(a⊗b))

k = otimes(g, h)
@test dom(k) == dom(g)⊗dom(h)
@test codom(k) == codom(g)⊗codom(h)
@test (nv(apex(k)), ne(apex(k))) == (8, 6)

# Hypergraph category
#--------------------

@test compose(create(a), mcopy(a)) == dunit(a)
@test compose(mmerge(a), delete(a)) == dcounit(a)

@test dom(dagger(g)) == codom(g)
@test codom(dagger(g)) == dom(g)

# Structured cospan of attributed C-sets
########################################

# Non-attributed boundary
#------------------------

const OpenWeightedGraphOb, OpenWeightedGraph = OpenACSetTypes(WeightedGraph, :V)

g0 = WeightedGraph{Float64}(2)
add_edge!(g0, 1, 2, weight=1.5)
g = OpenWeightedGraph{Float64}(g0, FinFunction([1],2), FinFunction([2],2))

h0 = WeightedGraph{Float64}(3)
add_edges!(h0, [1,1], [2,3], weight=[1.0,2.0])
h = OpenWeightedGraph{Float64}(h0, FinFunction([1],3), FinFunction([2,3],3))
@test dom.(legs(h)) == [WeightedGraph{Float64}(1), WeightedGraph{Float64}(2)]
@test feet(h) == [FinSet(1), FinSet(2)]
@test dom(h) == OpenWeightedGraphOb{Float64}(FinSet(1))
@test codom(h) == OpenWeightedGraphOb{Float64}(FinSet(2))

k = compose(g, h)
k0 = apex(k)
@test src(k0) == [1,2,2]
@test tgt(k0) == [2,3,4]
@test subpart(k0, :weight) == [1.5, 1.0, 2.0]

# Attributed boundary
#--------------------

@present TheoryVELabeledGraph <: TheoryGraph begin
  Label::AttrType
  vlabel::Attr(V,Label)
  elabel::Attr(E,Label)
end

@acset_type VELabeledGraph(TheoryVELabeledGraph, index=[:src,:tgt]) <: AbstractGraph
const VELGraph = VELabeledGraph

const OpenVELGraphOb, OpenVELGraph = OpenACSetTypes(VELGraph, :V)

g0 = VELGraph{Symbol}()
add_vertices!(g0, 2, vlabel=[:x,:y])
add_edge!(g0, 1, 2, elabel=:f)
g = OpenVELGraph{Symbol}(g0, FinFunction([1],2), FinFunction([2],2))

h0 = VELGraph{Symbol}()
add_vertices!(h0, 3, vlabel=[:y,:w,:z])
add_edges!(h0, [1,1], [2,3], elabel=[:g,:h])
h = OpenVELGraph{Symbol}(h0, FinFunction([1],3), FinFunction([2,3],3))
lfoot, rfoot = feet(h)
@test nparts(rfoot, :V) == 2
@test subpart(rfoot, :vlabel) == [:w,:z]
@test dom(h) == OpenVELGraphOb{Symbol}(FinSet(1), vlabel=:y)
@test codom(h) == OpenVELGraphOb{Symbol}(FinSet(2), vlabel=[:w,:z])

# Category
#---------

k = compose(g, h)
k0 = apex(k)
@test src(k0) == [1,2,2]
@test tgt(k0) == [2,3,4]
@test subpart(k0, :vlabel) == [:x, :y, :w, :z]
@test subpart(k0, :elabel) == [:f, :g, :h]

# Incompatible attributes.
set_subpart!(h0, 1, :vlabel, :y′)
h = OpenVELGraph{Symbol}(h0, FinFunction([1],3), FinFunction([2,3],3))
@test_throws ErrorException compose(g, h)

# Symmetric monoidal category
#----------------------------

k = otimes(g, h)
@test dom(k) == dom(g)⊗dom(h)
@test codom(k) == codom(g)⊗codom(h)
@test (nv(apex(k)), ne(apex(k))) == (5, 3)

a = OpenVELGraphOb{Symbol}(FinSet(2), vlabel=[:u,:v])
b = OpenVELGraphOb{Symbol}(FinSet(3), vlabel=[:x,:y,:z])
@test dom(braid(a, b)) == a⊗b
@test codom(braid(a, b)) == b⊗a

# Structured cospan of multiply attributed C-sets
#################################################

# Attributed boundary
#--------------------

@present TheoryMultiplyAttributedGraph <: TheoryGraph begin
  Length::AttrType
  Size::AttrType
  Label::AttrType
  vlabel::Attr(V,Label)
  elabel::Attr(E,Label)
  vsize::Attr(V, Size)
  elength::Attr(E, Length)
end

@acset_type MAGraph(TheoryMultiplyAttributedGraph, index=[:src,:tgt]) <: AbstractGraph
const OpenMAGraphOb, OpenMAGraph = OpenACSetTypes(MAGraph, :V)

g0 = @acset MAGraph{Int, Float64, Symbol} begin
  V = 5
  E = 3
  vlabel = [:a, :b, :c, :d, :e]
  elabel = [:a, :b, :c]
  vsize = [0.0, 1.0, 3.0, 3.0, 4.0]
  elength = [1,2,3]
  src = [1, 2, 3]
  tgt = [3, 4, 5]
end
g = OpenMAGraph{Int, Float64, Symbol}(g0, FinFunction([1],5), FinFunction([3,4],5))

g′ = OpenMAGraph{Int, Float64, Symbol}(g0, FinFunction([1],5), FinFunction([3],5), FinFunction([4],5))
@test bundle_legs(g′, [1, (2,3)]) == g

end

max_forks_repo_name: epatters/CompCat


text:
using Base.Test
using DynamicFactorModels
using DataFrames

## Some test data
data = readtable("/home/joi/Documents/Konstanz/Masterarbeit/data/1959-2014_normalized.csv")
data_matrix = reduce(hcat, [convert(Array{Float64}, col) for col in data.columns[2:size(data.columns, 1)]])
#ids = map(string, names(data)[2:end])
#titles = series_titles(ids)  # TODO: does not work at the moment because ICU.jl and with it Requests.jl seems to be broken
T = size(data_matrix, 1) - 4  # we include 4 lags
y = data_matrix[:,1]  # TODO: this is not something we actually want to predict...
data_matrix = data_matrix[:, 2:end]
lag1 = lag_vector(y)
lag2 = lag_vector(lag1)
lag3 = lag_vector(lag2)
lag4 = lag_vector(lag3)
y = y[5:end]
w = hcat(ones(T), array(lag1[5:end]), array(lag2[5:end]), array(lag3[5:end]), array(lag4[5:end]))
x = data_matrix[5:end, :]
T, N = size(x)

# test residuals from factor equation
model = DynamicFactorModel(y,w,x)
sum_squared_residuals = sum(model.factor_residuals.^2)


text:
using Zygote, Test, LinearAlgebra
zygotefun2(a, b) = sigmoid_fast(a + b) + LoopVectorization.relu(a * b) * tanh(a - b)
zygotefun3(a, b, c) = log(tanh_fast(a * b) + 1 + exp(b - c))
@testset "Zygote" begin
  @show @__LINE__

  for T ∈ (Float32, Float64)
    x = randn(T, 217)
    y = randn(T, 217)
    z = randn(T, 217)
    w = randn(T, 217)
    # Test 1-arg and `tanh_fast` vs `tanh`
    gref1 = gradient(x -> dot(w, map(tanh, x)), x)
    glv_1 = gradient(x -> dot(w, vmap(tanh_fast, x)), x)

    @test only(gref1) ≈ only(glv_1)

    # Test 2 arguments
    # Test 3 arguments
    gref2 = gradient(xyz -> dot(w, map(zygotefun2, xyz...)), (x, y))
    g_lv2 = gradient(xyz -> dot(w, vmap(zygotefun2, xyz...)), (x, y))

    @test all(map(≈, only(gref2), only(g_lv2)))

    gref3 = gradient(xyz -> dot(w, map(zygotefun3, xyz...)), (x, y, z))
    g_lv3 = gradient(xyz -> dot(w, vmap(zygotefun3, xyz...)), (x, y, z))

    @test all(map(≈, only(gref3), only(g_lv3)))
  end
end


text:
# Needing this file is a bit unfortunate,
# but seeing https://github.com/GabrielBB/xvfb-action/issues/3
# makes me loose all hope to get the string interpolation right :D
using Pkg; Pkg.test(coverage=true)