text:
"""
Store paths' information built by `learn_paths` or `build_paths`
"""
struct Result_Path_Info_Struct
    ngrams_ind::Array
    num_tolerance::Int64
    support::Float64
end

"""
Store gold paths' information including indices and indices' support and total
support. It can be used to evaluate how low the threshold needs to be set in
order to find most of the correct paths or if set very low, all of the correct paths.
"""
struct Gold_Path_Info_Struct
    ngrams_ind::Array
    ngrams_ind_support::Array
    support::Float64
end

"""
Store threshold and tolerance proportional for each timestep.
"""
struct Threshold_Stat_Struct
    threshold::Float64
    threshold_prop::Array
    tolerance::Float64
    tolerance_prop::Array
    n_rows::Int64
    timestep::Int64
end

"""
A sequence finding algorithm using discrimination learning to predict, for a given
word, which n-grams are best supported for a given position in the sequence of n-grams.
"""
function learn_paths end

"""
The build_paths function constructs paths by only considering those n-grams that are
close to the target. It first takes the predicted c-hat vector and finds the
closest n neighbors in the C matrix. Then it selects all n-grams of these neighbors,
and constructs all valid paths with those n-grams. The path producing the best
correlation with the target semantic vector (through synthesis by analysis) is selected.
"""
function build_paths end

"""
    learn_paths(data_train, data_val, C_train, S_val, F_train, Chat_val, A, i2f, f2i)

A sequence finding algorithm using discrimination learning to predict, for a given
word, which n-grams are best supported for a given position in the sequence of n-grams.

# Obligatory Arguments
- `data::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `C_train::Union{SparseMatrixCSC, Matrix}`: the C matrix for training dataset
- `S_val::Union{SparseMatrixCSC, Matrix}`: the S matrix for validation dataset
- `F_train::Union{SparseMatrixCSC, Matrix}`: the F matrix for training dataset
- `Chat_val::Union{SparseMatrixCSC, Matrix}`: the Chat matrix for validation dataset
- `A::SparseMatrixCSC`: the adjacency matrix
- `i2f::Dict`: the dictionary returning features given indices
- `f2i::Dict`: the dictionary returning indices given features

# Optional Arguments
- `gold_ind::Union{Nothing, Vector}=nothing`: gold paths' indices
- `Shat_val::Union{Nothing, Matrix}=nothing`: the Shat matrix for the validation dataset
- `check_gold_path::Bool=false`: if true, return a list of support values for the gold path; this information is returned as second output value
- `max_t::Int64=15`: maximum timestep
- `max_can::Int64=10`: maximum number of candidates to consider
- `threshold::Float64=0.1`:the value set for the support such that if the support of an n-gram is higher than this value, the n-gram will be taking into consideration
- `is_tolerant::Bool=false`: if true, select a specified number (given by `max_tolerance`) of n-grams whose supports are below threshold but above a second tolerance threshold to be added to the path
- `tolerance::Float64=(-1000.0)`: the value set for the second threshold (in tolerant mode) such that if the support for an n-gram is in between this value and the threshold and the max_tolerance number has not been reached, then allow this n-gram to be added to the path
- `max_tolerance::Int64=4`: maximum number of n-grams allowed in a path
- `grams::Int64=3`: the number n of grams that make up an n-gram
- `tokenized::Bool=false`: if true, the dataset target is tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator token
- `keep_sep::Bool=false`:if true, keep separators in cues
- `target_col::Union{String, :Symbol}=:Words`: the column name for target strings
- `start_end_token::Union{String, Char}="#"`: start and end token in boundary cues
- `issparse::Symbol=:auto`: control of whether output of Mt matrix is a dense matrix or a sparse matrix
- `sparse_ratio::Float64=0.05`: the ratio to decide whether a matrix is sparse
- `if_pca::Bool=false`: turn on to enable pca mode
- `pca_eval_M::Matrix=nothing`: pass original F for pca mode
- `activation::Function=nothing`: the activation function you want to pass
- `ignore_nan::Bool=true`: whether to ignore NaN when compare correlations, otherwise NaN will be selected as the max correlation value
- `check_threshold_stat::Bool=false`: if true, return a threshold and torlerance proportion for each timestep
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# basic usage without tokenization
res = JudiLing.learn_paths(
latin,
latin,
cue_obj.C,
S,
F,
Chat,
A,
cue_obj.i2f,
max_t=max_t,
max_can=10,
grams=3,
threshold=0.1,
tokenized=false,
keep_sep=false,
target_col=:Word,
verbose=true)

# basic usage with tokenization
res = JudiLing.learn_paths(
french,
french,
cue_obj.C,
S,
F,
Chat,
A,
cue_obj.i2f,
max_t=max_t,
max_can=10,
grams=3,
threshold=0.1,
tokenized=true,
sep_token="-",
keep_sep=true,
target_col=:Syllables,
verbose=true)

# basic usage for validation data
res_val = JudiLing.learn_paths(
latin_train,
latin_val,
cue_obj_train.C,
S_val,
F_train,
Chat_val,
A,
cue_obj_train.i2f,
max_t=max_t,
max_can=10,
grams=3,
threshold=0.1,
tokenized=false,
keep_sep=false,
target_col=:Word,
verbose=true)

# turn on tolerance mode
res_val = JudiLing.learn_paths(
...
threshold=0.1,
is_tolerant=true,
tolerance=-0.1,
max_tolerance=4,
...)

# turn on check gold paths mode
res_train, gpi_train = JudiLing.learn_paths(
...
gold_ind=cue_obj_train.gold_ind,
Shat_val=Shat_train,
check_gold_path=true,
...)

res_val, gpi_val = JudiLing.learn_paths(
...
gold_ind=cue_obj_val.gold_ind,
Shat_val=Shat_val,
check_gold_path=true,
...)

# control over sparsity
res_val = JudiLing.learn_paths(
...
issparse=:auto,
sparse_ratio=0.05,
...)

# pca mode
res_learn = JudiLing.learn_paths(
korean,
korean,
Array(Cpcat),
S,
F,
ChatPCA,
A,
cue_obj.i2f,
cue_obj.f2i,
check_gold_path=false,
gold_ind=cue_obj.gold_ind,
Shat_val=Shat,
max_t=max_t,
max_can=10,
grams=3,
threshold=0.1,
tokenized=true,
sep_token="_",
keep_sep=true,
target_col=:Verb_syll,
if_pca=true,
pca_eval_M=Fo,
verbose=true);

```
"""
function learn_paths(
    data_train,
    data_val,
    C_train,
    S_val,
    F_train,
    Chat_val,
    A,
    i2f,
    f2i;
    gold_ind = nothing,
    Shat_val = nothing,
    check_gold_path = false,
    max_t = 15,
    max_can = 10,
    threshold = 0.1,
    is_tolerant = false,
    tolerance = (-1000.0),
    max_tolerance = 3,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    target_col = "Words",
    start_end_token = "#",
    issparse = :auto,
    sparse_ratio = 0.05,
    if_pca = false,
    pca_eval_M = nothing,
    activation = nothing,
    ignore_nan = true,
    check_threshold_stat = false,
    verbose = false
)

    # initialize queues for storing paths
    n_val = size(data_val, 1)
    # working_q = Array{Queue{Array{Int64,1}},1}(undef, n_val)
    working_q = Vector{Queue{Tuple{Vector{Int64},Int64}}}(undef, n_val)
    # res = Array{Array{Array{Int64}},1}(undef, n_val)
    res = Vector{Vector{Tuple{Vector{Int64},Int64}}}(undef, n_val)
    for j = 1:n_val
        res[j] = Tuple{Vector{Int64},Int64}[]
    end

    # store thr proportion for each timestep and each utterance
    thr_stat = zeros(Float64, (n_val, max_t))
    tlr_stat = zeros(Float64, (n_val, max_t))

    # # initialize gold_path_info supports
    if check_gold_path && !isnothing(gold_ind)
        gold_path_info_supports =
            Vector{Vector{Float64}}(undef, size(data_val, 1))
        for i = 1:n_val
            gi = gold_ind[i]
            gold_path_info_supports[i] = Vector{Float64}(undef, length(gi))
        end
    end

    verbose && println("Making fac C")
    fac_C_train = make_transform_fac(C_train)

    C_dim = size(Chat_val, 2)

    for i = 1:max_t
        verbose && println("="^10)
        verbose && println("Timestep $i")
        verbose &&
            i >= 2 &&
            println("average $(mean(length.(working_q))) of paths currently")

        verbose && println("Calculating Yt...")
        Yt_train = make_Yt_matrix(
            i,
            data_train,
            f2i,
            grams = grams,
            target_col = target_col,
            tokenized = tokenized,
            sep_token = sep_token,
            keep_sep = keep_sep,
            start_end_token = start_end_token
        )

        verbose && println("Calculating Mt...")
        Mt_train = make_transform_matrix(
            fac_C_train,
            C_train,
            Yt_train,
            output_format = issparse,
            sparse_ratio = sparse_ratio,
            verbose = verbose,
        )
        # Mt = sparse(Mt)
        # verbose && println("Sparsity of Mt: $(length(Mt.nzval)/Mt.m/Mt.n)")

        verbose && println("Calculating Ythat...")
        # Ythat_val = Matrix(Chat_val)*Matrix(Mt_train)
        Ythat_val = Chat_val * Mt_train
        if is_truly_sparse(Ythat_val, verbose = verbose)
            Ythat_val = sparse(Ythat_val)
        end

        # apply activation to Yt hat
        if !isnothing(activation)
            Ythat_val = activation.(Ythat_val)
        end
        # collect supports for gold path each timestep
        if check_gold_path && !isnothing(gold_ind)
            for j = 1:size(data_val, 1)
                gi = gold_ind[j]
                if i <= length(gi)
                    gold_path_info_supports[j][i] = Ythat_val[j, gi[i]]
                end
            end
        end

        verbose && println("Finding paths...")
        iter = 1:n_val
        if verbose
            pb = Progress(n_val)
        end

        for j in iter
            # collect all n-grams which has greater support than the threshold
            candidates_t = findall(x -> x > threshold, Ythat_val[j, :])
            candidates_t_tlr =
                findall(x -> x > tolerance && x <= threshold, Ythat_val[j, :])

            # calculate threshold stat
            if check_threshold_stat
                prop_candidates_t = length(candidates_t) / C_dim
                prop_candidates_t_tlr = length(candidates_t_tlr) / C_dim
                thr_stat[j, i] = prop_candidates_t
                tlr_stat[j, i] = prop_candidates_t_tlr
            end

            # for timestep 2 and after 2
            if isassigned(working_q, j)
                tmp_working_q = Queue{Tuple{Vector{Int64},Int64}}()
                while !isempty(working_q[j])
                    a = dequeue!(working_q[j]) ## a = [11] A[11,5] == 1 # candidates = [1, 5, 7]

                    for c in candidates_t ## c = 5 # a = [11, 1, 5, 7] # a = [11, 1] [11, 5] [11, 7]
                        # if a n-grams is attachable then append it
                        if isattachable(a[1], c, A)
                            a_copy = deepcopy(a[1])
                            push!(a_copy, c)
                            # if the path is complete then move it to result list
                            if iscomplete(
                                a_copy,
                                i2f,
                                tokenized = tokenized,
                                sep_token = sep_token,
                                start_end_token = start_end_token,
                            )
                                push!(res[j], (a_copy, a[2]))
                            else
                                # otherwise enqueue it to the next timestep
                                enqueue!(tmp_working_q, (a_copy, a[2]))
                            end
                        end
                    end

                    if is_tolerant && a[2] < max_tolerance
                        for c in candidates_t_tlr ## c = 5 # a = [11, 1, 5, 7] # a = [11, 1] [11, 5] [11, 7]
                            # if a n-grams is attachable then append it
                            if isattachable(a[1], c, A)
                                a_copy = deepcopy(a[1])
                                push!(a_copy, c)
                                # if the path is complete then move it to result list
                                if iscomplete(
                                    a_copy,
                                    i2f,
                                    tokenized = tokenized,
                                    sep_token = sep_token,
                                    start_end_token = start_end_token,
                                )
                                    push!(res[j], (a_copy, a[2] + 1))
                                else
                                    # otherwise enqueue it to the next timestep
                                    enqueue!(tmp_working_q, (a_copy, a[2] + 1))
                                end
                            end
                        end
                    end

                end

                # refresh queue for the next timestep
                working_q[j] = tmp_working_q
                # for timestep 1
            else
                working_q[j] = Queue{Tuple{Vector{Int64},Int64}}()
                for c in candidates_t
                    # check whether a n-gram is a start n-gram
                    if isstart(
                        c,
                        i2f,
                        tokenized = tokenized,
                        sep_token = sep_token,
                        start_end_token = start_end_token
                    )
                        a = Int64[]
                        push!(a, c)
                        # check whether this n-gram is both start and complete
                        if iscomplete(
                            a,
                            i2f,
                            tokenized = tokenized,
                            sep_token = sep_token,
                            start_end_token = start_end_token,
                        )
                            push!(res[j], (a, 0))
                        else
                            enqueue!(working_q[j], (a, 0))
                        end
                    end
                end

                if is_tolerant && 0 < max_tolerance
                    for c in candidates_t_tlr
                        # check whether a n-gram is a start n-gram
                        if isstart(
                            c,
                            i2f,
                            tokenized = tokenized,
                            sep_token = sep_token,
                            start_end_token = start_end_token,
                        )
                            a = Int64[]
                            push!(a, c)
                            # check whether this n-gram is both start and complete
                            if iscomplete(
                                a,
                                i2f,
                                tokenized = tokenized,
                                sep_token = sep_token,
                                start_end_token = start_end_token,
                            )
                                push!(res[j], (a, 1))
                            else
                                enqueue!(working_q[j], (a, 1))
                            end
                        end
                    end
                end
            end
            if verbose
                ProgressMeter.next!(pb)
            end
        end
    end

    ts = Threshold_Stat_Struct(threshold, thr_stat, tolerance, tlr_stat, n_val, max_t)

    verbose && println("Evaluating paths...")
    res =
        eval_can(res, S_val, F_train, i2f, max_can, if_pca, pca_eval_M, ignore_nan, verbose)

    # initialize gpi
    if check_gold_path
        if isnothing(gold_ind)
            throw(ErrorException("gold_ind is nothing! Perhaps you forgot to pass gold_ind as an argument."))
        end

        if isnothing(Shat_val)
            throw(ErrorException("Shat_val is nothing! Perhaps you forgot to pass Shat_val as an argument."))
        end

        gpi =
            Vector{Gold_Path_Info_Struct}(undef, size(data_val, 1))

        # calculate all shat correlation with S
        Scors = [cor(Shat_val[i, :], S_val[i, :]) for i = 1:n_val]

        for i = 1:size(data_val, 1)
            gpi[i] = Gold_Path_Info_Struct(
                gold_ind[i],
                gold_path_info_supports[i],
                Scors[i],
            )
        end
    end

    if check_gold_path
        if check_threshold_stat
            return res, gpi, ts
        else
            return res, gpi
        end
    else
        if check_threshold_stat
            return res, ts
        else
            return res
        end
    end
