text:
using ClutteredEnvPathOpt
using LinearAlgebra
using Test
using Pipe
using Plots
using JuMP
using Gurobi

################################# solve_deits ########################################
# Create obstacles
num_obs = 3;
seed = 10;
obstacles = ClutteredEnvPathOpt.gen_field(num_obs, seed = seed)
plot();
ClutteredEnvPathOpt.plot_field(obstacles)
# points = ClutteredEnvPathOpt.find_points(obstacles)
# ClutteredEnvPathOpt.plot_points(points)
display(plot!())

# Set parameters
N = 30  # number of steps
f1 = [0.0, 0.1, 0.0]  # initial footstep pose 1
f2 = [0.0, 0.0, 0.0]  # initial footstep pose 2
goal = [1.0, 1.0, 0.0]  # goal pose
Q_g = 10 * Matrix{Float64}(I, 3, 3)  # weight between final footstep and goal pose
Q_r = [1 0 0; 0 1 0; 0 0 0] # Matrix{Float64}(I, 3, 3)  # weight between footsteps
q_t = -.05  # weight for trimming unused steps
# TODO: Assure second footstep is within reachability given first

# Optional named arguments
d1 = 0.2
d2 = 0.2
p1 = [0. 0.07]
p2 = [0, -0.27]
delta_x_y_max = 0.1  # max stride norm in space
delta_θ_max = pi/4  # max difference in θ
L = 5  # number of pieces in p.w.l approx. of sine/cosine

# Compute optimal path
x1, y1, θ1, t1, stats1 = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method="merged")
x = copy(x1); y = copy(y1); θ = copy(θ1); t = copy(t1);
term_status, r_solve_time, rel_gap, simplex_iters, r_node_count, num_vertices, merged_size, full_size, num_free_faces, num_free_face_ineq, method_used = stats1

x2, y2, θ2, t2, stats2 = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method="full")
term_status2, r_solve_time2, rel_gap2, simplex_iters2, r_node_count2, num_vertices2, merged_size2, full_size2, num_free_faces2, num_free_face_ineq2, method_used2 = stats2

x3, y3, θ3, t3, stats3 = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method="bigM")
term_status3, r_solve_time3, rel_gap3, simplex_iters3, r_node_count3, num_vertices3, merged_size3, full_size3, num_free_faces3, num_free_face_ineq3, method_used3 = stats3

seed_range = 1:3
num_obs_range = 1:4
time_diffs = []
for seed in seed_range
    for num_obs in num_obs_range
        println("\nOn test Seed = $seed, Num_Obs = $num_obs")
        obstacles = ClutteredEnvPathOpt.gen_field(num_obs, seed = seed)
        _,_,_,_,stats = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method = "merged")
        _,_,_,_,stats2 = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method="full")
        _,_,_,_,stats3 = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method="bigM")
        term_status, solve_time, rel_gap, simplex_iters, node_count, num_vertices, merged_size, full_size, num_free_faces, num_free_face_ineq, method_used = stats
        term_status2, solve_time2, rel_gap2, simplex_iters2, node_count2, num_vertices2, merged_size2, full_size2, num_free_faces2, num_free_face_ineq2, method_used2 = stats2
        term_status3, solve_time3, rel_gap3, simplex_iters3, node_count3, num_vertices3, merged_size3, full_size3, num_free_faces3, num_free_face_ineq3, method_used3 = stats3
        # Merged biclique times vs bigM
        t_diff = solve_time3 - solve_time
        push!(time_diffs, t_diff)
    end
end

# Trim excess steps
num_to_trim = length(filter(tj -> tj > 0.5, t[3:end]))
if num_to_trim % 2 == 0
    x = vcat(x[1:2], x[num_to_trim + 3 : end]);
    y = vcat(y[1:2], y[num_to_trim + 3 : end]);
    θ = vcat(θ[1:2], θ[num_to_trim + 3 : end]);
else
    x = vcat(x[1], x[num_to_trim + 3 : end]);
    y = vcat(y[1], y[num_to_trim + 3 : end]);
    θ = vcat(θ[1], θ[num_to_trim + 3 : end]);
end

num_to_trim2 = length(filter(tj -> tj > 0.5, t2[3:end]))
if num_to_trim2 % 2 == 0
    x2 = vcat(x2[1:2], x2[num_to_trim2 + 3 : end]);
    y2 = vcat(y2[1:2], y2[num_to_trim2 + 3 : end]);
    θ2 = vcat(θ2[1:2], θ2[num_to_trim2 + 3 : end]);
else
    x2 = vcat(x2[1], x2[num_to_trim2 + 3 : end]);
    y2 = vcat(y2[1], y2[num_to_trim2 + 3 : end]);
    θ2 = vcat(θ2[1], θ2[num_to_trim2 + 3 : end]);
end

num_to_trim3 = length(filter(tj -> tj > 0.5, t3[3:end]))
if num_to_trim3 % 2 == 0
    x3 = vcat(x3[1:2], x3[num_to_trim3 + 3 : end]);
    y3 = vcat(y3[1:2], y3[num_to_trim3 + 3 : end]);
    θ3 = vcat(θ3[1:2], θ3[num_to_trim3 + 3 : end]);
else
    x3 = vcat(x3[1], x3[num_to_trim3 + 3 : end]);
    y3 = vcat(y3[1], y3[num_to_trim3 + 3 : end]);
    θ3 = vcat(θ3[1], θ3[num_to_trim3 + 3 : end]);
end

# Plot footstep plan
plot_steps(obstacles, x, y, θ)
plot_steps(obstacles, x2, y2, θ2)
plot_steps(obstacles, x3, y3, θ3)
png("Path Seed $seed Num Obs $num_obs")

# Plot intersections of circles
plot_circles(x, y, θ)
plot_circles(x2, y2, θ2)
plot_circles(x3, y3, θ3)

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

num_obs = 4;
seed = 3;
obstacles, points, g, obstacle_faces, free_faces = ClutteredEnvPathOpt.plot_new(num_obs, "Obstacles Seed $seed Num Obs $num_obs", seed = seed)
skeleton = LabeledGraph(g)
all_faces = union(obstacle_faces, free_faces)

# Plot obstacles
plot(title="Obstacles");
ClutteredEnvPathOpt.plot_field(obstacles)
# points = ClutteredEnvPathOpt.find_points(obstacles)
ClutteredEnvPathOpt.plot_points(points, vertices=skeleton.labels)
# ClutteredEnvPathOpt.plot_lines(obstacles)
# ClutteredEnvPathOpt.plot_borders()
# ClutteredEnvPathOpt.plot_intersections(obstacles);
png("Obstacles Seed $seed Num Obs $num_obs")
display(plot!())

# Plot faces
ClutteredEnvPathOpt.plot_faces(obstacle_faces, points, plot_name = "Obstacle Faces", col = "dodgerblue")
png("Obstacle Faces Seed $seed Num Obs $num_obs")

ClutteredEnvPathOpt.plot_faces(free_faces, points, plot_name = "Free Faces")
png("Free Faces Seed $seed Num Obs $num_obs")

ClutteredEnvPathOpt.plot_faces(obstacle_faces, points, plot_name = "All Faces", col = "red")
ClutteredEnvPathOpt.plot_faces(free_faces, points, plot_name = "All Faces", col = "green", new_plot = false)
png("All Faces Seed $seed Num Obs $num_obs")

# Plot faces individually
ClutteredEnvPathOpt.plot_faces(obstacle_faces, points, col = "dodgerblue", individually = true)
ClutteredEnvPathOpt.plot_faces(free_faces, points, individually = true)

# Plot edges
ClutteredEnvPathOpt.plot_edges(skeleton, points, plot_name = "Skeleton", vertices=skeleton.labels)
png("Skeleton Seed $seed Num Obs $num_obs")

feg_S = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces))
ClutteredEnvPathOpt.plot_edges(feg_S, points, plot_name = "Finite Element Graph of S",col="black", vertices=feg_S.labels)
png("FEG of S Seed $seed Num Obs $num_obs")

# feg_obs = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(obstacle_faces))
# ClutteredEnvPathOpt.plot_edges(feg_obs, points, plot_name = "Finite Element Graph Obstacle Faces")
# png("FEG Obstacles Seed $seed Num Obs $num_obs")

# feg_all = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(all_faces))
# ClutteredEnvPathOpt.plot_edges(feg_all, points, plot_name = "Finite Element Graph")
# png("FEG All Seed $seed Num Obs $num_obs")

# # Duplicate tests
# obstacles = ClutteredEnvPathOpt.gen_field(num_obs, seed = seed)
# points, mapped, inside_quant = ClutteredEnvPathOpt.find_intersections(obstacles)
# neighbors = ClutteredEnvPathOpt._find_neighbors(points, mapped)
# points = ClutteredEnvPathOpt.find_points(obstacles)

# dup, dup_ind = face_duplicates(obstacle_faces)
# dup, dup_ind = face_duplicates(free_faces)
# dup, dup_ind = face_duplicates(all_faces)
# plot_intersections_indiv(points)
# dup, dup_ind = point_duplicates(points)
# dup, dup_ind = points_in_hs_duplicates(mapped)
# dup, dup_ind = hs_duplicates(mapped)
# points_in_mapped_ordered(obstacles, mapped)
# plot_neighbors(obstacles, neighbors, points)
# list_neighbors(neighbors, points)
# suspect = suspect_neighbors(neighbors)

# cover3 = ClutteredEnvPathOpt._find_biclique_cover_visual(skeleton, free_faces, points)
# tree = ClutteredEnvPathOpt.find_biclique_cover_as_tree(skeleton, free_faces)

# Free Faces Cover
cover = ClutteredEnvPathOpt.find_biclique_cover(skeleton, free_faces)
(valid_cover, size_diff, missing_edges, missing_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces)), cover)

feg_free = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces))
merged_cover = ClutteredEnvPathOpt.biclique_merger(cover, feg_free)
(valid_cover, size_diff, missing_edges, missing_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces)), merged_cover)

ClutteredEnvPathOpt.plot_edges(missing_edges_lg, points, plot_name = "Missing Edges",vertices=missing_edges_lg.labels)
png("Missing Edges Seed $seed Num Obs $num_obs Free")

# file_name = "Biclique Cover Debug Output Seed $seed Num Obs $num_obs Free Unedited.txt"
# f = open(file_name, "w")
# write(f, "Plot 1\n")
# flush(f)
# close(f)
# cover2 = ClutteredEnvPathOpt.find_biclique_cover_debug(skeleton, free_faces, points, obstacles, file_name)

# To produce individual FEGs/Skeletons in the recursion process
(C,A,B), skeleton_ac, faces_ac, skeleton_bc, faces_bc = ClutteredEnvPathOpt.find_biclique_cover_one_iter(skeleton, free_faces)
feg_S_ac = ClutteredEnvPathOpt._find_finite_element_graph(skeleton_ac, ClutteredEnvPathOpt._find_face_pairs(faces_ac))
ClutteredEnvPathOpt.plot_edges(feg_S_ac, points, plot_name = "Finite Element Graph of S_ac",vertices=feg_S_ac.labels,col="black")
# plot!(title="",axis=([], false))
png("FEG of S_ac Seed $seed Num Obs $num_obs")
feg_S_bc = ClutteredEnvPathOpt._find_finite_element_graph(skeleton_bc, ClutteredEnvPathOpt._find_face_pairs(faces_bc))
ClutteredEnvPathOpt.plot_edges(feg_S_bc, points, plot_name = "Finite Element Graph of S_bc", vertices=feg_S_bc.labels, col="black")
# plot!(title="",axis=([], false))
png("FEG of S_bc Seed $seed Num Obs $num_obs")

# Biclique Cover Validity Tests-----------------------------------------------------------------------
seed_range = 1:100
num_obs_range = 1:4
# fails, tuples = biclique_cover_validity_tests(seed_range, num_obs_range, faces = "all")
fails, tuples = biclique_cover_validity_tests(seed_range, num_obs_range, faces = "free", with_plots=true)

percent_decreases, tuples, c_sizes, mc_sizes = biclique_cover_merger_stats(seed_range, num_obs_range)
percent_decreases_hp, tuples_hp, c_sizes_hp, mc_sizes_hp = biclique_cover_merger_stats(seed_range, num_obs_range, file_name="Biclique Cover Merging Stats HP Partition.txt", partition="hp")
size_diff_cover = c_sizes_hp - c_sizes
size_diff_merged_cover = mc_sizes_hp - mc_sizes

