<jupyter_start><jupyter_text>Load data<jupyter_code># @time df_spatial, gene_names = Baysor.load_df("../run_results/spacejam2/allen_sm_fish/no_dapi/segmentation.csv");
# df_spatial[!, :x] = round.(Int, 10 .* (df_spatial.x .- minimum(df_spatial.x)));
# df_spatial[!, :y] = round.(Int, 10 .* (df_spatial.y .- minimum(df_spatial.y)));
# length(gene_names)
@time df_spatial, gene_names = Baysor.load_df("../run_results/merfish_moffit/segmentation.csv");
length(gene_names)<jupyter_output><empty_output><jupyter_text>Molecule clustering  Baysor<jupyter_code>bench_df = @where(df_spatial, :x .< -3300, :y .< -3300) |> deepcopy;
gn_bench = gene_names;

# confidence_nn_id = Baysor.default_param_value(:confidence_nn_id, 10);
confidence_nn_id = Baysor.default_param_value(:confidence_nn_id, 50);
@show confidence_nn_id
size(bench_df, 1)
bench_clust = BenchmarkGroup();
B.append_confidence!(bench_df, nn_id=confidence_nn_id);
bench_clust["confidence"] = @benchmarkable B.append_confidence!($bench_df, nn_id=$confidence_nn_id);
adjacent_points, adjacent_weights = B.build_molecule_graph(bench_df, filter=false);
bench_clust["mol_graph"] = @benchmarkable B.build_molecule_graph($bench_df, filter=false);
for cl in [2, 4, 6, 8, 10]
    bench_clust["clust_$cl"] = @benchmarkable B.cluster_molecules_on_mrf($bench_df.gene, $adjacent_points, $adjacent_weights, $bench_df.confidence; 
        n_clusters=$cl, max_iters=5000, n_iters_without_update=100, verbose=false);
end
bench_clust_res = run(bench_clust)
bench_res_df = vcat([DataFrame("Key" => k, "Mean time, sec" => mean(v.times) ./ 1e9, "Std time, sec" => std(v.times) ./ 1e9, 
            "Num. samples" => length(v.times)) for (k,v) in bench_clust_res]...)<jupyter_output><empty_output><jupyter_text>Leiden<jupyter_code>using RCall
nm_bench = B.neighborhood_count_matrix(bench_df, 50, normalize=false);
size(nm_bench)
R"""
library(pagoda2)
library(conos)
library(microbenchmark)

cm <- as($nm_bench, "dgCMatrix")
rownames(cm) <- $gn_bench
colnames(cm) <- paste0("c", 1:ncol(cm))

getClusters <- function(cm, verbose=FALSE) {
    p2 <- Pagoda2$new(cm, trim=5, n.cores=1, verbose=FALSE, log.scale=verbose)
    p2$calculatePcaReduction(nPcs=50, odgenes=rownames(cm), maxit=1000, verbose=verbose, var.scale=FALSE)
    p2$makeKnnGraph(k=30, type="PCA", center=T, distance="cosine", weight.type="none", verbose=verbose)
    p2$getKnnClusters(method=conos::leiden.community, type="PCA", name="leiden", resolution=1.0)
    
    return(p2$clusters$PCA$leiden)
}

b <- microbenchmark(
    "clustering" = {getClusters(cm)},
    times=5,
    control=list(warmup=1)
)
"""<jupyter_output><empty_output><jupyter_text>Aggregate<jupyter_code>leiden_times = rcopy(R"b").time;
bench_res_df
df1 = hcat(DataFrame("Method" => "MRF", "Num. clusters" => 2:2:10), bench_res_df[[3, 1, 5, 4, 2],2:end]);
df2 = vcat(df1, DataFrame("Method" => "Leiden", "Num. clusters" => "Any", "Mean time, sec" => mean(leiden_times) / 1e9, 
        "Std time, sec" => std(leiden_times) / 1e9, "Num. samples" => 5));