end

"""
    learn_paths(data, cue_obj, S_val, F_train, Chat_val)

A high-level wrapper function for learn_paths with much less control. It aims
for users who is very new to JudiLing and learn_paths function.

# Obligatory Arguments
- `data::DataFrame`: the training dataset
- `cue_obj::Cue_Matrix_Struct`: the C matrix object containing all information with C
- `S_val::Union{SparseMatrixCSC, Matrix}`: the S matrix for validation dataset
- `F_train::Union{SparseMatrixCSC, Matrix}`: the F matrix for training dataset
- `Chat_val::Union{SparseMatrixCSC, Matrix}`: the Chat matrix for validation dataset

# Optional Arguments
- `Shat_val::Union{Nothing, Matrix}=nothing`: the Shat matrix for the validation dataset
- `check_gold_path::Bool=false`: if true, return a list of support values for the gold path; this information is returned as second output value
- `threshold::Float64=0.1`:the value set for the support such that if the support of an n-gram is higher than this value, the n-gram will be taking into consideration
- `is_tolerant::Bool=false`: if true, select a specified number (given by `max_tolerance`) of n-grams whose supports are below threshold but above a second tolerance threshold to be added to the path
- `tolerance::Float64=(-1000.0)`: the value set for the second threshold (in tolerant mode) such that if the support for an n-gram is in between this value and the threshold and the max_tolerance number has not been reached, then allow this n-gram to be added to the path
- `max_tolerance::Int64=4`: maximum number of n-grams allowed in a path
- `activation::Function=nothing`: the activation function you want to pass
- `ignore_nan::Bool=true`: whether to ignore NaN when compare correlations, otherwise NaN will be selected as the max correlation value
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
res = learn_paths(latin, cue_obj, S, F, Chat)
```
"""
function learn_paths(
    data,
    cue_obj,
    S_val,
    F_train,
    Chat_val;
    Shat_val = nothing,
    check_gold_path = false,
    threshold = 0.1,
    is_tolerant = false,
    tolerance = (-1000.0),
    max_tolerance = 3,
    activation = nothing,
    ignore_nan = true,
    verbose = true)

    max_t = JudiLing.cal_max_timestep(data, cue_obj.target_col,
        tokenized=cue_obj.tokenized, sep_token=cue_obj.sep_token)

    learn_paths(
        data,
        data,
        cue_obj.C,
        S_val,
        F_train,
        Chat_val,
        cue_obj.A,
        cue_obj.i2f,
        cue_obj.f2i;
        gold_ind = cue_obj.gold_ind,
        Shat_val = Shat_val,
        check_gold_path = check_gold_path,
        max_t = max_t,
        max_can = 10,
        threshold = threshold,
        is_tolerant = is_tolerant,
        tolerance = tolerance,
        max_tolerance = max_tolerance,
        grams = cue_obj.grams,
        tokenized = cue_obj.tokenized,
        sep_token = cue_obj.sep_token,
        keep_sep = cue_obj.keep_sep,
        target_col = cue_obj.target_col,
        activation = activation,
        ignore_nan = ignore_nan,
        verbose = verbose,
    )
end