# TODO: WRITE STATS TO FILE

# Biclique validity tester function
function biclique_cover_validity_tests(seed_range, num_obs_range; faces::String="free", with_plots::Bool=false)
    fails = 0
    tuples = []
    for seed = seed_range
        for num_obs = num_obs_range
            println("On test Seed = $seed, Num_Obs = $num_obs")
            obstacles, points, g, obstacle_faces, free_faces = ClutteredEnvPathOpt.plot_new(num_obs,"Obstacles Seed $seed Num Obs $num_obs",seed=seed); # optional seed argument, default is 11
            skeleton = LabeledGraph(g)
            all_faces = union(obstacle_faces, free_faces)
            if faces == "all"
                cover = ClutteredEnvPathOpt.find_biclique_cover(skeleton, all_faces);
                feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(all_faces));
                merged_cover = ClutteredEnvPathOpt.biclique_merger(cover, feg)
                # (valid_cov, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, cover);
                (valid_cov, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, merged_cover);
                println("Difference: $(length(cover)-length(merged_cover)) | Merged: $(length(merged_cover)) | Original: $(length(cover))")
            else
                cover = ClutteredEnvPathOpt.find_biclique_cover(skeleton, free_faces);
                feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces));
                merged_cover = ClutteredEnvPathOpt.biclique_merger(cover, feg)
                # (valid_cov, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, cover);
                (valid_cov, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, merged_cover);
                println("Difference: $(length(cover)-length(merged_cover)) | Merged: $(length(merged_cover)) | Original: $(length(cover))")
            end
            if size_diff != 0
                fails += 1
                push!(tuples, (seed, num_obs))
                println("Test Failed")

                if with_plots
                    plot(title="Obstacles");
                    ClutteredEnvPathOpt.plot_field(obstacles)
                    # points = ClutteredEnvPathOpt.find_points(obstacles)
                    ClutteredEnvPathOpt.plot_points(points)
                    # ClutteredEnvPathOpt.plot_lines(obstacles)
                    # ClutteredEnvPathOpt.plot_borders()
                    # ClutteredEnvPathOpt.plot_intersections(obstacles);
                    display(plot!())
                    png("Obstacles Seed $seed Num Obs $num_obs")

                    ClutteredEnvPathOpt.plot_faces(free_faces, points, plot_name = "Free Faces")
                    png("Free Faces Seed $seed Num Obs $num_obs")

                    ClutteredEnvPathOpt.plot_edges(skeleton, points, plot_name = "Skeleton")
                    png("Skeleton Seed $seed Num Obs $num_obs")

                    if size_diff > 0
                        ClutteredEnvPathOpt.plot_edges(miss_edges_lg, points, plot_name = "Missing Edges")
                        if faces == "all"
                            png("Missing Edges All Faces Used Seed $seed Num Obs $num_obs")
                        else
                            png("Missing Edges Free Faces Used Seed $seed Num Obs $num_obs")
                        end
                    end
                end
            end
        end
    end
    return (fails, tuples)
end

# Biclique cover merger function for stats
function biclique_cover_merger_stats(seed_range, num_obs_range; file_name="Biclique Cover Merging Stats.txt", partition="CDT")
    # file_name = "Biclique Cover Merging Stats.txt" #Seed Range = $seed_range, Num Obs Range = $num_obs_range.txt"
    f = open(file_name, "w")
    write(f, "Seed Range = $seed_range, Num Obs Range = $num_obs_range\n")
    flush(f)
    close(f)

    f = open(file_name, "a")
    write(f, "Seed\tNum_obs\tFull_cover\tMerged_cover\tDecrease\tPercent_decrease\n")
    fails = 0
    tuples = []
    percent_vec = []
    cover_vec = []
    merged_cover_vec = []
    for seed = seed_range
        for num_obs = num_obs_range
            # f = open(file_name, "a")
            # write(f, "\nTest Seed = $seed, Num_Obs = $num_obs\n")
            write(f, "$seed\t$num_obs")
            # flush(f)
            # close(f)

            # println("On test Seed = $seed, Num_Obs = $num_obs")
            obstacles, points, g, obstacle_faces, free_faces = ClutteredEnvPathOpt.plot_new(num_obs,"Obstacles Seed $seed Num Obs $num_obs",seed=seed;partition=partition); # optional seed argument, default is 11
            skeleton = LabeledGraph(g)
            all_faces = union(obstacle_faces, free_faces)

            cover = ClutteredEnvPathOpt.find_biclique_cover(skeleton, free_faces);
            feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, ClutteredEnvPathOpt._find_face_pairs(free_faces));
            merged_cover = ClutteredEnvPathOpt.biclique_merger(cover, feg)
            (valid_cover, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, cover);
            (valid_cover_merged, size_diff, miss_edges, miss_edges_lg) = ClutteredEnvPathOpt._is_valid_biclique_cover_diff(feg, merged_cover);
            # println("Difference: $(length(cover)-length(merged_cover)) | Merged: $(length(merged_cover)) | Original: $(length(cover))")

            # if valid_cover && valid_cover_merged
            #     write(f,"Difference: $(length(cover)-length(merged_cover)) | Merged: $(length(merged_cover)) | Original: $(length(cover))")
            #     percent = round( ((length(cover)-length(merged_cover)) / length(cover)) * 100, digits = 2)
            #     write(f," | Percent Decrease: $percent\n")
            #     push!(percent_vec, percent)
            # else
            #     fails += 1
            #     push!(tuples, (seed, num_obs))
            #     if !valid_cover_merged
            #         write(f,"Merged cover invalid.\n")
            #     else
            #         write(f,"Original cover invalid.\n")
            #     end
            # end
            if valid_cover && valid_cover_merged
                percent = round( ((length(cover)-length(merged_cover)) / length(cover)) * 100, digits = 2)
                write(f,"\t$(length(cover))\t$(length(merged_cover))\t$(length(cover)-length(merged_cover))\t$percent\n")
                push!(percent_vec, percent)
                push!(cover_vec, length(cover))
                push!(merged_cover_vec, length(merged_cover))
            elseif valid_cover
                percent = round( ((length(cover)-0) / length(cover)) * 100, digits = 2)
                write(f,"\t$(length(cover))\t0\t$(length(cover)-0)\t$percent\n")
                push!(percent_vec, percent)
                push!(cover_vec, -1)
                push!(merged_cover_vec, -1)

                fails += 1
                push!(tuples, (seed, num_obs))
                # if !valid_cover_merged
                #     write(f,"Merged cover invalid.\n")
                # else
                #     write(f,"Original cover invalid.\n")
                # end
            end
        end
    end

    # write(f, "\nTotal Fails: $fails.\nTuples: $tuples")
    flush(f)
    close(f)

    return percent_vec, tuples, cover_vec, merged_cover_vec
end

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

# method <- "merged" for the compact biclique cover, "full" for the original biclique cover, "bigM" for big-M constraints
function solve_time_stats(seed_range, num_obs_range; method="merged", file_name="Solve Time Stats.txt")
    # file_name = "Solve Time Stats.txt" #Seed Range = $seed_range, Num Obs Range = $num_obs_range.txt"
    # f = open(file_name, "a")
    # write(f, "Method = $method, Seed Range = $seed_range, Num Obs Range = $num_obs_range\n")
    # flush(f)
    # close(f)

    # Set parameters
    N = 30  # number of steps
    f1 = [0.0, 0.1, 0.0]  # initial footstep pose 1
    f2 = [0.0, 0.0, 0.0]  # initial footstep pose 2
    goal = [1, 1, 0]  # goal pose
    Q_g = 10*Matrix{Float64}(I, 3, 3)  # weight between final footstep and goal pose
    Q_r = [1 0 0; 0 1 0; 0 0 0] # Matrix{Float64}(I, 3, 3)  # weight between footsteps
    q_t = -.05  # weight for trimming unused steps
    # Optional named arguments
    # d1 = 0.2 <- radius of reference foot circle
    # d2 = 0.2 <- radius of moving foot circle
    # p1 = [0, 0.07] <- center of reference foot circle
    # p2 = [0, -0.27] <- center of moving foot circle
    # delta_x_y_max = 0.10  # max stride norm in space
    # delta_θ_max = pi/4  # max difference in θ
    # L = 5  # number of pieces in p.w.l approx. of sine/cosine

    f = open(file_name, "a")
    # write(f, "Seed\tNum_obs\tTime\n")
    write(f, "Seed\tNum_obs\tSolve_time\tRel_gap\tSimplex_iterations\tNodes_explored\tNum_vertices\tMerged_size\tFull_size\tNum_free_faces\tNum_free_face_inequalities\n")
    times = zeros(length(seed_range) * length(num_obs_range))
    i = 1
    for seed = seed_range
        for num_obs = num_obs_range
            # f = open(file_name, "a")
            # write(f, "Seed = $seed, Num_Obs = $num_obs\n")
            # flush(f)
            # close(f)
            # println("On test Seed = $seed, Num_Obs = $num_obs")

            # Create obstacles
            obstacles = ClutteredEnvPathOpt.gen_field(num_obs, seed = seed)

            x, y, θ, t, stats = solve_deits(obstacles, N, f1, f2, goal, Q_g, Q_r, q_t, method=method)
            term_status, r_solve_time, rel_gap, simplex_iters, r_node_count, num_vertices, merged_size, full_size, num_free_faces, num_free_face_ineq, method_used = stats

            if method_used == method
                times[i] = r_solve_time # time
            else
                times[i] = -1
            end
            write(f, "$seed\t$num_obs\t$(times[i])\t$rel_gap\t$simplex_iters\t$r_node_count\t$num_vertices\t$merged_size\t$full_size\t$num_free_faces\t$num_free_face_ineq\n")
            flush(f)

            i += 1

            # # Trim excess steps
            # num_to_trim = length(filter(tj -> tj > 0.5, t[3:end]))
            # if num_to_trim % 2 == 0
            #     x = vcat(x[1:2], x[num_to_trim + 3 : end]);
            #     y = vcat(y[1:2], y[num_to_trim + 3 : end]);
            #     θ = vcat(θ[1:2], θ[num_to_trim + 3 : end]);
            # else
            #     x = vcat(x[1], x[num_to_trim + 3 : end]);
            #     y = vcat(y[1], y[num_to_trim + 3 : end]);
            #     θ = vcat(θ[1], θ[num_to_trim + 3 : end]);
            # end
            # plot_steps(obstacles, x, y, θ)
            # plot_circles(x, y, θ, p1=p1, p2=p2)

        end
    end

    # write(f, "\nEND")
    # flush(f)
    close(f)

    return times

end

seed_range = 21:30 # done 1:20
num_obs_range = 1:4
methods = ("merged", "full", "bigM")
times = Dict()
# method = "merged"
for method in methods
    file_name = "Solve Time Stats Method $method.txt"
    f = open(file_name, "a")   # write or append appropriately
    write(f, "Seed Range = $seed_range, Num Obs Range = $num_obs_range\nMethod = $method\n")
    flush(f)
    close(f)
    time = solve_time_stats(seed_range, num_obs_range, method=method, file_name=file_name)
    times[method] = time
end

file_name = "Solve Time Stats Comparison.txt" #Seed Range = $seed_range, Num Obs Range = $num_obs_range.txt"
f = open(file_name, "a")    # write or append appropriately
times_m_b = times["bigM"] .- times["merged"]
times_m_f = times["full"] .- times["merged"]
times_f_b = times["bigM"] .- times["full"]
write(f, "Seed Range = $seed_range, Num Obs Range = $num_obs_range\n")
write(f, "Merged_vs_BigM\tMerged_vs_Full\tFull_vs_BigM\n")
for i = 1:length(times_m_b)
    write(f, "$(times_m_b[i])\t$(times_m_f[i])\t$(times_f_b[i])\n")
end
flush(f)
close(f)

import DelimitedFiles
data_merged = DelimitedFiles.readdlm("Solve Time Stats Method merged.txt", '\t')
matrix_merged = Matrix{Float64}(data_merged[4:end,:])
data_bigM = DelimitedFiles.readdlm("Solve Time Stats Method bigM.txt", '\t')
matrix_bigM = Matrix{Float64}(data_bigM[4:end,:])
header = Vector{String}(data_merged[3,:])

