text:
KWARGS_DEFAULT = Dict([
        (:train_sample_size, 0),
        (:val_sample_size, 0),
        (:extension, ".csv"),
        (:val_ratio, 0.0),
        (:n_grams_target_col, :Word),
        (:n_grams_tokenized, false),
        (:n_grams_sep_token, nothing),
        (:grams, 3),
        (:n_grams_keep_sep, false),
        (:start_end_token, "#"),
        (:path_sep_token, ":"),
        (:random_seed, 314),
        (:sd_base_mean, 1),
        (:sd_inflection_mean, 1),
        (:sd_base, 4),
        (:sd_inflection, 4),
        (:isdeep, true),
        (:add_noise, true),
        (:sd_noise, 1),
        (:normalized, false),
        (:if_combined, false),
        (:learn_mode, :cholesky),
        (:method, :additive),
        (:shift, 0.02),
        (:multiplier, 1.01),
        (:output_format, :auto),
        (:sparse_ratio, 0.05),
        (:wh_freq, nothing),
        (:init_weights, nothing),
        (:eta, 0.1),
        (:n_epochs, 1),
        (:max_t, 0),
        (:A, nothing),
        (:A_mode, :combined),
        (:max_can, 10),
        (:threshold_train, 0.1),
        (:is_tolerant_train, false),
        (:tolerance_train, -0.1),
        (:max_tolerance_train, 2),
        (:threshold_val, 0.1),
        (:is_tolerant_val, false),
        (:tolerance_val, -0.1),
        (:max_tolerance_val, 2),
        (:issparse, false),
        (:n_neighbors_train, 10),
        (:n_neighbors_val, 20),
        (:output_dir, "out"),
        (:verbose, false)
        ])

KWARGS_DUMP = Dict()

"""
    test_combo(test_mode;kwargs...)

A wrapper function for a full model for a specific combination of parameters. A detailed introduction is in [Test Combo Introduction](@ref)

# Obligatory Arguments
- `test_mode::Symbol`: which test mode, currently supports :train_only, :pre_split, :careful_split and :random_split.

# Optional Arguments
- `train_sample_size::Int64=0`: the desired number of training data
- `val_sample_size::Int64=0`: the desired number of validation data
- `val_ratio::Float64=0.0`: the desired portion of validation data, if works only if :val_sample_size is 0.0.
- `extension::String=".csv"`: the extension for data n_features_inflections
- `n_grams_target_col::Union{String, Symbol}=:Word`: the column name for target strings
- `n_grams_tokenized::Boolean=false`: if true, the dataset target is assumed to be tokenized
- `n_grams_sep_token::String=nothing`: separator
- `grams::Int64=3`: the number of grams for cues
- `n_grams_keep_sep::Boolean=false`: if true, keep separators in cues
- `start_end_token::String=":"`: start and end token in boundary cues
- `path_sep_token::String=":"`: path separator in the assembled path
- `random_seed::Int64=314`: the random seed
- `sd_base_mean::Int64=1`: the sd mean of base features
- `sd_inflection_mean::Int64=1`: the sd mean of inflectional features
- `sd_base::Int64=4`: the sd of base features
- `sd_inflection::Int64=4`: the sd of inflectional features
- `isdeep::Boolean=true`: if true, mean of each feature is also randomized
- `add_noise::Boolean=true`: if true, add additional Gaussian noise
- `sd_noise::Int64=1`: the sd of the Gaussian noise
- `normalized::Boolean=false`: if true, most of the values range between 1 and -1, it may slightly exceed between 1 or -1 depending on the sd
- `if_combined::Boolean=false`: if true, then features are combined with both training and validation data
- `learn_mode::Int64=:cholesky`: which learning mode, currently supports :cholesky and :wh
- `method::Int64=:additive`: whether :additive or :multiplicative decomposition is required
- `shift::Int64=0.02`: shift value for :additive decomposition
- `multiplier::Int64=1.01`: multiplier value for :multiplicative decomposition
- `output_format::Int64=:auto`: to force output format to dense(:dense) or sparse(:sparse), make it auto(:auto) to determined by the program
- `sparse_ratio::Int64=0.05`: the ratio to decide whether a matrix is sparse
- `wh_freq::Vector=nothing`: the learning sequence
- `init_weights::Matrix=nothing`: the initial weights
- `eta::Float64=0.1`: the learning rate
- `n_epochs::Int64=1`: the number of epochs to be trained
- `max_t::Int64=0`: the number of epochs to be trained
- `A::Matrix=nothing`: the number of epochs to be trained
- `A_mode::Symbol=:combined`: the adjacency matrix mode, currently supports :combined or :train_only
- `max_can::Int64=10`: the max number of candidate path to keep in the output
- `threshold_train::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 for training data
- `is_tolerant_train::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  for training data
- `tolerance_train::Float64=-0.1`: 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  for training data
- `max_tolerance_train::Int64=2`: maximum number of n-grams allowed in a path for training data
- `threshold_val::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 for validation data
- `is_tolerant_val::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 for validation data
- `tolerance_val::Float64=-0.1`: 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 for validation data
- `max_tolerance_val::Int64=2`: maximum number of n-grams allowed in a path for validation data
- `n_neighbors_train::Int64=10`: the top n form neighbors to be considered for training data
- `n_neighbors_val::Int64=20`: the top n form neighbors to be considered for validation data
- `issparse::Bool=false`: if true, keep sparse matrix format when learning paths
- `output_dir::String="out"`: the output directory
- `verbose::Bool=false`: if true, more information will be printed
"""
function test_combo(test_mode; kwargs...)
    verbose = get_kwarg(kwargs, :verbose, required=false)
    train_sample_size = get_kwarg(kwargs, :train_sample_size, required=false)
    val_sample_size = get_kwarg(kwargs, :val_sample_size, required=false)
    val_ratio = get_kwarg(kwargs, :val_ratio, required=false)
    n_grams_target_col = get_kwarg(kwargs, :n_grams_target_col, required=false)
    n_grams_tokenized = get_kwarg(kwargs, :n_grams_tokenized, required=false)
    n_grams_sep_token = get_kwarg(kwargs, :n_grams_sep_token, required=false)
    grams = get_kwarg(kwargs, :grams, required=false)
    n_grams_keep_sep = get_kwarg(kwargs, :n_grams_keep_sep, required=false)
    start_end_token = get_kwarg(kwargs, :start_end_token, required=false)
    path_sep_token = get_kwarg(kwargs, :path_sep_token, required=false)
    random_seed = get_kwarg(kwargs, :random_seed, required=false)
    if_combined = get_kwarg(kwargs, :if_combined, required=false)

    verbose && println("="^20)
    verbose && println("Preparing datasets...")
    verbose && println("="^20)

    # split and load data
    if test_mode == :train_only
        data_path = get_kwarg(kwargs, :data_path, required=true)

        data_train, data_val = loading_data_train_only(data_path,
            train_sample_size = train_sample_size,
            val_sample_size = val_sample_size)
    elseif test_mode == :pre_split
        data_path = get_kwarg(kwargs, :data_path, required=true)
        data_prefix = get_kwarg(kwargs, :data_prefix, required=true)
        extension = get_kwarg(kwargs, :extension, required=false)

        data_train, data_val = loading_data_pre_split(
            data_path, data_prefix,
            train_sample_size = train_sample_size,
            val_sample_size = val_sample_size, extension=extension)

    elseif test_mode == :random_split
        data_path = get_kwarg(kwargs, :data_path, required=true)
        data_output_dir = get_kwarg(kwargs, :data_output_dir, required=true)
        data_prefix = get_kwarg(kwargs, :data_prefix, required=true)

        data_train, data_val = loading_data_randomly_split(
            data_path,
            data_output_dir,
            data_prefix,
            train_sample_size = train_sample_size,
            val_sample_size = val_sample_size,
            val_ratio = val_ratio,
            random_seed = random_seed,
            verbose = verbose)

    elseif test_mode == :careful_split
        data_path = get_kwarg(kwargs, :data_path, required=true)
        data_prefix = get_kwarg(kwargs, :data_prefix, required=true)
        data_output_dir = get_kwarg(kwargs, :data_output_dir, required=true)
        n_features_columns = get_kwarg(kwargs, :n_features_columns, required=true)

        data_train, data_val = loading_data_careful_split(
            data_path,
            data_prefix,
            data_output_dir,
            n_features_columns,
            train_sample_size = train_sample_size,
            val_sample_size = val_sample_size,
            val_ratio=val_ratio,
            n_grams_target_col=n_grams_target_col,
            n_grams_tokenized=n_grams_tokenized,
            n_grams_sep_token=n_grams_sep_token,
            grams = grams,
            n_grams_keep_sep= n_grams_keep_sep,
            start_end_token = start_end_token,
            random_seed = random_seed,
            verbose=verbose)
    else
        throw(ArgumentError("test_mode is incorrect, using :train_only," *
            " :pre_split, :careful_split or :random_split"))
    end

    check_used_token(
        data_train,
        n_grams_target_col,
        start_end_token,
        "start_end_token",
    )

    check_used_token(
        data_val,
        n_grams_target_col,
        path_sep_token,
        "path_sep_token",
    )

    check_used_token(
        data_train,
        n_grams_target_col,
        start_end_token,
        "start_end_token",
    )

    check_used_token(
        data_val,
        n_grams_target_col,
        path_sep_token,
        "path_sep_token",
    )

    n_features_base = get_kwarg(kwargs, :n_features_base, required=true)
    n_features_inflections = get_kwarg(kwargs, :n_features_inflections, required=true)

    sd_base_mean = get_kwarg(kwargs, :sd_base_mean, required=false)
    sd_inflection_mean = get_kwarg(kwargs, :sd_inflection_mean, required=false)
    sd_base = get_kwarg(kwargs, :sd_base, required=false)
    sd_inflection = get_kwarg(kwargs, :sd_inflection, required=false)
    isdeep = get_kwarg(kwargs, :isdeep, required=false)
    add_noise = get_kwarg(kwargs, :add_noise, required=false)
    sd_noise = get_kwarg(kwargs, :sd_noise, required=false)
    normalized = get_kwarg(kwargs, :normalized, required=false)

    # make cue matrix/matrices
    # make semantic matrix/matrices
    verbose && println("="^20)
    verbose && println("Making cue matrix object...")
    verbose && println("="^20)

    cue_obj_train, cue_obj_val = make_cue_train_val(
        data_train,
        data_val,
        grams,
        n_grams_target_col,
        n_grams_tokenized,
        n_grams_sep_token,
        n_grams_keep_sep,
        start_end_token,
        if_combined,
        verbose)

    verbose && println("="^20)
    verbose && println("Making S matrix...")
    verbose && println("="^20)

    n_features = size(cue_obj_train.C, 2)
    S_train, S_val = make_S_train_val(data_train, data_val,
        n_features_base, n_features_inflections,
        n_features, sd_base_mean, sd_inflection_mean, sd_base,
        sd_inflection, isdeep, add_noise, sd_noise,
        normalized, if_combined, verbose)

    # temporary fix, S_val is not valid using :train_only
    # add noise don't apply to both S
    if test_mode == :train_only
        if val_sample_size == 0
            val_sample_size = 2
        end
        S_val = S_train[1:val_sample_size, :]
    end

    verbose && println("="^20)
    verbose && println("Learning transformation mapping F and G...")
    verbose && println("="^20)

    learn_mode = get_kwarg(kwargs, :learn_mode, required=false)

    # cholesky params
    method = get_kwarg(kwargs, :method, required=false)
    shift = get_kwarg(kwargs, :shift, required=false)
    multiplier = get_kwarg(kwargs, :multiplier, required=false)
    output_format = get_kwarg(kwargs, :output_format, required=false)
    sparse_ratio = get_kwarg(kwargs, :sparse_ratio, required=false)

    # wh params
    wh_freq = get_kwarg(kwargs, :wh_freq, required=false)
    init_weights = get_kwarg(kwargs, :init_weights, required=false)
    eta = get_kwarg(kwargs, :eta, required=false)
    n_epochs = get_kwarg(kwargs, :n_epochs, required=false)

    if learn_mode == :cholesky
        F_train = make_transform_matrix(
            cue_obj_train.C,
            S_train,
            method = method,
            shift = shift,
            multiplier = multiplier,
            output_format = output_format,
            sparse_ratio = sparse_ratio,
            verbose = verbose,
        )

        G_train = make_transform_matrix(
            S_train,
            cue_obj_train.C,
            method = method,
            shift = shift,
            multiplier = multiplier,
            output_format = output_format,
            sparse_ratio = sparse_ratio,
            verbose = verbose,
        )
    elseif learn_mode == :wh
        learn_seq = make_learn_seq(wh_freq)

        F_train = wh_learn(
            cue_obj_train.C,
            S_train;
            eta = eta,
            n_epochs = n_epochs,
            weights = init_weights,
            learn_seq = learn_seq,
            verbose = verbose,
        )

        G_train = wh_learn(
            S_train,
            cue_obj_train.C;
            eta = eta,
            n_epochs = n_epochs,
            weights = init_weights,
            learn_seq = learn_seq,
            verbose = verbose,
        )
    else
        throw(ArgumentError("learn_mode is incorrect, using :cholesky," *
            ":wh"))
    end

    verbose && println("="^20)
    verbose && println("Predicting S and C...")
    verbose && println("="^20)

    Shat_train = cue_obj_train.C * F_train
    Shat_val = cue_obj_val.C * F_train
    Chat_train = S_train * G_train
    Chat_val = S_val * G_train

    # calculate max timestep
    max_t = get_kwarg(kwargs, :max_t, required=false)

    if max_t == 0
        max_t = cal_max_timestep(data_train, data_val,
            n_grams_target_col, tokenized = n_grams_tokenized,
            sep_token = n_grams_sep_token)
    end

    # choose A
    A = get_kwarg(kwargs, :A, required=false)
    A_mode = get_kwarg(kwargs, :A_mode, required=false)

    if isnothing(A)
        if A_mode == :combined
            A = make_combined_adjacency_matrix(
                data_train,
                data_val;
                grams = grams,
                target_col = n_grams_target_col,
                tokenized = n_grams_tokenized,
                sep_token = n_grams_sep_token,
                keep_sep = n_grams_keep_sep,
                start_end_token = start_end_token,
                verbose = verbose
            )
        elseif A_mode == :train_only
            A = cue_obj_train.A
        else
            throw(ArgumentError("A_mode $A_mode is not supported!" *
                "Please choose from :combined or :train_only"))
        end
    end

    verbose && println("="^20)
    verbose && println("Performing path-finding algorithms...")
    verbose && println("="^20)

    max_can = get_kwarg(kwargs, :max_can, required=false)
    threshold_train = get_kwarg(kwargs, :threshold_train, required=false)
    is_tolerant_train = get_kwarg(kwargs, :is_tolerant_train, required=false)
    tolerance_train = get_kwarg(kwargs, :tolerance_train, required=false)
    max_tolerance_train = get_kwarg(kwargs, :max_tolerance_train, required=false)
    threshold_val = get_kwarg(kwargs, :threshold_val, required=false)
    is_tolerant_val = get_kwarg(kwargs, :is_tolerant_val, required=false)
    tolerance_val = get_kwarg(kwargs, :tolerance_val, required=false)
    max_tolerance_val = get_kwarg(kwargs, :max_tolerance_val, required=false)
    issparse = get_kwarg(kwargs, :issparse, required=false)
    sparse_ratio = get_kwarg(kwargs, :sparse_ratio, required=false)

    res_learn_train, gpi_learn_train = learn_paths(
        data_train,
        data_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_train.gold_ind,
        Shat_val = Shat_train,
        check_gold_path = true,
        max_t = max_t,
        max_can = max_can,
        threshold = threshold_train,
        is_tolerant = is_tolerant_train,
        tolerance = tolerance_train,
        max_tolerance = max_tolerance_train,
        grams = grams,
        tokenized = n_grams_tokenized,
        sep_token = n_grams_sep_token,
        keep_sep = n_grams_keep_sep,
        target_col = n_grams_target_col,
        start_end_token = start_end_token,
        issparse = issparse,
        sparse_ratio = sparse_ratio,
        verbose = verbose)

    res_learn_val, gpi_learn_val = learn_paths(
        data_train,
        data_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_train.i2f,
        cue_obj_train.f2i,
        gold_ind = cue_obj_val.gold_ind,
        Shat_val = Shat_val,
        check_gold_path = true,
        max_t = max_t,
        max_can = max_can,
        threshold = threshold_val,
        is_tolerant = is_tolerant_val,
        tolerance = tolerance_val,
        max_tolerance = max_tolerance_val,
        grams = grams,
        tokenized = n_grams_tokenized,
        sep_token = n_grams_sep_token,
        keep_sep = n_grams_keep_sep,
        target_col = n_grams_target_col,
        start_end_token = start_end_token,
        issparse = issparse,
        sparse_ratio = sparse_ratio,
        verbose = verbose)

    n_neighbors_train = get_kwarg(kwargs, :n_neighbors_train, required=false)
    n_neighbors_val = get_kwarg(kwargs, :n_neighbors_val, required=false)

    rC_train = cor(Chat_train, Matrix(cue_obj_train.C), dims = 2)
    rC_val = cor(Chat_val, Matrix(cue_obj_train.C), dims = 2)

    res_build_train = build_paths(
        data_train,
        cue_obj_train.C,
        S_train,
        F_train,
        Chat_train,
        A,
        cue_obj_train.i2f,
        cue_obj_train.gold_ind;
        rC = rC_train,
        max_t = max_t,
        max_can = max_can,
        n_neighbors = n_neighbors_train,
        grams = grams,
        tokenized = n_grams_tokenized,
        sep_token = n_grams_sep_token,
        target_col = n_grams_target_col,
        verbose = verbose,
    )

    res_build_val = build_paths(
        data_val,
        cue_obj_train.C,
        S_val,
        F_train,
        Chat_val,
        A,
        cue_obj_val.i2f,
        cue_obj_train.gold_ind;
        rC = rC_val,
        max_t = max_t,
        max_can = max_can,
        n_neighbors = n_neighbors_val,
        grams = grams,
        tokenized = n_grams_tokenized,
        sep_token = n_grams_sep_token,
        target_col = n_grams_target_col,
        verbose = verbose,
    )

    verbose && println("="^20)
    verbose && println("Evaluating results...")
    verbose && println("="^20)

    acc_Chat_train = eval_SC(Chat_train, cue_obj_train.C)
    acc_Shat_train = eval_SC(Shat_train, S_train)
    acc_Shat_train_homo = eval_SC(Shat_train, S_train, data_train, n_grams_target_col)

    acc_Chat_val = eval_SC(Chat_val, cue_obj_val.C)
    acc_Chat_val_tv = eval_SC(Chat_val, cue_obj_val.C, cue_obj_train.C)

    acc_Shat_val = eval_SC(Shat_val, S_val)
    acc_Shat_val_tv = eval_SC(Shat_val, S_val, S_train)

    acc_Shat_val_homo = eval_SC(Shat_val, S_val, data_val, n_grams_target_col)
    acc_Shat_val_homo_tv = eval_SC(Shat_val, S_val, S_train, data_val, data_train, n_grams_target_col)

    acc_learn_train = JudiLing.eval_acc(
        res_learn_train,
        cue_obj_train.gold_ind,
        verbose=verbose
    )

    acc_learn_val = JudiLing.eval_acc(
        res_learn_val,
        cue_obj_val.gold_ind,
        verbose=verbose
    )

    acc_build_train = JudiLing.eval_acc(
        res_build_train,
        cue_obj_train.gold_ind,
        verbose=verbose
    )

    acc_build_val = JudiLing.eval_acc(
        res_build_val,
        cue_obj_val.gold_ind,
        verbose=verbose
    )

    verbose && println("="^20)
    verbose && println("Saving outputs...")
    verbose && println("="^20)

    output_dir = get_kwarg(kwargs, :output_dir, required=false)

    # write acc output
    mkpath(output_dir)

    save_S_matrix(S_train, joinpath(".", output_dir, "S_train.csv"),
        data_train, n_grams_target_col)
    save_S_matrix(S_val, joinpath(".", output_dir, "S_val.csv"), data_val,
        n_grams_target_col)

    accio = open(joinpath(output_dir, "acc.out"), "w")
    println(accio, "Acc for Chat train: $acc_Chat_train")
    println(accio, "Acc for Shat train: $acc_Shat_train")
    println(accio, "Acc for Shat train homophones: $acc_Shat_train_homo")
    println(accio, "Acc for Chat val: $acc_Chat_val")
    println(accio, "Acc for Chat val for both train and val: $acc_Chat_val_tv")
    println(accio, "Acc for Shat val: $acc_Shat_val")
    println(accio, "Acc for Acc for Shat val for both train and val: $acc_Shat_val_tv")
    println(accio, "Acc for Shat val homophones: $acc_Shat_val_homo")
    println(accio, "Acc for Shat val homophones for both train and val: $acc_Shat_val_homo_tv")
    println(accio, "Acc for learn_path train: $acc_learn_train")
    println(accio, "Acc for learn_path val: $acc_learn_val")
    println(accio, "Acc for build_path train: $acc_build_train")
    println(accio, "Acc for build_path val: $acc_build_val")
    close(accio)

    # write params into a file
    params_io = open(joinpath(output_dir, "params.out"), "w")
    for (k,v) in KWARGS_DUMP
        println(params_io, "$k: $v")
    end
    close(params_io)

    write2csv(
        res_learn_train,
        data_train,
        cue_obj_train,
        cue_obj_train,
        "learn_res_train.csv",
        grams=grams,
        tokenized=n_grams_tokenized,
        sep_token=n_grams_sep_token,
        start_end_token=start_end_token,
        output_sep_token=n_grams_sep_token,
        path_sep_token=path_sep_token,
        target_col=n_grams_target_col,
        root_dir=".",
        output_dir=output_dir
        )

    write2csv(
        gpi_learn_train,
        "learn_gpi_train.csv",
        root_dir=".",
        output_dir=output_dir
        )

    write2csv(
        res_learn_val,
        data_val,
        cue_obj_train,
        cue_obj_val,
        "learn_res_val.csv",
        grams=grams,
        tokenized=n_grams_tokenized,
        sep_token=n_grams_sep_token,
        start_end_token=start_end_token,
        output_sep_token=n_grams_sep_token,
        path_sep_token=path_sep_token,
        target_col=n_grams_target_col,
        root_dir=".",
        output_dir=output_dir
        )

    write2csv(
        gpi_learn_val,
        "learn_gpi_val.csv",
        root_dir=".",
        output_dir=output_dir
        )

    write2csv(
        res_build_train,
        data_train,
        cue_obj_train,
        cue_obj_train,
        "build_res_train.csv",
        grams=grams,
        tokenized=n_grams_tokenized,
        sep_token=n_grams_sep_token,
        start_end_token=start_end_token,
        output_sep_token=n_grams_sep_token,
        path_sep_token=path_sep_token,
        target_col=n_grams_target_col,
        root_dir=".",
        output_dir=output_dir
        )

    write2csv(
        res_build_val,
        data_val,
        cue_obj_train,
        cue_obj_val,
        "build_res_val.csv",
        grams=grams,
        tokenized=n_grams_tokenized,
        sep_token=n_grams_sep_token,
        start_end_token=start_end_token,
        output_sep_token=n_grams_sep_token,
        path_sep_token=path_sep_token,
        target_col=n_grams_target_col,
        root_dir=".",
        output_dir=output_dir
        )