"""
    build_paths(data_val, C_train, S_val, F_train, Chat_val, A, i2f, C_train_ind)

The build_paths function constructs paths by only considering those n-grams that are
close to the target. It first takes the predicted c-hat vector and finds the
closest n neighbors in the C matrix. Then it selects all n-grams of these neighbors,
and constructs all valid paths with those n-grams. The path producing the best
correlation with the target semantic vector (through synthesis by analysis) is selected.

# Obligatory Arguments
- `data::DataFrame`: the training dataset
- `data_val::DataFrame`: the validation dataset
- `C_train::SparseMatrixCSC`: the C matrix for the training dataset
- `S_val::Union{SparseMatrixCSC, Matrix}`: the S matrix for the validation dataset
- `F_train::Union{SparseMatrixCSC, Matrix}`: the F matrix for the training dataset
- `Chat_val::Matrix`: the Chat matrix for the validation dataset
- `A::SparseMatrixCSC`: the adjacency matrix
- `i2f::Dict`: the dictionary returning features given indices
- `C_train_ind::Array`: the gold paths' indices for the training dataset

# Optional Arguments
- `rC::Union{Nothing, Matrix}=nothing`: correlation Matrix of C and Chat, specify to save computing time
- `max_t::Int64=15`: maximum number of timesteps
- `max_can::Int64=10`: maximum number of candidates to consider
- `n_neighbors::Int64=10`: the top n form neighbors to be considered
- `grams::Int64=3`: the number n of grams that make up n-grams
- `tokenized::Bool=false`: if true, the dataset target is tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator
- `target_col::Union{String, :Symbol}=:Words`: the column name for target strings
- `if_pca::Bool=false`: turn on to enable pca mode
- `pca_eval_M::Matrix=nothing`: pass original F for pca mode
- `verbose::Bool=false`: if true, more information will be printed

# Examples
```julia
# training dataset
JudiLing.build_paths(
    latin_train,
    cue_obj_train.C,
    S_train,
    F_train,
    Chat_train,
    A,
    cue_obj_train.i2f,
    cue_obj_train.gold_ind,
    max_t=max_t,
    n_neighbors=10,
    verbose=false
    )

# validation dataset
JudiLing.build_paths(
    latin_val,
    cue_obj_train.C,
    S_val,
    F_train,
    Chat_val,
    A,
    cue_obj_train.i2f,
    cue_obj_train.gold_ind,
    max_t=max_t,
    n_neighbors=10,
    verbose=false
    )

# pca mode
res_build = JudiLing.build_paths(
    korean,
    Array(Cpcat),
    S,
    F,
    ChatPCA,
    A,
    cue_obj.i2f,
    cue_obj.gold_ind,
    max_t=max_t,
    if_pca=true,
    pca_eval_M=Fo,
    n_neighbors=3,
    verbose=true
    )
```
"""
function build_paths(
    data_val,
    C_train,
    S_val,
    F_train,
    Chat_val,
    A,
    i2f,
    C_train_ind;
    rC = nothing,
    max_t = 15,
    max_can = 10,
    n_neighbors = 10,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    target_col = :Words,
    start_end_token = "#",
    if_pca = false,
    pca_eval_M = nothing,
    ignore_nan = true,
    verbose = false,
)
    # initialize queues for storing paths
    n_val = size(data_val, 1)
    # working_q = Array{Queue{Array{Int64,1}},1}(undef, n_val)
    # res = Array{Array{Array{Int64}},1}(undef, n_val)
    res = Vector{Vector{Tuple{Vector{Int64},Int64}}}(undef, n_val)

    for j = 1:n_val
        res[j] = Tuple{Vector{Int64},Int64}[]
    end

    verbose && println("Finding all top features..")
    # findall features indices for all utterances
    isnothing(rC) && begin
        rC = cor(Chat_val, Matrix(C_train), dims = 2)
    end
    top_indices = find_top_feature_indices(
        # C_train,
        # Chat_val,
        rC,
        C_train_ind,
        n_neighbors = n_neighbors,
        verbose = verbose,
    )

    # verbose && println("="^10)
    # verbose && println("Timestep $i")

    verbose && println("Finding paths...")
    iter = 1:n_val
    if verbose
        pb = Progress(n_val)
    end
    for j in iter
        candidates_t = top_indices[j]

        # timestep 1
        working_q = Queue{Array{Int64,1}}()
        for c in candidates_t
            # check whether a n-gram is a start n-gram
            if isstart(c,
                i2f,
                tokenized = tokenized,
                sep_token = sep_token,
                start_end_token = start_end_token,
            )

                a = Int64[]
                push!(a, c)
                # check whether this n-gram is both start and complete
                if iscomplete(
                    a,
                    i2f,
                    tokenized = tokenized,
                    sep_token = sep_token,
                    start_end_token = start_end_token,
                )
                    push!(res[j], (a, 0))
                else
                    enqueue!(working_q, a)
                end
            end
        end

        for i = 2:max_t
            tmp_working_q = Queue{Array{Int64,1}}()
            while !isempty(working_q)
                a = dequeue!(working_q) ## a = [11] A[11,5] == 1 # candidates = [1, 5, 7]
                for c in candidates_t ## c = 5 # a = [11, 1, 5, 7] # a = [11, 1] [11, 5] [11, 7]
                    # if a n-grams is attachable then append it
                    if isattachable(a, c, A)
                        a_copy = deepcopy(a)
                        push!(a_copy, c)
                        # if the path is complete then move it to result list
                        if iscomplete(
                            a_copy,
                            i2f,
                            tokenized = tokenized,
                            sep_token = sep_token,
                            start_end_token = start_end_token,
                        )
                            push!(res[j], (a_copy, 0))
                        else
                            # otherwise enqueue it to the next timestep
                            enqueue!(tmp_working_q, a_copy)
                        end
                    end
                end
            end

            # refresh queue for the next timestep
            working_q = tmp_working_q
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    verbose && println("Evaluating paths...")
    eval_can(res, S_val, F_train, i2f, max_can, if_pca, pca_eval_M, ignore_nan, verbose)
end

"""
    learn_paths_rpi(data_train, data_val, C_train, S_val, F_train, Chat_val, A, i2f, f2i)

Calculate learn_paths with results indices supports as well.
"""
function learn_paths_rpi(
    data_train,
    data_val,
    C_train,
    S_val,
    F_train,
    Chat_val,
    A,
    i2f,
    f2i;
    gold_ind = nothing,
    Shat_val = nothing,
    check_gold_path = false,
    max_t = 15,
    max_can = 10,
    threshold = 0.1,
    is_tolerant = false,
    tolerance = (-1000.0),
    max_tolerance = 3,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    target_col = "Words",
    start_end_token = "#",
    issparse = :auto,
    sparse_ratio = 0.05,
    if_pca = false,
    pca_eval_M = nothing,
    activation = nothing,
    ignore_nan = true,
    check_threshold_stat = false,
    verbose = false
    )

    res = learn_paths(
        data_train,
        data_val,
        C_train,
        S_val,
        F_train,
        Chat_val,
        A,
        i2f,
        f2i,
        gold_ind = gold_ind,
        Shat_val = Shat_val,
        check_gold_path = check_gold_path,
        max_t = max_t,
        max_can = max_can,
        threshold = threshold,
        is_tolerant = is_tolerant,
        tolerance = tolerance,
        max_tolerance = max_tolerance,
        grams = grams,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        target_col = target_col,
        start_end_token = start_end_token,
        issparse = issparse,
        sparse_ratio = sparse_ratio,
        if_pca = if_pca,
        pca_eval_M = pca_eval_M,
        activation = activation,
        ignore_nan = ignore_nan,
        check_threshold_stat = check_threshold_stat,
        verbose = verbose
    )

    if check_gold_path
        res = res[1]
        gpi = res[2]
    end

    n = size(res)
    ngrams_ind = make_ngrams_ind(res, n)
    Shat = zeros(Float64, size(S_val))

    for i in 1:n[1]
        ci = ngrams_ind[i]
        Shat[i,:] = sum(F_train[ci, :], dims = 1)
    end

    tmp, rpi = learn_paths(
        data_train,
        data_val,
        C_train,
        S_val,
        F_train,
        Chat_val,
        A,
        i2f,
        f2i,
        gold_ind = ngrams_ind,
        Shat_val = Shat,
        check_gold_path = true,
        max_t = max_t,
        max_can = 1,
        threshold = 1,
        is_tolerant = false,
        tolerance = 1,
        max_tolerance = 1,
        grams = grams,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        target_col = target_col,
        start_end_token = start_end_token,
        issparse = issparse,
        sparse_ratio = sparse_ratio,
        if_pca = if_pca,
        pca_eval_M = pca_eval_M,
        activation = activation,
        ignore_nan = ignore_nan,
        check_threshold_stat = check_threshold_stat,
        verbose = false
    )

    if check_gold_path
        return res, gpi, rpi
    else
        return res, rpi
    end

end

"""
    eval_can(candidates, S, F, i2f, max_can, if_pca, pca_eval_M)

Calculate for each candidate path the correlation between predicted semantic
vector and the gold standard semantic vector, and select as target for production
the path with the highest correlation.
"""
function eval_can(
    candidates,
    S,
    F,
    i2f,
    max_can,
    if_pca,
    pca_eval_M,
    ignore_nan = true,
    verbose = false,
)

    verbose &&
        println("average $(mean(length.(candidates))) of paths to evaluate")

    res_l = Array{Array{Result_Path_Info_Struct,1},1}(undef, size(S, 1))
    iter = 1:size(S, 1)
    if verbose
        pb = Progress(size(S, 1))
    end

    if if_pca
        F = pca_eval_M
    end

    Threads.@threads for i in iter
        tid = Threads.threadid()
        res = Result_Path_Info_Struct[]
        if size(candidates[i], 1) > 0
            for (ci, n) in candidates[i] # ci = [1,3,4]
                Shat = sum(F[ci, :], dims = 1)
                Scor = cor(Shat[1, :], S[i, :])
                push!(res, Result_Path_Info_Struct(ci, n, Scor))
            end
        end

        if ignore_nan
            res = filter(x -> !isnan(x.support), res)
        end

        # we collect only top x candidates from the top
        res_l[i] = collect(Iterators.take(
            sort!(res, by = x -> x.support, rev = true),
            max_can,
        ))
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    res_l
end

"""
    find_top_feature_indices(rC, C_train_ind)

Find all indices for the n-grams of the top n closest neighbors of
a given target.
"""
function find_top_feature_indices(
    rC,
    C_train_ind;
    n_neighbors = 10,
    verbose = false,
)

    # collect num of val data
    n_val = size(rC, 1)

    # calculate correlation matrix
    # rC = cor(Chat_val, Matrix(C_train), dims=2)
    # display(rC)

    # initialize features list for all candidates
    features_all = Vector{Vector{Int64}}(undef, n_val)

    # create iter for progress bar
    verbose && println("finding all n_neighbors features...")
    iter = 1:n_val
    if verbose
        pb = Progress(n_val)
    end

    # find all features of n_neighbors
    for i in iter
        features = [
            C_train_ind[j]
            for j in sortperm(rC[i, :], rev = true)[1:n_neighbors]
        ]
        features_all[i] = unique(collect(Iterators.flatten(features)))
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    features_all
end

"""
    make_ngrams_ind(res, n)

Construct ngrams indices.
"""
function make_ngrams_ind(
    res,
    n
)
    ngrams_ind = Array{Array{Int, 1}, 1}(undef, n)
    for i = 1:n[1]
        if isempty(res[i])
            ngrams_ind[i] = []
        else
            ngrams_ind[i] = res[i][1].ngrams_ind
        end
    end
    ngrams_ind
end


text:
using SemiseparableMatrices, BandedMatrices, LinearAlgebra, Test
import SemiseparableMatrices: LowRankMatrix

@testset "LowRankMatrix" begin
    @test LowRankMatrix(Zeros(5,5), 2) == Zeros(5,5)
end

@testset "SemiseparableMatrix" begin
    n = 10
    v = rand(n, 2); w = rand(n, 2); y = rand(n, 3); z = rand(n, 3);
    U = ApplyMatrix(*, v, w'); L = ApplyMatrix(*, y, z')
    @test SemiseparableMatrix(L, U, 1, 2) isa SemiseparableMatrix{Float64}
    @test Matrix(SemiseparableMatrix(L, U, 1, 2)) ≈ tril(L, -1) + triu(U, 2)
    @test_throws ArgumentError SemiseparableMatrix(L, U, 11, 2)
    @test_throws ArgumentError SemiseparableMatrix(L, U, 1, 11)

    y = rand(n+1, 3); z = rand(n+1, 3)
    U = LowRankMatrix(v, w'); L = LowRankMatrix(y, z')
    @test_throws DimensionMismatch SemiseparableMatrix(L, U, 1, 2)
end

ext: jl


text:
# module for tap(ing) the simulation to achieve playaback at a later time
module tap

dir = @__DIR__
push!(LOAD_PATH, dir * "/../sim")
push!(LOAD_PATH, dir * "/../vis")
using sim
using vis

include(dir * "/../update_info.jl")

