text:
#=
Functions for generating a single replication of a particular experiment
Dependencies: util.jl, sgd.jl, reformulation.jl, random_forests_po.jl, validation_set.jl
=#

using Parameters, DataFrames, Statistics, Random, Gurobi

## Type for strong the results of an experiment replication
@with_kw mutable struct replication_results
    SPOplus_spoloss_test::Union{Float64, Missing} = missing # SPO loss of SPO+ SGD on test set
    LS_spoloss_test::Union{Float64, Missing} = missing # SPO loss of LS SGD on test set
    SSVM_spoloss_test::Union{Float64, Missing} = missing # SPO loss of SSVM-Hamming SGD on test set
    RF_spoloss_test::Union{Float64, Missing} = missing # SPO loss of random forests on test set
    Absolute_spoloss_test::Union{Float64, Missing} = missing # SPO loss of ell_1 on test set
    Huber_spoloss_test::Union{Float64, Missing} = missing  # SPO loss of Huber on test set
    Baseline_spoloss_test::Union{Float64, Missing} = missing # SPO loss of baseline that returns mean of training data, on test set
    SPOplus_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of SPO+ SGD on test set
    LS_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of LS SGD on test set
    SSVM_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of SSVM SGD on test set
    RF_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of RF on test set
    Absolute_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of ell_1 on test set
    Huber_hammingloss_test::Union{Float64, Missing} = missing  # Hamming loss of Huber on test set
    Baseline_hammingloss_test::Union{Float64, Missing} = missing # Hamming loss of baseline on test set
    zstar_avg_test::Union{Float64, Missing} = missing # Average of z^*(c_i) on testing set
end

"""
Make a blank data frame for storing the results of the experiment.
"""
function make_blank_complete_df()
    results_df = DataFrame(
        grid_dim = Int64[],
        n_train = Int64[],
        n_holdout = Int64[],
        n_test = Int64[],
        p_features = Int64[],
        polykernel_degree = Int64[],
        polykernel_noise_half_width = Float64[],
        SPOplus_spoloss_test = Union{Float64, Missing}[], # SPO loss of SPO+ SGD on test set
        LS_spoloss_test = Union{Float64, Missing}[],  # SPO loss of LS SGD on test set
        SSVM_spoloss_test = Union{Float64, Missing}[], # SPO loss of SSVM-Hamming SGD on test set
        RF_spoloss_test = Union{Float64, Missing}[], # SPO loss of random forests on test set
        Absolute_spoloss_test = Union{Float64, Missing}[], # SPO loss of ell_1 on test set
        Huber_spoloss_test = Union{Float64, Missing}[],  # SPO loss of Huber on test set
        Baseline_spoloss_test = Union{Float64, Missing}[], # SPO loss of baseline that returns mean of training data, on test set
        SPOplus_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of SPO+ SGD on test set
        LS_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of LS SGD on test set
        SSVM_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of SSVM SGD on test set
        RF_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of RF on test set
        Absolute_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of ell_1 on test set
        Huber_hammingloss_test = Union{Float64, Missing}[],  # Hamming loss of Huber on test set
        Baseline_hammingloss_test = Union{Float64, Missing}[], # Hamming loss of baseline on test set
        zstar_avg_test = Union{Float64, Missing}[]) # Average of z^*(c_i) on testing set

    return results_df
end

"""
Make a new data frame that contains one row of results for an experiment instance that was just run.
"""
function build_complete_row(grid_dim,
    n_train, n_holdout, n_test,
    p_features, polykernel_degree, polykernel_noise_half_width, results_struct)

    df = make_blank_complete_df()

    push!(df, [
        grid_dim,
        n_train,
        n_holdout,
        n_test,
        p_features,
        polykernel_degree,
        polykernel_noise_half_width,
        results_struct.SPOplus_spoloss_test,
        results_struct.LS_spoloss_test,
        results_struct.SSVM_spoloss_test,
        results_struct.RF_spoloss_test,
        results_struct.Absolute_spoloss_test,
        results_struct.Huber_spoloss_test,
        results_struct.Baseline_spoloss_test,
        results_struct.SPOplus_hammingloss_test,
        results_struct.LS_hammingloss_test,
        results_struct.SSVM_hammingloss_test,
        results_struct.RF_hammingloss_test,
        results_struct.Absolute_hammingloss_test,
        results_struct.Huber_hammingloss_test,
        results_struct.Baseline_hammingloss_test,
        results_struct.zstar_avg_test])

    return df
end

"""
Run one instance of the shortest path experiment. Uses the reformulation approach.
`num_lambda` is the number of lambdas used on the grid for each method with regularization.
Returns a replication_results struct.
"""
function shortest_path_replication(grid_dim,
    n_train, n_holdout, n_test,
    p_features, polykernel_degree, polykernel_noise_half_width;
    num_lambda = 10, lambda_max = missing, lambda_min_ratio = 0.0001, regularization = :ridge,
    gurobiEnvOracle = Gurobi.Env(), gurobiEnvReform = Gurobi.Env(), different_validation_losses = false,
    include_rf = true)

    # Get oracle
    sources, destinations = convert_grid_to_list(grid_dim, grid_dim)
    d_feasibleregion = length(sources)
    sp_oracle = sp_flow_jump_setup(sources, destinations, 1, grid_dim^2; gurobiEnv = gurobiEnvOracle)
    sp_graph = shortest_path_graph(sources = sources, destinations = destinations,
        start_node = 1, end_node = grid_dim^2, acyclic = true)

    # Generate Data
    B_true = rand(Bernoulli(0.5), d_feasibleregion, p_features)

    (X_train, c_train) = generate_poly_kernel_data_simple(B_true, n_train, polykernel_degree, polykernel_noise_half_width)
    (X_validation, c_validation) = generate_poly_kernel_data_simple(B_true, n_holdout, polykernel_degree, polykernel_noise_half_width)
    (X_test, c_test) = generate_poly_kernel_data_simple(B_true, n_test, polykernel_degree, polykernel_noise_half_width)

    # Add intercept in the first row of X
    X_train = vcat(ones(1,n_train), X_train)
    X_validation = vcat(ones(1,n_holdout), X_validation)
    X_test = vcat(ones(1,n_test), X_test)

    # Get Hamming labels
    (z_train, w_train) = oracle_dataset(c_train, sp_oracle)
    (z_validation, w_validation) = oracle_dataset(c_validation, sp_oracle)
    (z_test, w_test) = oracle_dataset(c_test, sp_oracle)

    c_ham_train = ones(d_feasibleregion, n_train) - w_train
    c_ham_validation = ones(d_feasibleregion, n_holdout) - w_validation
    c_ham_test = ones(d_feasibleregion, n_test) - w_test

    # Put train + validation together
    X_both = hcat(X_train, X_validation)
    c_both = hcat(c_train, c_validation)
    c_ham_both = hcat(c_ham_train, c_ham_validation)
    train_ind = collect(1:n_train)
    validation_ind = collect((n_train + 1):(n_train + n_holdout))

    # Set validation losses
    if different_validation_losses
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :least_squares_loss
        ssvm_val_loss = :hamming_loss
        absolute_val_loss = :absolute_loss
        huber_val_loss = :huber_loss
    else
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :spo_loss
        ssvm_val_loss = :spo_loss
        absolute_val_loss = :spo_loss
        huber_val_loss = :spo_loss
    end

    ### Algorithms ###

    # SPO+
    best_B_SPOplus, best_lambda_SPOplus = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :sp_spo_plus_reform, validation_loss = spo_plus_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :SPO_plus))

    # Least squares
    best_B_leastSquares, best_lambda_leastSquares = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :leastSquares))

    # SSVM Hamming
    best_B_SSVM, best_lambda_SSVM = validation_set_alg(X_both, c_ham_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :sp_spo_plus_reform, validation_loss = ssvm_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :SSVM_hamming))

    # RF
    if include_rf
        rf_mods = train_random_forests_po(X_train, c_train;
            rf_alg_parms = rf_parms())
    end

    # Absolute
    best_B_Absolute, best_lambda_Absolute = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = absolute_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization,
            po_loss_function = :absolute, gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :Absolute))

    # Huber
    best_B_fake, best_lambda_fake = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(lambda_max = 0.0001, regularization = regularization,
            num_lambda = 2, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber_LS_fake))

    fake_ls_list = abs.(vec(c_train - best_B_fake*X_train))
    delta_from_fake_ls = median(fake_ls_list)

    best_B_Huber, best_lambda_Huber = validation_set_alg(X_both, c_both, sp_oracle; sp_graph = sp_graph,
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = huber_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = regularization, po_loss_function = :huber,
            huber_delta = delta_from_fake_ls, lambda_min_ratio = lambda_min_ratio, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber))

    # Baseline
    c_bar_train = mean(c_train, dims=2)

    ### Populate final results ###
    final_results = replication_results()

    final_results.SPOplus_spoloss_test = spo_loss(best_B_SPOplus, X_test, c_test, sp_oracle)
    final_results.LS_spoloss_test = spo_loss(best_B_leastSquares, X_test, c_test, sp_oracle)
    final_results.SSVM_spoloss_test = spo_loss(best_B_SSVM, X_test, c_test, sp_oracle)

    if include_rf
        rf_preds_test = predict_random_forests_po(rf_mods, X_test)
        final_results.RF_spoloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), rf_preds_test, c_test, sp_oracle)
    else
        final_results.RF_spoloss_test = missing
    end

    final_results.Absolute_spoloss_test = spo_loss(best_B_Absolute, X_test, c_test, sp_oracle)
    final_results.Huber_spoloss_test = spo_loss(best_B_Huber, X_test, c_test, sp_oracle)

    c_bar_test_preds = zeros(d_feasibleregion, n_test)
    for i = 1:n_test
        c_bar_test_preds[:, i] = c_bar_train
    end
    final_results.Baseline_spoloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), c_bar_test_preds, c_test, sp_oracle)

    # Now Hamming
    final_results.SPOplus_hammingloss_test = spo_loss(best_B_SPOplus, X_test, c_ham_test, sp_oracle)
    final_results.LS_hammingloss_test = spo_loss(best_B_leastSquares, X_test, c_ham_test, sp_oracle)
    final_results.SSVM_hammingloss_test = spo_loss(best_B_SSVM, X_test, c_ham_test, sp_oracle)

    if include_rf
        final_results.RF_hammingloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), rf_preds_test, c_ham_test, sp_oracle)
    else
        final_results.RF_hammingloss_test = missing
    end

    final_results.Absolute_hammingloss_test = spo_loss(best_B_Absolute, X_test, c_ham_test, sp_oracle)
    final_results.Huber_hammingloss_test = spo_loss(best_B_Huber, X_test, c_ham_test, sp_oracle)
    final_results.Baseline_hammingloss_test = spo_loss(Matrix(1.0I, d_feasibleregion, d_feasibleregion), c_bar_test_preds, c_ham_test, sp_oracle)

    final_results.zstar_avg_test = mean(z_test)

    return final_results
end

"""
Run multiple replications of the shortest path experiment and return a data frame.
`rng_seed` and `num_trials` are integers
`n_train_vec`, `polykernel_degree_vec`, and `polykernel_noise_half_width_vec` are
vectors so that we try all combinations of these values.
"""
function shortest_path_multiple_replications(rng_seed, num_trials, grid_dim,
    n_train_vec, n_test,
    p_features, polykernel_degree_vec, polykernel_noise_half_width_vec;
    num_lambda = 10, lambda_max = missing, lambda_min_ratio = 0.0001, holdout_percent = 0.25, regularization = :ridge,
    different_validation_losses = false,
    include_rf = true)

    Random.seed!(rng_seed)

    big_results_df = make_blank_complete_df()

    envOracle = setup_gurobi_env(method_type = :default, use_time_limit = false)
    envReform = setup_gurobi_env(method_type = :method3, use_time_limit = false)

    for n_train in n_train_vec
        for polykernel_degree in polykernel_degree_vec
            for polykernel_noise_half_width in polykernel_noise_half_width_vec
                println("Moving on to n_train = $n_train, polykernel_degree = $polykernel_degree, polykernel_noise_half_width = $polykernel_noise_half_width")
                for trial = 1:num_trials
                    println("Current trial is $trial")
                    n_holdout = round(Int, holdout_percent*n_train)

                    current_results = shortest_path_replication(grid_dim,
                        n_train, n_holdout, n_test,
                        p_features, polykernel_degree, polykernel_noise_half_width;
                        num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio, regularization = regularization,
                        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform,
                        different_validation_losses = different_validation_losses,
                        include_rf = include_rf)

                    current_df_row = build_complete_row(grid_dim,
                        n_train, n_holdout, n_test,
                        p_features, polykernel_degree, polykernel_noise_half_width, current_results)

                    big_results_df = vcat(big_results_df, current_df_row)
                end
            end
        end
    end

    return big_results_df