end

function loading_data_train_only(
    data_path;
    train_sample_size = 0,
    val_sample_size = 0)

    data = DataFrame(CSV.File(data_path))

    if train_sample_size != 0
        data = data[1:train_sample_size, :]
    end

    if val_sample_size == 0
        val_sample_size = 2
    end

    data, data[1:val_sample_size, :]
end

function loading_data_pre_split(
    data_output_dir,
    data_prefix;
    train_sample_size = 0,
    val_sample_size = 0,
    extension=".csv")

    train_path = joinpath(data_output_dir, data_prefix * "_train" * extension)
    val_path = joinpath(data_output_dir, data_prefix * "_val" * extension)
    data_train = DataFrame(CSV.File(train_path))
    data_val = DataFrame(CSV.File(val_path))

    if train_sample_size != 0
        data_train = data_train[1:train_sample_size, :]
    end

    if val_sample_size != 0
        data_val = data_val[1:val_sample_size, :]
    end

    data_train, data_val
end

function loading_data_randomly_split(
    data_path,
    output_dir_path,
    data_prefix;
    train_sample_size = 0,
    val_sample_size = 0,
    val_ratio = 0.0,
    random_seed = 314,
    verbose = false)
    verbose && println("Spliting data...")

    train_val_random_split(
        data_path,
        output_dir_path,
        data_prefix,
        train_sample_size = train_sample_size,
        val_sample_size = val_sample_size,
        val_ratio = val_ratio,
        random_seed = random_seed,
        verbose = verbose,
        )

    # load data
    verbose && println("Loading CSV...")
    loading_data_pre_split(output_dir_path, data_prefix)
end

function loading_data_careful_split(
    data_path,
    data_prefix,
    output_dir_path,
    n_features_columns;
    train_sample_size = 0,
    val_sample_size = 0,
    val_ratio = 0.0,
    n_grams_target_col = :Word,
    n_grams_tokenized = false,
    n_grams_sep_token = nothing,
    grams = 3,
    n_grams_keep_sep = false,
    start_end_token = "#",
    random_seed = 314,
    verbose = false)

    verbose && println("Spliting data...")
    train_val_careful_split(
        data_path,
        output_dir_path,
        data_prefix,
        n_features_columns,
        train_sample_size = train_sample_size,
        val_sample_size = val_sample_size,
        val_ratio = val_ratio,
        n_grams_target_col = n_grams_target_col,
        n_grams_tokenized = n_grams_tokenized,
        n_grams_sep_token = n_grams_sep_token,
        grams = grams,
        n_grams_keep_sep = n_grams_keep_sep,
        start_end_token = start_end_token,
        random_seed = random_seed,
        verbose = verbose,
    )

    # load data
    verbose && println("Loading CSV...")
    loading_data_pre_split(output_dir_path, data_prefix)
end

function make_cue_train_only(data, grams, target_col, tokenized, sep_token,
    keep_sep, start_end_token, verbose)

    make_cue_matrix(
        data,
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose
        )
end

function make_cue_train_val(data_train, data_val, grams, target_col, tokenized,
    sep_token, keep_sep, start_end_token, if_combined, verbose)

    if if_combined
        cue_obj_train, cue_obj_val = make_combined_cue_matrix(
            data_train,
            data_val,
            grams = grams,
            target_col = target_col,
            tokenized = tokenized,
            sep_token = sep_token,
            keep_sep = keep_sep,
            start_end_token = start_end_token,
            verbose = verbose
            )
    else
        cue_obj_train, cue_obj_val = make_cue_matrix(
            data_train,
            data_val,
            grams = grams,
            target_col = target_col,
            tokenized = tokenized,
            sep_token = sep_token,
            keep_sep = keep_sep,
            start_end_token = start_end_token,
            verbose = verbose
            )
    end
    cue_obj_train, cue_obj_val
end

function make_S_train_only(data, n_features_base, n_features_inflections,
    ncol, sd_base_mean, sd_inflection_mean, sd_base, sd_inflection,
    isdeep, add_noise, sd_noise, normalized, verbose)
    verbose && println("Making S matrix...")
    make_S_matrix(
        data,
        n_features_base,
        n_features_inflections,
        ncol = ncol,
        sd_base_mean = sd_base_mean,
        sd_inflection_mean = sd_inflection_mean,
        sd_base = sd_base,
        sd_inflection = sd_inflection,
        isdeep = isdeep,
        add_noise = add_noise,
        sd_noise = sd_noise,
        normalized = normalized
    )
end