mutable struct TapeFrame
  X::Array{Array{Float64, 1}, 1}
  DIAG::Array{Tuple, 1}
  CV::Array{Array{Float64, 1}, 1}

  TapeFrame() = new(Array{Float64, 1}[], Tuple[], Array{Float64, 1}[])
end

mutable struct Tape
  world::World
  frames::Array{TapeFrame, 1}
  info_objs::Array{Union{vis.InfoBox, Nothing}, 1}
  frame_last::Int

  function Tape(world)
    info_objs = Array{Union{vis.InfoBox, Nothing}, 1}(undef,
                                                      length(world.agents))
    for i in 1:length(info_objs)
      info_objs[i] = nothing
    end
    new(world, TapeFrame[], info_objs, 0)
  end
end

function add_frame!(tape::Tape)
  push!(tape.frames, TapeFrame())
  tape.frame_last += 1
end

function add_to_frame!(tape::Tape, x::Array{Float64, 1}, diag::Tuple,
                       cv::Array{Float64, 1}=[NaN, NaN])
  push!(tape.frames[tape.frame_last].X, x)
  push!(tape.frames[tape.frame_last].DIAG, diag)
  push!(tape.frames[tape.frame_last].CV, cv)
end

function visualize_frame(tape::Tape, f::Int)
  @assert f > 0 && f <= length(tape.frames)

  to_vis = []
  if tape.world.road.road != nothing
    push!(to_vis, tape.world.road.road)
  end
  i = 1
  for agent in tape.world.agents
    agent_x = tape.frames[f].X[i]
    diag = tape.frames[f].DIAG[i]
    cv = tape.frames[f].CV[i]

    if tape.info_objs[i] != nothing
      update_info_from_tape(tape.info_objs[i], agent_x, agent.id, diag, cv)
      push!(to_vis, tape.info_objs[i])
    end
    (x, y, sx, sy, dx, u) = diag
    ## get additional information
    dp = dx[3]
    th = atan(sy, sx)
    th_car = atan(dp, agent_x[2]) + atan(sy, sx)
    agent_is_braking = dx[1] * u[1] < 0
    if agent.car != nothing
      ## update the car
      vis.car_lights!(agent.car, agent_is_braking)
      agent.car.T = (vis.scale_mat(tape.world.vis_scaling,
                                   tape.world.vis_scaling) *
                     vis.translate_mat(x, y) *
                     vis.rotate_mat(th_car - pi / 2))
      push!(to_vis, agent.car)
    end
    i += 1
  end

  return to_vis
end

end


text:
@testset "Visualization Unobserved Components" begin
    air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
    log_air_passengers = log.(air_passengers.passengers)
    model = UnobservedComponents(log_air_passengers; trend = "local linear trend", seasonal = "stochastic 12")
    fit!(model)
    kf = kalman_filter(model)
    ks = kalman_smoother(model)

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, kf)
    @test length(rec) == 4
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, ks)
    @test length(rec) == 4

    finland_fatalities = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
    log_finland_fatalities = log.(finland_fatalities.ff)
    model = UnobservedComponents(log_finland_fatalities; trend = "local linear trend")
    fit!(model)

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, kf)
    @test length(rec) == 3
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, ks)
    @test length(rec) == 3
end


text:
module Agents

using Requires
using Scratch
using Distributed
using DataStructures
using LightGraphs
using DataFrames
using Random
import ProgressMeter

import Base.length
import LinearAlgebra

# Core structures of Agents.jl
include("core/agents.jl")
include("core/model.jl")
include("core/space_interaction_API.jl")

# Existing spaces
include("spaces/nothing.jl")
include("spaces/graph.jl")
include("spaces/grid.jl")
include("spaces/discrete.jl")
include("spaces/continuous.jl")
include("spaces/openstreetmap.jl")
include("spaces/utilities.jl")

# Stepping and data collection functionality
include("simulations/step.jl")
include("simulations/collect.jl")
include("simulations/paramscan.jl")
include("simulations/sample.jl")
include("simulations/ensemblerun.jl")

# Other features that exist in submodules
include("submodules/pathfinding/all_pathfinders.jl")
include("submodules/schedulers.jl")
include("submodules/io/AgentsIO.jl")

function __init__()
    # Plot recipes
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("visualization/plot-recipes.jl")
    end
    # Workaround for Documenter.jl, so we don't need to include
    # heavy dependencies to build documentation
    @require Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" begin
        include("visualization/plot-recipes.jl")
    end
    # Update message:
    display_update = true
    version_number = "4.4"
    update_name = "update_v$(version_number)"

    if display_update
        # Get scratch space for this package
        versions_dir = @get_scratch!("versions")
        if !isfile(joinpath(versions_dir, update_name))
            printstyled(
                stdout,
                """
                \nUpdate message: Agents v$(version_number)
                Please see the changelog online. Some key features:

                * Agent data can be loaded from and saved to CSV files using `populate_from_csv!` and `dump_to_csv`
                * Support for saving and loading entire models using `save_checkpoint` and `load_checkpoint`
                """;
                color = :light_magenta,
            )
            touch(joinpath(versions_dir, update_name))
        end
    end
end

# Deprecations, that will be removed in future versions
include("deprecated.jl")

# Predefined models
include("models/Models.jl")
export Models

end # module


text:
function interreduce(ctx::SigPolynomialΓ{I, M},
                     G::Basis{I, M},
                     H::Syz{I, M};
                     verbose = false) where {I, M}

    unit = Base.one(ctx.po.mo)
    to_reduce = Set([(unit, (i, g)) for (i, Gi) in G for (g, lm) in Gi])

    symbolic_pp_timed = @timed symbolic_pp!(ctx, to_reduce, G, H, are_pairs = false, interreduction_step = true)
    mons = symbolic_pp_timed.value
    interred_mat = f5matrix(ctx, mons, collect(to_reduce), interreduction_step = true)
    mat_size = (length(rows(interred_mat)), length(interred_mat.tbl))

    if verbose
        println("symbolic pp for interreduction matrix took $(symbolic_pp_timed.time) seconds.")
        println("size of interreduction matrix: $(mat_size)")
        println("reducing...")
    end

    red = @timed reduction!(interred_mat, interreduction_step = true)
    red_time = red.time
    arit_operations_interreduction, _ = red.value

    if verbose
        println("reduction of interreduction matrix took $(red_time) seconds.")
    end

    G_new = new_basis(ctx, G)
    for (sig, row) in interred_mat.sigs_rows
        (m, (pos_key, t)) = sig
        m != unit && continue
        isempty(row) && continue
        p = unindexpolynomial(interred_mat.tbl, row)
        ctx(sig[2], p)
        push!(G_new[pos_key], (t, leadingmonomial(p)))
    end

    G_new, arit_operations_interreduction, mat_size
end

struct InterredOrder{I, M, J, T, SΓ <: SigPolynomialΓ}<:Base.Order.Ordering
    ctx::SΓ
    sigs_lms::Dict{MonSigPair{I, M}, J}
end
function interredorder(ctx::SigPolynomialΓ{I, M, T},
                       sigs_lms::Dict{MonSigPair{I, M}, J}) where {I, M, J, T}
    InterredOrder{I, M, J, T, typeof(ctx)}(ctx, sigs_lms)
end

function Base.Order.lt(ord::InterredOrder{I, M},
                       sig1::MonSigPair{I, M},
                       sig2::MonSigPair{I, M}) where {I, M}

    ctx = ord.ctx
    unit = Base.one(ctx.po.mo)
    sig1[1] == unit && sig2[1] != unit && return false
    sig1[1] != unit && sig2[1] == unit && return true
    ord.sigs_lms[sig1] > ord.sigs_lms[sig2]
end


text:
include("core.jl")
include("utils.jl")
abstract type AbstractCubicModel <: AbstractHelmholtzModel end
abstract type AbstractAlphaT end
struct PengRobinson{Mixture,AlphaT} <: AbstractCubicModel
    criticalTemperature::Array{Float64,1}
    criticalPressure::Array{Float64,1}
    criticalVolume::Array{Float64,1}
    acentricFactor::Array{Float64,1}
    delta1::Float64
    delta2::Float64
    constanta::Float64
    constantb::Float64
    alphat::AlphaT
    mixture::Mixture
end

function PengRobinson(data,alphat::abstractAlphaT,mixture::AbstractMixture)
    Tc=disallowmissing(getproperty(data,:criticalTemperature))
    Pc=disallowmissing(getproperty(data,:criticalPressure))
    Vc=disallowmissing(getproperty(data,:criticalVolume))
    omega=disallowmissing(getproperty(data,:acentricFactor))
    return Cubic(Tc,Pc,Vc,omega,
    1+sqrt(2),1-sqrt(2),0.45724,0.07780,
    alphat,mixture)
end

struct Soave <: abstractAlphaT end


text:
using Landscapes
using Documenter

DocMeta.setdocmeta!(Landscapes, :DocTestSetup, :(using Landscapes); recursive=true)

makedocs(;
    modules=[Landscapes],
    authors="Rafael Schouten <rafaelschouten@gmail.com>",
    repo="https://github.com/rafaqz/Landscapes.jl/blob/{commit}{path}#{line}",
    sitename="Landscapes.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rafaqz.github.io/Landscapes.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rafaqz/Landscapes.jl",
    devbranch="main",
)


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,DifferentialEquations,JLD2,DataFrames,StatsPlots,FileIO,MAT,RecursiveArrayTools
import KenyaCoV
using LinearAlgebra:eigen
using Statistics: median, quantile

gr()#Plotting frontend

"""
Load age structured data, define initial state and declare the KenyaCoV problem
"""

u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
@load "data/agemixingmatrix_china.jld2" M_China
@load "data/reporting_rate_for_MERS_like_scenario.jld2" reporting_rate

u0[KenyaCoV.ind_nairobi_as,5,4] = 5#Five initial infecteds in Nairobi in the 20-24 age group
prob = KenyaCoV.create_KenyaCoV_non_neg_prob(u0,(0.,365.),P)
"""
Declare the treatment/isolation rates considered, and a callback that limits intervention
"""
treatment_rates = [(0.,1),(0.,0.75),(0.,0.5),(1/7.,1.),(1/7,0.75),(1/7,0.5),(1/3.5,1.),(1/3.5,0.75),(1/3.5,0.5)]

function isolation_limit(u,t,integrator) # Isolation can only continue until 10000 detected cases
  integrator.p.isolating_detecteds && sum(u[:,:,8]) > 1e4