df2[:, 3:4] .= round.(df2[:, 3:4], digits=2);
df2
CSV.write("plots/clustering_profiling.csv", df2)<jupyter_output><empty_output><jupyter_text>Color embedding<jupyter_code>@time neighb_cm = B.neighborhood_count_matrix(df_spatial, 40);
@time color_transformation = B.gene_composition_transformation(neighb_cm, df_spatial.confidence; sample_size=20000, spread=2.0, min_dist=0.1);
@time color_emb = B.transform(color_transformation, neighb_cm);
bench_emb = BenchmarkGroup();
bench_emb["neighborhood_count_matrix_40"] = @benchmarkable B.neighborhood_count_matrix($df_spatial, 40)
bench_emb["gene_composition_transformation_20k"] = @benchmarkable B.gene_composition_transformation(neighb_cm, df_spatial.confidence; 
    sample_size=20000, spread=2.0, min_dist=0.1)
bench_emb["transform"] = @benchmarkable B.transform(color_transformation, neighb_cm)
bench_emb_res = run(bench_emb)
bench_df = vcat([DataFrame("Key" => k, "Mean time, sec" => mean(v.times) ./ 1e9, "Std time, sec" => std(v.times) ./ 1e9, 
            "Num. samples" => length(v.times)) for (k,v) in bench_emb_res]...)<jupyter_output><empty_output><jupyter_text>Segmentation<jupyter_code>bench_segmentation = BenchmarkGroup();
@time df_spatial, gene_names = B.load_df("../run_results/iss_hippo/ca1_no_prior/segmentation.csv");
df_spatial[!, :cell_dapi] = df_spatial.parent_id;
dapi_arr = Float16.(Images.load("/home/vpetukhov/data/spatal/iss/hippocampus/CA1/Viktor/CA1DapiBoundaries_4-3_right.tif"));
iss = Dict(:df => df_spatial, :gene_names => gene_names, :name => "ISS", :dapi_arr => dapi_arr);
B.append_confidence!(df_spatial, (args["prior_segmentation"]===nothing ? nothing : df_spatial.prior_segmentation), nn_id=confidence_nn_id, prior_confidence=args["prior-segmentation-confidence"])
adjacent_points, adjacent_weights = build_molecule_graph(df_spatial, filter=false)[1:2];

mol_clusts = cluster_molecules_on_mrf(df_spatial.gene, adjacent_points, adjacent_weights, df_spatial.confidence;
            n_clusters=args["n-clusters"], weights_pre_adjusted=true)

df_spatial[!, :cluster] = mol_clusts.assignment;

bm_data_arr = initial_distribution_arr(df_spatial; n_frames=args["n-frames"], scale=args["scale"], scale_std=args["scale-std"],
        n_cells_init=args["num-cells-init"], prior_seg_confidence=args["prior-segmentation-confidence"],
        min_molecules_per_cell=args["min-molecules-per-cell"], confidence_nn_id=0);

bm_data = run_bmm_parallel!(bm_data_arr, args["iters"], new_component_frac=args["new-component-fraction"], new_component_weight=args["new-component-weight"],
                            min_molecules_per_cell=args["min-molecules-per-cell"], assignment_history_depth=history_depth);
cur_df = deepcopy(iss[:df]);
bm_data = B.initial_distribution_arr(cur_df; n_frames=1, scale=14, scale_std="25%", min_molecules_per_cell=3)[1];
@time B.bmm!(bm_data, n_iters=350, new_component_frac=0.3, min_molecules_per_cell=3, assignment_history_depth=30, log_step=100);
cur_df[!, :cell] = B.estimate_assignment_by_history(bm_data)[1];

B.plot_comparison_for_cell(cur_df, B.val_range(cur_df.x), B.val_range(cur_df.y), nothing, iss[:dapi_arr];
    ms=2.0, bandwidth=5.0, size_mult=0.25, plot_raw_dapi=false)<jupyter_output><empty_output><jupyter_text>Full run  Run<jupyter_code>using ProgressMeter
dataset_paths = "/home/vpetukhov/spatial/Benchmarking/run_results/" .* 
    ["iss_hippo/ca1_no_prior", "merfish_moffit", "osm_fish", "star_map/vis_1020_cl0", "spacejam2/allen_sm_fish/no_dapi"];
param_dumps = dataset_paths .* "/segmentation_params.dump";

dataset_names = ["iss", "merfish", "osm_fish", "starmap_1020", "allen_smfish"];

param_strings = [open(p) do f readlines(f)[1][16:end-1] end for p in param_dumps];
baysor_path = "/home/vpetukhov/local/bin/baysor";
for i in 2:length(param_strings)
# for i in 2:2
    dataset = dataset_names[i]
    params = split(param_strings[i], ' ')

    out_path = expanduser("/home/vpetukhov/spatial/Benchmarking/run_results/profiling/$dataset/")
    mkpath(out_path)
    cmd = `/usr/bin/time -f '%e %U %P %M %t %K' -o ./profiling_output/$dataset.prof -a $baysor_path run --debug -o $out_path $params`;
