From f1d8f758f08299947ed6f9052dde183d6c760131 Mon Sep 17 00:00:00 2001 From: schochastics Date: Fri, 15 Dec 2023 21:36:23 +0100 Subject: [PATCH] cleaning up and new lintr --- CRAN-SUBMISSION | 3 + DESCRIPTION | 4 +- NEWS.md | 4 + R/bundle_force.R | 56 ++-- R/bundle_hammer.R | 76 +++-- R/bundle_path.R | 158 +++++----- R/bundle_stub.R | 372 +++++++++++----------- R/convert_edges.R | 60 ++-- R/data.R | 17 - R/flow_tnss.R | 624 ++++++++++++++++++------------------- R/metro_multicriteria.R | 78 +++-- cran-comments.md | 7 +- data/cali2010.rda | Bin 1214 -> 1229 bytes data/metro_berlin.rda | Bin 11357 -> 11366 bytes data/us_flights.rda | Bin 27473 -> 27502 bytes man/edge_bundle_force.Rd | 11 +- man/edge_bundle_path.Rd | 9 +- man/edge_bundle_stub.Rd | 26 +- man/edgebundle-package.Rd | 2 +- man/metro_multicriteria.Rd | 4 +- man/tnss_dummies.Rd | 4 +- man/tnss_smooth.Rd | 8 +- man/tnss_tree.Rd | 6 +- 23 files changed, 769 insertions(+), 760 deletions(-) create mode 100644 CRAN-SUBMISSION diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION new file mode 100644 index 0000000..958d534 --- /dev/null +++ b/CRAN-SUBMISSION @@ -0,0 +1,3 @@ +Version: 0.4.1 +Date: 2022-11-22 06:42:52 UTC +SHA: 9196928bc4eda31af1e25c30dc6e71912b87cf58 diff --git a/DESCRIPTION b/DESCRIPTION index 843e41c..763e8c1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: edgebundle Title: Algorithms for Bundling Edges in Networks and Visualizing Flow and Metro Maps -Version: 0.4.1 +Version: 0.4.2 Authors@R: person(given = "David", family = "Schoch", @@ -19,7 +19,7 @@ Config/testthat/edition: 2 Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 LinkingTo: Rcpp Imports: diff --git a/NEWS.md b/NEWS.md index a67deee..f9090cb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# edgebundle 0.4.2 + +* clean up codebase + # edgebundle 0.4.1 * added package level docs diff --git a/R/bundle_force.R b/R/bundle_force.R index 8268546..9ea25ec 100644 --- a/R/bundle_force.R +++ b/R/bundle_force.R @@ -20,39 +20,49 @@ #' Holten, Danny, and Jarke J. Van Wijk. "Force-Directed Edge Bundling for Graph Visualization." Computer Graphics Forum (Blackwell Publishing Ltd) 28, no. 3 (2009): 983-990. #' @examples #' library(igraph) -#' g <- graph_from_edgelist(matrix(c(1,12,2,11,3,10,4,9,5,8,6,7),ncol = 2,byrow = TRUE),FALSE) -#' xy <- cbind(c(rep(0,6),rep(1,6)),c(1:6,1:6)) -#' edge_bundle_force(g,xy) +#' g <- graph_from_edgelist( +#' matrix(c( +#' 1, 12, 2, 11, 3, 10, +#' 4, 9, 5, 8, 6, 7 +#' ), ncol = 2, byrow = TRUE), FALSE +#' ) +#' xy <- cbind(c(rep(0, 6), rep(1, 6)), c(1:6, 1:6)) +#' edge_bundle_force(g, xy) #' @export -edge_bundle_force <- function(object,xy, K=1,C=6,P=1,S=0.04, - P_rate=2,I=50,I_rate=2/3, +edge_bundle_force <- function(object, xy, K = 1, C = 6, P = 1, S = 0.04, + P_rate = 2, I = 50, I_rate = 2 / 3, compatibility_threshold = 0.6, eps = 1e-8) { - #initialize matrix with coordinates - edges_xy <- convert_edges(object,xy) - m <- nrow(edges_xy) + # initialize matrix with coordinates + edges_xy <- convert_edges(object, xy) + m <- nrow(edges_xy) - #initialize edge subdivision list - elist <- unname(lapply(split(edges_xy, rep(1:nrow(edges_xy), ncol(edges_xy))), - function(y) matrix(y,2,2,byrow = TRUE))) + # initialize edge subdivision list + elist <- unname(lapply( + split(edges_xy, rep(seq_len(nrow(edges_xy)), ncol(edges_xy))), + function(y) matrix(y, 2, 2, byrow = TRUE) + )) - #main force bundling routine - elist <- force_bundle_iter(edges_xy,elist,K,C,P,P_rate, - S,I, I_rate,compatibility_threshold, eps) + # main force bundling routine + elist <- force_bundle_iter( + edges_xy, elist, K, C, P, P_rate, + S, I, I_rate, compatibility_threshold, eps + ) - # assemble data frame - segments <- nrow(elist[[1]]) + # assemble data frame + segments <- nrow(elist[[1]]) - idx <- seq(0, 1, length.out = segments) - data_bundle <- as.data.frame(cbind( - do.call("rbind",elist), - rep(idx,m), - rep(1:m,each=segments))) + idx <- seq(0, 1, length.out = segments) + data_bundle <- as.data.frame(cbind( + do.call("rbind", elist), + rep(idx, m), + rep(1:m, each = segments) + )) - names(data_bundle) <- c("x","y","index","group") + names(data_bundle) <- c("x", "y", "index", "group") - data_bundle + data_bundle } #' @importFrom Rcpp sourceCpp diff --git a/R/bundle_hammer.R b/R/bundle_hammer.R index 16d5054..69794d2 100644 --- a/R/bundle_hammer.R +++ b/R/bundle_hammer.R @@ -14,40 +14,38 @@ #' @seealso [edge_bundle_force],[edge_bundle_stub], [edge_bundle_path] #' @export -edge_bundle_hammer <- function(object,xy,bw=0.05,decay=0.7){ - if (!requireNamespace('reticulate', quietly = TRUE)) { - stop('The `reticulate` package is required for this functionality') - } - if(any(class(object)=="igraph")){ - if (!requireNamespace('igraph', quietly = TRUE)) { - stop('The `igraph` package is required for this functionality') +edge_bundle_hammer <- function(object, xy, bw = 0.05, decay = 0.7) { + if (!requireNamespace("reticulate", quietly = TRUE)) { + stop("The `reticulate` package is required for this functionality") } - nodes <- data.frame(name=paste0("node",0:(igraph::vcount(object)-1)),x=xy[,1],y=xy[,2]) - el <- igraph::get.edgelist(object,names = FALSE) - el1 <- data.frame(source=el[,1]-1,target=el[,2]-1) - - } else if(any(class(object)=="tbl_graph")){ - if (!requireNamespace('tidygraph', quietly = TRUE)) { - stop('The `tidygraph` package is required for this functionality') + if (any(class(object) == "igraph")) { + if (!requireNamespace("igraph", quietly = TRUE)) { + stop("The `igraph` package is required for this functionality") + } + nodes <- data.frame(name = paste0("node", 0:(igraph::vcount(object) - 1)), x = xy[, 1], y = xy[, 2]) + el <- igraph::get.edgelist(object, names = FALSE) + el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1) + } else if (any(class(object) == "tbl_graph")) { + if (!requireNamespace("tidygraph", quietly = TRUE)) { + stop("The `tidygraph` package is required for this functionality") + } + object <- tidygraph::as.igraph(object) + nodes <- data.frame(name = paste0("node", 0:(igraph::vcount(object) - 1)), x = xy[, 1], y = xy[, 2]) + el <- igraph::get.edgelist(object, names = FALSE) + el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1) + } else if (any(class(object) == "network")) { + nodes <- data.frame(name = paste0("node", 0:(network::get.network.attribute(object, "n") - 1)), x = xy[, 1], y = xy[, 2]) + el <- network::as.edgelist(object) + el1 <- data.frame(source = el[, 1] - 1, target = el[, 2] - 1) + } else { + stop("only `igraph`, `network` or `tbl_graph` objects supported.") } - object <- tidygraph::as.igraph(object) - nodes <- data.frame(name=paste0("node",0:(igraph::vcount(object)-1)),x=xy[,1],y=xy[,2]) - el <- igraph::get.edgelist(object,names = FALSE) - el1 <- data.frame(source=el[,1]-1,target=el[,2]-1) - - } else if(any(class(object)=="network")){ - nodes <- data.frame(name=paste0("node",0:(network::get.network.attribute(object,"n")-1)),x=xy[,1],y=xy[,2]) - el <- network::as.edgelist(object) - el1 <- data.frame(source=el[,1]-1,target=el[,2]-1) - } else{ - stop("only `igraph`, `network` or `tbl_graph` objects supported.") - } - data_bundle <- shader_env$datashader_bundling$hammer_bundle(nodes,el1,initial_bandwidth = bw,decay = decay) - data_bundle$group <- is.na(data_bundle$y)+0 - data_bundle$group <- cumsum(data_bundle$group)+1 - data_bundle <- data_bundle[!is.na(data_bundle$y),] - data_bundle$index <- unlist(sapply(table(data_bundle$group),function(x) seq(0,1,length.out=x))) - data_bundle[,c("x","y","index","group")] + data_bundle <- shader_env$datashader_bundling$hammer_bundle(nodes, el1, initial_bandwidth = bw, decay = decay) + data_bundle$group <- is.na(data_bundle$y) + 0 + data_bundle$group <- cumsum(data_bundle$group) + 1 + data_bundle <- data_bundle[!is.na(data_bundle$y), ] + data_bundle$index <- unlist(sapply(table(data_bundle$group), function(x) seq(0, 1, length.out = x))) + data_bundle[, c("x", "y", "index", "group")] } #' @title install python dependencies for hammer bundling @@ -60,17 +58,17 @@ edge_bundle_hammer <- function(object,xy,bw=0.05,decay=0.7){ #' @export #' install_bundle_py <- function(method = "auto", conda = "auto") { - if (!requireNamespace('reticulate', quietly = TRUE)) { - stop('The `reticulate` package is required for this functionality') - } - reticulate::py_install("datashader", method = method, conda = conda, pip = TRUE) - reticulate::py_install("scikit-image", method = method, conda = conda, pip = TRUE) + if (!requireNamespace("reticulate", quietly = TRUE)) { + stop("The `reticulate` package is required for this functionality") + } + reticulate::py_install("datashader", method = method, conda = conda, pip = TRUE) + reticulate::py_install("scikit-image", method = method, conda = conda, pip = TRUE) } # Environment for globals shader_env <- new.env(parent = emptyenv()) .onLoad <- function(libname, pkgname) { - reticulate::configure_environment(pkgname) - assign("datashader_bundling", reticulate::import("datashader.bundling", delay_load = TRUE), shader_env) + reticulate::configure_environment(pkgname) + assign("datashader_bundling", reticulate::import("datashader.bundling", delay_load = TRUE), shader_env) } diff --git a/R/bundle_path.R b/R/bundle_path.R index 617cbfa..16cfd84 100644 --- a/R/bundle_path.R +++ b/R/bundle_path.R @@ -14,93 +14,99 @@ #' Wallinger, M., Archambault, D., Auber, D., Nollenburg, M., & Peltonen, J. (2021). Edge-Path Bundling: A Less Ambiguous Edge Bundling Approach. IEEE Transactions on Visualization and Computer Graphics. #' @examples #' library(igraph) -#' g <- graph_from_edgelist(matrix(c(1,2,1,6,1,4,2,3,3,4,4,5,5,6),ncol = 2,byrow = TRUE),FALSE) -#' xy <- cbind(c(0,10,25,40,50,50),c(0,15,25,15,0,-10)) -#' edge_bundle_path(g,xy) +#' g <- graph_from_edgelist(matrix(c( +#' 1, 2, 1, 6, +#' 1, 4, 2, 3, 3, 4, 4, 5, 5, 6 +#' ), ncol = 2, byrow = TRUE), FALSE) +#' xy <- cbind(c(0, 10, 25, 40, 50, 50), c(0, 15, 25, 15, 0, -10)) +#' edge_bundle_path(g, xy) #' @export -edge_bundle_path <- function(g,xy,max_distortion = 2,weight_fac = 2,segments = 20){ - # preprocess - if(!igraph::is.igraph(g)){ - stop("edge_bundle_path requires the input graph to be an ingraph object") - } - m <- igraph::ecount(g) - lock <- rep(FALSE,m) - skip <- rep(FALSE,m) - - el <- igraph::get.edgelist(g,names = FALSE) - exy <- cbind(xy[el[,1],1],xy[el[,1],2], - xy[el[,2],1],xy[el[,2],2]) - elen <- sqrt((exy[,1]-exy[,3])^2+(exy[,2]-exy[,4])^2) - weights <- elen^weight_fac - sortedEdges <- order(weights,decreasing = TRUE) - cpoints <- vector("list",m) - #iterate - for(e in sortedEdges){ - s <- el[e,1] - t <- el[e,2] - cpoints[[e]] <- xy[c(s,t),] - if(lock[e]){ - next() - } - skip[e] <- TRUE - g1 <- igraph::delete.edges(g,which(skip)) - sp_verts <- suppressWarnings(igraph::shortest_paths(g1,s,t,weights = weights[!skip])$vpath[[1]]) - if(length(sp_verts)<2){ - skip[e] <- FALSE - next +edge_bundle_path <- function(g, xy, max_distortion = 2, weight_fac = 2, segments = 20) { + # preprocess + if (!igraph::is.igraph(g)) { + stop("edge_bundle_path requires the input graph to be an ingraph object") } - sp_len <- path_length(sp_verts,xy) - if(sp_len >= max_distortion * elen[e]){ - skip[e] <- FALSE - next + m <- igraph::ecount(g) + lock <- rep(FALSE, m) + skip <- rep(FALSE, m) + + el <- igraph::get.edgelist(g, names = FALSE) + exy <- cbind( + xy[el[, 1], 1], xy[el[, 1], 2], + xy[el[, 2], 1], xy[el[, 2], 2] + ) + elen <- sqrt((exy[, 1] - exy[, 3])^2 + (exy[, 2] - exy[, 4])^2) + weights <- elen^weight_fac + sortedEdges <- order(weights, decreasing = TRUE) + cpoints <- vector("list", m) + # iterate + for (e in sortedEdges) { + s <- el[e, 1] + t <- el[e, 2] + cpoints[[e]] <- xy[c(s, t), ] + if (lock[e]) { + next() + } + skip[e] <- TRUE + g1 <- igraph::delete.edges(g, which(skip)) + sp_verts <- suppressWarnings(igraph::shortest_paths(g1, s, t, weights = weights[!skip])$vpath[[1]]) + if (length(sp_verts) < 2) { + skip[e] <- FALSE + next + } + sp_len <- path_length(sp_verts, xy) + if (sp_len >= max_distortion * elen[e]) { + skip[e] <- FALSE + next + } + lock[igraph::get.edge.ids(g, rep(as.integer(sp_verts), each = 2)[-c(1, 2 * length(sp_verts))])] <- TRUE + cpoints[[e]] <- xy[sp_verts, ] } - lock[igraph::get.edge.ids(g,rep(as.integer(sp_verts),each=2)[-c(1,2*length(sp_verts))])] <- TRUE - cpoints[[e]] <- xy[sp_verts,] - } - cpoints_bezier <- lapply(cpoints,approximateBezier,n=segments) + cpoints_bezier <- lapply(cpoints, approximateBezier, n = segments) - idx <- seq(0, 1, length.out = segments) - data_bundle <- as.data.frame(cbind( - do.call("rbind",cpoints_bezier), - rep(idx,m), - rep(1:m,each=segments))) + idx <- seq(0, 1, length.out = segments) + data_bundle <- as.data.frame(cbind( + do.call("rbind", cpoints_bezier), + rep(idx, m), + rep(1:m, each = segments) + )) - names(data_bundle) <- c("x","y","index","group") + names(data_bundle) <- c("x", "y", "index", "group") - data_bundle + data_bundle } -path_length <- function(verts,xy){ - plen <- 0 - for(i in 1:(length(verts)-1)){ - plen <- plen + sqrt((xy[i,1]-xy[i+1,1])^2+(xy[i,2]-xy[i+1,2])^2) - } - plen +path_length <- function(verts, xy) { + plen <- 0 + for (i in 1:(length(verts) - 1)) { + plen <- plen + sqrt((xy[i, 1] - xy[i + 1, 1])^2 + (xy[i, 2] - xy[i + 1, 2])^2) + } + plen } approximateBezier <- function(points, n) { - pnrow <- nrow(points)-1 - tseq <- seq(0,1,length.out=n) - if(pnrow==1){ - bezier <- cbind( - tseq*points[1,1]+(1-tseq)*points[2,1], - tseq*points[1,2]+(1-tseq)*points[2,2] - ) - } - binoms <- choose(pnrow,seq(0,pnrow)) - bezier <- matrix(0,length(tseq),2) - b <- 1 - for(t in tseq){ - p <- c(0,0) - for(i in 0:pnrow){ - tpi <- (1-t)^(pnrow-i) - coeff <- tpi*t^i - p[1] <- p[1] + binoms[i+1] * coeff * points[i+1,1] - p[2] <- p[2] + binoms[i+1] * coeff * points[i+1,2] + pnrow <- nrow(points) - 1 + tseq <- seq(0, 1, length.out = n) + if (pnrow == 1) { + bezier <- cbind( + tseq * points[1, 1] + (1 - tseq) * points[2, 1], + tseq * points[1, 2] + (1 - tseq) * points[2, 2] + ) + } + binoms <- choose(pnrow, seq(0, pnrow)) + bezier <- matrix(0, length(tseq), 2) + b <- 1 + for (t in tseq) { + p <- c(0, 0) + for (i in 0:pnrow) { + tpi <- (1 - t)^(pnrow - i) + coeff <- tpi * t^i + p[1] <- p[1] + binoms[i + 1] * coeff * points[i + 1, 1] + p[2] <- p[2] + binoms[i + 1] * coeff * points[i + 1, 2] + } + bezier[b, ] <- p + b <- b + 1 } - bezier[b,] <- p - b <- b+1 - } - bezier + bezier } diff --git a/R/bundle_stub.R b/R/bundle_stub.R index f3405f1..09ec43a 100644 --- a/R/bundle_stub.R +++ b/R/bundle_stub.R @@ -15,235 +15,235 @@ #' Nocaj, Arlind, and Ulrik Brandes. "Stub bundling and confluent spirals for geographic networks." International Symposium on Graph Drawing. Springer, Cham, 2013. #' @examples #' library(igraph) -#' g <- graph.star(10,"undirected") +#' g <- graph.star(10, "undirected") #' #' xy <- matrix(c( -#' 0,0, -#' cos(90*pi/180),sin(90*pi/180), -#' cos(80*pi/180),sin(80*pi/180), -#' cos(70*pi/180),sin(70*pi/180), -#' cos(330*pi/180),sin(330*pi/180), -#' cos(320*pi/180),sin(320*pi/180), -#' cos(310*pi/180),sin(310*pi/180), -#' cos(210*pi/180),sin(210*pi/180), -#' cos(200*pi/180),sin(200*pi/180), -#' cos(190*pi/180),sin(190*pi/180) -#'),ncol=2,byrow=TRUE) +#' 0, 0, +#' cos(90 * pi / 180), sin(90 * pi / 180), +#' cos(80 * pi / 180), sin(80 * pi / 180), +#' cos(70 * pi / 180), sin(70 * pi / 180), +#' cos(330 * pi / 180), sin(330 * pi / 180), +#' cos(320 * pi / 180), sin(320 * pi / 180), +#' cos(310 * pi / 180), sin(310 * pi / 180), +#' cos(210 * pi / 180), sin(210 * pi / 180), +#' cos(200 * pi / 180), sin(200 * pi / 180), +#' cos(190 * pi / 180), sin(190 * pi / 180) +#' ), ncol = 2, byrow = TRUE) #' -#' edge_bundle_stub(g,xy) +#' edge_bundle_stub(g, xy) #' # use ggforce::geom_bezier for plotting #' @export -edge_bundle_stub <- function(object,xy,alpha = 11,beta = 75,gamma = 40,t = 0.5,tshift = 0.5){ - if(any(class(object)=="igraph")){ - if (!requireNamespace('igraph', quietly = TRUE)) { - stop('The `igraph` package is required for this functionality') +edge_bundle_stub <- function(object, xy, alpha = 11, beta = 75, gamma = 40, t = 0.5, tshift = 0.5) { + if (any(class(object) == "igraph")) { + if (!requireNamespace("igraph", quietly = TRUE)) { + stop("The `igraph` package is required for this functionality") + } + el <- igraph::get.edgelist(object, FALSE) + adj <- igraph::get.adjlist(object, "all") + } else if (any(class(object) == "tbl_graph")) { + if (!requireNamespace("tidygraph", quietly = TRUE)) { + stop("The `tidygraph` package is required for this functionality") + } + object <- tidygraph::as.igraph(object) + el <- igraph::get.edgelist(object, FALSE) + adj <- igraph::get.adjlist(object, "all") + } else if (any(class(object) == "network")) { + el <- network::as.edgelist(object) + + stop("`network` objects not supported. Convert object to `igraph` or `tbl_graph` first.") + } else { + stop("only `igraph` or `tbl_graph` objects supported.") } - el <- igraph::get.edgelist(object,FALSE) - adj <- igraph::get.adjlist(object,"all") - } else if(any(class(object)=="tbl_graph")){ - if (!requireNamespace('tidygraph', quietly = TRUE)) { - stop('The `tidygraph` package is required for this functionality') + idx <- seq(0, 1, length.out = 8) # 4 control points per stub/2 stubs + + B <- compute_bundle_list(xy, adj, gamma, alpha) + + res <- vector("list", nrow(el)) + for (e in seq_len(nrow(el))) { + v <- el[e, 1] + w <- el[e, 2] + idw <- which(adj[[w]] == v) + Bw <- adj[[w]][which(B[[w]] == B[[w]][idw])] + idv <- which(adj[[v]] == w) + Bv <- adj[[v]][which(B[[v]] == B[[v]][idv])] + pv <- xy[v, ] + pw <- xy[w, ] + tst <- as.data.frame(control_points(xy, pv, pw, Bv, Bw, tshift, beta, t)) + tst$idx <- idx + tst$grp <- c(rep(paste0(e, ".", 1), 4), rep(paste0(e, ".", 2), 4)) + rownames(tst) <- NULL + res[[e]] <- tst } - object <- tidygraph::as.igraph(object) - el <- igraph::get.edgelist(object,FALSE) - adj <- igraph::get.adjlist(object,"all") - - } else if(any(class(object)=="network")){ - el <- network::as.edgelist(object) - - stop("`network` objects not supported. Convert object to `igraph` or `tbl_graph` first.") - } else{ - stop("only `igraph` or `tbl_graph` objects supported.") - } - - idx <- seq(0, 1, length.out = 8) #4 control points per stub/2 stubs - - B <- compute_bundle_list(xy,adj,gamma,alpha) - - res <- vector("list",nrow(el)) - for(e in 1:nrow(el)){ - v <- el[e,1] - w <- el[e,2] - idw <- which(adj[[w]]==v) - Bw <- adj[[w]][which(B[[w]]==B[[w]][idw])] - idv <- which(adj[[v]]==w) - Bv <- adj[[v]][which(B[[v]]==B[[v]][idv])] - pv <- xy[v,] - pw <- xy[w,] - tst <- as.data.frame(control_points(xy,pv,pw,Bv,Bw,tshift,beta,t)) - tst$idx <- idx - tst$grp <- c(rep(paste0(e,".",1),4),rep(paste0(e,".",2),4)) - rownames(tst) <- NULL - res[[e]] <- tst - } - data_bundle <- do.call(rbind,res) - names(data_bundle) <- c("x","y","index","group") - - data_bundle + data_bundle <- do.call(rbind, res) + names(data_bundle) <- c("x", "y", "index", "group") + + data_bundle } -euclidean_dist <- function(p,q){ - sqrt(sum((p-q)^2)) +euclidean_dist <- function(p, q) { + sqrt(sum((p - q)^2)) } -angle_edge <- function(P,Q){ - pvec <- c(P[3]-P[1],P[4]-P[2]) - qvec <- c(Q[3]-Q[1],Q[4]-Q[2]) - # dot_pq <- sum(pvec*qvec) - # mag_pq <- sum(pvec^2)*sum(qvec^2) - # acos(dot_pq/mag_pq)*180/pi - - m1 <- pvec[2]/pvec[1] - m2 <- qvec[2]/qvec[1] - a <- atan((m1 - m2 ) / (1+ m1*m2)) - aa <- c(a*180/pi,-a*180/pi) - aa[aa<0] <- 360 + aa[aa<0] - aa <- aa*pi/180 - min(aa) +angle_edge <- function(P, Q) { + pvec <- c(P[3] - P[1], P[4] - P[2]) + qvec <- c(Q[3] - Q[1], Q[4] - Q[2]) + # dot_pq <- sum(pvec*qvec) + # mag_pq <- sum(pvec^2)*sum(qvec^2) + # acos(dot_pq/mag_pq)*180/pi + + m1 <- pvec[2] / pvec[1] + m2 <- qvec[2] / qvec[1] + a <- atan((m1 - m2) / (1 + m1 * m2)) + aa <- c(a * 180 / pi, -a * 180 / pi) + aa[aa < 0] <- 360 + aa[aa < 0] + aa <- aa * pi / 180 + min(aa) } -angle_edge_vec <- Vectorize( angle_edge ) - -bundle_edges <- function(edges_xy,gamma,alpha){ - if(nrow(edges_xy)==1){ - return(1) - # return(t(c(1,1))) - } - dat <- as.data.frame(edges_xy) - dat[["V3"]] <- dat[["V3"]] - dat[["V1"]] - dat[["V4"]] <- dat[["V4"]] - dat[["V2"]] - dat[["V1"]] <- dat[["V2"]] <- 0 - elen <- sqrt(dat[["V3"]]^2+dat[["V4"]]^2) - dat[["x"]] <- dat[["V3"]]/elen - dat[["y"]] <- dat[["V4"]]/elen - dat[["angle"]] <- atan2(dat[["y"]],dat[["x"]])*180/pi - dat[["angle"]] <- ifelse(dat[["angle"]]<0,360+dat[["angle"]],dat[["angle"]]) - aord <- order(dat[["angle"]]) - angle_vec <- dat[["angle"]][aord] - w <- pmin(abs(angle_vec-angle_vec[c(2:length(angle_vec),1)]), - 360+angle_vec[c(2:length(angle_vec),1)]-angle_vec) - - if(length(w)==1){ - if(w>alpha){ - return(c(1,2)) - } else{ - return(c(1,1)) +angle_edge_vec <- Vectorize(angle_edge) + +bundle_edges <- function(edges_xy, gamma, alpha) { + if (nrow(edges_xy) == 1) { + return(1) + # return(t(c(1,1))) } - } - - start <- which.min(w) - bundles <- rep(0,length(aord)) - if(start!=1){ - w <- w[c(start:length(w),1:(start-1))] - aord <- aord[c(start:length(w),1:(start-1))] - } - bundles[aord[1]] <- 1 - bundles[aord[2]] <- 1 - cur <- 1 - for(i in 2:length(w)-1){ - if(w[i] alpha) { + return(c(1, 2)) + } else { + return(c(1, 1)) + } } - } - # pa <- graph.ring(nrow(edges_xy),directed = FALSE) - # E(pa)$weight <- 360-w - # bundles <- cluster_louvain(pa,weights = E(pa)$weight)$membership[order(aord)] + start <- which.min(w) + bundles <- rep(0, length(aord)) + if (start != 1) { + w <- w[c(start:length(w), 1:(start - 1))] + aord <- aord[c(start:length(w), 1:(start - 1))] + } + bundles[aord[1]] <- 1 + bundles[aord[2]] <- 1 + cur <- 1 + for (i in 2:length(w) - 1) { + if (w[i] < alpha && sum(bundles[bundles == cur]) < gamma) { + bundles[aord[i + 1]] <- cur + } else { + cur <- cur + 1 + bundles[aord[i + 1]] <- cur + } + } - bundles + # pa <- graph.ring(nrow(edges_xy),directed = FALSE) + # E(pa)$weight <- 360-w + # bundles <- cluster_louvain(pa,weights = E(pa)$weight)$membership[order(aord)] + + bundles } -compute_bundle_list <- function(xy,adj,gamma,alpha){ - bundles_lst <-lapply(1:length(adj),function(i){ - if(length(adj[[i]])>1){ - edges_xy <- (cbind(matrix(xy[i,],nrow=length(adj[[i]]),ncol=2,byrow=TRUE),xy[adj[[i]],])) - }else{ - edges_xy <- t((c(xy[i,],xy[adj[[i]],]))) - } - bundle_edges(edges_xy,gamma,alpha) - }) - bundles_lst +compute_bundle_list <- function(xy, adj, gamma, alpha) { + bundles_lst <- lapply(seq_len(length(adj)), function(i) { + if (length(adj[[i]]) > 1) { + edges_xy <- (cbind(matrix(xy[i, ], nrow = length(adj[[i]]), ncol = 2, byrow = TRUE), xy[adj[[i]], ])) + } else { + edges_xy <- t((c(xy[i, ], xy[adj[[i]], ]))) + } + bundle_edges(edges_xy, gamma, alpha) + }) + bundles_lst } -weighted_midpoint <- function(pv,pw,Bv,Bw,tshift){ - 0.5*(pv+pw)+(Bv/(Bv+Bw)-0.5)*tshift*(pw-pv) +weighted_midpoint <- function(pv, pw, Bv, Bw, tshift) { + 0.5 * (pv + pw) + (Bv / (Bv + Bw) - 0.5) * tshift * (pw - pv) } -centroid <- function(Bv,xy){ - if(length(Bv)>1){ - colMeans(xy[Bv,]) - } else{ - xy[Bv,] - } +centroid <- function(Bv, xy) { + if (length(Bv) > 1) { + colMeans(xy[Bv, ]) + } else { + xy[Bv, ] + } } -control_points <- function(xy,pv,pw,Bv,Bw,tshift,beta,t){ - #pv2 is the point on c such that angle(pv,pv2,pm)=beta - beta <- beta*pi/180 +control_points <- function(xy, pv, pw, Bv, Bw, tshift, beta, t) { + # pv2 is the point on c such that angle(pv,pv2,pm)=beta + beta <- beta * pi / 180 - #weighted midpoints - pm <- weighted_midpoint(pv,pw,length(Bv),length(Bw),tshift) - cv <- centroid(Bv,xy) - cw <- centroid(Bw,xy) + # weighted midpoints + pm <- weighted_midpoint(pv, pw, length(Bv), length(Bw), tshift) + cv <- centroid(Bv, xy) + cw <- centroid(Bw, xy) - #line between pv and cv - slope_v <- (cv[2] - pv[2])/(cv[1] - pv[1]) - incep_v <- cv[2] - slope_v * cv[1] + # line between pv and cv + slope_v <- (cv[2] - pv[2]) / (cv[1] - pv[1]) + incep_v <- cv[2] - slope_v * cv[1] - #line between pw and cw - slope_w <- (cw[2] - pw[2])/(cw[1] - pw[1]) - incep_w <- cw[2] - slope_w * cw[1] + # line between pw and cw + slope_w <- (cw[2] - pw[2]) / (cw[1] - pw[1]) + incep_w <- cw[2] - slope_w * cw[1] - # control points v ----------------------------------------------------------- - #angle pm-pv-cv - alpha <- angle_edge(c(pv,pm),c(pv,cv))*pi/180 - #remaining angle in triangle - gamma <- pi - alpha - beta + # control points v ----------------------------------------------------------- + # angle pm-pv-cv + alpha <- angle_edge(c(pv, pm), c(pv, cv)) * pi / 180 + # remaining angle in triangle + gamma <- pi - alpha - beta - #length of side opposite of beta - b <- euclidean_dist(pv,pm) + # length of side opposite of beta + b <- euclidean_dist(pv, pm) - l_pvcv <- sin(gamma) * b/sin(beta) + l_pvcv <- sin(gamma) * b / sin(beta) - pv2x <- c(pv[1] + l_pvcv/(sqrt(1+slope_v^2)),pv[1] - l_pvcv/(sqrt(1+slope_v^2))) - pv2y <- slope_v * pv2x + incep_v + pv2x <- c(pv[1] + l_pvcv / (sqrt(1 + slope_v^2)), pv[1] - l_pvcv / (sqrt(1 + slope_v^2))) + pv2y <- slope_v * pv2x + incep_v - idx <- which.min(c(euclidean_dist(pm,c(pv2x[1],pv2y[1])),euclidean_dist(pm,c(pv2x[2],pv2y[2])))) - pv2 <- c(pv2x[idx],pv2y[idx]) + idx <- which.min(c(euclidean_dist(pm, c(pv2x[1], pv2y[1])), euclidean_dist(pm, c(pv2x[2], pv2y[2])))) + pv2 <- c(pv2x[idx], pv2y[idx]) - pv1 <- pv + t * (pv2 - pv) + pv1 <- pv + t * (pv2 - pv) - # control points w ----------------------------------------------------------- - #angle pm-pv-cv - alpha <- angle_edge(c(pw,pm),c(pw,cw))*pi/180 - #remaining angle in triangle - gamma <- pi - alpha - beta + # control points w ----------------------------------------------------------- + # angle pm-pv-cv + alpha <- angle_edge(c(pw, pm), c(pw, cw)) * pi / 180 + # remaining angle in triangle + gamma <- pi - alpha - beta - #length of side opposite of beta - b <- euclidean_dist(pw,pm) + # length of side opposite of beta + b <- euclidean_dist(pw, pm) - l_pwcw <- sin(gamma) * b/sin(beta) + l_pwcw <- sin(gamma) * b / sin(beta) - pw2x <- c(pw[1] + l_pwcw/(sqrt(1+slope_w^2)),pw[1] - l_pwcw/(sqrt(1+slope_w^2))) - pw2y <- slope_w * pw2x + incep_w + pw2x <- c(pw[1] + l_pwcw / (sqrt(1 + slope_w^2)), pw[1] - l_pwcw / (sqrt(1 + slope_w^2))) + pw2y <- slope_w * pw2x + incep_w - idx <- which.min(c(euclidean_dist(pm,c(pw2x[1],pw2y[1])),euclidean_dist(pm,c(pw2x[2],pw2y[2])))) - pw2 <- c(pw2x[idx],pw2y[idx]) + idx <- which.min(c(euclidean_dist(pm, c(pw2x[1], pw2y[1])), euclidean_dist(pm, c(pw2x[2], pw2y[2])))) + pw2 <- c(pw2x[idx], pw2y[idx]) - pw1 <- pw + t * (pw2 - pw) + pw1 <- pw + t * (pw2 - pw) - # control point both - pvwm <- 0.5*(pv2 + pw2) + # control point both + pvwm <- 0.5 * (pv2 + pw2) - rbind( - rbind(pv,pv1,pv2,pvwm), - # rbind(pvwm,pw2,pw1,pw) - rbind(pw,pw1,pw2,pvwm) - ) + rbind( + rbind(pv, pv1, pv2, pvwm), + # rbind(pvwm,pw2,pw1,pw) + rbind(pw, pw1, pw2, pvwm) + ) } diff --git a/R/convert_edges.R b/R/convert_edges.R index 4afe78f..71335b0 100644 --- a/R/convert_edges.R +++ b/R/convert_edges.R @@ -6,54 +6,54 @@ #' @author David Schoch #' @export #' -convert_edges <- function(object,coords) UseMethod("convert_edges") +convert_edges <- function(object, coords) UseMethod("convert_edges") #' @rdname convert_edges #' @method convert_edges default #' @export -convert_edges.default <- function(object,coords){ - stop("don't know how to handle class ", dQuote(data.class(object))) +convert_edges.default <- function(object, coords) { + stop("don't know how to handle class ", dQuote(data.class(object))) } #' @rdname convert_edges #' @method convert_edges igraph #' @export -convert_edges.igraph <- function(object,coords){ - if (!requireNamespace('igraph', quietly = TRUE)) { - stop('The `igraph` package is required for this functionality') - } - el <- igraph::as_edgelist(object,names = FALSE) - if(igraph::vcount(object)!=nrow(coords)){ - stop('number of rows in `coords` does not match number of vertices') - } - edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2]) - edges_xy +convert_edges.igraph <- function(object, coords) { + if (!requireNamespace("igraph", quietly = TRUE)) { + stop("The `igraph` package is required for this functionality") + } + el <- igraph::as_edgelist(object, names = FALSE) + if (igraph::vcount(object) != nrow(coords)) { + stop("number of rows in `coords` does not match number of vertices") + } + edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2]) + edges_xy } #' @rdname convert_edges #' @method convert_edges network #' @export -convert_edges.network <- function(object,coords){ - if (!requireNamespace('network', quietly = TRUE)) { - stop('The `network` package is required for this functionality') - } - el <- network::as.edgelist(object) - if(network::get.network.attribute(object,"n")!=nrow(coords)){ - stop('number of rows in `coords` does not match number of vertices') - } - edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2]) - edges_xy +convert_edges.network <- function(object, coords) { + if (!requireNamespace("network", quietly = TRUE)) { + stop("The `network` package is required for this functionality") + } + el <- network::as.edgelist(object) + if (network::get.network.attribute(object, "n") != nrow(coords)) { + stop("number of rows in `coords` does not match number of vertices") + } + edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2]) + edges_xy } #' @rdname convert_edges #' @method convert_edges tbl_graph #' @export -convert_edges.tbl_graph <- function(object,coords){ - if (!requireNamespace('tidygraph', quietly = TRUE)) { - stop('The `tidygraph` package is required for this functionality') - } - el <- as.matrix(tidygraph::as_tibble(object,"edges")) - edges_xy <- cbind(coords[el[,1],1],coords[el[,1],2],coords[el[,2],1],coords[el[,2],2]) - edges_xy +convert_edges.tbl_graph <- function(object, coords) { + if (!requireNamespace("tidygraph", quietly = TRUE)) { + stop("The `tidygraph` package is required for this functionality") + } + el <- as.matrix(tidygraph::as_tibble(object, "edges")) + edges_xy <- cbind(coords[el[, 1], 1], coords[el[, 1], 2], coords[el[, 2], 1], coords[el[, 2], 2]) + edges_xy } diff --git a/R/data.R b/R/data.R index 9e60e67..022fcaa 100644 --- a/R/data.R +++ b/R/data.R @@ -26,20 +26,3 @@ #' @references #' Kujala, Rainer, et al. "A collection of public transport network data sets for 25 cities." Scientific data 5 (2018): 180089. "metro_berlin" - -# fl <- list.files("~/Documents/data/migration/",full.names = TRUE,pattern = "xls") -# map(fl,function(f){ -# df <- readxl::read_xls(f) -# names(df)[1] <- "first" -# idx <- which(df$first%in%state.name) -# idy <- which(df[6,]%in%state.name) -# tbl <- bind_cols(state.name,df[idx,idy]) -# names(tbl) <- c("from",state.name) -# tbl <- tbl %>% -# gather("to","weight",Alabama:Wyoming) %>% -# mutate(weight=as.numeric(weight)) %>% -# dplyr::filter(!is.na(weight)) %>% -# mutate(year=parse_number(f)) -# }) -> tbl_lst -# -# us_migration <- as.data.frame(do.call(rbind,tbl_lst)) diff --git a/R/flow_tnss.R b/R/flow_tnss.R index 2e702c7..e473a4c 100644 --- a/R/flow_tnss.R +++ b/R/flow_tnss.R @@ -33,54 +33,54 @@ tnss_dummies <- function(xy, root, ndiag = 50, ngrid = 50, nrand = 50) { - n <- nrow(xy) - verts <- 1:n - leafs <- setdiff(verts, root) - dat <- matrix(0, 0, 2) - - # circular points around leafs - if (circ) { - angles <- seq(0.01, 0.99 * 2 * pi, length.out = ncirc) - r <- rcirc - xy_circle <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] + r * cos(angles), xy[x, 2] + r * sin(angles)))) - dat <- rbind(dat, xy_circle) - } - - # points on line from source to leafs - if (line) { - tseq <- seq(0.2, 0.9, length.out = nline) - xy_lines <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] * tseq + xy[root, 1] * (1 - tseq), xy[x, 2] * tseq + xy[root, 2] * (1 - tseq)))) - dat <- rbind(dat, xy_lines) - } - - # diagonals through space - if (diag) { - pts_tr <- c(max(xy[, 1]), max(xy[, 2])) - pts_br <- c(max(xy[, 1]), min(xy[, 2])) - pts_bl <- c(min(xy[, 1]), min(xy[, 2])) - pts_tl <- c(max(xy[, 1]), max(xy[, 2])) - pts_extra <- rbind(pts_tr, pts_br, pts_bl, pts_tl) - tseq <- seq(0.1, 0.9, length.out = ndiag) - xy_extra <- do.call(rbind, lapply(1:4, function(x) cbind(pts_extra[x, 1] * tseq + xy[root, 1] * (1 - tseq), pts_extra[x, 2] * tseq + xy[root, 2] * (1 - tseq)))) - dat <- rbind(dat, xy_extra) - } - - # create an equidistant grid - if (grid) { - xdiff <- seq(min(xy[, 1]), max(xy[, 1]), length.out = ngrid) - ydiff <- seq(min(xy[, 2]), max(xy[, 2]), length.out = ngrid) - xy_grid <- as.matrix(expand.grid(xdiff, ydiff)) - colnames(xy_grid) <- NULL - dat <- rbind(dat, xy_grid) - } - - # some random points - if (rand) { - xy_rand <- cbind(stats::runif(nrand, min(xy[, 1]), max(xy[, 1])), stats::runif(50, min(xy[, 2]), max(xy[, 2]))) - dat <- rbind(dat, xy_rand) - } - - dat[!duplicated(dat), ] + n <- nrow(xy) + verts <- 1:n + leafs <- setdiff(verts, root) + dat <- matrix(0, 0, 2) + + # circular points around leafs + if (circ) { + angles <- seq(0.01, 0.99 * 2 * pi, length.out = ncirc) + r <- rcirc + xy_circle <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] + r * cos(angles), xy[x, 2] + r * sin(angles)))) + dat <- rbind(dat, xy_circle) + } + + # points on line from source to leafs + if (line) { + tseq <- seq(0.2, 0.9, length.out = nline) + xy_lines <- do.call(rbind, lapply(leafs, function(x) cbind(xy[x, 1] * tseq + xy[root, 1] * (1 - tseq), xy[x, 2] * tseq + xy[root, 2] * (1 - tseq)))) + dat <- rbind(dat, xy_lines) + } + + # diagonals through space + if (diag) { + pts_tr <- c(max(xy[, 1]), max(xy[, 2])) + pts_br <- c(max(xy[, 1]), min(xy[, 2])) + pts_bl <- c(min(xy[, 1]), min(xy[, 2])) + pts_tl <- c(max(xy[, 1]), max(xy[, 2])) + pts_extra <- rbind(pts_tr, pts_br, pts_bl, pts_tl) + tseq <- seq(0.1, 0.9, length.out = ndiag) + xy_extra <- do.call(rbind, lapply(1:4, function(x) cbind(pts_extra[x, 1] * tseq + xy[root, 1] * (1 - tseq), pts_extra[x, 2] * tseq + xy[root, 2] * (1 - tseq)))) + dat <- rbind(dat, xy_extra) + } + + # create an equidistant grid + if (grid) { + xdiff <- seq(min(xy[, 1]), max(xy[, 1]), length.out = ngrid) + ydiff <- seq(min(xy[, 2]), max(xy[, 2]), length.out = ngrid) + xy_grid <- as.matrix(expand.grid(xdiff, ydiff)) + colnames(xy_grid) <- NULL + dat <- rbind(dat, xy_grid) + } + + # some random points + if (rand) { + xy_rand <- cbind(stats::runif(nrand, min(xy[, 1]), max(xy[, 1])), stats::runif(50, min(xy[, 2]), max(xy[, 2]))) + dat <- rbind(dat, xy_rand) + } + + dat[!duplicated(dat), ] } #' @title Create Steiner tree from real and dummy points @@ -104,121 +104,121 @@ tnss_dummies <- function(xy, root, #' @export tnss_tree <- function(g, xy, xydummy, root, gamma = 0.9, epsilon = 0.3, elen = Inf, order = "random") { - xymesh <- rbind(xy, xydummy) - - n <- nrow(xy) - verts <- 1:n - leafs <- setdiff(verts, root) - - # triangulate points - tria <- interp::tri.mesh(xymesh[, 1], xymesh[, 2], duplicate = "remove") - - # create network - g1 <- igraph::graph_from_edgelist(rbind(tria$trlist[, 1:2], tria$trlist[, 2:3]), F) - g1 <- igraph::simplify(g1) - igraph::V(g1)$tnss <- "dummy" - igraph::V(g1)$tnss[seq_len(nrow(xy))] <- "real" - igraph::V(g1)$x <- xymesh[, 1] - igraph::V(g1)$y <- xymesh[, 2] - - # delete edges that are too long - el <- igraph::get.edgelist(g1, FALSE) - edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2]) - dist <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2)) - idx <- which(dist > elen) - if (length(idx) != 0) { - g1 <- igraph::delete.edges(g1, idx) - } - - # edge weights from distances - el <- igraph::get.edgelist(g1, FALSE) - edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2]) - igraph::E(g1)$weight <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2)) - - # calculate all shortest paths to eliminate dummy nodes - sp_nodes <- vector("list", length(leafs)) - sp_edges <- vector("list", length(leafs)) - k <- 0 - - g2 <- igraph::as.directed(g1, "mutual") - ide <- which(igraph::get.edgelist(g2, FALSE)[, 1] %in% leafs) - g2 <- igraph::delete.edges(g2, ide) - - - dist2root <- sqrt((xy[root, 1] - xy[leafs, 1])^2 + (xy[root, 2] - xy[leafs, 2])^2) - # minew <- min(igraph::E(g2)$weight) - # leafs_order <- leafs[order(dist2root,decreasing = TRUE)] - if (order == "near") { - leafs <- leafs[order(dist2root, decreasing = FALSE)] - } else if (order == "far") { - leafs <- leafs[order(dist2root, decreasing = TRUE)] - } else if (order == "weight") { - leafs <- leafs - } else { - leafs <- sample(leafs) - } - for (dest in leafs) { - k <- k + 1 - sp_list <- igraph::shortest_paths(g2, from = root, to = dest, weights = igraph::E(g2)$weight, output = "both") - sp_nodes[[k]] <- unlist(sp_list$vpath[[1]]) - sp_edges[[k]] <- unlist(sp_list$epath[[1]]) - - igraph::E(g2)$weight[sp_edges[[k]]] <- gamma * igraph::E(g2)$weight[sp_edges[[k]]] #+0.01*minew - } - del_nodes <- unique(unlist(sp_nodes)) - del_edges <- unique(unlist(sp_edges)) - - g3 <- igraph::delete.edges(g2, which(!((1:igraph::ecount(g2)) %in% del_edges))) - idx <- which(!igraph::V(g3) %in% del_nodes) - g3 <- igraph::delete.vertices(g3, idx) - - # straighten edges - xymesh1 <- xymesh[-idx, ] - g4 <- igraph::as.undirected(g3) - g4 <- igraph::delete_edge_attr(g4, "weight") - deg <- igraph::degree(g4) - del2 <- c() - for (dest in leafs) { - sp <- unlist(igraph::shortest_paths(g4, from = root, to = dest)$vpath) - keep <- which(duplicated( - rbind(xymesh1[sp, ], visvalingam(xymesh1[sp, ], epsilon = epsilon)), - fromLast = TRUE - )) - del <- sp[-keep] - del <- del[deg[del] == 2] - del2 <- c(del2, del) - } - - if (is.null(igraph::V(g4)$name)) { - igraph::V(g4)$name <- paste0("dummy_", 1:igraph::vcount(g4)) - } - if (!is.null(igraph::V(g)$name)) { - igraph::V(g4)$name[1:n] <- igraph::V(g)$name - } - del2_name <- igraph::V(g4)$name[unique(del2)] - g5 <- g4 - for (v in del2_name) { - ni <- igraph::neighborhood(g5, 1, v, mindist = 1) - g5 <- igraph::add.edges(g5, unlist(igraph::neighborhood(g5, 1, v, mindist = 1))) - g5 <- igraph::delete_vertices(g5, v) - } - - # calculate flow from edge weight - gfinal <- g5 - igraph::E(gfinal)$flow <- 0 - el <- igraph::get.edgelist(g, FALSE) - for (dest in leafs) { - ide <- which(el[, 1] == dest | el[, 2] == dest) - w <- igraph::E(g)$weight[ide] - sp <- igraph::shortest_paths(gfinal, root, dest, weights = NA, output = "epath") - - igraph::E(gfinal)$flow[unlist(sp$epath)] <- igraph::E(gfinal)$flow[unlist(sp$epath)] + w - } - gfinal$name <- "approx steiner tree" - igraph::V(gfinal)$tnss[root] <- "root" - igraph::V(gfinal)$tnss[leafs] <- "leaf" - class(gfinal) <- c("steiner_tree", class(gfinal)) - gfinal + xymesh <- rbind(xy, xydummy) + + n <- nrow(xy) + verts <- 1:n + leafs <- setdiff(verts, root) + + # triangulate points + tria <- interp::tri.mesh(xymesh[, 1], xymesh[, 2], duplicate = "remove") + + # create network + g1 <- igraph::graph_from_edgelist(rbind(tria$trlist[, 1:2], tria$trlist[, 2:3]), F) + g1 <- igraph::simplify(g1) + igraph::V(g1)$tnss <- "dummy" + igraph::V(g1)$tnss[seq_len(nrow(xy))] <- "real" + igraph::V(g1)$x <- xymesh[, 1] + igraph::V(g1)$y <- xymesh[, 2] + + # delete edges that are too long + el <- igraph::get.edgelist(g1, FALSE) + edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2]) + dist <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2)) + idx <- which(dist > elen) + if (length(idx) != 0) { + g1 <- igraph::delete.edges(g1, idx) + } + + # edge weights from distances + el <- igraph::get.edgelist(g1, FALSE) + edges_xy <- cbind(xymesh[el[, 1], 1], xymesh[el[, 1], 2], xymesh[el[, 2], 1], xymesh[el[, 2], 2]) + igraph::E(g1)$weight <- apply(edges_xy, 1, function(x) sqrt((x[1] - x[3])^2 + (x[2] - x[4])^2)) + + # calculate all shortest paths to eliminate dummy nodes + sp_nodes <- vector("list", length(leafs)) + sp_edges <- vector("list", length(leafs)) + k <- 0 + + g2 <- igraph::as.directed(g1, "mutual") + ide <- which(igraph::get.edgelist(g2, FALSE)[, 1] %in% leafs) + g2 <- igraph::delete.edges(g2, ide) + + + dist2root <- sqrt((xy[root, 1] - xy[leafs, 1])^2 + (xy[root, 2] - xy[leafs, 2])^2) + # minew <- min(igraph::E(g2)$weight) + # leafs_order <- leafs[order(dist2root,decreasing = TRUE)] + if (order == "near") { + leafs <- leafs[order(dist2root, decreasing = FALSE)] + } else if (order == "far") { + leafs <- leafs[order(dist2root, decreasing = TRUE)] + } else if (order == "weight") { + leafs <- leafs + } else { + leafs <- sample(leafs) + } + for (dest in leafs) { + k <- k + 1 + sp_list <- igraph::shortest_paths(g2, from = root, to = dest, weights = igraph::E(g2)$weight, output = "both") + sp_nodes[[k]] <- unlist(sp_list$vpath[[1]]) + sp_edges[[k]] <- unlist(sp_list$epath[[1]]) + + igraph::E(g2)$weight[sp_edges[[k]]] <- gamma * igraph::E(g2)$weight[sp_edges[[k]]] #+0.01*minew + } + del_nodes <- unique(unlist(sp_nodes)) + del_edges <- unique(unlist(sp_edges)) + + g3 <- igraph::delete.edges(g2, which(!((1:igraph::ecount(g2)) %in% del_edges))) + idx <- which(!igraph::V(g3) %in% del_nodes) + g3 <- igraph::delete.vertices(g3, idx) + + # straighten edges + xymesh1 <- xymesh[-idx, ] + g4 <- igraph::as.undirected(g3) + g4 <- igraph::delete_edge_attr(g4, "weight") + deg <- igraph::degree(g4) + del2 <- c() + for (dest in leafs) { + sp <- unlist(igraph::shortest_paths(g4, from = root, to = dest)$vpath) + keep <- which(duplicated( + rbind(xymesh1[sp, ], visvalingam(xymesh1[sp, ], epsilon = epsilon)), + fromLast = TRUE + )) + del <- sp[-keep] + del <- del[deg[del] == 2] + del2 <- c(del2, del) + } + + if (is.null(igraph::V(g4)$name)) { + igraph::V(g4)$name <- paste0("dummy_", 1:igraph::vcount(g4)) + } + if (!is.null(igraph::V(g)$name)) { + igraph::V(g4)$name[1:n] <- igraph::V(g)$name + } + del2_name <- igraph::V(g4)$name[unique(del2)] + g5 <- g4 + for (v in del2_name) { + ni <- igraph::neighborhood(g5, 1, v, mindist = 1) + g5 <- igraph::add.edges(g5, unlist(igraph::neighborhood(g5, 1, v, mindist = 1))) + g5 <- igraph::delete_vertices(g5, v) + } + + # calculate flow from edge weight + gfinal <- g5 + igraph::E(gfinal)$flow <- 0 + el <- igraph::get.edgelist(g, FALSE) + for (dest in leafs) { + ide <- which(el[, 1] == dest | el[, 2] == dest) + w <- igraph::E(g)$weight[ide] + sp <- igraph::shortest_paths(gfinal, root, dest, weights = NA, output = "epath") + + igraph::E(gfinal)$flow[unlist(sp$epath)] <- igraph::E(gfinal)$flow[unlist(sp$epath)] + w + } + gfinal$name <- "approx steiner tree" + igraph::V(gfinal)$tnss[root] <- "root" + igraph::V(gfinal)$tnss[leafs] <- "leaf" + class(gfinal) <- c("steiner_tree", class(gfinal)) + gfinal } #' @title Smooth a Steiner tree @@ -237,141 +237,141 @@ tnss_tree <- function(g, xy, xydummy, root, gamma = 0.9, epsilon = 0.3, elen = I #' @export tnss_smooth <- function(g, bw = 3, n = 10) { - if (!"steiner_tree" %in% class(g)) { - stop("g must be a steiner tree created with tnss_tree().") - } + if (!"steiner_tree" %in% class(g)) { + stop("g must be a steiner tree created with tnss_tree().") + } - root <- which(igraph::V(g)$tnss == "root") - leafs <- which(igraph::V(g)$tnss == "leaf") + root <- which(igraph::V(g)$tnss == "root") + leafs <- which(igraph::V(g)$tnss == "leaf") - el <- igraph::get.edgelist(g, names = FALSE) - xy <- cbind(igraph::V(g)$x, igraph::V(g)$y) + el <- igraph::get.edgelist(g, names = FALSE) + xy <- cbind(igraph::V(g)$x, igraph::V(g)$y) - ord <- order(igraph::distances(g, root, leafs), decreasing = TRUE) - res <- matrix(0, 0, 4) + ord <- order(igraph::distances(g, root, leafs), decreasing = TRUE) + res <- matrix(0, 0, 4) - for (v in ord) { - dest <- leafs[v] - sp <- unlist(igraph::shortest_paths(g, root, dest, output = "epath")$epath[[1]]) + for (v in ord) { + dest <- leafs[v] + sp <- unlist(igraph::shortest_paths(g, root, dest, output = "epath")$epath[[1]]) - path <- el[sp, , drop = FALSE] - edges_xy <- cbind(xy[path[, 1], 1], xy[path[, 1], 2], xy[path[, 2], 1], xy[path[, 2], 2]) - pathxy <- matrix(c(t(edges_xy)), ncol = 2, byrow = TRUE) - pathxy <- pathxy[!duplicated(pathxy), ] - flow <- igraph::E(g)$flow[sp] - flow <- c(flow, flow[length(flow)]) - pathxy <- cbind(pathxy, flow) + path <- el[sp, , drop = FALSE] + edges_xy <- cbind(xy[path[, 1], 1], xy[path[, 1], 2], xy[path[, 2], 1], xy[path[, 2], 2]) + pathxy <- matrix(c(t(edges_xy)), ncol = 2, byrow = TRUE) + pathxy <- pathxy[!duplicated(pathxy), ] + flow <- igraph::E(g)$flow[sp] + flow <- c(flow, flow[length(flow)]) + pathxy <- cbind(pathxy, flow) - if (nrow(pathxy) == 2) { - pathsm <- pathxy - } else { - pathsm <- edge_ksmooth(pathxy, bandwidth = bw, n = n) - } + if (nrow(pathxy) == 2) { + pathsm <- pathxy + } else { + pathsm <- edge_ksmooth(pathxy, bandwidth = bw, n = n) + } - res <- rbind(res, cbind(pathsm, dest)) - } - df <- as.data.frame(res, row.names = NA) - names(df) <- c("x", "y", "flow", "destination") - df + res <- rbind(res, cbind(pathsm, dest)) + } + df <- as.data.frame(res, row.names = NA) + names(df) <- c("x", "y", "flow", "destination") + df } # helpers ---- DouglasPeucker <- function(points, epsilon) { - dmax <- 0 - index <- 0 - end <- nrow(points) - ResultList <- numeric(0) - if (end < 3) { - return(ResultList <- rbind(ResultList, points)) - } - for (i in 2:(end - 1)) { - d <- ShortestDistance(points[i, ], line = rbind(points[1, ], points[end, ])) - if (d > dmax) { - index <- i - dmax <- d + dmax <- 0 + index <- 0 + end <- nrow(points) + ResultList <- numeric(0) + if (end < 3) { + return(ResultList <- rbind(ResultList, points)) + } + for (i in 2:(end - 1)) { + d <- ShortestDistance(points[i, ], line = rbind(points[1, ], points[end, ])) + if (d > dmax) { + index <- i + dmax <- d + } } - } - # if dmax is greater than epsilon recursively apply - if (dmax > epsilon) { - recResults1 <- DouglasPeucker(points[1:index, ], epsilon) - recResults2 <- DouglasPeucker(points[index:end, ], epsilon) - ResultList <- rbind(ResultList, recResults1, recResults2) - } else { - ResultList <- rbind(ResultList, points[1, ], points[end, ]) - } - ResultList <- as.matrix(ResultList[!duplicated(ResultList), ]) - colnames(ResultList) <- c("x", "p") - return(ResultList) + # if dmax is greater than epsilon recursively apply + if (dmax > epsilon) { + recResults1 <- DouglasPeucker(points[1:index, ], epsilon) + recResults2 <- DouglasPeucker(points[index:end, ], epsilon) + ResultList <- rbind(ResultList, recResults1, recResults2) + } else { + ResultList <- rbind(ResultList, points[1, ], points[end, ]) + } + ResultList <- as.matrix(ResultList[!duplicated(ResultList), ]) + colnames(ResultList) <- c("x", "p") + return(ResultList) } ShortestDistance <- function(p, line) { - x1 <- line[1, 1] - y1 <- line[1, 2] - x2 <- line[2, 1] - y2 <- line[2, 2] - x0 <- p[1] - y0 <- p[2] - d <- abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt((y2 - y1)^2 + (x2 - x1)^2) - return(as.numeric(d)) + x1 <- line[1, 1] + y1 <- line[1, 2] + x2 <- line[2, 1] + y2 <- line[2, 2] + x0 <- p[1] + y0 <- p[2] + d <- abs((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt((y2 - y1)^2 + (x2 - x1)^2) + return(as.numeric(d)) } point_distance <- function(x) { - d <- diff(x) - sqrt(d[, 1]^2 + d[, 2]^2) + d <- diff(x) + sqrt(d[, 1]^2 + d[, 2]^2) } seq_multiple <- function(start, end, n) { - f <- function(x, y, z) { - sq <- seq(from = x, to = y, length.out = z) - sq[-1] - } - sq_mult <- mapply(f, start, end, n, SIMPLIFY = FALSE) - c(start[1], do.call(c, sq_mult)) + f <- function(x, y, z) { + sq <- seq(from = x, to = y, length.out = z) + sq[-1] + } + sq_mult <- mapply(f, start, end, n, SIMPLIFY = FALSE) + c(start[1], do.call(c, sq_mult)) } densify <- function(x, n = 10L) { - stopifnot(is.matrix(x), ncol(x) == 3, nrow(x) > 1) - n_pts <- nrow(x) - n_vec <- rep(n + 1, n_pts - 1) - x_dense <- seq_multiple(start = x[1:(n_pts - 1), 1], end = x[2:n_pts, 1], n = n_vec) - y_dense <- seq_multiple(start = x[1:(n_pts - 1), 2], end = x[2:n_pts, 2], n = n_vec) - f_dense <- seq_multiple(start = x[1:(n_pts - 1), 3], end = x[2:n_pts, 3], n = n_vec) - # f_dense <- c(rep(x[-nrow(x),3],each=n)) - # f_dense <- c(f_dense,f_dense[length(f_dense)]) - - cbind(x_dense, y_dense, f_dense) + stopifnot(is.matrix(x), ncol(x) == 3, nrow(x) > 1) + n_pts <- nrow(x) + n_vec <- rep(n + 1, n_pts - 1) + x_dense <- seq_multiple(start = x[1:(n_pts - 1), 1], end = x[2:n_pts, 1], n = n_vec) + y_dense <- seq_multiple(start = x[1:(n_pts - 1), 2], end = x[2:n_pts, 2], n = n_vec) + f_dense <- seq_multiple(start = x[1:(n_pts - 1), 3], end = x[2:n_pts, 3], n = n_vec) + # f_dense <- c(rep(x[-nrow(x),3],each=n)) + # f_dense <- c(f_dense,f_dense[length(f_dense)]) + + cbind(x_dense, y_dense, f_dense) } edge_ksmooth <- function(x, smoothness = 1, bandwidth = 2, n = 10L) { - pad <- list( - start = rbind(x[1, ], 2 * x[1, ] - x[2, ]), - end = rbind(x[nrow(x), ], 2 * x[nrow(x), ] - x[nrow(x) - 1, ]) - ) - pad$start[2, 3] <- pad$start[1, 3] - - pad$start <- densify(pad$start[2:1, ], n = n) - pad$start[nrow(pad$start), 3] <- pad$start[nrow(pad$start) - 1, 3] - pad$end <- densify(pad$end, n = n) - - x_dense <- densify(x, n = n) - - n_pts <- nrow(x_dense) - x_pad <- rbind(pad$start, x_dense, pad$end) - d <- c(0, cumsum(point_distance(x_pad))) - x_smooth <- stats::ksmooth(d, x_pad[, 1], - n.points = length(d), - kernel = "normal", bandwidth = bandwidth - ) - y_smooth <- stats::ksmooth(d, x_pad[, 2], - n.points = length(d), - kernel = "normal", bandwidth = bandwidth - ) - keep_rows <- (x_smooth$x >= d[nrow(pad$start) + 1]) & - (x_smooth$x <= d[nrow(pad$start) + n_pts]) - x_new <- cbind(x_smooth$y, y_smooth$y, x_pad[, 3])[keep_rows, ] - x_new[1, ] <- pad$start[nrow(pad$start), ] - x_new[nrow(x_new), ] <- pad$end[1, ] - return(x_new) + pad <- list( + start = rbind(x[1, ], 2 * x[1, ] - x[2, ]), + end = rbind(x[nrow(x), ], 2 * x[nrow(x), ] - x[nrow(x) - 1, ]) + ) + pad$start[2, 3] <- pad$start[1, 3] + + pad$start <- densify(pad$start[2:1, ], n = n) + pad$start[nrow(pad$start), 3] <- pad$start[nrow(pad$start) - 1, 3] + pad$end <- densify(pad$end, n = n) + + x_dense <- densify(x, n = n) + + n_pts <- nrow(x_dense) + x_pad <- rbind(pad$start, x_dense, pad$end) + d <- c(0, cumsum(point_distance(x_pad))) + x_smooth <- stats::ksmooth(d, x_pad[, 1], + n.points = length(d), + kernel = "normal", bandwidth = bandwidth + ) + y_smooth <- stats::ksmooth(d, x_pad[, 2], + n.points = length(d), + kernel = "normal", bandwidth = bandwidth + ) + keep_rows <- (x_smooth$x >= d[nrow(pad$start) + 1]) & + (x_smooth$x <= d[nrow(pad$start) + n_pts]) + x_new <- cbind(x_smooth$y, y_smooth$y, x_pad[, 3])[keep_rows, ] + x_new[1, ] <- pad$start[nrow(pad$start), ] + x_new[nrow(x_new), ] <- pad$end[1, ] + return(x_new) } @@ -379,45 +379,45 @@ edge_ksmooth <- function(x, smoothness = 1, bandwidth = 2, n = 10L) { # credit to @coolbutuseless #' @importFrom utils head tail visvalingam <- function(points, epsilon) { - x <- points[, 1] - y <- points[, 2] - n <- max(c(floor(length(x) * epsilon), 2)) - # Sanity Check - stopifnot(length(x) == length(y)) - stopifnot(n <= length(x) && n >= 2) - if (length(x) == 2) { - return(list(x = x, y = y)) - } - - # Remove points - for (i in seq_len(length(x) - n)) { - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # Find areas - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - x1 <- head(x, -2) - x2 <- head(x[-1], -1) - x3 <- tail(x, -2) - - y1 <- head(y, -2) - y2 <- head(y[-1], -1) - y3 <- tail(y, -2) - - a0 <- x1 - x2 - a1 <- x3 - x2 - - b0 <- y1 - y2 - b1 <- y3 - y2 - - tri_areas <- abs(a0 * b1 - a1 * b0) # / 2 - - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # Find minimim area triangle and remove its centre vertex from all pts - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - min_tri <- which.min(tri_areas) - - x <- x[-(min_tri + 1)] - y <- y[-(min_tri + 1)] - } - - cbind(x, y) + x <- points[, 1] + y <- points[, 2] + n <- max(c(floor(length(x) * epsilon), 2)) + # Sanity Check + stopifnot(length(x) == length(y)) + stopifnot(n <= length(x) && n >= 2) + if (length(x) == 2) { + return(list(x = x, y = y)) + } + + # Remove points + for (i in seq_len(length(x) - n)) { + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # Find areas + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + x1 <- head(x, -2) + x2 <- head(x[-1], -1) + x3 <- tail(x, -2) + + y1 <- head(y, -2) + y2 <- head(y[-1], -1) + y3 <- tail(y, -2) + + a0 <- x1 - x2 + a1 <- x3 - x2 + + b0 <- y1 - y2 + b1 <- y3 - y2 + + tri_areas <- abs(a0 * b1 - a1 * b0) # / 2 + + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # Find minimim area triangle and remove its centre vertex from all pts + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + min_tri <- which.min(tri_areas) + + x <- x[-(min_tri + 1)] + y <- y[-(min_tri + 1)] + } + + cbind(x, y) } diff --git a/R/metro_multicriteria.R b/R/metro_multicriteria.R index ae7fd71..f9ecc76 100644 --- a/R/metro_multicriteria.R +++ b/R/metro_multicriteria.R @@ -21,60 +21,58 @@ #' # the algorithm has problems with parallel edges #' library(igraph) #' g <- simplify(metro_berlin) -#' xy <- cbind(V(g)$lon,V(g)$lat)*100 +#' xy <- cbind(V(g)$lon, V(g)$lat) * 100 #' #' # the algorithm is not very stable. try playing with the parameters -#' xy_new <- metro_multicriteria(g,xy,l = 2,gr = 0.5,w = c(100,100,1,1,100),bsize = 35) +#' xy_new <- metro_multicriteria(g, xy, l = 2, gr = 0.5, w = c(100, 100, 1, 1, 100), bsize = 35) #' @export -metro_multicriteria <- function(object,xy,l = 2,gr = 0.0025,w = rep(1,5),bsize = 5){ - n <- igraph::vcount(object) - # lg <- l*gr - adj <- as_adj_list1(object) - adj <- lapply(adj, function(x) x-1) - adj_deg2 <- adj[unlist(lapply(adj,length))==2] - el <- igraph::get.edgelist(object,FALSE) - 1 +metro_multicriteria <- function(object, xy, l = 2, gr = 0.0025, w = rep(1, 5), bsize = 5) { + adj <- as_adj_list1(object) + adj <- lapply(adj, function(x) x - 1) + adj_deg2 <- adj[unlist(lapply(adj, length)) == 2] + el <- igraph::get.edgelist(object, FALSE) - 1 - xy <- snap_to_grid(xy,gr) + xy <- snap_to_grid(xy, gr) - bbox <- station_bbox(xy,bsize,gr) + bbox <- station_bbox(xy, bsize, gr) - xy_new <- layout_as_metro_iter(adj,el,adj_deg2,xy,bbox,l,gr,w,bsize) - xy_new + xy_new <- layout_as_metro_iter(adj, el, adj_deg2, xy, bbox, l, gr, w, bsize) + xy_new } -#helper ---- -snap_to_grid <- function(xy,gr){ - xmin <- min(xy[,1]) - xmax <- max(xy[,1]) - ymin <- min(xy[,2]) - ymax <- max(xy[,2]) +# helper ---- +snap_to_grid <- function(xy, gr) { + xmin <- min(xy[, 1]) + xmax <- max(xy[, 1]) + ymin <- min(xy[, 2]) + ymax <- max(xy[, 2]) - deltax <- seq(xmin-4*gr,xmax+4*gr,by=gr) - deltay <- seq(ymin-4*gr,ymax+4*gr,by=gr) + deltax <- seq(xmin - 4 * gr, xmax + 4 * gr, by = gr) + deltay <- seq(ymin - 4 * gr, ymax + 4 * gr, by = gr) - xdiff <- outer(xy[,1],deltax,function(x,y) abs(x-y)) - ydiff <- outer(xy[,2],deltay,function(x,y) abs(x-y)) + xdiff <- outer(xy[, 1], deltax, function(x, y) abs(x - y)) + ydiff <- outer(xy[, 2], deltay, function(x, y) abs(x - y)) - xy_new <- cbind(deltax[apply(xdiff,1,which.min)],deltay[apply(ydiff,1,which.min)]) - dups <- duplicated(xy_new) - while(any(dups)){ - xy_new[which(dups),] <- xy_new[which(dups),] + c(sample(c(1,-1),1)*gr,sample(c(1,-1),1)*gr) + xy_new <- cbind(deltax[apply(xdiff, 1, which.min)], deltay[apply(ydiff, 1, which.min)]) dups <- duplicated(xy_new) - } - xy_new + while (any(dups)) { + xy_new[which(dups), ] <- xy_new[which(dups), ] + c(sample(c(1, -1), 1) * gr, sample(c(1, -1), 1) * gr) + dups <- duplicated(xy_new) + } + xy_new } -station_bbox <- function(xy,bsize,gr){ - cbind(xy - bsize * gr,xy + bsize * gr) +station_bbox <- function(xy, bsize, gr) { + cbind(xy - bsize * gr, xy + bsize * gr) } -as_adj_list1 <- function(g){ - n <- igraph::vcount(g) - lapply(1:n,function(i){ - x <- g[[i]][[1]] - attr(x,"env") <- NULL - attr(x,"graph") <- NULL - class(x) <- NULL - x - }) +as_adj_list1 <- function(g) { + n <- igraph::vcount(g) + lapply(1:n, function(i) { + x <- g[[i]][[1]] + attr(x, "env") <- NULL + attr(x, "graph") <- NULL + class(x) <- NULL + x + }) } diff --git a/cran-comments.md b/cran-comments.md index 50dc3af..3c4a871 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,11 +1,10 @@ -# Update to 0.4.1 +# Update to 0.4.2 -* Fixed broken URL in README -* bug fixes and non user facing changes +* code clean and maintenance ## Test environments -* Ubuntu 20.04: R 4.2.2 +* Ubuntu 22.04: R 4.3.2 * win-builder (devel and release) ## R CMD check results diff --git a/data/cali2010.rda b/data/cali2010.rda index bd75d45993a7cd8626a03a83f3847ee2710ba625..58bf33d43db58556d432083e2fc839bf1356dde1 100644 GIT binary patch literal 1229 zcmV;;1TyiR34#={-Pxsp@S#Q}mtW zqs=8mgdW3g?L1D~J*?MzZL8aM+HZxV_nr4;|6F!gakX*SA2YDldhQB|c-|(n$kBZR zfI%3C5MUl~^TvSsoliIr0k1A2%eV2|vSLL=}! z$`Qzzh+kMxqKY5Q-d`1nf97H@GZu?p7qFE_jUvR20r9iaXeHrV7f3xWX9ZF3>s*CJ zlx?Z0V6d+1Pn~!-9hqdUXsIXKvOE+w>Rd_mv3J9&^ilkA-N{ z9F&NKsU(!85eSV@B1DN2BuJ4WM2R15D1<1KrJ}e>f-Xr!iAp6TqEeKl5UBQ!Qh>~Q z#a*R=`>ZTf7bGPph=+wM95t(3C5O1N^Sdr$`o}qWMT}}NrGaxHS26Cfs<~GtR;7iB zd4-9AjhMu~xn8a3g~Q}@Dd=mgaa=VVS1I2(YZj$)gs7rbaJU=yHS<1J1?8`Sv_4N+ zn8h7;O_a=GF}tmsDV)S*e@4vqmoVCvyPo87+@||xtH9N(8s#HybLhPOPf=>GE-yAq z64G>jC5%g1?c2WZH!ofCyUpHiBXq{dhGTTAY`G?CGn6X)FMc-7;HIMkzLl<8T8mhh zo7)+KaKOsnF~0s52MM!sBkN&WX0z9sOlwtL645@RLB3kNWGI&1h}%9lVZv<}*6kOg z8Jgdy{Ojs1I#WhP9G!OURj`s$p<{zn3dWKozLLe|7|B$s)Rr^ZS2UGezIm#$v|3C4 z70fNuMYd#1!0@#V6Gh!iz-P5yQK;Fnv>Q{5xob^ohVe$4QralEq)?)nD7h~wDz_zy zFSjdM7Oq)wYc+vYVvCV0o0bI^GZFn&VOo^5TP8*o$}gB#C|b~JP*sttsnlpF+vRMU ztmdozylhY6*m#)+;C9CRaffVi2cjY#dQp7ny;U8z`G3V)>1Q?HI;a}yL&oGN*LfXf rP<;w0kEJNsOURR9jj-Cg#D2cRvkE37sUMD2|KjdQrwS4qa}FQ?P*ggs literal 1214 zcmV;v1VQ^kT4*^jL0KkKSr;?#mH-RKfB*mgOvF`J|NsC0-|+wM-||$Dh2#R!S*`qy zOb&3AKm*VPUH||9001DIBLr!n(*T&52rx`Qz>I*zFpLS4002w~G|P3XF&HCF2ABlIz(Imy1_WdV zA%tK|m;eG`L8gWyCIo4MX^6oZXf(hkCISo-5HKSkF$^ODWWWFu096qPCQT-q8e}z3 z8kvzNlPRXtDdSYenKEhWo|=Pd06k4MrXdkM1PfbtOcoL2=DkLEM zExWS%Uh6%s|E|+)?LF4{@id1(`=7m#qsBM{>Z z1CBnx&>ndEBhVay=#Eh3j$r1EVfM$QJu&JHAVucI}D2XsAs#ZH}r&VE%@2*}-^3{c`9p&|<<@JiMB()}TRC~&=v7!jv zA|Vn)G9?5-LM0NE5Je^;0%bDLN+55M1WXbAl88igovIPZltTK#g%nWv9rf{2JZz+1 zWi1x$J*28UNfsn%4~?FaK`n)$bc37B;VM22ZTC+47K^xbXYOp&xLQ|xgAMOf+`c!Q zABP7^ad6X-hXn>I^XeEZ+i%tk?WPKh6#L4AZ$}-ze6QxTd4rUZ5Va(dl%gRK!V)A% zks?Hi5+q2G)Jh==B`I!~2~b7JD3K{dq?Afhl%f?L;r5gUQ`su)EDr63ih|^Xr4bO= zO2(Sit)j!SvGh9aW%Z8Q@{1VMVM_w$Syr>Q#;WC9*;L0_VwU*+qbm& zo%Z)L5xiq%p_JY#+pbBNjjEAX;(QIeOjA*T;7ZoAT8mhho@I=ox?p5*7~O0v4kIS28r4@sv`zCtz*@XzluK~LZQmR9IE>}>xy#Xv%y3kGh4mJl38O7XGiAM1 zY$TMZSo-xatZ5=k=`3DRjFn2QNn<^cb4gXp=N(o?i%EZBa|?LUZP{sfTMa`*(RY&Y zneA5QH5sgp2HH6*mbdFvG}Ld>mU4@VMG7gGMagbbt9DqT`*yXKwQ|dgS&Ry+6kLg1 z-LNRRl#lAG3e=^d+A`QzD86A_p=&{@K~`H;Q>4&QxX9WyTUM+6yi8Bx^1r(*ZijGi zMSQ2+n{ diff --git a/data/metro_berlin.rda b/data/metro_berlin.rda index 4e25389671fd4fe479cfb25a6f889062d9548707..2e6396a80f9674109b4defc72fed758ff50d07df 100644 GIT binary patch literal 11366 zcmbWchgTEZ_daYF1cXF73I-Akp?AR^LP`~sc&-s5y?6cZ)Y|mb; z{~@tg`9Cc({)faK2^8|cN98Q#Jz0D9Bv=Vj>kl-#dOj0;2)XC#NxSD6)AP(kIRn4^ z@S=$9G54MnkwbBJ_xxj)xG!wao-29goD=KHcg`ePKRW%|B)9Y zM;q>RF%jz;%Y_4&{(J^q{Pr}J0Iz4Qa0&qw7$By)nog&wVL_@A1O$r?hGTFv3>V(O zfX4$!91D$!;lKd^!7V^w!89x?9{}LY0-O9)e^Ut?EZl|#1vtpqSX(F$%s^urVmQ!z z0F8jCz;ltNVgLsQ2J!(82m!?+;8>ceU77>}s{%vDbC3iK4bD6QL)kFJ(XfMHEQ$kW z#fxL%B%Ub^gJQiyaS(8{X-qyE7Z2y*np#;+v4yrwB2N`hWa8162u(CKBDuW4(^L{( zT47qCfoW(0zh?9nrio*b(=RZ11d&t`P?(m_^NNQNxD4)0hqyX$Z;KKa9xEdG@442EA>?c*3HZ5rmT(O9BdPf=p8&F7>1Z(HBG>geRFl(AWBGy>N?Pv29r$FAOn&x5T1K4Uj1N9 z5t(ia=h`x=R8U&dM-WwB+>>swNErOypb{iZfM*6m+A3TnG{y6nS{Zl^c(V(_3E?iT zV=2^ZAhAv$vjqusMlJ$I*^B2E!c1pqt~Pa!rgTh6CM-_hGF~e^eNcj2nyb2qDaJUi z@viJ^n!UszB(yt`*&t7biyOx38PZNnAqpQ-h8>Yi13SUk?wZ`jubHZn8TLuctkH=! zvn$$+hMgtmZuWxNNBFk4m73aI4<%7fKl_Wgi>2&PfqQwv&{L0*VEaU+=N=;l5v9)I zgPhpFV~{|L!u1lZj7)DKS;DLXuBY^|2alHLW-3!CEj!t)SP0z7dubml3oWfal$p|7 zemOX0umYi*sV+pM(7X0qr0qW!Xrv}$vTxOj60bYLkx6Kdwmrn-xM$r<24z)H$Zj`^ z9VrA?EM>oXM!^+pFYpJ?4B0~_FfOrIS7r%QbZx^D)8-McD|;(98ff*B9$qqT=1|NU zhqA64MCb?59Nl$fr%>(0k6B75dpLm?`(=dMD4pw5utXN&@sQLof>&OcHJF>-*jFZ< zbt7d6uYN0am@8o!<>-VV zpLujsL!wDGP5JPV57_J*S})XLFQp$}bII}-1{lmaV z5yWb0PGodn=ZWi*-9EhTFD4}_4ov4yOX=7r#2rjHqV9}!9nz{F&%Q4cCP7SB9JZR( zx0FF+5`@=h>~v%bXoAuAX%+vatw@s$wNHOW;_|xTfML2G=foDjUNAr`5T`d2VEPZH{SN zuMNP0fvKvg<^2dmZ&yD*i$ERa~-!JkVEjW>rFInq8S(rRz!8HyQ0G= zM~vp4d0R?~`v?)xQ7Vqtbu7UZt_d&@0?sdtdCbn3VoPVt3mrPHubi@2)`rkkf)fZ)r{)${D;5nj zbuij&SXe38%MuCvikBp-fvX1PAr#SkRT(sx3I?l)LqjoG6*L@Rq^g*L4Fi!Z%MJ;- ztbrF_t8^V6p+2QQDs!M5aZoIr5rw8r_dHU?w|dxV(wAz-pZkvEj&cDHp|Z zb4!($Jv9(MBWaMloih;!OM%o_pj@qxV$B^=zmJv}{5$=sB254V5uh*3$D5+5F-TC3 zpy|?VN-X1)>Zz^a8%g8%d@4%I0;O~j@}7efqoQn4iV$hp*;5a%ATvBuMnA$*0#6C2 ztXi1n={G2TJ$~|LD+U}ByC&6I>V&kZMQREV1f^mpr64X+`?jP~53SHsr9l#!xe|(q zcqS4at-f2}?kuci{Q;4v7;}RY^oYw908&?QG34f_5Lc`q6i5wXn+o)Fq99^DECE4* zn81HZ=TGP~KswJEQ*Q^olY3}{672<00;s6U->%{23B7d=E2keMtiKRea%sjYTdu)l zgwcDCuh=L;Ty5?twV#R(_2q3%M8B@(uIUH3w{<|1KiuBDi1>R1w)$93ePcX)c~V)7 zDjaxB@qFQy@imv%m8ZD3eq8l3bbE0@=xg2mWp4HP{d4-sw)6uiNWS-y18MQ9+mdnb z-^h?3SIyxwI4t4l7^?{bOEK+o8>&pdwMDz&Fngx*=lPCTcUFNL0w|M5Ru_R(Xbi>% zV?)Cn;2dBfaWNP)S_Q>mp!iooE6A!Nr={+tpZ9~ApbojhE3wbc*CgT_rq0AViOtCs z4xKF=dUAFxxwD?-gxwbP`0Y!psFZ4Ta(`>NTPSx1vmP7L0)I+(%-tNBd{?LoUgg3G z^ErioxZiUqZ=-DgEX(iK)oglioGLrFY=;XE@v>})a&_3b0eg4+h|IExQvW$!w@^Bn zPP?Ggf0m6LIA3!uyb`fGYXTTcj*K7QA9YLAbxK}Hh+j)oSwTO|yn4rEA5?|AiG2Yp z47Z&L32j$$)YwjjE{`tkpU^-Bm2BLouBa!;tizikXpMpw8 zYN;X3hhI^Wn%*ZJ96SB&u*Vm1)QQy?ziw9e(7!n>&6ArimUo`W$ODlNK6-X}o>*9@ z4%v1?ovX1mt3AtRhN@fsP`mS>)CFa9>vzT96ZXsL0=uHtzyGP;^r`CC-<%J#dlzYT z=@UoP;CrU;z{2sDq>)FmXO4x1JWmhR=&wSe$?0bqRA$N3@9F(e2JUwgF6mb~(ixKS zyK7uvs{1MS7v^YQ?Q10UVX;ZEIoX{f5q{}Q?da#zE!099vBI_FGX9d2BKlS0RBlwbfUZg}xv$5>(rZunPkcBwpYV0yZ|FKTX30i|43_ zW8hWMUhzdlUKEGqm9DKps&B4k(&IU1P=J8XXUoHgW)a9t z$Khx^htLjZaJ4vzYG|WVV}2iUg`=SuGl^AjzM0p6AO@r&ftG|F;=u0}C^n|zQ~9?? zSYTv<2nTEmMitv&ok;);r^_S)dEBpFR5_fvq>M)H$J6O?gW3C4Qs#%%V|j=;&$noo9D5nqR2o=Bv^Dw`G>Aeg#ynL zFEf;q9SO}w8xoz*S8tUdCvyrP;6FROjEYhkho}9yaq_2i!MwuKvcRYS)uh2(2m#B_lJt&$ZX&bdxPp!k#~|q=ofmIl7)* zO>riVk{#B|B(=18RfKq$G%gJ}6wew|hBdFP3AbFS*LA4*k+pTad94eDJfAw^>?1lc zn$5oB_cQO`uC6ZF?jv>@C80$l?70?YHm@RJs4`^k{0#}uR0WBuqqzcO9|`-Yaf|up z0vj=?bGl`?jkZ5mF()Wr&tB)HN%b5)8*b=3;dIf=Qy_&QD5wZVT0p4))XyDtE#;hv zchu0ykxU5ZVmpKygE^V|20UdHD-5V(V3z_UE>amF2uKku<^hg@Gz=V_DKS2@zbag6 z(We{Hsbrsqau>LsxV}3`q{{BYiUFzMUaGjaRT@~s-2LB!jdJ~W``ZoM=a&VN#&gfJ zm-K#mZ~B5tbS)awV46P$mQB50egLCz?iW*0r$On*xrc8pA!)vH!!z=! zlw5y!hm!g;i)#=LGKK)IjWD-1R%%t0WyIK6d^~cJffZv|(2yt$1hO~A!-A^K?tJ-Z zf4!U8kJ}vCJ?R%4%{*o2Q%mE}B5Tx7vC$AI$Porzj0J)~QBMK1oHWYgV`&OQ(E@$2 z9Rh)3E_i7>$GR%lFEN**MRN&zuC{c7)l+Q{Eb6+cEN)7Y>Vz*(=upSx zIqfz22?sESn}A`kY027a#FWFrDk6#b2JSum0Av3g-#?zT zWO_e6>IUiOiJE6O5FgYdZoCZESPqjg9vU(6BWziM8<+(A!x4Z~%RB;@Vj0@v(Dq24 ziklyb55I7ZU)NgTh}^KJcea%M7nv|qT$;X+KsvG$Gwp%PtvJD>JVvNtA6I}zV5i6r z>&%-`HBuf(FBuVkG;vAaFB??b*r%?Xmwl8z3EznxR zHF<q{wVn6%6 zzxB+>g)<%1v%i)nXwyJgrSX_kPSe?tgz*u(+3+F!)aQYT^DmzT{@nX}v}X9t!yJmt zKOaZ0>UP$em-!YC{Cu+f@jQLtm04_wJQ^4%XGJz6*hX`|I?^TS-;DRm~Lj zrVGbv!we7gojx=_U~DI_r6E#t`dLozne>LRSCd2I54&DT>5b)>+zs^;nKXQEG5^G+ z%&^kp=aV<<^;j@Cn$IgRa2Pa3 z5%=E!2FBO`9R4idI;cDtzaeGteuAeY|?JcwM*;#h+CBY<;&qlg1AP=LX2Vp%vo z#Q=kH_$4dIOPoIuOhewoLZ34DZ7l;p|JUGRRWLXh^!414e=irTXPr9m_h{4hgc38r zx~l73$BV#XqpDYu0Xy3lj*aLzskgm-=o@aQwxXclp?P-n&F$yuQE@*iTGawdW`~d0 z@Pazi-+OmY{zeD+Wsi)z-HO<}yPbsHa@=aG>u|X(*}3aS9G1UU`s#-GUtj&k*Uw-2 zioW_Rzr1n0Z>27`w>D(1${{>2K^@9(0q?Ua&*FM~q`gGfH`?g`motDrg2g=Vn zp+L$H&g#`G_(`V|zhw;##>XfL*I!HBTD+>tq9liJXnx;33K=>bk{kTtL`y|W{BCmi zrY7>+pwaUxuN=xVmq_!;;Gc+}IT7ED?GOx@3KGXR13&Qc_+3^W5I$W3{mA({gEue7 z%=J1F+N*DUD*Oej<{a*_?1LSDl4xuU14A)-k zWy@Nd!XHSde{sBN#b~6VW7zmxhs7&|_<;>kivckI=b608F^Biz2a|(FK5rMQ&5FCA za#g=8*`_V>`2J(Z^Z6q%T29Hy`^WBx1{uZR`=aiXuN}V)$ub1A{ZLD*ui7@sH$J8bz2Qls~w~PZYj9oXW7FSpY0W zTntByVS|REXcvKd{DvDQ0Ps@_q;iSB5<~@KK|{gBFb)l-V1R{-RuP~faa8_^%A@kd zP-r-auTCzMnhIy|zhOZISOWZ6K@3C?bCL!H11K1j4=XKz6f|6r?TImCPSJP^Z%N}H;g2#bSOz!+Af!yO* zO#mPqr~;d(JoG|;CraZUIw2u~C>R1tpC%=mBw!Z?^XXWH-`^!pn{CfcM7e5~fk$u+ z=3;bf@AY;`kOt&bZB+XO`aW$~~OG9xG=h&uJIUYV08Y!`6`XLpCLLOfg%at@TEQ+wJCmnwY z`#U?k{ST}0*GL*hA*VF|-r&W_$AotX>Xw4r#oww2K)a0iy6<$fCLE?>bARDN1hPQh$E4xTt3UPBb)iPzE;pM{j(p1E z<@wov_LkQ2dh|ChZcvUNX)5gKkg;rYa=OM+@P~#jPV2Q1ko>TOGO*CxOQI--f+sww z0DhKrLXnRu7b5Y$Cf*_S*Vtwcxt48E-$Ugrkyz-R@qNP6iA+^7|29S1(Oq$K*-s-V z_6(dIY&o94^d#0TJwWKkmYq0z=5&Pudf7hj_MyZdDIZMCENK7L2XLy2;>Um8V=+~g zCrb@a9k&U%ydPg*52RG6p!CtVG~3xk;GD6={UUUsCz)}6%9>Y5Au(rW^7~(4O>n$0 zf>b}jbwFz5Rkag=Kx(p;w#$%03=s#V(leyfnQo+YZYd`v?NE4vDe>z4ea$eT&I-H1 zOi!m<@ul2W7h)%|w^UuH$4%drm|>?`h9}q7@v=(m8`>~5@~>evE3|u=lvbGoX=P~H zP3IFavq-I?sLU?YOXG$nA`NLTj6g7PRcdJy9&*1u8a>C48R7aFfVs!leo#5s+}ARM z9Di5?I6vi$1+9{v|Hdc#?sPw08a)Hymlu<@+)(%oaBg+@RjiG-;ElI6UpT02Z9??d ziht-=USs61!l&J5aVEF#T<=^@cC0wMerghPaz!C*S3(H_r`x2rSXu#c-(J3HTN?y* zza3qw8)*cmE|PWs6l12QwOPxc`nUU{3va>%maBa(+-VMI^= ziHFE3J@#ZCANi+qxF5tu%SXZ0LKfu5$JB`I-0-GJV<7ct7d!?PQJ(tx9G99HywI<7 zt~~z%`=&z!yWy7l?&O7sUwy#Kp~{29#HBmOJ(B7H^Sql@e$#hOeT{4>bS1AW$Tj@l zviN~w5y9oTbC~=(gqgjq6q(#ie?kNON&55$BNoKz6fk!smG(g^J~2ri4)@f2pye2y z9-99r1Gf+G^%>b1vroXKP8jHWn-V+3>+pM2Y6}M>8@Um!Hb6aH`O8LkxJ&h_Ph!Wu z>4dmy_Rq>^WA*jxt&?_7mWFuR@@E`9oGhvsRV}(1$KW04(yYUJiHTV(&2 zS$U9`P&(~W3R*Vi>EFbj#8kq!>!R3H;NAyD9kLf~ z2IO{1PX(iJ5>YCO^I$Bd&&dZ_O0@Q?%|74jp$Ib`$*TbFFd+ zmr#Cbnqkg&v>e3Vt@Wesquei>4CigS>fV_!x9PAm>vU`cpJ zD0?XfEiIAxM(8L;B<{0Jr#ZQ+L4mx15=0|?!7g>fPN#LF0I8IY5qS(v5cWYP#EN@= z!$xOt8xxQ8(mP*C_S3tIDb12z&HB>nat{lqGh;l-fk*P3eV1ku0T%URLSAlpiC)Lt zgsY#Z-e$5!=e4udG_3l1XeR+^u!QAdqw~FF?ZD1{&%&D=RMivN>b_l>jL(DR*E*z! zi-bK%EGlayE&mL!G1;O)OTXKSFI3H1-AlE2Woeo`0|U%92{|w; zN#e5xo$Gw6*#3|SI6-L9OM~FGYEfm&p`gm+4c6m-wo+sC53b83SPdjtS*6)jWDIHQ zfGwl2rfH>$M4=@TSGL^Sd&a(NH=tR;s!a;~!S&|lGy4qoihcz27R2=v_u-iFP^N;B za$Cn(a=_m?7QWLXYjtk{#TGF%#Aju}gZjO_ddc^4|DITh+rI&zyhXKg(<+Lt@wacd zcMR3ILLD9Hj9A%_d}0N~RZpm+9kF?1+QkI*unFB;D{ixPbW7cTp56iGzVDS1vAW4U zo;rhpxIYSje^iT9v=sYT2}_D%c}k_a&NhsEa;W}FQworl9xO7qww$lP=-?U()-#p$ zIEkv6LhUe9u{D~907gGv>~;tnD?rE?ks2_B#mi%BbrPFuKkGNUHre%*=n(*dn0d5VFeN+ zXPpuX1UB2WD^@}HvoV=IwDi}nPRajAUHl)UKx6pb4*G6GTCX=scf!BOT~pl=E3OT5 zsbiDL39VErvY%64HZnlB%RsU0VR*~^kMxf)U9~EN7H1tG?K+=G_zb*{SEffp2UfGc zaLFhQh~x3z^s*b#V;q4|;XIDItkp9&*)Q~jDO4x9_K&`^SDLn@mZnOl_k+j+$r5>c z_O-G!g(+=;6$SrGN6P`^1-AgKnDh0uBP>WPjdBrR1v&sv87N1bIkf`VNRSlCt5vPSv-9(za9wBD{^$f z%*1P|xsVF*`9v-USNG}jWv}glcDN-Z<9;bL`B|zL^jyQPNa=_p`jw_ug&=swv1(pU zml~Ef=@i)Ym2B*^hggpy^>(lvsJ_S; zmCEBU9CeoM8KyHJsi8x%=rGQB zF3#6u8%2B&2bjpng^7hGg1Iccd8%n%(9eYW0vfz7%;3d9qW0s)iA0?(z!F~g$`ap7 zER1l3#FPRotcHYbMU)J@aiB!2QbJNkn$J2s3or6CO$CX?ZLnf-53`RD&rUlI81zWoyb!CI6D|~WO z8seap*Ik=tx(dQfXTWWaHjRxu6#_UnpX`{3NDPZrU07rds_{3gN_8Z)OAuA2&v=Ac z44nWA*miKB5DC)af!vq9jSbC(B-ZwJb1Xp`f8BV%Ru48SiQBmtW29))ic#8sOb9ky zTI6g~8*{6KJ1oEC(g2ehF_>hbuXAF}a#^0#k0qLV6d6w!4}e4$YM0~)&xT76!TU;H zb{*5IchBQPV?jR176kJf&0G~&CN1%~8@f~{9tLhhT|zPFjCMyX4Ybp1e2at90aKHf zz%HhEVNr>Snsq8&jxjTBF3_J=ZR^$4)R=`Y*1Ig%Chz06UnZL4qui@65^P&qCG2U| z{B_u?RCYj^B0L+Y0_p@O<@}8o=*{E*fO;&NAFk(l`$Ye=EmA7*ypq*TjL1u>f*a~R zQpEn~r4H5|i^cb{j|@6^dD;7zbi?~g6-2hNK=SDjE{8C?QJ1Vda$+ z4k<2hQ{m*U?^g4hz-*^L(dmX&x|8X`BmYPW*gR>#a2gS)e#Q;o|EABTuWr9wO6Uix z`iKn^zrY`Ylt}d3%7^$uYk>1!++Q%LCDW!bLOnd#k`K7ruL%Qa)e~p(wO?1-R8;BO z88ajh_Tg#!<|k~y@ZKd%vfJae2SI8 zV55=Cc`oFtly1!QD;?`|NI1EU$F{V_#DW3FkW3qPSseDWdm3Lwumg;gc?dW$7tYdp z<6jTOnn2jt$EjIW>R~MvMXQ}^5NHvGU%&|FnVyU|PYrDoXYYA`>X#N|^mLoKFcobs ztdyadn70n2i%{6 zxuUi!#lb%Y#r}Gpk9ca+iQbXp%6I*p)qhTY_$Kn;!D8$-Q(7dkool3*x<-?Kp|okd z`Qx@>2T3XLW_NMHJlSN930cEk0DZ^nGL(Mjc|Pv{&3(fAn;e9n1^bJ#a!)i^e4{Bf zcx%Sv14Z$|5YwT1pJ3g`Y# z{YFI2I7I(zdw9pNG1O0{$x#3O6Z6pIe*;9KJ^S!~X!)N_(OL_aga4S?lZT>U-xz$v^?}67XzPl}7JumoL{Z~i1LJ%Q{_ge@t| zq-+&;jj2VNxnwLlnEX8Va`Z`TZ%y<}`sqW=tnEStQ?oCVIq~u)4+C2ULcPtfQx~TC zt2RDOK5?&`zwE}0kgybf6?>c>5)mT2S=b^!4g`AyC2=a>?5%&|Gb8hW6R>!BzY}(G zv}U5F^7N>SAXPgJDj<=Frm~+=@cD{h?S4-t21dg)K}1NczxqQhl}@7=3;Mr zk@vHS<;TaW$OEwD%S}AP8?$DMs~RUtwa!P(y)rbLPjwBy!SK{mQ zu`mku_Iyp0cj57uZ++s~Ki}-^>isicUjAh~t0Z*Mfzp3ocd$4>vv}Bv`SSE(YKPBU z@7Y7`wYv$v%g6m!q!q-%)AJqum)xv&)rbGRdVj3s&BK}y-|?L^_SJ|76?aTyPbl^0 z&qXVocrJlty}p1%C0*~irCQ4nS^g>0;>sr zsL6?fpYl1hIdb3kvpX$OhHc?3LG_$(OLed_jmE`_hCMuByIg1+4XZ?9v1qbulk{9(=n%= z+xOUCvu@Rr5a`2HW$G}Ntk1IWlE5uXFhMmYuMQyAZuuOmcZUS@^D@;CYG|4&NV8vu zaa~+nTvq~yXSkWt5LHbq9kfLU0tkeBQ8lGlPgYDatjJvz=MbS0hpc6I-}f%!gOdCicTw`SepjucfoRzfpSJ;R7>E&mc zo#Xd0j{IGQU$RqCehW>rKYm-cHNEXyGDv@SI_{9foK>Mq6bCfVh)Vwx8N{8Qa51@V zc$gcUe3vvUM@bG|P!L>rvc4*n?MC`!C8ApG8t~1F?6N_QS6E+nH2GaAb-ZrJcmCVh zOX7@U@J2QTWL`>8n=?qi`Lp%A)4QTG&xeDR%YQsVCRo4Qdg2jZry7D^x9*x7fAl2S z`rFzBXza{>=uYy%zYbHS=daZJb^rRFwUcISD7O88+7Xykw3W=|jkE7X-1rrC^v$V> zhy&)wfBSr07#95$X}A&8oH=&p8E2?EEp(s?gi1p|Em4ssOmuaQ>O(?b4bILer1s=@ za9Y60ZkU1?^*aqFdvNxU&Q{`6*I*Ns0m#64rMqFz84(OZh>+x4h^_uABfJv1_5E#X zp-*Df5`#T$aKqzP5^V7*YtUgu{Ptnx-IB6ZO9!Dgg)U{vno-^b#A)+GzdJSB-j$ zo8t}Gi! zy}YxK;#+e5n|<-j%+T3YuRwvRc<834^HTPGukbGzMQYW8P2E{)O+*?g>f>vl_x$sy+OL*w7x8a+2y z>0bNdXpdH!oQW$G<^rp7${Rte(^sU>K0u8eY`;NJ;3|C=|e9%NgY(Hb}- o<sZUY_?4_`TQlp6kFl9PaZS_xgVB&p8M?PXohNc$fW`y`NvLm!nhl=o%q+lRWmU--;Elydr+beS{CivK95SqkfDc~P!-J(3aL zr_=T^&GVz@J{=v;k4IwuIeZQo`fBxeEQ3&Nv{6MYoPn|&=o?vf4qWd!QBFU5Mz91k zLy?5qv1rI*3=PUPi%%A0pAdJlt_+~w$H#$S=5_SuSTI21u&r?gMrtr=|m8udTbtsD?2ndQJ>5P1A8K;Q}#-mk;SVmq7hz1}5 zKpbk$0x|Mf_!JO>1t4PqQ$bBAXT6b)0U-G}>wHWB4glpK^J76Q8sr%Z17)zRafnzz z$Og()F=r5=47RBNv@AtH6=y1`N~{qlmQ{N96iNc#aUc#2s#?tj;)Iwfa44%;QZm01 zXeSa`wjet~Kw~V_A*BFAu@JH;Aq4bG!f{AWC7$3=L!uOB#ugHsoG=zlDJC0PL~!*- zP~xQE4n(eXK9vg;I1pKIJJvKifM(jvDRQu@M3%q_rtPNf1cJ4LU7Rhnn#hR-v&F$Z zlYJa|f1!mCo<#!}CVNNcbDG)Iis~F!RhTKwj=mlz0moPnaIR9StThlZ#qyjj%#SHF7{RMxEeKYDDRZn{IV7!PXA+W(?d3C)H7+N@okgk*3r6$p`{L4O4m;FA2io zDbX+-6i1@Lr9d^3pg6SMeLyYwusD`0$z)(TUjtxv)$}G-Pje^zeiKVI1&cF%R+bW% z?A9GZ%F$edX+QJFQqwb~AY*D)=cMU~JXd@K!Zh=$Hwb6QhT#(Xyq6@H?xqLqc z?RvMRv7S@Na8>07SOdsn7u6Jkdwv}khLP%Q<9ZP`AZ0BD{wUD;PGccc9EhM1+AYSJ zsxHL3HZ2S{-C}`3f@@udd7&9?7C;msZa%-7i_&7-0Lu8CSjluj$+8n&hu4Yt#~^%t zQ@l_*KbK4kRCEJi9(o8)DNEyoxVyAkM}B2M(P7)vJ_iOTw&+ZLb%T3~8`1p)Mq0=Q z?n-oMs4q6o3C ztc2vQ?j;dfIeYv#JjxAI*XeDicf3ndQzzatSJ&;J1F?w^khIP8Ig^WIN*VWR$=J|J zUdjxxQGq~?kFY=anHQ+xM#&8zD&NJ)E5RQXW}T~<_EEKAK&R+27WcK>rr0O`=2mC7 zcM%5ea=lza)vFtEty}xx#Wwd!3mKG~Q2`;ozTN8@Tny!A)ptfKUD!m{=moK_*mEdj z<8sRfj}n|(SB!~+XS#*VaiqJ9PKl??14S~NoL%S6t-zFzOzr|UdMSZ+ecEXmcTA37 z-V+W+3RY>Ev#yCw&h>$4!NEd&~1yCGWo zD=dfVd0z*S-rd@Y`6SfnmE;S4E&6d-tM2{yX0EA6`&$i~OyF=$P4ojVV6$QKLa}18sm%)*y zQ0+I$iLufpnk{uBws2Ik_qug}P&%x|c{z)0pDaoov1Qi(ryc~QQis1^RZktah$qmt|rAgpUdF-WV>@x=FZN63Xts5bU_$#Q2MAB&b5L8 zKxxSI0?ayx1uBlrr!wS=Gm8#Z_ zq{TTrcowRBnXR&Eeo5?DGE`CYoJ38Em}s=^2X*TjDDm>WQ45DZ(Qulwm>54X2gxFd zMMP{lRhlL|{d=?5av}yO<8MrhR(#YhRkc_=+74sewnZz7k;!6+ye5Pqe|I!%d(9lt zPJOZ%T)OOwA&5!M@*{=!J>*XT^G72j7&mm&GD+#}F*AV+RB^}*gAv;KS4l>GAX6#Ds|wORPvGKxa1RJ5*j%YVAa?Z1T+f1^upcc zsM6*4stNuG?Vei8cETZy`7afOx*JxM3yr2@%xiHhs+i~}FuwqmiZBQn4An=KQ3Fbn zk>`;t*#{R{e6b89KQiLQFH)f2#~x$*wX-4LLSIG4-^aT}28%{fPANUX974h{qR8O( z&rb>cb!OSV!u_GLHX&zrP_Ub6j#)5zEk&>3Q?;zAu9r;($HVU_C==m!J>V79t9W z+$BUUSx#+FdbxNk>yjL)&r7P!$hF;gck0c-d!CNRAG+{)EGxJ+8~i}-#hHF;o6d{3P~KdMRa7!LjEpXOEx+k(tBtAuBRsyUAol2vfJs;FbNhj1&#M#>Pt=K$ALmC zjO>iA(HVqq4?mDo>->|1Z9!b`W+%^@8U6@V?Tk!CMUpTm=ebX;D%Z!r?IPC=Rb$Qe zR%OK&dE4lq?J9v2wkgEiK~>dC(yTZqS8%2P=FrG0z~G?_T;PFz2_l!Yi>c^WCplC` zC6LvmBrx%4D5rq83aXZBmP)~ZM>(cju_+e&RW&66B&SC8P2qt8Od$gYSi4%t5a|I| zsy{m-qz{-g;?-sWYu+CG{PfaZh$K?Gd1Zg8*rNucPnGg(EV-#~Rm{ql zZ$n!hS{__(={X%oTvLEAIC!i3sb{@B%RExwYj}DrlAU*_#3j*)Y+|uYHm0p!q!yIJP*S9r0H|G@*?~)36?WNrdre=xO@&i0w+)fS%N$ z1Tr!#Rc~;3T2E=aUTlk5Y=n12Yj~phO@f1NrPz|f*THAbp8gfacfU^V-M@Z$Mi$U>^Hi8dVflshXWNZm7m?L;Y%a;Rj0}XT6CE`2{zPan!f z<_hB5?AFBXnwqL7vLn*wL$-bHM)W!>rA?&gQ_a%8g&!<$8hP4WZD6+Vi~PaeD%%{3 z-KrbvJr&1XCLP+#GG6v~It{NM3xBWL2iDrt4w4B{3c6%12$p?4OB8E-9RQX)KBRKb zW6!tXSmkQAnxhA0^V0WYx5myFe&g>}1c5FoLJ;QYL&>1%gjH(%OFr!A)8IK-dx%b( zG*VNPh7{L->ExGzZAbL6tfdsK`DuzEF_5C-ktitQq!D(zIO^jW`ToueZhn5SM%{FG z)WsOj&riQzr$#*k#01b{aA*{&rQH=g9mcn=erjU+$a7r%pcxezA}p|_;PoQxM(yol zat2CV2kP0B145D|mEvw}$2}M=eOMiy%0KWVlP|GHaqH=cZn4rBgtX!_#djXA`XA10dXyu9ER}oRdn|tMkcl-mRzPMy5`+PR&|*+?8WM8qOoD*e5gl?&pldi&$M5t! z?wLT)QKUz5F-Y%JyJCH^?1?f6cv=yUBZ(s*5Gn&iQ*!4|kdXv|kZ6HZHAo1AW)7j! zA<-;>L%Abdg^Gm?$)6)O<@`)3(V6kw3Ojw9k2%AJj|IWSLbHs>AF-_SkBi7ZJh=bR zPnN@Hc zcwKe-`*25joYU}Pmi=odQ>)7DkE8vO_@**8fdJX2NR|Q!wwSehys3cJ(=!<2=qEr@ z)jG90;!=lTzLTeXv$FrFDkhc<^Oh7-r=+E{ePo+>U-)geNimJp9$*BTbzqCRQXVsp>t_TsWmjVX|mGFGv`y-J*H@=2Vq&?(?t!CrV=xZWr@v z6V*2=hlo|0&>LJ2XUBn}*>%O%9J~YTJ|4+dLFO@7cskaig_!Rvfapv1#%(>hI`Z0V zj%@2N;ryoK)r|c0_f=AJ{#Pc}?B8M5p1!kw zOyIfS(kD((iSYW^@jL#OhIwJOjb!HJ;CNiuD=D3^ER%Z?HwGeYHowaqUwmc83`*Sk z?d-a$K<+3twqH3vHGi(E>YnZSS^KcY$X3+vM1xNK*;xjXEdU5$@E9~&1&KlPPLM1# z2E^c975fi{M)RlzJd%#V((ycaO~?=g(r=Mut69=1;SG&2G7Vj zjmkl?L1sh-6aZla00^DeFrCLx@bEYrjJHU+j@-|WXsoUS-#+X z_3Z4>wa%6K-&fq*!&{ZtCQQGdEN=;)gVk-QFYQbaw(T1SPjNdgUvzwjVXB^micZt0>{6Sz_Az|-xLgR&n$=9fn-y;8vjqDbL4_f!1 zO>K>oOz5t9nCTxsxbMgEzdg0e)j#}P_Evl+<%{3FIoMwvB0s2*v+K09EFaZd6nN!! zztPp_uHk0U3z@I$f17sBK2NEYV7SOh-P%o%AIT{9%l7u`s+|Uu6F4YlhuY764Z}g$ z&`Ej@&oZdC!BN9$L)}g)9nG^hl7*-7oIp{9*Mdch|GX*&5ctos(RmBbo*aGmYDh5# z4PoRcH*EuAZvY`)P^lP5?Lm3#H2beY|DOSH=pEw2tKq!e41oM=V z0#ozH9&V^~8@%IjGL~gar=xr>vWjK%3D1c+yb|g^Pf3KEf*zJ*DFS95bRO8oGw5_4 zZ-wT$A%kUU;zu^M{NqjI6-+GLA3AS}mSmHQsgXlADHuQ&gzEuLeo{M>$6!Ml=4c4- zIE_yNdrPXKPw!!NUt(j1eu~s{gt{rnpbybG-ZO_=_s9g|aCeHW7UQ~(s*L^*oK4uN zLId!+m8H_fmYip|<5PL{F3`g+hs5V`c}nF?9u4|WJv79BRT~X}&ybM}^M?fu-PS$Y+d4m5SmRUdKFHgiId=6#d^AR4(!Q zc*MvN07M~lA-qMTc@6&t=7L2|s3~qiz?|mWVNA~tl{^476StUQOI%V%CQm^|=%IIx zi02J{q#|$v7#jGOWG;3}JamQw?lXHAEciq=Jn+LguGFx=EfrTKs_6^8(0(EQRME%d z#xwGU2Q)RrT$1jj9IzrxLOC2M0GeX%7W(k=cAU+#RO6H{V+?aM>J|6w)z=}DbpR)dfxMTGuKq=b6 z$Flv)8Flag-l-W@H6Gea4nMqU`aT(nM93#lv2iCei^KHm?M{{he`jY`6-;SgDq)NZ zhOW57MN!TAAqO~0GwC9YNdxt%^1F-X?cvR$Injv~KmE>gz%2uR7CPvr(KiKU%fG-x z&3NS*iiMt^q_uFsTtu$-BdCUUUwH=W}QX3yrvcl4Bl`$2ynyY~4I4_Q>Rv{m6^^ zTk~S2Y#E&_ynO)bix1}15Z-L>r|KD99^18PagGLQ^u>5cBN&V1({cs1x`+O37CFyT{Wo`T!g6&92ZgO!`>ng zcw8@q92HO;oZw0(N72YOx~0TI!D|98H}abZl)*9>30v1mi11>j{O~9gEQ*ncZ%$A< z&@z(SaX!cKu#)59NjN9sMgr zul@z|tV7-wX?W&6wjW_=!)@gJpq)9CXbDZ^eksz8ZFfytpS2~vs(ylW{#4G*x=DH3PR0tIMCj zg+55bDD~0JP(SAR8j;DZY%!Ve?I>73^=#ETR(rTx&GA-3=Yi??nA+Txsu#&FE}?z? zY9DqM($bP91~a7oacFm#O79Y2T+bV_JeYxM$-N!{h@AD378Z=b4^YnZozY;;Y$v0% znq{>5Jmq@Ld(LP}qdR~JMdnXVw2m^3ApIQkP9?30dsM)bQzd!TI$F;?*&nJUpoC>8 z{HoP_HK5g_S#d11g#Ta`>s-t1Gl=$X;tW5b+la8g>X4d#Cvc;6e(lt81{aGP|Sj;@>Z#=0A9TlH#= zj)1xo5b>ytj(==ZR|wK5b}}1w=u{eJ^Cvn&_3qNoXjTepiQ~$73@fP+q@#w>&>FJG zYZ%v*p14+cqR*wm*Q%86A;$LaQcaGkJD(}5Vk(v+jW+Re@=EmURCMEN~+I8Q`<&UQ|lUHz|Enq4Gmv2m9%T1xoA{+ z=X-zRmZ#j%pXodiDrsKL3_Kp?H6&KX9~qtVf|yg}31SrCAo;)A^{kwx!yYaNaMBZw zJGtn>y_@ZJ*5mB&+qKl%p^{g{FD%qy8zOtC2~myE>%~6O*M<1KpML68y}!ObGo6FV z`p$&g4dq^}t7DZ0MCys_pwA<-)aFz2@bCcRk%GI0$xxgOH(pT;F@;-*J?c3PQUb{;0X`gJ}FDI#=1rHD9F=)uNqbj$Qn@o&mGg=Vq{WD@FARF09-R8C~Dm5C8 zvA2L16D!Hl7=TVdmKl^|3v85xu=Ve+zmz948+-B^bA7BkQIdAFNVc%A}dOuZ0J#xm0cx*~-n`u-Uxw?s|eOoKg1IfNHiH z(R|1V2aS$wnP#18=$`9Gn(&^V{8JDB=Eo@(dRhlBd^;+V){c$WTpDM=K=(csh99R)bq^aa6_Pfc=5XDjLPODCm{|v92+tQ{0KzScxyd zatRoA1h+n;f2)vzWtwR6IiVFS!F}do!8@|96|ZK8>KgkX7L3rrr#Ig#&WzF`zKwQb zzNX3JW~Jc}Z5cknkJT}h@Few(qDuJ$#c4Z96NBvQ6^XiOg;taHJqZV8UDUm)Q+Vj6 zpnOE3zp(4D^x^o8>$)33A?+dF1~!0oh==(c{$Z5~GFI_d-%Yc9W#2)lCJ!Cq)pbku zG`L?&%DzEMr{f|j?Uu7Ysk>%#VhGM*nueyH^D4uC-zVRPJJgqc*TD~yLM zEhvOyXB0+L*QX2KY#JYPSj?s%yHEz_jzlWX@E(g|H_)mdJuVa??Mw?amF)tkE=0N& z?Y+@s58sYb@GJr^y{G$`92$H$jY|&n zcEXX7FKro`nlKCjneQznm5a%Ba>eZILM!9=E6K2C@Cx8d zDuBQx@P8N{fRt-$RxEny9l;@4*Pv7=9T#xTqC*pWA=b>QFFyu>6I=PDEKa!CuqJCi z4MP&gQZ>M5_3|ov%5BuuAw6t=mglJdtj!tCG#Q(oTT0kLf{)On7RF6qIcDMhVr)y< z`x{!xC;x?h|15k-#`+1}d`IO(Ph>?xLYdJ>e{?WhHSpTPu^FX)_L14g?iiyEoXp_A zw>cTsu6CWbEG^*q!LJuETFJACOV*sa_0DhW!Z{O0m_)e>?GHOYeNC%0bZGXg)+Xo4S`t07 z>h6Z0-JI0eyH-~+yOgz63U-#^f!C4eI~1$yFK&8gqz3fqtc15r%d0&;Bfs@%bZjHa zZlcHUgR*yuR}Z6_l&GklwmvUGCJnN!YmjW%b8I`3&h4#zAR8AfWdB$ zw|0-VEl_%PyE6NB96wqOlz?Cpa5E$BIT;r5*df_-CrM`&Kb|5Vz`X;5aPRnMg{HiW z3`p-tx)w4O8GQVcT`bFF)(suC4vk~5TLu0~wE?ir>hFUoQjym#x%;8%?+i;{ z|LD$qc2tw$dY*BtS#)Ud(P`?Qm&}cFW33m9;rapAvmk#FR*qbw`FEPqcRi^VRfC85 z*iN+U0juvyTX~^3J_WBxkS=&%v3WE?H}FWA&$e9tT)~PhQF@E!8%G~$du*U+8#A4y zlx)Vb+l*`ScYn0#0DAZXJ(@FEJ=-=u97jVBymXW;xw%aKRBr#;-&Jn+Z|l1Cua2sN ze_$d%{XAEl)M9;o9-dE%O&Xh)J*^gc0o zoOXUUC`|d=())H&zwE!`CD(;Y!bV&&fUxB6G)%k)J2PU&H4$9da7p&{o-PhMS2 z9p7IRkmOR}%KVs{>?d;bYxdIJsnM5-#mAESUj2A__TW<8<1?FyjC5X6-x8i-;`q)7Dr?!%uc>N} z{gST<|0d^8ij#&tXFm_)yE|mryq58**4^Mn?B>#+_BRLD4IeVY^hW&rBxas|oO7g& zz13V+n91HK$=R_tytt4p`%`z8FFq_-KBBkYqFrBCLASB`#(_@Po5xZ{0SV%;<(tX+ z$0O?A^q(!0Lz(3pN_WU1$A(VG?Y_O~oA61#Ve#{=p{#=y0I8LI9;NRvHkxDe(;_0` z&d*^RgBxwPVslsCEc~gsGG9^gC5v7fv?%M-YdBcsFzxL4@~xVO--|5|2g-g+Br<>ZP*XrHD<13SV;Lj0!2=#H-ApjSbGZk zIV*FA^Lh06{c}p4Sz@1L(bUhDi(Og@DiigAefo2PYh#TU%I`&ntOTlGr+i7bnwifo z_^hwb4opms`1}6X80Sl-tkL|H?F6lvgg-MfA>|P%-D_^|p}lV(jMcrgEM84aITE;2 zr98cI(r)AD$*a156e#LnU%~v*0V_{_bvLDsN7!_%d^R08Q>k*-+vveeU9W74zaRF; z6-?8VL+||gg&*~qGv@|_8iwT5%%>#Jcm8?{n)1GsO*6|4k7=Abxd3YI36;D`B{`5E z{yy8fcW0^peQij=9lx|$C&e9S>2-N_TBShfi^xBv16wN-ZBFEhQbYamP+MrB6*;ed zw+l|_eq~)^-C4K7= z9;%oENJ|`u5EzGarzbNB&{>mu@{-yHOLe zNE{Wpy=}34ZAk&wP%YWWT^(H*s^3i}w{+nKdGR~<@_NSZ-gzY-y}rBZVx_MbeyE?! z-<36~X``Kup}3LkYx7ATn|2qENq6W0Hgy|y*1xU@t11pFZ9;#wU#PX~s~n(B^M8@+ zovhWR-ub|2xFcefyQp6tZuSi0T-T`{*?Q)$@UMD0nHaGttMm_X#Br3VCsiA!w;8EV?9)l`d4Iv;LsjZy&oQ5Z z+Z}J?3WE~mTHj=DoE5j|J|SPdRT8uw(=cawju!ZwI)R@S$C=%^$JFaslwRK2J^29R zL{2@ok}(zT*V#O~`cu8R`R(YEUTVYHm?IK%mkV7Z*kT=dyBa$`zZ2T4O@-6O5XC-}IY(JA>(>tq-BKO7E#-ka;c$G73DZV<6O{;9XdgWGX- zaD@B{oaDs0_WqXo*LD7R=k_AQ7sJ8I70a23c&j(Ro_faCYu~lrI@Uht_&7V*>igzr z^JO_QJ)-r)z1i)B56qR!h2OQ$>t_-Znv8S4z8sq>Gdu6?=^>z#_xRGF=ERH6zkQHW42FT2Uo*G?8!wMiUHH>vEN96|i z=sIqVKQW(>WsW;8 ztSE($r-D|KS*;5DFYe+7vJ)aowLfQRluTtBCf(IIb+6xXffF|OCZH{RdvV>Y_MX&UVCkJrE8e?I ztn-u7#@!W6$|E#|xbu#s>?- z%)a?7*U*}4PoLPEYYpuFB)>W?ZsEkEr0y2?@5;}JtWL^3p1JVaImF}V?32-s3#0cU zT8(q-!^J%_b1D9%iF>7ccYGI@6jT~?6Th!qAq^Os?46gK{E~D#sGXMf>8$l@6k|zP kJ2)_b(!0{*?wd4O({%21Z7AR0=CNceaR^_T$kxUG0}iR-xc~qF diff --git a/data/us_flights.rda b/data/us_flights.rda index d5cc1e8ffbf9bc089bb1467cf6b5590b1e9c137f..caab26884c6d86d2cd2d0c21c6399b4ccaaf85c8 100644 GIT binary patch literal 27502 zcma&Odmxiv{5ak{l(OZnn=x%;?w3?5Y}imX8#A{MX+vd^ODaNhEo*a68PjZYH!M;K zxh68pE!|K;Ar#g7`}F<(^Zn=h`Td@4kMo@8ob!5}*X=pa>-9RU$bsetmTp?kZtGV- zZ#Ozzzy9C<*YPg3{^^GQ;n=e2#)dnnt#Vs8Yye9Cjt84g7OG<(Y}l|-6OF%NM=!)~ zklBy`y8+z#zYh8TU48(ZW&W?)|5p*LNl(_qZP>Iy8Uo+KK;4mIYu-tMV3YJWIBeMR z^#9fCUp=ps>$tJR8kAHLtFN`$`L(sMXdi-t-;mhML=<2&JTMTx`q?P{UOHFMJcm_^VKaS|5Jt0 zh<#hzm24Oax4vXi2u~;6RhxqS9GFt(g!bA{Y z%LoLEfPpcRAvZ(_D4d8waZns7W-&YE3l0t*;gWcAI20F~423yqvOr**JR0nTgX7>3 zC>-ibz;fC6xMZje9?RzJw!jey7)>rXdG}|c6&M`Hk#}MzL*Y*Bcn@|d1}BX}IR)d{ z3RD=6h?B;~bI>vzLK223W5DHMxl9fk1OlNrR>|_%WElGvr4R?UA=@ip!MI2?hVV=K znLL8Tl44<1N@I8_TR6V3CLRwavU0FvB_v^!Gogx#dayt8bR1C;O3=ljw#2m-61(Mz zs0tVZM=2X910Ynm|0<}Sf-2fMKrV9((b-IKN?piA&{5?$bax6kxkjGo3rq4L zvx#B$?fp~};wKcC&f`H!z3rNX=;dThHUW_U7See{Y$?`g;Am@#B6JG|Ap>S$?66KO z*)0`!^c2g8zK|p=s+g$M1!jm7b|?_7TRaR0nymK>@RIceFgB4G2`=0z!_C})LGy{I z;^bl`RtqSK+XdXr0z4m|F_Vc2Z6N~E%p;s3mHJL97D_|F zj8aV|M;cX#mVs(|Ft{lM+`4vmC)NehhG!>ynB)_k*mxZhTb+X9Xe*G@v96L6NT+^p z63dH%t(3<3@)Qs_6j8S`2|u*83mTfdLmsGu`b-&Q9Ds%bb((d2+xDfalZHCF^X2}7B50(4Fh3dhn#R*P+5cH1QHAOkYVQ%85OaoSxl z9tJIs#n(2)XMpeUaBwt6m$*nETeY%Vxk(;Orz#>5ODBfng=@)(bi)Mt(<*2ejNP(G zXEK^als3LA(gOwX9psyW=kb_0BEV}Hn}8+Iz(ijcCod=X4&1bdM2 zH-X=p(CH!irEL$st4e1(r0&!d*lgv4h+ukZ%IS+NS4b+xqe<8X7#b|zqtH53sX@?^GsU$|%oScKvxfN?OZLEdeiH!J}I zX49EvSib8ZaivQXWS^x3m^3(ckVRtSdC6cHyDepF5>ZYt2QKqz613 zfGVisCZ@MPEM;%P!{3UEn!IEfvmDE@@1BtlfDCZiE;>$CT3aWHN{k}a6fpC6@)`&K zbWDL%g(^nlp|%P!S0&iR$&@yHWfKGc6`TxC<`4>>dkBd8pwn+s!mI5bGHwPXQ;2Tlhc!$jA*a&okFW9yD3(i&k84}$iH03ZObd< zDI4d49*<0)q}XhByM_UAvPa0P9`pxy<1!*+b%ZX^n9SDpmJX3`J9@ImL8i)3zrs*Y zPtQ=#P~U*t+NEEuJgl6f?4Yges;^&u*5j#&rrKJaBLz7~Y5<=auhI%$ zJA-znxH3NR+S|J$b=sXZ)TEPQb=m|7Ue}Ql%9L*kDceopa6Z6DS)zv|Ws*~abW&!q z)sjpD!>kLkhV}%MVk323G@IOrh1n zodZXOib6H5VH&05Rj_KAP4wn~B8)I7TO~M;$m{VJx~T|iRNdtDsU(Yam$e^+c(Pcsob`kT2b%Nb(7*p;VMD5u$mrh}PA+O#*ov_>tceY9qz{msi@a?g<8W@HW3BrVUa zmTF+?P@*M9lLd}0y&3s0eJDbeG8-yo#3;8c4_7isF*H$8@HY*1xf?hbG~8=rb9c`U zc#UbVqPt}srFOV8pgJ$Wk0K7v5u_2_O{tjd4iU3m>~Cm1!YlJP5UTP_l!g_Bs%2W@ zy1=?1ot4hC{8s;5K|9}oSCt|3H*p;iI4n=`178P|-jET-BmC}ea!`s(ZTFl@uQ0Hx zuGb~ck=@?n&&|vw4FO6-xEOZNUROFw19at;Lb!b-v1pX$*r3mA~OWl#Dak5zy%NF|qv=+4Uro_AvZqeod zVR5V2_bkQ4b)d^HDAUMoHqA78GwuP8%&P(8Leq$?vs0iTwq%`GFa zra`asg7X!1-D>`6#MCKq9Z4LVhfX7=C?dqfNZufY$1Zc>EfHx%DrGRRc33qnFR1EL ziLgXW5qOcB#OKR{vP|4v+=<*~v4@z{R2Ab#^s}(JNf=dgb`l40i-IV*C8fOJor$Gf zF@-y7iWQR6h{;G?K@mZ$L#Z?&P`F|oubapd@r1l87oMPwT}CAHRr#uAHN>PyUX$3% zD?mt6ErX`0@)d%5oe5$+p#u^S!&kU`F(Ar(BEZWxKm$IJJLqpb>{tV*1$6bA+W z07{vkn;;MwJkrZfgQ}*M)VQVjE)j)PJ>%iD0GqtPp>rl~z39}sn;Ici+yVQL;c8=L zB}7vO*(VTbtel8#}3gp&vPqLiQIH zj323My!_!vy@H@B_s8rDmCPSyTle9YK37%`8Qw1+5y>3cDZ3Qf6<$`u{OVa_`7Ct* zy`hF6oqA+1)f8*KD>?0u{C|ySg(W{y$mxdDKMh$_Azqlx9(pvbz#9xp3+OND@|Y05 z$=xZtt+%{G9FjZ7$F)))$$s>FbTKcy@y5*9wnkiMm}5Aj^s;PtUZCU>SQ4DqJ`8|P zO$*!G&iACbby4eIYTuF2$n&MGm$Kcb8gs|Z=B^lIPYK7eEIX(TY1_@4FYSG3G1QRm zeB+aeYH&^Vi1~#tG`G3Y#G%kJiTcZnBeTJ`wwBN2D-qh?);p9e^o6LsFe#s^Q=M(^ zOmcbtU@*v3ZRJlN0O+Pq3-Kw%AfMQ~v?pZzQF||M-~u@CHbJO>??0=PMTFzZa^gxk3m0 zA2~iRZV?U8JfYK=$CvL>6F(d0EvJ`OQ_OTzihF4aw0+RM+kVwu+^1C&G+{iPuEq-* z_+oNTYhi%G8MvG$tAr8PjszjiFk*efds=h0=>Fi*>60_}w*Bh9Fuv_i-$hTQjr|w* zuKeo1=ppu~m+rWD?>taU9DjCkoGJ)1%xA4a+d@6gL6PALZJ~LYjmeV@y;oE^$r5kc zSU_X)^x?9wX_Gc;L)l0AWJ5VEdCIihk4M{k&Y$)dSiX7LO zlfq0|XaMCBHqlEX+fzF?T{Z8s2EpdzQBMg@IpYX#N#wOvgn< zv$CRymqGb1Hxow=>#9*c8y?xEhJ5A!Xh00uvjD?r_6S##g7Sd%3d>uLdm_R;~f-pSzFn( zj`QZSSESQ;gH+*IU1RpNsJyh-hawULl{7{Pn|!Gvr!7SdUDzAG&=um5+4x|*A>BPx5^FyB?l5>7WQI)q z@`Q3^=sCPGNS9ijH)Q%eV(9F-hO}EMm4iWMX^KT5<8^xzTgnnAD?c82v2nMTx&r{#`Xj`=|Qys?3q{OW7j`Q~UX<*scr8s`J*K4Yi{CX;G+Fa;wI@jn}IeJ(S9N*9*drQ5tqnnm7R402A^Wa@o-L5$DuGoS% zXM=N8PvO&UsXsj@lHK;)vi59<`_4~wie~A44~A&D+btfvId`~c+B|>C&{Wy2JhbeV zx(h$G3SsQZ3NSCbEX?VutW(P0rm$+z;lbM5-cpX821t8YHj#E!9Hbj$lDEr<7@6N8 zYzrP^ihAm3p4n4%<$jjslQoi&S-?szdHSV}{Zmx>bcJX(-MyfPCNZ7Ai5?P4hD-)B zl@xgWwEP$5#qS$Zcj|_w<+NTLL0%zD9WIZyDF19SoK|!8v-zo;%$LD#f>VAECPi42 zvG&GXNyuK8yZIl@y7MIpLx5VVKj4QHWse<^s0*e|=Y_qru|x0{+^ciK!OMGXqN84g zj43QiLH1btg(8G;X&u3KbCy}h;cvpz(Tu!+p&=jzuPtI4F zyirA7QAsN2d)5Zlj=Nh3I);E}>^49>>&wDJ+o@5nskN$eOP*rm(WyV>X8{ZKwP4)5 zEZp?hOPccD#Kt?VBXz;FD2YGk$tAUVcvD!}IrUmM7{*{M>%oYfnM!W_G%-glb^f8~b^hUD$?i~8!mN+|Mq5(N7Uduukgo;^{w zj~f5^@UELazRA`3L(T<~GeP=Qre$GM+g0j_841rXjkZ5e9cxJOkXW|Ta(b_rWKWrl z?y~l=c@FRMzP#7fvgfjbU!eJHyHkG0nnj8vJ= zxxDYO&mA)_tH~V^C84@-5jG{S;m}Fb91j-)W2UCWwYbYa&Tr^)x8Y;`9m$J@Uas=b zk2LGUODscA2?UVeje5PRf3v)^S@q4C_)fIL!7L+U+C|4woKBehaViXx7bLURtIHri z>i8Yhfh}Dup0z`bGz)SMm4inHo=?ose+eB-J6d+KQ%i1yg;2993A}4}o}+e<=RwxQ zs{XBm5B6@#Iowp9GPC&ff%RD>%@6kqtWT>ybg7Li^*-m?elPA+rQ?yvOpI2}*_KDI z*lw4C?=!vKl)0XhQkm)IA`SfA9Z*{y)?t0ul3U0 z^g|&q~#F!GMPTISBXdpE6M4ZrlLYWUtc>k_|$D@B;~UE&WEnUDEz(SV(ngM(^P-w z&#EP|JTGNU(~YHPR9U52junFHdjS;JObe@8-Q)nbQ-ge|Al_n~HWjb>>LHo%6??DY z%KQSLqeYj$)Sq_wd?anOC#1k{)J)Idg}oaR)$}5pnRm{2X+O35Tq7}6I&9~z-+YN# z$glb+PK>OCBVH*Bp!9dN8Jas3T01LmHL@x*q|5V0E@3$&ry)gy5*M+s>>!7dtzej= z$~g;`8Ngsq8&gBg?^&ODC22$8u=l)O>saybuU@^95KP|R_Enjj)9Ql_P1+x&e(GNN3W#1HF(K&eK3wg%ENK)%1ntdzmi+Jg5nf>|j%c9ho!@5P5k%KtB z(Eg~Y=&>5Fvc~Al?V_NZzP(cx-I{fh^9S?)ZE@5E@YU#m7r@pduMR*QHZg3``>?%dl&Vdx47h~$mOcV16hBJv}4sE>5 z$U+F!3wvR27wIVoxq*QF4gRa`li*7EboQwHqFu;4>Ho&MX?{tvW217jZ?< zeOA+=^N(B7zV>*ZJN6<|G+8_1ptBn~X*hT+;fLT-fmMY{#$5g(dX=ltV}Adp@}h#t z23__v{o~QP6z-z(eTOR*>hLq+GW$DuGlw0w<0MBFR(-!5GyRn|Y2l8x#leF{C@K}C zo)|+fYm>o%M`=rK=bBuc+oH^p`@>WO^_Bh>H{pZ6k!3}7=M9hgX*~8nlFr5#aYq`0 zB##0rEV@T4I22D?kv$~62t8HnS};i)H`?!8VURgvb|JKV)p1*Y!hGsOM*EDx+1nL*4vCWv_+f?xL4 zwSSknJ9G(y`PdW>)eFITS;8Y=tjYbBTNn`)A!vDAuexSa9ONLAQD}R|3ZI1e+6uv; z5BSD2R$?`qI|*zUt1}WNP1F&42hg#n#gV7)palMig$xjw=Z0+-OWXHC(Yxy+lkDeO z1-MTca*X8T@8vLK+_g}U(f|%lY9YjQk3&vBQYWMj#wE1Y-ENd66pgDN(9CkOhp*lo zFB+>EyjpCOSNr(G!N3N+$B4n0lW5SP2`Jp&09N#tQOT@JZF8MK?nPzPpJU6hal4ER zP*sAfZSmP-CZ6tZA51{^QNQ)ud>I6RZ01m57&FYJi5np*DRMl!)oGPO3O^B0*|kO9 z<6p@2?{TRz8gWcyaKSIx54h*{-(N8^uKtutW?fdWs>r2VB`CY@0If-H1ruRKFCnQR zn-}a-?01zWxbS^4zW3<0Q;Wx&1;1LucHIqO8}EZ2NT=T|cq|4BW*(XT1c3#f zc=-b~asn%TaB*vfknD>9r zbzSRRE$`JE(TQ35{W%M|_y#=Zt@EcgX82bHif+C5cwtr{diCkmg!0fAm>;fh+lbzY zfAsQpN}rz20Ljr0WHP0;Aekr#M14(P4oU!+As4Iv`?lr%HAT6)jlOFca`6_?Es@gU z7qIglL)-RjQI}i#q>^T;WbU*1zc&Gga%?ZIRFHm)KgX%}WYbYL)_uQw%bJ$q+CFNQ zC2%QqkQBnfY9pvrI>4(KwMFNy4*dSI<}vCO9vh6adCE&Q4ejv(+5e#>GJhN5=v z6em|yRvf?gS(jXes>xOwwpUgj7C3N}vVp)shd-hRJ4M@LLPPpa`4R(QN&dPzem*|- zCJ70#eml&mn_AEKc;SdEctLVyRYjKA+Xt5v(SE0QZ!@~ZECGRFrLdS3Eq(C2LtLW_6gqC>sCD_J%4)rp6}x8?W- zOm8wyM|($)j6STH=Vt#kOV3ELB1?}eiut53ASdZV^Q%JFO_*Ke)|QvVcyCfSQKz=5 z+}n>Ny4sJ6O-r!evqRm_g2?Yymgt-%i||X}6cu|QcqtGkb21`hsNl7!EkNLczj>1j zLd5Ly^{+H2^YSLUvMVcg*-`eebZ9}aj8<|F@Z?^8cBS{ zqqT`IBe(#MCH7$TDvd~4EK)nW)b|Wd`QC@rS^*|W=e(Cb7Y9wUN`@j-b~gV`zaL}srOPpr&nVtnFe#Dpec zHB`k4P-QUNeg>uogh%CZpk)L%nUl{vV-Kg$`I}Ny5j?CD&q}O71Tz4I>cu2!cS2)h zTx>kIx!b6Vd$w4o)73XQ#ex3n-B@pPc=2aKcK-es@5q@am##r}Klt1n%8mEI%F$5= zTW2SxVjA-o?!6vtto9ssZnSv`8uiy$y}#n-mK`Uju>_7sX{;UC>ZnviT-~MmCgftw zVlnCnGVd|**B80iHNE$g$7(+(@uYzt^gTDSFTYh_IVP&``g3cJ`$tj{!~Six7bREr zZT$NaMT;mOaXq<`5nl+yDlGn38AX|v{m1H{A1^cbrT4qZHDY~rq}h7p1o{T3GJEwl zKcGRc8n%AujZRE*Qq0X+9D7}X5JQ`FfvrtPv2Irqhq<%N@SsVz^xk0(wI-m%YPKYPFAj4Kj3>Wrz88^#fb?&nv|t+TV3@Ljo{Cg~zu zN1M+@6-PqMb#q=H3CYPWO1|VUMOjL|WdCjcJo|Y5ov{6#L&vhZt}m~ar?dpJ@ee+CzGk}`bxweXru z*d>QltNkOJ6k^+PrWCY)(6tL0WgcUb;}YEt-0o+Zx0|D1Z2TsFjd5aQPt3)EZOZci+q^dcuOY{3_FYRYeA70AJ6BS1% zd5z0nIW{Rzh!g3*CO8v~pBh)MF$+p#sm)g=laI731zp*G5%^wy^<^)!3?)?uhousPZ_;Auy|ecmUsRA`mC?wlhT)g zwP}zv4UEzs?ca0TAI8KD9;vs&NyQbE9pn|~*?&HR>)5hs&!W#> zm$0Q1u8b9r)wdP3g2k;#9Q?rmsd#(4V&hWmy1`r~cS8S|KksqrT)uImeTGYm z;df^b_E5)FqXJrbbU~JRpZOfmW7UGVclFqGmIOAHzf^GXo}11tK~LB@7uu`8FGBvF zTbqdgFPaqZPWo>$_@=8zd)S0&Z(~G-y4NeS(~deV$Cv}vtxZ#WxBMa4km=-YacIG_ z-xoRq)fe!$eoeeuwtaO~cr<-X+s+e;J-F_2xpCV!k1NPKQ9-nJk~esB1RVuth7Lw; z+AThqp{WQ(F{;p~nejNdeOWjhCIoHe!RR-Cfe3N*^~=x+j^8giToM5yuNhCiLj-}r z5KV?AkY1Ajj%Q@e9Q>4ovtE?bfkG8x7WPxx;-6$lLtJ607!2WlG7(5AL1hh`%Nppr zreT@c6JbA_{-ewP_>Ni9{p4Z&$85=mvY7SJ$LH&bD=8`$+I-GLFOx4sHIFt;di^fC z*dLPraPYCO*Pm~F*M#53l%7N{6Cyf9GwI_W8d7F>3trmkD+u2AvAb=^h=<9`A>W%8 zh@=sgeT7*maDAnh(DrJ*t`VO>S&j4Pymkc7VPT80$q^>xF`CoJF8A74}Lu&2S7wDY9B6x(D*Y9<_qnFhP z)rJGD?cM&2J|FJ@Ly5$QLFS25b;W3vN9RDhk!t@jL77hhVwF zp0+%#iV#_?wJfHnDCzS=F$s||p{YAE<^A;$na0IJn{D_d2xqc~Ql(c}QLS%KrhAS% zqHEzx%AV*cMYmSpKX%u|cgPP;Swd3yN`?}}Y6-B!*W;c!Lue6ErqDgmsVJhC(ay=9 zpdmLRG0OI8M-Pyg;tq9dN0;EIz^N%fiWwhICbNeq$w(hBub@n(J;XsrtTfQt%m)@U zoR~}R3?(^^)sBixLD7>h%=nCsrf-z6TtYT$~$Ue{yv7YE~3`F~xKJ3D%4IE&b0#4%WEm*441Y*$?^7 z$(QUBR=NI)&6~F98Kk<9m)+lo7k~5^LGJ!JKhZhHW!>60r%Ys4XX3luP;rW^K&sc* z+#U;#QTkY=<4DtQ#Q|LDwOu!PD-!wK8}oxUjxVSl?aMv^?yDU&#i~Dc#}&O-&J(7; zW8>E{-@}8Gk8iqt;~Vkn_9yo5khm`{y13CSbDjE-D%YU7A*gLg_~*I)(MDg-Mn~Pj zxtC+_^0qyz43hnR8;HUl?y`Nk$kf}Y+snOo=tGO~>{ni0gTtPeLHk^`1_Y1Yg1fwL zz`EbOScQBBS}w18b$#>FJ^rkk?U%|EgENPBJ_`KrRME!^`7Z{>8qXE~+*^-{@XBN* zjaBYA;&$>rT0UGiuL}C6anHurL|Rqwz}>2FSe9?NuBwNNd_y{%k{J8MN_}_PmGjpY zm);nUAF{|C)=x4S{ti>9|LxOVajUv1GPng^r+kZiUw)v~~*6csF+^Y|N-(QIQmcAbTW$a1e+h%v3{f=3WyFs<@p^}9>+4hsKWQ&YjIH!CrH-eR2y*Mqa(4FRl_lz_RG!4;{W2 zsDF@;xF8q1`Ea0vPpY4kycFrVvXxary}R=_?Y*B2znMKe^6`XA!tGFp+>A3;`hrDu z-*c{7&-LjS$N59{ioT@LyvzPl1B(Z)tvniEHF)iRY)W3^qv?UM&rdeaoZ7soHkHSi zHV;W*6JK3cgnd@h(|AC;xBcX2eID?%|GaN2{UOosTP1KWD;^!XntK(#jitjP+hY}qlhF!4^py6m;%~oi zi}*@oFNz)3r5YdNilI-mma=Y1>;5e8(-nS*@l3}l(u0|faX!VagkY6f*2PT3f)~P~ zunygaFKUiP6D|_%;d8PsEuNjWp&hd_M~VA@!aRQd)JE5bwxL%JZJ0l`*#~*kM>ekg z*$c?7qSly?gA0GZtt$7&{3(Z|O;?>-FqCLpa!$_uOgZx7?TJOC-9o?gAe)Ft%is}x zKbUt?t0&_53-^=Rw;3}-C&CuIoU6SfD}=aDO}1wS!qT7BdfB4&inI2uP*mqBOYAQI zx$3~5xtZb^-HrYGMR$-l=N9r?}HkI7{B?AIh#%A5^XFzQZ}q zSD$`xA*cS)bE68!*~5-J)rRt8n%2VW1nb$W%b+#c{zqokm(%(~<y& z0`(vdm!YL~X8Lv=Fd67}jIWZ`gWJRVisrO{Rt_0lSDpaS%Md zg`Vt^ya$&B(&Q{+d8O6C{UAaNo!-HOlz;h_(nH0;bpv5pZ{@1=y7%9MS(QRzBkD4+ z2NX>N6Lsw~)2gi(jDUjeP>jgdLKFwCVSq56$i%^z+2EvTb*2ab0Xe%hTr(^U$zpA` zP=ByAeEyn-W6jl}ZPWz^y=Ul(Hht*NGTu@{*xcxXqtZJIh-b=zV#O1Mma8V+_f?1$ zm%A=0HgM9|vkh>UL1*72-!%N^5s)lfy*3yG59r4yhf4`>8&Jha`(A(1X%kn=Dc%zH zm;|GDaOt%1(rKIjeU|HwxJAG8R(%|D@XBqxTvo*FewGiBj|Otgao`LD490-Mth|d0 z(b$7UVkinibh4BU!)oG)mW0Q6ClpQ{!KHg)4<(L?07lK-T^xEdHd{g;-A#rh7DiQTq zWqmsSr{2Y1+6TVI)Oclvco}tsf4lcxE7Lg2V$&8W zwP|Fr=jc|egyly;eiPSsuAlikk1#)TG9m<#_M^$vv-4Wd(z>pL``M}TB;(MhvC-kW z?;FETZ9en%m>Ii>0CV?`QN2-oUA(FN&5?GIU*w$Xqw;tur^98*pO*$s$GmIKLN@!m zZhPMpa}$|ug@a$c61ezguWfCxx60auju`j9r+r5a``O?vu-!ku982(j+d$QTynL(# z8aGK>BBj}G+GHV$(q)?;uh74J>3dURZ$jj^!U5a5^t|FbcIsB~SPYnPn()b>dDvAJ zSda4}c3B5~^_m+im$gov7HG?zJQ+*Bq9qxgF5(uvqehx3;$jpB#=t;fP#6{m zw}L=nni-lLIA$Y^3BlNbt?(E)21tNrsKc=+HXf=8{G_39loc42gasNi*-!|ci2-9l zaJdX1VVrJ-;%wmnN({&1?U)=e6RKGZf#B&_V9N!Q17j95=}C4V_!c1h8>OhCzt{~?C|hYxV47=?k`u?Z+FklhX>Dx;v93;=7ncs7u< zjAt?!c#JxHqb4v(kY?e(nj|ni4j*3(!JvwP7J#RhX$fOm;Vm^W@PD35O=dg*7Yesy z;5I=ut^N;Xte~1GfDmzjLl_Q{V%Wt&@JUcPpaw`UmjYlwVR67sW+5CPf)t$@4+bA( zVx5>+fS`ra|FB2`W+?^W%uE6>gy4ZIyCmQ*3TtNt5cgkCI|c@eu>xkV$-u$r^ndIE z0mEqm6b1v{Wu`iS9X-Anr1=0Kc|43M4+UU>nGp4V;{v1sGVV1WFbh!`psgS{Mj8r& z0}bqwuwcePU>f*%4hRlF1Goo3vy#S0W96V=C>#Sv0k{DSD~@L{tZ*<(O*RH-06>Ms zGXNz1q0G@N#;60s00l4;paeiA6xIpFVBP_^E{*+11v`cl)+z~vfB~YwECld{1FQx! zQBW8-1DJ$5{2zgCdE8;2WM;sDE-)+n0W#AH1`LV;fiu8NX-yQ?1q21;5d!@uPq9D` z5d0(qKp3=z0|J5V01<=20Q8Gv;p%V{FwUlbtpFIDO;8x)B(Ti`g(`;NApnXHJo65- zQ0^ao@>qzxT^t_Ry#riAAyUAN2Ux%`g73r^1Jk1eqy8V3(ZLxT;U@t$W39j#IaDzq z#=w-5FdzWTK_=Y_uMWOLXRr~SZKs*U4k23qyciQX_sPpGWFJwx)wA(tRS%hlcEuHNn_1VXD^oEdzQ{KI{%M86e z0Wjct92Bm}WLkhv(Ry+FihvNJrBqd@cRV*;Mob@lp1;#~J z`5^NR|JhFo^dZK0v$T)3QT}b^?#u31_W8t)BDhl~Iharn#*lGd{l%xk(dsPv)D!77 z-t7R!vARY$;p(A~_d@Gp;acAs?9Z)wU)o3zr%TAXQfhluRMIq3$Il}QPP_0UYO5wh ztwI<5`dy_>XU<;NF6a*6R^BrB$Xz~fyAH0(unnC&acSa7(%Sj(Z^^vd?}%O37JhGk z8mTRREwp0rnkxS9Huu-e2&z(k9N(=^HhV?pn9=jgeX{SlS&0!;{NNd{XtTsCQp?-9 zI#KaaH^VQ;{W%nu@ptl+P_HWdo7c>^{FAk5(MznD-x+11sQ*a&{vR(GExQN0e|fx8 zeOKFBg~~)@ktph6z+clX#|RegGfXTFh~$OXhkOIcZlRADPc)4bWN?4iW1Nio zL{CZb!;G0ZZb4lIJW_WtHn^|Z9*O+933PktDNG0AlL|d_=Zm=h=z7j|(cBmf{jqqo zL1!C3MCAP3-azW+@@=@4N&oNDK`lRZZL@)t^`F}xE)@)q zo1eyG)4zPlo0Pv|dvTsI*99FvTFe>Oc!V2chQ9p#DgQ1bz#KLLk()apJIR5{;{bd7 zZykw++RVMbXJ-Cw?{;ZWPqfW)XB2BUpO)Wl{U|KGGrPH#U+_Wa(pZAha8T-$JnoA} zqLxCBc3WRRj@jBnc>g3JbIE&_3OP|ajhjC)PLa96}{Z=2b{iGu)jy#3}N7*x{n*qUh7^w z;4ookM<4J=KRRK7Ppcjp?!loyhK}WYzQw=3yjb9#{2xW(nCo;SQMK(r>a}Fa`wGq< zz2(d&ETT92c5>Dn^4g=@7U&8K*Mg;ic>-xo=oQ@uoql|YayyaV7voFjL8qHm8SJGe zh`~th0+20(4S4^XvKY*!yB(2-?ksv?5XZ7W$hf5^eo@u8T|zVdPOpe&5SE}nj;=)- zfB5$uzInAxKwn?^_R$M5JBaR6u&s#Rezr8}M<15o4y5v2@apyO^~H^HJ0dDh3E%2- zTQW^e)Cl|5zfX+47pPWHKJn-g7rtD#e=D%E^4@0~F!p%pi;Iv`aJQhm&`M($1p`h1 zr(Qo7Qu?d!L(JwsNLPz;Ot1Ni!yRx z?;uy|_aXDjm3i78c#L0YQIt=NmHV!nGG+0PrVE~~vl8*pl7q0YPMq|K&so30y!~?r z=!}P!?;fu`oPs_M;GWTVB^UcAIj{0>g6wNI_K%8)w@Ckb$L@=-qx{Pn!^aCKxP1Y&u35v~3b=%{&{mR=y*9NxYG5{-8W0Uqnai`dL1I=yGg&ScGwtJ2 z>9i%HRk0PsgHZJ3%FAC*qO6}r?B0ktwiPe@x}T%B5!JL=QUq5Xw!LvV9SA0zbNG^L zN1r~hu?KcypH-q*vG{h?Rj2|uJ@h~&I~E)=rnGf&@2bkuUb8bZUFJLk^U|nH=VyRc zJT&P&8z!)w)Hk|EwTSea&(3zKs@4exHRC{$)38cwakCY7-g_kiTl5P1Jsib4*C8&_!T4MGE4!DGmVX>KDuI(Q;1ylJ zJ8zG^AP)rP$tYm)7v;lv=6d+nmg##RFFsuk+L3pQcgEjnV#Ud60QLBh(*t^vZ3$-N z@%AeFR~FkJux?u{>^@oaO6~Pmke+4i=i{ngH$Nar-(8LwioJT@bKY~m_wD>Xbj-|) z|2!fvD=~eUH-w9MUy$)-w;hhAHlY3(6)l&3_tul+Wy>=AOYB>1L!*f}nGCRwQ+ofX zfRy#1_SQp1#v+h{V8;TpCOcNQwqEaVxnl2m_gT1%MNxq$l%7=__HAF9)b%4`xz)bE zC*1sg&nZk7ho}uLfZyh|JcsXkup;nVD_91eAe$w}%tD2eOMjr^WEFt<)>y3{bW8!10tXLkFwAFdls2Mtniz)L5${kJlvT&<<47r^g>BbZp-} zVNqAX^8@4NGbdtU&GAWc(e;0Fx_8SW8GzMLqlNt2B@=PLawX+kDOLv+Ej%!kHDLdz z55?9!$waYFd$A!~LCbv}v;RF_@4f@ESpV-r%+sk82ezkoN9n-m(Jy}1>-kp6&3cm; z?ic^;#5@E4`C&UJ0PTz~B(V2af;1U@(Y`X82fk@0*^xwOJmD9gh$Iy?>vOCv6`@I2 zHTO9&7SkH|pLe(**l9-Va(WO@ug9OQssgl|jG3%;(T^pjlGX>;>&`4!h#5q;PWcT=bpY`(c!X|)uYy50rgLKN*OucpQeQBW0 zGEP^AqAPbkbUF^ojp*)!b?ZGoY7yTojD(5)=p=?Bnpr)%zeji4bPPC&s*M!RpYUAl zs&Tuyd_Vh|5Az&!C$P*%O9K%V;B|_FfYJMiRtR<+r1T{4MuLO!2sreQVpfuE2lU63 z+ddz*P-t|!PNn-(^d&XC#d$=f%G9_@G289jsBt2rY*SbkB6B%6SlO$ySVXYL*@)kA za#ezY{Lzy<4#j=7!XDL18q&}0Xhut00slb2!Y!tR#s6X&pdl&PnDE$a+qU7@h^}*N zJceC4a6ZSbRHUe(=EPDn6@&j*N9P{T{oMZF}*{ib`+-y%+zON~H0wQLtRO;ydD3T$H= z&?^NWI_%7BEB)Hf35Xw14?02gFJ~`F_T9N-BJBwhvW>O1P#HKg(mI#XOC8NTHk_En zDmUwXt@HlFv9_QO7xruyqm+)n`{yM>j~*)y^}4@wJd?V1`OU%s3EG(zT}!V#;enUF zQ3Qp;wX&?jJN)Do=y=(~s+UXE0nODAtja)+(CXPH!PV|&`~!4|_)$xAGb|-lGa0$f`<* zs4KN3h@>qJX{ZqnlB~v9)O32nYmGfCHn0GGAp^=P?!8&mtC0zNMGklX>cSmPAjsh7 z_hCjZ;6K$1KezhD4#Vm3Vct!69fEX$d4qKR8Mg-3iAy`QK$`4BM){a0xHr^6fyhyLN+m5@^2X6Xl z`%DJG+Ae{18#=9OBT8_~{<#j0)M`0q_>~V>(pBZ1IdvSdD2krvNScZTPw|Tntm+*x z4ckmm{C8L>=J>gvdI*jd#kL45kU@58XZLE^cK-ODj!TrVD&CVMdL93UkX%K6lOa(o zc^Veblo)?I4^`?)CwPPVw-2b{5^c)VNlXS4aF${jRQ_PesjC@frzMd9b+RWZ;dur8 zy5n^cs+`AQm`QXTIyJ9`HL;GAya@a_Z`IpBoqNUz5J`Ds7+4Wz>GIwqH!nNx=-%Nf z3rnXYGo5UPvm=f?ioF4Gp;V#HiYelI z$kp~nwryg2+|fPwELpglmIW4Ge39RJ)T@#0do2hjf%1@IdLIRyS_E*h^he`! zbf@viD|ZjV$fF0bpvByli_=PF|40gP|42ukFYNy{rH8#SK$s2}9F{Hi>z|N+yYz+r zJL!#CHynj>!GR88{!_(SE1snAMDY)QylQ!m$@r@>W|;l9a1@$W>J;{E8bjKfCYxbpriP40E}keyHsoz+*Ld zPd6T%<}KqI1kjPQ7ixmql~AgFup2<^x3!+PBW6+P1BGi}(vsX=dyg*uxZaP;)`bui`j zl6>B#l8q05~#U3jRy?g;2s!AY9gV_yt z33C004^F)=hSx^=v<`HGrqm1!wPXnlRBj)VN$e+)PUbY}FONe%BPs63Y!`eO9%hvj znBuC@gED5bCn2E8D%t+#umaDlqAPQm=bJv#(zI=|6LPNj&S%N@9UIP~?B5&hA&V~Z zmF}D~LxVH?tqH2wlZgi-U$#;>$z}}4lUGo!os`xhHFJQkX9L9q|6c^-5Ltk7XAX#9 zfKaH415`O87Rm)gHXH#Ch-554{el7q2+lm&Eu|^4YGz`~eKf7+2 z-tEz{O@=(zX}uGBEXnOC+?N@3u&afkQ&8!ZhOQtd@5__H!cKZ@@81;uc8)JTgHsaE z_-*;=rlh#jus3U1yE|a=6;81pCw*gzPWNp`w~Kr=2#4<8`hd}VxN&IkPWyFT6mLlMx}ODymmqD*x6ZVQfZ-*x<0Ob|RONuH zB4|TAAW!~IhdZB+*|h7V{j{5#;9)|F2MW`B7RPz>b|dNwOPzLOk)iG(33>-g2<=)a z_!PX0>CNOR(*PE60vy8yq&GmrQ}E&(fC7pXpo3Wfx+FB62=q4GEi#GW0gB_?_{y(X zt6F&4-LL$x*jSNW{Dg+`e>=WPg)WY$7iGu3W%c3~?8-C9oMg_NcVQ$8hi&YKSy#I72f#&JlueDKt%$+G8nb|vPOrdEB^>FKs;b9$2M35GfPQrg6a zbzWDK;=xVY%LglXG5rH~6JCtb{UndWXiVl@)v^IVS{#2G9!Ra=6L1bZk^I8mL0q zB9Y8cWzu-t1845X!`IslnV%rXa+871m3t#+q+~jDEIcFK5*Qc2NQcu()Bh!BgzP`+ zX8l#h;@xxi21Q1cux>Am%7+GrU+%@t${TL!LQW7v~-{y>lyuQ(QbN`;L6j zkWUtSrSZV~k+V-xb|<<>OS*JQ5%lxXc1;wqIjnaTZn`M9TJLLl*Ns=>UiOa zX=P~Kxx=7f@C{)Wp6T1VlA~>Zz6^$8uP|=L-(MkS9$LxlUzS%Q0#m{AVPjzcOE<2> z9?|D=gnf~qK0o;S?6`e6! zdqY={gX&ARUyPnqe&q^J-7+eFO9Fv*tMYz5p2*uf(h!KY`qiJ>90@zTqjH<=dCn6D z`|FjzHZh4hE;)&c2z_pqQahyXe_548B~C{sRXv0hvJmH=f?zn9CvbUyp^8nuFt^>y zomKk|KO2=Tj1=)J(AhQci*sKO9nLPom$ z0Y0}eHio02cDfRX$}K0EB`imQ;H)C6CVokkE0V&oCVAb{M3bnbHqCnNRd~#UI0C1s zT-8+1B$qOhFg!I1n==9iA&yA%H*BMH<-g$58a|E#nrC2U+C2$2C8USFgeRQY**NGk2U4Gd zGvSn|gRiv3L0eFE+QlX}X`jS}Rv}3wDY*v6AXi6)hkF$v%Kb@2xP-!PvQ>G&taR$y znfQ-J^LaE2Gpz$=d-c?&RnbMdo*t(whj`A~fOU|EzB=lC3B0AjZEZ~dLv@pc`6o-U zaTHa2*cgIeEOsgXxGgTGCh zpLxK5ubi4Ge}wVj7{$U@<-zo@q6-NapP!GGzv$1@%_h`o+oI(JuR`N#q=`oYRn^d^ zw8l*P;P3xCJ+hN+#v8qLeQ61L0yYW0{cZ34T_MmNubDFVq~z%a0=nH5KHpAua3M*l z;ZhQCrwieEjU-74XVbOER%@3V9FIwzfwy*%6Cu`om{XB%vE+Uo49F4Pj^?Y`<6sD& zmW?wAL2>A;1ds@>kW#)TnVC5}Kv1O+(9|Jtu_O<2I^A9a+^a@r73OD|qqH*BjHTM( zE>^6qX5PufCo>OcN>HUZMnFGD+nSKUL8=qbvn8?XQs6|3swxZgA*4j7=yk$DLaMSN zBBRuum&+yqtz<0L0Mk??EvnREXSiE453m>3 z;i={ueqXJlBvLIlV%Mj!vhND=9G8-htQcY;3ZlS%t0JHSknijz?0J znQ+7Ph_8+shAYi1dR+Ull+V6h`3AjYlakojm3MTUP6NR+!)MgI@rGR|n9Gx5B{kFB zMo$C{@(D9<~ufnAE{KgyHQ)yAcCnm>70dN4;`W7?!4#QZMN1J-_C6 zxb|=AOY3b^Zg~qf@K1*PqDr?XTu;RROH32#b4EsuH_7Dt8>?kva=#^Y>5H~Ct(6!* z(~uDVtE@1{jwQ!RoG(m!NjqD}X%WSO*DiOX0S_7rG=ptPWvEg4D=CwnUPLx$g@1r| z@BqZ9R-+krS{B9pq^Tiajk=K2FJ0jEHtRaZ%RvJvg%lsX*`a`iS4XW$m~OLDwVVAv z2fdsH7pEjL7gCh}I6Qg&@B2si4Wb4Z!vzq?s)I*vH7Ku)kwk8?Za55aAh(?8sg$&3 zS+v-^)n@Wmj8WHUO;6uw!x`!e9{wY|99@Q^IX|PBepLR5$h5DWOSQOhNkuduLDwwe zqEFlmb4PsI>iu6HNS;bHtDg1&Oitovji zswO$_R4ftm3mVmN8Suaq*V`raM>5mF?hZYUr)zE-n zd4A8wEYSShx2hr(kU?xeXB|8G@JmBujhvXTO7m9@(+H?-zylZe2|*Z z^vAV5TGnUlx5{KjwS#u}+X&arpBMk6oZl|F;c(<@DqW6N^wgSr9YDg;-syEO?M8r7_Wi#T{!ZS(3?d1jx|U3r_Ym}f(4F} z)|j%3PeP8;$M4)tPeZWzqz*)WT!q$Q6%#57n&3vVbCcml|1}bhgq2l9NiFNoaxW;E zcVL%^QTv$VMz9e?+@{KBqwp3QEUtYHV)`~5cVcs?He1l2iM4J3r{Gdi`OgSH*fc2Z z7e#y>gdXfsHzO_OV-<#h9$KkAvj>t(h3ZCi8y-^z}Iu zuoe}3nV^{qG&+4zb_i1u066@}rg!_=PM&3}WrOZt&HE+Kc6Gf4iy z|4|4W>_Ay-JSWHw@=Bd?74@S8ALC1C(t4;IMh{^q?*QAJm4Q~Tl>mGp4;~p5-+gP` zM4WH_TUEXk9)(ZfC{J^DSO3i`JHr1#V?tCVORVDma98_pWV=*#Ygtf?N}&km;Nw-R zd&@-3=9h<0QU{|b@s|_})8?Rt2-EM}a8L0!=HU$Gg?>~dnfwPj2luK&Lkw{NcE^ReYfVqp-YoMeQkcaaCSn$#7uFU3 z%`j`Qs>;7vG0{YVNU>fn-k?{jGo>2p%H4KhZitHiTHKc00@dl+vi@t$$}Xe=u~MQ42q0sX7Zo!g^J!BAm9{ z!z~rj3H3ZiIC_JF~m0Pij)-Qa5*|nAT?jf53=hvyP6wh7vk@@y0`u*l z94`GQQE{7lZ-FI0qO!yt;ockIn&ps(w3A)P6w?&;tx3f>>Y#X?dFnbBb&Fold0ag_B(|>CaS4hfB{V{$lqda&_FvB3|-VF5CQSFP9+zVjLTl3(a z?vLla&U+ckBw7OirutjVb;2jPdz^*CZ&7u7+Kx>hcTVz8K9U;4HqjNmn7aG;*p4>l z{=$GE_)-u@`JbnW@EUgq8sCFefGJKwQ;W<6tdJtzA|*K@>V_!ieyxty;%2`sK_&Se__=|IGpTZ4c(d|^t7Qf;wO&9bif?b)=7uD+ zhKk}}hKh#F>&MD9rA(sgf`Bj#yFm*%R=9LP|K7{<=5un&L36VGpsli9IPk|CgY7j| zGCX71ixZ)7V{U>mKe(7MCa^kxf)#RE%GEf~cW%siQOqBhP`j~=8I$WNgU;KTlgAO{ zla&-~{;lsbx}zwb#6HA*T*chlnnWn zgIM8JD1ED$399jt-T%%J8TT$BJvr{JD)|#No8mX~=578MbbC!WDN(QUp6Rx&jG z23PL7HCdz6N~qC!Ur!IrrQH5I!(1;?$Y#z!lbg&<3t93Hl$!HEo2sBWqT(-b1S+-E z4`SE^Y%5208QGg&OmP=BDk{J-j-~vUXzU!53QsEELMc3C5!qaF=1-8w57?2Nw2cq! z`y4F5&?#^j5bBeE|ACg6+W=r?blaq?)i}UBmU#ElN-S|_1)c{n0XQnQ34?&CL%oh+ z5u*g!8(rV7X^v)|HWCV7C9GaPiu^$Mg9mSxJwrwgUt+G>T4s-^Q~SxE%M9X9HZcso zKO8IE=%}CW32zU!)GhAuJXb8g9NQ4y(QxJhmrs@td>ia#AVK_PCrHYF6Ot%5Rd)Rs zaW?nhTiAjWtM+R>*k1H&ew~io$oPZ8H5XqZU6Q-usaS_QY(fBuhp`mD^NwKVcp|Gr z;ch?|`-nS{5@2%S5m;drc*=l0!Q5mDzru<{#&DEmXXt&NGW(e}OfQjLC-aCI`YM=? z`D$creB=*9co$4J^=0V~+;3!~FrrF|xhx~ni8K9*zsuGO=?XVl%FF#``|`xQ0A3Yp zNFyjbXNui5IUcIn6)?ir^Wx7SW(#NFjkCwA0)}`U0sWb31LyPkXS=>z2b(2`M9Y}d zci41=&-XQ<^j(@ceK+Ca-vG?#cB-!~oBVbrKWvKR4pyasGK0&$t5C&9+D(e*3E#_*q5;P<< zBrtLLEe2O_XI4>)N1m5NpdBoJ(|;!zIt{y!Gkwd|$2p6mh{Aro#he!*GyETNbuVCZ zDWB)$RLnEY=cMO+#^V3HF8}CZ(`u{BMH$7H{HPWbU&u8&$N9vuiy(jG*HQU}^>vrF*^a;I1cfA6ELczb#3c{@X?49>P#5wpJ zyfvRocs_R?4e>ji#=&!gQdl}q#j-OGO+M7(=IAC@{9PMz&7(#Z`t<)8d2TVeP@#@r zyZ9aHC%sm}&ASd%zN*d$eqE`&l=yIqT^7D0qrE(4_`+?sG{WX7`P#RJ#d`Ey!F3sl zN988VL!<|}XeB9Y2IXU+PE`gk_~l%Nk|GoNNaV#vUK-0{ddSSI{Jo@62=T7&pC~0^E)ECD%xiZC+rWpllLM)>Qc5$_ z*!SzBLU*!5Gz=5$?yFzYzuRkg)xPT)evd(_@_UxyL4n_}l=dB3@>bQZ*Xq7BU%Txt z68}^4s>R21NeJEypE&<9de^I3cF!Gu533U%mktQiz9QIRsV*dno|nH?lB)Q3tux#YIYUFA zXs-}EdXg8Nly6^5tQRIlM2U-$ZxWTzN_Y04n7T_^H`8o$HiU{IIzEygJAYDf+hL0Z zLT8PN!M@F&wHVR36!2o@>c1oG${BPb6MYHczhz-mKiU3`+7e*n))n~Ga@B9MwmB!^ZL zNv_zF=KR$U87#YfyknMnAnPnv{-4$a0&%4kT6Wrq9z%k=!78KLK@f7D*c14hC#fII z(hN}&P-fsF0p(~6!6_0~lz};Gk9lkn3JHL!03J-H|IdJC;N2$Fa4q-986#1D{NDIH zD0|rHLI~WuocTqsB69Y7kM?kmdx<*!yL!~Ur3XJjUS$tAJF^oEg=Wi)`D>E(p-=dO z_7NXwVysl-DCZBP^;BchI`9F`{!Yo-aY^Bn;!lOYJM>Z&ugcBn3gk_S!iDdLi-rQ# z6f^SYo9pB6XWk}EDo)BjT!$|ahLNe7e+^w@PptHIFBW}nUfsM9Px7H27kXRheasc% z*lY9oYQ|Eo^TMs!Q@2ME!!z%QPyOtqw9i(l42Y-OY_e6W&dHg5e1%UuRM()RX-E6H zvi5~i^I1Q$;=5Y!BonP&BDlMj5ef^8EPMfaN4Riuefve{l=fDo@6q2D8s3#B-=6Kl z-yeHf=&QVOE7WDxPNnNxfqXw{6FJ}itE0dnar9O-lX3Y}GNM$sq`GI3;cQn&9MLIk zQw+UTur4Y1{B$WL?9$z$r(s1cjK?8|1OY#7FewA4%!jkqDWvx8(#!cCBTA5N(?Mh@-O8-f$&GUIk~5slx9WK@P4Uml3Pw&$@or2AE{3Ekge7V9Yu~5U6(uLLgQ|B(RMtd_emo~l>51@-X zonp=U%a#I5-Q>1M4TSCp^m%(nY;A@n|Jg}>sLg+jYn(1UEc#YgWnV(1*h*}LhkbtZ zo0pz#9$UFbYob!C9m0_BGIj-HFMpCZHM|xb5`KP>So$O*wi@)CzWjxrF=8g9V2Ft( zthG(a-WYF>cyDyMnTnt! z8b1e9v*b`WdtM5%8NnTz=Vz~;$MJ@hp9Qz$pXq+5a_szvN~*^L)~!SDr-cbYNec7_Jpteb}kC)NWCB z!R4G+a)9*7?&BS`hXlh8Lv3pz8f!z(lbqM=tW}5q13!U7L_qMhEn%npa%yooN&Ug6 zM~sj9Qcmo;qe*@_O;o3)YZ79`^VGvf|L@N@X>t1uzT#7b#wqWM8yRO%SMHb7qxu0G(2Tbaqjt+l*y}E`CWu>@Ks1??^Dt@gSIM4y50Qu zu>5*^(mW|N_PBf4PmpuH8fijqx-7)B-kzLa9~x8FdbaLek|V}#p*_UCCP=|8tYjEj z6PUZV$}#(e*!Yu(L&af#uhm%Z9MSP@=nfp`TtZYHfmxL_UulgLBo0QTM7*))e2VSP zMQh{Qv#H18o8C5NGz>NdmIjU@1MRG>bq{G3a_qlc^JThCBo1G0#4IG5 zl?y644{ng+v}aF_4lmSI`b-J!Jah-=HxVWcl$)r6LjJZOX5%210wMz0v7>77I09&r!) z%RWR~SSQzj%%+VbW9zUytaZ^I7$WxA=UV*vAM=m9FOG3Ps~EZl+1Y3o_#f?l=hCX0 zws@)=WAtg4vP$UjYB;v$XO?N7-&ZJ@NXoPq^jF6oA7li#ae{vc2Jq3rDKbQX*#0}a z!g2qGTW__0@XDpk_cop**=*Q%2$TNFALMTgnY_0>sOWOK7j%f}WldkmkmB*3dAL5rp(&@=W1Sft4!2JLIA1+B*HGH2MIZuRRx=jG1?$>jgISAmD zAV91G3j`cQ(x(qdxN@#GbLZ)Cr*TGSaRf;yEpPem0B`kB>y^PqXC$W4*~u|tNTMvO zelzvZR$G_r1*hWPjOtaXFPhVmb&1JR<{H+`1*BXJuV^FD zTACM*x~tgeaU8Px$$mKW*Y{e>?n-q13(?mB)Xb!hRqthj2-7)jw|}=sTumx&C0(rw z>jdYwMqsKO30ITK1FK0V^Q(i1i`L#H)?d@{)in14*eE}fr+=jyIQ4ibUV2Z@eoJ;2 zG!&ObW=`wTpZ@pj!f)(knfgdLR`c-A-=Sdt_oSNS99(Iwc&%HvWg3~8|@J_Qk~?+#_|=W z_fH38R6S{B2)b%x&?4{ot?^NIHK3>P%FeFZInBy$WR-t<7S7lC7O>Km-B>+dBY=eZ zE;i#*&XfUfy%-tSudo}d)2{}vtAxhzT>FfOqd@kuu$dPP4O2|0>UEXvA~1%}lP zEAXt+G5N8^je!oXEvaN@>Tq8+a?W?h(7px{NuL$+($pE_OjWpR8@D+v2#Oj)KbTfRW}ZSxgr;;- z^2bW3Q|Uni&b|^sqw}x4UioGh1@-)|vm~21`^ru0z_SNHp^CL{u6Ms~*?$XJUoy&@ zqo=;Eg2XDVcvedG$_0Z)Im1zld{rn`5evqI%qm+iZm(k1_|<)7G|T@+PRi`N@-iB5 zSNtsuzu)~u)^?mWn%LTyd0o0~`O@?JH`J?HIe(1aFl-o}Ne%8gT#I@&Z6JuN9X@^~ zAqAoPv>{R$n=BYjEesRaUgI1Knd^@hH?H>2*$-2jH)5&2a~ccfw$!(UROSInmA8WH zuC;|lR0QqwdA#n=PCo^`&1AQcSm%u|s(R#s}7HaJ!KR`>LM|K8`j&vWsB=j?ONUcdEQYwvUR-fLMCgUt=B+%%ls z;!+B3{p%3%=l}EbHs$_5|1Y25f1Jerw^&l_pB+|`cuAR4|NO)KpCa&X*e<#l9@yW( z)?_-Na88_m{yl|Ymj3Sz|Nnk7r4cBI-9M85z(ELDru0qqOmzebj6>br;q#Bi|5fE~ zJmT5yP^TI`AqJyk=Bl)A(Bfr7%m6i3mVggXxKIn&%I zDMy7o2P=9Jq(Gm!71q~xL-B3P!_trwt@PGwdy_V+f_Z$WK6F*v4R^3`j<*;!Pg*e1 z>fdC1qn_nv(NA`&;)>uOOXH3k=YV4&-xgpgA8vQrZGXtC%R^*NLglfEDaARAeQ=SJ zo>1&z5^5JGRlWG6j#{Z*q^N9$wLKQugTadx{}QI|cHxMk$XE=?3i#qca6~yQ3ABrl z1s=8wk?{~9;~qdzgu>murCw+{3dcc;+NrZUowD%w!*OgeRXBpd;IcK79&fTs3@irPCyHZotjjPl zBu+Au7^j|yMi(=%Oei>B8ph0|qS-na+Z8=MD4dWigHDo5W>w&=>2`clVtg5xS&k>0 zi+X$fCp?XeCgCILnopp}q*7cN5`=5ogxVA2V1P?#V}P=(+pI3c^&NBIMylfwXUB@V-hDtvlEqrF0Lb6CaO9C*rAC!B&!Hnr3|B z6F#&T#%^CtWHQ3sD=p$aaUd*MMx2MfpDRqgR%u=yhoA&l)5?wLgF!S z6g7T-Cbh%5%nJ{K?<+1x0ky-jq`mO6Re+Pb=wQfrJW^a1iiHUf*#+OIvQR1=Bc{%F ziZ8>fc(h@}7*He#PRw$tD&>99oHNk%@L(1*Ia*jeiH_#N?AXb&+OkkwIbJiVlnEG_ zd0lfWt&tR4CJ)F4GgdLtOt$ug(fi*B1O^wz#x`nZ>FJeXa+1+(I)QK9kt3Id^FpVP zwb4#n$xB)GJE?e5Cpc5aE=viyvF54O$I_EQx38k?AEIGs%`EI~398Vg6a(|%Ci?(H2==^tQnH<{6Ipl>8i=m^woX=Djph$(D0vZ0{A-cWq=(E2>GUiN^bh&F_)vPtK9k^&_yH;qjR%e9|^3HImuU z=wf4vCAc#E{0LoLT|6~-d?Y=i<4Jc**OPYksLk#^X^%ERCJo$WBWwLr+aV=oUud$M zd=3szQ-dQ{zr|A^ANYpm{iE&%v>jVpC42m&bW(FV`u)7U1NC72Xl#31M_aE)OJnL4 zHUP^_R$s5Ws#=@m@9*WW4eM^}>h;IOM@L3S#xs2kp)|-mvag@BiQwxS=vUgqZIEk_ zt?iv~-_)2)@~THXVSM78YEG}epRa*FG)EARce6zgc=Tf^QsO4`RP8Q5g}FH%xDS&W zpBYVeA+C05VAs?RM#e;D#>B+?sY&fr+o>ilE)E&Bku`NzQxjKH=rXl&R)j~#_`#X& z-RuOJ)z$ga@XSbZ3BuoAHKUscrzt)PFq}}@)UY{XklV7=S#E!Dx!x$JSI;bE^cFI@ zi$oZpdS9E6S1j2Rwe6Tn7hK{RW(BwLkR|O)EF;B0$`%1fx>&B2cey>NCNEym@BCn} zVgj9ZzHLI!FfU?OySU!*nPa3M>XA|QfW~Zb$6(l$3wAraZK5XEu^E*~!58<3m@bQ% zhBc$|yNkkm!xSUfqecExx?ye#5=HpZmp-8*A*K~~cIwoU1#T(zJT()RpFJdg6T6yP zYaU-quA2(X!%$I);Rj@S0Kv7s!2(5%JpZ=I8q<`~#fG4! z`H7rKqbbwe)b^B4#fVwD%>Yeoo@ZX~6E5%W;vC%7J!Q5p?3rR{*HrD)tDeWHd$1^7cM8_gb zZIH=91qTQ1lKv@Ux3E1E+~FW-Zg1G#R35%Qu*L)*I#HqF2bp(P#522x$y@=mYlK`n z;_XiyA=l-(Ot2@aYo|0!X&@H4R-wF$+hrqU3cbmEqp{qIB1eIrZ*b`_xk^tbq(+B0 ztS4|c6|%dgOo(Bo$%>*ZHZjBxQ;+9%3;q2G!-Vd>u)Cr2LhfW8QO9%-&%coyqHj>| z5XLUw z^Tc;-iwCpUFp(wvlC~-BFqVmn5um6azCKLW#6_FH!{g_97%n$SK;#YJ2k;d<413t0 zY#LV6*_At*hv)JmNo)9$#u^V{m|K_&dmB%~lU6%7Qt>Gp3KGttVUDNqx$IOvmz-Nt zR%fVCQiSIgg-&P+$tEr^flXAZGmg)FBE)cu@HRXh7w$ZFl-o(<)=VY|@Z@sfuYSM< zI?r9hM@HrfeaLkzZE`uj*2K{n*Cx~#xV!MIY-sp8a(#shU`t=0zdJ=3=5T;%KWCs8 zHcYPXud6YwC<<<$G$QJ_&~}u*kyl8qVuDRx@U&5G`8>PJ-72$z zFcJc72%S*aC!XaMz-?=6@Ms(H_aO+$1{?~WIqVOdD(_;94zMSTi4OA$#t&B5%y1`brc7sQf~#v;#`2^3wJCzG30nB0$ww>gUBY3_aQFI| zQeNn0O6pL>6yA5NcB*S;s7`Torkz|KbQBHa-JEJzdG=t~=ulzfIZw5b(pmkc2D|A7 z-PvIa+k-m8)X#1!!zGl4Glj!PQ`@^JBehB{?jy$17kMvB+xrN;*PF@jOBeLbUSu+t zr-Gm(%{I4U-?!SAoga4#cg~qMMIZfKW^*U;hP%I@Wc?E{O>^OgVHQPz738o-n-)Oa zk&7$Of`5d3;yMa4F7O>k4E^qafHY@r*r8K$j9^yAt48vO} zW;*Y4#%PKX5q|z1Jl_{&Q$>mG>!Xy07c}2?{1dW4K5}49Dtx9XVt_kzj+mNIbZxR; z!6CGZvF#f@;7ut%?=BIa2JW0ot?6?1ed6)9UI89zDv=7JJVM7>h|vy5@7`@RN)oR`$Ree)FG_P~aLgj-GFd0ow@u1gMSrvsPr1HKc6L-HS5 zJ|&x1Dq=DPRHaS#g6YO(t^Dr>FTH3}W0oV2!fB|-ghh9aqT)+uuA21W`)DKQ73U|O z(54b5+CplT%+M2@m*GmA-(HeAFKJT}a+br!lc73hmP6Hs`I{@YxTU;^pD!$2&D*GX z-_0LWj_Q^@xpZcyTic~G;yFC#?9tCsWfl*2KM$LhlNf1=cn0tGrYhO0c2dD(=ag)r zLnhO<@J`C3O6$?4ZMRNJiBwtqj`5BXY0>z_I7i?p?j7=I+Yx!(y9)^zG@jwK zQu8@Vx8cvI`K(_()V%mgXC~EyHB0ZI5??qoSu2(71)UO}r>kh?yxB^`Oyx`3L8S(r zKEj+*_Q7uN3)8z(KLx)w$~yvjFk7D{o-o(l;jL* z{Lkm=?1dJ*5lol#^ioX2E*sR*6h8Z<4f@cOp39WX)lFaeE?W7)D0k(u$)j`b5xm?% za$~B?KxH4fu>vPh8s+88ndItR9Ks)eaJCc#_aBNhEjXjIw5#cYZoa+#nP}ZC55n~*Sw3Amxg?mzAFPP58K|@O?DQtQqTJIeqg0tpCp3FL z089~fS*LOlfD`-fk*CUr=rK-EgU-w~us= z{L5N8=oqxSeC*ux5zts;*lnlAZc>PW&czioiTLVSfqU1SAk=7R;;7@T^W)(P`@*Gm zS~dBE7`yM4c^Y_B;d`!)ucOej7S<@R65e+X~4NA?99 zPL(_}7;Z?rH*_>yyev?0-VfPF9Vag9MC6qP2&Ugt#{!Jh ztGMTm)lIB znYWSaT z)tIA!CfTN+X~+@5*r>@6OAf>xq0eSZP4*mW3|Dh27^CGbcTLpL z^0Uo*f=tUtkLbLh=4UB`xpcj?!{hY$m94X?Oj^-E{45elVz0jA$G&nK@=PY?!~h;_~jwqKx--?w6sB zl=?uR$xb#UG+!SkQL2rJpS7esU7CXjLz)7wp6{MEbBws<`67I36yV*YHAFBe5%IHf zkr=jieasNRxKTc1o*H5LbAURo6w#d4Fq=H8Q2cWpc4^T|YXz~-`$K0L+ zMp|&_UNHMii+gK6aM`?_=CNl_ZN%KsN*c9$+Hf+BQbtuM52xW@kgEi;<^yJ;rK3PU ze_7`dU~R`h9(9h>AMVjL(I4hf#|Tq@rY0>~7aMXa&j%c`2c(&IMq@aXW}4+#hw9R~ z=?K5$up^-O44*!KAqJy9Aq&^&^_dCndW}hqs1Vc}v@68EF*xkH;Gsc5x3P)nLDuOi>s8?v9-r zXbwc}?{3ntI`1Q0eoHg)GCH91#>LQ2`gtLD&Moye@M!z{j74u-59n7Y7W&||`m9T? z0FxFDBSyjBPn8AU6LaGP=5$cHZ7Y-C+sqo2uSU=7^lvRK%9TW4h+Gyw^Quw;>&ZqN zrH>OgmcJ)$9Pk>vdRr!PS>;h%R@sfRWp}7D70cf(nHD0s6hB+x1zP;x!NI{j0*vk> z1tA++MWo45{tbH`&Jm2K&wkz$<=^SpyX zzaK%l#qWM3S)VyvB21-%!wMY*nj6+(WnNBqbg%(MzLPBe%Lir7F!US=%popHb))8- zkJkDqLzG3WTY16?3qh)*e#rW;d~DnceQk0p#~M=z%g!1IY^bA;^3yN|1EkdZc1P~6 z_Z8p@CwEH`ZU^NEETm3#^3l>;F~&=UAvQE@Vi7 z<4Ta_lNYxg^8IpmJMZJWRe#jHG;?N1S4;gWuJ(*^Pn|{#A>YM--0|r|t=mI&>6 z=M{Ij)Nu#@_(KQPBN>@mySBj5oU?@*8M(-$OQLMIQaNPqfKx5bAPLEcbm`K*v)n7z z<51%>nLzg!Xy;0+FAYQjor0`6&4XYPlj4mNTc-*lT`NL!ThBSOWtMMwBvI}u=oI)% z66Y-|eRDo&=y-JUvJ+xX-47Ces%bPHC}Z?12=2S`Sj8HV#;t7EIQ5!V-70(DRQ9#c zuiY=M25}b)GNAS6wq<3+TYBlT7&`1wqDoq|3hsmV8 zo3AgHiC1)LDOSS$uTq}Yy zj>lLyzSz7d&NeduC9n|=VF<~J%#9uK$c)25Mcq``2cOOZ`ioZ-zC_-9SkoeOXjl71N+%B z73hJ(wKr@k2tt7!x}KltW*`prBgR#^luPk;Ej1n*5hR=JczP=C0ijM@CRUK`R@HaQ zBW_e{B<0p&KG{V}CW#-|YA?|BZ->GC~X?!;4ecB}rTC?S`>mS@Ab-4L$){t->-yTNdy zSVMZml;xCeZNR8~(>TXd*n zm(Ly5Gp9Yf;5m$3I396CaO`;+v)W?rK%g-4g_PBTMFTQ=koH>YTGy%V~w$Px9n;?IUxYo_@LbBlQiibiRU${ZQZHP zx7Fuzl8{c-gyZu)kq_3=%~{>w(RE(pqv{hJjh?1UL1uy5V!cRH>zVC*zsQRw?j(5~VvvEx)Q6Mfc4 z6`Tay_E*oY7lr;tViVb}&<}^9xLSO4lPdfq{zKabI-wW;jtzn0YOOX1@z`UAVshZg6aSfD(<9=p|EqV@ z9{xk^wc3o)$Tg->UR~?b;owJl?V2Oe!AQ`NSq|LZ09Nz{T+OUc>2!UQm{zDBHt{hI z#yp+WkI|;6ztLN7h;|!IJ6Qgu`7e2wZ-%qt5Us`ou5!p3E&5$6Uzm1vNw1Z` zZF8Ti>c!rF34y3eo-}|-inuWn6A_5hi%(8v(GjW$)g8WTzvKqX#3g0CF_LkXSD#$z zx%B1ChOYZg(OAj2fSPe7LEM4gZ?g-H=P;8kuI1OMTR^djqJ4J1qb#Ua5reQ2NwmU} zBt=r|pxIHGl9Dxkhn;&klom&jh*u5*0r`=Wn3x!cux7l9dR}vFx8fywq$X6)WX`gE zwcXombE%p`Ca8FHbhPu4=^)xpc+60ZfhJ1>yQPi6m3s)jxMtym*8=yDJ6*?J3r@e+ zCGruqIdT*B^70eq4jj208xpoVK(ik`Pm?mC>5%40_&^vvKv&z(*Vo=8K0d}z+?=wb z!^hVPgWtrKCsp&SvW4Eh7<#ndj?Q@v+Kyy-X_yyUH;~KThHJ@>jdCI*HPlOcy*++t zbNX1^w%%ZUeLY}yZ&@G*B|EjTVw3dCfn^9{7k68%%eZEIpO^u^g5=TOk|k@%DB96c zrO9|EwwC#Vka#$-;oSMDrog#p=UDn*?rc-!Sof^IAqj40LwRxBEj*DY)+xu*2Nx;c z%u3S6M3Y^(fD6cYdnECXyVm=7d9c?Mx-^_cDdy2g*TGM&!hi~W4mSA{8XE5#8Dq0s zPG}F0rq0cigsXEKa(bD9V9aLOID=LpD+6alu+USV$2BsqIPMqYj# zXP{prYa}GL@KdQ>NM4qn0gTUVVY_M$5V&q~*`s{DN+0wnG^tu2O5)vhDA91=sUV!& zaN}m*m1AP)_sDc&V@#y8!YHw8Z(^SZPi*@j(;qzw5n{9Ez~1%><>EQnS-$wx_tlz> zbmZn93V4+#KZ4hG5L1&9!urr1IgGBX_hRv20uxwjE6lWsDs5p&v*lT($#;{WFfOpv zC#n-zD_h}gSm;vhE$W3mk@Fj!2Crg|GY zb_)Z?AVDfbZX(yM$7WXzt2C(^#%FVjIby{+7+{@I4@%b_$zilG2y#OIfL8%^B*{c- z=B~pIlGZ0X7utS0>-=$TkF1*YzJ~Yy=zaLto@c9ggqoOBoVDtnBgk#NM%tl?CC_J# z(CMsW%+T(x^!3MO#QqM+7!@lR!usU&)||02sjBHb>si%*G9rng#`lJ*Uf9JQU~g5} z5`TT4ue80n3Cy9_N9sYrn_sBeDSX!yUKb{=sOxG_67>C{5B$$wa|gR9v3y<82d|8M zCF3OeOEeDuHmh~|{-uZHau{fqV~;&MPg}C$>)q`n4_WLJU;AFh-o0~ zhoPQ~G2K)&VF>z(uMVFwh3U!k{9s$L__9E8nU@k9oCc`w-vf z)6y=yus8uezVPA1t?hz}_ESGZMK%|S*^9^I79Jmb-!?}*9JgO>V_fAk0ijTKb?>xg z$YAu5uP1%4m{#Eat_S#hEs$6kf27;FX@B4cNAl7lLgnAOaVz(^ zeCeC+PIiaLqC56R%;@X=cxy0OoJ7oehxn9+NC<+D$HRNi~4AJcI( zfo2lfaWww>cjG}7f}N^yxdOlW`oO#RBAeQQ;h9j@*DGH#UmGXtt*?MnPdgPkDg>%7 z-UjLL8gDF8*1@Oozb~&o9+@)lIsUwoazgy~7oGSIH5n%_s0X&?J&rmn6WMxW?!U9w zXCJ?P_vkYOqLPMEN%&s5^L6Pj_3K8g|73ULyn|1osLzKqJAOf54qm_0kA7}WdfpJZ zkh}Ie%KM+SmzI%D9+TRU%Li^Yx)!XB7epmB&n-j=T(9wJFzhc`m!Gv>?z#M8L-*$! z=nMkJ)N-=h^L(mr|_9)*CBiGmA0BA!;Kh!yattuyy{8v9j<7(|1LM zFZXKa*}iyW^b%fp>C5`FMsyu%XL#%>+Ce$d;tTg{|7*QG^q=fz1;@%Gx1O{6Q?KUc z7wPY!=xu|~vZZhIkK?(mjh~f_n~^>EN!Sa_oH!u_97W8Y5DdLx!aYz+pl%H^&r)%p)YYfNo=RIF0}WmuiM zUh&LlV@zi8Yb&%UEm`A2_=|ZYk=YIztygI!MzWJ^VZ!p)u{Vvhq?qp6;H~0=E zL|E3HpLV%0Z?U&0=;ZjFuS_Rq5wmDW@_N}W`eDANdRCI=H710P_Oxd+k};A*c4Fem zco{obdb73Xs>{jv-SS!1s@8VMvfUVvDg=Uppir7H3=9Ro&fb42OZueFj&>}@bvAQe zn_&OLi2#}cok6F&F71mRQ5a81Ke*@;(fM+xm)j|~Icgz8s?mD|UQ z@~0m|m#1$lM@AwSh;<9i8=h|GTShH}CA+zr>wRD@x8^y4?_z4Q&yH{7?AwKB3L$koVI z52_vLKBGfO&1C+ETwt;CzURv*Hnd#+lxOr`j|$^w-9PK9@U#ExY5zrC zSQc_lwt3+iv(@eOy|p$&X2{+B>x%;u9Y^Dy>Yut<-?ZCet9MXaueTS&m7RO8?bl*GI|5$*{$Qa&qPglK2B|D(H3%KN=&;uc{~|e@ zc~1@r{&-JLLRU@9^K>Gy@WHuPiKmzE4S7wDe`JIbBjqOy;nL6EVjTX0Cr=|`w!=l&Yq6Z2!`C3&S`Mec3O3(d)edqw3QUNqA$-Yc8O_6APt zISkJ^oAFk~rYB42$cGXz)-1A8YFS13OV)KM%XP&d_{@OTH1R(ZjsFu}s@Qb%7W%j9l{G^&C=Z}u%^F8(mC zG9NViQ(!)SGO+Molt!dWn^@lg!@B5Z>x#R1k%!1Csa@nXr zC45sa<%8ujsg5fT>+tOaUgwq4ZxS`SQ5@w@9eCT_xS9f0?znveY-K|Ef^A?Curt%- z+4w|r1Ge&F`mJ)~lTc`NsA@G~{Ah<&6|Q;o9Q^jDU&X{-_Bpef`FlL;&RBHTnFQ;k zi)wc~miQxlDf1)A(#Vl3wCSG8B0Ng#alQ*PCB{!8R6I1cuTfV+nS21ZN}G` zJ$FhkDk4A7hJG=dvZrC8=DZhIdD@tlmr%43dvaU$KMQG5=AW+?upNZhK_k$UqdRY& zmjE%>q$x*N)^6qy^EsZSRTQSw1SS+O|)Bgl?-#GBcMm=mEwe%7|{|4Oc!+^ z;^%HqB<$5Gl&X56{1DrfkpzX~@otM(&URv`Y#f4-{Z_5eZly3WiGGR#1XK&?l zM%h@8R~WsgswD4xOSDrw7LBR1md2bKvBSd>e+|OolXKN^X!Q)b`VaL&1P88WplSRj z9|L3NfazBan2nkckn{QJHfmXTcGgZy)d%ZtA=jN9Yp>hyG+A@d>qb`P`q(KKu5B^+5l%Om!r2#sxxp2JG%XeQj=@oJ4Nwv0^NsJFc20DPn+cC8MCx*=eZ;SLA*e0HqHs7j&_QB>}Z|Dw6f))6A|18u?kB}J2p4AEUNB7BHJ2| zk!)dFNuxoySSbA321pwJ4Tl%K(K`EB@CjQdm_9KeeB3XPZ8)&0naEr5H2Me~ZsY$a z-T{hKc|FzYH?7~SQS|a^cG4gBcP%B>s+oFcw)eigv9UIxY}9peS!53wqbsHT+Blp~ z`U2X2>_1bj=JQ`X*UGPbi9gc7y5X$ni=lk3Ct`1v_mCUNS?vkR; zYTkFIzOt12Jk+YOJ3BN#K?0KMAxX$af69A!^o_;nN@Jwh!`zK?_5K0H-@LajtsLn* zry^SX^9zWc2nBEHr7XU&th-UMs`BE`>2C%6C*khTCu=6+py2H15qj{L5lx0}&-|nP zjXfb?yAGe5u`#SM}a;wvk)yxSS`MQ*#CA>X~F+CnFe z{62a2SJ>Y1$AVvBvPlpSTpf+FbWiL0DfZ}GmQo4MKvYAyD$sRPFxMvNq>89Z%i`w@ zDZ@OPry)J#`k9}U=>2y$b7Um?|6!M-)TBKw$U-1_}j2WFpX5C>)IvVPFB|#dP(!I66=QW+pOGa3}lTz3Y^5o6$2O`5Uiw~T`Z6n@;A$X5eI04MgP?j0$|1h-NzGtda2IF2KVhJZ5=jGIgb zLz3;}^p|k7Q(_!YZ{q*pW-?;2K*&Er8o^OlXCjK3Xbc1jL)+bCW&#Q`5gb4om?)4qfPi`cI(8KUSq@AQz`_581AYJ))h-c$2+(Yejk5;uNn&7$7>FneFsU{66f+S5 zoQQ${m1qrS0P6j%^?yxv7>(G)`G*O>5J9P1ff;}s11dw+0WpVTv8tkWacFCRXi=cW zhybiZ0mn(a$;^PN!XZ!y1h|eP3Bev_#)_g<|Aj;TiUff(0Sf=;*ct*;`=CCAR?d6# z_{Dq8pSsv2Xd{MHcAedc&B6-EiQM*NIR;w?E#&H==j1dCwK+*JE}baA3VXHD!gZ^i zJ7nxPhJK!13V;2y4Mnk(n2g`EE9qqFZ=_0pOCmdc>T;m-w?lWf-c0lV)XR+(EOtis zb|AqWt0Tm${+g=awFgt?DLWtH>D%e8RKoLj97Wt?jOq_K;ZE>>&%VAc`++Fk+Wjnw zX-dEHd=d76(-)qbzFv2q&>q5Qti`zUeQ~Pn@# z2%qYmQJPlx$0T87AoPC0k2L%9lV-`Tc|pH(^1Eh2u66T%@6+oJ&yR4w79`^n?mGT5 zcn}uTj|KwT=P)haLEP>)#q6|ubU~nz97y~i+Ax1kvOCds;siVKhWJ;e<&a(bo%(A> zdc30zU5=0a9Qx%qRR!I8K61iftGnyT!b!>X184B#mPgl@)Ym@_B7os<1UpHZhQn)e z{K$85R-_RW`+AD`em~Zyd%G*PZ?5Fn%PTD)%jYY7pmJ#(WO=en4eOZfE#a~H9gqXp znWsMw#f8-OpcB(2)FT$75ly4d*h7SuktyHTPlzL9UXiQn!*bJZDJ-}Y(i~ZzW#xvf zzB0bk5z(Jw3#T4EOT2*TJX-2;fOsdZ7j{FS@$<;;3pUTv4s4zvJ=>v%RJ-!}n4JG3 zLGN{5_rV@SMG{VX$EwOhvz)f~|;?#Y8+N%LD(;g)XF)p{v4S+*S!m5Mm?Pu}If z<5g#`LnCsA2ZdcH+}i#o9nrryk!~82A=BO6zx2FJU;29rR&lOB6h31p$ztwrFV4B zrM5uIU>SGCHtypgsdkK>$X6Dp^hwKRy2wiA^t~RxMq6t5z&n-Vui%q+oy+W3?(xPg zf^gY?fPj(=0eRKyE;2>8jP@tFb3f%^T+y3jCA(Ayuk$)Dzj(Xyx?=ofOw0k$0m9ZC z80fZI{W7MS>rPuiOjFH1*FStnK6Zr3m%Lnd*fI)%?BEx#T3!`T z6Gl%h9j&zLK04FR|^jjbP<=5ksdhz=9 zw0@gnoXP7Q9%eNGP`yXqmn!a{@}KluOa8g_ecJWW<^A?TaM7a($JL}`x9Ts8MsUEN zmD{nOLDR66UUS>?)csc7MoMvEAdv{Yq*VqZ3HJB7#|BcZz(o;X({AMg^O3|Vx3 zD%p5hAn2j(Cz&U4`VBLb8mes_(D5LovgRJ+2ayBAj!4D8JLe}G&R8VO;uX$GTzLb1 zbg^Q%nDC*=W&rd%DEw#IuZ_y9Q+j(7Oc#M%!B6fi+m(|#EAI1)AO-u(uB+dHiF)Fc zIT_Y+F)iaI!qV?gI%758cD2yCy%p=&(haW4h3DzHH{$3sXrpa5!Y>Hl=9Q#zd zE=^pQi8^4VJA+KG!hf{f;dN__Rr-V5y0^<7k7(WdU{%`dawjoemG%_2DsA6={C37F zI6LFYm&4y)pYDd*W-Nin^^VP(mBTh}Lsk0!Jaa?ZK~8eG>m4^V7X+8n{l8#89$WaT zCletCSx77zw@j|Wb!=8S@j8KkopU||9wZWGvb6h3I#w`KUE0nvHBrPVugyJfd`2+w zn(4+%e!buQ+qR7mXI}NIZRh7V`0p=^XX?Maym1VCxc6jc9bP8jq}x~FNz3goB{%ro zQ3`{$KV-lV0#Z>H?}5q-K_JBg4cBk=3trgIw23q4EN-`465;owZ>D$C+z(La`G(n} z7D&Zm@lV)sS`I!GDs0`e_jS2prjEhCr-3(p)GyH+otlMTKtJJ{w1S;)x!0-gqvV7g zmZ~BulV!~Z-aE~hNywF69sNepX}Xx&=QyDQ#GS4TRq39S_srcc4wLqE`hCYTchJi_ zNE>e>YtjcIjV3QojRkkSsWG^_j0ma^-Taq?eZnH*3E%q^FG9?#R=M$ZafnBB+iPDu zKV4s6xc+5k3wL>mLm6$wdT!A#FI3xL4QkTgisOIw?a?^r**WU9pv=Z1ntYaR!bocDPT^C04pBrS3@Kti2ytt&j zF_+#(@+Ee4^XqJYMd*bB_R|A|?JUgH zZjn#A!T2g4owmQvz>Y2`y)R0;(6s|}pz*{{a9f&!9e>}~Lb!$Q&nKs^zb^fg^ugo5 z-SWFo|5XGWXhU}()VIghnK#5-M8mMR3pVxsupazo-Ao2jKj~MG=sZ8CJXIEA>KHG^ zfgcsK$Wpq^CHYd6r>IXhQpEGJc|HL~Z#JBqh7hg2P7f04w*Gf&qt5zg7e}2%X&t&X zXbDNTJ9jpP9_4TV^$ojTkR(k1*RNd%zj6rd1x^?Fnl9?seu>+1B#wT??u$Z2d-4~N z^)F|lsL466?@27FP3)pgTmIbmQsK0L^@1FuO>q-(iJ7Fh2wfqmqZL|gX;+v1_w5@| z%B<40&GoQdKYC`CQgvm2qisD!*MkjRoS}62S*v^|vJd_G)9^|A#@3)$RK#5CIgrj` zi?(6L+U*8Em7(vw$AsO7T4#@4s)E0%L}J##Nr!(eENlK&2iSjegX27MVb(Ie7DG0nDVV*j&A38(|WTUDMDM^GbVDp!k&L@=ilu-BU_|E8 z=Z!TV_;P*egGhSppW!!O&ClJx)fLdx4@YVDc`L-NpFI$)tS$zZb%G?ccY%FzfIGcH(n;Q;d(rAGMd6xTBwrXKh~GQb?T=a$iwCY-}*(W*A zovzif4<*_CnomlMZr+wgoQ#GgE-2puB4qbkja3A5>)vX6yyI94#XsFUmc(5^QWw4% zLL|HSvS9X6uI7&bak+_^5qAAqfGpIrSIjgagwPGBYq$=hlkqw(xEvJcwqGUg`ZuR ze5O9|nlG)M=){0v|7N?ZS>i>wNg_#{9l&@L!7jrf<%Sq8P!+2Er+qaRG$e%_q;fb& zFY_}Z%xFM({PvN-dWDn%HQk#J-_q1Ifs8Z2RQe+ulssFpV#BIJDA}zMz_pG)fqGX9WB4~ za)8rfYyCmPaOAa{go;^4Au$nl31)kz>Le&4)i6yp?<-B2=7Mz2y2Z6adYxi&opaTG zr$yGtAHh7SJd?PA>i0TxMhKFTDoD&C(pn>c$oJiB&#m5K-^5@}oUD-N_=acVKM6s6 zhbI*qob({wx(fy5*gXlA2J=|Y8I3xt<)r#T%V&zvsO=@e@Z)@C=IAAD{aNIBcSd;J z6_0rt1rfIa9y}!n^M9+bq-Q}y)yJb+QpOaq1D)beRGehME0;6hV z5|-E#SfReTpy2Nl8sj#fva@N$T%K1uUUBf?Vi81SE?oWC`zvx}bjWE*A2azfg}~ ze8)N(FIFbUE!rv*HR^S#>KXFeQ>8G0;dA0c!y_-?Us>WH6*{!XsJRRr`dvlXs^hL2YStQ}n$$3+*!@H~>GsXc?CX=3XSNe7-QudfGgyP(J+;C(>qGml z?_;1^C9V!Q-i>R1@4n;}aAeQ2k<@t0(F7b3#`Vlf%xQf$QuClPg9_p6Oem4B!YPGZ zYPHAsA9g?&;4>SZ~Qlv;->G&g)0B2uXYZW&c#aVEtKDZ+QE!E}6tJ}8<- zV`lwG6@oXQ;uB%S20JSk<}+%X9y<%Sx+LY4I$wlVdhe<`a)9JsRqTkI-p=8F^HCF= z^m^qKBH3pfdt6xf{&&`{_Z0plf4w(wm zrOn8WS9!36LLTKo1lD@kd~p~ed1GA)vMR!7*bjS(JN-_rC)ynE$RfH=u_`@WXc&S5 zUZ&XBReF$3yVsp5CueQOfDmG*6(~oUT_;|@W`GX!K&P34V#$>4! zjp;J{&%PwOu|TCYL(KA3D3Y~*l&~*4(o`3TwY_m|Mj`lQpS(dY-1g>d6V>+SYJ#cU zda@jlw}8t4cco@&0o!c_r5^l6RiX*9L;JU4;OlEgyW8xru28{r5`DqiVxCB<6~G-(|p#A?wMsZH-3mEQ6= zAbQT<4>dyhPsonyjQX{a0A*bmYf?C@a-0Oi7+{L6q#Z|*qiNm zK+ejMlV6_}bEl}Dg-2yyKP@iraki)us&M^~86egDr-K2~8V1k@A%K1dXoXA}Kybsr zsHuR?CIiTXwEv49QfYt&NJr&fwr+#=i1D~kgFM^Ym?%~ujcql4`Ndh;Qzv>|E)!ml zpFJkzsFY<~*kG`!s3HOz%Z(Yl5bsoZmPhY06>HqTUA8Xp!P{J_m4=ccZ+nmWhzK|h zdr*f}d;KOGu#hIKnAa4Ue7Dn}^Tapp1~U2m_S!dM*0MQYtI?+`??&av=t(E{GX-nh zxlHM~Kd#V>4M<%@4TX>0a&mrWkAHs6y}i<2A;M~=5&FZTDkq{8 z70a+1RKUQZNXO*X5bc>p<*6)HrgjCty~oO{79; zyuX}TG*R!7kmE7oS>xlyNLx63mIof3afm2~dNR1qrlv=OkSt=gS^6vv?*Wlb6?O(! zI>7Y6R&ff%BLGM!+W_R|E{CK`OI_D%!?szwo<~J$ZqXoo(^V0l_U!_H&%O_%MR}@S zf9pP9$ypuX`sETyWpH6zye!QwfH{+yQK_~r5Pr<^K$uE7LP=Yhf=EAJ<)m5yWlXa zGV!FX<#EoPGF!ys)7lpORhQ3X>tBa4{jji4SOwf6wg?hAqFNGRnH63husbry-xxdl&rTCts^Hy~5a5 zg(O#$>!dcz2I7E2V1kvJ%j6TE9_yKC^Au!Tk)=$y(N zGP9uDYX6GguNb(UuwzHNSB8#s(&SJ4PYWGnGNs0{J>0yuk3z3-`Bg!Ay)L3C?=_^A z*6C$^R1zB(Sn3V0B2FXMg+dfznJ0qRATmCE@CZe|Nv7f{pf#9I@&D2gC4M>ea6R^X zAiZDWM(}30c&r@FZ0o159y!j@R$fsgUT{nPj;f3wZ==67Vr^pdgCoS^!s(9tf9+8M zNmL1u&VM5(q4my)dQD?#df`mlhlZorf13m=05o^@<)i5U2!7Q8D60#Kw-@i>Xr4 zOln*Z`hTrfQ%dk})Uo>GwUW5Ytxb2QdHfH*y0_H)_8&iuH2I8Tq!4h_N70;_=m+d# z#X<)}eD^82w;2;nH$e0zI}{eiHowF!kr%^Ru+oJ|1wDA8rLb>;H&1vM&(dDkq>2?^K6Sw%<9$nyJ>GH9YTSvH9&b6EnZc%(hzgTvEOx?R^hr?#eAdA& zM406k!Zr<*(3<+-1MX@U-reLZ>6F}7u~jeL9RT7Gw{t$BTH+U#DwvU?V$pHdw+&{& z85F{VS{IvbMM%tWx_R>RRkbX3@& z8_Zy1426#hUFo8?mbyl!E!*5Kvobu|Ls6;{k!}}8##ymVq}v}RrgG%~R8@^Zd?!H6 z!lTTA93Bl1#_K_Nh$k&{at7A-_9$>!fYXrGN_6DXW56ZtQ#WVh!9*dy&J!_yvt;RC zGYJ&1ofWO9-~Y+uoOE(}FScT|zdI1tmZxW0g`uVW zpprM&g7?_TtjbaC#fabm!ygL>tFS8Kg~bPSO?~50MT9rl8GPFcUyqaPcCOnh= zBVr3i5iyg$DFT5s-?3iDThmrU3i?b;8%FZaTf8`bs($Fbn9g=M@npN`ujn$1pKok2 zibo%sW7{d|7#aLHi@hNx)(IKY^*snGz8}gp97TL^FvEqWuBGcw6gkVGAuYD%<6_Sq z!E3DI7p6Ep2INW^Lp=tKiwCs@x=U>Yu_0R(p{rj?1{)l09g`y&o{X=Of6m~}*pYS? z!**J5Uz_kAUz;f%kJz@n`)?3GZow*=ge(3xC+x0u`VE%kY|HblNvNSUsd^b-8o-*6 zxb@2}xOj(0+k;G>;w>DU@Ck_^QbW3L&IgT@e8(Am3}2%n$FB(EFr#&g`~+L z%>K^2f+Jpgfzj4&O`)@b+bvIeRV+@uJ*26;8{*1wjU10L`x8YU-1$Tl_>R@;Fq^*G z=}UoIiitu;pX%@AV;uhUxfj+LHDM&5EPi$UD0l7mqEyZpOFt^j{cpP>a4b&4K-efv~}Dmz}qCyXeAw2X1$X1ek{iXJS*gru>D|^ zKPC0I_cycmzuQ5Q1d~}rEG=zEME9XN3jYbN|H2cYGJtL~8paIQuZ^42-S%sH$x(&Bi=la`iNh@N>RHo~yLpGW%6*T> zK~dn;29zF*=AX(F?hyCmaGA6MX5l<_E{FBY6@|)Y<--(&`Wz&gX7;Lu3uA)}h;tbZxs*gzgVz+VQebP62!svjtIz|oc_ea!a5C)DSp6HQHphvJZ9 zS}>H8qxT@Nlq}BpjD#VX4 ziO`D!Q}M`BGXwCwB?OF06rTkNK?t|F`S#wVWYCYNIQ6MZO=pqXy6I$pmWeK zhzs=<36%?oP=$WkJQXr$R`Ny#X+60V=nFgGvswVMZ#{KZQfZ3OrN!vdJWdJSp#;ZL zLLH@^tzaPn{aDBZ^F*RhK8q@~e9v6!I0@DLQz6dzJ1za zQPSk5L@3Lf8Zo_NjkbC*s;xl|7J#Fq3W~^$`DIDif{qdH)cCB@J%k{d%k;9 zF@G^64tjDgKdbnS(I~7X?W}1X=9yPwWvyboQACcOq`T2UD=O>=9rCABpwWo@^Q!?0 z6)3+99p%89v)a6dwBtv$gqXJnwS?JP$8NmFe`UNzE{lF4B(auqN+by7#b|Q|8!=nU z$}|r!DOM;R6b?X}S1N!8piy_O8vYbqA)wcb_Tv8ik$g@6xzQfC@aKS((jRxh2tV{o zqpdjoZ&xYOIvS>bUzMRe;ZiLrNfd(;A;@1D?IXD^+43sJpG1;ZFuXM1QvEM z(W0_6$*srZY!|%LdLnal^{^^ke;<6`0d}8<4Swenz5q2lT0_`!WKpM)#GC#$mM7 zeK$Q|{49$2qz|5VCwsPoKvf~FN2e}U-3CYDHs7l%66wV+i25{{X6~oo@qwffKjk<)pIP=e$!Gas`8oF*Qi}y z#8>35tQ36i)sg)zyStT%zhpBfuzMSwj;@l`KZk%hQtCP(gE~(g(|!{%qm`LNZ@0y^ zHBt9FgO1cG%ziEpLb8E@yQgpQXlu4bW7>nDi+-LGoq$KY_(A+r-c!k}=|a$JjKq>z z_|B68^nsYiLEi(ug)zZuYn%!2rfhcdPZXzP;-lCS% z#8Js?1GYwxb*43=9I>3l8Pu3A$1BuVwqJcNH1*^L3OWTBOgdx+{kvAFdmTUBc$B~L zLIWOjfGRw^3y!4f&%=^m8wJZizrh0BG%r6vrShi(hBlor26gp1X+&NH>&_goA4*Ol zq~xLRFp|3nU$DN?J{xpOaZgv?JhIDd|46Dxk?u7_Y_CzWMcS_c@;)sksTjH^6?*WJ zt{V;6>O_U!>2)v{Fy8j&tH;kUwmo^%_w=IvUnO`7}|Yn1}LRbQBYvT<)dcL3PUe_802qH*})Ytovq`~sTb zmfMMMFn}hTvIaRV5cC7~UJ~Agtl8VntJK%u~GeL9A6OqzJk5!?SM&$2yGS?(;{q3!+ zRR)e6gBkstbsDfql3yAfAkj(GD>4uPb(khB9vp@KXq0T0DiD=64vR<&E}Ub6P=x}! zVij%XrGOD+n_*ha*mg_O6c9@KzswQ<>kMnd?Zm+^PO8|&d~mLk z$o4<{b2wi~>e8h0=UJ(gO2w>w<>8)tVp2-9-9+kaFFYwVWV%vK2s7Wb{RRDsIG-2R zhhv7!+J<48t#(S~m-~pL&~L)H!+pX0*+9Q!d$`HhJm`D-?MC8??ICN%>)<2Jlyyt3 z;}3AJwT=WCzHQKCs7hufr6L}>n+8+Gs|%)?Ai_yml+y9hmG7CR3YGyqG!K_rhF#k_&= zt7Ob_j0ev#iN1q=G+(Zrp9w)7N0{MYa5Q%bnSmS^f2PI_gCja5Bl4djvef4kB^ded zy;RTv?BVX>KbzpFt9=r#EUzy&z6JAM)afqog9k+{>Y34Cohs!9&+>EN6IvERslKE@ z{C16eSebr($$s#N_R+1(Ak%i#PWf$@KgZjMFXmMzCa0Q+;Gf{1Lcv%q{kEXRFlTl% zgSVrWxN18pIf^bLEq8-gE#PC-HYEncvxZeJt2V_*PKWy{CJ^oJ{~;_7m*?M2>;M!oDwM8A5 z`wA6yck3zGZKAF_xBt$KPV=Sfk3-5;urU@HMi~&jOl`Ebq3XHy#8~Y<1S8o03tl0l zPkmu53;yQlqf3S-W?2FOsO8VGBc)+|Vl0z}=hc7vevOn0FUP;~wEfg*f*h&A5{C=C!CsbCRty&UTMJIc z&qkwNC_K7Uu|)@~gB||HIknzyiI4#j94KSWVV&{#{z~Uu0@$4bz8F)Xx>_`9HSeFTHbnP!fPukSoygltn zoztJx`0l79RmT>+E!d1F9-Q8QvvUn4LLRw~t%remGhj(5{SAuh4HIn+-*KISKlBja zPM=rdJAJ4gN}jqJ3r$g1FdOj zZeS!d{uazB8w%?@bJ`$o=^N9s(X@+RqjNm9GGzQrmpmi zw*{}!h#Qd>U&UPWfZLO6j-yer^uF zc}4a{WHM-BQ%ER8(r4n<-Bg=A@yPP9ZV%h{DxUtptRl38(0nK93uCzbwCj?c{+B=# zPMy~{dQb(l>-wkK-VUWFF$Iy)8y-)WE(TNOaq}Rc$}OZGVHXO#7tJ4F<0RsQ74X-3 zQ6loXZlD5)-<+yJr+%*@D{_7I48fX~Z!g{+XJEl367TPt2wS>hvUhHF^}cuJc1?xTU|c3Q&L{6=Ww2?@0Wopo zB}oaB?$->^|L(2dOd(wDMt40kgIslgR@rSYrY1b6q7bssY?(DfzZGP^^*MUJwzgy)9}TmUW{%H~}SU zl@m%OX~uucpgv;?Wz(cuhUDoF?jxcqX@nNgqwA!vqh&`m{RI+;{W>Fh#Tfe7B&uqmj%p_jm zd1RlZGbFL34O`k|1fEF0wO_Na&E#@iw*1h7_zo>JxX<~$^QY^luePc!dwerwjnq*Z zIRb$LPQl-AZS$9H1GltoVXMux4Jpr%wMUhP7ne&Tj<3MwOj*{uA4-?T;sV^2G-@o> zmzaZ$k^nbc$<+=Ey3NhOrzGE`o+2EQV->{#ZNI`3c)Z>}*WaxE!T@GapXD^ff3oIUsOLuD@7!m>f6Q*FIdt6D z?DHZe5=i7hx$w+`4=1;zmhNkQ&s+3==`+!LyUSKP{=rIq>w<3fPYIQIjf^UTau0R} z$8qrNyn1;?XB@I-&_6MNG}lnm4Y;z(X({{7r&dZ5QvP=C;Z>wWpF07mhs;)1W z&GuJt1s%~3@8qy!;GJHn-xFuFg-_2!E|ZjUjny%AE;TBq|NJ|m%V`f5II2(#$ZM=0 z=D~_bepn)ZSjq*<4|F{GLHLNTRtQ>B)PUQFS)v&91?oZrADg0q_@u!fUl51y z*SILxIUbC0g>jRh%<>W)Q}2g^Z?e|?^1d7czHm+pan1Ip$5+Hom-hN_!WVUbjAzZDyl&a*w&k?>ifGOz zleYS-E1zw#nag<`**T^`zTb`1iyfb=Q)tId}`=tN3)dcgjKhm}Un63k4KRgFgY@Umn**%FKI?t0sMs zt0F+2YWgj4k!{66|2Ag&eSheFD}y)lf%V4^f1Z1Dy|2eP*mW<3_VrhhCQd%sO&*La znQ6Mi!M!t(RB8%Z^U16WW;YRu0Fa{+2|!rFT!e*#OGGtQ3GSL#%rjjWHjdAaux)Z* z2BI2O30ucq0`CCpQe_YN(%iW>$juz)=_c<_nnz9VqT%L5?x)LJr5?fbYk>kpt-W!zf8O{K`#4SqhvN-oL&5HEE&UmB^tfuh~(iY4!v^)UQqUy1Gon#vEo>~73i?B_qJuzzR1 z(Xu-t*C*!NXQ)qj#ZnY?-*AiLKjfRAGoj