end
function affect_isolation_limit!(integrator)
  integrator.p.τ = 0
  integrator.p.isolating_detecteds = false
end
cb_iso_limit = DiscreteCallback(isolation_limit,affect_isolation_limit!)

"""
SCENARIO A1:
* Age-specific susceptibilties calculated from Chinese case data
* Asymptomatics are not infectious
* 20% of infections are symptomatics
"""
P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.
P.δ = 0.2
P.γ = 1/2.5
P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = rand(KenyaCoV.d_R₀)*P.γ/(P.δ + P.ϵ*(1-P.δ))

results_A1 = KenyaCoV.run_scenario(P,prob,10,treatment_rates,cb_iso_limit)

@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_A1.jld2") results_A1

"""
SCENARIO B1:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 50% as infectious as symptomatics
* 20% of infections are symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
sus_matrix = repeat(P.χ,1,KenyaCoV.n_a)
eigs, = eigen(sus_matrix.*M_China)
R₀_scale = Real(eigs[end])
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.5
P.δ = 0.2
P.γ = 1/2.5
P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = rand(KenyaCoV.d_R₀)*P.γ/(P.δ + P.ϵ*(1-P.δ))
results_B1 = KenyaCoV.run_scenario(P,prob,1000,treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_B1.jld2") results_B1
results_B1 = 0

"""
SCENARIO A2:
* No age-dependent susceptibility
* Asymptomatics are not infectious
* 80% of infections are symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")

P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.
P.δ = 0.8
P.γ = 1/2.5
P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = rand(KenyaCoV.d_R₀)*P.γ/(P.δ + P.ϵ*(1-P.δ))

results_A2 = KenyaCoV.run_scenario(P,prob,1000,treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_A2.jld2") results_A2

"""
SCENARIO B2:
* No age-specific susceptibilties --- its disease difference
* Asymptomatics are 50% as infectious as symptomatics
* 80% of infections are symptomatics
"""
u0,P,P_dest = KenyaCoV.model_ingredients_from_data("data/data_for_age_structuredmodel.jld2",
                                            "data/flight_numbers.csv",
                                            "data/projected_global_prevelance.csv")
#Redistribute susceptibility JUST to rescale infectiousness so we get the correct R₀/r
P.χ = ones(KenyaCoV.n_a)
sus_matrix = repeat(P.χ,1,KenyaCoV.n_a)
eigs, = eigen((P.β/P.γ)*sus_matrix.*P.M)
R₀_scale = Real(eigs[end])
P.χ = ones(KenyaCoV.n_a)/R₀_scale

P.dt = 0.25;
P.ext_inf_rate = 0.;
P.ϵ = 0.5
P.δ = 0.8
P.γ = 1/2.5
P.σ = 1/rand(KenyaCoV.d_incubation)
P.β = rand(KenyaCoV.d_R₀)*P.γ/(P.δ + P.ϵ*(1-P.δ))
results_B2 = KenyaCoV.run_scenario(P,prob,1000,treatment_rates,cb_iso_limit)
@save joinpath(homedir(),"Github/KenyaCoVOutputs/results_B2.jld2") results_B2
results_B2 = 0


text:
module Day10

export get_inputs, get_solution1, get_solution2

## Input getting
function read_input(f_name)
    sorted = sort(parse.(Int, readlines(joinpath(@__DIR__, f_name))))
    return join(string.(diff([0; sorted; sorted[end]+3])))
end

function get_inputs()
    test_input1 = test_input2 = read_input.(("test_input11.txt", "test_input12.txt"))
    test_output1 = (5*7, 22*10)
    test_output2 = (8, 19208)
    data = read_input("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
# Part 1
get_solution1(data::Tuple) = get_solution1.(data)
get_solution1(data) = count(==('1'), data) * count(==('3'), data)

# Part 2
# Geez. This took me so long to figure out.
# Since diffs are either 1 or 3, we can count the length of each run of 1-diffs by splitting by the number 3
get_solution2(data::Tuple) = get_solution2.(data)
get_solution2(data) = prod(binomial.(length.(split(data, '3', keepempty=false)), 2) .+ 1)

end


text:
include("../src/SpinMonteCarlo.jl")

using SpinMonteCarlo

const ups = [0,1,2]
const Qs = [2,3,4]
const Ls = [8,16,24]
const MCS = 8192
const Therm = MCS >> 3

params = Dict{String, Any}[]
for update in ups
    for Q in Qs
        const Tc = 1.0/log1p(sqrt(Q))
        const Ts = Tc*linspace(0.85, 1.15, 11)
        for L in Ls
            for T in Ts
                push!(params,
                  Parameter("Model"=>Potts, "Lattice"=>square_lattice,
                            "Q"=>Q, "L"=>L, "T"=>T, "J"=>1.0,
                            "MCS"=>MCS, "Thermalization"=>Therm,
                            "Update Method"=> (update==0 ? local_update! :
                                               update==1 ? SW_update! : Wolff_update!),
                            "update"=>update,
                            "Verbose"=>true,
                           ))
            end
        end
    end
end

obs = map(runMC, params)

const pnames = ["update", "Q", "L", "T"]
const onames = ["Magnetization",
                "|Magnetization|",
                "Magnetization^2",
                "Magnetization^4",
                "Binder Ratio",
                "Susceptibility",
                "Connected Susceptibility",
                "Energy",
                "Energy^2",
                "Specific Heat",
                "MCS per Second",
                "Time per MCS",
               ]

const io = open("res-potts.dat", "w")
i=1
for pname in pnames
    println(io, "# \$$i : $pname")
    i+=1
end
for oname in onames
    println(io, "# \$$i, $(i+1): $oname")
    i+=2
end

for (p,o) in zip(params, obs)
    for pname in pnames
        print(io, p[pname], " ")
    end
    for oname in onames
        @printf(io, "%.15f %.15f ", mean(o[oname]), stderror(o[oname]))
    end
    println(io)
end


text:
function test_physicalmodeequation(scenario)
    @unpack ea, tx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)

    me = PhysicalModeEquation(ea, tx.frequency, waveguide)
    @test me isa LMP.ModeEquation
    @test me isa PhysicalModeEquation

    me2 = PhysicalModeEquation(tx.frequency, waveguide)
    @test me2 isa PhysicalModeEquation
    @test iszero(me2.ea.θ)

    @test LMP.setea(ea, me2) == me
    @test LMP.setea(ea.θ, me2) == me
end

function test_wmatrix(scenario)
    @unpack ea, tx, bfield, species = scenario

    C = ea.cosθ
    L = @SMatrix [C 0 -C 0;
                  0 -1 0 -1;
                  0 -C 0 C;
                  1 0 1 0]

    M = LMP.susceptibility(80e3, tx.frequency, bfield, species)
    T = LMP.tmatrix(ea, M)
    W = LMP.wmatrix(ea, T)

    @test [W[1] W[3]; W[2] W[4]] ≈ 2*(L\T)*L
end

function test_wmatrix_deriv(scenario)
    @unpack tx, bfield, species = scenario

    M = LMP.susceptibility(80e3, tx.frequency, bfield, species)

    for i = 1:4
        for j = 1:4
            Wfcn(θ) = (ea = EigenAngle(θ); T = LMP.tmatrix(ea, M); LMP.wmatrix(ea, T)[i][j])
            dWref = FiniteDiff.finite_difference_derivative(Wfcn, θs, Val{:central})
            dW(θ) = (ea = EigenAngle(θ); T = LMP.tmatrix(ea, M);
                     dT = LMP.dtmatrix(ea, M); LMP.dwmatrix(ea, T, dT)[i][j])

            @test maxabsdiff(dW.(θs), dWref) < 1e-3
        end
    end
end

function test_dRdz(scenario)
    @unpack ea, tx, bfield, species, ground = scenario

    params = LMPParams()
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)

    Mtop = LMP.susceptibility(params.topheight, me; params=params)
    Rtop = LMP.bookerreflection(ea, Mtop)

    # sharply bounded R from bookerreflection satisfies dR/dz = 0
    @test isapprox(LMP.dRdz(Rtop, me, params.topheight), zeros(2, 2); atol=1e-15)

    # Compare default and optional argument
    M = LMP.susceptibility(72e3, me; params=params)
    R = LMP.bookerreflection(ea, M)
    dR1 = @inferred LMP.dRdz(R, me, 72e3)

    Mfcn = z -> LMP.susceptibility(z, me; params=params)
    dR2 = @inferred LMP.dRdz(R, me, 72e3, Mfcn)

    Mfcn2 = LMP.susceptibilityspline(me; params=params)
    dR3 = @inferred LMP.dRdz(R, me, 72e3, Mfcn2)

    @test dR1 == dR2  # identical functions Mfcn
    @test isapprox(dR1, dR3, atol=1e-12)  # interpolating spline
end

function test_dRdθdz(scenario)
    @unpack ea, tx, bfield, species, ground = scenario

    params = LMPParams()
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)

    Mtop = LMP.susceptibility(params.topheight, me)
    z = params.topheight - 500

    for i = 1:4
        dRfcn(θ) = (ea = EigenAngle(θ); Rtop = LMP.bookerreflection(ea, Mtop);
            me = LMP.setea(ea, me); LMP.dRdz(Rtop, me, z)[i])
        dRdθref = FiniteDiff.finite_difference_derivative(dRfcn, θs, Val{:central})
        dRdθtmp(θ) = (ea = EigenAngle(θ); me = LMP.setea(ea, me);
            (Rtop, dRdθtop) = LMP.bookerreflection(ea, Mtop, LMP.Dθ());
            RdRdθtop = vcat(Rtop, dRdθtop);
            LMP.dRdθdz(RdRdθtop, (me, params), z))
        dR(θ) = dRdθtmp(θ)[SVector(1,2),:][i]
        dRdθ(θ) = dRdθtmp(θ)[SVector(3,4),:][i]

        @test dRfcn.(θs) ≈ dR.(θs)
        @test maxabsdiff(dRdθ.(θs), dRdθref) < 1e-6
    end
end

function test_integratedreflection_vertical(scenario)
    @unpack tx, ground, bfield, species = scenario

    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(tx.frequency, waveguide)

    R = LMP.integratedreflection(me)

    @test R[1,2] ≈ R[2,1]

    # Let's also check with interpolation susceptibilityfcn here
    Mfcn = z -> LMP.susceptibility(z, me)
    R2 = @inferred LMP.integratedreflection(me; susceptibilityfcn=Mfcn)

    Mfcn2 = LMP.susceptibilityspline(me)
    R3 = @inferred LMP.integratedreflection(me; susceptibilityfcn=Mfcn2)

    @test R2 == R
    @test maxabsdiff(R, R3) < 1e-6