end

### Portfolio stuff ###

"""
Run one instance of the portfolio experiment. Uses the SGD approach.
`num_lambda` is the number of lambdas used on the grid for each method with regularization.
Returns a replication_results struct.
"""
function portfolio_replication(num_assets, num_factors, n_train, n_holdout, n_test, p_features,
    n_sigmoid_polydegree, noise_multiplier_tau;
    num_lambda = 10, lambda_max = :practical, lambda_min_ratio = 0.0001,
    gurobiEnvOracle = Gurobi.Env(), gurobiEnvReform = Gurobi.Env(), different_validation_losses = false,
    data_type = :sigmoid)

    # Get oracle
    factor_size = 0.0025*noise_multiplier_tau
    sigma_noise = 0.01*noise_multiplier_tau

    L_mat = 2*factor_size*rand(num_assets, num_factors) - factor_size*ones(num_assets, num_factors)
    Sigma = L_mat*L_mat' + (sigma_noise^2)*Matrix(1.0I, num_assets, num_assets)

    avg_vec = ones(num_assets)/num_assets
    gamma = avg_vec'*Sigma*avg_vec
    gamma = 2.25*gamma

    port_oracle = portfolio_simplex_jump_setup(Sigma, gamma; gurobiEnv = gurobiEnvOracle)

    # Generate Data
    return_intercept = 0
    if data_type == :rbf_kernel
        sigma_kernel = sqrt(p_features)
        # Rademacher instead of Bernoulli
        B_true = 2*rand(Bernoulli(0.5), num_assets, p_features) - ones(num_assets, p_features)

        (X_train, c_train) = generate_rbf_kernel_returns_data(B_true, n_train, sigma_kernel, L_mat, sigma_noise, return_intercept)
        (X_validation, c_validation) = generate_rbf_kernel_returns_data(B_true, n_holdout, sigma_kernel, L_mat, sigma_noise, return_intercept)
        (X_test, c_test) = generate_rbf_kernel_returns_data(B_true, n_test, sigma_kernel, L_mat, sigma_noise, return_intercept)
    elseif data_type == :sigmoid
        bernoulli_prob = 0.5

        (X_all, c_all) = generate_sigmoid_returns_data(p_features, num_assets, n_sigmoid_polydegree, n_train + n_holdout + n_test, bernoulli_prob,
        	L_mat, sigma_noise, return_intercept)
        c_all = c_all ./ n_sigmoid_polydegree

        X_train = X_all[:, 1:n_train]
        c_train = c_all[:, 1:n_train]
        X_validation = X_all[:, (n_train + 1):(n_train + n_holdout)]
        c_validation = c_all[:, (n_train + 1):(n_train + n_holdout)]
        X_test = X_all[:, (n_train + n_holdout + 1):(n_train + n_holdout + n_test)]
        c_test = c_all[:, (n_train + n_holdout + 1):(n_train + n_holdout + n_test)]
    elseif data_type == :poly_kernel
        #### zzz need to do properly tomorrow
        B_true = rand(Bernoulli(0.5), num_assets, p_features)
        #B_true = rand(num_assets, p_features)
        polykernel_degree = n_sigmoid_polydegree

        (X_train, c_train) = generate_poly_kernel_returns_data(B_true, n_train, polykernel_degree,
        	L_mat, sigma_noise, return_intercept)
        (X_validation, c_validation) = generate_poly_kernel_returns_data(B_true, n_holdout, polykernel_degree,
        	L_mat, sigma_noise, return_intercept)
        (X_test, c_test) = generate_poly_kernel_returns_data(B_true, n_test, polykernel_degree,
        	L_mat, sigma_noise, return_intercept)
    else
        error("Enter a valid data generation type")
    end

    #println("Train returns are:")
    #println(-c_train[:,1])

    # Add intercept in the first row of X
    X_train = vcat(ones(1,n_train), X_train)
    X_validation = vcat(ones(1,n_holdout), X_validation)
    X_test = vcat(ones(1,n_test), X_test)

    # Put train + validation together
    X_both = hcat(X_train, X_validation)
    c_both = hcat(c_train, c_validation)
    train_ind = collect(1:n_train)
    validation_ind = collect((n_train + 1):(n_train + n_holdout))

    # Set validation losses
    if different_validation_losses
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :least_squares_loss
        ssvm_val_loss = :hamming_loss
        absolute_val_loss = :absolute_loss
        huber_val_loss = :huber_loss
    else
        spo_plus_val_loss = :spo_loss
        ls_val_loss = :spo_loss
        ssvm_val_loss = :spo_loss
        absolute_val_loss = :spo_loss
        huber_val_loss = :spo_loss
    end

    ### Algorithms ###

    # SPO+
    best_B_SPOplus, best_lambda_SPOplus = validation_set_alg(X_both, c_both, port_oracle;
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :spo_plus_sgd, resolve_sgd = false, validation_loss = spo_plus_val_loss),
        path_alg_parms = sgd_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio,
            obj_accuracy = 0.001, step_type = :long_dynamic, verbose = true))

    # Least squares
    best_B_leastSquares, best_lambda_leastSquares = validation_set_alg(X_both, c_both, port_oracle;
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = :ridge,
            gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :leastSquares))

    # RF
    rf_mods = train_random_forests_po(X_train, c_train;
        rf_alg_parms = rf_parms())

    # Absolute
    best_B_Absolute, best_lambda_Absolute = validation_set_alg(X_both, c_both, port_oracle;
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = absolute_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = :ridge,
            po_loss_function = :absolute, gurobiEnv = gurobiEnvReform, lambda_min_ratio = lambda_min_ratio, algorithm_type = :Absolute))

    # Huber
    best_B_fake, best_lambda_fake = validation_set_alg(X_both, c_both, port_oracle;
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = ls_val_loss),
        path_alg_parms = reformulation_path_parms(lambda_max = 0.0001, regularization = :ridge,
            num_lambda = 2, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber_LS_fake))

    fake_ls_list = abs.(vec(c_train - best_B_fake*X_train))
    delta_from_fake_ls = median(fake_ls_list)

    best_B_Huber, best_lambda_Huber = validation_set_alg(X_both, c_both, port_oracle;
        train_ind = train_ind, validation_ind = validation_ind,
        val_alg_parms = val_parms(algorithm_type = :ls_jump, validation_loss = huber_val_loss),
        path_alg_parms = reformulation_path_parms(num_lambda = num_lambda, lambda_max = lambda_max, regularization = :ridge, po_loss_function = :huber,
            huber_delta = delta_from_fake_ls, lambda_min_ratio = lambda_min_ratio, gurobiEnv = gurobiEnvReform, algorithm_type = :Huber))

    # Baseline
    c_bar_train = mean(c_train, dims=2)

    ### Populate final results ###
    final_results = replication_results()

    final_results.SPOplus_spoloss_test = spo_loss(best_B_SPOplus, X_test, c_test, port_oracle)
    final_results.LS_spoloss_test = spo_loss(best_B_leastSquares, X_test, c_test, port_oracle)

    rf_preds_test = predict_random_forests_po(rf_mods, X_test)
    final_results.RF_spoloss_test = spo_loss(Matrix(1.0I, num_assets, num_assets), rf_preds_test, c_test, port_oracle)

    final_results.Absolute_spoloss_test = spo_loss(best_B_Absolute, X_test, c_test, port_oracle)
    final_results.Huber_spoloss_test = spo_loss(best_B_Huber, X_test, c_test, port_oracle)

    c_bar_test_preds = zeros(num_assets, n_test)
    for i = 1:n_test
        c_bar_test_preds[:, i] = c_bar_train
    end
    final_results.Baseline_spoloss_test = spo_loss(Matrix(1.0I, num_assets, num_assets), c_bar_test_preds, c_test, port_oracle)

    (z_test, w_test) = oracle_dataset(c_test, port_oracle)
    final_results.zstar_avg_test = mean(z_test)

    return final_results
end

function portfolio_multiple_replications(rng_seed, num_trials,
    num_assets, num_factors, n_train_vec, n_test, p_features,
        n_sigmoid_polydegree_vec, noise_multiplier_tau_vec;
        num_lambda = 10, lambda_max = :practical, lambda_min_ratio = 0.0001, holdout_percent = 0.25,
        different_validation_losses = false, data_type = :poly_kernel)

    Random.seed!(rng_seed)

    big_results_df = make_blank_complete_df()

    envOracle = setup_gurobi_env(method_type = :default, use_time_limit = false)
    envReform = setup_gurobi_env(method_type = :method3, use_time_limit = false)

    for n_train in n_train_vec
        for n_sigmoid_polydegree in n_sigmoid_polydegree_vec
            for noise_multiplier_tau in noise_multiplier_tau_vec
                println("Moving on to n_train = $n_train, n_sigmoid_polydegree = $n_sigmoid_polydegree, noise_multiplier_tau = $noise_multiplier_tau")
                for trial = 1:num_trials
                    println("Current trial is $trial")
                    n_holdout = round(Int, holdout_percent*n_train)

                    current_results = portfolio_replication(num_assets, num_factors, n_train, n_holdout, n_test, p_features,
                        n_sigmoid_polydegree, noise_multiplier_tau;
                        num_lambda = num_lambda, lambda_max = lambda_max, lambda_min_ratio = lambda_min_ratio,
                        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform, different_validation_losses = different_validation_losses,
                        data_type = data_type)

                    current_df_row = build_complete_row(num_assets,
                        n_train, n_holdout, n_test,
                        p_features, n_sigmoid_polydegree, noise_multiplier_tau, current_results)

                    big_results_df = vcat(big_results_df, current_df_row)
                end
            end
        end
    end

    return big_results_df
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
include("..\\libs\\FunctionsRobot.jl")

function zakrasit!(r::Robot)
    sides = moveAndReturnDirections!(r)
    marksArea!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end

ext: jl


text:
include("bandit_utils.jl")
using LinearAlgebra, PyPlot, MultivariateStats, BSON

### load model ###
using BSON: @load
@load "models/bandit.bson" m

##### lets do some evaluation ####

test_as = zeros(T, 0)
for i = 1:100
    L, ys, rews, as, ps, hs, states = run_episode(m, ps = [0.20; 0.80], hidden = true)
    global test_as = [test_as as]
end