function make_S_train_val(data_train, data_val,
    n_features_base, n_features_inflections,
    ncol, sd_base_mean, sd_inflection_mean, sd_base, sd_inflection,
    isdeep, add_noise, sd_noise, normalized, if_combined, verbose)
    verbose && println("Making S matrix...")
    if if_combined
        make_combined_S_matrix(
            data_train,
            data_val,
            n_features_base,
            n_features_inflections,
            ncol = ncol,
            sd_base_mean = sd_base_mean,
            sd_inflection_mean = sd_inflection_mean,
            sd_base = sd_base,
            sd_inflection = sd_inflection,
            isdeep = isdeep,
            add_noise = add_noise,
            sd_noise = sd_noise,
            normalized = normalized
        )
    else
        make_S_matrix(
            data_train,
            data_val,
            n_features_base,
            n_features_inflections,
            ncol = ncol,
            sd_base_mean = sd_base_mean,
            sd_inflection_mean = sd_inflection_mean,
            sd_base = sd_base,
            sd_inflection = sd_inflection,
            isdeep = isdeep,
            add_noise = add_noise,
            sd_noise = sd_noise,
            normalized = normalized
        )
    end
end

function get_kwarg(kwargs, kw; required=false)
    if !haskey(kwargs, kw)
        if required
            throw(ArgumentError("$kw is not specified!"))
        else
            kwarg = get_default_kwargs(kw)
        end
    else
        kwarg = kwargs[kw]
    end

    # dump args
    KWARGS_DUMP[kw] = kwarg

    return kwarg
end

function get_default_kwargs(kw)
    if !haskey(KWARGS_DEFAULT, kw)
        throw(ArgumentError("$kw don't have default value!"))
    else
        return KWARGS_DEFAULT[kw]
    end
end


text:
export solve

"""
    default_operator(system::InitialValueProblem)

Return the default continous post operator for the initial value problem of a
discrete or continuous system.

### Input

- `system` -- an initial value problem wrapping a mathematical system (continuous or
              discrete) and a set of initial states

### Output

A continuous post operator with default options.
"""
function default_operator(system::InitialValueProblem)
    if isaffine(system)
        op = BFFPSV18()
    else
        error("no default reachability algorithm available for system of " *
              "type $(typeof(system))")
    end
    return op
end

"""
    solve(system, options)  or  solve(system, :key1 => val1, [...], keyK => valK)

Solves a reachability problem s.t. the given options.
If some options are not defined, we may fall back to default values.

### Input

- `system`    -- a (discrete or continuoues) system specification
- `options`   -- algorithm options for solving the problem
- `algorithm` -- (optional, default: dispatched on the system's type) the
                 reachability algorithm for the computation

### Output

A solution object whose content depends on the input options.

### Notes

To see all available input options, see
`keys(Reachability.available_keywords.dict)`.
"""
function solve(system::InitialValueProblem,
               options::Options;
               op::ContinuousPost=default_operator(system))
    solve!(system, copy(options), op=op)
end

solve(system::AbstractSystem, options::Pair{Symbol,<:Any}...) =
    solve(system, Options(Dict{Symbol,Any}(options)))

function solve!(problem::InitialValueProblem,
                options::Options;
                op::ContinuousPost=default_operator(problem)
               )

    @assert statedim(problem) == dim(problem.x0) "the state-space dimension should match the " *
    "dimension of the initial states, but they are of size $(statedim(problem)) and $(dim(problem.x0)) respectively"

    # normalize system to a canonical form if needed
    problem = IVP(normalize(problem.s), problem.x0)

    # initialize the algorithm-specific options
    options = init!(op, problem, options)

    # compute a coordinate transformation if needed
    problem, options = transform(problem, options)

    # run the continuous-post operator
    res = post(op, problem, options)

    # undo a coordinate transformation
    res = backtransform(res, options)
end

"""
    solve(system::InitialValueProblem{<:HybridSystem},
          options::Options)

Interface to reachability algorithms for a hybrid system PWA dynamics.

### Input

- `system`  -- hybrid system
- `options` -- options for solving the problem
"""
function solve(system::InitialValueProblem{<:HybridSystem, <:LazySet},
               options::Options,
               opC::ContinuousPost=BFFPSV18(),
               opD::DiscretePost=LazyDiscretePost())
    return solve!(distribute_initial_set(system), copy(options), opC, opD)
end

# if the initial states are distributed, no need to call distribute_initial_set(system)
function solve(system::InitialValueProblem{<:HybridSystem,
                                           <:Vector{<:Tuple{Int64,<:LazySet}}},
               options::Options,
               opC::ContinuousPost=BFFPSV18(),
               opD::DiscretePost=LazyDiscretePost())
    return solve!(system, copy(options), opC, opD)
end

function solve!(system::InitialValueProblem{<:HybridSystem,
                                            <:Vector{<:Tuple{Int64,<:LazySet{N}}}},
               options_input::Options,
               opC::ContinuousPost,
               opD::DiscretePost
              ) where N<:Real
    # update global variable
    global discrete_post_operator = opD

    HS = system.s
    init_sets = system.x0
    options = init!(opD, HS, options_input)
    time_horizon = options[:T]
    max_jumps = options[:max_jumps]
    property = options[:mode] == "check" ? options[:property] : nothing

    # waiting_list entries:
    # - (discrete) location
    # - (set of) continuous-time reach sets
    # - number of previous jumps
    waiting_list = Vector{Tuple{Int, AbstractReachSet{LazySet{N}}, Int}}()

    for (loc_id, x0) in init_sets
        loc = HS.modes[loc_id]
        source_invariant = loc.X

        if x0 ⊆ source_invariant
            loc_x0set = x0
        else
            loc_x0set = intersection(source_invariant, x0)
        end

        if !isempty(loc_x0set)
            aux = ReachSet{LazySet{N}}(loc_x0set, zero(N), zero(N))
            push!(waiting_list, (loc_id, aux, 0))
        end
    end

    # passed_list maps the (discrete) location to the (set of) continuous-time
    # reach sets
    passed_list = options[:fixpoint_check] == :none ?
        nothing :
        Vector{Vector{AbstractReachSet{LazySet{N}}}}(undef, nstates(HS))

    flowpipes = Vector{Flowpipe}()
    while (!isempty(waiting_list))
        loc_id, X0, jumps = pop!(waiting_list)
        loc = HS.modes[loc_id]
        source_invariant = loc.X

        # compute reach tube
        options_copy = copy(options)
        options_copy.dict[:T] = time_horizon - time_start(X0)
        options_copy.dict[:project_reachset] = false
        if property != nothing
            # TODO temporary hack, to be resolved in #447
            options_copy[:mode] = "reach"
        end

        if opC isa BFFPS19
            opC.options.specified[:HS] = HS
            opC.options.specified[:loc_id] = loc_id
            opC.options.specified[:opD] = opD
        end

        reach_tube = solve!(IVP(loc, set(X0)), options_copy, op=opC)
        reach_sets = reach_tube.flowpipes[1].reachsets

        # get the property for the current location
        property_loc = property isa Dict ?
                       get(property, loc_id, nothing) :
                       property

        # add the very first initial approximation
        if passed_list != nothing &&
                (!isassigned(passed_list, loc_id) || isempty(passed_list[loc_id]))
            reach_set = reach_sets[1]
            # TODO For lazy X0 the fixpoint check is likely to fail, so we
            # currently ignore that. In general, we want to add an
            # *underapproximation*, which we currently do not support.
            X = set(reach_set)
            if !(X isa CartesianProductArray) || !(array(X)[1] isa CH)
                Xoa = overapproximate(X)
                ti, tf = time_start(reach_set), time_end(reach_set)
                passed_list[loc_id] = [ReachSet{LazySet{N}}(Xoa, ti, tf)]
            end
        end

        Rsets = tube⋂inv(opD, reach_sets, loc.X, [time_start(X0), time_end(X0)])

        if property_loc != nothing
            if opD isa DecomposedDiscretePost
                temp_vars = opD.options[:temp_vars]
                n_lowdim = length(temp_vars)
                n = dim(set(X0))
                property_loc_lowdim = project(property_loc, temp_vars)
                for (i, reach_set) in enumerate(reach_sets)
                    X = set(reach_set)
                    if (dim(X) == n_lowdim && n_lowdim < n)
                        if !property_loc_lowdim(X)
                            return CheckSolution(false, i, options)
                        end
                    elseif !property_loc(X)
                        return CheckSolution(false, i, options)
                    end
                end
            else
                for (i, reach_set) in enumerate(Rsets)
                    if !property_loc(set(reach_set))
                        return CheckSolution(false, i, options)
                    end
                end
            end
        end

        if jumps == max_jumps
            continue
        end
        post(opD, HS, waiting_list, passed_list, loc_id, Rsets, jumps, options)
        # create vector with concrete set type (needed by ReachSolution)
        push!(flowpipes, Flowpipe([rs for rs in Rsets]))
    end
    if options[:mode] == "check"
        return CheckSolution(true, -1, options)
    end

    # Projection
    if options[:project_reachset] || options[:projection_matrix] != nothing
        info("Projection...")
        RsetsProj = @timing project(ReachSolution(flowpipes, options))
    else
        RsetsProj = ReachSolution(flowpipes, options)
    end
    return RsetsProj
end

ext: jl


text:
module MPolyFact

using Hecke
import Hecke: Nemo, @vprint, @hassert, @vtime, rational_reconstruction, set_precision!,
              set_precision, find_morphism

import Nemo: shift_left, shift_right
import Base: *

export factor_absolute

add_verbose_scope(:AbsFact)
add_assert_scope(:AbsFact)

function Hecke.norm(f::MPolyElem{nf_elem})
  Kx = parent(f)
  K = base_ring(Kx)
  n = nvars(Kx)
  Qx, x = PolynomialRing(QQ, [String(x) for x= symbols(Kx)], cached = false)
  Qxy, y = PolynomialRing(Qx, "y", cached = false)
  gg = [MPolyBuildCtx(Qx) for i=1:degree(K)]
  for (c, e) = zip(coefficients(f), exponent_vectors(f))
    for i=0:degree(K)-1
      d = coeff(c, i)
      if !iszero(d)
        push_term!(gg[i+1], d, e)
      end
    end
  end
  g = Qxy(map(finish, gg))
  return resultant(g, K.pol(y))
end

#= AbsFact: strategy
 f in Q[x,y] = Q[x][y], not necc. monic, irre. over Q

 find evaluation x <- x+ly s.th. g = f(0, y) is square-free
 find p s.th g mod p is square-free and such that the lcm of the factors
   is small
 compute the roots up to a bound B over the F_q splitting field
   as power series
 compute the coeffs. of f (as elements in Q[x]) as power series

 TODO: factor f over F_p[[t]], then compute the roots of the
       factors in F_q[[t]]. Each factor over F_p[[t]] should
       correspong to a Frobenius orbit of roots.
       This way the degrees are smaller....
       Done.

 TODO: refactor: a HenselCtx for fmpq_mpoly (or fmpz_mpoly) factored
       in F_p[[t]][x]
       a RootCtx for poly in F_p[[t]][x] to find roots in F_q[[t]]
       Done.

       See what Dan did in this case.
=#

function set_precision(f::PolyElem{T}, n::Int) where {T <: SeriesElem}
  g = parent(f)()
  for i=0:length(f)
    setcoeff!(g, i, set_precision(coeff(f, i), n))
  end
  return g
end

function set_precision!(f::PolyElem{T}, n::Int) where {T <: SeriesElem}
  for i=0:length(f)
    setcoeff!(f, i, set_precision!(coeff(f, i), n))
  end
  return f
end

mutable struct HenselCtxFqRelSeries{T}
  f :: fmpz_mpoly # bivariate
  n :: Int # number of factors
  lf :: Vector{PolyElem{T}} # T should be nmod_rel_series or fq_nmod_rel_series
  cf :: Vector{PolyElem{T}} # the cofactors for lifting
  t :: Int # shift, not used, so might be wrong.

  function HenselCtxFqRelSeries(f::fmpz_mpoly, lf::Vector{<:PolyElem{S}}, lg::Vector{<:PolyElem{S}}, n::Int, s::Int = 0) where {S <: Union{Nemo.FinFieldElem, Nemo.nmod}}
    @assert ngens(parent(f)) == 2
    k = base_ring(lf[1])
    R, t = PowerSeriesRing(k, 10, "t", cached = false) #, model = :capped_absolute)
    Rx, x = PolynomialRing(R, cached = false)
    r = new{typeof(t)}()
    r.f = f
    r.n = n
    r.t = s
    r.lf = [set_precision(s(x), 1) for s = lf]
    r.cf = [set_precision(s(x), 1) for s = lg]
    return r
  end

  function HenselCtxFqRelSeries(f::fmpz_mpoly, lf::Vector{<:PolyElem{<:SeriesElem{qadic}}}, lc::Vector{<:PolyElem{<:SeriesElem{qadic}}}, n::Int, s::Int = 0)
    @assert ngens(parent(f)) == 2
    r = new{elem_type(base_ring(lf[1]))}()
    r.f = f
    r.n = n
    r.t = s
    r.lf = lf
    r.cf = lc
    return r
  end

  function HenselCtxFqRelSeries(f::fmpz_mpoly, p::Int, s::Int = 0)
    k = GF(p, cached = false)
    return HenselCtxFqRelSeries(f, k, s)
  end

  function HenselCtxFqRelSeries(f::fmpz_mpoly, k::FinField, s::Int = 0)
    kt, t = PolynomialRing(k, cached = false)
    g = evaluate(f, [t, kt(-s)])
    issquarefree(g) || return nothing
    @assert issquarefree(g)

    lf = collect(keys(factor(g).fac))
    return HenselCtxFqRelSeries(f, lf, s)
  end

  function HenselCtxFqRelSeries(f::fmpz_mpoly, lf::Array{<:PolyElem{<:SeriesElem{<:FinFieldElem}}}, s::Int = 0)
    k, mk = ResidueField(base_ring(lf[1]))
    kt, t = PolynomialRing(k, cached = false)
    return HenselCtxFqRelSeries(f, [map_coefficients(mk, x, parent = kt) for x = lf], s)
  end

  function HenselCtxFqRelSeries(f::fmpz_mpoly, lf::Array{<:PolyElem{<:FinFieldElem}}, s::Int = 0)
    n = length(lf)
    lg = typeof(lf[1])[]
    i = 1
    while i < length(lf)
      g, a, b = gcdx(lf[i], lf[i+1])
      @assert isone(g)
      push!(lg, a)
      push!(lg, b)
      push!(lf, lf[i] * lf[i+1])
      i += 2
    end
    k = base_ring(lf[1])
    if isa(k, Nemo.GaloisField)
      p = Int(characteristic(k))
      k = quo(ZZ, p)[1]
      kt, t = PolynomialRing(k, cached = false)
      lf = [map_coefficients(k, x, parent = kt) for x = lf]
      lg = [map_coefficients(k, x, parent = kt) for x = lg]
    end

    return HenselCtxFqRelSeries(f, lf, lg, n, s)
  end