end

function test_integratedreflection_deriv(scenario)
    @unpack tx, ground, bfield, species = scenario
    freq = tx.frequency

    # 1e-10 is hardcoded in Dθ form of `integratedreflection`
    # This function is a test of the Dθ form - the lower default tolerance in the non-Dθ
    # form would cause tests to fail if not lowered to the same threshold
    ip = IntegrationParams(tolerance=1e-10)
    params = LMPParams(integrationparams=ip)
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(freq, waveguide)

    Rref(θ) = LMP.integratedreflection(LMP.setea(θ, me); params=params)
    RdR(θ) = LMP.integratedreflection(LMP.setea(θ, me), LMP.Dθ(); params=params)

    Rs = Vector{SMatrix{2,2,ComplexF64,4}}(undef, length(θs))
    dRs = similar(Rs)
    Rrefs = similar(Rs)
    dRrefs = similar(Rs)

    Threads.@threads for i in eachindex(θs)
        v = RdR(θs[i])
        Rs[i] = v[SVector(1,2),:]
        dRs[i] = v[SVector(3,4),:]
        Rrefs[i] = Rref(θs[i])
        dRrefs[i] = FiniteDiff.finite_difference_derivative(Rref, θs[i], Val{:central})
    end

    for i = 1:4
        R = getindex.(Rs, i)
        dR = getindex.(dRs, i)
        Rr = getindex.(Rrefs, i)
        dRr = getindex.(dRrefs, i)

        # maxabsdiff criteria doesn't capture range of R and dR so rtol is used
        @test isapprox(R, Rr; rtol=1e-5)
        @test isapprox(dR, dRr; rtol=1e-3)
    end
end

function test_fresnelreflection(scenario)
    @unpack ea, tx, bfield, species, ground = scenario

    # PEC ground
    pec_ground = LMP.Ground(1, 1e12)
    vertical_ea = LMP.EigenAngle(π/2)
    Rg = LMP.fresnelreflection(vertical_ea, pec_ground, Frequency(24e3))
    @test isapprox(abs.(Rg), I; atol=1e-7)

    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)
    @test LMP.fresnelreflection(ea, ground, tx.frequency) == LMP.fresnelreflection(me)
end

function test_fresnelreflection_deriv(scenario)
    @unpack ea, tx, bfield, species, ground = scenario
    freq = tx.frequency

    dRgtmp(θ) = (ea = EigenAngle(θ); LMP.fresnelreflection(ea, ground, freq, LMP.Dθ()))
    for i = 1:4
        Rgref(θ) = (ea = EigenAngle(θ); LMP.fresnelreflection(ea, ground, freq)[i])
        dRgref = FiniteDiff.finite_difference_derivative(Rgref, θs, Val{:central})
        Rg(θ) = dRgtmp(θ)[1][i]
        dRg(θ) = dRgtmp(θ)[2][i]

        @test Rgref.(θs) ≈ Rg.(θs)
        @test maxabsdiff(dRg.(θs), dRgref) < 1e-6
    end

    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)
    Rg1 = LMP.fresnelreflection(ea, ground, freq, LMP.Dθ())
    Rg2 = LMP.fresnelreflection(me, LMP.Dθ())
    @test Rg1 == Rg2
end

function test_solvemodalequation(scenario)
    @unpack ea, tx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)

    f = LMP.solvemodalequation(me)
    me2 = PhysicalModeEquation(EigenAngle(0.0+0.0im), tx.frequency, waveguide)
    f2 = LMP.solvemodalequation(ea.θ, me2)
    @test f == f2

    # Test with specified susceptibilityfcn
    Mfcn = z -> LMP.susceptibility(z, me)
    f3 = @inferred LMP.solvemodalequation(me; susceptibilityfcn=Mfcn)

    Mfcn2 = LMP.susceptibilityspline(me)
    f4 = @inferred LMP.solvemodalequation(me; susceptibilityfcn=Mfcn2)

    @test f == f3
    @test maxabsdiff(f, f4) < 1e-3
end

function test_modalequation_resonant(scenario)
    @unpack ea, tx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    me = PhysicalModeEquation(ea, tx.frequency, waveguide)

    R = @SMatrix [1 0; 0 1]
    Rg = @SMatrix [1 0; 0 -1]
    @test isapprox(LMP.modalequation(R, Rg), 0; atol=1e-15)

    R = LMP.integratedreflection(me)
    Rg = LMP.fresnelreflection(ea, ground, tx.frequency)
    f = LMP.modalequation(R, Rg)
    @test isroot(f; atol=1e-4)  # this test is a little cyclic

    @test f == LMP.solvemodalequation(me)

    θ = 1.5 - 0.02im  # not resonant
    @test abs(LMP.solvemodalequation(θ, me)) > 1e-3
end

function test_modalequation_deriv(scenario)
    @unpack ea, tx, ground, bfield, species = scenario
    freq = tx.frequency

    waveguide = HomogeneousWaveguide(bfield, species, ground)
    modeequation = PhysicalModeEquation(ea, freq, waveguide)

    # `solvedmodalequation` calls the `Dθ` form of `integratedreflection`, which uses a
    # hard coded tolerance of 1e-10. Tests will fail if comparing to the lower default
    # tolerance used by `solvemodalequation`
    ip = IntegrationParams(tolerance=1e-10)
    params = LMPParams(integrationparams=ip)

    dFref = FiniteDiff.finite_difference_derivative(θ->LMP.solvemodalequation(θ, modeequation; params=params),
        θs, Val{:central})

    dFs = Vector{ComplexF64}(undef, length(θs))
    Threads.@threads for i in eachindex(θs)
        dFdθ, R, Rg = LMP.solvedmodalequation(θs[i], modeequation)
        dFs[i] = dFdθ
    end

    @test isapprox(dFs, dFref; rtol=1e-3)
end

function test_findmodes(scenario)
    @unpack tx, bfield, species, ground = scenario
    waveguide = HomogeneousWaveguide(bfield, species, ground)
    modeequation = PhysicalModeEquation(tx.frequency, waveguide)

    origcoords = LMP.defaultmesh(tx.frequency)

    # params = LMPParams(grpfparams=LMP.GRPFParams(100000, 1e-6, true))
    params = LMPParams()
    modes = @inferred findmodes(modeequation, origcoords; params=params)
    modes2 = @inferred findmodes(modeequation, origcoords; params=LMPParams(params; approxsusceptibility=true))

    @test length(modes) == length(modes2)

    # 5e-4 needed to accomodate multiplespecies_scenario. Otherwise they satisfy 1e-5
    @test all(maxabsdiff(modes[i].θ, modes2[i].θ) < 5e-4 for i in 1:length(modes))

    for m in modes
        f = LMP.solvemodalequation(m, modeequation; params=params)
        isroot(f) || return f
        @test isroot(f)
    end

    # return modes
end

########
# `TEST_MODES` must be filled for tests below!

function test_roots(scenario)
    x = 0.0005
    z = complex(x, -x)

    # isroot Real and Complex
    @test isroot(x)
    @test isroot(x; atol=1e-6) == false
    @test isroot(z)
    @test isroot(z; atol=1e-6) == false

    z2 = complex(1, 0)
    z3 = complex(0, 1)
    @test isroot(z2) == false
    @test isroot(z3) == false
end

# function evalroot(root, scenario)
#     @unpack tx, bfield, species, ground = scenario
#     waveguide = HomogeneousWaveguide(bfield, species, ground)
#
#     modeequation = PhysicalModeEquation(tx.frequency, waveguide)
#     LMP.solvemodalequation(EigenAngle(root), modeequation)
# end

@testset "modefinder.jl" begin
    @info "Testing modefinder"

    for scn in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
        multiplespecies_scenario)

        test_physicalmodeequation(scn)  # just test with one scenario?

        test_wmatrix(scn)
        test_wmatrix_deriv(scn)

        test_dRdz(scn)
        test_dRdθdz(scn)

        test_integratedreflection_deriv(scn)

        test_fresnelreflection(scn)
        test_fresnelreflection_deriv(scn)

        test_solvemodalequation(scn)
        test_modalequation_deriv(scn)

        test_findmodes(scn)
    end

    # Fill in TEST_MODES
    @info "  Mode finding..."
    for scenario in (verticalB_scenario, resonant_scenario, nonresonant_scenario,
        multiplespecies_scenario)
        if !haskey(TEST_MODES, scenario)
            TEST_MODES[scenario] = findroots(scenario)
        end
    end

    for scn in (resonant_scenario, )
        test_roots(scn)
        test_modalequation_resonant(scn)
    end
    for scn in (verticalB_scenario, )
        test_integratedreflection_vertical(scn)
    end
end


text:
module rdma
	export libibverbs, librdmacm

	module libibverbs
		using sys

		@pkgconf begin
			using libc.POSIX, linux

			c`-libverbs`
			c"""
				#include <infiniband/verbs.h>
			"""jiw
		end

		function __init__()
			ENV["IBV_FORK_SAFE"] = "1"
		end
	end

	module librdmacm
		using sys

		@pkgconf begin
			using libc.POSIX, linux, ..libibverbs

			c`-lrdmacm`
			c"""
				#include <rdma/rdma_verbs.h>
			"""ji
		end

		function __init__()
			ENV["RDMAV_FORK_SAFE"] = "1"
		end
	end
end