file_name = "Solve Stats Comparison.txt" #Seed Range = $seed_range, Num Obs Range = $num_obs_range.txt"
f = open(file_name, "w")    # write or append appropriately
write(f, "S_T_m\tS_T_b\tS_T_diff\tS_I_m\tS_I_b\tS_I_diff\tN_C_m\tN_C_b\tN_C_diff\n")
S_T_m = matrix_merged[:,3]
S_T_b = matrix_bigM[:,3]
S_T_diff = S_T_b .- S_T_m
S_I_m = matrix_merged[:,5]
S_I_b = matrix_bigM[:,5]
S_I_diff = S_I_b .- S_I_m
N_C_m = matrix_merged[:,6]
N_C_b = matrix_bigM[:,6]
N_C_diff = N_C_b .- N_C_m
for i = 1:length(S_T_m)
    write(f, "$(S_T_m[i])\t$(S_T_b[i])\t$(S_T_diff[i])\t")
    write(f, "$(S_I_m[i])\t$(S_I_b[i])\t$(S_I_diff[i])\t")
    write(f, "$(N_C_m[i])\t$(N_C_b[i])\t$(N_C_diff[i])\n")
end
flush(f)
close(f)

data_summary = DelimitedFiles.readdlm("Solve Stats Comparison.txt", '\t')
matrix_summary = Matrix{Float64}(data_summary[2:end,:])
header_summary = Vector{String}(data_summary[1,:])

# Debugging Functions ---------------------------------------------------------------

# Check if free_faces has any duplicates---------------------------------------------
function face_duplicates(faces)
    dup = 0
    dup_ind = []
    col_faces = collect(faces)
    for i in 1:length(col_faces)
        for j = (i+1):length(col_faces)
            sort_i = sort(col_faces[i])
            sort_j = sort(col_faces[j])
            if sort_i == sort_j
                dup += 1
                push!(dup_ind, (i,j))
            end
        end
    end
    return dup, dup_ind
end

# Check if points has any duplicates---------------------------------------------
function point_duplicates(points)
    dup = 0
    dup_ind = []
    tol = Rational(1.e-13)
    for i in 1:(length(points))
        for j in (i+1):(length(points))
            #if (points[i].first == points[j].first) && (points[i].second == points[j].second)
            if abs(points[i].first - points[j].first) < tol && abs(points[i].second - points[j].second) < tol
                dup += 1
                push!(dup_ind, (i,j))
            end
        end
    end
    # for (i,j) in dup_ind
    #     println("Points $i and $j: \n ($(points[i])) \n ($(points[j]))")
    # end
    return dup, dup_ind
end

# Check if halfspaces in mapped have any duplicated points within them ---------
function points_in_hs_duplicates(mapped)
    dup = 0
    dup_ind = []
    for (l,hs) in enumerate(mapped)
        for i in 1:(length(hs))
            for j in (i+1):(length(hs))
                if (hs[i].first == hs[j].first) && (hs[i].second == hs[j].second)
                    dup += 1
                    push!(dup_ind, (l,i,j))
                end
            end
        end
    end
    return dup, dup_ind
end

# Check if any halfspaces in mapped are duplicates-----------------------------
function hs_duplicates(mapped)
    dup = 0
    dup_ind = []
    for i = 1:length(mapped)
        for j = (i+1):length(mapped)
            if mapped[i] == mapped[j]
                dup += 1
                push!(dup_ind, (i,j))
            end
        end
    end
    return dup, dup_ind
end

# Check to see if points are in order in mapped---------------------------------
function points_in_mapped_ordered(obstacles, mapped)
    for hs in mapped
        plot(legend=false,xlims=(-0.05,1.05),ylims=(-0.05,1.05));
        ClutteredEnvPathOpt.plot_field(obstacles)
        ClutteredEnvPathOpt.plot_lines(obstacles)
        for point in hs
            display(scatter!([point.first], [point.second]))
        end
    end
end

# Plot points individually----------------------------------------------------------
function plot_intersections_indiv(intersections)
    x = map(point -> point[1], intersections)
    y = map(point -> point[2], intersections)
    Plots.plot(legend=false,xlims=(-0.05,1.05),ylims=(-0.05,1.05))
    for i = 1:length(intersections)
        # Plots.plot(legend=false,xlims=(-0.05,1.05),ylims=(-0.05,1.05))
        display(Plots.scatter!([x[i]],[y[i]],title="Point $i:\n$(x[i])\n$(y[i])"))
    end
end

# Check if neighbors are correct----------------------------------------------------
function plot_neighbors(obstacles, neighbors, points)
    for (i,point) in enumerate(points)
        plot();
        ClutteredEnvPathOpt.plot_field(obstacles)
        ClutteredEnvPathOpt.plot_lines(obstacles)
        display(scatter!([point.first],[point.second],title="Point $i"))
        x_i_neighbors = map(node -> points[node].first, neighbors[i])
        y_i_neighbors = map(node -> points[node].second, neighbors[i])
        display(scatter!(x_i_neighbors,y_i_neighbors, title="Point $i. Neighbors $(neighbors[i])"))
        #println("Point $i has neighbors $(neighbors[i]).")
    end
end

function list_neighbors(neighbors, points)
    for i in 1:length(points)
        println("Point $i has neighbors $(neighbors[i]).")
    end
end

function duplicates(n)
    dup = 0
    for i in 1:length(n)
        for j in (i+1):length(n)
            if n[i] == n[j]
                dup += 1
            end
        end
    end
    return dup != 0
end

function suspect_neighbors(neighbors)
    suspect = []
    for n in 1:length(neighbors)
        if length(neighbors[n]) > 4 || duplicates(neighbors[n])
        #if duplicates(neighbors[n])
            push!(suspect, n)
        end
    end
    for point in suspect
        println("Point $point has neighbors $(neighbors[point])")
    end
    return suspect
end

# # Face finder ----------------------------------------------------------------------------
# t_seed = 7
# t_num_obs = 4
# obstacles = ClutteredEnvPathOpt.gen_field(t_num_obs, seed = t_seed) # optional seed
# construct_graph_debug(obstacles)


text:
function gxbprint(io::IO, x, name = "", level::libgb.GxB_Print_Level = libgb.GxB_SUMMARY)
    #if x isa AbstractGBType
    #    libgb.GxB_Type_fprint(x, name, level, C_NULL)
    #elseif x isa libgb.GrB_UnaryOp
    #    libgb.GxB_UnaryOp_fprint(x, name, level, C_NULL)
    #elseif x isa libgb.GrB_BinaryOp
    #    libgb.GxB_BinaryOp_fprint(x, name, level, C_NULL)
    #elseif x isa libgb.GrB_Monoid
    #    libgb.GxB_Monoid_fprint(x, name, level, C_NULL)
    #elseif x isa SelectUnion
    #    libgb.GxB_SelectOp_fprint(x, name, level, C_NULL)
    #elseif x isa libgb.GrB_Semiring
    #    libgb.GxB_Semiring_fprint(x, name, level, C_NULL)
    #elseif x isa AbstractDescriptor
    #    libgb.GxB_Descriptor_fprint(x, name, level, C_NULL)
    #elseif x isa GBVector
    #    if level == libgb.GxB_SUMMARY
    #        libgb.GxB_Matrix_fprint(x, name, libgb.GxB_SHORT, C_NULL)
    #    else
    #        libgb.GxB_Matrix_fprint(x, name, level, C_NULL)
    #    end
    #elseif x isa GBMatrix
    #    if level == libgb.GxB_SUMMARY
    #        libgb.GxB_Matrix_fprint(x, name, libgb.GxB_SHORT, C_NULL)
    #    else
    #        libgb.GxB_Vector_fprint(x, name, level, C_NULL)
    #    end
    #elseif x isa GBScalar
    #    libgb.GxB_Scalar_fprint(x, name, libgb.GxB_SHORT, C_NULL)
    #end
    print(io, gxbstring(x, name, level))
end

function gxbstring(x, name = "", level::libgb.GxB_Print_Level = libgb.GxB_SUMMARY)
    str = mktemp() do _, f
        cf = Libc.FILE(f)
        if x isa AbstractGBType
            libgb.GxB_Type_fprint(x, name, level, cf)
        elseif x isa libgb.GrB_UnaryOp
            libgb.GxB_UnaryOp_fprint(x, name, level, cf)
        elseif x isa libgb.GrB_BinaryOp
            libgb.GxB_BinaryOp_fprint(x, name, level, cf)
        elseif x isa libgb.GrB_Monoid
            libgb.GxB_Monoid_fprint(x, name, level, cf)
        elseif x isa SelectUnion
            libgb.GxB_SelectOp_fprint(x, name, level, cf)
        elseif x isa libgb.GrB_Semiring
            libgb.GxB_Semiring_fprint(x, name, level, cf)
        elseif x isa AbstractDescriptor
            libgb.GxB_Descriptor_fprint(x, name, level, cf)
        elseif x isa GBVector
            if level == libgb.GxB_SUMMARY
                libgb.GxB_Matrix_fprint(x, name, libgb.GxB_SHORT, cf)
            else
                libgb.GxB_Matrix_fprint(x, name, level, cf)
            end
        elseif x isa GBMatrix
            if level == libgb.GxB_SUMMARY
                libgb.GxB_Matrix_fprint(x, name, libgb.GxB_SHORT, cf)
            else
                libgb.GxB_Vector_fprint(x, name, level, cf)
            end
        elseif x isa GBScalar
            libgb.GxB_Scalar_fprint(x, name, libgb.GxB_SHORT, cf)
        end
        flush(f)
        seekstart(f)
        x = read(f, String)
        close(cf)
        x
    end
    return str
end

function getdocstring(x)
    str = gxbstring(x)
    return strip(split(str, '\n')[2])
end

ext: jl


text:
__precompile__()

module Tools

#################### utilities
function discount_sum{DT<:AbstractFloat}(x::Vector{DT}, gamma::DT, terminal::DT)
   y = similar(x)
   runsum = terminal
   for t=length(x):-1:1
      runsum = x[t] + gamma*runsum
      y[t] = runsum
   end
   return y
end

function discount_sum!{DT<:AbstractFloat}(y::Matrix{DT}, x::Matrix{DT},
                                          t::Int, gamma::DT, terminal::DT)
   runsum = terminal
   for i=size(x, 1):-1:1
      runsum = x[i,t] + gamma*runsum
      y[i,t] = runsum
   end
end

function discount_sum!{T<:AbstractFloat,
                       DT<:AbstractFloat}(y::Matrix{T}, x::Matrix{DT},
                                          t::Int, gamma::DT, terminal::DT)
   runsum = terminal
   for i=size(x, 1):-1:1
      runsum = x[i,t] + gamma*runsum
      y[i,t] = runsum
   end
end

# reward: T x numT
function compute_returns!{DT<:AbstractFloat}(returns::Vector{DT},
                                             rewards::Matrix{Float64},
                                             gamma::Float64)
   T, numT = size(rewards)
   mat_ret = reshape(returns, T, numT)
   for t=1:numT
      discount_sum!(mat_ret, rewards, t, gamma, 0.0)
   end
end

function compute_advantages(returns::Vector{Float64},
                            baseline::Vector{Float64},
                            obs::Matrix{Float64},
                            gamma::Float64)
   return returns - baseline
end

function computeGAEadvantages!{DT<:AbstractFloat}(adv::Vector{DT},
                                                  baseline::Vector,
                                                  rewards::Matrix{Float64},
                                                  obs::Matrix{Float64},
                                                  gae::Float64,
                                                  gamma::Float64)
   T, numT = size(rewards)
   tdsum = Array{DT}(size(rewards))
   mat_adv = reshape(adv, T, numT)

   nthread = min(numT, Threads.nthreads()) # cant have more threads than data
   Threads.@threads for tid=1:nthread # WHY IS MT SO MUCH SLOWER??
      thread_range = Distributed.splitrange(numT, nthread)[tid]

      for t=thread_range
         base = view(baseline, (t-1)*T+1:t*T)
         for k=1:T-1
            tdsum[k,t] = rewards[k,t] + gamma * base[k+1] - base[k]
         end
         tdsum[T,t] = rewards[T,t] - base[T]
         #mat_adv[:,t] = discount_sum(tdsum[:,t], DT(gae*gamma), DT(0.0))
         discount_sum!(mat_adv, tdsum, t, DT(gae*gamma), DT(0.0))
      end
   end
end