### plot 'suboptimal pulls' ###
figure(figsize = (6, 4))
imshow(1 .- test_as', cmap = "Greys", aspect = "auto")
xlabel("trial", fontsize = 12)
ylabel("episode", fontsize = 12)
savefig("figs/bandit_suboptimal_pulls.png", bbox_inches = "tight")
close()

### run for different bandit tasks ###
bandits = [[0.99; 0.01], [0.75; 0.25], [0.5;0.5], [0.25; 0.75], [0.01; 0.99]]
all_hs = zeros(Nhidden, 0)
for ps = bandits
    println(ps)
    L, ys, rews, as, ps, hs = run_episode(m, ps = ps, hidden = true)
    global all_hs = [all_hs hs]
    println(mean(rews))
end

### fit PCA across hidden states and reshape to episodes ###
M = MultivariateStats.fit(PCA, all_hs; maxoutdim=3)
Z = transform(M, all_hs)
Z = (Z .- minimum(Z, dims = 2)) ./ (maximum(Z, dims = 2) - minimum(Z, dims = 2))
zs = [Z[:, (i-1)*T+1:i*T] for i = 1:length(bandits)]

### plot figure ###
figure(figsize = (10, 2.))
axs = 150 .+ (1:5)
for i = 1:length(zs)
    subplot(axs[i])
    scatter(zs[i][1,:], zs[i][2,:], c = 1:T, cmap = "coolwarm", s = 5)
    xticks([])
    yticks([])
    xlim(-0.05, 1.05 )
    ylim(-0.05, 1.05 )
    xlabel("PC 1")
    if i == 1
        ylabel("PC 2")
    end
    title(L"$p_L = $"*string(bandits[i][1]))
end
show()
savefig("figs/bandit_dynamics_pca.png", bbox_inches = "tight")
close()

max_stars_repo_name: KrisJensen/metaRL

max_issues_repo_name: KrisJensen/metaRL

max_forks_repo_name: KrisJensen/metaRL


text:
module Backend

export contacts, chats, update_contact_list, get_key, set_key, get_state, set_state

contacts = Dict()
chats = Dict()

function update_contact_list(message)
    id = message["message"]["from"]["id"]
    if !(id in  keys(contacts))
        contacts[id] = Dict("state"=> 0, "message"=>message)
    end
end

function get_key(id, key)
    return contacts[id][key]
end

function set_key(id, key, value)
    contacts[id][key] = value
end

function get_state(id)
    return get_key(id,"state")
end

function set_state(id, state)
    set_key(id,"state", state)
end

end


text:
function label(colors)

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ThermalCleanEnergy <: ThermalGen
        name::String
        available::Bool
        status::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits::Union{Nothing, Min_Max}
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        operation_cost::OperationalCost
        base_power::Float64
        time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        prime_mover::PrimeMovers
        fuel::ThermalFuels
        services::Vector{Service}
        time_at_status::Float64
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

Data Structure for thermal generation technologies.

# Arguments
- `name::String`
- `available::Bool`
- `status::Bool`
- `bus::Bus`
- `active_power::Float64`, validation range: `active_power_limits`, action if invalid: `warn`
- `reactive_power::Float64`, validation range: `reactive_power_limits`, action if invalid: `warn`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(0, nothing)`, action if invalid: `warn`
- `reactive_power_limits::Union{Nothing, Min_Max}`
- `ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute, validation range: `(0, nothing)`, action if invalid: `error`
- `operation_cost::OperationalCost`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `fuel::ThermalFuels`: Prime mover fuel according to EIA 923
- `services::Vector{Service}`: Services that this device contributes to
- `time_at_status::Float64`
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ThermalCleanEnergy <: ThermalGen
    name::String
    available::Bool
    status::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits::Union{Nothing, Min_Max}
    "ramp up and ramp down limits in MW (in component base per unit) per minute"
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    operation_cost::OperationalCost
    "Base power of the unit in MVA"
    base_power::Float64
    "Minimum up and Minimum down time limits in hours"
    time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    "Prime mover fuel according to EIA 923"
    fuel::ThermalFuels
    "Services that this device contributes to"
    services::Vector{Service}
    time_at_status::Float64
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ThermalCleanEnergy(name, available, status, bus, active_power, reactive_power, rating, active_power_limits, reactive_power_limits, ramp_limits, operation_cost, base_power, time_limits=nothing, prime_mover=PrimeMovers.OT, fuel=ThermalFuels.OTHER, services=Device[], time_at_status=INFINITE_TIME, dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    ThermalCleanEnergy(name, available, status, bus, active_power, reactive_power, rating, active_power_limits, reactive_power_limits, ramp_limits, operation_cost, base_power, time_limits, prime_mover, fuel, services, time_at_status, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function ThermalCleanEnergy(; name, available, status, bus, active_power, reactive_power, rating, active_power_limits, reactive_power_limits, ramp_limits, operation_cost, base_power, time_limits=nothing, prime_mover=PrimeMovers.OT, fuel=ThermalFuels.OTHER, services=Device[], time_at_status=INFINITE_TIME, dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    ThermalCleanEnergy(name, available, status, bus, active_power, reactive_power, rating, active_power_limits, reactive_power_limits, ramp_limits, operation_cost, base_power, time_limits, prime_mover, fuel, services, time_at_status, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function ThermalCleanEnergy(::Nothing)
    ThermalCleanEnergy(;
        name="init",
        available=false,
        status=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        rating=0.0,
        active_power_limits=(min=0.0, max=0.0),
        reactive_power_limits=nothing,
        ramp_limits=nothing,
        operation_cost=ThreePartCost(nothing),
        base_power=0.0,
        time_limits=nothing,
        prime_mover=PrimeMovers.OT,
        fuel=ThermalFuels.OTHER,
        services=Device[],
        time_at_status=INFINITE_TIME,
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`ThermalCleanEnergy`](@ref) `name`."""
PSY.get_name(value::ThermalCleanEnergy) = value.name
"""Get [`ThermalCleanEnergy`](@ref) `available`."""
PSY.get_available(value::ThermalCleanEnergy) = value.available
"""Get [`ThermalCleanEnergy`](@ref) `status`."""
PSY.get_status(value::ThermalCleanEnergy) = value.status
"""Get [`ThermalCleanEnergy`](@ref) `bus`."""
PSY.get_bus(value::ThermalCleanEnergy) = value.bus
"""Get [`ThermalCleanEnergy`](@ref) `active_power`."""
PSY.get_active_power(value::ThermalCleanEnergy) = get_value(value, value.active_power)
"""Get [`ThermalCleanEnergy`](@ref) `reactive_power`."""
PSY.get_reactive_power(value::ThermalCleanEnergy) = get_value(value, value.reactive_power)
"""Get [`ThermalCleanEnergy`](@ref) `rating`."""
PSY.get_rating(value::ThermalCleanEnergy) = get_value(value, value.rating)
"""Get [`ThermalCleanEnergy`](@ref) `active_power_limits`."""
PSY.get_active_power_limits(value::ThermalCleanEnergy) = get_value(value, value.active_power_limits)
"""Get [`ThermalCleanEnergy`](@ref) `reactive_power_limits`."""
PSY.get_reactive_power_limits(value::ThermalCleanEnergy) = get_value(value, value.reactive_power_limits)
"""Get [`ThermalCleanEnergy`](@ref) `ramp_limits`."""
PSY.get_ramp_limits(value::ThermalCleanEnergy) = get_value(value, value.ramp_limits)
"""Get [`ThermalCleanEnergy`](@ref) `operation_cost`."""
PSY.get_operation_cost(value::ThermalCleanEnergy) = value.operation_cost
"""Get [`ThermalCleanEnergy`](@ref) `base_power`."""
PSY.get_base_power(value::ThermalCleanEnergy) = value.base_power
"""Get [`ThermalCleanEnergy`](@ref) `time_limits`."""
PSY.get_time_limits(value::ThermalCleanEnergy) = value.time_limits
"""Get [`ThermalCleanEnergy`](@ref) `prime_mover`."""
PSY.get_prime_mover(value::ThermalCleanEnergy) = value.prime_mover
"""Get [`ThermalCleanEnergy`](@ref) `fuel`."""
PSY.get_fuel(value::ThermalCleanEnergy) = value.fuel
"""Get [`ThermalCleanEnergy`](@ref) `services`."""
PSY.get_services(value::ThermalCleanEnergy) = value.services
"""Get [`ThermalCleanEnergy`](@ref) `time_at_status`."""
PSY.get_time_at_status(value::ThermalCleanEnergy) = value.time_at_status
"""Get [`ThermalCleanEnergy`](@ref) `dynamic_injector`."""
PSY.get_dynamic_injector(value::ThermalCleanEnergy) = value.dynamic_injector
"""Get [`ThermalCleanEnergy`](@ref) `ext`."""
PSY.get_ext(value::ThermalCleanEnergy) = value.ext
"""Get [`ThermalCleanEnergy`](@ref) `time_series_container`."""
PSY.get_time_series_container(value::ThermalCleanEnergy) = value.time_series_container
"""Get [`ThermalCleanEnergy`](@ref) `internal`."""
PSY.get_internal(value::ThermalCleanEnergy) = value.internal

"""Set [`ThermalCleanEnergy`](@ref) `name`."""
PSY.set_name!(value::ThermalCleanEnergy, val) = value.name = val
"""Set [`ThermalCleanEnergy`](@ref) `available`."""
PSY.set_available!(value::ThermalCleanEnergy, val) = value.available = val
"""Set [`ThermalCleanEnergy`](@ref) `status`."""
PSY.set_status!(value::ThermalCleanEnergy, val) = value.status = val
"""Set [`ThermalCleanEnergy`](@ref) `bus`."""
PSY.set_bus!(value::ThermalCleanEnergy, val) = value.bus = val
"""Set [`ThermalCleanEnergy`](@ref) `active_power`."""
PSY.set_active_power!(value::ThermalCleanEnergy, val) = value.active_power = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `reactive_power`."""
PSY.set_reactive_power!(value::ThermalCleanEnergy, val) = value.reactive_power = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `rating`."""
PSY.set_rating!(value::ThermalCleanEnergy, val) = value.rating = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `active_power_limits`."""
PSY.set_active_power_limits!(value::ThermalCleanEnergy, val) = value.active_power_limits = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `reactive_power_limits`."""
PSY.set_reactive_power_limits!(value::ThermalCleanEnergy, val) = value.reactive_power_limits = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `ramp_limits`."""
PSY.set_ramp_limits!(value::ThermalCleanEnergy, val) = value.ramp_limits = set_value(value, val)
"""Set [`ThermalCleanEnergy`](@ref) `operation_cost`."""
PSY.set_operation_cost!(value::ThermalCleanEnergy, val) = value.operation_cost = val
"""Set [`ThermalCleanEnergy`](@ref) `base_power`."""
PSY.set_base_power!(value::ThermalCleanEnergy, val) = value.base_power = val
"""Set [`ThermalCleanEnergy`](@ref) `time_limits`."""
PSY.set_time_limits!(value::ThermalCleanEnergy, val) = value.time_limits = val
"""Set [`ThermalCleanEnergy`](@ref) `prime_mover`."""
PSY.set_prime_mover!(value::ThermalCleanEnergy, val) = value.prime_mover = val
"""Set [`ThermalCleanEnergy`](@ref) `fuel`."""
PSY.set_fuel!(value::ThermalCleanEnergy, val) = value.fuel = val
"""Set [`ThermalCleanEnergy`](@ref) `services`."""
PSY.set_services!(value::ThermalCleanEnergy, val) = value.services = val
"""Set [`ThermalCleanEnergy`](@ref) `time_at_status`."""
PSY.set_time_at_status!(value::ThermalCleanEnergy, val) = value.time_at_status = val
"""Set [`ThermalCleanEnergy`](@ref) `ext`."""
PSY.set_ext!(value::ThermalCleanEnergy, val) = value.ext = val
"""Set [`ThermalCleanEnergy`](@ref) `time_series_container`."""
PSY.set_time_series_container!(value::ThermalCleanEnergy, val) = value.time_series_container = val

max_stars_repo_name: sourabhdalvi/PowerSystemExtensions

max_issues_repo_name: sourabhdalvi/PowerSystemExtensions

max_forks_repo_name: sourabhdalvi/PowerSystemExtensions


text:
function FieldAbsolutePVCoordinates(arg0::FieldAbsoluteDate, arg1::FieldAbsolutePVCoordinates, arg2::FieldAbsolutePVCoordinates)
    return FieldAbsolutePVCoordinates((FieldAbsoluteDate, FieldAbsolutePVCoordinates, FieldAbsolutePVCoordinates), arg0, arg1, arg2)
end

function FieldAbsolutePVCoordinates(arg0::FieldAbsoluteDate, arg1::RealFieldElement, arg2::FieldAbsolutePVCoordinates)
    return FieldAbsolutePVCoordinates((FieldAbsoluteDate, RealFieldElement, FieldAbsolutePVCoordinates), arg0, arg1, arg2)
end

function FieldAbsolutePVCoordinates(arg0::FieldAbsoluteDate, arg1::RealFieldElement, arg2::FieldAbsolutePVCoordinates, arg3::RealFieldElement, arg4::FieldAbsolutePVCoordinates)
    return FieldAbsolutePVCoordinates((FieldAbsoluteDate, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4)
end

function FieldAbsolutePVCoordinates(arg0::FieldAbsoluteDate, arg1::RealFieldElement, arg2::FieldAbsolutePVCoordinates, arg3::RealFieldElement, arg4::FieldAbsolutePVCoordinates, arg5::RealFieldElement, arg6::FieldAbsolutePVCoordinates)
    return FieldAbsolutePVCoordinates((FieldAbsoluteDate, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function FieldAbsolutePVCoordinates(arg0::FieldAbsoluteDate, arg1::RealFieldElement, arg2::FieldAbsolutePVCoordinates, arg3::RealFieldElement, arg4::FieldAbsolutePVCoordinates, arg5::RealFieldElement, arg6::FieldAbsolutePVCoordinates, arg7::RealFieldElement, arg8::FieldAbsolutePVCoordinates)
    return FieldAbsolutePVCoordinates((FieldAbsoluteDate, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates, RealFieldElement, FieldAbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function FieldAbsolutePVCoordinates(arg0::Frame, arg1::FieldAbsoluteDate, arg2::FieldPVCoordinates)
    return FieldAbsolutePVCoordinates((Frame, FieldAbsoluteDate, FieldPVCoordinates), arg0, arg1, arg2)
end

function FieldAbsolutePVCoordinates(arg0::Frame, arg1::FieldAbsoluteDate, arg2::FieldVector3D)
    return FieldAbsolutePVCoordinates((Frame, FieldAbsoluteDate, FieldVector3D), arg0, arg1, arg2)
end

function FieldAbsolutePVCoordinates(arg0::Frame, arg1::FieldAbsoluteDate, arg2::FieldVector3D, arg3::FieldVector3D)
    return FieldAbsolutePVCoordinates((Frame, FieldAbsoluteDate, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3)
end

function FieldAbsolutePVCoordinates(arg0::Frame, arg1::FieldAbsoluteDate, arg2::FieldVector3D, arg3::FieldVector3D, arg4::FieldVector3D)
    return FieldAbsolutePVCoordinates((Frame, FieldAbsoluteDate, FieldVector3D, FieldVector3D, FieldVector3D), arg0, arg1, arg2, arg3, arg4)
end

function FieldAbsolutePVCoordinates(arg0::Frame, arg1::TimeStampedFieldPVCoordinates)
    return FieldAbsolutePVCoordinates((Frame, TimeStampedFieldPVCoordinates), arg0, arg1)
end

function cross_product(obj::FieldPVCoordinates, arg0::FieldPVCoordinates)
    return jcall(obj, "crossProduct", FieldPVCoordinates, (FieldPVCoordinates,), arg0)
end

function estimate_velocity(::Type{FieldPVCoordinates}, arg0::FieldVector3D, arg1::FieldVector3D, arg2::jdouble)
    return jcall(FieldPVCoordinates, "estimateVelocity", FieldVector3D, (FieldVector3D, FieldVector3D, jdouble), arg0, arg1, arg2)
end

function get_acceleration(obj::FieldPVCoordinates)
    return jcall(obj, "getAcceleration", FieldVector3D, ())
end

function get_angular_velocity(obj::FieldPVCoordinates)
    return jcall(obj, "getAngularVelocity", FieldVector3D, ())
end

function get_date(obj::TimeStampedFieldPVCoordinates)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_frame(obj::FieldAbsolutePVCoordinates)
    return jcall(obj, "getFrame", Frame, ())
end

function get_momentum(obj::FieldPVCoordinates)
    return jcall(obj, "getMomentum", FieldVector3D, ())
end

function get_position(obj::FieldPVCoordinates)
    return jcall(obj, "getPosition", FieldVector3D, ())
end

function get_pv_coordinates(obj::FieldAbsolutePVCoordinates)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, ())
end

function get_pv_coordinates(obj::FieldAbsolutePVCoordinates, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::FieldAbsolutePVCoordinates, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (Frame,), arg0)
end

function get_velocity(obj::FieldPVCoordinates)
    return jcall(obj, "getVelocity", FieldVector3D, ())
end

function get_zero(::Type{FieldPVCoordinates}, arg0::Field)
    return jcall(FieldPVCoordinates, "getZero", FieldPVCoordinates, (Field,), arg0)
end

function interpolate(::Type{FieldAbsolutePVCoordinates}, arg0::Frame, arg1::FieldAbsoluteDate, arg2::CartesianDerivativesFilter, arg3::Stream)
    return jcall(FieldAbsolutePVCoordinates, "interpolate", FieldAbsolutePVCoordinates, (Frame, FieldAbsoluteDate, CartesianDerivativesFilter, Stream), arg0, arg1, arg2, arg3)
end

function interpolate(::Type{TimeStampedFieldPVCoordinates}, arg0::FieldAbsoluteDate, arg1::CartesianDerivativesFilter, arg2::Collection)
    return jcall(TimeStampedFieldPVCoordinates, "interpolate", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, CartesianDerivativesFilter, Collection), arg0, arg1, arg2)
end

function interpolate(::Type{TimeStampedFieldPVCoordinates}, arg0::FieldAbsoluteDate, arg1::CartesianDerivativesFilter, arg2::Stream)
    return jcall(TimeStampedFieldPVCoordinates, "interpolate", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, CartesianDerivativesFilter, Stream), arg0, arg1, arg2)
end

function interpolate(obj::FieldAbsolutePVCoordinates, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldAbsolutePVCoordinates, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function negate(obj::FieldPVCoordinates)
    return jcall(obj, "negate", FieldPVCoordinates, ())
end

function normalize(obj::FieldPVCoordinates)
    return jcall(obj, "normalize", FieldPVCoordinates, ())
end

function shifted_by(obj::FieldAbsolutePVCoordinates, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldAbsolutePVCoordinates, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldAbsolutePVCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldAbsolutePVCoordinates, (jdouble,), arg0)
end

function to_absolute_pv_coordinates(obj::FieldAbsolutePVCoordinates)
    return jcall(obj, "toAbsolutePVCoordinates", AbsolutePVCoordinates, ())
end

function to_derivative_structure_pv(obj::FieldPVCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructurePV", FieldPVCoordinates, (jint,), arg0)
end

function to_derivative_structure_vector(obj::FieldPVCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructureVector", FieldVector3D, (jint,), arg0)
end

function to_pv_coordinates(obj::FieldPVCoordinates)
    return jcall(obj, "toPVCoordinates", PVCoordinates, ())
end

function to_string(obj::TimeStampedFieldPVCoordinates)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::TimeStampedFieldPVCoordinates, arg0::TimeScale)
    return jcall(obj, "toString", JString, (TimeScale,), arg0)
end

function to_taylor_provider(obj::FieldAbsolutePVCoordinates)
    return jcall(obj, "toTaylorProvider", FieldPVCoordinatesProvider, ())
end

function to_time_stamped_pv_coordinates(obj::TimeStampedFieldPVCoordinates)
    return jcall(obj, "toTimeStampedPVCoordinates", TimeStampedPVCoordinates, ())
end

function to_univariate_derivative1_vector(obj::FieldPVCoordinates)
    return jcall(obj, "toUnivariateDerivative1Vector", FieldVector3D, ())
end

function to_univariate_derivative1pv(obj::FieldPVCoordinates)
    return jcall(obj, "toUnivariateDerivative1PV", FieldPVCoordinates, ())
end

function to_univariate_derivative2_vector(obj::FieldPVCoordinates)
    return jcall(obj, "toUnivariateDerivative2Vector", FieldVector3D, ())
end

function to_univariate_derivative2pv(obj::FieldPVCoordinates)
    return jcall(obj, "toUnivariateDerivative2PV", FieldPVCoordinates, ())
end


text:
using PkgBenchmark

group_target = PkgBenchmark.readresults(joinpath(@__DIR__, "result-target.json"))
group_baseline = PkgBenchmark.readresults(joinpath(@__DIR__, "result-baseline.json"))

judgement = judge(group_target, group_baseline)

using Markdown
let io = IOBuffer()
    PkgBenchmark.export_markdown(io, judgement)
    seekstart(io)
    md = Markdown.parse(io)
    display(md)
end


text:
# initialize the attractor
n = 100
n = 10000
s = [10., 28., 8/3]
u0 = rand(3,1)
u_trj = lorenz63(u0, s, n)
# initialize a 3D plot with 1 empty series
plt = path3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
                xlab = "x", ylab = "y", zlab = "z",
                title = "Lorenz Attractor", marker = 1)

# build an animated gif, saving every 10th frame
u_trj = permutedims(u_trj,[3,2,1])
@gif for i=1:n
		push!(plt, u_trj[i,1,:], u_trj[i,2,:], u_trj[i,3,:])
end every 10

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
"""
    AirBase(id::String)

A wrapper for an [Airtable Base](https://support.airtable.com/hc/en-us/articles/202576419-Introduction-to-Airtable-bases),
containing its unique identifier.
The ID can be identified from the URL of the base (the part right after `airtable.com`),
or by clicking `HELP -> API documentation` from within your base.
"""
struct AirBase
    id::String
end

"""
    id(::AirBase)

Accessor function for a base identifier.
"""
id(bs::AirBase) = bs.id

"""
    path(::AirBase)

Accessor function for the API path to the base.
This is just the base ID preceded by the API version number (only `v0` for now).

```jldoctest
julia> b = AirBase("appphImnhJO8AXmmo");

julia> path(b)
"/v0/appphImnhJO8AXmmo"
```
"""
path(bs::AirBase) = joinpath("/", API_VERSION, id(bs))
Base.show(io::IO, bs::AirBase) = println(io, "Airtable Base '$(id(bs))'")

"""
    AirTable(id::String, ::AirBase)

A wrapper for an [Airtable Table](https://support.airtable.com/hc/en-us/articles/360021333094-Getting-started-tables-records-and-fields),
containing its unique identifier or name, and parent [`AirBase`](@ref).
"""
struct AirTable
    id::String
    base::AirBase
end

"""
    id(::AirTable)

Accessor function for a table identifier.
"""
id(tab::AirTable) = tab.id

"""
    base(::AirTable)

Accessor function for the parent base of a table
"""
base(tab::AirTable) = tab.base

"""
    path(::AirTable)

Accessor function for the API path to the table.
This is just the table ID preceded by the [path to the base](@ref `path(::AirBase)`)

```jldoctest
julia> b = AirBase("appphImnhJO8AXmmo");

julia> tab = AirTable("Table 1", b);

julia> path(tab)
"/v0/appphImnhJO8AXmmo/Table 1"
```
"""
path(tab::AirTable) = joinpath(path(base(tab)), id(tab))
Base.show(io::IO, tab::AirTable) = print(io, "AirTable(\"$(id(tab))\")")

"""
    AirRecord(id::String, table::AirTable, fields::NamedTuple)

A wrapper for an [Airtable Record](https://support.airtable.com/hc/en-us/articles/360021333094-Getting-started-tables-records-and-fields),
containing its unique identifier, parent [`AirTable`](@ref),
and values for any stored fields in a NamedTuple.

Typically, you won't crete these on your own, but they will be returned from API queries.

Field values can be accessed using [`getindex`](@ref `Base.getindex(::AirRecord, ::Symbol)`).
"""
struct AirRecord
    id::String
    table::AirTable
    fields::NamedTuple

    function AirRecord(id, table, fields)
        new(id, table, NamedTuple(fields))
    end
end

AirRecord(; id, table, fields=Dict()) = AirRecord(id, table, fields)

"""
    id(::AirRecord)

Accessor function for a record identifier.
"""
id(rec::AirRecord) = rec.id

"""
    table(::AirRecord)

Accessor function the the parent [`AirTable`](@ref) of a record.
"""
table(rec::AirRecord) = rec.table

"""
    base(::AirRecord)

Accessor function the the parent [`AirBase`](@ref) of a record.
"""
base(rec::AirRecord) = base(table(rec))

"""
    base(::AirRecord)

Accessor function the the fields within an [`AirRecord`](@ref).
Individual values can also be accessed with [`getindex`](@ref `Base.getindex(::AirRecord, ::Symbol)`).
"""
fields(rec::AirRecord) = rec.fields

"""
    path(::AirTable)

Accessor function for the API path to the table.
This is just the table ID preceded by the [path to the table](@ref `path(::AirTable)`)
"""
path(rec::AirRecord) = joinpath(path(table(rec)), id(rec))

JSON3.write(rec::AirRecord) = string("""{ "id": "$(id(rec))", "fields": """, JSON3.write(fields(rec)), "}")
AirRecord(rec::JSON3.Object, tab::AirTable) = AirRecord(rec.id, tab, rec.fields)

"""
    Airtable.query(cred::Credential, baseid, tablename; query_kwargs...)

Shorthand for a "GET" request that handles continuation and rate-limiting.

The Airtable API will return a maximum of 100 records per requests,
and only allows 5 requests / sec.
This function uses the `offset` field returned as part of a requst
that does not contain all possible records to make additional requests
after pausing 0.21 seconds in between.

Required arguments:

- `cred`: an [`Airtable.Credential`](@ref) containing your API key
- `baseid`: the endpoint of your Airtable base. See https://airtable.com/api for details
- `tablename`: the name of the table in your base (eg `"Table 1"`)

Query parameters are in the form of keyword arguments,
eg `filterByFormla = "NOT({Name} = '')", maxRecords=2`.
See Airtable API reference for more information.

If you know the exact record id, pass that as a fourth positional argument
"""
function query(cred::Credential, baseid, tablename; query_kwargs...)
    tab = AirTable(tablename, AirBase(baseid))
    resp = get(cred, path(tab); query_kwargs...)
    records = map(rec-> AirRecord(rec.id, tab, rec.fields), resp.records)

    if haskey(resp, :offset)
        prog = ProgressUnknown("Making additional requests:", spinner=true)
        while haskey(resp, :offset)
            ProgressMeter.next!(prog)
            resp = get(cred, path(tab); offset=resp.offset, query_kwargs...)
            append!(records, map(rec-> AirRecord(rec.id, tab, rec.fields), resp.records))
        end
        ProgressMeter.finish!(prog)
    end

    return records
end

query(cred::Credential, tab::AirTable; query_kwargs...) = query(cred, id(base(tab)), id(tab); query_kwargs...)
query(tab::AirTable; query_kwargs...) = query(Credential(), tab; query_kwargs...)

get(cred::Credential, rec::AirRecord) = AirRecord(get(cred, path(rec)), table(rec))
get(rec::AirRecord) = get(Credential(), rec)

_extract_records(tab::AirTable, resp) = map(rec-> AirRecord(rec.id, tab, rec.fields), resp.records)
_extract_record(tab::AirTable, rec) = AirRecord(rec.id, tab, rec.fields)

function post!(cred::Credential, tab::AirTable, rec::AirRecord)
    resp = post!(cred, path(tab), ["Content-Type" => "application/json"], JSON3.write(rec))
    return _extract_record(tab, resp)
end

function post!(cred::Credential, tab::AirTable, recs::IO)
    resp = post!(cred, path(tab), ["Content-Type" => "application/json"], JSON3.write(JSON3.read(recs)))
    return _extract_records(tab, resp)
end

function post!(cred::Credential, tab::AirTable, rec::NamedTuple)
    resp = post!(cred, path(tab), ["Content-Type" => "application/json"], JSON3.write((; fields=rec)))
    return _extract_record(tab, resp)
end

function post!(cred::Credential, tab::AirTable, recs::Vector{<:NamedTuple})
    resps = AirRecord[]
    for recpart in Iterators.partition(recs, 10)
        topost = (; records = [(; fields = nt) for nt in recpart])
        resp = post!(cred, path(tab), ["Content-Type" => "application/json"], JSON3.write(topost))
        append!(resps, _extract_records(tab, resp))
    end
    return resps
end

post!(tab::AirTable, rec::AirRecord)           = post!(Credential(), tab, rec)
post!(tab::AirTable, recs::IO)                 = post!(Credential(), tab, recs)
post!(tab::AirTable, rec::NamedTuple)          = post!(Credential(), tab, rec)
post!(tab::AirTable, recs::Vector{<:NamedTuple}) = post!(Credential(), tab, recs)

delete!(cred::Credential, rec::AirRecord) = delete!(cred, path(rec))
delete!(rec::AirRecord) = delete!(Credential(), rec)

patch!(cred::Credential, rec::AirRecord) = _extract_record(table(rec), patch!(cred, path(rec), ["Content-Type" => "application/json"], JSON3.write(rec)))
patch!(cred::Credential, rec::AirRecord, fields::NamedTuple) = _extract_record(table(rec), patch!(cred, path(rec), ["Content-Type" => "application/json"],  JSON3.write((; fields))))

function patch!(cred::Credential, tab::AirTable, recs::Vector{<:AirRecord})
    resps = AirRecord[]
    for recpart in Iterators.partition(recs, 10)
        resp = patch!(cred, path(tab), ["Content-Type" => "application/json"], JSON3.write((; records = [
            (; id=id(rec), fields=fields(rec)) for rec in recpart
        ])))
        append!(resps, _extract_records(tab, resp))
    end
    return resps
end

function patch!(cred::Credential, tab::AirTable, recs::Vector{<:AirRecord}, fields::Vector{<:NamedTuple})
    length(recs) == length(fields) || throw(ArgumentError("Lengths of records vector and fields vector must be the same"))
    patch!(cred, tab, [AirRecord(Airtable.id(rec), tab, fs) for (rec, fs) in zip(recs, fields)])
end

patch!(rec::AirRecord) = patch!(Credential(), rec)
patch!(rec::AirRecord, fields::NamedTuple) = patch!(Credential(), rec, fields)
patch!(tab::AirTable, recs::Vector{<:AirRecord}) = patch!(Credential(), tab, recs)
patch!(tab::AirTable, recs::Vector{<:AirRecord}, fields::Vector{<:NamedTuple}) = patch!(Credential(), tab, recs, fields)

Base.getindex(rec::AirRecord, k::Symbol) = fields(rec)[k]
Base.keys(rec::AirRecord) = keys(fields(rec))
Base.haskey(rec::AirRecord, key::Symbol) = key ∈ keys(rec)
Base.get(rec::AirRecord, key, default) = haskey(rec, key) ? rec[key] : default


text:
using WAV
using DSP

f, fs, nbits, opt = wavread("audio.wav")

f_down = fill(0.0, round(Int, floor(length(f)/32)))
for i = 1:32:length(f)
    f_down[round(Int, floor(i / 32))] = f[i]
end

rate = int(fs)
impulse_train = fill(0.0, rate * 3)
impulse_train[1] = 1.0
impulse_train[rate + 1] = 1.0
impulse_train[2 * rate + 1] = 1.0
sum(conv(f_down[1:rate*3], impulse_train))

max_stars_repo_name: broluwo/ciett

max_issues_repo_name: broluwo/ciett

max_forks_repo_name: broluwo/ciett


text:
#
# Qubit site type
#

# Define Qubit space in terms of
# S=1/2 space, but use different
# defaults for QN names

"""
    space(::SiteType"Qubit";
          conserve_qns = false,
          conserve_parity = conserve_qns,
          conserve_number = false,
          qnname_parity = "Parity",
          qnname_number = "Number")

Create the Hilbert space for a site of type "Qubit".

Optionally specify the conserved symmetries and their quantum number labels.
"""
space(::SiteType"Qubit";
      conserve_qns=false,
      conserve_parity=conserve_qns,
      conserve_number=false,
      qnname_parity="Parity",
      qnname_number="Number") =
  space(SiteType("S=1/2");
        conserve_qns=conserve_qns,
        conserve_sz=conserve_number,
        conserve_szparity=conserve_parity,
        qnname_sz=qnname_number,
        qnname_szparity=qnname_parity)

state(::SiteType"Qubit", ::StateName"0") = 1

state(::SiteType"Qubit", ::StateName"1") = 2

# Use S=1/2 definition of any operators
# called using Qubit SiteType
op(o::OpName, ::SiteType"Qubit";kwargs...) =
  op(o, SiteType("S=1/2");kwargs...)

#
# 1-Qubit gates
#

op(::OpName"σx",t::SiteType"S=1/2") =
  op("X",t)

op(::OpName"σ1",t::SiteType"S=1/2") =
  op("X",t)

op(::OpName"√NOT",::SiteType"S=1/2") =
  [(1+im)/2 (1-im)/2
   (1-im)/2 (1+im)/2]

op(::OpName"√X",t::SiteType"S=1/2") =
  op("√NOT",t)

op(::OpName"σy",t::SiteType"S=1/2") =
  op("Y",t)

op(::OpName"σ2",t::SiteType"S=1/2") =
  op("Y",t)

op(::OpName"iσy",t::SiteType"S=1/2") =
  op("iY",t)

op(::OpName"iσ2",t::SiteType"S=1/2") =
  op("iY",t)

op(::OpName"σz",t::SiteType"S=1/2") =
  op("Z",t)

op(::OpName"σ3",t::SiteType"S=1/2") =
  op("Z",t)

op(::OpName"H",::SiteType"S=1/2") =
  [1/sqrt(2) 1/sqrt(2)
   1/sqrt(2) -1/sqrt(2)]

# Rϕ with ϕ = π/2
op(::OpName"Phase",::SiteType"S=1/2") =
  [1  0
   0 im]

op(::OpName"P",t::SiteType"S=1/2") =
  op("Phase",t)

op(::OpName"S",t::SiteType"S=1/2") =
  op("Phase",t)

## Rϕ with ϕ = π/4
op(::OpName"π/8",::SiteType"S=1/2") =
  [1  0
   0  1/sqrt(2) + im/sqrt(2)]

op(::OpName"T",t::SiteType"S=1/2") =
  op("π/8",t)

# Rotation around X-axis
op(::OpName"Rx",::SiteType"S=1/2"; θ::Number) =
  [    cos(θ/2)  -im*sin(θ/2)
   -im*sin(θ/2)      cos(θ/2)]

# Rotation around Y-axis
op(::OpName"Ry",::SiteType"S=1/2"; θ::Number) =
  [cos(θ/2) -sin(θ/2)
   sin(θ/2)  cos(θ/2)]

# Rotation around Z-axis
op(::OpName"Rz",::SiteType"S=1/2"; ϕ::Number) =
  [1         0
   0 exp(im*ϕ)]

# Rotation around generic axis n̂
op(::OpName"Rn",::SiteType"S=1/2";
    θ::Real, ϕ::Real, λ::Real) =
  [          cos(θ/2)    -exp(im*λ)*sin(θ/2)
   exp(im*ϕ)*sin(θ/2) exp(im*(ϕ+λ))*cos(θ/2)]

op(::OpName"Rn̂",t::SiteType"S=1/2"; kwargs...) =
  op("Rn",t; kwargs...)

#
# 2-Qubit gates
#

op(::OpName"CNOT",::SiteType"S=1/2") =
  [1 0 0 0
   0 1 0 0
   0 0 0 1
   0 0 1 0]

op(::OpName"CX",t::SiteType"S=1/2") =
  op("CNOT",t)

op(::OpName"CY",::SiteType"S=1/2") =
  [1 0  0   0
   0 1  0   0
   0 0  0 -im
   0 0 im   0]

op(::OpName"CZ",::SiteType"S=1/2") =
  [1 0 0  0
   0 1 0  0
   0 0 1  0
   0 0 0 -1]

op(::OpName"SWAP",::SiteType"S=1/2") =
  [1 0 0 0
   0 0 1 0
   0 1 0 0
   0 0 0 1]

op(::OpName"√SWAP",::SiteType"S=1/2") =
  [1        0        0 0
   0 (1+im)/2 (1-im)/2 0
   0 (1-im)/2 (1+im)/2 0
   0        0        0 1]

# Ising (XX) coupling gate
op(::OpName"XX",::SiteType"S=1/2"; ϕ::Number) =
  [    cos(ϕ)          0          0 -im*sin(ϕ)
            0     cos(ϕ) -im*sin(ϕ)          0
            0 -im*sin(ϕ)     cos(ϕ)          0
   -im*sin(ϕ)          0          0     cos(ϕ)]

# Ising (YY) coupling gate
op(::OpName"YY",::SiteType"S=1/2"; ϕ::Number) =
  [    cos(ϕ)          0          0  im*sin(ϕ)
            0     cos(ϕ) -im*sin(ϕ)          0
            0 -im*sin(ϕ)     cos(ϕ)          0
    im*sin(ϕ)          0          0     cos(ϕ)]

# Ising (ZZ) coupling gate
op(::OpName"ZZ",::SiteType"S=1/2"; ϕ::Number) =
  [    exp(-im*ϕ)         0         0          0
            0     exp(im*ϕ)         0          0
            0             0 exp(im*ϕ)          0
            0             0         0 exp(-im*ϕ)]

#
# 3-Qubit gates
#

op(::OpName"Toffoli",::SiteType"S=1/2") =
  [1 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0
   0 0 1 0 0 0 0 0
   0 0 0 1 0 0 0 0
   0 0 0 0 1 0 0 0
   0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 1
   0 0 0 0 0 0 1 0]

op(::OpName"CCNOT",t::SiteType"S=1/2") =
  op("Toffoli",t)

op(::OpName"CCX",t::SiteType"S=1/2") =
  op("Toffoli",t)

op(::OpName"TOFF",t::SiteType"S=1/2") =
  op("Toffoli",t)

op(::OpName"Fredkin",::SiteType"S=1/2") =
  [1 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0
   0 0 1 0 0 0 0 0
   0 0 0 1 0 0 0 0
   0 0 0 0 1 0 0 0
   0 0 0 0 0 0 1 0
   0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 1]

op(::OpName"CSWAP",t::SiteType"S=1/2") =
  op("Fredkin",t)

op(::OpName"CS",t::SiteType"S=1/2") =
  op("Fredkin",t)

#
# 4-Qubit gates
#

op(::OpName"CCCNOT",::SiteType"S=1/2") =
  [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]


text:
# tests with change direction update
# if done twice, does it leave the same network?
# Claudia March 2015

# Case G
include("../examples/case_g_example.jl");
net0 = deepcopy(net);
printEdges(net)
printNodes(net)
node = searchHybridNode(net);
node[1].number
success,newnode = changeDirectionUpdate!(net,node[1]);
printEdges(net)

success,node = changeDirectionUpdate!(net,newnode);
printEdges(net)
printEdges(net0)
# not the same, bad diamond I found in between

# Case H
include("../examples/case_h_example.jl");
net0 = deepcopy(net);
printEdges(net)
printNodes(net)
node = searchHybridNode(net);
node[1].number
success,newnode = changeDirectionUpdate!(net,node[1]);
printEdges(net)

success,node = changeDirectionUpdate!(net,newnode);
printEdges(net)
printEdges(net0)
# not the same, bad diamond II found in between

# Case J
include("../examples/case_j_example.jl");
net0 = deepcopy(net);
printEdges(net)
printNodes(net)
node = searchHybridNode(net);
node[1].number
success,newnode = changeDirectionUpdate!(net,node[1]);
printEdges(net)

success,node = changeDirectionUpdate!(net,newnode);
printEdges(net)
printEdges(net0)


text:
module Day03
export day03

const TESTINPUT = """
00100
11110
10110
10111
10101
01111
00111
11100
10000
11001
00010
01010
"""

function parseLine(line)
    # work around string being a scalar for broadcasting
    return parse.(Bool, Iterators.take(line, length(line)))
end

function readInput(inputfile="inputs/input03.txt")
    return vcat((parseLine(line)' for line in eachline(inputfile))...)
end

function day03(inputdata = readInput("inputs/input03.txt"))
    return (day03Part1(inputdata), day03Part2(inputdata))
end

function day03Part1(values)
    resultbits = sum(values; dims = 1) .> div(size(values, 1), 2)
    γvalue = bits_to_int(resultbits)
    ϵvalue = bits_to_int(.!(resultbits))
    return γvalue*ϵvalue
end

function day03Part2(values)
    oxygen_rating = compute_diagnostic(values, >=)
    c02scrub_rating = compute_diagnostic(values, <)
    return oxygen_rating*c02scrub_rating
end

function compute_diagnostic(values, cmp)
    mask = fill!(BitArray(undef, size(values, 1)), true)
    col = 1
    while col <= size(values, 2) && sum(mask) > 1
        updatemask!(mask, values[:, col], cmp)
        col += 1
    end
    @assert sum(mask) == 1
    row = argmax(mask)
    return bits_to_int(values[row,:])
end

function updatemask!(mask, colvalues, cmp)
    colsum = sum(colvalues[mask])
    selectVal = cmp(2*colsum, sum(mask)) ? 1 : 0
    mask .&= (colvalues .== selectVal)
end

bits_to_int(something) = bits_to_int(convert(BitArray, something))
bits_to_int(bitarr::BitArray) = reinterpret(Int64, reverse(bitarr).chunks[1]) # only works for very small bitarrays

end #module


text:
#! /usr/bin/env julia

using MeshSlicer
using Base.Test
using ImmutableArrays

# Since we use relative paths, you need to test from MeshSlicer/test

# setup test handlers
test_handler(r::Test.Success) = println("Test success: $(r.expr)")
test_handler(r::Test.Failure) = error("Test fail: $(r.expr)")
test_handler(r::Test.Error) = rethrow(r)

# getmesh/getface test
println("Testing PolygonMesh...")

# get the first face from cube.stl
testFace = Face(Vector3[Vector3(10.0,0.0,10.0),
                      Vector3(10.0,0.0,0.0),
                      Vector3(10.0,10.0,10.0)],
                      Vector3(1.0,0.0,0.0))

testBounds = Bounds{Float64}(10.0,10.0,10.0,0,0,0)

resultmesh = PolygonMesh("./data/cube.stl")

Test.with_handler(test_handler) do
    @test resultmesh.faces[end].vertices == testFace.vertices
    @test resultmesh.faces[end].normal == testFace.normal
end

println("Testing Bounds...")
Test.with_handler(test_handler) do
    @test resultmesh.bounds.xmax == testBounds.xmax
    @test resultmesh.bounds.ymax == testBounds.ymax
    @test resultmesh.bounds.zmax == testBounds.zmax
    @test resultmesh.bounds.xmin == testBounds.xmin
    @test resultmesh.bounds.ymin == testBounds.ymin
    @test resultmesh.bounds.zmin == testBounds.zmin
    @test resultmesh.bounds == testBounds
end

# Test bounds on translated cube
# translate([1,2,3])cube([4,5,6]);
testBounds2 = Bounds(4.0,6.0,9.0,1.0,2.0,4.0)
updateFace = Face(Vector3[[5.0,3.0,4.0],
                      [2.0,7.0,4.0],
                      [2.0,3.0,3.0]],
                      Vector3(-1.0,0.0,0.0)) # contains correct bounds for cube
update!(testBounds2, updateFace)
resultmesh2 = PolygonMesh("./data/translated_cube.stl")
Test.with_handler(test_handler) do
    @test resultmesh2.bounds.xmax == testBounds2.xmax
    @test resultmesh2.bounds.ymax == testBounds2.ymax
    @test resultmesh2.bounds.zmax == testBounds2.zmax
    @test resultmesh2.bounds.xmin == testBounds2.xmin
    @test resultmesh2.bounds.ymin == testBounds2.ymin
    @test resultmesh2.bounds.zmin == testBounds2.zmin
    @test resultmesh2.bounds == testBounds2
end

# Test patio detection
println("Tesing patio detection...")
Test.with_handler(test_handler) do
    @test resultmesh2.patios == [3.0, 9.0]
    @test resultmesh.patios == [0.0, 10.0]
end

# Test LineSegment
testpoint1 = Vector3[[0.0,0.0,0.0],[0.0,1.0,0.0],[0.0,1.0,1.0]]
normal = Vector3(1.0, 1.0, 1.0)
l1 = LineSegment(testpoint1[3],testpoint1[1],testpoint1[2], 0.5, normal)
l2 = LineSegment(testpoint1[3],testpoint1[2],testpoint1[1], 0.5, normal)
l3 = LineSegment(testpoint1[3],testpoint1[2],testpoint1[1], 0.0, normal)

testpoint2 = Vector3[[0.0,0.0,0.0],[1.0,1.0,1.0],[0.0,1.0,0.0]]
l4 = LineSegment(testpoint2[2],testpoint2[1],testpoint2[3], 0.5, normal)
l5 = LineSegment(testpoint2[2],testpoint2[3],testpoint2[1], 0.5, normal)

l6 = LineSegment(testFace, 5.0)
l7 = LineSegment(testpoint2[2],testpoint2[3],testpoint2[1], 1.0, normal)

Test.with_handler(test_handler) do
    @test l1.start == Vector2(0.0, 0.5)
    @test l1.finish == Vector2(0.0, 1.0)
    @test l2.start == Vector2(0.0, 1.0)
    @test l2.finish == Vector2(0.0, 0.5)
    @test l3.start == Vector2(0.0, 1.0)
    @test l3.finish == Vector2(0.0, 0.0)
    @test l4.start == Vector2(0.5, 0.5)
    @test l4.finish == Vector2(0.5, 1.0)
    @test l5.start == Vector2(0.5, 1.0)
    @test l6.start == Vector2(10.0, 0.0)
    @test l6.finish == Vector2(10.0, 5.0)
    @test l6.normal == Vector3(1.0,0.0,0.0)
    @test l7.start == Vector2(1.0, 1.0)
    @test l7.finish == Vector2(1.0, 1.0)
end

# test binary STL
println("Testing Binary STL import...")
binarymesh = PolygonMesh("./data/cube_binary.stl")
Test.with_handler(test_handler) do
    @test binarymesh.bounds == Bounds(10.0,10.0,10.0,0.0,0.0,0.0)
end

println("Testing slicing...")
slice1 = MeshSlice(binarymesh, [0.0,5.0,10.0])
Test.with_handler(test_handler) do
    for slice in slice1
        println("hi")
        for poly in slice.polygons
            println("hi again")
            n = length(poly.segments)
            for i = 1:n-1
                @test poly.segments[i].finish == poly.segments[i+1].start
            end
        end
    end
end

# test mesh rotation
println("Testing Mesh Rotation...")
r1 = rotate(1, 1, 0, 0, 0, 0, 1, 0, 0, pi/2)
r2 = rotate(1, 1, 0, 0, 0, 0, 0, 1, 0, pi/2)
r3 = rotate(1, 1, 0, 0, 0, 0, 0, 0, 1, pi/2)
rr1 = [1.0, 0, 1.0]
rr2 = [0.0, 1.0, -1.0]
rr3 = [-1.0, 1.0, 0.0]
Test.with_handler(test_handler) do
    for i = 1:3
        @test isapprox(r1[i] ,rr1[i])
        @test isapprox(r2[i] ,rr2[i])
        @test isapprox(r3[i] ,rr3[i])
    end
end

rotate!(binarymesh, 45.0, [1.0,0.0,0.0], [0.0,0.0,0.0])
Test.with_handler(test_handler) do
    @test binarymesh.bounds == Bounds(10.0,10.0,13.762255133518483,0.0,-8.509035245341185,0.0)
end
rotate!(binarymesh, 45.0, [0.0,1.0,0.0], [0.0,0.0,0.0])
Test.with_handler(test_handler) do
    @test binarymesh.bounds == Bounds(16.96357128682594,10.0,13.762255133518483,0.0,-8.509035245341185,-8.509035245341185)
end


text:
function jacvec(f::Function,x::Array{T,1},v::Array{T,1}) where {T <: AbstractFloat}

    ep = 1e-10
    jvec = imag(f(x+ep*im*v))/ep

    return jvec

end

function jacvec!(jvec::Array{T,1},f::Function,x::Array{T,1},v::AbstractArray{T,1}) where {T <: AbstractFloat}

    ep = 1e-10
    jvec .= imag(f(x+ep*im*v))/ep

end

function jacvec_inplace!(jvec::Array{T,1},f::Function,x::Array{T,1},v::AbstractArray{T,1}) where {T <: AbstractFloat}

    n = length(x)
    ffk = Array{Complex{T},1}(undef,n)

    ep = 1e-10
    f(ffk,x+ep*im*v)

    jvec .= imag(ffk)/ep

end


text:
import VimBindings: TextObject, a_word
@testset "a word textobject" begin
end


text:
using MakieRecipes
using Literate, MakieCore, CairoMakie
using Test

cd(@__DIR__)

@testset "Examples" begin
    literatedir = joinpath(@__DIR__, "..", "docs", "src", "literate")
    ispath("test_examples") || mkpath("test_examples")
    cd("test_examples") do
        @test try
            include(joinpath(literatedir, "examples.jl")) # execute the source file directly
            true
        catch
            false
        end
    end
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Beam(ndim, nst, nxe, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Beam(2, 1, 100, 1, :x, Line(2, 1), false),
  :properties => [1.0;],
  :x_coords => 0.0:0.01:1.0,
  :support => [
    (1, [0 1]),
    (101, [0 0])
    ],
  :limit => 100,
  :tol => 0.00001
)

data |> display
println()

@time m = p46(data)
println()

println("\nThe buckling load = $(m[1])\n")
buckling_df = DataFrame(
  translation = m[4][m[5][1,:]+1],
  rotation = m[4][m[5][2,:]+1]
)
display(buckling_df)
println()

if VERSION.minor < 6
  using Plots
  gr(size=(400,600))

  p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
  titles = ["p46.1a translation", "p46.1a rotation"]
  p[1] = plot(
    convert(Array, buckling_df[:translation]),
    ylim=(-0.1, 0.3), xlabel="node",
    ylabel="y translation [m]", color=:blue,
    marker=(:circle,1,0.1,stroke(1,:black)),
    title=titles[1], leg=false)
  p[2] = plot(
    convert(Array, buckling_df[:rotation]),
    ylim=(-1.0, 1.0), xlabel="node",
    ylabel="rotation [radians]", color=:red,
    marker=(:circle,1,0.1,stroke(1,:black)),
    title=titles[2], leg=false)

  plot(p..., layout=(2, 1))
  savefig(ProjDir*"/Ex46.1a.png")

end


text:
function gausslobatto(n)
    # Gauss-Legendre-Lobatto Quadrature Nodes and Weights
    if n == 1
        error("Lobatto undefined for n = 1.")
    elseif n == 2
        [-1.0, 1.0], [1.0, 1.0]
    elseif n == 3
        [-1.0, 0.0, 1.0], [1.0 / 3, 4.0 / 3, 1.0 / 3]
    else
        # Compute via GaussJacobi:
        x, w = gaussjacobi(n - 2, 1.0, 1.0)
        @inbounds for i in 1:length(x)
            w[i] = w[i] / (1 - x[i]^2)
        end
        pushfirst!(x, -1.0)
        push!(x, 1.0)
        pushfirst!(w, 2 / (n * (n - 1)))
        push!(w, 2 / (n * (n - 1)))
        x, w
    end
end


text:
module test_infer_geometry_based_bonds!

using Test, Xtals, XtalsPyTools, Graphs

function visual_check(xtal::String)
    c = Crystal(xtal)
    c = replicate(c, (2, 2, 2))
    strip_numbers_from_atom_labels!(c)
    infer_geometry_based_bonds!(c, true) # must
    write_xyz(c, "temp/c.xyz")
    write_bond_information(c, "temp/$(c.name)")
    @info c.name * " see .vtk and .xyz to visually check bonds"
end

if !isdir("temp")
    mkdir("temp")
end

@testset "scipy dep" begin
    scipy = rc[:scipy]
    rc[:scipy] = nothing
    @test_throws ErrorException infer_geometry_based_bonds!(Crystal("IRMOF-1.cif"), true)
    rc[:scipy] = scipy
end

@testset "NiPyC2 Tests" begin
    # NiPyC2 bonding
    c = Crystal("NiPyC2_relax.cif")
    strip_numbers_from_atom_labels!(c)
    c = replicate(c, (4, 4, 4))
    bonding_rules = [BondingRule(:H, :*, 1.2),
                     BondingRule(:Ni, :O, 2.5),
                     BondingRule(:Ni, :N, 2.5),
                     BondingRule(:*, :*, 1.9)]
    infer_bonds!(c, true, bonding_rules=bonding_rules)
    c1 = deepcopy(c)
    conn_comps = connected_components(c.bonds)

    @test length(conn_comps) == 2 # interpenetrated

    c_red = getindex(c, conn_comps[1])
    c_blue = getindex(c, conn_comps[2])

    @test ne(c_red.bonds) + ne(c_blue.bonds) == ne(c.bonds)

    remove_bonds!(c)
    infer_geometry_based_bonds!(c, true)

    @test length(conn_comps) == 2 # interpenetrated

    # debugging outputs (delete me)
    write_bond_information.([c, c1], ["temp/c", "temp/c1"])
    write_xyz(c, "temp/c")

    @test c1.bonds == c.bonds # consistency between two different bonding schemes
end

@testset "FIQCEN Tests" begin
    # FIQCEN bonding
    c = Crystal("FIQCEN_clean.cif")
    strip_numbers_from_atom_labels!(c)
    infer_geometry_based_bonds!(c, true)

    @test length(connected_components(c.bonds)) == 1 # not interpenetrated

    @test c.atoms.species[neighbors(c.bonds, 1)] == [:Cu, :O, :O, :O, :O]

    visual_check("FIQCEN_clean.cif")

    # reduce covalant radius to see Cu-Cu bond disappear
    covalent_radii = rc[:covalent_radii]
    covalent_radii[:Cu] = 1.15
    remove_bonds!(c)

    @test ne(c.bonds) == 0

    infer_geometry_based_bonds!(c, true, covalent_radii=covalent_radii)

    @test c.atoms.species[neighbors(c.bonds, 1)] == [:O, :O, :O, :O]
end

end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
using RandomVariates
using Gadfly
import Cairo, Fontconfig

N = 1000000
A = [1 .8; .8 1]

struct Dist
    title
    x
    xmin
    xmax
end

function plot_dist(d::Dist)
    plot(x=d.x, Stat.density, Geom.polygon(fill=true, preserve_order=true), Guide.title(d.title), Coord.cartesian(xmin=d.xmin, xmax=d.xmax))
end

function plot_hist(d::Dist)
    plot(x=d.x, Geom.histogram(density=true), Guide.title(d.title), Coord.cartesian(xmin=d.xmin, xmax=d.xmax))
end

uniform = Dist("Standard Uniform", uniform_rng(0,1,N), 0, 1)
normal = Dist("Standard Normal", normal_rng(0,1,N), nothing, nothing)
bernoulli = Dist("Bernoulli(.34)", bernoulli_rng(.34, N), nothing, nothing)
beta = Dist("Beta(2,8)", beta_rng(2,8, N), 0, 1)
binomial = Dist("Binomial(.43, 10)", binomial_rng(.43, 10, N), nothing, nothing)
erlang = Dist("Erlang(7, 4)", erlang_rng(7, 4, N), 0, nothing)
expon = Dist("Exponential(4)", expon_rng(4, N), 0, nothing)
gamma = Dist("Gamma(6,4)", gamma_rng(6,4,N), 0, nothing)
geom = Dist("Geometric(.2)", geometric_rng(.2, N), 0, nothing)
neg_binom = Dist("Negative Binomial(.25, 20)", conv_neg_binomial_rng(.25, 20, N), 0, nothing)
poisson = Dist("Poission(10)", poisson_rng(10, N), 0, nothing)
triang = Dist("Triangular(1,5,2)", triag_rng(1,5,2,N), nothing, nothing)
weibull = Dist("Weibull(1,5)", weibull_rng(1,5,N), 0, nothing)
mv_norm = Dist("Multivariate Normal([0,2], Σ)", mv_normal_rng([0,2], A, Int(N/1000)), nothing, nothing)

for d in [uniform, normal, expon, erlang, gamma, triang, weibull, beta]
    out = plot_dist(d)
    out |> SVG("images/$(d.title).svg", 5inch, 5inch)
    out |> PNG("images/$(d.title).png", 5inch, 5inch)
end

for d in [bernoulli, binomial, geom, neg_binom, poisson]
    out = plot_hist(d)
    out |> SVG("images/$(d.title).svg", 5inch, 5inch)
    out |> PNG("images/$(d.title).png", 5inch, 5inch)
end

out = plot(x=mv_norm.x[:, 1], y=mv_norm.x[:, 2], Geom.density2d(levels=10), Geom.point, alpha=[0.4], Guide.title(mv_norm.title))
out |> SVG("images/$(mv_norm.title).svg", 5inch, 5inch)
out |> PNG("images/$(mv_norm.title).png", 5inch, 5inch)


text:
module InteractiveUtilsPlus

export @show_sexpr, @show_tree, print_tree, show_expr, @show_expr, show_Sexpr, @show_Sexpr

using InteractiveUtils: subtypes
using AbstractTrees

"""
    @show_sexpr(expr, linenums=false)

shows the lisp style S-expression of `expr` and prints the line number nodes if `linenums` is true.

## Example

```julia
julia> @show_sexpr 2x+1
(:call, :+, (:call, :*, 2, :x), 1)
```
"""
macro show_sexpr(expr, linenums=false)
    linenums || Base.remove_linenums!(expr)
    :(Meta.show_sexpr($(QuoteNode(expr))))
end

AbstractTrees.printnode(io::IO, expr::Expr) = show(io, expr.head)

"""
    @show_tree(expr, maxdepth=10, linenums=false)

shows the tree form of the expression `expr` with maxdepth and prints the line number nodes if `linenums` is true.

## Example
```julia
julia> @show_tree 2x+1
:call
├─ :+
├─ :call
│  ├─ :*
│  ├─ 2
│  └─ :x
└─ 1
```
"""
macro show_tree(expr, maxdepth=10, linenums=false)
    linenums || Base.remove_linenums!(expr)
    :(print_tree($(QuoteNode(expr)), $(esc(maxdepth))))
end

AbstractTrees.children(T::Type) = subtypes(T)

"""
    print_tree(T::Type, maxdepth=5; kwargs...)

prints the subtree of the type `T`.

## Example
```julia
julia> print_tree(AbstractRange)
AbstractRange
├─ LinRange
├─ OrdinalRange
│  ├─ AbstractUnitRange
│  │  ├─ IdentityUnitRange
│  │  ├─ OneTo
│  │  ├─ Slice
│  │  └─ UnitRange
│  └─ StepRange
└─ StepRangeLen
```
"""
print_tree

"""
    show_expr([io::IO,], ex)

shows expression `ex` as a Julia style expression.

# Examples
```julia
julia> show_expr(:(f(x, g(y, z))))
Expr(:call, :f, :x,
    Expr(:call, :g, :y, :z))
```
"""
show_expr(ex; indent=expr_indent, head="Expr") =
    show_expr(stdout, ex; indent, head)
show_expr(io::IO, ex; indent=expr_indent, head="Expr") =
    show_expr(io, ex, 0; indent, head)
show_expr(io::IO, ex, inner; indent=expr_indent, head="Expr") =
    show(io, ex)

const expr_indent = 4

function show_expr(io::IO, ex::QuoteNode, inner; indent=expr_indent, head="Expr")
    print(io, "(:quote, ")
    show_expr(io, ex.value, inner+indent; indent)
    print(io, ')')
end
function show_expr(io::IO, ex::Expr, inner; indent=expr_indent, head="Expr")
    print(io, (iszero(inner) ? "" : "\n"), " "^inner, "$head(")
    show_expr(io, ex.head, inner+indent; indent, head)
    for arg in ex.args
        print(io, ", ")
        show_expr(io, arg, inner+indent; indent, head)
    end
    if isempty(ex.args)
        print(io, ",)")
    else
        print(io, ')')
    end
end

"""
    @show_expr(expr, linenums=false)

shows the Juia style expression of `expr` and prints the line number nodes if `linenums` is true.

## Example

```julia
julia> @show_expr 2x+1
Expr(:call, :+,
    Expr(:call, :*, 2, :x), 1)
```
"""
macro show_expr(expr, linenums=false)
    linenums || Base.remove_linenums!(expr)
    :(show_expr($(QuoteNode(expr))))
end

"""
    show_Sexpr([io::IO,], ex)

Show expression `ex` as a lisp style expression.

Remark: The indentation is different from `@show_sexpr`.

# Examples
```julia
julia> show_Sexpr(:(f(x, g(y, z))))
Expr(:call, :f, :x,
    Expr(:call, :g, :y, :z))
```
"""
show_Sexpr(io::IO, ex; indent=expr_indent) = show_expr(io, ex; indent, head="")
show_Sexpr(ex; indent=expr_indent) = show_Sexpr(stdout, ex; indent)

"""
    @show_Sexpr(expr, linenums=false)

shows the lisp style S-expression of `expr` and prints the line number nodes if `linenums` is true.

Remark: The indentation is different from `@show_sexpr`.

## Example

```julia
julia> @show_Sexpr 2x+1
(:call, :+,
    (:call, :*, 2, :x), 1)
```
"""
macro show_Sexpr(expr, linenums=false)
    linenums || Base.remove_linenums!(expr)
    :(show_Sexpr($(QuoteNode(expr))))
end

end


text:
################################################################################
#
#   bgjt.jl : functions for the BGJT algorithm
#
################################################################################

# Flint functions

function rrefMod!(M::Nemo.fmpz_mat, n::Nemo.fmpz)
#    perm::Array{Int, 1} = collect(1:rows(M))
    perm = C_NULL
    rank = ccall((:fmpz_mat_rref_mod, :libflint), Cint, (Ptr{Int}, Ptr{fmpz_mat},
                                                  Ptr{fmpz}), perm, &M, &n)
    return M, rank, perm
end

function rrefMod!(M::Nemo.fmpz_mat, d::Integer)
    n = Nemo.fmpz(d)
#    perm::Array{Int, 1} = collect(1:rows(M))
    perm = C_NULL
    rank = ccall((:fmpz_mat_rref_mod, :libflint), Cint, (Ptr{Int}, Ptr{fmpz_mat},
                                                  Ptr{fmpz}), perm, &M, &n)
    return M, rank, perm
end

function rrefMod(M::Nemo.fmpz_mat, n::Nemo.fmpz)
    N = deepcopy(M)
    M, rank, perm = rrefMod!(N, n)
    return N, rank, perm
end

function rrefMod(M::Nemo.fmpz_mat, n::Integer)
    N = deepcopy(M)
    M, rank, perm = rrefMod!(N, n)
    return N, rank, perm
end

# Julia functions

export subMatrix
"""
    subMatrix(M::MatElem, nrow::Integer, ncol::Integer)

Return the submatrix consisting of the first `nrow` rows and `ncol` colums of
`M`.
"""
function subMatrix(M::MatElem, nrow::Integer, ncol::Integer)

    # We create the matrix subspace
    S = MatrixSpace(parent(M[1,1]), nrow, ncol)

    # We create a new matrix
    sub = S()

    # And we copy the entries of the given matrix
    for j in 1:ncol
        for i in 1:nrow
            sub[i, j] = M[i, j]
        end
    end
    return sub
end

export pivots
"""
    pivots(M::MatElem)

Find `r` pivots in the row reduced echelon form matrix `M`.
"""
function pivots(M::MatElem, r::Int = 0)

    # There is as many pivots as the rank of
    if r == 0
        r = rank(M)
    end

    # We set an array with the coordinates of the pivots
    piv = Array(Int, r)
    c = 1

    # And we find these coordinates
    for i in 1:r
        while M[i, c] == 0
            c += 1
        end
        piv[i] = c
    end
    return piv
end


text:
type SGNHTHMCState <: SamplerState
  x::Array{Float64}
  p::Array{Float64}
  zeta::Array{Float64,1}

  niters::Int64
  stepsize::Float64
  mass::Array{Float64,1}
  D::Array{Float64,1}
  Best::Array{Float64,1}
  function SGNHTHMCState(x::Array{Float64};stepsize = 0.001, p=:none, mass=[1.0],niters=10, D=[1.0], Best=[0.0], zeta=[0.0])
    if isa(mass,Number)
      mass = mass * ones(length(x))
    end
    if p == :none
      p = sqrt(mass).*randn(length(x))
    end
    new(x,p,zeta,niters,stepsize,mass,D,Best)
  end
end

function sample!(s::SGNHTHMCState, llik, sgrad)
  D = s.D
  niters = s.niters
  stepsize = s.stepsize
  Best = s.Best
  zeta = s.zeta
  m = s.mass

  for iter=1:s.niters

    p_grad = zeta.*s.p./ m
    n = sqrt(stepsize.*(2D.-stepsize.*Best)).*randn(length(s.x))
    s.p[:] += stepsize.*(sgrad(s.x)-p_grad) + n
    s.x[:] += stepsize.*s.p./m
    zeta[:] += stepsize.* (s.p'*(s.p./(length(s.x)*m.^2)) - mean(1./m))
  end
  s
end


text:
using SearchLight

try
  SearchLight.Configuration.load()

  if SearchLight.config.db_config_settings["adapter"] !== nothing
    eval(Meta.parse("using SearchLight$(SearchLight.config.db_config_settings["adapter"])"))
    SearchLight.connect()

    @eval begin
      using Genie.Renderer.Json

      """
      $TYPEDSIGNATURES
      """
      function Genie.Renderer.Json.JSON3.StructTypes.StructType(::Type{T}) where {T<:SearchLight.AbstractModel}
        Genie.Renderer.Json.JSON3.StructTypes.Struct()
      end

      """
      $TYPEDSIGNATURES
      """
      function Genie.Renderer.Json.JSON3.StructTypes.StructType(::Type{SearchLight.DbId})
        Genie.Renderer.Json.JSON3.StructTypes.Struct()
      end
    end
  end
catch ex
  @error ex
end


text:
using LinearAlgebra
using FastLapackInterface

struct LyapdWs
    AA::Matrix{Float64}
    AAtemp::Matrix{Float64}
    AA2::Matrix{Float64}
    BB::Matrix{Float64}
    temp1::Matrix{Float64}
    XX::Vector{Float64}
    nonstationary_variables::Vector{Bool}
    nonstationary_trends::Vector{Bool}
    dgees_ws::DgeesWs
    linsolve_ws1::LinSolveWs
    linsolve_ws2::LinSolveWs
    function LyapdWs(n::Int64)
        AA = Matrix{Float64}(undef, n, n)
        AAtemp = Matrix{Float64}(undef, n, n)
        AA2 = Matrix{Float64}(undef, 2*n, 2*n)
        BB = similar(AA)
        temp1 = Matrix{Float64}(undef, n, n)
        XX = Vector{Float64}(undef, 2*n)
        nonstationary_variables = Vector{Bool}(undef, n)
        nonstationary_trends = Vector{Bool}(undef, n)
        dgees_ws = DgeesWs(n)
        linsolve_ws1 = LinSolveWs(n)
        linsolve_ws2 = LinSolveWs(2*n)
        new(AA, AAtemp, AA2, BB, temp1, XX, nonstationary_variables,
            nonstationary_trends, dgees_ws, linsolve_ws1,
            linsolve_ws2)
    end
end

function solve_one_row!(X::AbstractMatrix{Float64},
                        A::AbstractMatrix{Float64},
                        B::AbstractMatrix{Float64},
                        n::Int64, row::Int64, ws::LyapdWs)
    α = A[row, row]
    vA = view(A, 1:row, 1:row)
    vAA = view(ws.AAtemp, 1:row, 1:row)
    copy!(vAA, vA)
    lmul!(α, vAA)
    vAA .-= I(row)
    #exploiting symmetry to have row vectors
    vX = view(X, 1:row, row)
    vB = view(B, 1:row, row)
    vX .= .-vB
    linsolve_core!(vAA, vX, ws.linsolve_ws1)
    i = row
    j = n*(row - 1) + 1
    while i < row*n
        X[i] = X[j]
        i += n
        j += 1
    end
end

function solve_two_rows!(X::AbstractMatrix{Float64},
                         A::AbstractMatrix{Float64},
                         B::AbstractMatrix{Float64},
                         n::Int64, row::Int64, ws::LyapdWs)
    α11, α21, α12, α22 = A[(row - 1):row, (row - 1):row]
    l1 = 1
    l2 = 1
    while l1 <= row
        k1 = 1
        k2 = 1
        while k1 <= row
            ws.AA2[k2, l2] = α11*A[k1, l1]
            ws.AA2[k2  + 1, l2] = α21*A[k1, l1]
            ws.AA2[k2, l2 + 1] = α12*A[k1, l1]
            ws.AA2[k2 + 1, l2 + 1] = α22*A[k1, l1]
            k1 += 1
            k2 += 2
        end
        ws.XX[l2] = -B[row - 1, l1]
        ws.XX[l2 + 1] = -B[row, l1]
        l1 += 1
        l2 += 2
    end
    vAA = view(ws.AA2, 1:2*row, 1:2*row)
    vAA .-= I(2*row)
    vXX = view(ws.XX, 1:2*row)
    linsolve_core!(vAA, vXX, ws.linsolve_ws2)
    i = row - 1
    j = n*(row - 2) + 1
    m = 1
    while i <= row*n
        X[i] = vXX[m]
        X[i + 1] = vXX[m+1]
        X[j] = X[i]
        X[j + n] = X[i+1]
        i += n
        j += 1
        m += 2
    end
end

"""
    function lyapd(Σ, A, B, ws) solves equation Σ - A*Σ*A' = B
"""
function extended_lyapd!(Σ::AbstractMatrix{Float64},
                         A::AbstractMatrix{Float64},
                         B::AbstractMatrix{Float64},
                         ws::LyapdWs)
    n = size(A, 1)
    copy!(ws.AA, A)
    dgees!(ws.dgees_ws, ws.AA, >, 1-1e-6)
    mul!(ws.temp1, transpose(ws.dgees_ws.vs), B)
    mul!(ws.BB, ws.temp1, ws.dgees_ws.vs)

    extended_lyapd_core!(Σ, ws.AA, ws.BB, ws)
    mul!(ws.temp1, ws.dgees_ws.vs, Σ)
    mul!(Σ, ws.temp1, ws.dgees_ws.vs')

    for i = 1:n
        if ws.nonstationary_trends[i]
            for j = 1:n
                if abs(ws.dgees_ws.vs[j, i]) > 1e-10
                    ws.nonstationary_variables[j] = true
                end
            end
        end
    end

    for i = 1:n
        if  ws.nonstationary_variables[i]
            for j = i:n
                Σ[j, i] = NaN
                Σ[i, j] = NaN
            end
        else
            for j = i:n
                if  ws.nonstationary_variables[j]
                    Σ[j, i] = NaN
                    Σ[i, j] = NaN
                end
            end
        end
    end
end

function extended_lyapd_core!(Σ::AbstractMatrix{Float64},
                              A::AbstractMatrix{Float64},
                              B::AbstractMatrix{Float64},
                              ws::LyapdWs)
    fill!(Σ, 0.0)
    fill!(ws.nonstationary_trends, false)
    fill!(ws.nonstationary_variables, false)
    n = size(A, 1)
    row = n
    while row >= 1
        if row == 1 || A[row, row - 1] == 0
            if A[row, row] > 1 - 1e-6
                ws.nonstationary_trends[row] = true
            else
                solve_one_row!(Σ, A, B, n, row, ws)
                vB = view(B, 1:row - 1, 1:row - 1)
                vtemp = view(ws.temp1, 1:row - 1, 1)
                vΣ = view(Σ, 1:row -1 , row)
                vA = view(A, 1:row - 1, 1:row - 1)
                mul!(vtemp, vA, vΣ)
                vA = view(A, 1:row - 1, row)
                mul!(vB, vtemp, vA', 1.0, 1.0)
                vΣ = view(Σ, 1:row, row)
                vA = view(A, 1:row - 1, 1:row)
                mul!(vtemp, vA, vΣ)
                vA = view(A, 1:row - 1, row)
                vB .= vB .+ vA.*vtemp'
            end
            row -= 1
        else
            a = A[row, row]
            if a*a + A[row, row - 1]*A[row - 1, row] > 1 - 2e-6
                ws.nonstationary_trends[row] = true
                ws.nonstationary_trends[row - 1] = true
            else
                solve_two_rows!(Σ, A, B, n, row, ws)
                vB = view(B, 1:row - 2, 1:row - 2)
                vtemp = view(ws.temp1, 1:row - 2, 1:2)
                vΣ = view(Σ, 1:row - 2, row - 1:row)
                vA = view(A, 1:row - 2, 1:row - 2)
                mul!(vtemp, vA, vΣ)
                vA = view(A, 1:row - 2, row - 1:row)
                mul!(vB, vtemp, vA', 1.0, 1.0)
                vΣ = view(Σ, 1:row, row - 1:row)
                vA = view(A, 1:row - 2, 1:row)
                mul!(vtemp, vA, vΣ)
                vA = view(A, 1:row - 2, row - 1:row)
                mul!(vB, vA, vtemp', 1.0, 1.0)
            end
            row -= 2
        end
    end
end

function is_stationary(ws::LyapdWs)
    return !any(ws.nonstationary_trends)
end