end

function Hecke.precision(H::HenselCtxFqRelSeries{<:Generic.RelSeries{qadic}})
  return precision(coeff(coeff(H.lf[1], 0), 0)), precision(coeff(H.lf[1], 0))
end

function shift_coeff_left!(f::PolyElem{<:SeriesElem}, n::Int)
  for i=0:length(f)
    setcoeff!(f, i, shift_left(coeff(f, i), n))
  end
end

function shift_coeff_right!(f::PolyElem{<:SeriesElem}, n::Int)
  for i=0:length(f)
    setcoeff!(f, i, shift_right(coeff(f, i), n))
  end
end

function shift_coeff_left(f::PolyElem{<:SeriesElem}, n::Int)
  g = parent(f)()
  for i=0:length(f)
    setcoeff!(g, i, shift_left(coeff(f, i), n))
  end
  return g
end

function shift_coeff_right(f::PolyElem{<:SeriesElem}, n::Int)
  g = parent(f)()
  for i=0:length(f)
    setcoeff!(g, i, shift_right(coeff(f, i), n))
  end
  return g
end

function lift(C::HenselCtxFqRelSeries{<:SeriesElem})
  St = parent(C.lf[1])
  S = base_ring(C.lf[1])

  C.lf[1]
  pr = precision(coeff(C.lf[1], 0))
  N2 = 2*pr

  S.prec_max = N2+1

  i = length(C.lf)
  if i == 1
    C.lf[1] = evaluate(C.f, [gen(St), St(gen(S)-C.t)])
    C.lf[1] *= inv(leading_coefficient(C.lf[1]))
    return
  end
  j = i-1
  while j > 0
    if i==length(C.lf)
      f = evaluate(C.f, [gen(St), St(gen(S)-C.t)])
      f *= inv(leading_coefficient(f))
    else
      @assert ismonic(C.lf[i])
      f = set_precision(C.lf[i], N2)
      @assert ismonic(C.lf[i])
    end
    @assert ismonic(f)

    #formulae and names from the Flint doc
    h = C.lf[j]
    g = C.lf[j-1]
    b = C.cf[j]
    a = C.cf[j-1]
    set_precision!(h, N2)
    set_precision!(g, N2)
    set_precision!(a, N2)
    set_precision!(b, N2)

    fgh = shift_coeff_right(f-g*h, pr)

    G = shift_coeff_left(rem(fgh*b, g), pr)+g
    H = shift_coeff_left(rem(fgh*a, h), pr)+h

    t = shift_coeff_right(1-a*G-b*H, pr)

    B = shift_coeff_left(rem(t*b, g), pr)+b
    A = shift_coeff_left(rem(t*a, h), pr)+a
    if i < length(C.lf)
      @assert ismonic(G)
      @assert ismonic(H)
      C.lf[i] = G*H
      @assert ismonic(C.lf[i])
    end
    @assert ismonic(G)
    @assert ismonic(H)
    C.lf[j-1] = G
    C.lf[j] = H
    C.cf[j-1] = A
    C.cf[j] = B
    i -= 1
    j -= 2
  end
end

function _set_precision(f::PolyElem{<:SeriesElem{qadic}}, n::Int)
  g = deepcopy(f)
  return _set_precision!(g, n)
end

function _set_precision!(f::PolyElem{<:SeriesElem{qadic}}, n::Int)
  for i=0:length(f)
    c = coeff(f, i)
    for j=0:pol_length(c)
      setprecision!(polcoeff(c, j), n)
    end
  end
  return f
end
# TODO: bad names...
function _shift_coeff_left(f::PolyElem{<:SeriesElem{qadic}}, n::Int)
  g = parent(f)()
  for i = 0:length(f)
    setcoeff!(g, i, map_coefficients(x -> shift_left(x, n), coeff(f, i), parent = base_ring(f)))
  end
  return g
end

function check_qadic(a::qadic)
  v = a.val
  a.val = 0
  f = Hecke.lift(Hecke.Globals.Zx, a)
  p = prime(parent(a))^a.N
  a.val = v
  @assert all(x->abs(x) < p, coefficients(f))
end

function shift_right(a::qadic, n::Int)
  b = deepcopy(a)
  b.val -= n
  return b
end

function shift_left(a::qadic, n::Int)
  b = deepcopy(a)
  b.val += n
  return b
end

function _shift_coeff_right(f::PolyElem{<:SeriesElem{qadic}}, n::Int)
  g = parent(f)()
  for i = 0:length(f)
    @assert all(y -> valuation(polcoeff(coeff(f, i), y)) >= n, 0:pol_length(coeff(f, i)))
    setcoeff!(g, i, map_coefficients(x -> shift_right(x, n), coeff(f, i), parent = base_ring(f)))
  end
  return g
end

mutable struct Preinv{T, S <: PolyElem{T}}
  f::S
  n::Int
  fi::S
  function Preinv(f::PolyElem)
    r = new{elem_type(base_ring(f)), typeof(f)}()
    r.f = reverse(f)
    @assert degree(f) == degree(r.f)
    r.fi = one(parent(f))
    r.n = 1
    return r
  end
end

preinv(f::PolyElem) = Preinv(f)

function lift(P::Preinv)
  f = truncate(P.f, 2*P.n)
  P.fi = truncate(P.fi*truncate((parent(f)(2)-f*P.fi), 2*P.n), 2*P.n)
  P.n *= 2
end
# von zur Gathen: Modern Computer Algebra, p 243:
#  9.1. Division with remainder using Newton iteration
function Base.rem(g::PolyElem, P::Preinv)
  if degree(g) < degree(P.f)
    return g
  end
  if degree(g) == degree(P.f)
    return g - leading_coefficient(g)*reverse(P.f)
  end

  n = degree(g) # a in the book
  m = degree(P.f) # b in the book
  @assert m <= n

  gr = reverse(g)

  while P.n < n-m+1
    lift(P)
  end
  q = truncate(gr*P.fi, n-m+1)
  q = reverse(q, n-m+1)
  r = g-q*reverse(P.f)
  r = truncate(r, m)
  @hassert :AbsFact 2 r == rem(g, reverse(P.f))
  return r
end

function check_data(f::PolyElem{<:SeriesElem{qadic}})
  for c = coefficients(f)
    for i=1:pol_length(c)
      check_qadic(polcoeff(c, i))
    end
  end
end

function lift_q(C::HenselCtxFqRelSeries{<:SeriesElem{qadic}})
  St = parent(C.lf[1])
  S = base_ring(C.lf[1])
  Q = base_ring(S)

  pr = precision(coeff(coeff(C.lf[1], 0), 0))
  N2 = 2*pr

  setprecision!(Q, N2)

  i = length(C.lf)
  @assert i > 1
  @assert all(ismonic, C.lf[1:end-1])
  j = i-1
  while j > 0
    if i==length(C.lf)
      f = evaluate(map_coefficients(Q, C.f), [gen(St), St(gen(S))])
      f *= inv(leading_coefficient(f))
    else
#      f = _set_precision(C.lf[i], N2)
      f = C.lf[i]
      @assert precision(coeff(coeff(f, 0), 0)) >= N2
      @assert ismonic(C.lf[i])
    end
    @assert ismonic(f)
    #formulae and names from the Flint doc
    h = C.lf[j]
    g = C.lf[j-1]
    b = C.cf[j]
    a = C.cf[j-1]
    h = _set_precision(h, N2)
    g = _set_precision(g, N2)
    a = _set_precision(a, N2)
    b = _set_precision(b, N2)

    fgh = _shift_coeff_right(f-g*h, pr)

    @assert ismonic(g)
    @assert !iszero(constant_coefficient(g))
    @assert ismonic(h)
    @assert !iszero(constant_coefficient(h))
    gi = preinv(g)
    hi = preinv(h)

    G = _shift_coeff_left(rem(fgh*b, gi), pr)+g
    @assert ismonic(G)
    H = _shift_coeff_left(rem(fgh*a, hi), pr)+h
    @assert ismonic(H)

    t = _shift_coeff_right(1-a*G-b*H, pr)

    B = _shift_coeff_left(rem(t*b, gi), pr)+b
    A = _shift_coeff_left(rem(t*a, hi), pr)+a
#    check_data(A)
#    check_data(B)

    C.lf[j-1] = G
    C.lf[j] = H
    C.cf[j-1] = A
    C.cf[j] = B
    i -= 1
    j -= 2
  end
end

mutable struct RootCtxSingle{T}
  f::PolyElem{T}
  R::T  # the root
  o::T  # inv(f'(R)) for the double lifting.

  function RootCtxSingle(f::PolyElem{S}, K::FqNmodFiniteField) where {S <: SeriesElem}
    #not used I think
    RR,  = PowerSeriesRing(K, max_precision(R), string(var(R)), cached = false) #can't get the modell
    return RootCtxSingle(f, RR)
  end

  function RootCtxSingle(f::PolyElem{<:SeriesElem{T}}, r::T) where {T}
    R = base_ring(parent(f))
    k, mk = ResidueField(R)
    g = map_coefficients(mk, f)
    # should be zero-ish, but if T is acb, this is difficult.
    isexact_type(T) && @assert iszero(g(r))
    o = inv(derivative(g)(r))
    return new{elem_type(R)}(f, R([r], 1, 1, 0), R([o], 1, 1, 0))
  end

  function RootCtxSingle(f::PolyElem{S}, RR::FqNmodRelSeriesRing) where {S <: SeriesElem}
    K = base_ring(RR)
    R = base_ring(f) # should be a series ring
    r = new{elem_type(RR)}()
    r.f = map_coefficients(x->map_coefficients(K, x, parent = RR), f)
    k, mk = ResidueField(R)
    _, mK = ResidueField(RR)
    g = map_coefficients(mk, f)
    @vtime :AbsFact 2 rt = Nemo.any_root(map_coefficients(K, g))
    r.R = preimage(mK, rt)
    g = map_coefficients(K, g)
    @vtime :AbsFact 2 r.o = preimage(mK, inv(derivative(g)(r.R)))
    return r
  end
end

"""
Computes the roots of `f` as power series over the complex numbers up to
power series precision `pr` and complex precision `prec`. The roots
are given as power series around `r`, the 2nd argument.

Not very bright algorithm, `f` has to be square-free and `r` cannot be a singularity.
"""
function analytic_roots(f::fmpz_mpoly, r::Integer, pr::Int = 10; prec::Int = 100, max_roots = degree(f, 2))
  return analytic_roots(f, fmpz(r), pr, prec = prec, max_roots = max_roots)
end

function analytic_roots(f::fmpz_mpoly, r::fmpz, pr::Int = 10; prec::Int = 100, max_roots = degree(f, 2))
  @assert ngens(parent(f)) == 2
  g = evaluate(f, [Hecke.Globals.Zx(r), gen(Hecke.Globals.Zx)])
  @assert issquarefree(g)
  C = ComplexField(prec)
  rt = Hecke.roots(g, C)[1:max_roots]
  @assert all(x->parent(x) == C, rt)
  Cs, s = PowerSeriesRing(C, pr+2, "s", cached = false)
  Cst, t = PolynomialRing(Cs, cached = false)
  ff = evaluate(f, [Cst(s+C(r)), t])
  RT = []
  for x = rt
    R = RootCtxSingle(ff, x)
    for i=1:clog(pr, 2)
      lift(R)
    end
    push!(RT, R.R)
  end
  return RT
end

"""
Computes the roots of `f` as power series over some number field up to
power series precision `pr`. The roots
are given as power series around `r`, the 2nd argument.

Not very bright algorithm, `f` has to be square-free and `r` cannot be a singularity.
"""
function symbolic_roots(f::fmpz_mpoly, r::Integer, pr::Int = 10; max_roots::Int = degree(f, 2))
  return symbolic_roots(f, fmpz(r), pr; max_roots...)
end

function symbolic_roots(f::fmpz_mpoly, r::fmpz, pr::Int = 10; max_roots::Int = degree(f, 2))
  @assert ngens(parent(f)) == 2
  g = evaluate(f, [Hecke.Globals.Zx(r), gen(Hecke.Globals.Zx)])
  @assert issquarefree(g)
  lg = factor(g)
  rt = vcat([Hecke.roots(x, number_field(x)[1]) for x = keys(lg.fac)]...)
  rt = rt[1:min(length(rt), max_roots)]
  RT = []
  for i = 1:length(rt)
    x = rt[i]
    C = parent(x)
    if i>1 && C == parent(rt[i-1])
      Cs = parent(RT[i-1])
      s = gen(Cs)
    else
      Cs, s = PowerSeriesRing(C, pr+2, "s", cached = false)
    end
    Cst, t = PolynomialRing(Cs, cached = false)
    ff = evaluate(f, [Cst(s+C(r)), t])
    R = RootCtxSingle(ff, x)
    for i=1:clog(pr, 2)
      lift(R)
    end
    push!(RT, R.R)
  end
  return RT
end

"""
A single lifting step for the easy Newton iteration for a isolated simple root.
"""
function lift(R::RootCtxSingle)
  pr = precision(R.R)
  set_precision!(R.o, 2*pr)
  set_precision!(R.R, 2*pr)
  parent(R.R).prec_max = 2*pr+1
  R.o = R.o*(2-R.o*derivative(R.f)(R.R))
  R.R = R.R - R.f(R.R)*R.o
end