# solver algorithms
function cpu_cg_solve{T<:AbstractFloat}(fim::Matrix{T}, # ((n+1)*m + m) X T*numT
                                        vpg::Vector{T}, # ((n+1)*m + m)
                                        cg_iters::Integer=10,
                                        reg::Float64=1e-4,
                                        tol::Float64=1e-10)
   # Initialize cg variables
   r = copy(vpg)
   p = copy(vpg)
   x = zeros(T, size(vpg))  # I want x to be same shape as vpg but full of zeros
   rdr = dot(vpg, vpg)  # outputs a scalar
   z = fim*p
   #z += p*reg

   iters = 1
   for i=1:cg_iters
      v = rdr/dot(p, z)      # scalar

      x .+= v.*p
      r .-= v.*z

      rdr_new = dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      iters = i
      if rdr < tol
         break   # this should break the for loop
      end

      p .*= ratio
      p .+= r

      #println(norm(z))
      A_mul_B!(z, fim, p)
   end
   println("    Used $iters cojugate-gradient iterations.");
   return x
end

function hvp_cg_solve{DT<:AbstractFloat}(hvpfim,
                                         vpg::Vector{DT}, # ((n+1)*m + m)
                                         cg_iters::Integer=10,
                                         reg::Float64=1e-4,
                                         tol::Float64=1e-10)
   # Initialize cg variables
   r = copy(vpg)
   p = copy(vpg)
   x = zeros(DT, size(vpg))
   rdr = dot(vpg, vpg)  # outputs a scalar
   z = zeros(DT, size(vpg))
   hvpfim(z, p)

   iters = 1
   for i=1:cg_iters
      v = rdr/dot(p, z)      # scalar

      x += v*p
      r -= v*z
      #@. x += v*p
      #@. r -= v*z

      rdr_new = dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      iters = i
      if rdr < tol
         break   # this should break the for loop
      end

      @. p = r + ratio*p
      #p[:] .*= ratio
      #p[:] += r

      #println(norm(z))
      hvpfim(z, p)
   end
   println("    Used $iters HVP-cg iterations.");
   return x
end

end


text:
module EyeFix

export IVT, IDT, GMM, kmeans

include("IVT.jl")
include("IDT.jl")
include("HMM.jl")
include("GMM.jl")
include("kmeans.jl")
include("common.jl")

end # module


text:
function increaseindex!(h::Dict{K,Int}, key::K) where K
    index = Base.ht_keyindex2!(h, key)
    if index > 0
        h.age += 1
        #@inbounds v0 = h.vals[index]
        @inbounds h.keys[index] = key
        @inbounds h.vals[index] += 1
    else
        @inbounds Base._setindex!(h, 1, key, -index)
    end
    return
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
# Julia wrapper for header: /usr/include/libavfilter/buffersrc.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_buffersrc_buffer,
    av_buffersrc_write_frame

function av_buffersrc_buffer(s,buf)
    ccall((:av_buffersrc_buffer,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{AVFilterBufferRef}),s,buf)
end

function av_buffersrc_write_frame(s,frame)
    ccall((:av_buffersrc_write_frame,libavfilter),Cint,(Ptr{AVFilterContext},Ptr{AVFrame}),s,frame)
end


text:
InputFolder = './Images/FISH/';
OutputFolder = './Results/Images/FISH_nuc/';

@iA = '*_C00*.tif';

@fxg_mLapThrBinWatTiles [iA] > [L];
params.GRad = 7;
params.Thr = 0.02;
params.GaussianD = 2;
params.MinArea = 150;
/endf

@fxm_mClearBorders [L] > [L2];
params = [];
/endf

/show iA > L2;
/keep L2 > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
#=
Copyright 2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

const CT = ChrBase.CaseTables

function _check_char(ch, tab)
    cmp = ch & ~0x001f
    for (c, mask) in tab
        cmp < c && break
        cmp == c && return ((mask%Int64) << (62 - ((ch & 0x1f)<<1))) >> 62
    end
    0
end

function _calc_len(pnt, fin, len, mask, tab)
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0xf0
            if ch > 0x7f
                c16 = ch < 0xe0 ? get_utf8_2byte(pnt += 1, ch) : get_utf8_3byte(pnt += 2, ch)
                ((mask >>> (c16 >>> 9)) & 1) != 0 && (len += _check_char(c16, tab))
            end
            pnt += 1
        else
            pnt += 4 # no 4 byte changes size in UTF-8
        end
    end
    len
end

# These are more complex case folding functions, and maybe belong in a separate UTF8Str.jl package

# Note: these only check for cases in Unicode where 2 byte sequences
# could expand to 3 byte sequences.  In the standard Unicode tables,
# that is the only expansion that occurs for upper or lower case

function _lower_utf8(beg, off, len)
    # Note, the final length may be larger or smaller
    fin = beg + len
    pnt = beg + off
    # Calculate final length
    len = _calc_len(pnt, fin, len, CT.ct.siz_l_flg, CT.sizvecl)
    buf, out = _allocate(UInt8, len)
    unsafe_copyto!(out, beg, off)
    out += off
    pnt = beg + off
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0x80
            set_codeunit!(out, ch + (_is_upper_a(ch) << 5))
            out += 1
        elseif ch < 0xc4
            ch = (ch << 6) | (get_codeunit(pnt += 1) & 0x3f)
            out = output_utf8_2byte!(out, ch + (_is_upper_l(ch) << 5))
        elseif ch < 0xe0
            # 2 byte
            c16 = get_utf8_2byte(pnt += 1, ch)
            if _can_lower_bmp(c16)
                c16 = _lower_bmp(c16)
                if c16 < 0x80
                    set_codeunit!(out, c16%UInt8)
                    out += 1
                elseif c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_2byte!(out, c16)
            end
        elseif ch < 0xf0
            # 3 byte
            c16 = get_utf8_3byte(pnt += 2, ch)
            if _can_lower_bmp(c16)
                c16 = _lower_bmp(c16)
                if c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_3byte!(out, c16)
            end
        else
            # 4 byte
            c32 = get_utf8_4byte(pnt += 3, ch)
            c32 <= 0x1ffff && _can_lower_slp(c32) && (c32 = _lower_slp(c32))
            out = output_utf8_4byte!(out, c32)
        end
        pnt += 1
    end
    Str(UTF8CSE, buf)
end

function _upper_utf8(beg, off, len)
    fin = beg + len
    pnt = beg + off
    # Note, the final length may be larger or smaller
    len = _calc_len(pnt, fin, len, CT.ct.siz_u_flg, CT.sizvecu)

    buf, out = _allocate(UInt8, len)
    unsafe_copyto!(out, beg, off)
    out += off
    while pnt < fin
        ch = get_codeunit(pnt)
        if ch < 0x80
            set_codeunit!(out, ch - (_is_lower_a(ch)<<5))
            out += 1
        elseif ch < 0xc4
            ch = (ch << 6) | (get_codeunit(pnt += 1) & 0x3f)
            out = (ch == 0xdf
                   ? output_utf8_3byte!(out, 0x1e9e)
                   : output_utf8_2byte!(out, _upper_bmp(ch)))
        elseif ch < 0xe0
            # 2 byte
            c16 = get_utf8_2byte(pnt += 1, ch)
            if _can_upper_bmp(c16)
                c16 = _upper_bmp(c16)
                # Check if still 2 byte, could increase to 3 byte, or decrease to 1 byte
                if c16 < 0x80
                    set_codeunit!(out, c16%UInt8)
                    out += 1
                elseif c16 < 0x800
                    out = output_utf8_2byte!(out, c16)
                else
                    out = output_utf8_3byte!(out, c16)
                end
            else
                out = output_utf8_2byte!(out, c16)
            end
        elseif ch < 0xf0
            # 3 byte
            c16 = get_utf8_3byte(pnt += 2, ch)
            if _can_upper_bmp(c16)
                c16 = _upper_bmp(c16)
                # Check if still 3 byte, uppercase form could drop to 2 byte
                out = (c16 < 0x800) ? output_utf8_2byte!(out, c16) : output_utf8_3byte!(out, c16)
            else
                out = output_utf8_3byte!(out, c16)
            end
        else
            # 4 byte
            c32 = get_utf8_4byte(pnt += 3, ch)
            c32 <= 0x1ffff && _can_upper_slp(c32) && (c32 = _upper_slp(c32))
            out = output_utf8_4byte!(out, c32)
        end
        pnt += 1
    end
    Str(UTF8CSE, buf)
end

function lowercase(str::MaybeSub{S}) where {C<:UTF8CSE,S<:Str{C}}
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            (ch < 0x80
             ? _is_upper_a(ch)
             : (ch < 0xc4
                ? _is_upper_l((ch << 6) | (get_codeunit(pnt += 1) & 0x3f))
                : _can_lower_ch(ch >= 0xf0
                                ? get_utf8_4byte(pnt += 3, ch)
                                : (ch < 0xe0
                                   ? get_utf8_2byte(pnt += 1, ch)
                                   : get_utf8_3byte(pnt += 2, ch))%UInt32))) &&
                            return _lower_utf8(beg, prv-beg, ncodeunits(str))
            pnt += 1
        end
        str
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:UTF8CSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if ch <= 0x7f
            _is_upper_a(ch) || return str
            buf, out = _allocate(UInt8, len)
            len > 1 && unsafe_copyto!(out, pnt, len)
            set_codeunit!(out, ch + 0x20)
        elseif ch < 0xf0
            siz = 2 + (ch >= 0xe0)
            c16 = (ch < 0xe0 ? get_utf8_2byte(pnt + 1, ch) : get_utf8_3byte(pnt + 2, ch))
            _can_lower_bmp(c16) || return str
            # BMP characters might change size on lowercasing
            #print("lf: $c16")
            c16 = _lower_bmp(c16)
            ns = 1 + (c16 >= 0x80) + (c16 >= 0x800)
            buf, out = _allocate(UInt8, len - siz + ns)
            len > siz && unsafe_copyto!(out + ns, pnt + siz, len - siz)
            #println(" => $c16 len=$len ns=$ns siz=$siz out=$out pnt=$pnt $buf")
            (c16 < 0x80
             ? set_codeunit!(out, c16%UInt8)
             : (c16 < 0x800
                ? output_utf8_2byte!(out, c16)
                : output_utf8_3byte!(out, c16)))
            #println(" => $buf")
        else
            c32 = get_utf8_4byte(pnt + 3, ch)
            (c32 <= 0x1ffff && _can_lower_slp(c32)) || return str
            buf, out = _allocate(UInt8, len)
            len > 4 && unsafe_copyto!(out + 4, pnt + 4, len - 4)
            output_utf8_4byte!(out, _lower_slp(c32))
        end
        Str(C, buf)
    end
end

# Check if can be uppercased
@inline function _check_uppercase(ch, pnt)
    # ch < 0xc2 && return false (not needed, validated UTF-8 string)
    cont = get_codeunit(pnt)
    ch == 0xc3 ? ((cont > 0x9e) & (cont != 0xb7)) : (cont == 0xb5)
end

function uppercase(str::MaybeSub{S}) where {C<:UTF8CSE,S<:Str{C}}
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            ch = get_codeunit(pnt)
            prv = pnt
            (ch < 0x80
             ? _is_lower_a(ch)
             : (ch > 0xc3
                ? _can_upper_ch(ch >= 0xf0
                                ? get_utf8_4byte(pnt += 3, ch)
                                : (ch < 0xe0
                                   ? get_utf8_2byte(pnt += 1, ch)
                                   : get_utf8_3byte(pnt += 2, ch))%UInt32)
                : _check_uppercase(ch, pnt += 1))) &&
                    return _upper_utf8(beg, prv-beg, ncodeunits(str))
            pnt += 1
        end
        str
    end
end

