From bdeab7acf53d41f9b0eb2b5748e78db7e68a96f8 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 15 Apr 2021 18:47:48 +0100 Subject: [PATCH 01/18] minor changes for tnblast testing --- R/neuriteblast.r | 4 +++- workbook.R | 57 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 1 deletion(-) create mode 100644 workbook.R diff --git a/R/neuriteblast.r b/R/neuriteblast.r index ab83657..637cf8a 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -473,7 +473,9 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL # sqrt seems reasonable since if alpha1=alpha2=0.5 then # the scalefac will be 0.5 # zapsmall to make sure there are no tiny negative numbers etc. - scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) + #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) + maxlen = max(max(alphas1), max(alphas2)) + scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])^2)/maxlen) dps=dps*scalefac } } diff --git a/workbook.R b/workbook.R new file mode 100644 index 0000000..8765cb6 --- /dev/null +++ b/workbook.R @@ -0,0 +1,57 @@ +library(nat) +library(catmaid) + +conn <- catmaid_login(server="https://neuropil.janelia.org/tracing/fafb/v14/", + authname="fly", authpassword="superfly", + token = "c3990b12d59fb66d59107fb16e0540c4f3e91aa5") +catmaid_login(conn) + +neurons <- read_catmaid_selection("../tnblast/data/DA1s.json", readNeurons = TRUE) +neurons <- resample(prune_twigs(neurons/1e3, twig_length = 2), stepsize=1) + +dps <- nat::dotprops(neurons) + +nrn1 <- neurons[[1]] +nrn2 <- neurons[[2]] +nrn3 <- neurons[[3]] +nrn4 <- neurons[[4]] + +#somaid(nns) # resampled neurons don't have soma + +dps1 <- dps[[1]] +dps2 <- dps[[2]] +dps3 <- dps[[3]] +dps4 <- dps[[4]] + +gw1 <- as.ngraph(nrn1, weights=TRUE) +gw2 <- as.ngraph(nrn2, weights=TRUE) +gw3 <- as.ngraph(nrn3, weights=TRUE) +gw4 <- as.ngraph(nrn4, weights=TRUE) + +igraph::topo_sort(as.ngraph(nrn4)) + +max(igraph::diameter(gw1), igraph::diameter(gw2)) + +dst1 <- igraph::distances(gw1, v=as.character(rootpoints(nrn1))) + +get_dist_to_soma <- function(nrn) { + gw <- as.ngraph(nrn, weights=TRUE) + dst <- igraph::distances(gw, v=as.character(rootpoints(nrn))) + as.numeric(dst) +} + +dps1$alpha <- get_dist_to_soma(nrn1) +dps2$alpha <- get_dist_to_soma(nrn2) +dps3$alpha <- get_dist_to_soma(nrn3) +dps4$alpha <- get_dist_to_soma(nrn4) + +dps_list <- neuronlist(dps1, dps2, dps3, dps4) + +devtools::load_all() + +dps_aba1 <- nblast_allbyall(dps_list, normalisation = "raw") +dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "raw") +par(mfrow=c(1,2)) +heatmap(dps_aba1) +heatmap(dps_aba2) + From f50c18fc516b6af640d13ef2b02dbc21bc10db25 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 15 Apr 2021 20:42:09 +0100 Subject: [PATCH 02/18] simplified workbook --- workbook.R | 50 ++++++++++++++++++++------------------------------ 1 file changed, 20 insertions(+), 30 deletions(-) diff --git a/workbook.R b/workbook.R index 8765cb6..722fd9b 100644 --- a/workbook.R +++ b/workbook.R @@ -11,28 +11,9 @@ neurons <- resample(prune_twigs(neurons/1e3, twig_length = 2), stepsize=1) dps <- nat::dotprops(neurons) -nrn1 <- neurons[[1]] -nrn2 <- neurons[[2]] -nrn3 <- neurons[[3]] -nrn4 <- neurons[[4]] - -#somaid(nns) # resampled neurons don't have soma - -dps1 <- dps[[1]] -dps2 <- dps[[2]] -dps3 <- dps[[3]] -dps4 <- dps[[4]] - -gw1 <- as.ngraph(nrn1, weights=TRUE) -gw2 <- as.ngraph(nrn2, weights=TRUE) -gw3 <- as.ngraph(nrn3, weights=TRUE) -gw4 <- as.ngraph(nrn4, weights=TRUE) - -igraph::topo_sort(as.ngraph(nrn4)) - -max(igraph::diameter(gw1), igraph::diameter(gw2)) - -dst1 <- igraph::distances(gw1, v=as.character(rootpoints(nrn1))) +#somaid(neurons) # resampled neurons don't have soma +#igraph::topo_sort(as.ngraph(nrn4)) +#max(igraph::diameter(gw1), igraph::diameter(gw2)) get_dist_to_soma <- function(nrn) { gw <- as.ngraph(nrn, weights=TRUE) @@ -40,18 +21,27 @@ get_dist_to_soma <- function(nrn) { as.numeric(dst) } -dps1$alpha <- get_dist_to_soma(nrn1) -dps2$alpha <- get_dist_to_soma(nrn2) -dps3$alpha <- get_dist_to_soma(nrn3) -dps4$alpha <- get_dist_to_soma(nrn4) +make_topo_dotprops <- function(nrn) { + tdps <- nat::dotprops(nrn) + tdps$alpha <- get_dist_to_soma(nrn) + tdps +} -dps_list <- neuronlist(dps1, dps2, dps3, dps4) +dps_list <- nlapply(neurons, make_topo_dotprops) devtools::load_all() dps_aba1 <- nblast_allbyall(dps_list, normalisation = "raw") dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "raw") -par(mfrow=c(1,2)) -heatmap(dps_aba1) -heatmap(dps_aba2) +par(mfrow = c(1,2)) +image(dps_aba1) +title("NBLAST") +image(dps_aba2) +title("TNBLAST") + +par(mfrow = c(1,2)) +image(dps_aba1 / diag(dps_aba1)) +title("NBLAST") +image(dps_aba2 / diag(dps_aba2)) +title("TNBLAST") From 6505da6cd661fc1d1590ae2d456dee0c0dd8a285 Mon Sep 17 00:00:00 2001 From: dokato Date: Thu, 15 Apr 2021 20:43:23 +0100 Subject: [PATCH 03/18] Update workbook.R --- workbook.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/workbook.R b/workbook.R index 722fd9b..3ec8ddc 100644 --- a/workbook.R +++ b/workbook.R @@ -1,9 +1,6 @@ library(nat) library(catmaid) -conn <- catmaid_login(server="https://neuropil.janelia.org/tracing/fafb/v14/", - authname="fly", authpassword="superfly", - token = "c3990b12d59fb66d59107fb16e0540c4f3e91aa5") catmaid_login(conn) neurons <- read_catmaid_selection("../tnblast/data/DA1s.json", readNeurons = TRUE) From 8d0022a7b6b869250b63226db8b52f961508145b Mon Sep 17 00:00:00 2001 From: dokato Date: Thu, 15 Apr 2021 20:43:39 +0100 Subject: [PATCH 04/18] Update workbook.R --- workbook.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workbook.R b/workbook.R index 3ec8ddc..9b61403 100644 --- a/workbook.R +++ b/workbook.R @@ -1,7 +1,7 @@ library(nat) library(catmaid) -catmaid_login(conn) +catmaid_login() neurons <- read_catmaid_selection("../tnblast/data/DA1s.json", readNeurons = TRUE) neurons <- resample(prune_twigs(neurons/1e3, twig_length = 2), stepsize=1) From 1b15421cdf4313fca10bcb74d122eab15c5d7d50 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Fri, 16 Apr 2021 21:21:03 +0100 Subject: [PATCH 05/18] fix scalefac cost fun --- R/neuriteblast.r | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index 637cf8a..eeeb410 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -475,7 +475,7 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL # zapsmall to make sure there are no tiny negative numbers etc. #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) maxlen = max(max(alphas1), max(alphas2)) - scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])^2)/maxlen) + scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^2) dps=dps*scalefac } } From 3901bd3c061564f0d7454fa22dbca66a08fc12fe Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Fri, 16 Apr 2021 21:27:57 +0100 Subject: [PATCH 06/18] fix scalefac function --- R/neuriteblast.r | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index eeeb410..f546867 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -475,7 +475,7 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL # zapsmall to make sure there are no tiny negative numbers etc. #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) maxlen = max(max(alphas1), max(alphas2)) - scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^2) + scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^0.5) dps=dps*scalefac } } From fbd4504b08ce29fcdc22fd6642719e81ae5a2d6b Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Fri, 16 Apr 2021 21:58:35 +0100 Subject: [PATCH 07/18] adding Strahler Ord score --- R/neuriteblast.r | 17 +++++++++++++++-- workbook.R | 23 ++++++++++++++++++----- 2 files changed, 33 insertions(+), 7 deletions(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index f546867..d15ebe5 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -474,8 +474,21 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL # the scalefac will be 0.5 # zapsmall to make sure there are no tiny negative numbers etc. #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) - maxlen = max(max(alphas1), max(alphas2)) - scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^0.5) + if (is.numeric(alphas1) && is.numeric(alphas2)) { + maxlen = max(max(alphas1), max(alphas2)) + scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^0.5) + } else if (is.list(alphas1) && is.list(alphas2)) { + # distance from soma + maxlen = max(max(alphas1$distance), max(alphas2$distance)) + scalefac1 = zapsmall(1 - (abs(alphas1$distance[idxArray[,1]] - alphas2$distance[idxArray[,2]])/maxlen)^0.5) + # strenghten by Strahler Order + maxlen = max(max(alphas1$so), max(alphas2$so)) + scalefac2 = zapsmall(1 - (abs(alphas1$so[idxArray[,1]] - alphas2$so[idxArray[,2]])/maxlen)^0.5) + scalefac = scalefac1 * scalefac2 + } else{ + warning("Unknown alpha format") + scalefac = rep(1, length(dps)) + } dps=dps*scalefac } } diff --git a/workbook.R b/workbook.R index 9b61403..018adab 100644 --- a/workbook.R +++ b/workbook.R @@ -32,13 +32,26 @@ dps_aba1 <- nblast_allbyall(dps_list, normalisation = "raw") dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "raw") par(mfrow = c(1,2)) -image(dps_aba1) +image(dps_aba1 / diag(dps_aba1)) title("NBLAST") -image(dps_aba2) +image(dps_aba2 / diag(dps_aba2)) title("TNBLAST") +make_sotopo_dotprops <- function(nrn) { + tdps <- nat::dotprops(nrn) + tdps$alpha <- list() + tdps$alpha$distance <- get_dist_to_soma(nrn) + tdps$alpha$so <- strahler_order(nrn)$points + tdps +} + +dps_list2 <- nlapply(neurons, make_sotopo_dotprops) + +dps_aba21 <- nblast_allbyall(dps_list2, normalisation = "raw") +dps_aba22 <- nblast_allbyall(dps_list2, UseAlpha = T, normalisation = "raw") + par(mfrow = c(1,2)) -image(dps_aba1 / diag(dps_aba1)) +image(dps_aba21 / diag(dps_aba21)) title("NBLAST") -image(dps_aba2 / diag(dps_aba2)) -title("TNBLAST") +image(dps_aba22 / diag(dps_aba22)) +title("TNBLAST + SO") From 889bcbd0d4ba910a07aa5b1867bd3a1638b8f85b Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 29 Apr 2021 13:03:46 +0100 Subject: [PATCH 08/18] workbooks for a few datasets added --- R/neuriteblast.r | 4 +-- workbook.R | 59 +++++++++++++++++++----------------- workbook_fw.R | 51 +++++++++++++++++++++++++++++++ workbook_hemi.R | 13 ++++++++ workbook_swc.R | 79 ++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 177 insertions(+), 29 deletions(-) create mode 100644 workbook_fw.R create mode 100644 workbook_hemi.R create mode 100644 workbook_swc.R diff --git a/R/neuriteblast.r b/R/neuriteblast.r index d15ebe5..525722e 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -476,11 +476,11 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) if (is.numeric(alphas1) && is.numeric(alphas2)) { maxlen = max(max(alphas1), max(alphas2)) - scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^0.5) + scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^2) } else if (is.list(alphas1) && is.list(alphas2)) { # distance from soma maxlen = max(max(alphas1$distance), max(alphas2$distance)) - scalefac1 = zapsmall(1 - (abs(alphas1$distance[idxArray[,1]] - alphas2$distance[idxArray[,2]])/maxlen)^0.5) + scalefac1 = zapsmall(1 - (abs(alphas1$distance[idxArray[,1]] - alphas2$distance[idxArray[,2]])/maxlen)^2) # strenghten by Strahler Order maxlen = max(max(alphas1$so), max(alphas2$so)) scalefac2 = zapsmall(1 - (abs(alphas1$so[idxArray[,1]] - alphas2$so[idxArray[,2]])/maxlen)^0.5) diff --git a/workbook.R b/workbook.R index 018adab..014c46f 100644 --- a/workbook.R +++ b/workbook.R @@ -1,57 +1,62 @@ library(nat) library(catmaid) +library(dendroextras) + +devtools::load_all() catmaid_login() neurons <- read_catmaid_selection("../tnblast/data/DA1s.json", readNeurons = TRUE) neurons <- resample(prune_twigs(neurons/1e3, twig_length = 2), stepsize=1) -dps <- nat::dotprops(neurons) +plot3d(neurons, soma = TRUE) #somaid(neurons) # resampled neurons don't have soma #igraph::topo_sort(as.ngraph(nrn4)) #max(igraph::diameter(gw1), igraph::diameter(gw2)) -get_dist_to_soma <- function(nrn) { - gw <- as.ngraph(nrn, weights=TRUE) - dst <- igraph::distances(gw, v=as.character(rootpoints(nrn))) - as.numeric(dst) -} - -make_topo_dotprops <- function(nrn) { - tdps <- nat::dotprops(nrn) - tdps$alpha <- get_dist_to_soma(nrn) - tdps -} - dps_list <- nlapply(neurons, make_topo_dotprops) -devtools::load_all() - dps_aba1 <- nblast_allbyall(dps_list, normalisation = "raw") dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "raw") par(mfrow = c(1,2)) -image(dps_aba1 / diag(dps_aba1)) +image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n') title("NBLAST") -image(dps_aba2 / diag(dps_aba2)) +image(dps_aba2 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n') title("TNBLAST") -make_sotopo_dotprops <- function(nrn) { - tdps <- nat::dotprops(nrn) - tdps$alpha <- list() - tdps$alpha$distance <- get_dist_to_soma(nrn) - tdps$alpha$so <- strahler_order(nrn)$points - tdps -} - dps_list2 <- nlapply(neurons, make_sotopo_dotprops) dps_aba21 <- nblast_allbyall(dps_list2, normalisation = "raw") dps_aba22 <- nblast_allbyall(dps_list2, UseAlpha = T, normalisation = "raw") par(mfrow = c(1,2)) -image(dps_aba21 / diag(dps_aba21)) +image(dps_aba21 / diag(dps_aba21), zlim= c(0,1), yaxt='n', xaxt='n') title("NBLAST") -image(dps_aba22 / diag(dps_aba22)) +image(dps_aba22 / diag(dps_aba22), zlim= c(0,1), yaxt='n', xaxt='n') title("TNBLAST + SO") + +# ---------------------- artificial neurons +nrn3x <- neurons[[3]] +nrn3x$StartPoint <- endpoints(neurons[[3]])[length(endpoints(neurons[[3]]))] + +nrn3x <- neurons[[3]] +nrn3x$StartPoint <- endpoints(neurons[[3]])[106] +nrn4x <- neurons[[4]] +nrn4x$StartPoint <- endpoints(neurons[[4]])[length(endpoints(neurons[[4]]))] + +newnrns <- neuronlist(neurons[[1]], neurons[[2]], nrn3x, nrn4x) + +plot(newnrns,soma = T) + +dps_listT1 <- nlapply(newnrns, make_topo_dotprops) + +dps_aba1 <- nblast_allbyall(dps_listT1, normalisation = "raw") +dps_aba2 <- nblast_allbyall(dps_listT1, UseAlpha = T, normalisation = "raw") + +par(mfrow = c(1,2)) +image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n') +title("NBLAST") +image(dps_aba2 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n') +title("TNBLAST") diff --git a/workbook_fw.R b/workbook_fw.R new file mode 100644 index 0000000..dc57151 --- /dev/null +++ b/workbook_fw.R @@ -0,0 +1,51 @@ +library(nat) +library(fafbseg) +library(dendroextras) + +devtools::load_all() + +OUT_DIR <- "meshes" + +fw_pair <- flywire_fetch("https://globalv1.flywire-daf.com/nglstate/5478428614590464", return="text") +fw_scene <- flywire_fetch("https://globalv1.flywire-daf.com/nglstate/6012751403024384", return="text") + +fw_nids1 <- ngl_segments(ngl_decode_scene(fw_pair), as_character = TRUE) +fw_nids2 <- ngl_segments(ngl_decode_scene(fw_scene), as_character = TRUE) +#save_cloudvolume_meshes(fw_nids, savedir = OUT_DIR, format = 'obj') + +nids <- c(fw_nids1, fw_nids2) + +fw_nrns <- skeletor(nids, method = "wavefront") + +saveRDS(fw_nrns, "fw_neurons.rds") +fw_nrns <- readRDS("fw_neurons.rds") + +fw_nrns_so <- nlapply(fw_nrns, function(x) reroot_hairball(x, k.soma.search = 50, radius.soma.search = 2500)) + +fw_nrns_cl <- nlapply(fw_nrns_so/1e3, stitch_neurons_mst, OmitFailures=T) %>% + prune_twigs(twig_length=2, OmitFailures=T) + +fw_nrns_sp <- nlapply(fw_nrns_cl, function(x) simplify_neuron(x, n=2)) + +plot3d(fw_nrns_sp, soma = T) + +dps_list <- nlapply(fw_nrns_sp, make_topo_dotprops) + +dps_aba1 <- nblast_allbyall(dps_list, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "normalised") + +par(mfrow = c(1,2)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n') +title("NBLAST") +image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n') +title("TNBLAST") + +hckcs1 <- nhclust(scoremat=dps_aba1) +dev.off() +dkcs1 <- colour_clusters(hckcs1, k=3) +plot(dkcs1) + +hckcs2 <- nhclust(scoremat=dps_aba2) +dkcs2 <- colour_clusters(hckcs2, k=3) +plot(dkcs2) + diff --git a/workbook_hemi.R b/workbook_hemi.R new file mode 100644 index 0000000..23c4068 --- /dev/null +++ b/workbook_hemi.R @@ -0,0 +1,13 @@ +library(nat) +library(fafbseg) +library(hemibrainr) +library(dendroextras) + +SKEL_DIR <- "/Users/dominik/projects/nblast_scoring_data/datasets/hemibrain_mismatch" +hemi_nrns <- read.neurons(SKEL_DIR) +plot3d(hemi_nrns, soma=T) + +sc_pairs <- read.csv(paste(SKEL_DIR, "../hemibrain_mismatch_pairs.csv", sep="/")) + +names(hemi_nrns) <- nlapply(hemi_nrns, function(nn) strsplit(nn$NeuronName,"_")[[1]][[3]]) + diff --git a/workbook_swc.R b/workbook_swc.R new file mode 100644 index 0000000..8e7802e --- /dev/null +++ b/workbook_swc.R @@ -0,0 +1,79 @@ +library(nat) +library(fafbseg) +library(dendroextras) +library(elmr) +library('fossil') + +devtools::load_all() + +SKEL_DIR <- "/Users/dominik/projects/nblast_scoring_data/datasets/FAFB_RHS_uPNs" +nrns <- read.neurons(SKEL_DIR) +plot3d(nrns, soma=T) + +lin_type <- unlist(nlapply(nrns, function(n) strsplit(n$NeuronName, "PN_")[[1]][[1]])) + +nrns <- nrns/1e3 +nrns_s <- nlapply(nrns, function(x) elmr::unspike(x, threshold=5)) +nrns_s <- nlapply(nrns_s, function(x) smooth_neuron(x, sigma=1)) +nrns_cl <- prune_twigs(nrns_s, twig_length=5, OmitFailures=T) + +plot3d(nrns_cl, soma=T) + +dps_list <- nlapply(nrns_cl, make_topo_dotprops) + +plot3d(dps_list) + +dps_aba1 <- nblast_allbyall(dps_list, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(dps_list, UseAlpha = T, normalisation = "normalised") + +dps_list_so <- nlapply(nrns_cl, make_sotopo_dotprops) + +dps_aba3 <- nblast_allbyall(dps_list_so, UseAlpha = T, normalisation = "normalised") + +par(mfrow = c(1,2)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST") + +par(mfrow = c(1,2)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba3, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST+SO") +dev.off() + +hcpn1 <- nhclust(scoremat=dps_aba1) +hcpn2 <- nhclust(scoremat=dps_aba2) +hcpn3 <- nhclust(scoremat=dps_aba3) + +k = length(unique(lin_type)) + +dkcs1 = colour_clusters(hcpn1, k=k) +labels(dkcs1) <- lin_type +plot(dkcs1) + +plot3d(hcpn1, k=k, db=dps_list) + +dkcs2 = colour_clusters(hcpn2, k=k) +labels(dkcs2) <- lin_type +plot(dkcs2) + +dkcs3 = colour_clusters(hcpn3, k=k) +labels(dkcs2) <- lin_type +plot(dkcs3) + +cutk1 <- cutree(hcpn1, k = k) +rand.index(cutk1, encode_ordinal(lin_type)) + +cutk2 <- cutree(hcpn2, k = k) +rand.index(cutk2, encode_ordinal(lin_type)) + +cutk3 <- cutree(hcpn3, k = k) +rand.index(cutk3, encode_ordinal(lin_type)) + +dunn((dps_aba1-1)*-1, cutk1) + +dunn((dps_aba2-1)*-1, cutk2) + +dunn((dps_aba3-1)*-1, cutk3) From a0123ed0379eaebd17c4e378b8ccf9644302e34f Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Mon, 17 May 2021 16:10:41 +0100 Subject: [PATCH 09/18] topo testing 2 --- workbook.R | 10 +++++ workbook_hemi.R | 113 ++++++++++++++++++++++++++++++++++++++++++++++++ workbook_swc.R | 8 ++-- 3 files changed, 128 insertions(+), 3 deletions(-) diff --git a/workbook.R b/workbook.R index 014c46f..2438947 100644 --- a/workbook.R +++ b/workbook.R @@ -60,3 +60,13 @@ image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n') title("NBLAST") image(dps_aba2 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n') title("TNBLAST") + +dps_listT2 <- nlapply(newnrns, make_sotopo_dotprops) + +dps_aba3 <- nblast_allbyall(dps_listT2, UseAlpha = T, normalisation = "raw") + +par(mfrow = c(1,2)) +image(dps_aba1 / diag(dps_aba1), zlim= c(0,1), yaxt='n', xaxt='n') +title("NBLAST") +image(dps_aba3 / diag(dps_aba2), zlim= c(0,1), yaxt='n', xaxt='n') +title("TNBLAST + SO") diff --git a/workbook_hemi.R b/workbook_hemi.R index 23c4068..1bfd1ca 100644 --- a/workbook_hemi.R +++ b/workbook_hemi.R @@ -9,5 +9,118 @@ plot3d(hemi_nrns, soma=T) sc_pairs <- read.csv(paste(SKEL_DIR, "../hemibrain_mismatch_pairs.csv", sep="/")) +hemi_types <- unlist(lapply(hemi_nrns, function(nn) strsplit(nn$NeuronName,"_")[[1]][[2]])) names(hemi_nrns) <- nlapply(hemi_nrns, function(nn) strsplit(nn$NeuronName,"_")[[1]][[3]]) +hemi_nrns_raw <- hemi_nrns +hemi_nrns <- hemi_nrns/125 +hemi_nrns <- nlapply(hemi_nrns,stitch_neurons_mst) + +htst <- make_topo_dotprops(hemi_nrns[[1]], resample = 1, k=5) +plot3d(htst) + +library(doMC) +registerDoMC(cores=4) +hemi_dps <- nlapply(hemi_nrns, make_topo_dotprops) +hemi_dps_so <- nlapply(hemi_nrns, make_sotopo_dotprops) + +saveRDS(hemi_dps, "hemi_dotprops.rds") +hemi_dps <- readRDS("hemi_dotprops.rds") + +dps_aba1 <- nblast_allbyall(hemi_dps, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(hemi_dps, UseAlpha = T, normalisation = "normalised") +dps_aba3 <- nblast_allbyall(hemi_dps_so, UseAlpha = T, normalisation = "normalised") + +par(mfrow = c(1,3)) +image(dps_aba1, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("NBLAST") +image(dps_aba2, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST") +image(dps_aba3, zlim= c(-1,1), yaxt='n', xaxt='n', col = hcl.colors(12, "viridis", rev = TRUE)) +title("TNBLAST+SO") +dev.off() + +hcpn1 <- nhclust(scoremat=dps_aba1) +hcpn2 <- nhclust(scoremat=dps_aba2) +hcpn3 <- nhclust(scoremat=dps_aba3) + +plot3d(hemi_nrns["5813027132"], soma=T) +plot3d(hemi_nrns["2009184596"], soma=T, add=T, color="blue") + +val1 <- c() +val2 <- c() +val3 <- c() +for (i in 1:dim(sc_pairs)[[1]]) { + sc1 <- dps_aba1[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + sc2 <- dps_aba2[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + sc3 <- dps_aba3[as.character(sc_pairs[i,]$query_id), as.character(sc_pairs[i,]$target_id)] + val1 <- c(val1, sc1) + val2 <- c(val2, sc2) + val3 <- c(val3,sc3) +} + +ii<-5 +plot3d(hemi_nrns[as.character(sc_pairs[ii,]$query_id)], soma=T) +plot3d(hemi_nrns[as.character(sc_pairs[ii,]$target_id)], soma=T, add=T, color="blue") + + +tanglegram(dkcs1,dkcs2) + +library('fossil') +library('clValid') + +k = length(unique(hemi_types)) + +cutk1 <- cutree(hcpn1, k = k) +rand.index(cutk1, encode_ordinal(hemi_types)) + +cutk2 <- cutree(hcpn2, k = k) +rand.index(cutk2, encode_ordinal(hemi_types)) + +cutk3 <- cutree(hcpn3, k = k) +rand.index(cutk3, encode_ordinal(hemi_types)) + +# test on neurons represented in multiple types + +nrns_types <- data.frame(type=hemi_types) +multi_types <- nrns_types %>% group_by(type) %>% count() %>% filter(n>1) + +m_hemi_dps <- hemi_dps[nrns_types$type %in% multi_types$type] +m_hemi_dps_so <- hemi_dps_so[nrns_types$type %in% multi_types$type] + +nrns_types_m <- hemi_types[nrns_types$type %in% multi_types$type] + +dps_aba1 <- nblast_allbyall(m_hemi_dps, normalisation = "normalised") +dps_aba2 <- nblast_allbyall(m_hemi_dps, UseAlpha = T, normalisation = "normalised") +dps_aba3 <- nblast_allbyall(m_hemi_dps_so, UseAlpha = T, normalisation = "normalised") + +hc1 <- nhclust(scoremat=dps_aba1) +hc2 <- nhclust(scoremat=dps_aba2) +hc3 <- nhclust(scoremat=dps_aba3) + +k = length(multi_types$type) + +dkcs1 = colour_clusters(hc1, k=k) +labels(dkcs1) <- nrns_types_m +plot(dkcs1) +dkcs2 <- colour_clusters(hc2, k=k) +labels(dkcs2) <- nrns_types_m +plot(dkcs2) +dkcs3 <- colour_clusters(hc3, k=k) +labels(dkcs3) <- nrns_types_m +plot(dkcs3) + +cutk1 <- cutree(hc1, k = k) +rand.index(cutk1, encode_ordinal(nrns_types_m)) + +cutk2 <- cutree(hc2, k = k) +rand.index(cutk2, encode_ordinal(nrns_types_m)) + +cutk3 <- cutree(hc3, k = k) +rand.index(cutk3, encode_ordinal(nrns_types_m)) + +dunn((dps_aba1-1)*-1, cutk1) + +dunn((dps_aba2-1)*-1, cutk2) + +dunn((dps_aba3-1)*-1, cutk3) diff --git a/workbook_swc.R b/workbook_swc.R index 8e7802e..d616396 100644 --- a/workbook_swc.R +++ b/workbook_swc.R @@ -1,8 +1,10 @@ library(nat) library(fafbseg) library(dendroextras) +library(dendextend) library(elmr) library('fossil') +library('clValid') devtools::load_all() @@ -13,8 +15,8 @@ plot3d(nrns, soma=T) lin_type <- unlist(nlapply(nrns, function(n) strsplit(n$NeuronName, "PN_")[[1]][[1]])) nrns <- nrns/1e3 -nrns_s <- nlapply(nrns, function(x) elmr::unspike(x, threshold=5)) -nrns_s <- nlapply(nrns_s, function(x) smooth_neuron(x, sigma=1)) +nrns_s <- nlapply(nrns, function(x) elmr::unspike(x, threshold=4)) +nrns_s <- nlapply(nrns_s, function(x) smooth_neuron(x, sigma=0.8)) nrns_cl <- prune_twigs(nrns_s, twig_length=5, OmitFailures=T) plot3d(nrns_cl, soma=T) @@ -60,7 +62,7 @@ labels(dkcs2) <- lin_type plot(dkcs2) dkcs3 = colour_clusters(hcpn3, k=k) -labels(dkcs2) <- lin_type +labels(dkcs3) <- lin_type plot(dkcs3) cutk1 <- cutree(hcpn1, k = k) From 4cc24c608c9490ce5d6dfd0e65c5538f9712ded5 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Fri, 28 May 2021 13:53:15 +0100 Subject: [PATCH 10/18] utils for nblasting / clust --- utils.R | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 utils.R diff --git a/utils.R b/utils.R new file mode 100644 index 0000000..d24d136 --- /dev/null +++ b/utils.R @@ -0,0 +1,27 @@ +# utils + +get_dist_to_soma <- function(nrn) { + gw <- as.ngraph(nrn, weights=TRUE) + dst <- igraph::distances(gw, v = rootpoints(nrn)) + as.numeric(dst) +} + +make_topo_dotprops <- function(nrn, resample = 1, k = 5) { + tdps <- nat::dotprops(nrn, resample = resample, k = k, .parallel=TRUE) + tdps$alpha <- get_dist_to_soma(nrn) + tdps +} + +make_sotopo_dotprops <- function(nrn, resample = 1, k = 5) { + tdps <- nat::dotprops(nrn, resample = resample, k = k, .parallel=TRUE) + tdps$alpha <- list() + tdps$alpha$distance <- get_dist_to_soma(nrn) + so <- strahler_order(nrn) + tdps$alpha$so <- abs(so$points-max(so$points)) # normalizing so the main branch is always 0 + tdps +} + +encode_ordinal <- function(x, order = unique(x)) { + x <- as.numeric(factor(x, levels = order, exclude = NULL)) + x +} From 71638765a9d118e97dcff51a726f517f1cc6f464 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Tue, 14 Jun 2022 10:52:27 +0100 Subject: [PATCH 11/18] new version of topo nblast scoring --- R/neuriteblast.r | 49 ++++++++++++++++-------- man/WeightedNNBasedLinesetMatching.Rd | 12 +++++- man/nblast.Rd | 4 ++ tests/testthat/test-NBLAST1.r | 7 ++++ tests/testthat/test-NBLAST2.r | 11 ++++++ tests/testthat/testdata/testtopodps.rds | Bin 0 -> 45550 bytes 6 files changed, 67 insertions(+), 16 deletions(-) create mode 100644 tests/testthat/testdata/testtopodps.rds diff --git a/R/neuriteblast.r b/R/neuriteblast.r index bb23971..1060d54 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -46,6 +46,8 @@ #' @param UseAlpha whether to weight the similarity score for each matched #' segment to emphasise long range neurites rather then arbours (default: #' FALSE, see \bold{\code{UseAlpha}} section for details). +#' @param UseTopo whether to weight the similarity score for each matched +#' segment to relative topology of the neurons (default: FALSE). #' @param OmitFailures Whether to omit neurons for which \code{FUN} gives an #' error. The default value (\code{NA}) will result in \code{nblast} stopping #' with an error message the moment there is an error. For other values, see @@ -164,7 +166,7 @@ #' } nblast <- function(query, target=getOption("nat.default.neuronlist"), smat=NULL, sd=3, version=c(2, 1), normalised=FALSE, - UseAlpha=FALSE, OmitFailures=NA, ...) { + UseAlpha=FALSE, UseTopo = FALSE, OmitFailures=NA, ...) { version <- as.character(version) version <- match.arg(version, c('2', '1')) @@ -187,10 +189,11 @@ nblast <- function(query, target=getOption("nat.default.neuronlist"), } if(is.character(smat)) smat=get(smat) NeuriteBlast(query=query, target=target, NNDistFun=lodsby2dhist, smat=smat, - UseAlpha=UseAlpha, normalised=normalised, OmitFailures=OmitFailures, ...) + UseAlpha=UseAlpha, UseTopo=UseTopo, normalised=normalised, + OmitFailures=OmitFailures, ...) } else if(version == '1') { NeuriteBlast(query=query, target=target, NNDistFun=WeightedNNBasedLinesetDistFun, - UseAlpha=UseAlpha, sd=sd, normalised=normalised, + UseAlpha=UseAlpha, UseTopo=UseTopo, sd=sd, normalised=normalised, OmitFailures=OmitFailures, ...) } else { stop("Only NBLAST versions 1 and 2 are currently implemented. For more advanced control, see NeuriteBlast.") @@ -370,6 +373,8 @@ WeightedNNBasedLinesetMatching <- function(target, query, ...) { #' WeightedNNBasedLinesetMatching. These will be used to scale the dot #' products of the direction vectors for nearest neighbour pairs. #' @param UseAlpha Whether to scale dot product of tangent vectors (default=F) +#' @param UseTopo Whether to use topological information to scale dot products +#' of tangent vectors (default=F) #' @param ... extra arguments to pass to the distance function. #' @export #' @seealso \code{\link[nat]{dotprops}} @@ -380,12 +385,18 @@ WeightedNNBasedLinesetMatching <- function(target, query, ...) { #' segvals=WeightedNNBasedLinesetMatching(kcs20[[1]], kcs20[[2]], NNDistFun=list) #' names(segvals)=c("dist", "adotprod") #' pairs(segvals) -WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE, ...) { +WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE, + UseTopo=FALSE, + ...) { if(!"dotprops" %in% class(query)) query <- dotprops(query) if(UseAlpha) WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect, alphas1=target$alpha,alphas2=query$alpha,...) - else + else if(UseTopo) { + if(!"topo.dotprops" %in% class(query)) stop("query must be `topo.dotprops`") + WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect, + alphas1=target$topo,alphas2=query$topo,...) + } else WeightedNNBasedLinesetMatching(target$points,query$points,dvs1=target$vect,dvs2=query$vect,...) } @@ -396,10 +407,13 @@ WeightedNNBasedLinesetMatching.dotprops<-function(target, query, UseAlpha=FALSE, #' @rdname WeightedNNBasedLinesetMatching #' @importFrom nat dotprops WeightedNNBasedLinesetMatching.neuron<-function(target, query, UseAlpha=FALSE, + UseTopo=FALSE, OnlyClosestPoints=FALSE,...) { if(!"neuron" %in% class(query)) { target <- dotprops(target) - return(WeightedNNBasedLinesetMatching(target=target, query=query, UseAlpha=UseAlpha, OnlyClosestPoints=OnlyClosestPoints, ...)) + return(WeightedNNBasedLinesetMatching(target=target, query=query, + UseAlpha=UseAlpha, UseTopo=FALSE, + OnlyClosestPoints=OnlyClosestPoints, ...)) } if(UseAlpha) stop("UseAlpha is not yet implemented for neurons!") @@ -475,16 +489,21 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL # zapsmall to make sure there are no tiny negative numbers etc. #scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) if (is.numeric(alphas1) && is.numeric(alphas2)) { - maxlen = max(max(alphas1), max(alphas2)) - scalefac = zapsmall(1 - (abs(alphas1[idxArray[,1]] - alphas2[idxArray[,2]])/maxlen)^2) + scalefac=sqrt(zapsmall(alphas1[idxArray[,1]]*alphas2[idxArray[,2]])) } else if (is.list(alphas1) && is.list(alphas2)) { - # distance from soma - maxlen = max(max(alphas1$distance), max(alphas2$distance)) - scalefac1 = zapsmall(1 - (abs(alphas1$distance[idxArray[,1]] - alphas2$distance[idxArray[,2]])/maxlen)^2) - # strenghten by Strahler Order - maxlen = max(max(alphas1$so), max(alphas2$so)) - scalefac2 = zapsmall(1 - (abs(alphas1$so[idxArray[,1]] - alphas2$so[idxArray[,2]])/maxlen)^0.5) - scalefac = scalefac1 * scalefac2 + # use local neuron topology + if (any(names(alphas1) != names(alphas2))) + stop("Topo features not consistent between neurons.") + compare_feature <- function(feat) { + #' compares the features between points according to metric + #' 1 - ( |f_1 - f_2| / max(features))^2 + maxlen = max(max(alphas1[[feat]]), max(alphas2[[feat]])) + .Machine$double.eps + zapsmall(1 - (abs(alphas1[[feat]][idxArray[,1]] - alphas2[[feat]][idxArray[,2]])/maxlen)^2) + } + sim_scores <- sapply(names(alphas1), compare_feature, + USE.NAMES = TRUE, simplify = TRUE) + # Aggregate feature scores + scalefac = apply(sim_scores, 1, prod) } else{ warning("Unknown alpha format") scalefac = rep(1, length(dps)) diff --git a/man/WeightedNNBasedLinesetMatching.Rd b/man/WeightedNNBasedLinesetMatching.Rd index 55e0b78..21c5483 100644 --- a/man/WeightedNNBasedLinesetMatching.Rd +++ b/man/WeightedNNBasedLinesetMatching.Rd @@ -8,12 +8,19 @@ \usage{ WeightedNNBasedLinesetMatching(target, query, ...) -\method{WeightedNNBasedLinesetMatching}{dotprops}(target, query, UseAlpha = FALSE, ...) +\method{WeightedNNBasedLinesetMatching}{dotprops}( + target, + query, + UseAlpha = FALSE, + UseTopo = FALSE, + ... +) \method{WeightedNNBasedLinesetMatching}{neuron}( target, query, UseAlpha = FALSE, + UseTopo = FALSE, OnlyClosestPoints = FALSE, ... ) @@ -26,6 +33,9 @@ compare (must be of the same class)} \item{UseAlpha}{Whether to scale dot product of tangent vectors (default=F)} +\item{UseTopo}{Whether to use topological information to scale dot products +of tangent vectors (default=F)} + \item{OnlyClosestPoints}{Whether to restrict searches to the closest points in the target (default FALSE, only implemented for \code{dotprops}).} } diff --git a/man/nblast.Rd b/man/nblast.Rd index 74a7aba..f3fb1f6 100644 --- a/man/nblast.Rd +++ b/man/nblast.Rd @@ -12,6 +12,7 @@ nblast( version = c(2, 1), normalised = FALSE, UseAlpha = FALSE, + UseTopo = FALSE, OmitFailures = NA, ... ) @@ -38,6 +39,9 @@ query} segment to emphasise long range neurites rather then arbours (default: FALSE, see \bold{\code{UseAlpha}} section for details).} +\item{UseTopo}{whether to weight the similarity score for each matched +segment to relative topology of the neurons (default: FALSE).} + \item{OmitFailures}{Whether to omit neurons for which \code{FUN} gives an error. The default value (\code{NA}) will result in \code{nblast} stopping with an error message the moment there is an error. For other values, see diff --git a/tests/testthat/test-NBLAST1.r b/tests/testthat/test-NBLAST1.r index 3be6293..0e25d72 100644 --- a/tests/testthat/test-NBLAST1.r +++ b/tests/testthat/test-NBLAST1.r @@ -25,3 +25,10 @@ test_that("we can calculate normalised nblast v1 scores", { scores.norm=scale(scores, scale=c(scores[1,1],scores[2,2]), center = FALSE) expect_equivalent(nblast(testneurons[1:2], testneurons, version=1, normalised=TRUE), scores.norm) }) + +testtopodps <- readRDS('testdata/testtopodps.rds') + +test_that("topoNBLAST works correct", { + scores <- nblast_allbyall(testtopodps, version=1, UseTopo = TRUE) + expect_true(scores["18820","20262"] < 0.1) +}) diff --git a/tests/testthat/test-NBLAST2.r b/tests/testthat/test-NBLAST2.r index 0b5b849..bdea92b 100644 --- a/tests/testthat/test-NBLAST2.r +++ b/tests/testthat/test-NBLAST2.r @@ -100,3 +100,14 @@ test_that("we can handle all combinations of dotprops and neurons, both as neuro expect_is(nblast(Cell07PNs[[1]], testneurons[[1]]), 'numeric') expect_is(nblast(Cell07PNs[1:3], testneurons[[1]]), 'numeric') }) + + +testtopodps <- readRDS('testdata/testtopodps.rds') + +test_that("topoNBLAST works correct", { + scores <- nblast_allbyall(testtopodps, version=2, UseTopo = TRUE) + expect_true(scores["18820","20262"] > 0) + scores <- nblast_allbyall(testtopodps, version=2, + normalisation = "normalised", UseTopo = TRUE) + expect_true(scores["18820","20262"] > 0.9) +}) diff --git a/tests/testthat/testdata/testtopodps.rds b/tests/testthat/testdata/testtopodps.rds new file mode 100644 index 0000000000000000000000000000000000000000..d09bf61369b0e4862405adff995f7befbd43eaff GIT binary patch literal 45550 zcmV()K;OR~iwFP!000001MK~IJeA-3@PSgsCL)q1B8mp35M^BosSqg?GS8d6*=CUj zrBaCogiM)Arb4D=%q;Uf&*R=oNkZrD`|I<2&hK}g=b!U>UgtUIJm2+t-My`KuXU}p zuIpOU-uJ!JUBke@u$qC1k@%OH_>+P8(naO{GDM8~Pvn1a)dqrvVD_QQB_sDXnAPv@ zk9Cg%vtjLBXO4XVv+2*v$*+>Z?E9UL^OU_{w#fh9!HmfFWslt0`S)NpzLZX>I1FZ! zCDNr&)_~d9_73Ll!eGAYy_BHUelX{G%Egj+5zM8O(@)&61#{I9lY}&RFgH28()$6x z{Jxs0_ERA+&ldk$-Cz&q6IVCiVs-!v-r}?6Y^%ZI%xsM`pD|dt7HvvCwHqwLSjtlF z1%pL)$T_iIBE4%%{Yd|Qu=ue*bK5RSu;h8I|9r#-ERTFvAU#_TmKXCA?i+A|rQHWZ zQcojT22OJL_j6Etd06oQqp zz2(fs1+WUdXVz&d3RZ>mj4X#QVD-!D{L5i4uokQPOpbj4)~a2fWbfPuYmaxw)cWjaVy*{@B=j5!Ent1Dj~t@8h?1z^0{TI%oSG zuwBodRPreuY!xi1rRN;L*2%W${k=-CP2!n9JJ|)c-A|oIx1I$OXW5J8Q#*ikn#b$C z_d6hYU4D9EtOrPG>s30;F92yKFxhR}Q?O&Hw&cy%2RnurjSpYT5Nr_>pb3Hv;Kut?*JSJ31#1ABz#qK)tR@<=(&oRQK~|>lgfi>f~*E^CmY?$@||% zpDqWgiL*y~Y6DQO?s-y@u^p&7JE9)&_W||%2DO9l%7ChT=R(=WNT8lvAFaPv7pSMC zEY7Vu4^-I)kwuhXpdRYS4eq=^-GAMphxsc|#ic&XY#RfrNZ#%%oSi@wt68h0ZARR$ zbYMC|1*&AqsMYFfg0}ak<28UPEzm2Kv_iquho*MuZ_ug z^dnIBeGEMs7!OoINkcaYLLXc`Vb9jA1}aOFpqZaGIIhI3(^MqnGdoD^EoB481-X}f z#hl^EO|riFnmuSrH>RW-QIm zn9qV^vi?`@m2q&4*|^^R!gp{COe|u*5D$)i$!oQao(D&_quaTeo59iik?c}gDmd!g z=1!F*DJBENj_;E;G_;X25dJictza5cBnS;-ZD6vzEh=`>;Xp=Dmho|Bsd;;KzG+X0FGQiR`myV zgTs%IPM00_;Ltt3w&A)UIF!xG42S;$hd1vBQp|h7AJU(F zSBn*;@B!t;;=+k{B|v%2ReqrHJ5b)A+$R0v0#I;LVNNhJQ1bO^IaTREDfseGVAXA) zBwK#Ev-t>6p5DATE~@|(*X8bW95;b-#dnhZ!Wp0_RQxPs^=DvUT1Bh_j0{W+?0*&x zGiy5n>xt`6GW=^T`L9{nh%ENTcIH&#Z|=YOR)~x9Kh~XpP*W;8AQKw?n>{X6(R4* z;jOkG5xuv(GpCG-=zS4tKMI!+y)e}~*NXws=Xm0L-7X+{k8E#*kv*cPe-IeeBl?%s zlkn_IKcf3~g?{WbK(u4eSmSiE5$(9=nRv@hh<5H^hQzu{h{h~Fy31n|qC4CiFT6(7 zzh^)1S_3IW`}A;m3~~`|Vm`9<{w6fjn>`?PkPXo`_bfadWI=R&KlYq`Nr)DIN$z&U z5kxl$>UR(%^meOivtq0pnu+(z?Q4rbG`ZXiftw*{x;(BoS#SW+az^CE)ai&m63 zh8-h_zT>Vy{x)Akk88H@3I2ped!KgPIC39NNE!U7veQEJww87qgC~eyJ(p)5evZ(4 zcd({$!aw|H_kmgph3bsYdHyFG9P2-H8w#Pqv858_8)s13 z;UV4P+slXY(D?m_6EclY(8SEa*4=KK(Zrmnt63sHnz0oS$tfE_)9ceqVx8sD`1Oqn zaWTnggyoW9`gRsH;j^)8X7>g((Rb>`o6xwLEStxiVM`bLJb4Xe;p+p!eVq4Q{pG%Ua7p#Yk!y1tn7 z>jRp+o9`lej{!}NFLi~A`J<_xSNNjmBQ)**G4ljJ9Zh$xdNOw?9gXTe^UD$xM@?ft z_1*J}akS>kiWrGT>}hmR!JK^;O0?W`s;_|KFW2D;0gj58r^o+N5tp--R@Q+Z_Bd~# znK$o)6AzrfnR(v`*E*{u%^%Ul&7`GZU zbNAq;w=rHjE==L(&&u14#n0iEXhBkx8YeFI46CbSWW(hrO9Lv+8E{+jw*^JLJGlE| zXVx1(E?l3-_iE>&E)LtR!1m^OKTcZKW=yE6z@7f_OozG#aqsfh@BZ@HxaWn+xgE8Q zxQ8dYBC&l3uHMX4Tc5~@LmqJHmxLAKy6=XsJau((!&|O?s&*J|e)V#pPEHkfGZb+v zZc@Tskv}HyzgfbAtXZjC%B*-;yEslPsR<9**fB(x@ZoN%E@;&<;Qoe7&b$@5|H+wa zTjo6m=l}9MrpX7*&;~AHdNzkmfQ2@u56L<4M2ukiTJ?vmGxsr#yEOUJAyqt`dNQ_; z)e_Ha*isdx?~iA0o^dtOlEX8TLEjS&{Khk^zDE^LQt%X|A?B>m9z3lV(sd^3EuKo; z{!V*l8IQ%zXx|PAz(Z%XDX+5m@lauN@#HBJJoLUR_?1Nr9x5}!BHP08h-tXiZj~iG z`SpnN(~-M)%D+I^SL`^RdZTwW)#V+Y8oiyIkyL}GM3}z)ln%mEtKPjbTu(z&ferVx z`p=@Nr=>TJ1TLbfcRR+Tj$OtRcZ6P_eeH$EGQOu-r5{CO7Rhrt3@&Ii%Sq)&`T!nY z%k@0l@D?7tzPX>Z)ERdRw(>NHzsH{`Gv8Nci*ak3j@c1IZai@NJ7s2u6OUP}H9D~{ zjR(mceE!~_anG5yq7bfm+||Ev^lkcaJXpEnF!t*=9%|UU%5QrU?kORqRGof|KRv71 z__g^lF8lPW2pD?N%jRs0+OT<4;hJ|(deaTm%UNI6xfH=dx_;SI=0b_QKJ_HQ;9uYhhgUBlBr4KB+xp86hs94=>H5Pg~)4_Bnimp7X}fh$?c zd!<;qK<|h$%SGKm&`arV`u*e%=pTF(#AQVW{e+p5B2xljAmTKeC_e!Pukq;cF%!7D z`>wB#nh0Ed%^9%eOT~ZpgMWwrTf=`jmNEZxCgWuIGgfe8}GylEKLLpn*J5m3UXiTD58UKwixtN2O0IKH69CgAN@v{z9)Rg85vlY5I+%QVwA1F)PcXUnUeO|k z3rsGtG-Q^(1CxCck8XswfC(Fa--)0oFdh$laICHljGL~7sZX(kan=48rIb`KE*aZ; zDpv@MD=x|0dJzQ1^}hllhgO4eL$^;|x&jziKOAw5I1a`iLn+@SBEdML?Yw-pEf^=x zjekG46O1EXyFD&H4aSeJlXC7{1Y@t*rlD+mFn0QCX&eDyOgbf%GyGJ4#u*f1)g88gRxBgw)yyo#y;b&z8Lz69{iTrFbtfyu)z?hM@L~$_@jFw5i%0Im#SWceb^%IO) zWX+bzd0E}A(N3(Vu0^?sJ zZtmMIfywS;yVHfPfQfq6U<1Dim^itSUnF#c$?cHG_1iyz$wO)W8>=J0BqY*@#6sjt zU=2C<{tK97U)Vd_xd%+DLL>cnn!#jbQdGQv9ZWeyS>R_4m?~4ZySQ8cQ=cD-YF}4@ zX@17~v&o!b`n_*kUT7((^1@r37mHCEO!D1Eb7jwf-usC***Ih{-EF511 zgp!)T;%(|RRoiV~F)lm)QR5_7ifk(u9i)S$xpF^WK_ytmYDfJr%m&N;%p2|J9)i_Q zW%BBXXt28aLL{xM8?3_9?ya8B2CMD`Mt++Cu;v$|o`0?l)|VnT-PIZb>&N{kM{j9> zb;*xjhTo}RO~2}P;pJ1X;i+0NYY7D#Xo%gHJw1>3YQ+)AOR%#c-iH5ia>IF^deR;1W55a_XxOe15)3H z1;o?{c3ZXe`t|a_PU%A~TFVD^4kJAIP0zqCynXT2yYpaII>z9N}&tT8`h?H*V0d~Xc`|h}jfL;E{cb9dZfL+jYktfYc zU}w9#<<>+P8HaM1gc$=3!CNIUrHYSMB`B0wfEKq~Y)NKr)}lx5~|dWZSigx$F~=?DI|^ zToD2i^}SK2RU(kwM))K@djsiO0w3RWGmyL_l~29w1-V;{r-N&DB1i_aW@q||)! z&bKKbNiB|s&~5=qaz6QmJTs7l2iby^y@0eaKfbk@sCPAAR-c(a*v`{T_I$ksw!`=R z=QehNZM}-t_Jb*%3yM{TSy2}3AH&_Ey^IeykD^>&bl0NO1!=F$5$3-Ut4+Hho`nsGT z1XNkIofj`K0`;iG?N@cb2{I%X_znQ|WO3!KjQd2n)1<>ojf9+;ixbO4+&;}qmFFl> z?W{SDwd(`bgtu=*xC^Li>!QNn9|x+`{p&}yS%JEPO*LZ;A)nc6Vf}|raGWVB+03mC zj%{~W3)&}vW3I8Zp@}0nM*L16&_dvN%bw5oT_!l%lFEdibAh8)q${mo860Kv%#LZW zfTNI);TEqVa9o#HZ93Np4&NVs7yQTu4!t>H`RC=qq1yGl64*^OQjiJ z7g>NZtb<&VDu7aXyuwB9I8Z*EaJ$<58Ypi^ndI2N10@k2mKK`=C8mPCKH3{7&(BJR z^78}buG-?U=V?Hp{OH&po&=OD0VNyE^?{<4TBX$~02GOJ+TxE@fx`P>$@|hVATP@v z%;|_l4Z_q{4w5 zX>wb3A_d4FHCXxmh;}#WwO+iK3FMi_WvSE=psYG2ccIQ5C~Fwyw-<&2WxaWP;L%c` zY^d|x_F^xQ&g;;U)(RA6|J6@B!hphp7lU$`fwFF$x=Fh^P&fnCRXhR-y*O;ENG0sB z?U{Dp&YeIJ+&b~JYYk9DnT!MkxPT(ob98k!k-m?HE?qqf6!D3RNxiRtvPa_}shQBT zuyu$W&oiJ17;WB`=1<6TUPX8OO0*+hsm#}ssCWC-d*5q;{Bx@B^^pvs{V%?#(dvLa z@QbC@YL=jr5J&e*Adg)-bLdbikZB8%+F@IO{56X-&p8a_-dzD#(P|(!hE8WaBJ5L= zt1f%o1jrfm!l`O5Ais|I-p{!M$YI|_zVcTB`Dw`fNlxPaBi4q(glQl@+$Jrg`W(nF z)?4=KD*-u}nbd7$Nzm!A`3E*2#||o=nEytUegvQ^y$!b^_VP zna8)_0FYfg*QM_!;tx6Q-M@30sBhOU{(FSK#VKTX&94yg4_uoFzb1!dY70;M1Nj~; zJMr^FAiGLRZaAk5WD|!#%*IdD`zv@p+XKjo7o`0*GXYuF@9vID!at6@RDEh917t}t zo8MW){R3yCug=8*d6#mjrl%H=xmEMi-Vp6z{+ayqgA?JGaS;kQ7wkvh8tq>s{I{7o z+CAhK*f+eez8YB#_T^`AUl}vl=Z|{8j1t&qc7zhl!F3?c-wbNIyvsMVyW{w3r z`KvqB-MzqW6RUO77h51riVZxNu?A9Epi-tl3y@;y{kOA;b?`Pvv0%Fvkc|1xX5Rfw z><_jVVX;6U3C0=kHkAU>1|g4i!<3{( z1+V8Br-MS&!0WC==-n1M@N!)}li)fJUgq19(xV2!OQXEtRpmkOl5Kj?Q~w3L1Pews znK6JDr|6Dh-&NqnTyZF*Y6`pNv&$~Z;6 z0G?q=N3$LWfahJS^M~{0z>|8Q{OGA|;AznKxJOABJXIul8cABE&Za;b3 z$;=i!*F4XVRB`~1Z$*0re|-aw9ua~0(mL>{aD6oP?k#wvzt=rqECe25kFGgyNe2(# z_X&qfI>E!mCR6t%aL>DOR_%NfxTiGayj6Pv?lCeB{vJQTJ!CDfRhbdE zKkaAX4(BJB(fLYIhPbc%++sEh+%p>YALC8{_fNbUp+d*Nz4PTh);BxBotETpJo^be z7+yQM$_s)A_k(jEjL(6GKz{r9i9GNS%XAnUdJZ0v;XXet4ui)DyIS$|OW>h$FwIAn z0v?y54Ua@cfQKbx(HXM|@NlEGZIQnL9`~Ex4k<>1M+kRiuyYM~#11$|-nsxDDbEw_ z8=b%-i*Y32h#h#8Y&xaPP2{U3$Ja7#1&_wUL(lvOy*A$+=6fklkXN2j&X7na=N4Y$ z0gnbvXB0359<49xQVuE*$RDBH|)oHxtw}J^f)H!hddPvao^A1N_g34QV2y6t8 z!k>wLtJ%RLv+thMIih_@u4Cig8yEe+W6QhOZXpifzH&jyLdy-@hXfvv=KFwqWr+z3(OzH&*a64%(FO?h(ZlVWAHk2I(H+K4!3m=@pb><6y=PPM&ZFpFzwJjN3 z-xoJcpCf~7*h^`R${=v{-u84)fgHFRzq>7$TMVvB^~tI{2H+~TS(7VP6kIu6o(n3o zfXla*u7hJE;L`0zH#~L*T*{Jbz>OD!UN!9EXs88 z&JwsNjkJx(z6TdEh1w+NE^yhP^WYR^D>#2OooP<+2Ip4m(I&IC;QTS6Ea8MLIERKz z@K)7>v*)>ivA|?-){CBJTd)Raxp${`rT{o^ok?kU%Lh(N5rX?B{J^Q_`I`zuYj7$! z>v>?`3^+xsoDEuZ04JZhSBj<5;AFb_mQ7a$IGr*+EUQoePC`A>nUl}JiAibwhA;!5 z4!esN))M>K(r>qwGB|+xdc6Dhi%rD7wq{{nQ7N(ST*-VG7Yb8mKarB~5XA zfV%xT*0wGM$KSn72Gq;o*k5Ed<3I4mEXYr1+59HzYdcf3dkhqlFg z16t$YP*_Bl+hYa}ajBNnb@t%!NJHVO7_q-}VQ*d_p8*HG{7jK{Q*b!{^`PjHh5yNC z{do-k`FO#qKk5H(K2DH_=)+DV3%cN63L*Lma+Ra+K=cbOg<)}z5Z$=<>+*ayqSqXm z{m?0e=!0c_TY5eqddbb-jNbHAlRL8z zElR=DBHfCp$EHl>^jAc0+B?@x4M6nSH6~_?k%->@QBRfTh3Li6efZiCqRT|Qjm%L% zbPrDdy8av8ZzujY&zso%z$o@d|%BP!FZQMolzv6hm$~r<{eT6T# z_9A+j@!F#PlY|@*2|XTBL{GZ)=BwW+ME|mB`Vp4~A;(hfe)| zM)Xd0(x!XjSMy>{g$beOhI|F~w}_r*VrUUav`bpe_28ET z!oCa{i0=fVXTSfVS5<-Nb(Ydn8HAki7q3|E<|2CLa5(Gd8;IWcDAw!FSw!#4XtTU# zKVBdA?hE~SWs^y z?BCZs6?dYL(Bqj=xw#%fj|w072eJrzHz`>^`AM|@roW|$8>08z_A_^>LG+AprqWmH zh+Y^M9dYjoA?Ilg#gFi_=jYxZXlz1s**iZI(+U5qdl6?o`jqIm*wu1&J7RoT1Uh|M zM)douSf5xsAle(>VqyJt5p$W0=+)X%T%N@Ei2fS&?bmlizc&9=O>_d$G9$K0GrU5yhqEj(SxboSrrtlh zp3q~g%74ElA!_BF0OlW%bjj3y8i*d3$Bo zH$?xmeM^2y1EM>X250aMBU(@GIySV5m>;Nd{*8MHKO0C-=U7A7L9Vsd^%bJe`>9m4 z5aamAyG7^9dLlijDBzqmkw3e9~816#!+<=~)%Y^^XTPj?0h6#FeN;wi7 zkGQW{NsQz9fal&c!k&xd-JSbK2s;eRJt!L>{BmhWnrRB*R|}&#BF>1OT%dJ~@N;^9 z+|2AxZA4!Pxodxj7;p39HuUGjJn+@1BJu;_hlPzAaiIl-z7kIT)}J8Sm2b1pd77Ai zOh%XOi1v4?x)ojsMD(fwYOf+Op9Jh+AO7r#=&vVcVjYNnjPRP?c~VU1%l&cw@gzha zXh}${C;V6=s#hYR716Ih|2VkQh#0@0`9`wi3B7Q0aGDYN7`a?(ueFWvgQd>{$Bz-| z$8&FO=OO01p&zrjAJLy!s5^c;is+&DHY1QB+8tM%{DFe#KetwH|aljp4>%EaZegUv!OcjC_#k?S9-(eS4$p#>bLwQ!kYjlA#g08}}vIv^oG zjq7SGyuUosz}UHkLruaOzwhT!+)L`fnI@*9*<=^|(JZ@d3<={85`2d(+B$Km&`r-^ zYc2Fy?`}wySv;<*l~I=oTaELyKD_peJB$P8E>8Lv3FAbiGw*HeF5)cy9ezWK_Bc>9 z{h46mI#lR9byHqy3e`Ui+}^fz0Cnm-VV|ooNA34lg?&-%LnD`tm(6F%q4A{?$Z$yj zjl9c0@KIa>4JB@CRwm}HA#GXfFlA0OT-4rp)5Qc0sqxuv)%b+EG-zUcR(Db3%_zCB z*JY@r;6!9$g9NJKm{eyvatzgR-~pClCqU~X(8HqKbDVplbm1Yk#fs{>k!i({ESy69^-~bx; zX7$k+3Pls^EE<{5X`m^VUsZ>!ccZDtXH&F&LeP|}ckK+Z>v`k`o-G{a%ErZ_I-JoP+}$BgCtw`8BkBVnZZV8>l}(Avy+TCx-k^aS6U=aWH$ z#v$yd*fyYkRFE8ZfR3}%zL__QORy-O{;-TQWv$2|Cc=Yh)uppOAH0n6#dVTK>G?t@g z+i-n7n#l0salP(}rg}%d={T)H(*yTr%6GP*318i{k^$yu>;T)mYW7PsQuo4$r*H)g?3DQ>J?4sgqTX6r@Snn6(tDGtx@U3s)hs=y^>oygn|#~h z)Dmilbjb`L=Bc8?EkD#ASE2Hj{rkNd_P6pdwr~Anf4O(C=G~K$j{*p>)>bc{=vANv>+k8)DKFX-rj7La9D2-k2SElV7l%aI;<*LeFl(F@9guyHgz3%$r z)_`u|N4Bc(7LVM<`SHo_+oQa3`Nohv2|7kNC*XO_nr%)vE#%Ej?~lGXaqbxBg_qlK zko1?KgI=!qj--FLx@#fMytSI~g3%gW+;b@Esh2aZF^-M2bF;^Vy@xKSt3Jm?`@#lK z$%NrL*}56Vj|#Y6k(}sFQN%U0U24~TD{=khm@2jxnz&IlOkRHUHENA{yYk(+9e3_% zXl+s6h`KfwRXnF}MV(D5KRCtBaoYkhkFzS^Rs%K(k)wRLBlmo0{`?5;ZcDImKYtq! z@cL<}QpmXH=xT9?R&U(FF2;LRvP3V(g!DZ zwtt9Pm~M*fH{-?su+4LN+*{cQ_ZyV#-Nh-6hYnh-Va;K|L*Lm@L;eUJJng7tUL%D2 z&*e0_Is4-F?c#x#e$!B+m#9riMJ#Se?{7Miu83u>P} z8*oqMp;t%4lkh;8^5bpE_PAwn;QWBlCtMk+e_ z`?}M=@R*w2^|A6vJWiJHuGzl@kFN<8N;+qb$DZ%D*RI@&Cr9sCoRu5I)1G_ZBYAc_ z!*JlKlYTLttevXYuFJ;b4gZjXAX}u%O17Zb03Y}m14eCvIUQ=lm0eE z--X96S^u8h;>pqId z3hvuixToU5=~9&qN>6azxbLvcaVuQzic~f@E#P+ksI&1+yKwiW)O8BHNAU0_f5E*c zw&0PtDzC3sH{($znOpYnH{em`Z>PRnlo5If{4Skcho@(Cn60LyFl|+}hA4RoPjzSQ z?~`!C6Fvtu$3uzvIGp)u;H6$XCeX%m$&wR~x~Mg}clh9beW^HoyFZUVuH=gH_u}#$ zxhAVN=iy8Tc2BPga>P&_JnEw`)5s1m&D|C|IY!a(?UZrr)&weWl+AX zpL&HmK5sTz%R7lWg_D#GRra8^2W@_>l?$k0eSDCF;VM+z{_RHUw;&Yz=J0vX#C15E z(figJ{wZ9TnReZjth!ypYl1RC5y_Mbx2lEa#6Xm#+Bxn zOjLGd*2CK4HY!{7m|kp9kBe0nUdL~-!g-Oo@@IMyaHMv^aCSl}PBP{`s~$Crb8^X_ zmHi5F$rYWYHFy@+vAv?!8%p5%?3ceiSz2*@agMJC(+1qIt#szsonTzQy{Y5zi(UBB zLoc}=;{}{@_HmHZISF*5Ci-e{I0MevZFG58z6OqnFp}l@5{Waa?@X(Fal^G&PDK5d z55v`UTTh#ax#7H=I)VR}4qlww5{L(3k$*dx77}E3Deg6t>5*81$yBC2Q?<;=s(*K6*bw(b{dGq7a z##O==c7?bob38e0VjWK3lyRm^v>m6ldptPhHNi|ZZCW3;V033lfrSQyuY)htAQYR8Q%>jcfGU2xazf!D%5VYr{dlH+9WhWqg~ud@!Bxc5Z= zoa8$O-1IoRhN;R7*P8b-wrq{X?=wueSX=$@n?n7>w>gn0>)N`p8%7Jb5h*@V$?!q# zVan^zFV>@e+gCw}Vf=WAd0s}^*B%dRDOqm4TZl)lS{)t@F~K8pW<1^w3wQ_|?p`}H zh6h~Q`0Sld;$huv|<&j%XW zVm7!V>s!Zb<#V{X`c3M++C^OJ_0RxKB_b_P(-i02NL2Ai_tjh8U8rX$XzSSPBs93% z*Mg&s2MsZ0*Wadgq2X1Tkt#2`(BP8t@U#9DG$hX4IXUnO4RrBv_Zx+vmM_7V_a1$W znxb}e9@$}q8fod%{kK1(W>-t~{^uJ}C$n6%k&ao%JEr&I+A;}s6>T}1enbbgr0^?^_pC>aO(BZW{Pn2mjL?XI?0(d8azwW&M+;Y% z7d(}JF@#DkUe6a~KZ+=Vd**(zv*38~yB9qM*zx10mlK6?bQB$Mv(GI?4QE89NAoDF z<04+=Ms9K-&P%s!?Wt+Osm~>ESZwyiF@3Zyujy!9>UA&5EBh|4jb-)MvQ5P8K}D=r z6mQ~w7>e(_JdQ>To*Nx2I)}!zinV!s$I+zdR!X%37oNPXQquiZ8c!^aZuIKX##2H? zLI-Pe@zlOO?)Ypfo=oL)Jd~b?C(fr{9vzm)<6<8-%L^*v@qsVQ+f!%p$fu zBDXeKdD9*|@gd@|Wfm)%;OX}o@!gEZSo6~!3sUju8y5xEX5!qt-|u|< zq3SjG9zwL^;@)2U7*CI?UU%?J!c(O`Le0dj@YJJ>7z1)En(T<%n)c8OO~?c)n@!K6 z@nb7j`k9!C{gL?XBSz2AnD!Nt=fgW_{LQYM;&a|;vZ3^r{Ow9KJweZ%A5O>98>Bd! zC6&-rr??L1zGXc1)67vLI}cASJ-OMST1)tSYEIHm0X$P_P=0LgDyH8P?+N`NiRr1b z+=OZd5#nwns`><7`+taMDv zpgRsGzQHu}SNhb?{Fv^}8nF zx!dw+br7cMh`)+l*@dS@$0ArxZNyW8M_0D!KEsnUI?E}1qIgO?^;v$|6+ELq{>Z1m z0@Kc|Q**b^#k6K?`=QSVFfH6(T;k0KG*f9t*AErLQ-#3^;_=#9gu8>%`1g-ln^@XoIB)RIJFOi78Xu? zho2xiukcZ#X(Ncw4UeDWk%xn`GVaRG9x9faGHb@R1UU3q4fTM#KIhra1S)+FhCdh8 z;7rFOuQSaIaJKmx_jFA*oNH=sE}7H-)h0W^El*iNt;xD?cNPNm#>r?K))R2Pel}Hn-9J4VWut|?X117Kv(61)E0AQ**{8Ec+qf>C!zwdsL+ zFy8iubzrRjxy-&v5RAkAvit3;M=~TZ+4h(EQKC15TfrnOa%U881(&ovB zi>n0lt*=w8FZF}Dp6=AS=twYs`G>E-yo>Go>exE4;P34==-mt!SKh{57CrzLp-JC9 z%@=@0=L?T#Lz7_1|A!M`dF5&3rr|uW41JKs`1KlCc6tO<@8Sb1et{kPx95Y^6^0m* zc{*5y{$ZAZ6V3zM%NARXz10QVkiTp%V!u3`4kW(6)b;H;x?lpNV7IG> zT@`@TZggh6(hltS{!&*li>;#(?1B}zHeFo=yLN${mu~8TJ>OsItf<*3*Mj|vhc#*& zE`WX8iCs>;wLsqdmpW^T?oSB-`Ps1h_{t)Xo5OdR{Za+W8ucfKz2bl(#;EUhSqCU* z|46dWr5a)OQd%DeC2$o_}(;J9b#E-$MCII8?11stuK z!3l)ML2w!A~v&pH8h&tIzirBxhn zrNR@S`u$-bP!lYgqVtu3TJeWUKpi)J(7rJjoY?*_2b}ik_bzp|gOduJHMqY6oUH!R zkN+prWG^@+{AKy9;|{Yz;NrkXLc|31(EJN!HRJN!HRJN*A!_%E+e z{vUedlIf4Y|H+G%%>PvMzvJCXhDxw6t?)|EPzHONVTN{&E08%_6kDGa0a>zCIvIBY zS*m>Rb+w~FK2XHkzgiN=>ayPzo9+Pl!j_+><;;L=CscdU?g^0ZJTl-))dljyH4k01 z3xFK8=Fsb{eL&9jvGH8}7RYt|3`|D@fm|_jH>z(G$OY*ix!g*D{C-W`M2rrQ6ZIW- zubl((P5J5zo5X;8@!}QxvqC@?reZmAJ=iZhY;PB@1N*T{R`G34VBc^o!}Zoyu&+DE zq!&^I_Me3{I<5wTeeFcu;8G3Pm%I9oTULX8whD8A!y&Mbmn#l**97|zUL`$AVXzO# z3>~P91N#Rl#~u%V0ehbXgQu>?!QS<4_3a&1U{5u>_u_Uo*poHy%V;Wry=B*pUnU%2 zuP;_npYjRp&yS8tYV?CW=(t;19R>SC*Y3?WQovq-ODe)Q6YMvxuH`Aa4)#pmjZI^_ z!0y)upLohF*nM>hrllSPyAeh{zlzsj*O>S5+}jGU`?yh)i8z$m#S|2&OO}9Lz-*e5 zNDJ8cn7yn__yKn2)s&`fJHSpIC&}Ht3wC>63FQAe1a=GySMPUD0;%U{(Mnz-kV=Fq zk}n7WDW4)yuFy=dr*MV37f3jk-05~7NVzZ9_sYlsDTC$R;28%XMO)XoWlI7n#9HV! z*GnKhw|i~8ZxI}4$HS9NHsW+-srGQyMTjw%DvWfL2x)|==EUy0XT4c>qzmV z0p;sOLz74jpftN*?+mI3O6Kv#Y_~k1yy)BZ<*YbR-1!6EXiEa+3Wx58Kw`5lca^#< za~3Guf-Vn9Yy$F%(DZtyejpeB?vB>C1G4YOGDR*~Agk=giHaD=>q7MaKtR90=EM~V zzns3|ew*;eu-M4xNfxlT@a9W$86^B8Q#m-S7VH*F&uRAef?cWTn@N)oVE16LTRZG2 z*j-q^CCtAH>^4W3U0i#bu;=CHP2_Jtin(-nym>W{EDa~_`Kkd)>}<#9@3mn2RrW4} z{VTA|V|~KyI0?458Cf=+*#fqj55ns=xPa}>_}~U-MzEQZJ@;^h4{S2VV90p@Y+Qns z9-T-58|jnzFUvlF^{;15-;W1_bEgGGHZ~RDgw)r1FGUiEg`-OobrayknIsIZ z(Lh}kXAzp425PtHIx$T_pcd%yK9b-EYUHUpd08Q#`VKEWZj1*iHFex4S_7zN?arY) z2$~#8+t087RO59as$?A^o-OmU=?Zb*v^`$hfH=eqP0mJq2de!GjboAffa>u^YV(Oc zpx&No4E2fvs{d}$4Z>+aeL|JoNPYp-=bG>KA7}w;n3uPOSpiUAA9COy8VBlIwka2G zZJ=hUc^dDx18T{MqTGj`K&{SC>+47o1Mzsu;+0fRo0vzPYw!aMJt4 zQ`@-}oXn&?TB}|m4pVmHT_IzHUf0}z5b+zFNFOb?|9nB%#oplO_pgNB>6+YU@(8}@ zT`agp(0oMbbr(2!x3}efj0Go;Qt1oL-r(f;XmUF7CpcLTY#MVCBIKywYSsJ(POhKT zG*{^nc6DYCKWPh2kJ;_wS9^g|xai*6)4#y!EyF^|4K8rXxqf`aj{ZOSfcO-{|CVRK zPa(ReVSoj_<1gFMWRwM;LRmT*R(Afq?X31+UPMj9+{@#~GtuBjTEomCPc-Rwi^6^* z9L=;t{!+pdMC;}|@Kw;E$r9rBxgkTvgUrp zJ&Ncwe%W^tb%=g$?*^T1+YrrFq0QH85Ye|+yzVz*LiB^3=6W{o5nWn%G2Cb`qWwPg zoXhbxq7OfPR_3UU=v136UJ@6G`fAu?Hp?J-jo#AM=XVhOQyK5|`glaom0QPhUMsG#*@v?2txl9nf zTQe}Fb_mf&S1(5%`16eB^X$))-iY3(q`PDB&$G#bCf`nvA^MLB&MSO-5&dU8spHfH zqW_vnP?tGJdeE( z7$W*$zRrVvM7s?e4_5{pM|3$BHs@sxL>pSAa%sR9(N}MMmEeC9(K|hr8MK8EeGD%2 zpNb~>x43(BGm%bfOcJWOOVlsLAMU2}5C7T1pyAVt+-(;I{s*5FS^hqD5Oq!Jv$xlH z|D_jdd=@?_@zw@)l-t_dRqR4-e7OfX_lKhngTk0nw(F=VDmH*_9*wG*`$T&U@1y3| zY6nJk7NVw0gAeyFg`#rz0j<-XCvnks7$b3A-kc5N?N7KiL~`)%tqZv1O~@A`Q7X=_ z4p37$O-F?@vb{W#Zn$`I;sM{X1IA`4+oZ%Laq81=d{;kh!Wo~lM{`RSa8aC5!qyUfi4_oXVnf0JkmZrf%--!A*u=-}YNQ!i{T} zE|86CaJwd_K;-^H-1U5QmCWJaxcA1y`quRdXmIzb8e80e2lv(fjOu5^T?PxJZKZ*@ zS?p9J&(&!hdfDt~=i!@3*CY1fm#wFc0I^qpV$)Ae+(=JA8LO1&?l zW&3WN(Py@LL(1vDRK(*SzlI&{Q2t9IPN}23%zcWot79h`VrBjzo_1zce4)>Sr_*{J zlshHy*q%p%CX7RPLfB09d{6EdKb$&r7!TNV2{wuI;DH_NoHyo=7hC3KI^}1@TQQx{beY?;1kv8zKP)fWifP-} zESgtM;+eWfft`$(FwI>_^Fcx~rj;_3XdfBJv>>Ky0h@gBj2KCvIwcoR-?~G8e@quo zsR&j{|s$z9H7^^?p39_x_pxs;hX0Z)*ix zX(FDneW@9Ijf!V1oL$T;@O_-+Ezqau}7M^jY?X4n{@YJs-Uk4IDNX ziT7~-B^P?PPb72IcO^6(*@XHsZT=x{4HvjL^Pmcsd;Z)TVDS_WHa1A(*T-?M`8egm zi%{IB4b0vGDR{(fFhZqR6^~IFISpzR@VMt_slJE^o}{?Ty-9Y%Gk%+n2ES#%Gr99G z`?<^URK<=V{-Gs2LN|E2KED+YM$WipP3^`lF(-fR9WKXl7xqrb9-qZG*>k)qwaic| zqtXwSL=}{mxIyX2Oam%D=g+` zy#?#oWwh6#-b=@1kA-hRJ?{6Eo3SM77BfmQ*tQw92V7=$&hJBw>s*>f zfzt>ZO8gx3(?_}9EU`OZJVOnq!#rubLs65OLH=#^n<)E3wQDGnKu>I(w_FW=ff5z? z-ke*_f(uXAk_1m0;_46S0o?m?as76SvHB(w4&aJ>z2kxpDym^_m`Q6#txT$J9y+U0 z+W=3wG?GE>;dlDlN?lM_`<2K%cUIJ+D{P-IO8QIN|MeE;KfmKX#qfXpJIDWTy8g@e zpQrx$lmh+(?64sr;;c!lJn+8J|2~`+2|HriZ{_?JhFu*du8Kj6Ao4HLVDH3$nEsi^ zus_F8EWA?{B$}m7ZddPw!~f!DIMUFRy*>USNHrW{H@hYcGIi`;LleJ1_Fr6p;O1i~sukEDg^5i>Yu{_k4BfqNqfw?NUhx%%X?J5Y>8*3CW?1d7dzrK&Ru zKygqD<80dk6c-c6G}h}txt98)wrM?3yyXwY$vXkX%R75lpe|6ZpNV_jngx_we0t?- zV?gooaZOZ`Bko@p87nOz%F%xoL`)O;^4%ym&jICL+xaY!<3M@(X3McvSAg=e^pS5x zGEn00WzZgr0p)#E4EIK&pLujvpB189VZUo%_BWsuDIKx7;s=!MFh%9*IH07vS?^r# zB+8BL9mpZ%yb^iJOc4dj^F#IpNpzqDjM!`XiU8$)efHMgG@#rlT@Vmo28xHJZ zP#m>#PgyRRG@$dyU_F9K#>aCt2s9b z6p6>HV?-!G5mg-4y=@5;fooNlPZIjsw6p!d$vr?}J~g;}*A2)&zq&7EngMx&y2{E7 z0eNiqw2<9%AopFonQdD@kmb3bb2*T|$chSo(ggBUyvcgS6d;ef8;d*Y5am1)7#42> zg@K_|UwIW!*h=y~OFIK){jBMHb`W9bQ%lnqI|zCm9XdKsaKu1q^b=4vP-|D(&H{yl zfyu$Cm|%U)@8~(8@USJ7am5gR(sAa|_IjZ31zV4HD--F%Te3sV2)UUx@9+8og=MYy zn(tFUVRRnn|Ck8m*(ABl{sJ6h9!#Y&;q`HV9-+R=(NN0I;9( z;?4NE8tjXQ?Aitvz&`1k3iY}>*k3=|zkA*f>@6zL+}i@Mzr@En)gWd6rtOIvg!0yE5TP2a@U?;PCQ;>5Z*d2J6 z-D|iB?1aMdGh_(vbZ$HO_AJ=(Ha&8Ycm#Hvq|`zze8G-SrDuwz}rKpxix(obMEtUUyzHsxq_ zUjrcJ=DkT--wmYTeV+|WN`T}j`>Qc(07zOgQT*5QfOJSPY(nD~kajjreuQ`+ZF+Fa zxYHU)EM8_MQYXOnr_Qedt_iT6P2QYs-2)`%-1D-x!hj^QW=i%(G>{Uh~XmofY$z;;biN_lQ*@cK-&j%dwAHW2z5!6Y0laD}Dld z-ot^!E1+Z^c6dA}00sC+9G5s4J z0eMUImSh&fk57DMt`=SavbN=B=SN*YwtV$`o`n<0&W3w9N_~NRzv&FcAqL2=_&m}V z^9g@eInWTo1mu?qH!XBCfgDuG%3D1~NO4de?U?$7BKtXO$%eauD`h zAG^|L0;I`LGSn?VtyG5xp2IRm~X3nw)@Xk03{8rv8pK}=50l$m*Ry$aeH~(OYSjH z^tMEmJpR~U_Y?%TQ&7E*uSn|eX?^E*xOC+mT|HI`y)5Xj^}8BJtIwdbMFVRYrS9` z)JF!p1oPk8iKmG9AVVcgdk*Z>9Luk*R%wq=YtwmS6-VcIRkJ4p3aVxNjQjd?lm;zRIzR&V~ zMZxNXWZa#iH(C2^J6dw&%J$2aDS$9V-VfgT)=*AW2g{u<-YNC6l)oEc{f)?sV3J#jRGY^L1;$ z!b5amOY=Bb*i3T7khX$_R@T;s-#&rGQBi;WJ`S+pYFr`PMu7Q5Y}ZDP&tP8Jo}#wa z2h1ZqmHoKA!Q6M-;Mjx*G|dE)$-Ai1NEJ-- z#z%D@T7!u{UEf%%8%zvdpPO3Q1||nnj*ludfeD-1q;VJm{9kW*Jp;y!y=z~@O@L9q9<_8l7>quCtK#5M z1*7M)d3U>?gOPI_@9H!`Fw*1i${>$}(dp`PM;i?=lB&#gIETS#&+xi)s};aV_|amS z#t0Y*fy86xI|PTdL!+d?Xs`UED|h0+=$PIn?{|bewKYP$*0W$_vu&b(T@WI z4+3PsXvLvCON8jhTDJ%?U5=PB{R^*8|Q^BrbvR-Mv5Mg37?yUxa07m}V3NjmZ0^{h8)k_C^!8pA^AW!xO7*|V6Y#5+|as6gn_+&uXVP%%*&O0!k z>QE`^N&(|#Mv31R>|nBata#v~9GGmckb1Ji08AuvzKAJ~gNd?+Z1ZEnzLx@OZL}@G z#5B)^fw3A)EVH^=WA=fG$#c%@iL1cm^4uPds6;Tiaz#J#j4zlN=O0*gn4q1$kQd8m zFmVts*0Q(;CaxpP6=iB*a%*j5(QxT zsATKrf?Z&mb=rmfz$-BAXE-M=@D|Ki0yqah5#vc*g#UEy447T^*!=Ws2vD8ab8TXC zfNF9VXPmkSRINC>g_D%s+;pk_FgQ7rHx#)GcWF0~C+Rw!PkE)VeeCy*WpsG?cRyJw@RZnB&!v5nx zHCui}m6Qgm&2Ocb^161d*!kzbdiTkhhQ>&f= z)mFPup)Ve&jsq%2gG7B^b-t;0`+)j*>mq|1!B_c3yRYQ|HHGpjL%#>8g-Yus#yI8=T2^8k~03r=C?Q0;l7ZGF~4G!AZxl>PLey zIFavghu6%5)4g+hcGCI4DeCO&usIcQ%EI$)qL$$F`Det|1`;?82becAOo7wS)m3*( z+`)O{R=52uPr!Lk{?M~8Zs3gWZC(_r0B3C_&JeXz;B2?bBF1|LoNu49KHOpo&M$Kx ziQ`;lfi|1 z-M|-v)8MkXcsNMx4MDB*lZ?FJ!lCzehYUA3FNw;3^*RL3U-`Sv2Hyea!QIwa1JT?*-%mt?x{$)L7hTwGV)0fjZE#PEgutA9Z zB{-e@y0ZbtNQGL+92w`b1xH6KWN<>qE!Jk#PzFp z5#gWLZmkQxGYeE>7Phl{{eh~q@#4dTETHbCSSfFk0V;=>&#yEoa9lk0YEN7pIF6u~ z8)|C6u_;a@ERg{m@p)aAF(Gh#$7AeZ@eCY8UJR&A?*_*kJs+wP)4M@Eg*O}_vf#`%9#i70?W-JYLM`-8wCYgDa~GYuS~ zdD$PcHG;!^{+nNIm%zc{yPE!p88}?s8)@gF1rBGvWGay=!QrUi-H6*I;IJ#*?DKgx zaM-ZHh~G;MD2u$FJBMgM89KgqcAE!KY9If+xm}Prr+=`BpnU*Jj9TN}3}K)=Y@B*k zoDCFbuTHD2#JT+H$;gho#Ccw+_v`IU5>UijpI)li1r%-}2OrmK`<-_O=K+b!x z{C3+3Ap2=qyraDV@d8YRl??r)4_La-dm3m;~E$lMETnsjf$;qoYq`-!wG*uW^ zf^{E<^SNuU!8-kmozJ}ru)Z|k}qN@osl0xSND;T?I)U^$e=<=gibEZ=z8?pN6lmM(8^ z82Pn=;%i@+?rdLdcD z2F#3Pw5nFtgPEX0l-}M9Fdeix&HS(yOk z@}Xz!9b$cN|5-457mON~&fb0W6pVs@ocT&(10y48->Q5QFxvlO?Ds=wFk{ z3^RNeroYL7p{K&p+}#X+^3h6sIndm8B4>GPe#Dix>q z>2|m(?b>iW${7s0;Ao3^J{SaM*Lq}r{GWXJ)4=e*<@uEcG+mw%slCkgmpq7OXKK$B z;)G@*vR@_|e?v5F9n-#72N0d3fRi<{649fts+a|LBRb8$f4g_N1{91R#2*?iaHlGemdFVjbExf#^+w zW2KHn`KEQ0p`m(0-lt8O^MZ)(^IH`8G$Z=UMsHlY7twlRhFd>MBbxm=-E~5vi0jn=4-yvG* z$$O@|mk{lqbkvN^0;2bfn7iiPLG->T=Gw(?h(4Fvvoq`-q8E{9sp|$2J?+i1oS-eD zr2+|fnaE!)vcF4#0nr;u8oKUVA$nHXRwi8{y;y>KY~P>X zAnP71d+LMe^m&FbUlv54xytKi{)FJrFXQk|qMh<-a=jXH9k z=vR<#WG-R9iO#V@I)vPjM?tHjLJ9eHZ*)DLBYGQt@-A7Iu%GQhMfsoqdh&CWsGka= zOO8@CZV~!zlTO;`NwmK(`9g9h(f?4EgW(HAdQ8K!CGS7K4T<8L{`_WU?}xK(kJ=EO z^=l7rUm)ft_%hS2+r=t+&Z6Np|r z&>2z8is-MGe1gLY5Zy6TN%W!%qE}SK&;PDP^zJ4d3&D?w{)sGNJD!B-Z#njST1Dtl zH~Zcu(I!L}7f{fcI*RBP%6G)e2|YV*)#r3?MD!Yg{Ri8v|NI7>71O?Fh(0fT5nfDLE&I9F&;zYUWTe}XP+DOQskb0jNN$BOBcwB56qJNXTF6mFy zH+DZuiaLVmUuCXeD3&1fIvA#yyGXQa+;`YmhG>U|^k{_z;a6;n3J!!_K8toxh6uZK zo`@@~;`sAlhPDNjITPjH)%kt+^ILkeD;Cp){Dp=uFG~r1eqCOY>>%vd(3bS-*H=U@ zU^y5}evjy?LBjr*`w(sUQ~xo^21K7cKyyk8NA$Qe0?V?5|0R?=o(d$=hfA&umn9>5 z>${j=kr{|?x43VIpc|sU;%vTXS6VEh`SDeUOcVBNJaJkM_9FV!XNP5dVmyA+ymzMV zAYq@DwqpKRqQ5Si&QXYR3;9tgTtxgkug?cgS;7u-a`spr(JLdDsxk<_u}@8yxv&e- z@9(aSk|0=o$Vougi_mj`-!9b(Vm@n)TBS|2!`sYwuy{RT=a>w()hUEtBQKnZ?nCr4 z`Inkggq<8bJcW*XBl>|5d%+ySpI=*#G3}5-^a1~rw|5CWyE?Gor495`MGNc8&f+%^p5U4OC2#pzs<`iSrkuiedGN!+867s^H zeIdpw-8$;&lfyRUQn>=wCdR|_(rES6NklJy zt4MadPWbnM#mzm$IIiyP`u^GGAO7zgfmLQCoB zMPJz3o6+;A;>rDU3+Lidd)-|h{{c?a^C~VRL3j=7TIk$-iu4S1$-OoBpl65LeZR|5 zdDo-nVX-&1oa<1%+FI`U(JWNP`dnGy5GShe_?djQAQqK&RGsyyjKyMAjlb?$u0!vB zIS4rF%b;vIDIV46WE4BOubMyXAkJ1boK8q!#1U1pk3P(hkneKMI>-AZs9v4sbiw%n zYCbx7k1co&Y7F#m?o?+-O-EAYW*C)GyT3qtwOt=--e|kzSQm;ac6EC*T@XW0>LO-j zxW3}-6LXK;blh>>J+pF-{WZA$b%#izSR-y)j&)<@6~ayL$g(?EJ;n{k=U;Mvq2cn9 z^Fp`j7jQkV_Q=gJPTU}ul6a7kf$Lkw&OdkEfNOj&AF-?4fxSlaH;&dzpx5{R% z=Qw46ySeqB2Pewl&d2-3BV~8sW-j@UJTl*LVZE`c&J8*IA;%|wvxp^5dB$9?9D4{ zN}`Vk*%Q~=NT4(_dSovPDDX^mP}lG&d?SwP;Su4BUkzP&f4(mpz1r~PX8Mg3RG8tS ze)g^{Ds#Qn#ZgN_nae>(^b*?f3%v|kdR7gtwj(*yjy2=DkF??!K{Q+yDfF2oBp26) za~u|s%0cx~Gh*X-1vfa?Z@B5`gBw*ee{>v~#|`GqNjDDo;>K5PDw6XesO8kT>i0%l zadWY%p|olUZr0-qSN{mOY4JMG^LatsC3PY1VGtkgKb|!wKhS`NOK7?S784OU6s4j*a20*4*TS>UFp%MEPJJ=?gC5R`vV2c?B2mOf!&V7RQCHh0y3@ zg>(JJPri45i?dAchTKUo!P$-cW;3Q$_=B#&{GM$JIBC9u+XWZnq890kKfcxBQpXpD zkHpk*Q4+^ZZ3e*E#;pwdxHjSxZhD4#X%~vFW7M{a`h=27XGeJ_ub>wDb?Y zZ+0waD1K6y7S52`ji?)v-58(A;5%JQm1}Z%ap4Pv`2NIQxPHrzFDD}1aO*Y^{rwK< zxVL=GLg>II)Zc!P%}wb$9+>{V=f$RAJRoVY`t=@r)aQ`4K27-vYUj$|)9t+nH3sLN zw^>R?HCH2^F^tnt^}^KY&1a@iWyw0u2ZuzlMJM}}^`7&%HJa_yqzXUkzI&}kZ$S}t zUFWt{?^=hu+YZPu2A;yL*Ehdhb09(%GXN+tF|$?SWmuBpzx|{A3oQiALnbJG+jP@#uxy z?rg`V@R;iF)W;5o&{+LVv4gJ#(O6h4&%W@}c;t%UT?;!M+&?K0Uo)}~cP{Q5PRtp^ zT|61?^^uNvaC?K%Q$V>#kYa~T?T%^B4aw0za7U|9@cN&hf&C$wRyZ5cBovRAIS^~|a6{hVJO*`3Ty^ZS=@H>LRCF(z`eyMZDnD*s!!9R*D*YSkOA97YrSez-zb-dw zuxhW>{oaFG`RBfzxNnVG%PAg|J|on6-h@3UE5)>5Zxzl5%8e7g2+f#XY?**HO#&fn{^=kEk`A{NksH7wTwy`H|2%Dcl?5lw6-v62Ty(3n!q+$;Cj zXy8la&KXxV)PFvWwRrLu>g=q3^XADz)UJFWS0(<+Kg46R*3Z@WAN=JVJluY>N$7ky zo}99|bT{@krf2N(zxkdU)35q)yLJ3|{FFZ56slL~O*k^?b z9ijPkVtP<>2Kk6Lrr-Lea=jZcee~YB6P;|Bu9|mZaEWL~k@=58)eV^b`_g#V6-PW> zuCs-=|0teZ_9+&NypP9&^*@wNti|I`YuAlDios){V-gE&t$0%7rsIn+J3Mu^W{=4; zH9WCK(;OKX;n9QopA78>@W5WNU+*?lEizW`Zg$SQRMax@!_UOm}GD7I(F-f)FX(K4f9K#?b_~ zV%bN#MI~^HX5s6E#J9L>e}b{o;tt$DT*@(%VT-#S<$qacGK*V;D;v+5`r*$z!#a3+ zRB_V&@UdML^*H9$qic`Y>hTBlZujjsQgQwK*X=c$k8n9{ZmWee31_+fl(NfZ!%>Q} z@8jGU@EO$(r|u;N9G}qYc>hfy&e>Xizj!wTu50De(o?X)O{=ACuUUOV9dDUizG^$+ z0n-O>s840_&~D-7M$LzKutYDxbrZ2ZJ)i&PzWogDFHJW47%GpuCa)~5W+=p6DCp+^ zZ5;Oo#F>>>7ULeDMUd5;N? zS6p5$<0s?sOU>>#JG$}M_g{xdFGlfb!-w;@YBwG`W^m}iVg(+n8{EpPsf$P3>u4iI z3wXpoDU@N;6Feez@tgCPhj=);?eejnYlOdvKYf(ai-#B-9DBNc<00muJk3=n@W82j zj%n_9xIbOVbhG<39^4;n{9~g&9yAa=ow4aI9^5PN=B~6p9*8g(_hLxLJ;zsKYVUo+ z9U3tq9TRM}!n6nNi>3$emX63}b1`SMIhYRtC2ig=P zA6;BBt9Db)W*9eUu}-MqoY)?!A5-*WXa@o={H3RpDXzMba|3@HhL#uS%P63G1bupW^&* z-IC13n?`21gPG+`lXW(3(OSJcanT`wKK>q8IB4Tj4NK_LmSkba zdQM!zWuw@$_C3ba%(G8;B~Z?LKOYb0di1e(*nTvs0>6J;eaV@cic>OWS6X(Y;?m|K zPuHv6xYnh?mQ!jRH-5>zo5`{4Sz$qzd0OSqmRLF-3?F0N!+=M=EZ z3`adJQ5_vN!nqDT5i5S#I8*TA-nP?`_~~v(F|J}CTzqlI{8o#77?1chb-ZW5IeqK( z-*0lmMW(tzQO5#sap~aLN5y7bW*d2e*+~Ewe>-cxMNk7LJPP``@qqw-@xVNld9wft z)RF1h=VgSeE*4+DaE*kkv@?Y|52fNy4I=WPo&LB?IB(U5HQl&zO%`*g)EeB7{m`N6 zbTF>iWizK6nuqHpo>*0WUcoJ)RrPYb-*Bg3oB(C(bKK9RE2fvGga^jXXvmmm;vt<| zPu?vJ;odMRhR`wSQ5*J+)&THD5ZL6lf@i zOLrR7rfUEyk>j_&ulpMnzc5xSEPaX6(5f+sx30K6y)44oTNM{X+vk2FeG#*)w;cUQ^d+A?$Zh-=5FO z8g5N!Xfnyw@b(up+I_2-MP~z=7| zUE^XJ&frsu#?PHTv0o_+jUPVks8}Nayq(do_WX~>SvS$(9j20{FdNi+z@xTtg3eUJGJmM%nQO>c5 zrz-LMU55Rbc9J<{^BZnVyOpmj7`YkKGU;tk=)}H&`+MpuM-@yH7*;$$jl?sqTlYHE z*W#H7;n3{eet0JJyWr08AUxg7nPv9<7oHLu;LV;7!_yni*dFgJ$1_JO9OW+$V4CJ> z%RG4=Oy|wrC0(=$(P%dw+`n}Q)83>D3HcS`8Mn&`>s^0i+EJyT=N2)T=5GFki7^q; z_U5{jnLfufCEf-;t1FnMraya(#|_htCG1e6oWoO7yIFD)zv78O%@1ewN6_fb(+ydN z&f}q*$Cfo`#))+(WO3zAA|BctXPxxuJ|6arQ{b&N#}f|fB_+3$@l@4!#@$a7@zi3~ z8LHSAo>;1NF&6K@6Q|9_KD0&R@u^jOh2PHL3EfWD(7r=> zOhsR}czhQg`?egHm6L+U*FBPamrU%-@1*MoC|}0Icefj#OykGHcdlEkFP_0e9Yq~U zOtyHa_F(U?I!iq4z1_x)gAtFcTD5zvI6EGdp8B10@-ZIiln8Vab;2VW31bK9U*nPQ z2Ln3{$MG27%yvlyDjs_w|NiJw8y;7bNiN=(i^qSM{+5VRz!TYCeM+5mcp}-tRpd-H zo}>+K<|Wp@sdFz3c1&5~$&1Dx?uggn$s=jL9I-@ybhW28+^E2lrh$g>rhIrxElS*C zy)R)$`QGNx=Xf%typUV<1)4lP9Z%368+XrW;z>FazXjh4n(7U5Z}L(=Gk2L} zE;kvDjp4@nMKtxW}9tjDl{J{)fuBnYt925Yh}?ChtKrYYE?AVB6YQH$r?{J z-H|*XtB9vXEea-XzD84LSD*B|!;PkEd`=$PfzZr%$#_y+8KOW)7P8oxRy;+Lz2$IkG^ zdJIgXQ9Fakcj3R$$P@I+LY{_4w&f{2E4)aYgQ`~4_Nt(<^i9L-_6VV=&9^_SEgeA9 z{_6+Nv3)_)YckUxPY~ldLg%}I5=IlY5&W*U57D$`BYDSaJv4Ls{14tjK{S&VachYZ zji%#RQaQbw(Tti!?ATRjM2nD=^X$t;G>c+=tJ+nF)~U^;LKR1}A7{_LyXb>xEg4Vl z;JcW%FM+*og&or#v6a6wi^TN9qiil8Nr;{%94KpUfatr}m=&A+5Pi!p(|a=a5nbXI z+Zy*EM7t*R%KWYmnn|?zCL^7NW^8MNTZg&OOmvBb*cmlMyBqE@wJsXbr;@^Vf4+z5 zE1GV6f*%n5K!&Zj>Ul)-yOvhNeGSo0f3zM|*F`hWv)Ay z--ISO@ZS0e7Bnot(l&DV4H{hg^r9=J4Rxh(RvNaPK~2{+O^4|8m3s`YWjaeE+C{;s5wM zKmWe|ef|IX^WWyXS(L{hxl*9*|yo__?^CRuEEQx;erFGmflI2o9hF$+~9sjbuCa& zhrb{36Ka<2WO%uUDv0{tmgEBD4=aYZD^9Ni#=yKZjN*7#>GieT08V5uEJBNzG zPk`amSSHoi=fH^Nl}DWTEih8rd)A~@4vcQxUu-?<3r2ZoGY;@R2cz#-O%}C}gR!JU z_s?5XVCW{x`s6 zH13H%bsv}tANh?)_rTOB+8!<40@IM%8`a$cVA^Vtu+(w}%(z#WFNA1;nc9zL@u+)X z=9jyK&rXBcr*}g?Ua*4sD(%apmm^@VAU;HUwgTopZLH2mEWy06c)K$73z+{7dK0%R z1S}9`F}k(^EWB$uBfe|`i)`DiS^Uqz;s>oWd}}dSN~JEp3pfatF3%HR@rr`whZ@JV z4C}#i?qif>&m>qKcIZAd?Fv@p6}zvoQ(*PhOz9eB7_6qE->W})2-f0aT^dt&!P+KP z)^X-JSSP(W&G7XySdUL`VF z=-7jv&LFV0ObL*=Ap^GYink5FGlT7z;;mFZ4j_rLFYmEn1Cqt$w?MN&AjOT98Rjno zX>@4b%F7AtL@I^F9xH*Jxi-t&?4Mv4OJB2vF$wI3>Fg*<0_=tB^)K!V1$&dyNk8TH zU>|+#7WE1%*blVDkYm>YSupdLpMNxvjWr)n{+#~v6IQoMv=@NfU#T=U(+(6)`7p!r za-c{^e0@H26)5MElxayV*Wh3=XJ#Z=4-R+6cXmC$2@Y{XZ(#l|IF#LH zm|1rn97fz9Pk4udBh$;)cU)88D7+!Z@y#i4JXJR-WG@7crqYl1^1lSfTcS0bwsYVZ zz49ZbB>^0Zi#sFllE85=d8E|tBT!d8?()(U2kNdnkK%;i0rjMp{z>`$KsE7Vl`4J$ z)SJ%ax35nCHCpeN@|9Mg7VqM%`SU2{fZEk6qfT&QJSXDs<_J!LaHitNI5@%VW?|F& z;AB*7bmPq)aJp`G?#Ak;;Ph&fT0T<_I2GW54wb#&)YEg|nNTY@{rad@<=zd>d=JLb z`U=2V=A^Ud!47cN<-Bx{t`5$wUmC~dmcaQ%>C?0`vEZDgtEX7X?xT z!1aUmny1xE;M&5Smee2xuHU}i*v}pSZkvWJGuzw2?TEjd4aYTbyZE4NCTk4boJ>Nc zHY9`FQ~ssn8gt;5QFnXG@pN!&djHM!B^S8O`Odui;sfs7dLnI8p;c#}nj!Y_C1BLDW`K=EJ ziTkmj&}9e|u`TQ#^I||bwnFjbBIGLvpYa}u0m_BI786$=py<)eeI;)J#XPJ`=-_^! z*rjA|8MX(C)AYv|P0>K{pftH0B8mZ!V88&1N-_{c5R{fENpj9PGeb0jzJ2aF zzf<>}y6^q*>Q>#VSEs9{zUtna)!&M}cEKKi9wDk8`Cpr>YqZ0KkL`qP7>uH5f{ z9xH66uF(hd@FQjB0gf?o&dG$l8%zsT-g;s`WW)-8`fllF- zjukmV=({p?>nO3GsV({GLkiG!viCm@RRj9*gkzDK#6BoD@A_{Cpes!0QIGQgUGhfI zS^fHu@FxHR4kw1FJv zr2=!JiDwtIPTvc(@od$5`B^|)U?(5?ssgmJCokT(EGO0}E|&~tfi@GnYph%X=qoQ> z)NQ>+_;cW=3l}k8E2S!Gf)=$nL=mvs<+;BQSxJvB@4Hf3LtIRy0OFE__{ssZg=hpO@Q z6oSXyb<1@~1n=IX8+hV@)@&UrsH+Y%Tx%u!b_38dBsOG55PHNfIo+5Q0b2MuDLJeJ zG~eeBcZ!|>n*EwmD^&xaotamO=vM?9%I-gV${A=|xSl_>w*}|fd9xOF1kRPMu};Bt z;QY#^Wo^|>a1N6~?#EKV*&&~Fl2ioFdh+~V&Jov7s_QH?On|eZN!9H?3g9eaWZY0P z3eGEWkIZ)*rKP51J0#!QfdcE5vpo*d51}hBI&G7}2Q=LE+@B*~2 z8>s7#d4?WW1?pO%luVgkpl(QhRGhjGs9dXdjYnt$b)Bf>W8K#TO4@F=&Li&g&fujP z0hM=K{GhNcP*?w0tGp!|C_f7o7)&>yjP@m$j3S`4=`8tcs{^HouURZ_6;Ph9T>e&z z1eCzkm#4%ffMTI8BcFH+C~CzDa^HJ_B5t>jsrUgXTNbQO?RW*0l}W`>b4uX!eb-%O zy998WIXU5)Z3&dsO=5E~;y~Hy@7tK|0TkqhQ0X3^oL9cat~LY|-`Nq__l-b_k>L%A zJq(nGKO0vlQ~)J0pW{jw4Jet_PV1M-^6!lQX_N-kLwrN~ zMF@LZ*iIHDSQGv^s5Ke%0jL=lI|P>GfZDyH=)2h?aOT-e*33`_XJBqE>Gua`O7~Dw zMlU!=Wm_4Adx7)o7^A(bf`~XO5N)451qGMca4y=j@?b6>IA<=Gu3$KU^W6=~RcHe^`%LBN6%ct~ z{V_;3E&!ZW=g29$iTql}tJt^b1k|20voGF60X11x^l{B5p!$SvoHYmos>Ys^OyMs; z-7;VCxZ(&<=433kg)9*MLZLB3K|sk<&k)aA0GHiT+NN&9;If~!-YKI3T+}<_AMojcnn<5oGFkk-wG}&Wa|e4Rp7E!y}{z^ zMQ~wXA}3xZ>e*QQwy{2)+ZPISuMErm>_k02_fxMu z@D0$_PQKi06$x|^?ko?jBBEa1JHX(^K&#lkJzbgbN7zj3g!eI^>8t2_4vzzE{Zw5M zwE~=*AH}UZ&;-u+rEYIl=>TW5`Mi>mDR36H4DS=m0_u3nG_(+TnY;JIC(AuV9Iik7 z^lb=G&5DXFeFK3i(QKHzVh2#ZhPP6A@_i3Qqg%)yqhe;KXxUx!5=n9H&0h7x(T1$A+hsRh7Qrm=Wdi^vfhT zM#RKzG?N2I+GgP#TSIWv(TWQRn40gMV#j{K>*sVQr=+uerVEaQ<`0&|aux)A(5^FjJwjTsn z_1517+my4cq__QG`}pmAON$@aKKeb-7C%U!g>Rqc1+a}aVtlva1=|~@1;xa9YVd5td*p?EU%^(t z<;i!d8`!EixbiBn@tdF8C3*vF ztW*@DtSi7qr+AHc`a7_ZH4e;JiUONWz5P3WeFp2_D!YZ;zJqnQ_qIu@B3RdTlmFz8 zfOSSu%B2+#z&iAJv2DsTuqK^7=D|G)){5LF`(hq~^@{Sbb1tD^)u^EOS>hR3z0mj^ z_$e5yJR~k0TKfsCPUuOk@mLF18=8eCcgTR{Sild7(tTi=zkdP3AA)7r*!vSf>0oIQ z{iOGzKUm7WanLlo1eP4$K%&|JxlO#%I_@rz(=Ku;@3;o!OI2U?WbFd-d9$tVA(KEp z_Ne#VuAe~OG5RUkaW{}xwAJUo$^}x-;^w?l!a#a&L)z}33#5AoZp}$}0?B*oeoL4Y zkk0uo?kC3pN$nD&ss196_I$6>YkveJp~5rGa#27M_+rF*Ujd}GeXm4MS_6si0WVKp z5sRC7gB(rz>E}NAC>9W~q|LF%nioQC%!V zbx;b*aZV+!0x z)~p3`vbANnWE7Bd#CSzQp98t^2S>ryLLiSvhY93o0QvWEZ12d53g$kOR2n} zng=OhsU9|+8X*dn7tY^sN-hLTpN#V@kMDwI$h50v{Z_EN{pvPn_!qE@s2dbL>jjn% zx4WHs_Zlo;dg_iI`vI0XRGwN&tXDW+oqBl{EWdPxKH(4m%bIsF`*@APvf5WEc_yEj zn^%c@Ie}%>rsvDAihyN_o&=i*A-71s6qrdW~8ISZ3`mSl1{9mU*W-mKHjR z`||cUb`tB=Y>gUceZjJGUL{50I9SdPK9%Qq1XgP$ZrgLEfz>v@)9T`7U?taF?BAmV zRvO-%XFR3B>TI*gt<^4I<$7`TVBR8F-4xMrK41V=Y1;n$6@6gUCZ==fhYwhe&a01S=$5>o=R#QNNY(zz*Mm_+ZIvZn*%3|Gdf zH-o@1Fk~d`P9ZSWOM!9-7-l9K)dTnCBI`}|2Z8%$lT!!h ztiXNW+1F1u-vf6I6+!g#F}PcSi5Ia8L0F*uBXc-1A4A zO-YBqy*_BVQfDQ&536qbq#y?F-`+gU+PVuoR=KZ=c=Q-NcsG{RkBx%I#xnVVkZ0hr zp>y*w+fMKh_~~HW^A0@NXN`=H%!2!O2k!SbBf))4KR&QCAKV*#I@C6(f_wJP>p_3} zCtb0RRr?g>!9BnrR-mZ_+{vb4*YBi)yE-LAug(?RcO7uL@zNgL*N7zv$tj(|lhCD0agPIC3HtNeBPD}yUKe;hA(}LhSF5XxaR|T$Za)~iw;^10B zzc0A`61cvt=!f{WF9kL3C^aM3H<{_WuQI-=&cVE|UDok}3#X1iUW(Qv3ogZ2oxx z3lE9DE>>vA$Znw5e;!}=hzazM7uvK#G=ZMDCT!U29;iH0Q<0+#6eYiuyElV*{#EZ8Y$Zw$XD33jWFJAUE@mKz{ z1pZ)yh7R6AIbg$1(Nf@*2J7OQLAyteV14!FlEl6rV6And`8{p`>vgz}-k!;(W1HP;=FA0Lxvk z=8~#2fjr5E7eCqsj{5X(q z-EN(k{RSixcAqzb{Xmi&Nikc@0MfFLpMQF80E>^I77<^oz~Yj}tHA^=usFPQ;q>lN zuweW1bIY(9n74SWez!LZ%;US1l*jggIZb-Rc#a0;k#F*#3^4zY$B7V$i#m2eHD)oIbVIh;>%%6GK!3F$)i9+G_V%2j*Ew^tTIB(1%Za!E1VFkR2<@CA`we-hyh@urVii>{6X_=Q-R}9^eV*X2mpAie=s3aq zW&h#dfm%C*LV%4?Z z+vxumv2uSbNdC?mJ1Eq-BL(MlmJ&l<69lu!&ZzufQ ze=TIyJH#wXrrMd9ps75Hm9c*hn))**q4(fDVlh`l@s(*2`ez*RnIU-hu%(T&2t8O! zU#gC7A@CyfgjFC`=DO35S40zjkQ{=1B3_2h&hF?U>}SEq6L@lpu=8Q*mLHXf^`ukd zL)m(Qw^9r}>u>M%^AGZ@I*3@&37KoW9T2Ol~nOVomRhsC&vntlvj!R^gR|{W`1an)VQKUS+(OPDHE# z{$2b}2>UnipZ|SS6tUipaUNSm@NKVg|5(0+Sj~ek1}O^!-|Q_rG#?<=frvZ4OSXu` z;1KRSybiIhx^$&d3BOkEXWdgI__I!EN2s`@5V{5bN9A`1y}Dgq*;>^(BP7%nI-M8V|$@DiN@^B=mXto;Svu zi03-h@ex@fj>e4Nyt5|eqw@wOC1pfB&5QWO6Y*3PZkzO$u+!|az#c9!#F{ebSz0eh z*n{HdB|z{^llgYxTsvY73-KVmsep4zA#g{CUquYb7HjF{<$y82O@5Nqb@ce0^AVvW?vrZ3-3 z=wWwmA59gpu73y`7Qy#TN5cxroJD+pi88Lt7OyAPBM65(%+VHLu zh_zSh`qDnaj^)Ye=KDAadz~wJa8mnk->XnkNg?bpDxX~2`!_#4;F7s#M8p9@MEC>2 zo0;O-ZA$p7%kQM5K2axfxKC6)BCaQ2RjBbJ_Nx~%ELpys(C?0F{XT+E3hU6XBq9%U zx1M(~Hz)k=R4lUJjPPGwCBr@)v7#L0pWhuutT;2Hwg&_s=264bMXrQ?0Y9S;UMBL( zyiBQ=$S-}>WwA?yya3*ptz|YudM1v z!kVsy5OqgG#q~T&{$WZqPAo#Yh zshuM7H&##g<|1+301hTg5cww1!KRf=;C)ZYy%|J(>J2|!rcLmDcjvI9i#MTf!~lgv z#OoLHrZLyQeGh=Uqk#y)ucBqy5k&A~*(&jw67{j3L)*!dkn`Dp)q6=IE^3rXlD7yw z6K;D)f72lBIQ=EEIE~P2>e)&!WyI2HNvHe2{0sm36&3pUb!_ENjplz~6RN&R3p-*q zi|WG$#oh$>{n(u%r(9`!*@$6Tx^9L}h^qQvKbqB71>w&fYDmckV}9&7h*O~bXu z^tqI}aJjguX0#UI z+L)&&isH}X#=V+yPuAJtIw9LL$}K-}!}InQdHx@`S-0Y(h7b>~$Id%WJ^zSngD8FZ z>-uqRB8#iTSOGV23y<^h9mVwyqxTmMSK)fs@#8YjQczpJhKEf49PU^ZbpN3FTHGS# zF@62Ob6kC4xo$!6QJm(o{ABLiuQ)aC3;*-aFHm-Jc0kUluc&bEnB$%?3;crDFl#ev zz*#m|$G=IWq9VMfY2@Y))L1^Mc0MKy)&COR^U%^B3~uM< zY~P3*Ukp#L82f;`cnX!n+-LAW4I--sZbAdUM_Fc79Jrq=Ad5{(6Zc-~R*mdc$Gx5k zNA3%r$DJaxnaU}OxMqDoZEngNTr=LRcHJ`+H$OXcPWMYIZhhH#Eb8zW?({n3FX^)j zcmDBM6O?)!x5(Tq({As^E%CRM9^K}_t@86OnN6E6!|qqN)bRV1 zo3E9wWZ>M@3p&Z^llZOvds^;gC0uB-S=}sM7#BO~^v+2i#08;>H`C&Z@tfd%6FMq# z__^_ZnYk8DbpNinjP%YzeC;E*7{Biy{8&X|gO{u+GViT4ZVnv7HafT1Bd=4@m3y1% z7B~P$E?3`tarPyyJRCgWW&0i1rZZn{wT#5IzZ-mC`wL7ZY75(*9(pgmf^37a)dp)YXn4aJG zS`p>OE%&OGWa6qacB_NibWjT~)tg6r26xGbm7(gApIU2V!eGZq>u=ovW2x{XItzRHRQXyVal!JEzU zC(x*(WSu~*E*{*K%ddm7aZku<0r@KvxJw{$%yGIBb?>Nn@L8%D_YD0^4hWk--F`;D znrpt}&O@h`ifz?!>!CJ|qH70nZS@2zD)%%lSKsj`mBSgAbFg^Uc-iC6$zcW^@!K&z zoO!PK$YJc48UAy#{4$iCdCQRd+9gyRxTNbmE{6(y<@n|=N~1z%`NQ4&9-<<3M-HB> z?dXH_=<@i~LG)fcb@0%N_o#g3(_}7xdDQr^$O*Ck!7UM+FZH}#fm>EA9-cJt!;Rr7 z*8TCmxO8`o&eMygxGI9894B!Imuu=}z4lVV`DbN@YVY!+%myC(AvXxswk(q;D}={{RjyldpMuxrB`3rHc ztZ+oagVKLuUr6HC>VX*CMLwSY{YN~m9>1=7HJy&H49f|XAizEB+2IX)I`QbbA02p` zH)dr&v9yo4hgoW@>@yBlc)aL=>m`fbc)Wlayl|u%Po^)OJl8afr}DRvWX2QlRC4*Y zbLLm@RDW@k@XH*`^!)WyDeDDh>fMNCZBfL`yV5q2WGBq>dUfFRJsHfp!w$KRUtv~< z>G9~1=a@CZufVq;fLRl{?6xE|V*SAC#Qj~E_2a0RZ($^674B?uUhjpc%6sL_lAhy< zH;VdNb4z$kl1;NrUmcIu$4&jbvKNo*$UfEF{REHa)b}S@a^i85t)&m0QI9ZRH0aZg9c7MN%_}H=kP#3pM?-Qf(EI-_rGWw z!Gp%dIb1vY@o-q)AE%Xhc+Bw7cBTVi-w)Dr_l*KE^MymfhZ0ZBs=ik`npA~Zt_Aah z6E1kFM{Jlm!G|ZiQugYkXW)sJ2C=Dk9e7|VG*sb_JnrF5}{J4`BZ;Kxrr#=;v( z_``tDTxvuhF3d31-Bcrhimq)_y|8N$XS@?}NU{$^PBbyIi^th;d~)j)pBOLRHZ!&u z%EpF2JrEH&{pLIVjOI=Eal{knT~T@Vz)IZatgN9i&4D^6w|1HoZ^Qip1#e+XD<1S& zuGdxe8xMYRyJ{t_hI{+1N8pq=?zpsK#<5Kuw;kJ;AqQ-@L#iva(Vq=>1}^WnAqV63 zREvr@-frCal76Ubt{1m#wb#_fM!31tRjvEHBknbRNdG?NjR(s{lz*9iCh8r>hW6XZ zc;uP(yW19Nc+})w>Ewn(cvNE5z`DUzc+|10q(ajRkET5PQ+h%OkGgv7j?u}(BTb`& zl`oy~@Wx{+Ob)!n!>)O{f8O1}BhO3ypWOJ3NA}-6vhscj9=df>W} z^pJofLb%w#!CO*Y3T02W9<{pm2?wxmOwer##2K$6zpm{W#U=D_B0AhjxK7k(&&w4u zxaFPPaPCAs?!Wo+t4ypi?y&4^xVOa_)i-A)1>~>BrM7b>GY!1B%3HXhQGFCQtlH+- z79fi2Eb~4DpC-=9Vsadb(+GdcFwpW!kHDo(f<@L>PT?=*{Z>Y4leh)0#!OYK;KqpV zh#mSesIs`{upHo7%czX7t-1;4oCxYlWZNHQ3{dX!j z`_g^E-d}$4?a$TV`Jgkq?{Oc(+hAd9+H_bo0IS`i>9{zbq(+}4f z4e0UvsNkBm;d6;6s&HX;ak@6;CQdMm{57g0i@o+U)MOI_aj0=%eBtLn{OsZMnHA-1 zIIU7i{q$%r&XBiVsiDe-Gtbs|nH!DboQ*q#gE#|m(xl?SJ9VG%yW7Vu>?g2@KAjR$ zRE0}l(79JN)#ADu9<|+?=W+9o?W&8)M!0GB%Wc&~PjQ34=(grNyKs}x*ZMOh5xBvK zmG|i6YusqB6~~!eh}#@wXAG`f!@UnL*zaGq1^0(Vt)*@g!$ZB)WZkArJX~}8%w6x9 zzYtBB3e~b5Wc&kN(O9zSLxVjx(fHJclP9(~A?B(_Zf2|IFmwN7^`aarJo>)sdDloB z?i)^UZ;MUzP3 z&xmt>PN<+p;W8@Zvi9`gGDaCmOoqSaM^yINeox3PTU6E2^khLj9Th66BrV8wApC+S z=jerTl#-$M%>`1o5cl+ZbqxybVgdDISmnQgy@Vs!2?v(C4o*8>5qwO`v& zJw7?;J`s-EKcpv}6+4Q00=PGQ3=c*|3Za{KjU!mRG^&5e59^r}8C z-b5p=R=pZL!^HW%@kH#Zel*sxW~a%-x_HghZgRCO<=^J(g%RLu?^-V>ucv-&4&{*FhutH~wr% z4Z(x+kC*Lci$H^~lqM^V8>4~X6Hl@)`l0^SPqwbKw?QqMVx4g}HBmip+k}+a2zsq} z;>BBa3eG>hL3~>_Q8&w)tY6T>aL1O6`KR4;xNG4FHR!z`9)7<=C#WtAj~hg`%N8;* z^OV}C0*)a1g_El?&Z}Z(=;->LVQ!fDIbGxuKPPeC6<;y;>_AgJH#ikI*5k?F30aRn zb>K<8j>}{E33zgTgLrk&DLg67K_Z`a#}jth`&|;h;7RWXAvy|Q@D%Hy9gnRaW(Gfc zQqjH~vz(vodREbaS;Cx~X`2Z7Q8g=T1oaWqxYXz{&ne9ORDZc4$qq5!U9=w3yoH$R z>Q!wV?gZvPsXCa5nLDVOnO}qP)cd5~C$8)8gs698Rg)1OdPVZ+xvvoq>ec_wYo(!{}+LcX+}r(CcV;H6G2nIP_dm6%V&K zw{I(JN5ezc+ZuNU;h_s*4i<}Cc<}nMkXj)fJQ&p(sw6mt2UmtXP37;!14QN4;u67w zC~~dX)$@2r?7=~f)_6S3jJ)Ug=_VeIw2`@vE7M_&z`c=x*tPdtmLBL(^4iC=kZ_j($mNp3rL z$;DYbp>MTCj$Ig!mB>7_923RkD)l{M{_J?X$y(G-I2Mn;{go#hC5tDdoE9T=8}LM@ z?YZ8&8hBEHT->Ft)mUQ9G;_zJZ&*P+QRV{V;w z$Iy7R$CBC(KQz&u5GC4Fh^7M1Hox9WM=Y`T0!M;#5Hn}_&OGl*G*!SZac!qH8W&x^ ze68+TG?o}07}a8iM~&aGKR4KcN0Yc6(rl#A=neIC+gFm%sB$|;-W!BR4*0J9xZyV% zqW|H!)?9)H_s4d`{At32x0r(B7nRW{r8$dz=WR5RQTs`~#dyU45xNmrH$D)bOr1uvuh@q(?MhQ>Od7&x0HOd}Ya&dOcUyC(QIwvWdMLiCF8_4IHi}^b7vc<#&G+vEn$N%kmyU zEKaJ2;tFHLqzvwKo@zkMigC?l{@jQaFjPFqK8#pqQDf0-HS3qFjIAts3Wr%vuJvbaYIDDZLw{%@(;wU12nM!yW^PkAwVs*I0&&mJM0f; ze}XjX*)!<+FfBqHKKi}uCXZs(27y93ye}Dde`2270^w+cfW1v6dD(n_5i?D9aW`^^S zJ#1pd=5gOu1mXIwsJHinL2NMhTH;ZE*j*#OMF)IfZ_<6%mbj;|kMiAOYo#p6J>}5; z>6;Au@74A{;V1;f+nzS`Q5Oh7I> zHq?A<1<2lHWGb5oN75$uCXQT#qc0XZeyco&V==kqw%VPbe(g$ublxX89+L2?sV5jt zT=J+KyEO7c{N54+o522d#50>n1qXfVM$%GQ*@5 zbTp(d9NSGFx2QYWyFJ0sLp1)q0v{MoTvy=BRD{#E$>JtI72x#18h6j*7!T9K-y4lO_V4TPI0UYDN zL=kbUthNM`w-Th#iN>!cdmVtA-dg1!FayS>)tZ{Ux9n7}Mt{stM z2eX*+&%ciu!}$#t+3zt!;r#6_-Rn|f;ljFv2cLen!G&A&_{>LVzL~Fb|r| zIC*sv%ztN;ul~*k7M>olrlz!TN}RK};UHKQrJ4^{o&~Evf9{oW zXn-}U`xT-23aovG4|s#EWzDru+Dx#0@_31M#R+Wt)OfGl=mI;DtfL(PB(Njh z;ME#g4R&!o2MvXefL)Je(U4*q*l(NTjd`OA_UF7T*SFMy{aup($Ups($kB@P!(n@Z1*R-^>OST^q-% zn_GbrW@F2YBmt$`##UbV1yFg65ARmL1618phkzvk)UX`~;|}`)wOXWp%kg<|<`I;N zooNPVoumO1JOs|yTj$xt6T!Lc&uJM+51_3sU+d4)0kl(ZSIK##0_~b6q zXyCj6boLfQixp>pF812S^xzAiA7%AL9nS)~gXtMQU?0digkt3*MMI5 zeqtu}70~+wUI_yaxGdKm^;rr57t!6Nq2>YLa%8^nv#~6=kk3eNTh#zAL6Sv1F;BoH z-VeEGn}N$Gmml$aM8T!k@dj(-MsQ^_5teyU53VALaosn$!4=8BdL{Gx5j8-_{+3x7x4$iqxDvk3L3!p$Z@35+zl9=NuFAge85;tv@(i$1a8|jEIPJ) z1UJ=VJk*1~!OcvV)K}*RZoXW@!I{K7cBaSUv;w&0-q!bhK?1jSMVA=MGI0C3@wHGo z3*5J)P0o;H!2LiPYS_*IchiKf?VeKL?(J)7=Cu#p@9VD9mS%!`PRVygLO=IbhwetX zYH!7{EAvCZ!;_yZASVVMci9q#9e07pyOBp> zm4Cpa`Sa0?vv0uT>+3BbE)1R>OaC7K9{(Qy9{(Qy9{;oB|Lb@B|4X(jviPWC@f38yQ~Q^uKCXZ$xB{Vz|h_*cqn{;t@6 S`1@ac|GxlpQvIJT Date: Tue, 14 Jun 2022 10:54:12 +0100 Subject: [PATCH 12/18] added dk to descr --- DESCRIPTION | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3e9be63..bb4b724 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,7 +7,9 @@ Authors@R: c( person("Gregory", "Jefferis", email= "jefferis@gmail.com", role = c("aut"), comment = c(ORCID = "0000-0002-0587-9355")), person("James", "Manton", email="ajd.manton@googlemail.com", - role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9260-3156")) + role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9260-3156")), + person("Dominik", "Krzeminski", role = c("ctb"), + comment = c(ORCID = "0000-0003-4568-0583")) ) Description: Extends package 'nat' (NeuroAnatomy Toolbox) by providing a collection of NBLAST-related functions for neuronal morphology comparison (Costa et al. (2016) ). From 0f32a09815cebceb4617849df47e8a640378cef0 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Fri, 17 Jun 2022 14:18:38 +0100 Subject: [PATCH 13/18] another idea for topo scaling --- R/neuriteblast.r | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index 1060d54..59a631f 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -503,7 +503,8 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL sim_scores <- sapply(names(alphas1), compare_feature, USE.NAMES = TRUE, simplify = TRUE) # Aggregate feature scores - scalefac = apply(sim_scores, 1, prod) + scalefac = apply(sim_scores, 1, mean) + scalefac = scales::rescale(scalefac) } else{ warning("Unknown alpha format") scalefac = rep(1, length(dps)) From f8a6b316265d476a6a93f9c12bd1c72dca776866 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 23 Jun 2022 11:48:03 +0100 Subject: [PATCH 14/18] new fixes to NBLAST --- R/neuriteblast.r | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index 59a631f..b7b8ac4 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -335,7 +335,15 @@ NeuriteBlast <- function(query, target, targetBinds=NULL, normalised=FALSE, #' @importFrom stats dnorm -WeightedNNBasedLinesetDistFun<-function(nndists,dotproducts,sd,...){ +WeightedNNBasedLinesetDistFun<-function(nndists,dotproducts,sd, sd_dp=0.5,...){ + summaryfun=function(x) 1-mean(sqrt(x),na.rm=T) + sapply(sd,function(sd) summaryfun( + (dnorm(nndists,sd=sd)/dnorm(0,sd=sd)) * (dnorm(1-dotproducts,sd=sd_dp)/dnorm(0,sd=sd_dp)) + )) +} + +#' @importFrom stats dnorm +WeightedNNBasedLinesetDistFun.ORIG<-function(nndists,dotproducts,sd,...){ summaryfun=function(x) 1-mean(sqrt(x),na.rm=T) sapply(sd,function(sd) summaryfun(dnorm(nndists,sd=sd)*dotproducts/dnorm(0,sd=sd))) } @@ -481,7 +489,9 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL idxArray=idxArray[!targetdupes,,drop=FALSE] nntarget$nn.dists=nntarget$nn.dists[!targetdupes] } - dps=abs(dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],])) + dps=dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) + #dps=abs(dps) + dps[dps<0] = 0 if(!is.null(alphas1) && !is.null(alphas2)){ # for perfectly aligned points, alpha = 1, at worst alpha = 0 # sqrt seems reasonable since if alpha1=alpha2=0.5 then @@ -504,18 +514,31 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL USE.NAMES = TRUE, simplify = TRUE) # Aggregate feature scores scalefac = apply(sim_scores, 1, mean) - scalefac = scales::rescale(scalefac) + #scalefac = scales::rescale(scalefac) } else{ warning("Unknown alpha format") scalefac = rep(1, length(dps)) } + # dps = rep(1, length(dps)) dps=dps*scalefac } } + #if (!all(dps == 1)) + # dps <- runif(length(dps)) + dists <- disc_distance(nntarget$nn.idx, nntarget$nn.dists) - NNDistFun(as.vector(nntarget$nn.dists),dps,...) + NNDistFun(dists,dps,...) } +disc_distance <- function(nn_idx, nn_dist) { + tbl = table(nn_idx) + for (idx in names(tbl)) { + cnt = as.numeric(tbl[idx]) + idx = as.integer(idx) + nn_dist[nn_idx == idx,] = cnt * nn_dist[nn_idx == idx,] + } + nn_dist +} dotprod=function(a,b){ # expects 2 matrices with n cols each From 241c9de05748290fab3ce6e5cfe905064ffd533c Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 7 Jul 2022 15:01:54 +0100 Subject: [PATCH 15/18] distance discounting by log --- R/neuriteblast.r | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index b7b8ac4..ed4e575 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -526,7 +526,7 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL #if (!all(dps == 1)) # dps <- runif(length(dps)) dists <- disc_distance(nntarget$nn.idx, nntarget$nn.dists) - + #dists <- nntarget$nn.dists NNDistFun(dists,dps,...) } @@ -535,6 +535,7 @@ disc_distance <- function(nn_idx, nn_dist) { for (idx in names(tbl)) { cnt = as.numeric(tbl[idx]) idx = as.integer(idx) + cnt = 1 + log(cnt) nn_dist[nn_idx == idx,] = cnt * nn_dist[nn_idx == idx,] } nn_dist From 637e071e9da2d592b425bc3e14d1fd672ee8eeb9 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Sat, 30 Jul 2022 13:07:14 +0100 Subject: [PATCH 16/18] cleaning code from experimental --- R/neuriteblast.r | 19 +------------------ 1 file changed, 1 insertion(+), 18 deletions(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index ed4e575..b658924 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -490,7 +490,6 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL nntarget$nn.dists=nntarget$nn.dists[!targetdupes] } dps=dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) - #dps=abs(dps) dps[dps<0] = 0 if(!is.null(alphas1) && !is.null(alphas2)){ # for perfectly aligned points, alpha = 1, at worst alpha = 0 @@ -514,33 +513,17 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL USE.NAMES = TRUE, simplify = TRUE) # Aggregate feature scores scalefac = apply(sim_scores, 1, mean) - #scalefac = scales::rescale(scalefac) } else{ warning("Unknown alpha format") scalefac = rep(1, length(dps)) } - # dps = rep(1, length(dps)) dps=dps*scalefac } } - #if (!all(dps == 1)) - # dps <- runif(length(dps)) - dists <- disc_distance(nntarget$nn.idx, nntarget$nn.dists) - #dists <- nntarget$nn.dists + dists <- nntarget$nn.dists NNDistFun(dists,dps,...) } -disc_distance <- function(nn_idx, nn_dist) { - tbl = table(nn_idx) - for (idx in names(tbl)) { - cnt = as.numeric(tbl[idx]) - idx = as.integer(idx) - cnt = 1 + log(cnt) - nn_dist[nn_idx == idx,] = cnt * nn_dist[nn_idx == idx,] - } - nn_dist -} - dotprod=function(a,b){ # expects 2 matrices with n cols each c=a*b From 930c301bed84fc72d203fea2b864d2e75d15faa6 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 12 Jan 2023 16:30:38 +0000 Subject: [PATCH 17/18] cleaning of topoNBLAST --- R/neuriteblast.r | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/neuriteblast.r b/R/neuriteblast.r index b658924..a3c040f 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -489,8 +489,8 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL idxArray=idxArray[!targetdupes,,drop=FALSE] nntarget$nn.dists=nntarget$nn.dists[!targetdupes] } - dps=dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) - dps[dps<0] = 0 + dps = dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) + dps[dps < 0] = 0 if(!is.null(alphas1) && !is.null(alphas2)){ # for perfectly aligned points, alpha = 1, at worst alpha = 0 # sqrt seems reasonable since if alpha1=alpha2=0.5 then From 5c884dde9391185b1d670e05115d4d46981fab17 Mon Sep 17 00:00:00 2001 From: Dominik Krzeminski Date: Thu, 12 Jan 2023 16:53:02 +0000 Subject: [PATCH 18/18] cleaning topo NBLAST --- DESCRIPTION | 2 +- R/neuriteblast.r | 25 ++++++++++--------------- man/nblast.Rd | 4 ++-- 3 files changed, 13 insertions(+), 18 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index bb4b724..7ce476a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -35,6 +35,6 @@ Suggests: License: GPL-3 LazyData: yes VignetteBuilder: knitr -RoxygenNote: 7.1.1 +RoxygenNote: 7.2.1 Language: en-GB Encoding: UTF-8 diff --git a/R/neuriteblast.r b/R/neuriteblast.r index a3c040f..1a91137 100644 --- a/R/neuriteblast.r +++ b/R/neuriteblast.r @@ -46,8 +46,8 @@ #' @param UseAlpha whether to weight the similarity score for each matched #' segment to emphasise long range neurites rather then arbours (default: #' FALSE, see \bold{\code{UseAlpha}} section for details). -#' @param UseTopo whether to weight the similarity score for each matched -#' segment to relative topology of the neurons (default: FALSE). +#' @param UseTopo whether to use directed dotprops vectors (pointing towards +#' soma) and features in neuron nodes (default: FALSE). #' @param OmitFailures Whether to omit neurons for which \code{FUN} gives an #' error. The default value (\code{NA}) will result in \code{nblast} stopping #' with an error message the moment there is an error. For other values, see @@ -335,15 +335,7 @@ NeuriteBlast <- function(query, target, targetBinds=NULL, normalised=FALSE, #' @importFrom stats dnorm -WeightedNNBasedLinesetDistFun<-function(nndists,dotproducts,sd, sd_dp=0.5,...){ - summaryfun=function(x) 1-mean(sqrt(x),na.rm=T) - sapply(sd,function(sd) summaryfun( - (dnorm(nndists,sd=sd)/dnorm(0,sd=sd)) * (dnorm(1-dotproducts,sd=sd_dp)/dnorm(0,sd=sd_dp)) - )) -} - -#' @importFrom stats dnorm -WeightedNNBasedLinesetDistFun.ORIG<-function(nndists,dotproducts,sd,...){ +WeightedNNBasedLinesetDistFun<-function(nndists,dotproducts,sd,...){ summaryfun=function(x) 1-mean(sqrt(x),na.rm=T) sapply(sd,function(sd) summaryfun(dnorm(nndists,sd=sd)*dotproducts/dnorm(0,sd=sd))) } @@ -489,8 +481,12 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL idxArray=idxArray[!targetdupes,,drop=FALSE] nntarget$nn.dists=nntarget$nn.dists[!targetdupes] } - dps = dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) - dps[dps < 0] = 0 + if (!is.list(alphas1)) { + dps = abs(dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],])) + } else { + dps = dotprod(dvs1[idxArray[,1],],dvs2[idxArray[,2],]) + dps[dps < 0] = 0 + } if(!is.null(alphas1) && !is.null(alphas2)){ # for perfectly aligned points, alpha = 1, at worst alpha = 0 # sqrt seems reasonable since if alpha1=alpha2=0.5 then @@ -520,8 +516,7 @@ WeightedNNBasedLinesetMatching.default<-function(target,query,dvs1=NULL,dvs2=NUL dps=dps*scalefac } } - dists <- nntarget$nn.dists - NNDistFun(dists,dps,...) + NNDistFun(as.vector(nntarget$nn.dists),dps,...) } dotprod=function(a,b){ diff --git a/man/nblast.Rd b/man/nblast.Rd index f3fb1f6..0fcfd46 100644 --- a/man/nblast.Rd +++ b/man/nblast.Rd @@ -39,8 +39,8 @@ query} segment to emphasise long range neurites rather then arbours (default: FALSE, see \bold{\code{UseAlpha}} section for details).} -\item{UseTopo}{whether to weight the similarity score for each matched -segment to relative topology of the neurons (default: FALSE).} +\item{UseTopo}{whether to use directed dotprops vectors (pointing towards +soma) and features in neuron nodes (default: FALSE).} \item{OmitFailures}{Whether to omit neurons for which \code{FUN} gives an error. The default value (\code{NA}) will result in \code{nblast} stopping