mutable struct RootCtx
  f :: fmpq_mpoly
  H :: HenselCtxFqRelSeries{nmod_rel_series}
  R :: Vector{RootCtxSingle{fq_nmod_rel_series}}
  RP :: Vector{Vector{fq_nmod_rel_series}}
  all_R :: Vector{fq_nmod_rel_series}

  function RootCtx(f::fmpq_mpoly, p::Int, d::Int, t::Int = 0)
    r = new()
    r.f = f
    den = lcm(map(denominator, coefficients(f)))
    g = map_coefficients(numerator, den*f)
    @vtime :AbsFact 2 mu = HenselCtxFqRelSeries(g, p, t)
    mu === nothing && return mu
    r.H = mu
    r.R = RootCtxSingle{fq_nmod_rel_series}[]
    K = GF(p, d)[1]
    S, _ = PowerSeriesRing(K, 10, "s", cached = false)
    for i=1:r.H.n
      @vtime :AbsFact 2 push!(r.R, RootCtxSingle(r.H.lf[i], S))
    end
    return r
  end
end

"""
Computes `R[i]^j`, cached
"""
function root(R::RootCtx, i::Int, j::Int)
  if length(R.all_R) == 0 || precision(R.R[1].R) > precision(R.all_R[1])
    empty!(R.all_R)
    for i=1:R.H.n
      push!(R.all_R, copy(R.R[i].R))
      @hassert :AbsFact 2 iszero(R.R[i].f(R.all_R[end]))
      S = parent(R.all_R[end])
      for j=1:degree(R.R[i].f)-1
        push!(R.all_R, map_coefficients(frobenius, R.all_R[end], parent = S))
        @hassert :AbsFact 3 iszero(R.R[i].f(R.all_R[end]))
      end
    end
  end

  if length(R.RP) == 0 || precision(R.all_R[1]) != precision(R.RP[1][1])
    o = one(parent(R.all_R[1]))
    R.RP = [[set_precision(o, precision(R.all_R[1])) for x = R.all_R], copy(R.all_R)]
  end
  if length(R.RP) > j
    return (R.RP[j+1][i])
  end
  while length(R.RP) <= j+1
    push!(R.RP, R.RP[2] .* R.RP[end])
  end

  s = R.RP[j+1][i]
  return (s)
end

"""
Computes the roots of `f` in a suitable splitting field: a power series
over a finite field. The characteristic will be chosen, trying to have the
degree small. The initial precision if `2`.

Returns, not the roots, but the root context `RootCtx`, so the precision
can be increased (`newton_lift`, `roots`, `more_precision`).
"""
function roots(f::fmpq_mpoly, p_max::Int=2^15; pr::Int = 2)
  @assert nvars(parent(f)) == 2
  #requires f to be irred. over Q - which is not tested
  #requires f(x, 0) to be same degree and irred. - should be arranged by the absolute_bivariate_factorisation

  #f in Qxy
  Zx = Hecke.Globals.Zx
  f *= lcm([denominator(x) for x = coefficients(f)])
  ff = map_coefficients(ZZ, f)
  #TODO: 0 might not be a good evaluation point...
  #f needs to be irreducible over Q and g square-free
  g = evaluate(ff, [gen(Zx), Zx(0)])
  @assert degree(g) == degree(f, 1)

  d = typemax(Int)
  best_p = p_max
  pc = 0
  local p
  while true
    p_max = next_prime(p_max)
    gp = factor(g, GF(p_max, cached = false))
    if any(x->x>1, values(gp.fac))
      @vprint :AbsFact 1 "not squarefree mod $p_max\n"
      continue
    end
    e = lcm([degree(x) for x = keys(gp.fac)])
    if e < d
      d = e
      best_p = p_max
      pc = 0
    else
      pc += 1
    end

    if e == 1 || pc > 1.5 * degree(g)
      p = best_p
      @vprint :AbsFact 1 "using $best_p of degree $d\n"
      break
    end
  end

  R = RootCtx(f, p, d)
  R === nothing && return R, p_max

  @vprint :AbsFact 1 "need to seriously lift $(R.H.n) elements\n"
  R.all_R = typeof(R.R[1].R)[]
  R.RP = typeof(R.R[1].R)[]

  for i=1:pr
    more_precision(R)
  end

  return R, p_max
end

function more_precision(R::RootCtx)
  lift(R.H)
  S = base_ring(R.R[1].f)
  S.prec_max = 2*S.prec_max + 1
  T = parent(R.R[1].f)
  K = base_ring(S)
  for i=1:R.H.n
    R.R[i].f = map_coefficients(x->map_coefficients(K, x, parent = S), R.H.lf[i], parent = T)
    lift(R.R[i])
  end
end

"""
Internal use.
"""
function combination(RC::RootCtx)
  #R is a list of roots, ie. power series over F_q (finite field)
  root(RC, 1, 1)
  f = RC.f
  R = RC.all_R
  Ft = parent(R[1])
  t = gen(Ft)
  n = precision(R[1])
  @assert all(x->precision(x) == n, R)

  #ps = [[div(x^i % tn, td^i) for i = 1:n] for x = R]

  F = base_ring(Ft)
  k = degree(F)

  p = characteristic(F)
  Fp = GF(p, cached = false)

  #the paper
  # https://www.math.univ-toulouse.fr/~cheze/facto_abs.m
  #seems to say that a combination works iff the sum of the
  # deg 2 terms vanish (for almost all choices of specialisations)
  #this would be even easier
  #also: it would make the lifting above trivial (precision 2)
  #deal: if monic in x of degree n and deg_y(coeff of x^(n-1)) = r,
  # then for the 1st power sum we get r as a degree bound. In the paper
  # r == 1...
  # reasoning: sum of second highest terms of factor == second highest
  # of poly => bound
  # similar would be for other power sums - if I'd compute them.
  #all vanishing results in the paper do not depend on the assumption
  #D. Rupprecht / Journal of Symbolic Computation 37 (2004) 557–574
  # doi:10.1016/S0747-7171(02)00011-1
  # https://core.ac.uk/download/pdf/82425943.pdf
  #
  # This allows (should allow) to only use the 1st power sum (trace)
  # which makes denominators easier
  #
  #for non-monic: think about interaction with any_order
  #
  # instead of LLL use direct rref over GF(q)
  # for non-monics, the generalised equation order might be better than
  # trying to scale
  #
  j = 0
  nn = matrix(Fp, 0, length(R), [])

  d = degree(f, 2)+degree(f, 1)
  lc = leading_coefficient(f, 1)
  d += degree(lc, 2)

  ld = evaluate(map_coefficients(x->F(ZZ(x)), lc), [set_precision(Ft(0), n), set_precision(gen(Ft), n)])
  @assert precision(ld) >= n
  R = R .* ld

  pow = 1
  bad = 0
  stable = 0
  last_rank = length(R)
  while true
    j += 1
    while pow*d+j >= n
      @vprint :AbsFact 1 "need more precicsion: $n ($d, $pow, $j)\n"
      @vtime :AbsFact 2 more_precision(RC)
      n = precision(RC.R[1].R)
      if false && n > 170 #too small - but a safety valve
        error("too much n")
      end
    end
    root(RC, 1, 1)
    R = RC.all_R
    n = precision(R[1])
    ld = evaluate(map_coefficients(x->F(ZZ(x)), lc), [set_precision(Ft(0), n), set_precision(gen(Ft), n)])
    R = R .* ld
    @assert precision(R[1]) >= n

    mn = transpose(matrix([[Fp(coeff(coeff(x^pow, pow*d+j), lk)) for lk = 0:k-1] for x = R]))

    if false && iszero(mn)
      @vprint :AbsFact 2 "found zero column, disgarding\n"
      bad += 1
      if bad > max(2, div(length(R), 2))
        pow += 1
        d += degree(lc, 2)
        @vprint :AbsFact 1 "increasing power to $pow\n"
        j = 0
        bad = 0
      end
      continue
    else
      @vprint :AbsFact 2 "found non zero column\n"
    end

    nn = vcat(nn, mn)

    ke = kernel(nn)
    @vprint :AbsFact 2 "current kernel dimension: $(ke[1])\n"
    if last_rank == ke[1]
      bad += 1
      if bad > max(2, div(length(R), 2))
        pow += 1
        @vprint :AbsFact 2 "increasing power to $pow\n"
        j = 0
        bad = 0
        continue
      end
    else
      bad = 0
      stable = 0
      last_rank = ke[1]
    end
    if ke[1] == 0 || mod(length(R), ke[1]) != 0 || mod(total_degree(f), ke[1]) != 0
      continue
    end
    m = ke[2]
    z = transpose(m)*m
    if z != div(length(R), ke[1])
      @vprint :AbsFact 2 "not a equal size partition\n"
      continue
    end
    stable += 1
    if stable < max(5, degree(f, 1)/3)
      @vprint :AbsFact 2 "need confirmation...\n"
      continue
    end
    return transpose(m)
  end
end

function lift(R::PolyRing{S}, s::SeriesElem{S}) where {S}
  t = R()
  for x = 0:pol_length(s)
    setcoeff!(t, x, polcoeff(s, x))
  end
  return shift_left(t, valuation(s))
end

#TODO: possibly rename and export. This is used elsewhere
# Combine with the Hecke._meet into a POSet structure?
"""
Compute an array of arrays of indices (Int) indicating a partitioning of
the indices according to the values in `a`
"""
function block_system(a::Vector{T}) where {T}
  d = Dict{T, Vector{Int}}()
  for i=1:length(a)
    if haskey(d, a[i])
      push!(d[a[i]], i)
    else
      d[a[i]] = [i]
    end
  end
  return sort(collect(values(d)), lt = (a,b) -> isless(a[1], b[1]))
end

#= bounds:
  f in Z[x,y], g, h in C[x,y]
  H(f) = max abs value coeff of f
  gh = f
  then
  H(g) H(h) <= 2^(deg_x(f) + deg_y(g) - 2) ((deg_x(f)+1)(deg_y(f)+1))^(1/2) H(f)