function uppercase_first(str::MaybeSub{S}) where {C<:UTF8CSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if ch <= 0x7f
            _is_lower_a(ch) || return str
            buf, out = _allocate(UInt8, len)
            len > 1 && unsafe_copyto!(out, pnt, len)
            set_codeunit!(out, ch - 0x20)
        elseif ch < 0xf0
            siz = 2 + (ch >= 0xe0)
            c16 = (ch < 0xe0 ? get_utf8_2byte(pnt + 1, ch) : get_utf8_3byte(pnt + 2, ch))
            (tc = _titlecase(c16)) == c16 && return str
            # BMP characters might change size on titlecasing
            ns = 1 + (tc >= 0x80) + (tc >= 0x800)
            buf, out = _allocate(UInt8, len - siz + ns)
            len > siz && unsafe_copyto!(out + ns, pnt + siz, len - siz)
            (tc < 0x80
             ? set_codeunit!(out, tc%UInt8)
             : (tc < 0x800
                ? output_utf8_2byte!(out, tc)
                : output_utf8_3byte!(out, tc)))
        else
            c32 = get_utf8_4byte(pnt + 3, ch)
            (c32 <= 0x1ffff && _can_upper_slp(c32)) || return str
            buf, out = _allocate(UInt8, len)
            len > 4 && unsafe_copyto!(out + 4, pnt + 4, len - 4)
            output_utf8_4byte!(out, _upper_slp(c32))
        end
        Str(C, buf)
    end
end


text:
using Base.Iterators
using Distributions
using Statistics
using SparseArrays
using Random

export Grid, State, get_states, dims, num_states, neighbors, random_state,
  target_dynamics, RangingSensor, CoverageSensor, generate_observation,
  compute_likelihoods, transition_matrix, default_num_targets, continuous

abstract type StateSpace end
# methods
#
# get_states

# State and trajectory objects for convenience
const State = Tuple{Int64,Int64}
const Trajectory = Vector{State}

struct Grid <: StateSpace
  width::Int64
  height::Int64
  states::Array{State,2}
  transition_matrix::SparseMatrixCSC{Float64,Int64}

  # We will precompute some of the large objects that we use frequently
  function Grid(width, height)
    x = new(width, height, get_states(width, height),
            sparse([], [], Float64[], width*height, width*height)
           )
    x.transition_matrix .= generate_transition_matrix(x)

    x
  end
end

# Construct default grids according to the number of robots
const grid_cells_per_robot = 12.5
function Grid(;num_robots)
  grid_size = round(Int64, sqrt(grid_cells_per_robot * num_robots))
  Grid(grid_size, grid_size)
end

# default number of targets as a function of the number of robots
const targets_per_robot = 1
default_num_targets(num_robots) = round(Int64, num_robots * targets_per_robot)

get_states(width::Real, height::Real) = collect(product(1:width, 1:height))
get_states(g::Grid) = g.states
dims(g::Grid) = (g.width, g.height)
num_states(g::Grid) = length(g.states)

in_bounds(g::Grid, state) = all(x->in(x[1],1:x[2]), zip(state, dims(g)))

state_to_index(g::Grid, s) = s[1] + (g.width) * (s[2] - 1)
function index_to_state(g::Grid, x)
  b = div(x-1, g.width) + 1
  a = x - (b - 1) * g.width
  (a, b)
end

# Wrap states back into the grid
wrap(g::Grid, s) = mod.(s .- 1, dims(g)) .+ 1

# Produces the out neighbors of the transition graph based on a grid model
# We preallocate the array with the max number of neighbors
neighbors_buffer() = Vector{State}(undef, 5)
function neighbors(g::Grid, state; buffer=neighbors_buffer())
  dirs = ((1,0), (0,1))
  offsets = (-1, 1)

  resize!(buffer, 0)

  for dir in dirs, offset in offsets
    candidate = state .+ offset .* dir

    if in_bounds(g, candidate)
      push!(buffer, candidate)
    end
  end

  push!(buffer, state)

  buffer
end

# Returns a left stochastic matrix A so that posterior = A * prior
# This should only be run while initializing the grid
function generate_transition_matrix(g::Grid)
  rows = Int64[]
  columns = Int64[]
  weights = Float64[]

  # Push uniform transition probabilities for each state
  for state in get_states(g)
    ns = neighbors(g, state)

    source = state_to_index(g, state)
    weight = 1 / length(ns)

    for neighbor in ns
      dest = state_to_index(g, neighbor)

      push!(columns, source)
      push!(rows, dest)
      push!(weights, weight)
    end
  end

  size = length(get_states(g))
  sparse(rows, columns, weights, size, size)
end
transition_matrix(grid::Grid) = grid.transition_matrix

random_state(g::Grid; rng=Random.GLOBAL_RNG) = sample(rng, get_states(g))
function target_dynamics(g::Grid, s; rng=Random.GLOBAL_RNG,
                         buffer=neighbors_buffer())

  sample(rng, neighbors(g, s, buffer=buffer))
end

#################
# Ranging sensors
#################

# variance of observations is: constant + scaling * norm_squared
struct RangingSensor
  variance_constant::Float64
  variance_scaling_factor::Float64
  range_limit::Float64

  function RangingSensor(;variance_constant=0.25,
                          variance_scaling_factor=0.5,
                          range_limit=10)
    new(variance_constant, variance_scaling_factor, range_limit)
  end
end

saturate_range(r::RangingSensor, distance) = min(r.range_limit, distance)

# Sensor forward distribution variance and standard deviation
sensor_variance(r::RangingSensor, distance) =
  (r.variance_constant
   + r.variance_scaling_factor * saturate_range(r, distance)^2)

sensor_stddev(r::RangingSensor, distance) = sqrt(sensor_variance(r, distance))

# Mean value for the sensor distribution
# * Saturated distance
sensor_mean(g::Grid, r::RangingSensor, a, b) = saturate_range(r, norm(a .- b))

# sample a ranging observation
function generate_observation(g::Grid, r::RangingSensor, a, b;
                              rng=Random.GLOBAL_RNG)
  m = sensor_mean(g, r,a,b)
  m + randn(rng) * sensor_stddev(r, m)
end

const normal_lookup_table = NormalLookup(increment=0.001, max=4.0)

likelihoods_buffer(states) = Array{Float64}(undef, length(states))

# Computes (non-normalized) likelihoods of data
function compute_likelihoods(robot_state, target_states, sensor::RangingSensor,
                             grid::Grid,
                             range::Real;
                             buffer = likelihoods_buffer(target_states)
                            )

  resize!(buffer, length(target_states))

  for (ii, target_state) in enumerate(target_states)
    distance = sensor_mean(grid, sensor, robot_state, target_state)

    buffer[ii] = pdf(Normal(distance, sensor_stddev(sensor, distance)), range)
  end

  buffer
end

##################
# Coverage sensors
##################

struct CoverageSensor
  range_limit::Float64

  function CoverageSensor(; range_limit=2)
    new(range_limit)
  end
end

# produces an array of continuous ranges of states in the intput trajectory
# (as the robot may cros over the grid)
continuous(a::State, b::State) = sum(abs.(a .- b)) <= 1
function continuous(states::Trajectory)
  for ii = 2:length(states)
    if !continuous(states[ii-1], states[ii])
      return false
    end
  end

  true
end
function continuous(grid::Grid, states::Trajectory)
  for ii = 2:length(states)
    if ! in(states[ii], neighbors(grid, states[ii-1]))
      return false
    end
  end

  true
end

function continuous_ranges(states::Trajectory)
  ret = Vector{State}[]

  # Ends of ranges are states that are discontinuous with the following
  range_ends = findall(1:length(states)-1) do ii
    !continuous(states[ii], states[ii+1])
  end
  push!(range_ends, length(states))

  # Push the ranges of states
  start = 1
  for range_end in range_ends
    push!(ret, states[start:range_end])
    start = range_end + 1
  end

  ret
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using Test
using Arpack: eigs
using ExactDiagonalization.FED
using SparseArrays: SparseMatrixCSC
using QuantumLattices
using QuantumLattices.Mathematics.AlgebraOverFields: idtype

@testset "GBasis" begin
    basis = GBasis{UInt}(2)
    @test kind(basis) == kind(typeof(basis)) == :g
    @test idtype(basis) == idtype(typeof(basis)) == Int
    @test dimension(basis) == 2^2
    @test id(basis) == 2
    @test basis[4] == 3
    @test searchsortedfirst(basis, 3) == 4
    @test string(basis) == "GBasis(2):\n  0\n  1\n  10\n  11\n"
    @test repr(basis) == "GBasis(2)"
end

@testset "SparseMatrixCSC" begin
    lattice = Lattice("L2P", [Point(PID(1), [0.0]), Point(PID(2), [1.0])])
    config = Config{Fock{:f}}(pid->Fock{:f}(atom=pid.site%2, norbital=1, nspin=2, nnambu=2), lattice.pids)
    table = Table(config, usualfockindextotuple)
    t = Hopping(:t, 1.0, 1)
    U = Hubbard(:U, 1.0)
    μ = Onsite(:μ, -0.5)
    g = Generator((t, U, μ), Bonds(lattice), config, half=false, table=table)
    m = SparseMatrixCSC(expand(g), GBasis{UInt}(4), g.table)
    @test eigs(m, which=:SR, ritzvec=false)[1][1] ≈ -2.5615528128088303
end


text:
module BISC195Labs

export nwscore

include("needleman_wunch.jl")

end


text:
using Base.Test

import IJulia: helpcode, error_content

@test !haskey(Docs.keywords, :+)
@test "Base.Docs.@repl +" == helpcode("+")

@test haskey(Docs.keywords, :import)
@test """eval(:(Base.Docs.@repl \$(Symbol("import"))))""" == helpcode("import")

content = error_content(UndefVarError(:a))
@test "UndefVarError" == content["ename"]

@test haskey(IJulia.docdict("sum"), "text/plain")


text:
###############################
# Global state for annotation #
###############################

# TODO: Support multi-threaded usage by storing a global array with an instance
# of this annotation state for each thread
curr_delayed_task = DelayedTask()

function set_task(t::DelayedTask)
    global curr_delayed_task
    curr_delayed_task = t
end

function get_task()
    global curr_delayed_task
    curr_delayed_task
end

function finish_task()
    global curr_delayed_task
    curr_delayed_task = DelayedTask()
end

function get_pa()
    global curr_delayed_task
    last(curr_delayed_task.pa_union)
end

function finish_pa()
    global curr_delayed_task
    push!(curr_delayed_task.pa_union, PartitionAnnotation())
end

#################################################
# Setting up sample properties in beforecompute #
#################################################

# These functions are necessary because when we construct PT's to assign to
# futures for annotated code (especially code in BDF.jl), we rely on statistics
# of the samples of the futures to construct the PTs. But computing statistics
# is expensive and we only want to compute it for the keys/columns that are
# actually used for grouping/sorting/etc. So we propagate this through the
# `:groupingkeys` sample property of different futures (each future has a
# set of sample properties).
#
# These `keep_*` functions specify how `statistics` and `groupingkeys` sample
# properties should propagate. These functions are called inside of a
# `partitioned_using` so that they are called in a forward pass through all the
# tasks (code regions) detected and also in a backward pass to properly
# propagate the `groupingkeys` sample property. For example, a groupby may only
# happen in the very last task indicating that a particular column can be in
# the `groupingkeys` (the column used for the groupby) but we need to specify
# that that column could have been used for grouping the data throughout (if
# all tasks are amenable like that).
#
# Of course, columns can be dropped or
# renamed and so that makes things a bit challenging in properly propagating
# information about sample properties and that's why we have a lot of these
# keep_* functions.

function partitioned_using(handler::Function)
    global curr_delayed_task
    curr_delayed_task.partitioned_using_func = handler
end

function keep_all_sample_keys(participants::AbstractFuture...; drifted = false)
    # This can be use on an input and output that share the same keys that they
    # can be grouped by.

    # Take the union of discovered grouping keys. Grouping keys are discovered
    # through calls to `keep_keys` in functions where we know that we need
    # to allow grouping by certain keys such as sorting and join functions.
    groupingkeys = union([sample(p, :groupingkeys) for p in participants]...)
    statistics = if !drifted
        merge([sample(p, :statistics) for p in participants]...)
    else
        nothing
    end

    # Only allow keys that are actually in the keys of the participants
    # TODO: Maybe replace this with a function that doesn't require calling
    # a potentially expensive function to iterate through _all_ columns.
    for p in participants
        p_keys = intersect(groupingkeys, sample(p, :keys))
        setsample!(p, :groupingkeys, p_keys)
        if !drifted
            setsample!(
                p,
                :statistics,
                Dict(
                    key => statistic for
                    (key, statistic) in statistics if key in sample(p, :keys)
                ),
            )
        end
    end
end

function keep_all_sample_keys_renamed(old::AbstractFuture, new::AbstractFuture)
    for (old_key, new_key) in zip(sample(old, :keys), sample(new, :keys))
        if old_key in sample(old, :groupingkeys)
            setsample!(new, :groupingkeys, union(sample(new, :groupingkeys), [new_key]))
            setsample!(new, :statistics, new_key, sample(old, :statistics, old_key))
        end
        if new_key in sample(new, :groupingkeys)
            setsample!(old, :groupingkeys, union(sample(old, :groupingkeys), [old_key]))
            setsample!(old, :statistics, old_key, sample(new, :statistics, new_key))
        end
    end
end

function keep_sample_keys_named(
    participants::Pair{<:AbstractFuture,<:Any}...;
    drifted = false,
)
    # `participants` maps from futures to lists of key names such that all
    # participating futures have the same sample properties for the keys at
    # same indices in those lists
    participants = [
        (participant, Symbol.(to_vector(key_names))) for
        (participant, key_names) in participants
    ]
    nkeys = length(first(participants)[2])
    for i = 1:nkeys
        # Copy over allowed grouping keys
        for (p, keys) in participants
            p_key = keys[i]
            setsample!(p, :groupingkeys, union(sample(p, :groupingkeys), [p_key]))
        end

        # Copy over statistics if they haven't changed
        if !drifted
            key_statistics =
                merge([sample(p, :statistics, keys[i]) for (p, keys) in participants]...)
            for (p, keys) in participants
                setsample!(p, :statistics, key_statistics)
            end
        end
    end
end

# NOTE: For operations where join rate and data skew might be different, it is
# assumed that they are the same. For example, column vectors from the result
# of a join should have the same sample rate and the same data skew.

keep_sample_keys(keys, participants::AbstractFuture...; drifted = false) = begin
    keys = Symbol.(to_vector(keys))
    keep_sample_keys_named([p => keys for p in participants]..., drifted = drifted)
end

# This is useful for workloads that involve joins where the sample rate is
# diminished quadratically for each joinv
keep_sample_rate(fut::AbstractFuture, relative_to::AbstractFuture...) =
    setsample!(fut, :rate, prod([sample(r, :rate) for r in relative_to]))

###############################
# Using samples to assign PTs #
###############################

function partitioned_with(handler::Function)
    # NOTE: `partitioned_with` should mainly include calls to `partition`;
    # calls to `mut` should be made outside
    global curr_delayed_task
    curr_delayed_task.partitioned_with_func = handler
end

function pt(
    args::Union{AbstractFuture,PartitionType,PartitionTypeComposition,Vector}...;
    kwargs...,
)
    pa = get_pa()

    # Extract PT and args to assign the PT to from given arguments
    futs, sharedptype = args[1:end-1], last(args)

    if length(futs) > 1 && sharedptype isa Vector && length(sharedptype) > 1
        throw(
            ArgumentError(
                "Multiple partition types cannot be applied to multiple futures at once",
            ),
        )
    end

    # Assign PT to each future
    for fut in futs
        # Start new PA if this assignment would overwrite one for the current
        # PA. When we start a new PA, we append the old one the PA union for
        # the task being currently constructed.
        fut = convert(Future, fut)
        if fut.value_id in keys(pa.partitions.pt_stacks)
            finish_pa()
            pa = get_pa()
        end

        # Take a copy so that we can mutate
        ptype = deepcopy(sharedptype)

        # Put in PT
        if ptype isa Vector
            for pty in ptype
                pt(fut, PartitionTypeComposition([pty]); kwargs...)
            end
        elseif ptype isa PartitionType
            pt(fut, PartitionTypeComposition([ptype]); kwargs...)
        elseif ptype isa Vector{PartitionType} && length(ptype) == 1
            pt(fut, PartitionTypeComposition(ptype); kwargs...)
        elseif ptype isa PartitionTypeComposition
            # Handle constraints that have been delayed till PT assignment
            for pty in ptype.pts
                for (i, constraint) in enumerate(pty.constraints.constraints)
                    if constraint isa Function
                        pty.constraints.constraints[i] = constraint(fut)
                    end
                end
            end

            # Add to PAs information about how partitions are produced
            pa.partitions.pt_stacks[fut.value_id] = ptype

            # Add to PAs information about how partitioning is constrained
            for c in get(kwargs, :constrain, [])
                push!(pa.constraints.constraints, c)
            end

            # Handle `match`, `on` in keyword arguments
            if :match in keys(kwargs) && !isnothing(kwargs[:match])
                to_match_with = to_vector(kwargs[:match])
                if :on in keys(kwargs) && !isnothing(kwargs[:on])
                    for to_match_on in to_vector(get(kwargs, :on, []))
                        push!(
                            pa.constraints.constraints,
                            MatchOn(to_match_on, fut, to_match_with...),
                        )
                    end
                else
                    push!(pa.constraints.constraints, Match(fut, to_match_with...))
                end
            end

            if :cross in keys(kwargs) && !isnothing(kwargs[:cross])
                to_cross = to_vector(kwargs[:match])
                push!(pa.constraints.constraints, Cross(to_cross...))
            end

            # TODO: Implement support for other constraints in kwargs
        else
            throw(
                ArgumentError(
                    "Expected partition type (PT) or a composition or union of PTs",
                ),
            )
        end
    end
end

# NOTE: `mutated` should not be used inside of `partitioned_with` or
# `partitioned_using`

mutated(f::AbstractFuture) = mutated(f => f)
mutated(ff::Pair{<:AbstractFuture,<:AbstractFuture}) = mutated(first(ff), last(ff))

function mutated(old::AbstractFuture, new::AbstractFuture)
    global curr_delayed_task
    old = convert(Future, old)
    # if haskey(curr_delayed_task.value_names, old.value_id)
    #     finish_pa()
    curr_delayed_task.mutation[old] = convert(Future, new)
end

#################################################
# Macro for wrapping the code region to offload #
#################################################

function apply_mutation(mutation::Dict{Future,Future})
    for (old, new) in mutation
        if old != new
            # Apply the mutation by setting the value ID of the old future the
            # value ID of the new one. That way, the user can continue using
            # the old future as if it were mutated but it will be having a
            # different value ID.

            # Swap references in `futures_on_client` if either side of the
            # mutation is on the client
            futures_on_client = get_job().futures_on_client
            if old.value_id in keys(futures_on_client) &&
               new.value_id in keys(futures_on_client)
                futures_on_client[new.value_id], futures_on_client[old.value_id] =
                    futures_on_client[old.value_id], futures_on_client[new.value_id]
            elseif old.value_id in keys(futures_on_client)
                futures_on_client[new.value_id] = futures_on_client[old.value_id]
                delete!(futures_on_client, old.value_id)
            elseif new.value_id in keys(futures_on_client)
                futures_on_client[old.value_id] = futures_on_client[new.value_id]
                delete!(futures_on_client, new.value_id)
            end

            # Swap other fields of the `Future`s and their locations
            job_locations = get_job().locations
            old.value,
            old.value_id,
            old.mutated,
            old.stale,
            job_locations[old.value_id],
            new.value,
            new.value_id,
            new.mutated,
            new.stale,
            job_locations[new.value_id] = new.value,
            new.value_id,
            new.mutated,
            new.stale,
            job_locations[new.value_id],
            old.value,
            old.value_id,
            old.mutated,
            old.stale,
            job_locations[old.value_id]
        end
    end
end

invert(mutation::Dict{Future,Future}) = Dict(new => old for (old, new) in mutation)

function partitioned_using_modules(m...)
    global curr_delayed_task
    union!(curr_delayed_task.used_modules, m)
end

macro partitioned(ex...)
    res = quote end

    # Load in variables and code from macro
    variables = [esc(e) for e in ex[1:end-1]]
    variable_names = [string(e) for e in ex[1:end-1]]
    code = ex[end]

    # # Expand splatted variables
    # # if any(endswith(name, "...") for name in variable_names)
    # if true
    #     new_variables = []
    #     new_variable_names = []

    #     # Iterate through variables specified through the annotation
    #     for (variable, name) in zip(variables, variable_names)
    #         splatted = Meta.parse(name)
    #         println(__module__.eval(quote $variable end))
    #         # println(eval(quote $(variable) end))
    #     end

    #     # Iterate through variables specified through the annotation
    #     for (variable, name) in zip(variables, variable_names)
    #         # Variables with names ending with ... need to be expanded
    #         if endswith(name, "...")
    #             # Get information about the variable being splatted
    #             expanded_expr = Meta.parse(name[1:end-3])
    #             println(typeof(expanded_expr))
    #             println(:([$(esc(ex[2]))]))
    #             println(Base.eval(__module__, :([$(esc(ex[2]))])))
    #             expanded::Vector{AbstractFuture} = Base.eval(__module__, esc(expanded_expr))

    #             # Append to the variables and variable names to be loaded
    #             # into the string
    #             expanded_variables = []
    #             for (i, e) in enumerate(expanded)
    #                 # Appemnd to the variables and variable names to be used.
    #                 # We use a randomly-suffixed name for each variable in the
    #                 # expansion.
    #                 new_variable_name = name[1:end-3] * "_" * randstring(4)
    #                 new_variable = Meta.parse(new_variable_name)
    #                 push!(new_variable_names, new_variable_name)
    #                 push!(expanded_variables, new_variable)

    #                 # Then, in the code that this macro compiles to, we
    #                 # construct these randomly-suffixed variables to reference
    #                 # the appropriate future.
    #                 res = quote
    #                     $res
    #                     $new_variable = $expanded_expr[$i]
    #                 end
    #             end
    #             append!(new_variables, expanded_variables)

    #             # Append to the code so that the variable can be used as-is in
    #             # the code regiony
    #             code = quote
    #                 $expanded_expr = [$(expanded_variables...)]
    #                 $code
    #             end
    #         else
    #             push!(new_variables, variable)
    #             push!(new_variable_names, name)
    #         end
    #     end
    #     variables = new_variables
    #     variables_names = new_variable_names
    # end

    # Assign samples to variables used in annotated code
    assigning_samples = [
        quote
            unsplatted_future = unsplatted_futures[$i]
            # println(get_job().locations)
            # println(typeof(unsplatted_future))
            # unsplatted_variable_names = [$(variable_names...)]
            # println(unsplatted_variable_names[$i])
            $variable =
                unsplatted_future isa Tuple ? sample.(unsplatted_future) :
                sample(unsplatted_future)
        end for (i, variable) in enumerate(variables)
    ]

    # Re-assign futures to variables that were used in annotated code
    # TODO: Ensure that it is okay for different `quote...end` blocks to refer
    # to the same variable name. They shouldn't have different gensym-ed
    # variable names in the macro expansion.
    reassigning_futures = [
        quote
            $variable = unsplatted_futures[$i]
        end for (i, variable) in enumerate(variables)
    ]

    quote
        # $res

        # Convert arguments to `Future`s if they aren't already
        unsplatted_futures = [$(variables...)]
        # [unsplatted_future isa Base.Vector ? convert.(Future, unsplatted_future : convert(Future))]
        # for unsplatted_future in unsplatted_futures
        #     println(unsplatted_future)
        #     println(typeof(unsplatted_future))
        # end
        splatted_futures::Vector{Future} = []
        for unsplatted_future in unsplatted_futures
            if unsplatted_future isa Tuple
                for uf in unsplatted_future
                    push!(splatted_futures, convert(Future, uf))
                end
            else
                push!(splatted_futures, convert(Future, unsplatted_future))
            end
        end
        # splatted_futures::Vector{Future} = [
        #     (
        #         unsplatted_future isa Tuple ?
        #         Tuple([convert(Future, uf) for uf in unsplatted_future]) :
        #         convert(Future, unsplatted_future)
        #     ) for unsplatted_future in unsplatted_futures
        # ]
        # splatted_futures::Vector{AbstractFuture} = vcat(unsplatted_futures...) .|> x->convert.(Future,x)

        # TODO: Allow for any Julia object (even stuff that can't be converted
        # to `Future`s) to be passed into an @partitioned and by default have
        # it be replicated across all workers and batches
        # for pa in get_pa_union()
        #     for fut in futures
        #         if !(fut.value_id in keys(pa.partitions.pt_stacks))
        #             pa.partitions.pt_stacks[fut.value_id] = [Replicated()]
        #             # TODO: Add necessary default constraints so this
        #             # replication spans what everything else spans
        #         end
        #     end
        # end

        # If any sample computation fails, before we rethrow
        # the error (so that it is displayed in a notebook or crashes a
        # script/Lambda function that is running the job) we ensure that
        # we haven't recorded a faulty task or messed up the state in any way.

        # Fill in task with code and value names pulled using the macror
        unsplatted_variable_names = [$(variable_names...)]
        splatted_variable_names = []
        task = get_task()
        # Get code to initialize the unsplatted variable in the code region
        # TODO: Generate code in codegen that puts each code region in a
        # seperate function (where we reuse functions with the hash of the
        # function body) so that we don't have scoping-related bugs
        task.code = ""
        for (variable, unsplatted_variable_name) in
            zip(unsplatted_futures, unsplatted_variable_names)
            task.code *= "$unsplatted_variable_name = "
            if variable isa Tuple
                task.code *= "["
                for (i, v) in enumerate(variable)
                    splatted_variable_name = unsplatted_variable_name * "_$i"
                    push!(splatted_variable_names, splatted_variable_name)
                    task.code *= "$splatted_variable_name, "
                end
                task.code *= "]\n"
            else
                push!(splatted_variable_names, unsplatted_variable_name)
                task.code *= "$unsplatted_variable_name\n"
            end
        end
        task.code *= $(string(code))
        task.value_names = [
            (fut.value_id, var_name) for (fut, var_name) in
            # zip(futures, [$(variable_names...)])
            zip(splatted_futures, splatted_variable_names)
        ]
        # task = DelayedTask(
        #     ,
        #     ,
        #     # Dict(
        #     #     fut.value_id =>
        #     #         if get_mutated(fut.value_id)
        #     #             "MUT"
        #     #         else
        #     #             "CONST"
        #     #         end for fut in futures
        #     # ),
        #     Dict(),
        #     get_pa_union(),
        # )

        # Set `mutated` field of the `Future`s that have been mutated. This is
        # to ensure that future calls to `evaluate` on those `Future`s with
        # `mutated=true` and _only_ those `Future`s will result in an actual
        # evaluation
        for fut in splatted_futures
            if any((fut.value_id == m.value_id for m in values(task.mutation)))
                task.effects[fut.value_id] = "MUT"
            else
                task.effects[fut.value_id] = "CONST"
            end
        end

        # TODO: When mutating a value, ensure that the old future has a sample
        # of the old value and the new future of the new

        # Perform computation on samples
        # begin
        $(assigning_samples...)
        # y=10
        # x= (5,6,7,8,y,y)
        # # Store samples in variables
        # # let $(variables...) = [f isa Vector ? sample.(f) : sample(f) for f in unsplatted_futures]
        # # TODO: Generate code ro store smaples in varibles
        # samples = [(f isa Vector ? sample.(f) : sample(f)) for f in unsplatted_futures]
        # ($(variables...)) = samples

        # Run the actual code. We don't have to do any splatting here
        # because the variables already each contain either a single
        # future or a list of them.
        # TODO: Fix error
        # MethodError: Cannot `convert` an object of type Int32 to an object of type Symbol
        # Closest candidates are:
        #   convert(::Type{T}, ::PropertyDicts.PropertyDict) where T at /home/calebwin/.julia/packages/PropertyDicts/6Kqy6/src/PropertyDicts.jl:24
        #   convert(::Type{T}, ::LazyJSON.Object) where T at /home/calebwin/.julia/packages/LazyJSON/LpIGF/src/AbstractDict.jl:80
        #   convert(::Type{T}, ::T) where T at essentials.jl:205
        #   ...
        # Stacktrace:
        #   [1] merge(a::NamedTuple{(:dims,), Tuple{Colon}}, itr::Int32)
        #     @ Base ./namedtuple.jl:281
        try
            $(esc(code))
        catch
            # I
            finish_task()
            $(reassigning_futures...)
            rethrow()
        end

        # NOTE: We only update futures' state, record requests, update samples,
        # apply mutation _IF_ the sample computation succeeds. Regardless of
        # whether it succeeds, we make sure to clear the task (so that future
        # usage in the REPL or notebook can proceed, generating new tasks) and
        # reassign all variables.

        # Update mutated futures
        for fut in splatted_futures
            if any((fut.value_id == m.value_id for m in values(task.mutation)))
                fut.stale = true
                fut.mutated = true
            end
        end

        # Apply all delayed source and destination assignment. This will
        # perform any expensive sample collection that may require for example
        # an expensive scan of S3. This will record `RecordLocationRequest`s.
        for splatted_future in splatted_futures
            apply_sourced_or_destined_funcs(splatted_future)
        end

        # Record request to record task in backend's dependency graph and reset
        record_request(RecordTaskRequest(task))
        finish_task()

        # Move results from variables back into the samples. Also, update the
        # memory usage accordingly.
        # TODO: Determine if other sample properties need to be invalidated (or
        # updated) after modified by an annotated code region.
        for (f, value) in zip(unsplatted_futures, [$(variables...)])
            if f isa Tuple
                for (fe, ve) in zip(f, value)
                    setsample!(fe, ve)
                    setsample!(fe, :memory_usage, sample_memory_usage(ve))
                end
            else
                setsample!(f, value)
                setsample!(f, :memory_usage, sample_memory_usage(value))
            end
        end
        # end
        # end

        # This basically undoes `$(assigning_samples...)`.
        $(reassigning_futures...)

        # Make a call to `apply_mutation` to handle calls to `mut` like
        # `mutated(df, res)`
        apply_mutation(task.mutation)
    end
end

############################################################################
# Helper functions for compiling PAs to send as part of tasks in `compute` #
############################################################################

function duplicate_args(
    args::Vector{PartitionTypeReference},
    pa::PartitionAnnotation,
)::Vector{PartitionTypeReference}
    [(v, idx + div(length(pa.partitions.pt_stacks[v].pts), 2)) for
    # (v, idx + length(pa.partitions.pt_stacks[v].pts)) for
     (v, idx) in args]
end

function apply_default_constraints!(pa::PartitionAnnotation)
    # Add Cross constraints for all unconstrained PTs
    for (v, pt_stack) in pa.partitions.pt_stacks
        for i = 1:length(pt_stack.pts)
            # Check if in_cross_or_co
            in_cross_or_co = false
            for c in pa.constraints.constraints
                if (c.type == "CROSS" || c.type == "CO") && (v, i - 1) in c.args
                    in_cross_or_co = true
                elseif c.type == "CO_GROUP" && any((v, i - 1) in group for group in c.args)
                    in_cross_or_co = true
                end
            end

            # Add Cross constraint for those not constrained in any way
            if !in_cross_or_co
                push!(
                    pa.constraints.constraints,
                    PartitioningConstraintOverGroup("CROSS", [(v, i - 1)]),
                )
            end
        end
    end

    # Add Co constraint for all Cross-ed PTs

    # Find all co-partitioned PTs
    # inv: Every PT has been Cross-ed
    co_args = []
    co_group_args = []
    for c in pa.constraints.constraints
        if c.type == "CROSS"
            if length(c.args) == 1
                push!(co_args, c.args[1])
            elseif length(c.args) > 1
                push!(co_group_args, deepcopy(c.args))
            end
        end
    end
    if length(co_group_args) == 1 && length(co_args) > 0
        push!(co_group_args, [co_args[1]])
    end

    # Add constraints
    if length(co_args) > 0
        push!(pa.constraints.constraints, PartitioningConstraintOverGroup("CO", co_args))
    end
    if length(co_group_args) > 0
        push!(
            pa.constraints.constraints,
            PartitioningConstraintOverGroups("CO_GROUP", co_group_args),
        )
    end

    # TODO: Add MemoryUsage constraint by default which computes sample size or
    # defaults to 0
end

duplicated_constraints_for_batching(pc::PartitioningConstraints, pa::PartitionAnnotation) =
    begin
        # promote(
        #     Union{PartitioningConstraint, Function}[],
        new_pts = vcat(
            [
                if c.type == "CO" ||
                   c.type == "EQUAL" ||
                   c.type == "SEQUENTIAL" ||
                   c.type == "MATCH" ||
                   startswith(c.type, "MATCH_ON=")
                    [
                        deepcopy(c),
                        PartitioningConstraintOverGroup(c.type, duplicate_args(c.args, pa)),
                    ]
                elseif c.type == "CROSS" || startswith(c.type, "AT_MOST=")
                    # Note that with Cross constraints, the order of the
                    # arguments matters. But actually that doesnt matter.
                    # The scheduler will automaticcally ensure that the order
                    # of PTs in a PT stack is obeyed.
                    # ||
                    # c.type == "MATCH" || startswith(c.type, "MATCH_ON")
                    [
                        PartitioningConstraintOverGroup(
                            c.type,
                            [deepcopy(c.args); duplicate_args(c.args, pa)],
                        ),
                    ]
                elseif c.type == "CO_GROUP"
                    [
                        PartitioningConstraintOverGroups(
                            c.type,
                            [duplicate_args(group, pa) for group in c.args],
                        ),
                    ]
                elseif startswith(c.type, "SCALE_BY=")
                    # `ScaleBy` constraints are not duplicated. They must refer to
                    # only the first PT of the PT compositions they reference.
                    [c]
                else
                    []
                end for c in pc.constraints
            ]...,
        )
        # println()
        # println(new_pts)
        # )
        PartitioningConstraints(new_pts)
    end

function duplicate_for_batching!(pa::PartitionAnnotation)
    # Duplicate PT stacks
    for (v, pt_stack) in pa.partitions.pt_stacks
        # Copy over the PT stack
        second_half = deepcopy(pt_stack.pts)

        # # Don't duplicate parameters that have bang values
        # for pt in second_half
        #     for (k, v) in pt.parameters
        #         if v == "!"
        #             pop!(pt.parameters, k)
        #         end
        #     end
        # end

        # Append to form a compositions of PTs that is twic in length
        append!(pt_stack.pts, second_half)
    end

    # Duplicate annotation-level constraints for Co, Equal, Cross, AtMost, ScaleBy
    pa.constraints = duplicated_constraints_for_batching(pa.constraints, pa)
    # println(pa.constraints)

    # Add constraints for second half being Sequential and Match-ing the first
    # half
    for (v, pt_stack) in pa.partitions.pt_stacks
        for i = 1:div(length(pt_stack.pts), 2)
            dupi = i + div(length(pt_stack.pts), 2)

            # Add in `Sequential` and `Match` constraints for the duplicated
            # part of the PT composition
            push!(
                pa.constraints.constraints,
                PartitioningConstraintOverGroup("SEQUENTIAL", [(v, dupi - 1)]),
            )
            # Since we have duplicated all constraints, we don't need to
            # constrain both halves of the PT composition to match
            # push!(
            #     pa.constraints.constraints,
            #     PartitioningConstraintOverGroup("MATCH", [(v, i - 1), (v, dupi - 1)]),
            # )

            # Duplicate PT-level constraints
            pt_stack.pts[dupi].constraints =
                duplicated_constraints_for_batching(pt_stack.pts[dupi].constraints, pa)
        end
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using DAMMmodel
using Test

@testset "DAMMmodel Tests" begin
    # Testing DAMM.jl
	x = [18.0 0.35; 22.0 0.22] # Ind variables test
	p = (1e9, 64.0, 3.46e-8, 2.0e-3, 0.4, 0.0125, 1.0) # Parameters test
	@test typeof(DAMM(x, p)) == Vector{Float64} # DAMM test

	Ind_var = hcat([19.0, 22.0], [0.35, 0.22])
	@test lenght( DAMMfit(Ind_var, [2.2, 2.8], 0.4) ) == 7  # new DAMM test JF
end


text:
A = [
  1 2 3
  2 3 4
  3 4 5
]

B = [
  4 5 6
  5 6 7
  6 7 8
]

C = A + B

println(C)

max_stars_repo_name: mya-ake/matrix-playground

max_issues_repo_name: mya-ake/matrix-playground

max_forks_repo_name: mya-ake/matrix-playground


text:
abstract type Day10 <: AocDay end
day10 = Day10

function solve(day::Type{Day10})
    some_other_variable = 1
    println(day)
end


text:
using FixedEffects, Random, Statistics
Random.seed!(1234)
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
fes = [FixedEffect(id1), FixedEffect(id2)]
x = rand(N)

# simple problem
@time solve_residuals!(deepcopy(x), fes)
#   0.654833 seconds (1.99 k allocations: 390.841 MiB, 3.71% gc time)
@time solve_residuals!([x x x x], fes)
#   2.319452 seconds (7.71 k allocations: 620.016 MiB, 0.36% gc time)

# More complicated problem
N = 800000 # number of observations
M = 400000 # number of workers
O = 50000 # number of firms
Random.seed!(1234)
pid = rand(1:M, N)
fid = [rand(max(1, div(x, 8)-10):min(O, div(x, 8)+10)) for x in pid]
x = rand(N)
fes = [FixedEffect(pid), FixedEffect(fid)]
@time solve_residuals!(deepcopy(x), fes; double_precision = false)
# 2.239687 seconds (110.81 k allocations: 35.470 MiB)
@time solve_residuals!(deepcopy(x), fes; maxiter = 300)
# 3.681389 seconds (97.09 k allocations: 62.593 MiB)

@time solve_residuals!([x x x x], fes; double_precision = false)
# 5.253745 seconds (285.98 k allocations: 184.213 MiB, 0.55% gc time)

@time solve_residuals!([x x x x], fes; maxiter = 300)
# 9.889438 seconds (225.38 k allocations: 311.964 MiB, 0.33% gc time)

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects


text:
function toric_ed(A)
    d, n = size(A)
    @polyvar t[1:d] y[1:n] u[1:n]

    φ = map(j -> prod(i -> t[i]^A[i,j], 1:d), 1:n)
    Dφ = [differentiate(φ[i], t[j]) for i=1:n, j=1:d]

    F = [φ + y - u; Dφ' * y]

    # We sample a random t, find a normal vector on this point and then assemble
    # an u where we know the solution to.
    t_rand = randn(Float64, size(A, 1))
    null = nullspace(map(fi -> fi(t => t_rand), Dφ'))
    y_rand = null * randn(Float64, size(null, 2))
    u₀ = map(fi -> fi(t => t_rand), φ) + y_rand

    x₀ = complex.([t_rand; y_rand])

    F, u, u₀, x₀
end

@testset "Monodromy" begin

    @testset "monodromy_solve" begin
        F, p, p₀, x₀ = toric_ed([3 2 1 0; 0 1 2 3])

        # test that timeout works
        Random.seed!(51232)
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21, timeout=1e-12)
        @test length(result.solutions) < 21

        result = monodromy_solve(F, x₀, p₀, parameters=p,
                target_solutions_count=21,
                max_loops_no_progress=20)
        @test result.returncode == :success
        @test length(solutions(result)) == 21
        @test length(solutions(result, only_real = true)) >= 1
        @test result.statistics.ntrackedpaths ≥ 21
        @test result.statistics.nparametergenerations ≥ 1
        @test length(HC.UniquePoints(result.solutions).points) == 21
        @test isempty(sprint(show, result)) == false

        # test seed
        result2 = monodromy_solve(F, x₀, p₀, parameters=p,
                target_solutions_count=21,
                max_loops_no_progress=20,
                seed=result.seed)
        @test result2.statistics.ntrackedpaths == result.statistics.ntrackedpaths

        # test input of length > 1
        result = monodromy_solve(F, [x₀ for _ in 1:30], p₀, parameters=p)
        @test length(solutions(result)) == 21

        # test with false input
        result = monodromy_solve(F, [rand(6) for _ in 1:10], p₀, parameters=p)
        @test result.returncode == :invalid_startvalue
        result = monodromy_solve(F, [x₀, rand(6)], p₀, parameters=p)
        @test length(solutions(result)) == 21
        result = monodromy_solve(F, [x₀, rand(6)], p₀, parameters=p, check_startsolutions = false)
        @test result.returncode == :invalid_startvalue

        # different distance function
        result = monodromy_solve(F, x₀, p₀, parameters=p, distance = (x,y) -> 0.0)
        @test length(solutions(result)) == 1

        # Test stop heuristic using too high target_solutions_count
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=25)
        @test result.returncode == :heuristic_stop
        # Test stop heuristic with no target solutions count
        result = monodromy_solve(F, x₀, p₀, parameters=p)
        @test result.returncode == :heuristic_stop
        # Test stop heuristic with no target solutions count
        result = monodromy_solve(F, x₀, p₀, parameters=p, strategy=Triangle(useweights=true))
        @test result.returncode == :heuristic_stop

        # By group_actions=nothing we force that complex conjugation is not used.
        result2 = monodromy_solve(F, x₀, p₀, parameters=p,
                        target_solutions_count=21, complex_conjugation=false,
                        max_loops_no_progress=100)
        @test result2.returncode == :success

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
                done_callback=((_, _) -> true))
        @test length(result.solutions) == 2

        roots_of_unity(s) = begin
            t = cis(π*2/3)
            t² = t * t
            (vcat(t * s[1], t * s[2], s[3:end]),
             vcat(t² * s[1], t² * s[2], s[3:end]))
        end

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            max_loops_no_progress=100,
            equivalence_classes=false,
            group_action=roots_of_unity)
        @test length(result.solutions) == 21

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            max_loops_no_progress=100,
            complex_conjugation=false, # disable complex conjugation to test it as a group action.
            equivalence_classes=false,
            group_actions=(roots_of_unity, s -> (conj.(s),)))
        @test length(result.solutions) == 21
        @test length(solutions(result)) == 21
        @test length(real_solutions(result)) < 21
        test_treeviews(result)

        # group_actions as a vector
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            max_loops_no_progress=100,
            complex_conjugation=false,
            equivalence_classes=false,
            group_actions=[roots_of_unity, s -> (conj.(s),)])
        @test length(result.solutions) == 21

        # equivalence classes
        result = monodromy_solve(F, x₀, p₀, parameters=p,
            equivalence_classes=true,
            target_solutions_count=7,
            max_loops_no_progress=100,
            group_actions=roots_of_unity)
        @test length(result.solutions) == 7
        # Test that equivalence classes are on by default if we supply a group action
        result = monodromy_solve(F, x₀, p₀, parameters=p,
                            group_action=roots_of_unity,
                            max_loops_no_progress=20)
        @test length(result.solutions) == 7

        # Test affine tracking
        result = monodromy_solve(F, x₀, p₀, parameters=p, affine_tracking=true,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        max_loops_no_progress=200)
        @test length(result.solutions) == 7

        # AbstractSystem as input
        F_p = SPSystem(F; parameters=p)
        result = monodromy_solve(F_p, x₀, p₀, affine_tracking=true,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        max_loops_no_progress=200)
        @test length(result.solutions) == 7

        # test reuse strategies
        result = monodromy_solve(F_p, x₀, p₀,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        reuse_loops=:all,
                        max_loops_no_progress=200)
        @test length(result.solutions) == 7
        result = monodromy_solve(F_p, x₀, p₀,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        reuse_loops=:random,
                        max_loops_no_progress=200)
        @test length(result.solutions) == 7
        result = monodromy_solve(F_p, x₀, p₀,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        reuse_loops=:none,
                        max_loops_no_progress=200)
        @test length(result.solutions) == 7
    end

    @testset "Method of Moments" begin
        Random.seed!(130793)
        @polyvar a[1:3] x[1:3] s[1:3]

        f0 = a[1]+a[2]+a[3];
        f1 = a[1]*x[1]+a[2]*x[2]+a[3]*x[3];
        f2 = a[1]*(x[1]^2+s[1])+a[2]*(x[2]^2+s[2])+a[3]*(x[3]^2+s[3]);
        f3 = a[1]*(x[1]^3+3*s[1]*x[1])+a[2]*(x[2]^3+3*s[2]*x[2])+a[3]*(x[3]^3+3*s[3]*x[3]);
        f4 = a[1]*(x[1]^4+6*s[1]*x[1]^2+3*s[1]^2)+a[2]*(x[2]^4+6*s[2]*x[2]^2+3*s[2]^2)+
        	a[3]*(x[3]^4+6*s[3]*x[3]^2+3*s[3]^2);
        f5= a[1]*(x[1]^5+10*s[1]*x[1]^3+15*x[1]*s[1]^2)+a[2]*(x[2]^5+10*s[2]*x[2]^3+15*x[2]*s[2]^2)+a[3]*(x[3]^5+10*s[3]*x[3]^3+15*x[3]*s[3]^2);
        f6= a[1]*(x[1]^6+15*s[1]*x[1]^4+45*x[1]^2*s[1]^2+15*s[1]^3)+
        	a[2]*(x[2]^6+15*s[2]*x[2]^4+45*x[2]^2*s[2]^2+15*s[2]^3)+
        	a[3]*(x[3]^6+15*s[3]*x[3]^4+45*x[3]^2*s[3]^2+15*s[3]^3);
        f7= a[1]*(x[1]^7+21*s[1]*x[1]^5+105*x[1]^3*s[1]^2+105*x[1]*s[1]^3)+
        	a[2]*(x[2]^7+21*s[2]*x[2]^5+105*x[2]^3*s[2]^2+105*x[2]*s[2]^3)+
        	a[3]*(x[3]^7+21*s[3]*x[3]^5+105*x[3]^3*s[3]^2+105*x[3]*s[3]^3);
        f8= a[1]*(x[1]^8+28*s[1]*x[1]^6+210*x[1]^4*s[1]^2+420*x[1]^2*s[1]^3+105*s[1]^4)+
        	a[2]*(x[2]^8+28*s[2]*x[2]^6+210*x[2]^4*s[2]^2+420*x[2]^2*s[2]^3+105*s[2]^4)+
        	a[3]*(x[3]^8+28*s[3]*x[3]^6+210*x[3]^4*s[3]^2+420*x[3]^2*s[3]^3+105*s[3]^4)

        f = [f0, f1, f2, f3, f4, f5, f6, f7, f8]

        S₃ = SymmetricGroup(3)
        relabeling(v) = map(p -> (v[1:3][p]..., v[4:6][p]..., v[7:9][p]...), S₃)

        @polyvar p[1:9]

        function sample_moments(_=nothing)
        	a₀ = rand(3)
        	a₀ ./= sum(a₀)
        	x₀ = randn(3)
        	s₀ = rand(3) .* 2

        	y₀ = [a₀; x₀; s₀]
        	p₀ = [fᵢ(a => a₀, x => x₀, s => s₀) for fᵢ in f]
        	y₀, p₀
        end

        y₀, p₀ = sample_moments()

        R = monodromy_solve(f - p, y₀, p₀;
        			parameters=p, group_action=relabeling,
        			parameter_sampler=last ∘ sample_moments,
                    show_progress=false,
                    max_loops_no_progress=5)
        @test length(solutions(R)) ≤ 225

        R = monodromy_solve(f - p, y₀, p₀;
        			parameters=p, group_action=relabeling,
        			target_solutions_count=225,
                    affine_tracking=true,
        			parameter_sampler=last ∘ sample_moments,
                    show_progress=false,
                    max_loops_no_progress=20)
        @test length(solutions(R)) == 225
    end

    @testset "Trace Test" begin

    	@testset "Simple" begin
    		@polyvar x y a b c;
    		f = x^2+y^2-1;
    		l = a*x+b*y+c;
    		res = monodromy_solve([f,l], [-0.6-0.8im, -1.2+0.4im], [1,2,3]; parameters=[a,b,c])
    		@test verify_solution_completeness([f,l], res; parameters=[a,b,c])
    	end

    	@testset "ED Example" begin
    		# define f
    		@polyvar x y
    		f = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
    		# define new variables u₁, u₂ and λ₁
    		@polyvar u[1:2] λ[1:1]
    		# define the jacobian of F
    		J = differentiate([f], [x,y])
    		# J' defines the transpose of J
    		C = [[x,y] - u - J'*λ; f]

    		# create a random linear space
    		l = sum(randn(ComplexF64,3) .* [x,y,1])
    		# compute intersection points of the V(f) and V(l) and take the first solution
    		s₀ = first(solutions(solve([f, l])))
    		# The normal space is now generated by
    		J_s₀ = [p([x,y]=>s₀) for p in J]
    		# creat normal vector
    		λ₀ = cis(2π*rand()) / norm(J_s₀)
    		# x₀ - u₀ = J_s₀' * λ₀ => u = x -  J_s₀' * λ₀
    		u₀ = s₀ - transpose(J_s₀) * λ₀ |> vec

    		# This is not as realiable as I would like it to be, so just make a couple runs...
    		is_complete = false
    		for i in 1:5
    			res = monodromy_solve(C, [s₀;λ₀], u₀; parameters=u)
    			 if verify_solution_completeness(C, res; parameters=u)
    			 	is_complete = true
    				break
    			end
    		 end
    		 @test is_complete
    	 end
    end
end


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export _wsave

using PyPlot: Figure
import DrWatson: _wsave

_wsave(s, fig::Figure; dpi::Int=200) = fig.savefig(s, bbox_inches="tight", dpi=dpi)