#     cmd = `/usr/bin/time -f '%e %U %P %M %t %K' -o ./profiling_output/$dataset.prof -a $baysor_path run --debug --n-clusters=0 -o $out_path $params`;
    @show cmd
    
    println(dataset)
    @showprogress for ri in 1:5
        run(pipeline(cmd, stdout="./profiling_output/$dataset.log", stderr="./profiling_output/$dataset.err", append=true))
        run(pipeline(`echo -e \\n\\n\\n ----- RUN $ri ----- \\n\\n\\n`, stdout="./profiling_output/$dataset.log", append=true))
    end
end<jupyter_output><empty_output><jupyter_text>Summarize<jupyter_code>using DataFrames
using Statistics

printed_names = ["ISS", "MERFISH", "osmFISH", "STARmap 1020", "Allen smFISH"];

seg_results = dataset_paths .* "/segmentation.csv";
dataset_parameters = hcat([[size(df, 1), length(unique(df.gene))] for df in DataFrame!.(CSV.File.(seg_results))]...);
bench_vals = [hcat(split.(readlines("./profiling_output/$ds.prof"), ' ')...) for ds in dataset_names];
mem_vals = hcat([parse.(Float64, x[4,:]) / 1e6 for x in bench_vals]...);
cpu_vals = hcat([parse.(Float64, x[1,:]) / 60 for x in bench_vals]...);

bench_mat = round.(vcat(mean(cpu_vals, dims=1), std(cpu_vals, dims=1), mean(mem_vals, dims=1), std(mem_vals, dims=1))', digits=2);
bench_strs = [["$(r[i[1]]) ± $(r[i[2]])" for r in eachrow(bench_mat)] for i in ((1, 2), (3, 4))];
bench_df = DataFrame("Dataset" => printed_names, "Num. molecules" => dataset_parameters[1,:], "Num. genes" => dataset_parameters[2,:],
    "CPU time, min" => bench_strs[1], "Max RSS, GB" => bench_strs[2], "Num. samples" => 5)
CSV.write("./plots/segmentation_profiling.csv", bench_df)<jupyter_output><empty_output><jupyter_text>Parameter table<jupyter_code>import Pkg: TOML
using DataFrames
import CSV
data_paths = ["Allen smFISH" => "allen_smfish", "ISS" => "iss_hippo", "osmFISH" => "osmfish", "STARmap 1020" => "starmap_vis1020", "MERFISH Hypothalamus" => "merfish_moffit", "MERFISH Gut" => "merfish_membrane"];
prior_subfolders = ["No" => "baysor", "Paper" => "baysor_prior", "DAPI" => "baysor_dapi_prior", "Membrane" => "baysor_membrane_prior"];

p_keys = ["gene-composition-neigborhood", "scale", "prior-segmentation-confidence", "min-molecules-per-gene", "min-molecules-per-cell", "n-clusters", 
    "iters", "force-2d", "x-column", "y-column", "z-column", "gene-column", "prior_segmentation", "nuclei-genes", "cyto-genes"];
path_df = DataFrame([Dict(:Dataset => d, :Prior => pr, :Path => datadir("exp_pro", md, sd, "segmentation_params.dump")) for (d, md) in data_paths for (pr, sd) in prior_subfolders]);
path_df = path_df[isfile.(path_df.Path),:];

param_dicts = [OrderedDict(k => get(d, k, "NA") for k in p_keys) for d in TOML.parsefile.(path_df.Path)];
param_df = hcat(path_df[:,[:Dataset, :Prior]], vcat(DataFrame.(param_dicts)...))
CSV.write(plotsdir("parameters.csv"), param_df)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1734857601283594, 'fraction_numerical': 0.021460088247091856, 'mean_word_length': 6.532477341389728, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 77, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '369734', 'license': 'MIT', 'n_tokens_mistral': 4126, 'n_tokens_neox': 3836, 'n_words': 614, 'path': 'notebooks/benchmarking/performance_profiling.ipynb', 'repo_name': 'hiraksarkar/BaysorAnalysis', 'size': '9972'}