=#
"""
Internal use.
"""
function field(RC::RootCtx, m::MatElem)
  R = RC.all_R
  P = RC.f

  bnd = numerator(maximum(abs(x) for x = coefficients(RC.f))) * fmpz(2)^(degree(RC.f, 1) + degree(RC.f, 2)-2) * Hecke.root(fmpz((degree(RC.f, 1)+1)*(degree(RC.f, 2)+1)), 2)
  #all coeffs should be bounded by bnd...

  #we have roots, we need to combine roots for each row in m where the entry is pm 1
  #the coeffs then live is a number field, meaning that the elem sym functions or
  #the power sums will be needed
  #the field degree apears to be nrows(m)...

  #need primitive element, can use power sums up to #factors

  #we will ONLY find one factor, the others are Galois conjugate
  #the field here is not necc. normal

  #the leading_coeff of P needs to be monic.

  d_f = div(ncols(m), nrows(m))
  @assert ncols(m) == length(R)

  Ft = parent(R[1])
  F = base_ring(Ft)
  t = gen(Ft)
  d = precision(R[1])

  #TODO think, the bound might be too large...
  tf = divexact(total_degree(P), nrows(m)) + degree(P, 2)

  #TODO use Frobenius (again) to save on multiplications
  #TODO invest work in one factor only - need only powers of the roots involved there
  #     the other factor is then just a division away
  #     if complete orbits are combined, use the trace (pointwise) rather than powers
  @vprint :AbsFact 2 "combining: $([findall(x->!iszero(x), collect(m[i, :])) for i=1:nrows(m)])\n"
  k, mk = ResidueField(parent(R[1]))
  kt, t = PolynomialRing(k, cached = false)
  RR = map(mk, R)
  RP = [copy(RR)]
  for j=2:d_f
    push!(RP, RR .* RP[end])
  end
  @vtime :AbsFact 2 el = [[sum(RP[j][i] for i=1:ncols(m) if m[lj, i] != 0) for j=1:d_f] for lj=1:nrows(m)]

  #now find the degree where the coeffs actually live:
  k = 1
  for x = el
    for y = x
      d = degree(minpoly(y))
      k = lcm(d, k)
    end
  end

  @vprint :AbsFact 1 "target field has (local) degree $k\n"

  Qq = QadicField(characteristic(F), k, 1, cached = false)[1]
  Qqt = PolynomialRing(Qq, cached = false)[1]
  k, mk = ResidueField(Qq)

  phi = find_morphism(k, F) #avoids embed - which stores the info

  kt, t = PolynomialRing(k, cached = false)

  fl = [power_sums_to_polynomial(map(t->preimage(phi, t), x)) for x = el]
  fl = [map_coefficients(x->x, y, parent = kt) for y = fl]
  HH = HenselCtxFqRelSeries(RC.H.f, fl)
  while precision(coeff(HH.lf[1], 0)) < tf+2
    lift(HH)
  end

  kXY, (X, Y) = PolynomialRing(k, ["X", "Y"], cached = false)

  nl = []
  kS = PowerSeriesRing(k, tf+2, "s")[1]
  for x = HH.lf[1:HH.n]
    f = MPolyBuildCtx(kXY)
    lc = one(kt)
    cz = []
    z = x
    for i=0:degree(x)
      c = coeff(z, i)
      local fl, n, d = rational_reconstruction(c, parent = kt)
      @assert fl
      b = lcm(lc, d)
      n *= divexact(b, d)
      if b != lc
        cz = divexact(b, lc) .* cz
        lc = b
      end
      push!(cz, n)
    end
    for i=0:degree(z)
      c = cz[i+1]
      for j=0:degree(c)
        if !iszero(coeff(c, j))
          push_term!(f, coeff(c, j), [i,j])
        end
      end
    end
    push!(nl, finish(f))
  end
  @vprint :AbsFact 2 "now building the leading coeff...\n"
  lc = map(x->evaluate(leading_coefficient(x, 1), [0*t, t]), nl)

  for i = 1:length(lc)
    l = leading_coefficient(lc[i])
    lc[i] *= inv(l)
    nl[i] *= inv(l)
  end

  llc = coeff(P, 1)
  P *= inv(llc)

  _lc = evaluate(leading_coefficient(P, 1), [zero(Hecke.Globals.Qx), gen(Hecke.Globals.Qx)])
  _lc = Hecke.squarefree_part(_lc)
  local H, fa
  if !isone(_lc)
    ld = coprime_base(vcat(lc, [map_coefficients(k, _lc, parent = kt)]))
    if sum(map(degree, ld)) != degree(_lc) #TODO: this should(?) be caught earlier
      @vprint :AbsFact 1 "leading coeff not square-free mod p, bad prime\n"
      return nothing
    end

    fa = [[valuation(x, y) for y = ld] for x = lc]
    lc = _lc
    H = Hecke.HenselCtxQadic(map_coefficients(Qq, lc, parent = Qqt), ld)
  else
    @vprint :AbsFact 2 "is monic, no leading coefficient...\n"
    lc = _lc
    fa = []
  end

  el = nl

  @vprint :AbsFact 1 "locating primitive element...\n"
  #currently, we're in F_q[[t]], not Q_q[[t]], however, the primitivity
  #can already be decided over F_q as lifting can not enlarge the field

  #if no single coefficient is primtive, use block systems and sums of coeffs
  #to find a primitive one.
  all_bs = []
  pe_j = 0
  for i = 1:length(el[1])
    bs = block_system([coeff(x, i) for x = el])
    @vprint :AbsFact 3 "block system of coeff $i is $bs\n"
    @assert all(x->length(x) == length(bs[1]), bs)
    if length(bs[1]) == 1
      pe_j = i
      break
    end
    push!(all_bs, bs)
  end

  local pe, pow, used
  if pe_j == 0
    @vprint :AbsFact 2 "no single coefficient is primitive, having to to combinations\n"
    bs = [collect(1:length(el))]
    used = Int[]
    for i=1:length(el[1])
      cs = Hecke._meet(bs, all_bs[i])
      if length(cs) > length(bs)
        bs = cs
        push!(used, i)
      end
      if length(bs[1]) == 1
        break
      end
    end
    @vprint :AbsFact 2 "using coeffs $used to form primtive element\n"
    pow = 1
    while true
      pe = x -> (sum(coeff(x, t) for t = used)^pow)
      bs = block_system([pe(x) for x = el])
      if length(bs[1]) == 1
        @vprint :AbsFact 2 "using sum to the power $pow\n"
        break
      end
      pow += 1
    end
  else
    @vprint :AbsFact 2 "$(pe_j)-th coeff is primitive\n"
    pe = x -> coeff(x, pe_j)
    pow = 1
    used = [1]
  end

  @vprint :AbsFact 1 "hopefully $(length(el)) degree field\n"

  bnd = (length(el)*(length(used)*bnd))^(pow*length(el))
  @vprint :AbsFact 1 "power sums (coeffs of minpoly of field) should have at most $(clog(bnd, characteristic(F))) p-adic digits\n"

  #TODO: Think: Do we need to lift all n factors? In the end we're going
  #      to return el[1] and prod(el[2:n]) only.
  #      Problem: currently I need all conjugates of the coeffs, hence all
  #      the el's individually.

  SQq, _ = PowerSeriesRing(Qq, tf+2, "s", cached = false)
  SQqt, _ = PolynomialRing(SQq, cached = false)

  mc(f) = # PolyElem{SeriesElem{Fq}} -> PolyElem{SeriesElem{Qq}}
    map_coefficients(x->map_coefficients(y->setprecision(preimage(mk, y), 1), x, parent = SQq), f, parent = SQqt)

  HQ = HenselCtxFqRelSeries(HH.f, map(mc, HH.lf), map(mc, HH.cf), HH.n)

  QqXY, (X, Y) = PolynomialRing(Qq, 2, cached = false)

  pr = 1
  while true
    pr *= 2
    if pr > 800
      error("too bas")
    end
    @vprint :AbsFact 1  "using p-adic precision of $pr\n"

    setprecision!(Qq, pr+1)
    if length(fa) > 0
      H.f = map_coefficients(Qq, _lc, parent = Qqt)
      @vprint :AbsFact 2 "lifting leading coeff factorisation\n"
      @vtime :AbsFact 2 Hecke.lift(H, pr+1)
      fH = factor(H)
      lc = [prod(fH[i]^t[i] for i=1:length(t)) for t = fa]
    end

    @vprint :AbsFact 1 "lifting factors\n"
    @vtime :AbsFact 2 while precision(coeff(coeff(HQ.lf[1], 0), 0)) < pr+1
      lift_q(HQ)
    end

    if length(fa) > 0
      z = [lc[i](gen(SQq)) * HQ.lf[i] for i=1:HQ.n]
    else
      z = HQ.lf[1:HQ.n]
    end

    setprecision!(coeff(X, 1), pr+2)
    setprecision!(coeff(Y, 1), pr+2)
    el = [map_coefficients(q -> lift(Qqt, q)(Y), f)(X) for f = z]