text:
let
    n = 50
    p = 3
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);
    Wt,c = egrss.potrf(Ut,Vt,1e-4*ones(n))

    Lref = egrss.full_tril(Ut,Wt,c)
    nrmref = sum(Lref.^2,dims=1)'
    nrm = egrss.trnrms(Ut,Wt,c)
    @test ( isapprox(nrm,nrmref) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

const prob = DDEProblemLibrary.prob_dde_constant_1delay_long_ip

@testset for constrained in (false, true)
  sol = solve(prob, MethodOfSteps(Tsit5(); constrained = constrained))

  @test allunique(sol.t)
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
###############################################################################
#
# The primary module for computing the Riemann theta function.
#
# .. math::
#
#   Θ(z, ω) = Σ
    # \theta(z, \Omega) = \sum_{n \in \mathbb{Z}^g}
                        # e^{2 \pi i \left( \tfrac{1}{2} n \cdot \Omega n
                           # + n \cdot z \right)}
#
#
#  Original Authors
#  -------
#  * Chris Swierczewski (@cswiercz) - September 2012, July 2016
#  * Grady Williams (@gradyrw) - October 2012
#  * Jeremy Upsal (@jupsal) - July 2016
#
# References
# ----------
#
# .. [CRTF] B. Deconinck, M.  Heil, A. Bobenko, M. van Hoeij and M. Schmies,
#    Computing Riemann Theta Functions, Mathematics of Computation, 73, (2004),
#    1417-1442.
#
# .. [DLMF] B. Deconinck, Digital Library of Mathematics Functions - Riemann
#    Theta Functions, http://dlmf.nist.gov/21
#
# .. [SAGE] Computing Riemann theta functions in Sage with applications.
#    C. Swierczewski and B. Deconinck.Submitted for publication.  Available
#    online at
#    http://depts.washington.edu/bdecon/papers/pdfs/Swierczewski_Deconinck1.pdf
#
###############################################################################

"""
         oscillatory_part(zs::Vector{Vector{ComplexF64}},
                          Ω::Matrix{ComplexF64};
                          eps::Float64=1e-8,
                          derivs::Vector{Vector{ComplexF64}}=Vector{ComplexF64}[],
                          accuracy_radius::Float64=5.)::Vector{ComplexF64}

Return the value of the oscillatory part of the Riemann theta function for Ω and
all z in `zs` if `derivs` is empty, or the derivatives at all z in `zs` for the
given directional derivatives in `derivs`.

Parameters
----------
- `zs` : A vector of complex vectors at which to evaluate the Riemann theta function.
- `Omega` : A Riemann matrix.
- `eps` : (Default: 1e-8) The desired numerical accuracy.
- `derivs` : A vector of complex vectors giving a directional derivative.
- `accuracy_radius` : (Default: 5.) The radius from the g-dimensional origin
where the requested accuracy of the Riemann theta is guaranteed when computing
derivatives. Not used if no derivatives of theta are requested.

Returns
-------
- The value of the oscillatory part of the Riemann theta function at each point appearing in `z`.
"""
function oscillatory_part(zs::Vector{Vector{ComplexF64}},
                          Ω::Matrix{ComplexF64};
                          eps::Float64=1e-8,
                          derivs::Vector{Vector{ComplexF64}}=Vector{ComplexF64}[],
                          accuracy_radius::Float64=5.)
    # extract the requested information: the real part, inverse of the
    # imaginary part, and the cholesky decomposition of the imaginary part
    X = real.(Ω)
    Y = imag.(Ω)

    # In python version numpy.linalg.cholesky returns the lower triangular
    #  matrix, which is then transposed. Julia's cholesky returns the upper
    #  triangular matrix, hence no need to transpose.
    T = Matrix(cholesky(Y).U)

    Yinv = inv(Y)

    # compute the integer points over which we approximate the infinite sum to
    # the requested accuracy
    R = radius(eps, T, derivs, accuracy_radius)
    S = innerpoints(T, R)

    finite_sum(X, Yinv, T, zs, S, derivs)
end

"""
         exponential_part(zs::Vector{Vector{ComplexF64}},
                          Ω::Matrix{ComplexF64})::Vector{Float64}

Return the value of the exponential part of the Riemann theta function for Ω and
all z in `zs`.

Parameters
----------
- `zs` : A vector of complex vectors at which to evaluate the Riemann theta function.
- `Omega` : A Riemann matrix.

Returns
-------
The value of the exponential part of the Riemann theta function at
each point appearing in `zs`.

"""
function exponential_part(zs::Vector{Vector{ComplexF64}},
                          Ω::Matrix{ComplexF64})::Vector{Float64}
    # extract the imaginary parts of z and the inverse of the imaginary part
    # of Omega
    y = [ imag.(z) for z in zs ]
    Yinv = inv(imag.(Ω))

    # apply the quadratic form to each vector in z
    map(y -> π * dot(y, Yinv * y), y)
end

"""
     riemanntheta(zs::Vector{Vector{ComplexF64}},
                  Ω::Matrix{ComplexF64};
                  eps::Float64=1e-8,
                  derivs::Vector{Vector{ComplexF64}}=Vector{ComplexF64}[],
                  accuracy_radius::Float64=5.)::Vector{ComplexF64}

Return the value of the Riemann theta function for Ω and all z in `zs` if
`derivs` is empty, or the derivatives at all z in `zs` for the given directional
derivatives in `derivs`.

Parameters
----------
- `zs` : A vector of complex vectors at which to evaluate the Riemann theta function.
- `Omega` : A Riemann matrix.
- `eps` : (Default: 1e-8) The desired numerical accuracy.
- `derivs` : A vector of complex vectors giving a directional derivative.
- `accuracy_radius` : (Default: 5.) The radius from the g-dimensional origin
where the requested accuracy of the Riemann theta is guaranteed when computing
derivatives. Not used if no derivatives of theta are requested.

Returns
-------
The value (or derivative) of the Riemann theta function at each point in `zs`.
"""
function riemanntheta(zs::Vector{Vector{ComplexF64}},
                      Ω::Matrix{ComplexF64};
                      eps::Float64=1e-8,
                      derivs::Vector{Vector{ComplexF64}}=Vector{ComplexF64}[],
                      accuracy_radius::Float64=5.)::Vector{ComplexF64}

    u = exponential_part(zs, Ω)
    v = oscillatory_part(zs, Ω, eps=eps, derivs=derivs,
                         accuracy_radius=accuracy_radius)

    exp.(u) .* v
end


text:
using ConstrainedDynamics
using ConstrainedControl
using LinearAlgebra

path = "examples/examples_files/sawyer_arm.urdf"
mech = Mechanism(path, floating=false, g = 0.0)

setPosition!(mech,mech.eqconstraints["right_j0"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j1"],[-pi/2])
setPosition!(mech,mech.eqconstraints["right_j2"],[-pi/2])
setPosition!(mech,mech.eqconstraints["right_j3"],[pi/2])

xd=[mech.bodies[i].state.xc for i=1:length(mech.bodies)]
qd=[mech.bodies[i].state.qc for i=1:length(mech.bodies)]

setPosition!(mech,mech.eqconstraints["right_j0"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j1"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j2"],[0.0])
setPosition!(mech,mech.eqconstraints["right_j3"],[0.0])

Q = [diagm(ones(12))*1000.0 for i=1:7]
R = [ones(1,1) for i=1:7]

ConstrainedDynamics.deactivateConstraints!(mech)
# A, Bu, Bλ, G = linearsystem(mech, xd, [zeros(3) for i=1:length(xd)], qd, [zeros(3) for i=1:length(xd)], [[0.] for i=1:length(mech.eqconstraints)], getid.(mech.bodies), getid.(mech.eqconstraints))
lqr = LQR(mech, getid.(mech.bodies), getid.(mech.eqconstraints), Q, R, 20., xd=xd, qd=qd)
ConstrainedDynamics.activateConstraints!(mech)
@test true


text:
module OpenTelemetryAPI

include("context.jl")
include("propagator.jl")
include("trace.jl")

end # module


text:
using DiffEqBase
mutable struct DummySolution
  retcode
end

mutable struct DummyIntegrator{Alg, IIP, U, T} <: DiffEqBase.DEIntegrator{Alg, IIP, U, T}
  uprev::U
  tprev::T
  u::U
  t::T
  dt::T
  tdir
  tstops
  sol::DummySolution

  DummyIntegrator() = new{Bool,Bool,Vector{Float64},Float64}([0.0],0,[0.0],0,1,1,[],DummySolution(:Default))
end

function DiffEqBase.add_tstop!(integrator::DummyIntegrator,t)
  integrator.tdir * (t - integrator.t) < 0 && error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
  push!(integrator.tstops,t)
end

function DiffEqBase.step!(integrator::DummyIntegrator)
  t_next = integrator.t + integrator.dt
  if ! isempty(integrator.tstops) && integrator.tstops[1] < t_next
    t_next = integrator.tstops[1]
  end
  integrator.uprev .= integrator.u
  integrator.u[1] += 2 * (t_next - integrator.t)
  integrator.tprev = integrator.t
  integrator.t = t_next
end

function step_dt!(integrator, args...)
  t = integrator.t
  step!(integrator, args...)
  integrator.t - t
end

function DiffEqBase.done(integrator::DummyIntegrator)
  integrator.t > 10
end

integrator = DummyIntegrator()
@test step_dt!(integrator, 1.5) == 2
@test step_dt!(integrator, 1.5, true) == 1.5
@test_throws ErrorException step!(integrator, -1)

for (u, t) in tuples(DummyIntegrator())
  @test u[1] == 2*t
end
@test eltype(collect(tuples(DummyIntegrator()))) == Tuple{Vector{Float64},Float64}

for (uprev, tprev, u, t) in intervals(DummyIntegrator())
  @test u[1] - uprev[1] == 2
  @test t - tprev == 1
end
@test eltype(collect(intervals(DummyIntegrator()))) == Tuple{Vector{Float64},Float64,Vector{Float64},Float64}


text:
using IndirectTrajOpt
using DifferentialEquations
using StaticArrays
using LinearAlgebra

# Initialize BVP function
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
tspan = (0.0, 8.6404*24*3600/ps.crp.TU)

# Error producing initial conditions
y0 =   @SVector [-0.0194885115,
                 -0.0160334798,
                 0.0,
                 8.9188819237,
                 -4.081936888,
                 0.0,
                 1.0,
                 -18.20557373360214,
                 -16.683730836929307,
                 40.0,
                 0.04174822756326783,
                 -0.005163201648349852,
                 -0.03929039276763477,
                 0.12683528760241147]

# Set thrust type
cSc = ps.sp.isp*9.81*ps.crp.TU / (ps.crp.LU*1000.0)
λv = norm(view(y0, 11:13))
S = IndirectTrajOpt.computeS(y0, λv, cSc)
if S > ps.ϵ; ps.utype = 0
elseif S < -ps.ϵ; ps.utype = 2
else; ps.utype = 1; end

cb = VectorContinuousCallback(
            IndirectTrajOpt.cr3bpEomsCondition,
            IndirectTrajOpt.cr3bpEomsAffect!,
            IndirectTrajOpt.cr3bpEomsAffect!, 4;
            idxs = nothing,
            rootfind = DiffEqBase.LeftRootFind,
            interp_points = 10,
            abstol = 1e-14,
            reltol = 0.0,
            save_positions = (true, true))
ff = ODEFunction{false}(IndirectTrajOpt.cr3bpEomIndirect)
prob = ODEProblem(ff, y0, tspan, ps; callback=cb)

sol = solve(prob, Vern9(),reltol=1e-14,abstol=1e-14)


text:
function scheduling(jobs)
    perm = sortperm(jobs, by=x->first(x))
    ans = Vector{Char}()
    cfreeon = 0
    jfreeon = 0
    for job in jobs[perm]
        if first(job) >= cfreeon
            push!(ans, 'C')
            cfreeon = last(job)
        elseif first(job) >= jfreeon
            push!(ans, 'J')
            jfreeon = last(job)
        else
            return "IMPOSSIBLE"
        end
    end
    string(ans[invperm(perm)]...)
end

function solution()
    n_cases = parse(Int, readline())
    for x = 1:n_cases
        n_jobs = parse(Int, readline())
        jobs = Vector{Tuple{Int, Int}}()
        for i = 1:n_jobs
            push!(jobs, Tuple(parse(Int, t) for t in split(readline())))
        end
        println("Case #$x: $(scheduling(jobs))")
    end
end

solution()


text:
#+Documented_Item,Type
# A *central* structure containing documented item
#
struct Documented_Item
    _filename::String
    _tokenized::Tokenized
    _extracted_tag::Extract_Tag_Result
    _doc::Extracted_Comment
    _code::Union{Extracted_Function,
                 Extracted_Struct}
end

#+Documented_Item,Methods
#
filename(di::Documented_Item)::String = di._filename


text:
using Documenter
using Reproduce

makedocs(
    sitename = "Reproduce",
    format = Documenter.HTML(),
    modules = [Reproduce],
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Developing an Experiment"=>"manual/experiment.md"
            "Parallel Jobs"=>"manual/parallel.md"
        ],
        "Documentation" => [
            "Experiment"=>"docs/experiment.md",
            "Parallel"=>"docs/parallel.md",
            "Data Structure"=>"docs/parse.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

deploydocs(
    repo = "github.com/mkschleg/Reproduce.jl.git",
    devbranch = "master",
    versions = ["stable" => "v^"]
)


text:
push!(LOAD_PATH, abspath(@__DIR__,".."))
using Plots
using AttractorNetworksBase ; const A = AttractorNetworksBase
using Calculus
using BenchmarkTools, Cthulhu
using Random
using Test
using Statistics,LinearAlgebra
Random.seed!(0)

using Plots ; theme(:dark)
function plotvs(x::AbstractArray{<:Real},y::AbstractArray{<:Real})
  x,y=x[:],y[:]
  @info """
  The max differences between the two are $(extrema(x .-y ))
  """
  plt=plot()
  scatter!(plt,x,y;leg=false,ratio=1,color=:white)
  lm=xlims()
  plot!(plt,identity,range(lm...;length=3);linestyle=:dash,color=:yellow)
  return plt
end
function test_Jgradient(utest::Vector{R},ntw::A.BaseNetwork) where R
    Jalloc = A.JAlloc(ntw)
    gradW=similar(ntw.weights)
    gradu=similar(ntw.weights)
    J=similar(gradW)
    # analytic gradient
    A.jacobian!(J,gradW,gradu,utest,ntw,Jalloc)
    # numerical gradient weights
    grad_numW = similar(gradW)
    ntw2=copy(ntw)
    for ij in eachindex(ntw2.weights)
        function gradfun(w)
            ntw2.weights[ij] = w
            A.jacobian!(J,nothing,nothing,utest,ntw2,Jalloc)
            return J[ij]
        end
        grad_numW[ij] = Calculus.gradient(gradfun,ntw2.weights[ij])
    end
    # numerical gradient u
    grad_numu = similar(gradu)
    for lj in CartesianIndices(gradu)
        (j,l)=Tuple(lj)
        function gradfun(u)
            u2 = copy(utest)
            u2[l]=u
            A.jacobian!(J,nothing,nothing,u2,ntw,Jalloc)
            return J[lj]
        end
        grad_numu[lj] = Calculus.gradient(gradfun,utest[l])
    end
    return gradW,grad_numW,gradu,grad_numu
end

## linear stable dynamcs goes to zero
ne,ni = 59,33
ntot = ne+ni
ntw = A.BaseNetwork(ne,ni ; gfun=A.IOId{Float64}())
ntw.weights .= 0.3.*randn(ntot)
B = ntw.weights - I
ntw.external_input .= 0.0
ntw.membrane_taus .= 1.0
r_start = 30.0.*randn(ntot)

# check velocity
A.velocity(r_start,ntw)

u_end,r_end=A.run_network_to_convergence(ntw,r_start;veltol=1E-5)

vel_test=similar(u_end)
A.velocity!(vel_test,u_end,ntw.iofunction.(u_end),ntw)

t,rdyn,_=A.run_network(ntw,r_end,4.0)

plot(t,rdyn[4,:];leg=false)

##
fill!(ntw.weights,0.0)
u_start = randn(A.ndims(ntw))
r_start = ntw.iofunction.(u_start)

u_end,r_end=A.run_network_to_convergence(ntw,r_start)
@test all(isapprox.(u_end,ntw.external_input;atol=0.05))

ne,ni = 20,23
ntot = ne+ni
ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.02))
u_start = randn(A.ndims(ntw))
r_start = ntw.iofunction.(u_start)
u_end,r_end=A.run_network_to_convergence(ntw,r_start)
t,ut,rt=A.run_network(ntw,r_end,1.0)
@test all(isapprox.(u_end,ut[:,end];atol=0.05))

##

##
Wtest = A.make_wmat(80,20,A.WDalian)
heatmap(Wtest ; ratio=1)

##
pe,pi = A.populations(80,20)
mah = rand(100,100)

testee = pe*pe'
testie = pi*pe'
all( view(testee,pe,pe))
all( view(testie,pi,pe))

view(mah,pe*pe')

heatmap(testie ; ratio=1)
heatmap(pe*pi' ; ratio = 1)
heatmap(pi*pe' ; ratio = 1)

mah[pe,pe] .- mah[1:80,1:80]

##
ne,ni = 34,44
wmeanf = [ 0.333  -1.123
            12.122 -0.06432]
W = A.make_wmat(ne,ni,wmeanf;noautapses=true)

pe,pi = A.populations(ne,ni)

sum(W[pe,pe]; dims=2)
#
# whatevs = A.run_network(10.0 .* rand(ntot), 0.5,ntw)
#
# using Plots; theme(:dark)
# _ = let t=whatevs.t,
#   u12 = [ uu[1:2] for uu in whatevs.u]
#   u12=hcat(u12...) |> permutedims
#   plot(t,u12;leg=false,linewidth=4)
# end

##

ne,ni = 4,3
ntot = ne+ni
ntw = A.BaseNetwork(ne,ni ; gfun=A.IOQuad(0.02))

utest = randn(ntot)

grad_ana,grad_num = test_Jgradient_weights(utest,ntw)
grad_ana,grad_num = test_Jgradient_u(utest,ntw)

grad_ana ./ grad_num

grad_num

grad_ana


text:
M, ΔM = 0.0266, 0.0001
L, ΔL = 0.380, 0.001
L1,ΔL1= 0.025, 0.001
L2,ΔL2= 0.190, 0.001
d,Δd= 0.008,0.001

m,Δm = 0.0753, 0.0001
h,Δh = 0.02, 0.001
de,Δde = 0.025, 0.001
di,Δdi = 0.008, 0.001

# Barra

Ih = (1/12)*M*L^2 + (1/16)*M*de^2
ΔIh = (1/12)*(ΔM*L^2 + 2*M*L*ΔL) + (1/16)*(ΔM*d + 2*M*d*Δd)

Ih_eixo = Ih + M*(L2 - L1)^2
ΔIh_eixo = ΔIh + ΔM*(L2-L1)^2 + 2*M*(L2-L1)*(ΔL2+ΔL1)

# Cilindro

Icl = (1/16)*m*(de^2 + di^2) +  (1/12)*m*h^2
ΔIcl = (1/16)*(Δm*(de^2+di^2) + 2*m*(de*Δde+di*Δdi)) + (1/12)*(Δm*h^2 + 2*m*h*Δh)

Icl_cm = Icl + m*(L2 - h/2)^2
ΔIcl_cm = ΔIcl + Δm*(L2- h/2)^2 + 2*m*(L2-h/2)*(ΔL2+ Δh/2)

Icl_eixo = Icl + m*(L - L1 - h/2)^2
ΔIcl_eixo = ΔIcl + Δm*(L - L1 - h/2)^2 + 2*m*(L-L1-h/2)*(ΔL +ΔL1 + Δh/2)

println("\n Sistema no centro de massa\n")
println("======================")
println("Ih ± ΔIh =" , round(Ih,digits=8), " ± ", round(ΔIh,digits=8))
println("Icl_cm ± ΔIcl_cm =" , round(Icl_cm,digits=8), " ± ", round(ΔIcl_cm,digits=8))

println("\n Sistema no eixo \n ")
println("======================")
println("Ih_eixo ± ΔIh_eixo =" , round(Ih_eixo,digits=8), " ± ", round(ΔIh_eixo,digits=8))
println("Icl_eixo ± ΔIh_eixo =" , round(Icl_eixo,digits=8), " ± ",round(ΔIcl_eixo,digits=8))


text: include_folder(Main, @__FILE__, is_sorted=false, is_test=true)


text:
# This file was generated, do not modify it. # hide
using CategoricalArrays

v = categorical(["AA", "BB", "CC", "AA", "BB", "CC"])


text:
using Base.Test
# pointwiselayerstate.jl
import WaveletScattering: PointwiseLayerState
# inputlayer.jl
import WaveletScattering: InputLayer
# inputlayerstate.jl
import WaveletScattering: InputLayerState
# path.jl
import WaveletScattering: Path
# pointwiselayer.jl
import WaveletScattering: PointwiseLayer

backend = Mocha.CPUBackend()
data = map(Float32, randn(256, 2))
inputlayer = InputLayer(
        tops = [:signal],
        symbols = [:time, :chunk],
        data = data)
inputlayerstate = InputLayerState(backend, inputlayer)

bottoms = [:signal]
tops = [:identity]
pointwiselayer = PointwiseLayer(bottoms = bottoms, tops = tops)
pointwiselayerstate = Mocha.setup(backend, pointwiselayer,
    inputlayerstate.blobs, Mocha.Blob[Mocha.NullBlob()])

@test length(pointwiselayerstate.blobs) == 1
pointwisedata = pointwiselayerstate.blobs[1].nodes[Path()].data
@test pointwisedata ≈ data