#    # lift mod p^1 -> p^pr x^2+y^2+px+1 was bad I think
#    @vtime :AbsFact 1 ok, el = lift_prime_power(P*inv(coeff(P, 1)), el, [0], 1, pr)
#    ok || @vprint :AbsFact 1 "bad prime found, q-adic lifting failed\n"
#    ok || return nothing
#    @assert ok  # can fail but should fail for only finitely many p

    #to make things integral...
    fl = Qq(llc) .* el

    p = [coeff(sum(pe(x)^l for x = fl), 0) for l=1:length(el)]
    p = map(rational_reconstruction, p)

    if !all(x->x[1], p)
      @vprint :AbsFact 2 "reco failed (for poly), increasing p-adic precision\n"
      if 2*clog(bnd, prime(Qq)) < pr
        @vprint :AbsFact 2 "bad prime? too much precision and still no poly, so changing prime\n"
        return nothing
      end
      continue
    end

    p = [x[2]//x[3] for x = p]
    p = Hecke.power_sums_to_polynomial(p)
    if any(x->denominator(x)>1, coefficients(p))
      @vprint :AbsFact 2 "poly wrong, increasing p-adic precision\n"
      continue
    end

    k, a = number_field(p)

    @vprint :AbsFact 1  "using as number field: $k\n"

    m = transpose(matrix([[pe(x)^l for x = fl] for l=0:degree(k)-1]))
    kx, x = PolynomialRing(k, "x", cached = false)
    kX, _ = PolynomialRing(k, ["X", "Y"], cached = false)
    B = MPolyBuildCtx(kX)
    for j=1:length(el[1])
      n = transpose(matrix([[coeff(x, j)] for x = fl]))
      s = solve(m, transpose(n))
      @assert all(x->iszero(coeff(s[x, 1], 1)), 1:degree(k))
      s = [rational_reconstruction(coeff(s[i, 1], 0)) for i=1:degree(k)]
      if !all(x->x[1], s)
        break
      end
      push_term!(B, k([x[2]//x[3] for x = s]), exponent_vector(el[1], j))
    end
    q = finish(B)
    if length(q) < length(el[1])
      continue
    end
    b, r = divrem(map_coefficients(k, P, parent = kX), [q])
    if iszero(r)
      return q, b[1]
    end
    @vprint :AbsFact 2 "division failed, increasing precision\n"
  end
end

"""
Given an irreducible bivariate polynomial over `Q` compute the
absolute factorisation.

Returns two polynomials:
 - the (absolutely) irreducible factor over the smallest number field
 - the co-factor

All the other factors would be conjugates of the first, hence representing
them exactly requires the splitting field to be computed.
"""
function absolute_bivariate_factorisation(f::fmpq_mpoly)
  d = degree(f, 1)
  R = parent(f)
  x, y = gens(R)

  lf = factor(f)
  if length(lf.fac) > 1 || any(x->x>1, values(lf.fac))
    error("poly must be irreducible over Q")
  end

  if degree(f, 2) < d
    @vprint :AbsFact 1 "swapping variables to have smaller degree\n"
    f = evaluate(f, [y, x])
    a, ca = absolute_bivariate_factorisation(f)
    S = parent(a)
    X, Y = gens(S)
    return evaluate(a, [Y, X]), evaluate(ca, [Y, X])
  end

  if degree(f, 2) == d && !isone(leading_coefficient(f, 1)) && isone(leading_coefficient(f, 2))
    @vprint :AbsFact 1 "swapping variables to be monic\n"
    f = evaluate(f, [y, x])
    a, ca = absolute_bivariate_factorisation(f)
    S = parent(a)
    X, Y = gens(S)
    return evaluate(a, [Y, X]), evaluate(ca, [Y, X])
  end

  Qt, t = PolynomialRing(QQ, cached = false)

  if degree(f, 1) == 0 #constant in
    k, a = number_field(evaluate(f, [Qt(0), t]), cached = false)
    kXY, (X, Y) = PolynomialRing(k, ["X", "Y"], cached = false)
    b = Y-a
    return b, divexact(map_coefficients(k, f, parent = kXY), b)
  end

  s =-1
  while true
    s += 1
    @vprint :AbsFact 1 "substitution to $s\n"
    z = evaluate(f, [t, Qt(s)])
    if degree(z) == d && issquarefree(z)
      break
    end
  end
  ff = evaluate(f, [x, y+s])
  d = lcm([denominator(x) for x = coefficients(ff)])
  ff *= d
  @assert all(x->isone(denominator(x)), coefficients(ff))
  gg = ff

  p = 2^24
  local aa
  while true
    @vtime :AbsFact 1 r, p = roots(gg, p)
    if r === nothing
      continue
    end
    @vtime :AbsFact 1 z = combination(r)
    if nrows(z) == 1
      return f, one(parent(f))
    end

    @vtime :AbsFact 1 aa = field(r, z)
    aa !== nothing && break
  end
  a, b = aa
  S = parent(a)
  X, Y = gens(S)
  a = evaluate(a, [X, Y-s])
  b = evaluate(b, [X, Y-s])
  return a, b
end

# From Dan, the lifting in Qq[x,y]

function map_down(Rp, a, mKp :: Map, pr::Int)
    M = MPolyBuildCtx(Rp)
    pk = prime(base_ring(a))^pr
    for (c, v) in zip(coefficients(a), exponent_vectors(a))
        @assert valuation(c) >= pr
        q = divexact(c, pk)  #should be a shift
        push_term!(M, mKp(q), v)
    end
    return finish(M)
end

function map_up(R, a, mKp :: Map, pr::Int)
    Rp = parent(a)
    Kp = base_ring(Rp)
    M = MPolyBuildCtx(R)
    pk = prime(base_ring(R))^pr

    for (c, v) in zip(coefficients(a), exponent_vectors(a))
        d = preimage(mKp, c)
        push_term!(M, d*pk, v)
    end
    return finish(M)
end

#=
    supposed to have a = prod(fac) mod p^kstart
    furthermore, the fac's are supposed to be pairwise coprime univariate in
        Fq[gen(1)] when evaluated at gen(2) = alphas[1], ... gen(n) = alphas[n-1]

    try to lift to a factorization mod p^kstop  (or maybe its mod p^(kstop+1))
=#
function lift_prime_power(
    a::fmpq_mpoly,
    fac::Vector{Generic.MPoly{qadic}},
    alphas::Vector,
    kstart::Int,
    kstop::Int)

    if kstop <= kstart
        return
    end

    r = length(fac)
    R = parent(fac[1])
    n = nvars(R)
    ZZ = base_ring(R)
    Kp, mKp = ResidueField(ZZ)
    Rp, x = PolynomialRing(Kp, n, cached = false)

    minorvars = [i for i in 2:n]
    degs = [degree(a, i) for i in 2:n]

    md = [map_down(Rp, f, mKp, 0) for f in fac]
    ok, I = Hecke.AbstractAlgebra.MPolyFactor.pfracinit(md, 1, minorvars, alphas)
    @assert ok  # evaluation of fac's should be pairwise coprime

    a = map_coefficients(ZZ, a, parent = parent(fac[1]))

    for l in kstart:kstop
        error = a - prod(fac)

        for c in coefficients(error)
          if valuation(c) < l
            throw(AssertionError("factorization is not correct mod p^$l"))
          end
        end

        if iszero(error)
            break
        end
        t = map_down(Rp, error, mKp, l)
        ok, deltas = Hecke.AbstractAlgebra.MPolyFactor.pfrac(I, t, degs, true)
        if !ok
            return false, fac
        end
        for i in 1:r
            fac[i] += map_up(R, deltas[i], mKp, l)
        end
    end
    return true, fac
end

function example(k::AnticNumberField, d::Int, nt::Int, c::AbstractRange=-10:10)
  kx, (x, y) = PolynomialRing(k, 2, cached = false)
  f = kx()
  for i=1:nt
    f += rand(k, c)*x^rand(0:d)*y^rand(0:d)
  end
  return norm(f)
end

function Hecke.isirreducible(a::fmpq_mpoly)
  af = factor(a)
  return !(length(af.fac) > 1 || any(x->x>1, values(af.fac)))
end

# f is bivariate. return f(xvar, 0) where xvar is in the multivar ring R
function _yzero_image(R, f, xvar::Int)
  z = MPolyBuildCtx(R)
  zexps = zeros(Int, nvars(R))
  for (c, exps) in zip(coefficients(f), exponent_vectors(f))
    @assert length(exps) == 2
    if exps[2] == 0
      zexps[xvar] = exps[1]
      push_term!(z, base_ring(R)(c), zexps)
    end
  end
  return finish(z)
end

function absolute_multivariate_factorisation(a::fmpq_mpoly)

  Qxy, (x, y) = PolynomialRing(QQ, ["x", "y"], cached = false)

  R = parent(a)
  K = base_ring(R)

  alphas = zeros(ZZ, nvars(R))
  bi_sub = zeros(Qxy, nvars(R))

  @assert length(a) > 0

  unit = coeff(a, 1)
  if !isone(unit)
    a *= inv(unit)
  end

  degs = degrees(a)
  vars = Int[]    # variables that actually appear
  for v in 1:nvars(R)
    if degs[v] > 0
      push!(vars, v)
    end
  end

  if isempty(vars)
    @assert isone(a)
    return (unit, [])
  end

  sort!(vars, by = (v -> degs[v]), alg=InsertionSort)

  if degs[1] == 1
    # linear is irreducible by assumption
    return (unit, [a])
  elseif length(vars) == 1
    uni_sub = zeros(Hecke.Globals.Qx, nvars(R))
    uni_sub[vars[1]] = gen(Hecke.Globals.Qx)
    K1, alpha = number_field(evaluate(a, uni_sub), cached = false)
    R1 = PolynomialRing(K1, map(string, symbols(R)), ordering = ordering(R), cached = false)[1]
    A = map_coefficients(K1, a, parent = R1)
    x = gen(R1, vars[1])
    return (unit, [x - alpha, divexact(A, x - alpha)])
  elseif length(vars) == 2
    bi_sub[vars[1]] = x
    bi_sub[vars[2]] = y
    f, fbar = absolute_bivariate_factorisation(evaluate(a, bi_sub))
    K1 = base_ring(f)
    R1 = PolynomialRing(K1, map(string, symbols(R)), ordering = ordering(R), cached = false)[1]
    revsub = [gen(R1, vars[1]), gen(R1, vars[2])]
    return (unit, [evaluate(f, revsub), evaluate(fbar, revsub)])
  end

  maindeg = degree(a, vars[1])
  mainvar = vars[1]
  minorvars = vars[2:end]

  lcaf = factor_squarefree(Hecke.AbstractAlgebra.MPolyFactor.get_lc(a, mainvar))

  lcc_fails_remaining = 3
  bits = 1

@label next_alpha

  if bits > 1000
    error("too many iterations")
  end

  bi_sub[mainvar] = x
  for i in 1:length(minorvars)
    alphas[i] = rand_bits(ZZ, rand(1:bits))
    bi_sub[minorvars[i]] = Qxy(alphas[i])
  end

  uni_a = evaluate(a, bi_sub)

  if degree(uni_a, 1) != maindeg || !isirreducible(uni_a)
    bits += 1
    @goto next_alpha
  end

  bi_sub_degree = 1

  if bi_sub_degree > 1 + bits/8
    bits += 2
    @goto next_alpha
  end

  # The substitution down to univariate is good and produces an irreducible.
  # Therefore, whatever bivariate is generated, the primitive part will be
  # irreducible if the y = 0 image agrees with the univariate substitution.
  bi_sub[mainvar] = x
  for i in 1:length(minorvars)
    bi_sub[minorvars[i]] = Qxy(alphas[i])
    for j in 1:bi_sub_degree
      bi_sub[minorvars[i]] += Qxy(rand_bits(ZZ, rand(1:bits)))*y^j
    end
  end

  bi_a = evaluate(a, bi_sub)
  bi_a = Hecke.AbstractAlgebra.MPolyFactor.primitive_part(bi_a, 1)
  if degree(bi_a, 2) < 2
    if degree(bi_a, 2) == 1
      # a is abs irreducible
      return (unit, [a])
    end
    @goto next_alpha
  end

  f, fbar = absolute_bivariate_factorisation(bi_a)

  K1 = base_ring(f)
  R1 = PolynomialRing(K1, map(string, symbols(R)), ordering = ordering(R), cached = false)[1]
  f = _yzero_image(R1, f, mainvar)
  fbar = _yzero_image(R1, fbar, mainvar)

  if degree(f, mainvar) < 1 || degree(fbar, mainvar) < 1
    # a is abs irreducible
    return (unit, [a])
  end

  # map the stuff in Q to the number field
  A = map_coefficients(K1, a, parent = R1)
  lcAf = Fac{elem_type(R1)}()
  lcAf.unit = map_coefficients(K1, lcaf.unit, parent = R1)
  for i in lcaf.fac
    lcAf[map_coefficients(K1, i[1], parent = R1)] = i[2]
  end

  ok, divs = Hecke.AbstractAlgebra.MPolyFactor.lcc_kaltofen(
                       lcAf, A, mainvar, maindeg, minorvars, alphas, [f, fbar])
  if !ok
    lcc_fails_remaining -= 1
    if lcc_fails_remaining >= 0
      @goto next_alpha
    end
  end

  ok, fac = Hecke.AbstractAlgebra.MPolyFactor.hlift_with_lcc(
                                A, [f, fbar], divs, mainvar, minorvars, alphas)
  if !ok
    @goto next_alpha
  end

  @assert length(fac) == 2
  return (unit, map(Hecke.AbstractAlgebra.MPolyFactor.make_monic, fac))
end

"""
    factor_absolute(f::fmpq_mpoly)

Compute the factorisation of f in Q[X] over C, returns an array with first
component the leading coefficient, and then, for each irreducible factor over Q[X]
a tuple, containing
- an irreducible factor over a number field
- the product of all the other conjugate factors over this field
- the multiplicity.

Recall that for each irreducible over Q there is one (minimal) number field
(possibly Q) where this factor splits further. In this case all the factor over
this field are Galois-conjugate. In order to not create a huge splitting field,
only one factor and the product of the conjugates is returned.

# Examples
```julia
julia> Qx, (x, y) = PolynomialRing(QQ, ["x", "y"]);

julia> f = (x^3+5*y^3)*(x^2+2*y^2);

julia> z = factor_absolute(f)

3-element Vector{Any}:
                                                                            1
                AbstractAlgebra.Generic.MPoly{nf_elem}[x + _a*y, x - _a*y] => 1
 AbstractAlgebra.Generic.MPoly{nf_elem}[x + _a*y, x^2 - _a*x*y + _a^2*y^2] => 1

julia> z[2][1][1]
x + _a*y

julia> base_ring(ans)
Number field over Rational Field with defining polynomial x^2 + 2

julia> base_ring(z[3][1][1])
Number field over Rational Field with defining polynomial x^3 - 5
```
"""
function factor_absolute(a::fmpq_mpoly)
  result = Any[]
  @vprint :AbsFact 1 "factoring over QQ first...\n"
  @vtime :AbsFact 2 fa = factor(a)
  push!(result, fa.unit)
  for (p, e) in fa
    @vtime :AbsFact 2 unit, fp = absolute_multivariate_factorisation(p)
    result[1] *= unit
    push!(result, fp => e)
  end
  return result
end

end

using .MPolyFact
export factor_absolute

#application (for free)

function factor(f::Union{fmpq_mpoly, fmpz_mpoly}, C::AcbField)
  fa = factor_absolute(f)
  D = Dict{Generic.MPoly{acb}, Int}()
  Cx, x = PolynomialRing(C, map(String, symbols(parent(f))), cached = false)
  for i=2:length(fa)
    K = base_ring(fa[i][1][1])
    if K == FlintQQ
      D[map_coefficients(C, fa[i][1][1], parent = Cx)] = fa[i][2]
      continue
    end
    g = [MPolyBuildCtx(Cx) for i=1:degree(K)]
    for (c, e) = zip(coefficients(fa[i][1][1]), exponent_vectors(fa[i][1][1]))
      d = Hecke.conjugates(c, precision(C))
      for j=1:degree(K)
        push_term!(g[j], C(d[j]), e)
      end
    end
    for j=1:degree(K)
      D[finish(g[j])] = fa[i][2]
    end
  end
  return Fac(map_coefficients(C, fa[1], parent = Cx), D)
end

function factor(f::Union{fmpq_mpoly, fmpz_mpoly}, R::ArbField)
  fa = factor_absolute(f)
  D = Dict{Generic.MPoly{arb}, Int}()
  Rx, x = PolynomialRing(R, map(String, symbols(parent(f))), cached = false)
  C = ComplexField(precision(R))
  Cx, x = PolynomialRing(C, map(String, symbols(parent(f))), cached = false)

  for i=2:length(fa)
    K = base_ring(fa[i][1][1])
    if K == FlintQQ
      D[map_coefficients(R, fa[i][1][1], parent = Rx)] = fa[i][2]
      continue
    end
    g = [MPolyBuildCtx(Cx) for i=1:degree(K)]
    for (c, e) = zip(coefficients(fa[i][1][1]), exponent_vectors(fa[i][1][1]))
      d = Hecke.conjugates(c, precision(C))
      for j=1:degree(K)
        push_term!(g[j], C(d[j]), e)
      end
    end
    r,s = signature(K)
    for j=1:r
      D[map_coefficients(x->R(real(x)), finish(g[j]), parent = Rx)] = fa[i][2]
    end
    for j=1:s
      gg = finish(g[r+j])*finish(g[r+j+s])
      D[map_coefficients(x->R(real(x)), gg, parent = Rx)] = fa[i][2]
    end
  end
  return Fac(map_coefficients(R, fa[1], parent = Rx), D)
end

#= revised strategy until I actually understand s.th. better
 - "shift" to make f(o, y) square-free
 - find a prime s.th. f(0, y) is square-free with a small degree splitting field
 - compute roots in F_q[[x]] (finite field!)
   TODO: first factor f in F_p[[x]], then compute roots of the factors
         should be faster (see above, more comments)
 - use combine to find the combinations giving the proper factor
 - find the (small) field the factor lives over
 - lift the factor in the q-adic field
 - find the number field

 TODO:
  bounds on the precisions (poly prec)

example:

Qxy, (y, x) = PolynomialRing(QQ, ["y", "x"])
include("/home/fieker/Downloads/n60s3.m");
  #from  https://www.math.univ-toulouse.fr/~cheze/n60s3.m

  r = absolute_bivariate_factorisation(P)

  from the Rupprecht paper, but the 3rd is boring (probably wrong)

  y^9+3*x^5*y^6+5*x^4*y^5+3*x^10*y^3-3*x^6*y^3+5*x^9*y^2+x^15

  y^5*x^5+9*x^8*y^4-6*x^14*y^2-18*x^10*y^6-18*x^6*y^10+x^20+5*x^16*y^4+10*x^12*y^8+10*x^8*y^12+5*y^16*x^4+9*y^8*x^4-6*y^14*x^2+y^20

  -125685*x+151959*x^8+917230*x^6+8717398*y^5*x^5+5108544*x^8*y^4-1564434*x^5+7744756*x^5*y^6+306683*x^3*y^6+413268*x^4*y^6+9081976*x^6*y^6+1317780*x^6*y^5+76745*x^4*y^5-15797040*x^7*y^5+99348*x^3*y^5+4106178*x^6*y^4+2010995*x^4*y^4-11264228*x^7*y^4-12465712*x^5*y^4+40908*x^2*y^4+404227*x^3*y^4-9204694*x^7*y^3-49266*x^2*y^3-3500343*x^4*y^3+1512264*x^3*y^3+6405504*x^8*y^3+9879662*x^6*y^3-3821606*x^5*y^3-592704*x^9*y^3-8503779*x^5*y^2-783216*x^9*y^2+10608275*x^6*y^2+574917*x^2*y^2-10143*x*y^2+5943180*x^4*y^2-3295022*x^3*y^2+3452692*x^8*y^2-6432756*x^7*y^2-344988*x^9*y+67473*x*y+2548458*x^4*y-2646351*x^7*y+1059606*x^8*y-3698541*x^5*y-491400*x^2*y+430155*x^3*y+4011984*x^6*y+1530912*x^4+617526
=#


text:
include("env.jl")
include("resources.jl")
include("devtools.jl")
include("dash_creation.jl")
include("handlers.jl")
include("config_functional.jl")
include("utils.jl")
include("context.jl")
include("core.jl")
include("misc.jl")
include("callbacks.jl")
include("components_utils.jl")
include("table_format.jl")
include("reload_hash.jl")
#include("dev.jl")


text:
abstract type AbstractParameters end

#holds the main parameter object
mutable struct Parameters <: AbstractParameters
  σ²G::Float64
  ζ²G::Float64
  ζ²P::Float64
  ζGP::Float64
  SG::Float64
  SGP::Float64 #not an actual
end

#named constructor for testing purposes
function Parameters(;
  σ²G::Float64 = D_σ²G,
  ζ²G::Float64 = D_ζ²G,
  ζ²P::Float64 = D_ζ²P,
  ζGP::Float64 = D_ζGP,
  SG::Float64 = D_SG,
  SGP::Float64 = D_SGP)::Parameters

  return Parameters(σ²G, ζ²G, ζ²P, ζGP, SG, SGP)
end

function Base.print(Θ::Parameters)
  return  "σ²G = $(Θ.σ²G); ζ²G = $(Θ.ζ²G); ζ²P = $(Θ.ζ²P);
    ζGP = $(Θ.ζGP); SG = $(Θ.SG); SGP = $(Θ.SGP)"
end

Base.println(Θ::Parameters)::String = Base.print(Θ) * "\n"

copy2vector(Θ::Parameters)::Vector{Float64} = [Θ.σ²G, Θ.ζ²G, Θ.ζ²P, Θ.ζGP, Θ.SG, Θ.SGP]

Base.names(Θ::Parameters) = ["sigma2G", "zeta2G", "zeta2P", "zetaGP", "SG", "SGP"]

#holds the priors
struct Priors <: AbstractParameters
  θG::Float64
  δ²G::Float64
  Ψ::Matrix{Float64}
  ν::Float64
  θSG::Float64
  δ²SG::Float64
  θSGP::Float64
  δ²SGP::Float64
end

#named constructor for testing purposes
function Priors(;priortype::Symbol = :default)

  local θG::Float64
  local δ²G::Float64
  local Ψ::Matrix{Float64}
  local ν::Float64
  local θSG::Float64
  local δ²SG::Float64
  local θSGP::Float64
  local δ²SGP::Float64

  if priortype == :default || priortype == :informative
    θG = D_θG
    δ²G = D_δ²G
    Ψ = D_Ψ
    ν = D_ν
    θSG = D_θSG
    δ²SG = D_δ²SG
    θSGP = D_θSGP
    δ²SGP = D_δ²SGP
  elseif priortype == :uninformative
    θG = U_θG
    δ²G = U_δ²G
    Ψ = U_Ψ
    ν = U_ν
    θSG = U_θSG
    δ²SG = U_δ²SG
    θSGP = U_θSGP
    δ²SGP = U_δ²SGP
  else
    error("Prior type $priortype not recognized.")
  end

  return Priors(θG, δ²G, Ψ, ν, θSG, δ²SG, θSGP, δ²SGP)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using GLVisualize, Reactive, GeometryTypes
using GLWindow, GLAbstraction

if !isdefined(:runtests)
	window = glscreen()
end
const static_example = true

"""
functions to halve some rectangle
"""
xhalf(r)  = SimpleRectangle(r.x, r.y, r.w÷2, r.h)
xhalf2(r) = SimpleRectangle(r.w÷2, r.y, r.w÷2, r.h)

# create first screen with window as the parent screen
# and a different area.
screen2D = Screen(
	window, name=:screen2D,
	area=const_lift(xhalf2, window.area)
)
# create second screen with window as the parent screen
# and a different area.
screen3D = Screen(
	window, name=:screen3D,
	area=const_lift(xhalf, window.area)
)

# create something to look at!
bars = visualize(rand(Float32, 10,10))
points = visualize([rand(Point2f0)*1000 for i=1:50], scale=Vec2f0(40))

# _view them in different screens
_view(bars,   screen3D, camera=:perspective)
_view(points, screen2D, camera=:orthographic_pixel)

if !isdefined(:runtests)
	renderloop(window)
end


text:
using PkgEval
using Test

# determine the version to use
const version = get(ENV, "JULIA_VERSION", string(VERSION))
@info "Testing with Julia $version"
const julia = PkgEval.obtain_julia(version)::VersionNumber
@info "Resolved to Julia v$julia"
const install = PkgEval.prepare_julia(julia)

@testset "sandbox" begin
    PkgEval.prepare_runner()
    mktemp() do path, io
        PkgEval.run_sandboxed_julia(install, `-e 'print(1337)'`; stdout=io,
                                       tty=false, interactive=false)
        close(io)
        @test read(path, String) == "1337"
    end

    # print versioninfo so we can verify in CI logs that the correct version is used
    PkgEval.run_sandboxed_julia(install, `-e 'using InteractiveUtils; versioninfo()'`;
                                   tty=false, interactive=false)
end

const pkgnames = ["TimerOutputs", "Crayons", "Example", "Gtk"]

@testset "low-level interface" begin
    PkgEval.prepare_registry()

    pkgs = PkgEval.read_pkgs(pkgnames)

    # timeouts
    results = PkgEval.run([julia], pkgs; time_limit = 0.1)
    @test all(results.status .== :kill)
end

@testset "main entrypoint" begin
    results = PkgEval.run([julia], pkgnames)
    if !(version == "master" || version == "nightly")
        @test all(results.status .== :ok)
        for result in eachrow(results)
            @test occursin("Testing $(result.name) tests passed", result.log)
        end
    end
end

@testset "reporting" begin
    lts = v"1.0.5"
    stable = v"1.2.0"
    results = PkgEval.run([lts, stable], ["Example"])
    PkgEval.compare(results, lts, stable)
end

PkgEval.purge()
rm(install; recursive=true)


text:
using AlphaShapes
using Test

tol = 1e-6

@testset "Test Circum Sphere" begin
   p = [0. 1.;1. 0.;1. 1.]
   c,r = AlphaShapes.SimplexCircumSphere(p)
   @test (r .- 0.7071067811865476) < tol
   @test (c[1] .- 0.5) < tol
   @test (c[1].-0.5) < tol
end


text:
using Distributions
import Distributions.logpdf
using Base.Cartesian

export logpdf
export rand

export For
struct For{F,T,D,X}
    f :: F
    θ :: T
end

#########################################################
# T <: NTuple{N,J} where {J <: Integer}
#########################################################

For(f, θ::J...) where {J <: Integer} = For(f,θ)

function For(f::F, θ::T) where {F, N, J <: Integer, T <: NTuple{N,J}}
    d = f.(ones(Int, N)...)
    D = typeof(d)
    X = eltype(d)
    For{F, NTuple{N,J}, D, X}(f,θ)
end

@inline function logpdf(d::For{F,T,D,X1},xs::AbstractArray{X2,N}) where {F, N, J <: Integer, T <: NTuple{N,J}, D,  X1,  X2 <: X1}
    s = 0.0
    @inbounds @simd for θ in CartesianIndices(d.θ)
        s += logpdf(d.f(Tuple(θ)...), xs[θ])
    end
    s
end

function Base.rand(dist::For)
    map(CartesianIndices(dist.θ)) do I
        (rand ∘ dist.f)(Tuple(I)...)
    end
end

#########################################################
# T <: NTuple{N,J} where {J <: AbstractUnitRange}
#########################################################

For(f, θ::J...) where {J <: AbstractUnitRange} = For(f,θ)

function For(f::F, θ::T) where {F, N, J <: AbstractRange, T <: NTuple{N,J}}
    d = f.(ones(Int, N)...)
    D = typeof(d)
    X = eltype(d)
    For{F, NTuple{N,J}, D, X}(f,θ)
end

@inline function logpdf(d::For{F,T,D,X1},xs::AbstractArray{X2,N}) where {F, N, J <: AbstractRange,  T <: NTuple{N,J}, D, X1, X2 <: X1}
    s = 0.0
    @inbounds @simd for θ in CartesianIndices(d.θ)
        s += logpdf(d.f(Tuple(θ)...), xs[θ])
    end
    s
end

function Base.rand(dist::For{F,T}) where {F,  N, J <: AbstractRange, T <: NTuple{N,J}}
    map(CartesianIndices(dist.θ)) do I
        (rand ∘ dist.f)(Tuple(I)...)
    end
end

#########################################################
# T <: Base.Generator
#########################################################

function For(f::F, θ::T) where {F, T <: Base.Generator}
    d = f(θ.f(θ.iter[1]))
    D = typeof(d)
    X = eltype(d)
    For{F, T, D, X}(f,θ)
end

@inline function logpdf(d :: For{F,T}, x) where {F,T <: Base.Generator}
    s = 0.0
    for (θj, xj) in zip(d.θ, x)
        s += logpdf(d.f(θj), xj)
    end
    s
end

@inline function rand(d :: For{F,T,D,X}) where {F,T <: Base.Generator, D, X}
    rand.(Base.Generator(d.θ.f, d.θ.iter))
end

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


text:
# Kernel definition is taken from here: https://en.wikipedia.org/wiki/Kernel_(statistics)#Kernel_functions_in_common_use

abstract type CollocationKernel end
struct EpanechnikovKernel <: CollocationKernel end
struct UniformKernel <: CollocationKernel end
struct TriangularKernel <: CollocationKernel end
struct QuarticKernel <: CollocationKernel end
struct TriweightKernel <: CollocationKernel end
struct TricubeKernel <: CollocationKernel end
struct GaussianKernel <: CollocationKernel end
struct CosineKernel <: CollocationKernel end
struct LogisticKernel <: CollocationKernel end
struct SigmoidKernel <: CollocationKernel end
struct SilvermanKernel <: CollocationKernel end

function calckernel(::EpanechnikovKernel,t)
    if abs(t) > 1
        return 0
    else
        return 0.75*(1-t^2)
    end
end

function calckernel(::UniformKernel,t)
    if abs(t) > 1
        return 0
    else
        return 0.5
    end
end

function calckernel(::TriangularKernel,t)
    if abs(t) > 1
        return 0
    else
        return (1-abs(t))
    end
end

function calckernel(::QuarticKernel,t)
  if abs(t)>0
    return 0
  else
    return (15*(1-t^2)^2)/16
  end
end

function calckernel(::TriweightKernel,t)
  if abs(t)>0
    return 0
  else
    return (35*(1-t^2)^3)/32
  end
end

function calckernel(::TricubeKernel,t)
  if abs(t)>0
    return 0
  else
    return (70*(1-abs(t)^3)^3)/80
  end
end

function calckernel(::GaussianKernel,t)
  exp(-0.5*t^2)/(sqrt(2*π))
end

function calckernel(::CosineKernel,t)
  if abs(t)>0
    return 0
  else
    return (π*cos(π*t/2))/4
  end
end

function calckernel(::LogisticKernel,t)
  1/(exp(t)+2+exp(-t))
end

function calckernel(::SigmoidKernel,t)
  2/(π*(exp(t)+exp(-t)))
end

function calckernel(::SilvermanKernel,t)
  sin(abs(t)/2+π/4)*0.5*exp(-abs(t)/sqrt(2))
end


text:
struct VLAE{E,D,G,F,XD,ZD,V,XDIST<:Val} <: AbstractVLAE
    e::E
    d::D
    g::G # extracts latent variables
    f::F # concatenates latent vars with the rest
    xdim::XD # (h,w,c)
    zdim::ZD # scalar
    var::V # dense or conv last layer
    xdist::XDIST # gaussian or bernoulli
end
VLAE(e,d,g,f,xdim,zdim,var,xdist::Symbol=:gaussian) =
    VLAE(e,d,g,f,xdim,zdim,var,Val(xdist))
Flux.@functor VLAE
(m::VLAE)(x) = reconstruct(m, x)

"""
	VLAE(zdim::Int, ks, ncs, stride, datasize; layer_depth=1, var=:dense, activation="relu",
        xdist=:gaussian)

Basic VLAE constructor.
"""
function VLAE(zdim::Int, ks, ncs, strd, datasize; var=:dense, xdist=:gaussian, kwargs...)
    (xdist in [:gaussian, :bernoulli]) ? nothing :
        error("xdist must be either :gaussian or :bernoulli")

	# get encoder, decoder and latent extractors and reshapers
   	e,d,g,f = basic_model_constructor(zdim, ks, ncs, strd, datasize; var=var, xdist=xdist, kwargs...)

    return VLAE(e,d,g,f,datasize[1:3],zdim,var,xdist)
end

"""
	train_vlae(zdim, batchsize, ks, ncs, stride, nepochs, data, val_x, tst_x;
    	λ=0.0f0, epochsize = size(data,4), layer_depth=1, lr=0.001f0, var=:dense,
    	activation=activation)

"""
function train_vlae(zdim, batchsize, ks, ncs, strd, nepochs, data, val_x, tst_x;
    λ=0.0f0, epochsize = size(data,4), layer_depth=1, lr=0.001f0, var=:dense,
    activation=activation, xdist=:gaussian)
    gval_x = gpu(val_x[:,:,:,1:min(1000, size(val_x,4))])
    gtst_x = gpu(tst_x)

    model = gpu(VLAE(zdim, ks, ncs, strd, size(data), layer_depth=layer_depth, var=var,
        activation=activation, xdist=xdist))
    nl = length(model.e)

    ps = Flux.params(model)
    loss(x) = -elbo(model,gpu(x)) + λ*sum(l2, ps)
    opt = ADAM(lr)
    rdata = []
    hist = MVHistory()
    zs = [[] for _ in 1:nl]

    # train
    println("Training in progress...")
    for epoch in 1:nepochs
        local data_itr
        @suppress begin # suppress the batchsize warning
            data_itr = Flux.Data.DataLoader(data[:,:,:,sample(1:size(data,4), epochsize)], batchsize=batchsize)
        end
        Flux.train!(loss, ps, data_itr, opt)
        l = Flux.mean(map(x->loss(x), Flux.Data.DataLoader(val_x, batchsize=batchsize)))
        println("Epoch $(epoch)/$(nepochs), validation loss = $l")
        for i in 1:nl
            z = encode(model, gval_x, i)
            push!(zs[i], cpu(z))
        end
        push!(hist, :loss, epoch, l)
        push!(rdata, cpu(reconstruct(model, gval_x)))
    end

    return model, hist, rdata, zs
end

function logpdf(m::VLAE{E,D,G,F,X,Z,V,XD}, x, zs...) where {E,D,G,F,X,Z,V,XD<:Val{:gaussian}}
    μx, σx = _decoded_mu_var(m, zs...)
    _x = (m.var == :dense) ? vectorize(x) : x
    return  gauss_logpdf(_x, μx, σx)
end
function logpdf(m::VLAE{E,D,G,F,X,Z,V,XD}, x, zs...) where {E,D,G,F,X,Z,V,XD<:Val{:bernoulli}}
    _x = _decoder_out(m, zs...)
    return  bernoulli_logpdf(_x, x)
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
module Whirl
using LinearAlgebra
using StaticArrays
using Statistics
using Plots
import Base

include("num.jl")
include("error.jl")
include("debug.jl")
include("geometry.jl")
include("units.jl")
include("eos.jl")
include("kernel.jl")
include("initial.jl")
include("velocity.jl")
include("mass.jl")
include("energy.jl")
include("pressure.jl")
include("diffusion.jl")
include("boundary.jl")
include("fluid.jl")
include("stepping.jl")
include("simulation.jl")
include("experiments.jl")
include("plotting.jl")

include("test.jl")

end # module

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using Interpolations, Test

@testset "ScalingWithExtrapTests" begin

    xs = range(-5, stop=5, length=10)
    ys = map(sin, xs)

    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp, ::OnGrid) where {T,N,IT}
        for x in xs
            @test ≈(sut(x),sin(x),atol=sqrt(eps(sin(x))))
        end
        @test sut(-5) == sut(-5.1) == sut(-15.8) == sut(-Inf) == itp(1)
        @test sut(5) == sut(5.1) == sut(15.8) == sut(Inf) == itp[end]
    end

    function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT}, itp, ::OnCell) where {T,N,IT}
        halfcell = (xs[2] - xs[1]) / 2
        itps, axs = Interpolations.itpinfo(itp)
        for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
            @test sut(-x) == itp(Interpolations.lbound(axs[1], itps[1]))
            @test sut(x) == itp(Interpolations.ubound(axs[1], itps[1]))
        end
    end

    run_tests(sut, itp) = run_tests(sut, itp, Interpolations.degree(Interpolations.itpflag(itp)).bc.gt)

    for GT in (OnGrid, OnCell)
        itp = interpolate(ys, BSpline(Quadratic(Flat(GT()))))

        # Test extrapolating, then scaling
        eitp = extrapolate(itp, Flat())
        seitp = scale(eitp, xs)
        run_tests(seitp, itp)

        # Test scaling, then extrapolating
        sitp = scale(itp, xs)
        esitp = extrapolate(sitp, Flat())
        run_tests(esitp, itp)
    end

end