From 463afb07a5b69eb914d77787d7261610cbe0692a Mon Sep 17 00:00:00 2001 From: ctbus Date: Thu, 11 Aug 2022 10:20:52 -0400 Subject: [PATCH 1/4] Remove non-CI changes --- R/dist.R | 11 ----------- tests/testthat/test_distance_matrix.R | 6 ------ 2 files changed, 17 deletions(-) diff --git a/R/dist.R b/R/dist.R index 75d5a00..0515af7 100755 --- a/R/dist.R +++ b/R/dist.R @@ -80,17 +80,6 @@ dist_get <- function (d, idx1, idx2) { #' dist_subset(dm4, c("A", "B", "C")) #' dist_subset(dm4, c("D", "C", "B", "A")) dist_subset <- function (d, idx) { - m <- as.matrix(d) - l <- length(colnames(m)) - - if(is.logical(idx) & length(idx) > l) { - message("Logical vector too long for given distance matrix") - } else if (is.numeric(idx) & (any(idx < 0) | any(idx >= l + 1))) { - message(paste("Numeric vector inputs out of bounds: ", paste(idx[(idx < 0) | (idx >= l)], collapse = ", "))) - } else if (any(!(idx %in% colnames(m)))) { - message(paste("Character vector inputs don't match matrix column names: ", paste(idx[!(idx %in% colnames(m))], collapse = ", "))) - } - stats::as.dist(m[idx, idx]) } diff --git a/tests/testthat/test_distance_matrix.R b/tests/testthat/test_distance_matrix.R index 7186b3b..3acda24 100755 --- a/tests/testthat/test_distance_matrix.R +++ b/tests/testthat/test_distance_matrix.R @@ -66,12 +66,6 @@ test_that('dist_subset works with named vectors', { expect_equal(res, dm_124) }) -test_that('dist_subset provides informative errors', { - expect_error(expect_message(dist_subset(dm, c(T, F, T, F, T)), "Logical vector too long for given distance matrix")) - expect_error(expect_message(dist_subset(dm, c(-1, 1, 2, 3, 5)), "Numeric vector inputs out of bounds: -1, 5")) - expect_error(expect_message(dist_subset(dm, c("a", "b", "d", "DNE")), "Character vector inputs don't match matrix column names: DNE")) -}) - context("dist_groups") test_that("dist_groups labels groups correctly", { From 8d53177175ac172c311a6df7d4459dfeeed36de3 Mon Sep 17 00:00:00 2001 From: ctbus Date: Thu, 11 Aug 2022 10:22:47 -0400 Subject: [PATCH 2/4] Minor change --- R/dist.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/dist.R b/R/dist.R index 0515af7..0396c8b 100755 --- a/R/dist.R +++ b/R/dist.R @@ -80,7 +80,7 @@ dist_get <- function (d, idx1, idx2) { #' dist_subset(dm4, c("A", "B", "C")) #' dist_subset(dm4, c("D", "C", "B", "A")) dist_subset <- function (d, idx) { - stats::as.dist(m[idx, idx]) + stats::as.dist(as.matrix(d)[idx, idx]) } #' Create a data frame of distances between groups of items. From 1596abf40061661c249863491a2704b1fb7e1dda Mon Sep 17 00:00:00 2001 From: ctbus Date: Tue, 14 Nov 2023 14:34:40 -0500 Subject: [PATCH 3/4] Remove pkgdown workflow --- .github/workflows/pkgdown.yaml | 35 ---------------------------------- 1 file changed, 35 deletions(-) delete mode 100755 .github/workflows/pkgdown.yaml diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml deleted file mode 100755 index 63cbb18..0000000 --- a/.github/workflows/pkgdown.yaml +++ /dev/null @@ -1,35 +0,0 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/master/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help -on: - push: - branches: [main, master] - release: - types: [published] - workflow_dispatch: - -name: pkgdown - -jobs: - pkgdown: - runs-on: ubuntu-latest - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - steps: - - uses: actions/checkout@v2 - - - uses: r-lib/actions/setup-pandoc@v1 - - - uses: r-lib/actions/setup-r@v1 - with: - use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v1 - with: - extra-packages: pkgdown - needs: website - - - name: Deploy package - run: | - git config --local user.name "$GITHUB_ACTOR" - git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" - Rscript -e 'pkgdown::deploy_to_branch(new_process = FALSE)' From 106eeab2461a4fb859fd0258c299f875f58a0d76 Mon Sep 17 00:00:00 2001 From: ctbus Date: Tue, 14 Nov 2023 17:04:54 -0500 Subject: [PATCH 4/4] Update some docs and README badges --- .Rbuildignore | 16 +- DESCRIPTION | 5 +- NAMESPACE | 4 + R/long_format.R | 4 + README.Rmd | 9 +- README.md | 674 ++++++++++++++-------------- man/pivot_to_matrix.Rd | 4 +- tools/readme/centroid_example-1.png | Bin 21212 -> 4630 bytes 8 files changed, 363 insertions(+), 353 deletions(-) mode change 100644 => 100755 DESCRIPTION mode change 100644 => 100755 NAMESPACE mode change 100644 => 100755 R/long_format.R mode change 100644 => 100755 README.Rmd mode change 100644 => 100755 man/pivot_to_matrix.Rd mode change 100644 => 100755 tools/readme/centroid_example-1.png diff --git a/.Rbuildignore b/.Rbuildignore index bff01e8..c16beff 100755 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,7 +1,9 @@ -^.*\.Rproj$ -^\.Rproj\.user$ -^README\.Rmd$ -^\.travis\.yml$ -^cran-comments\.md$ -^CRAN-RELEASE$ -^\.github$ +^renv$ +^renv\.lock$ +^.*\.Rproj$ +^\.Rproj\.user$ +^README\.Rmd$ +^\.travis\.yml$ +^cran-comments\.md$ +^CRAN-RELEASE$ +^\.github$ diff --git a/DESCRIPTION b/DESCRIPTION old mode 100644 new mode 100755 index e18b3b6..eccbca8 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,10 +8,11 @@ Description: Functions to re-arrange, extract, and work with distances. License: GPL-3 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.1.1 +RoxygenNote: 7.2.3 Suggests: testthat, tibble, tidyr (>= 1.0.0), rlang, - future.apply + future.apply, + ggplot2 diff --git a/NAMESPACE b/NAMESPACE old mode 100644 new mode 100755 index 9a0bf69..7db2e05 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,3 +10,7 @@ export(dist_subset) export(dist_to_centroids) export(pivot_to_matrix) export(pivot_to_numeric_matrix) +importFrom(rlang,as_name) +importFrom(rlang,ensym) +importFrom(tibble,column_to_rownames) +importFrom(tidyr,pivot_wider) diff --git a/R/long_format.R b/R/long_format.R old mode 100644 new mode 100755 index e94cdec..77d945c --- a/R/long_format.R +++ b/R/long_format.R @@ -9,6 +9,10 @@ #' @param obs_col,feature_col,value_col The same as \code{rows_from}, #' \code{cols_from}, and \code{values_from}, respectively. #' +#' @importFrom tidyr pivot_wider +#' @importFrom tibble column_to_rownames +#' @importFrom rlang as_name ensym +#' #' @details #' The parameters \code{rows_from}, \code{cols_from}, and \code{values_from} #' should be provided as bare column names. diff --git a/README.Rmd b/README.Rmd old mode 100644 new mode 100755 index e0cb0f8..bb5c8f0 --- a/README.Rmd +++ b/README.Rmd @@ -2,7 +2,7 @@ output: github_document --- - + ```{r, echo = FALSE} knitr::opts_chunk$set( @@ -19,7 +19,10 @@ set.seed(0) This package provides useful functions for distance matrix objects in R. -[![Travis-CI Build Status](https://travis-ci.org/kylebittinger/usedist.svg?branch=master)](https://travis-ci.org/kylebittinger/usedist) + +[![R-CMD-check](https://github.com/kylebittinger/usedist/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kylebittinger/usedist/actions/workflows/R-CMD-check.yaml) +[![Codecov test coverage](https://codecov.io/gh/kylebittinger/usedist/branch/master/graph/badge.svg)](https://app.codecov.io/gh/kylebittinger/usedist?branch=master) + ## Installation @@ -235,4 +238,4 @@ parallel to save time. library(future.apply) future::plan(future::multisession) dist_make(data_matrix, rms_distance) -``` \ No newline at end of file +``` diff --git a/README.md b/README.md index bfb27e1..39a380b 100644 --- a/README.md +++ b/README.md @@ -1,339 +1,335 @@ - - - -# usedist - -This package provides useful functions for distance matrix objects in R. - -[![Travis-CI Build -Status](https://travis-ci.org/kylebittinger/usedist.svg?branch=master)](https://travis-ci.org/kylebittinger/usedist) - -## Installation - -You can install usedist from github with: - -``` r -# install.packages("devtools") -devtools::install_github("kylebittinger/usedist") -``` - -## Introduction to distance matrices in R - -In R, the `dist()` function is used to compute a distance matrix. But -the result you get back isn’t really a matrix, it’s a `"dist"` object. -Under the hood, the `"dist"` object is stored as a simple vector. When -it’s printed out, R knows how to make it look like a matrix. Let’s make -a distance object representing the distances between six rows of data. - -Here is our data matrix, `X`: - -``` r -X <- matrix(rnorm(30), nrow=6) -rownames(X) <- c("A", "B", "C", "D", "E", "F") -X -``` - - ## [,1] [,2] [,3] [,4] [,5] - ## A 1.2629543 -0.928567035 -1.1476570 0.4356833 -0.05710677 - ## B -0.3262334 -0.294720447 -0.2894616 -1.2375384 0.50360797 - ## C 1.3297993 -0.005767173 -0.2992151 -0.2242679 1.08576936 - ## D 1.2724293 2.404653389 -0.4115108 0.3773956 -0.69095384 - ## E 0.4146414 0.763593461 0.2522234 0.1333364 -1.28459935 - ## F -1.5399500 -0.799009249 -0.8919211 0.8041895 0.04672617 - -And here is our `"dist"` object, `d`, representing the distance between -rows of `X`: - -``` r -d <- dist(X) -d -``` - - ## A B C D E - ## B 2.603430 - ## C 1.821423 2.047355 - ## D 3.472394 3.727228 3.056922 - ## E 2.672239 2.653173 2.734967 2.069155 - ## F 2.843420 2.543180 3.369470 4.373791 3.129488 - -These `"dist"` objects are great, but R does not provide a set of -functions to work with them conveniently. That’s where the `usedist` -package comes in. - -## Working with “dist” objects - -The `usedist` package provides some basic functions for altering or -selecting distances from a `"dist"` object. - -``` r -library(usedist) -``` - -To start, we can make a new `"dist"` object, containing the distances -between rows B, C, F, and D. Our new object contains the rows *in the -order we specified*: - -``` r -dist_subset(d, c("B", "C", "F", "D")) -``` - - ## B C F - ## C 2.047355 - ## F 2.543180 3.369470 - ## D 3.727228 3.056922 4.373791 - -This is especially helpful when arranging a distance matrix to match a -data frame, for instance with the `adonis()` function in `vegan`. - -We can extract distances between specified pairs of rows. For example, -we’ll pull out the distances for rows A-to-D, B-to-E, and C-to-F. To -extract specific distance values, we use `dist_get()`. This function -takes two vectors of row labels: one vector for the rows of origin, and -another for the rows of destination. - -``` r -origin_row <- c("A", "B", "C") -destination_row <- c("D", "E", "F") -dist_get(d, origin_row, destination_row) -``` - - ## [1] 3.472394 2.653173 3.369470 - -If rows are arranged in groups, we might like to have a data frame -listing the distances alongside the groups for each pair of rows. The -`dist_groups()` function makes a data frame from the groups, and also -adds in a nice label that you might use for plots. - -``` r -item_groups <- rep(c("Control", "Treatment"), each=3) -dist_groups(d, item_groups) -``` - - ## Item1 Item2 Group1 Group2 Label Distance - ## 1 A B Control Control Within Control 2.603430 - ## 2 A C Control Control Within Control 1.821423 - ## 3 A D Control Treatment Between Control and Treatment 3.472394 - ## 4 A E Control Treatment Between Control and Treatment 2.672239 - ## 5 A F Control Treatment Between Control and Treatment 2.843420 - ## 6 B C Control Control Within Control 2.047355 - ## 7 B D Control Treatment Between Control and Treatment 3.727228 - ## 8 B E Control Treatment Between Control and Treatment 2.653173 - ## 9 B F Control Treatment Between Control and Treatment 2.543180 - ## 10 C D Control Treatment Between Control and Treatment 3.056922 - ## 11 C E Control Treatment Between Control and Treatment 2.734967 - ## 12 C F Control Treatment Between Control and Treatment 3.369470 - ## 13 D E Treatment Treatment Within Treatment 2.069155 - ## 14 D F Treatment Treatment Within Treatment 4.373791 - ## 15 E F Treatment Treatment Within Treatment 3.129488 - -You might have your own distance function that you’d like to use, beyond -the options available in `dist()` or `vegan::vegdist()`. For example, -the RMS distance is kind of like the Euclidean distance, but you take -the mean of the squared differences instead of the sum inside the square -root. Let’s define the distance function: - -``` r -rms_distance <- function (r1, r2) { - sqrt(mean((r2- r1) ^ 2)) -} -``` - -Then, we can pass it to `dist_make()` to create a new distance matrix of -RMS distances. - -``` r -dist_make(X, rms_distance) -``` - - ## A B C D E - ## B 1.1642895 - ## C 0.8145653 0.9156050 - ## D 1.5529017 1.6668670 1.3670972 - ## E 1.1950614 1.1865353 1.2231143 0.9253541 - ## F 1.2716160 1.1373449 1.5068729 1.9560190 1.3995495 - -## Centroid functions - -The `usedist` package contains functions for computing the distance to -group centroid positions. This is accomplished without finding the -location of the centroids themselves, though it is assumed that some -high-dimensional Euclidean space exists where the centroids can be -situated. References for the formulas used can be found in the function -documentation. - -To illustrate, let’s create a set of points in 2-dimensional space. Four -points will be centered around the origin, and four around the point (3, -0). - -``` r -pts <- data.frame( - x = c(-1, 0, 0, 1, 2, 3, 3, 4), - y = c(0, 1, -1, 0, 0, 1, -1, 0), - Item = LETTERS[1:8], - Group = rep(c("Control", "Treatment"), each=4)) - -library(ggplot2) -``` - -``` r -ggplot(pts, aes(x=x, y=y)) + - geom_point(aes(color=Group)) + - geom_text(aes(label=Item), hjust=1.5) + - coord_equal() -``` - -![](tools/readme/centroid_example-1.png) - -Our goal is to figure out distances for the group centroids using only -the distances between points. First, we need to put the data in matrix -format. - -``` r -pts_matrix <- as.matrix(pts[,c("x", "y")]) -rownames(pts_matrix) <- pts$Item -``` - -Now, we’ll compute the point-to-point distances with `dist()`. - -``` r -pts_distances <- dist(pts_matrix) -pts_distances -``` - - ## A B C D E F G - ## B 1.414214 - ## C 1.414214 2.000000 - ## D 2.000000 1.414214 1.414214 - ## E 3.000000 2.236068 2.236068 1.000000 - ## F 4.123106 3.000000 3.605551 2.236068 1.414214 - ## G 4.123106 3.605551 3.000000 2.236068 1.414214 2.000000 - ## H 5.000000 4.123106 4.123106 3.000000 2.000000 1.414214 1.414214 - -The function `dist_between_centroids()` will calculate the distance -between the centroids of the two groups. Here, we expect to get a -distance of 3. - -``` r -dist_between_centroids( - pts_distances, c("A", "B", "C", "D"), c("E", "F", "G", "H")) -``` - - ## [1] 3 - -The function is only using the distance matrix; it doesn’t know where -the individual points are in space. - -We can use another function, `dist_to_centroids()`, to calculate the -distance from each individual point to the group centroids. Again, this -works without knowing the point locations, only the distances between -points. In our example, the distances within the Control group and -within the Treatment group should all be equal to 1. - -``` r -dist_to_centroids(pts_distances, pts$Group) -``` - - ## Item CentroidGroup CentroidDistance - ## 1 A Control 1.000000 - ## 2 B Control 1.000000 - ## 3 C Control 1.000000 - ## 4 D Control 1.000000 - ## 5 E Control 2.000000 - ## 6 F Control 3.162278 - ## 7 G Control 3.162278 - ## 8 H Control 4.000000 - ## 9 A Treatment 4.000000 - ## 10 B Treatment 3.162278 - ## 11 C Treatment 3.162278 - ## 12 D Treatment 2.000000 - ## 13 E Treatment 1.000000 - ## 14 F Treatment 1.000000 - ## 15 G Treatment 1.000000 - ## 16 H Treatment 1.000000 - -You can use the Pythagorean theorem to check that the other distances -are correct. The distance between point “G” and the centroid for the -*Control* group should be sqrt(32 + 12) = sqrt(10) -= 3.162278. - -## Long format data - -Many times, the data is not stored as a matrix, but is represented in -“long” format as a data frame. In this case, one column of the data -frame gives the row label for the matrix, another indicates the column -label, and a third provides the value. To get a real data matrix, we -have to “pivot” the data frame and convert to matrix form. Because this -is such a common operation, `usedist` includes a convenience function, -`pivot_to_matrix()`. - -Here is an example of data in long format: - -``` r -data_long <- data.frame( - row_id = c("A", "A", "A", "B", "B", "C", "C"), - column_id = c("x", "y", "z", "x", "y", "y", "z"), - value = rpois(7, 12)) -data_long -``` - - ## row_id column_id value - ## 1 A x 11 - ## 2 A y 10 - ## 3 A z 10 - ## 4 B x 9 - ## 5 B y 7 - ## 6 C y 15 - ## 7 C z 10 - -The data table has no value for row “B” and column “z”. By default, a -value of 0 is filled in for missing combinations when we convert to -matrix format. Here is how we convert: - -``` r -data_matrix <- pivot_to_matrix(data_long, row_id, column_id, value) -data_matrix -``` - - ## x y z - ## A 11 10 10 - ## B 9 7 0 - ## C 0 15 10 - -Note that we provide bare column names in the call to -`pivot_to_matrix()`. This function requires some extra packages to be -installed. They are listed as suggestions for `usedist`. If the -additional packages are not installed on your system, you’ll get an -error message with the missing packages listed. - -The matrix format is what we need for distance calculations. If you want -to convert from long format and use a custom distance function, you can -combine `pivot_to_matrix()` with `dist_make()`: - -``` r -dist_make(data_matrix, rms_distance) -``` - - ## A B - ## B 6.137318 - ## C 6.976150 9.036961 - -## Parallelization - -Distance calculations can get computationally expensive with large -sample sizes. With the installation of future.apply package, you cna -compute the distances in parallel to save time. - -``` r -library(future.apply) -``` - - ## Loading required package: future - -``` r -future::plan(future::multisession) -dist_make(data_matrix, rms_distance) -``` - - ## A B - ## B 6.137318 - ## C 6.976150 9.036961 + + + +# usedist + +This package provides useful functions for distance matrix objects in R. + + + +[![R-CMD-check](https://github.com/kylebittinger/usedist/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/kylebittinger/usedist/actions/workflows/R-CMD-check.yaml) +[![Codecov test +coverage](https://codecov.io/gh/kylebittinger/usedist/branch/master/graph/badge.svg)](https://app.codecov.io/gh/kylebittinger/usedist?branch=master) + + +## Installation + +You can install usedist from github with: + +``` r +# install.packages("devtools") +devtools::install_github("kylebittinger/usedist") +``` + +## Introduction to distance matrices in R + +In R, the `dist()` function is used to compute a distance matrix. But +the result you get back isn’t really a matrix, it’s a `"dist"` object. +Under the hood, the `"dist"` object is stored as a simple vector. When +it’s printed out, R knows how to make it look like a matrix. Let’s make +a distance object representing the distances between six rows of data. + +Here is our data matrix, `X`: + +``` r +X <- matrix(rnorm(30), nrow=6) +rownames(X) <- c("A", "B", "C", "D", "E", "F") +X +``` + + ## [,1] [,2] [,3] [,4] [,5] + ## A 1.2629543 -0.928567035 -1.1476570 0.4356833 -0.05710677 + ## B -0.3262334 -0.294720447 -0.2894616 -1.2375384 0.50360797 + ## C 1.3297993 -0.005767173 -0.2992151 -0.2242679 1.08576936 + ## D 1.2724293 2.404653389 -0.4115108 0.3773956 -0.69095384 + ## E 0.4146414 0.763593461 0.2522234 0.1333364 -1.28459935 + ## F -1.5399500 -0.799009249 -0.8919211 0.8041895 0.04672617 + +And here is our `"dist"` object, `d`, representing the distance between +rows of `X`: + +``` r +d <- dist(X) +d +``` + + ## A B C D E + ## B 2.603430 + ## C 1.821423 2.047355 + ## D 3.472394 3.727228 3.056922 + ## E 2.672239 2.653173 2.734967 2.069155 + ## F 2.843420 2.543180 3.369470 4.373791 3.129488 + +These `"dist"` objects are great, but R does not provide a set of +functions to work with them conveniently. That’s where the `usedist` +package comes in. + +## Working with “dist” objects + +The `usedist` package provides some basic functions for altering or +selecting distances from a `"dist"` object. + +``` r +library(usedist) +``` + +To start, we can make a new `"dist"` object, containing the distances +between rows B, C, F, and D. Our new object contains the rows *in the +order we specified*: + +``` r +dist_subset(d, c("B", "C", "F", "D")) +``` + + ## B C F + ## C 2.047355 + ## F 2.543180 3.369470 + ## D 3.727228 3.056922 4.373791 + +This is especially helpful when arranging a distance matrix to match a +data frame, for instance with the `adonis()` function in `vegan`. + +We can extract distances between specified pairs of rows. For example, +we’ll pull out the distances for rows A-to-D, B-to-E, and C-to-F. To +extract specific distance values, we use `dist_get()`. This function +takes two vectors of row labels: one vector for the rows of origin, and +another for the rows of destination. + +``` r +origin_row <- c("A", "B", "C") +destination_row <- c("D", "E", "F") +dist_get(d, origin_row, destination_row) +``` + + ## [1] 3.472394 2.653173 3.369470 + +If rows are arranged in groups, we might like to have a data frame +listing the distances alongside the groups for each pair of rows. The +`dist_groups()` function makes a data frame from the groups, and also +adds in a nice label that you might use for plots. + +``` r +item_groups <- rep(c("Control", "Treatment"), each=3) +dist_groups(d, item_groups) +``` + + ## Item1 Item2 Group1 Group2 Label Distance + ## 1 A B Control Control Within Control 2.603430 + ## 2 A C Control Control Within Control 1.821423 + ## 3 A D Control Treatment Between Control and Treatment 3.472394 + ## 4 A E Control Treatment Between Control and Treatment 2.672239 + ## 5 A F Control Treatment Between Control and Treatment 2.843420 + ## 6 B C Control Control Within Control 2.047355 + ## 7 B D Control Treatment Between Control and Treatment 3.727228 + ## 8 B E Control Treatment Between Control and Treatment 2.653173 + ## 9 B F Control Treatment Between Control and Treatment 2.543180 + ## 10 C D Control Treatment Between Control and Treatment 3.056922 + ## 11 C E Control Treatment Between Control and Treatment 2.734967 + ## 12 C F Control Treatment Between Control and Treatment 3.369470 + ## 13 D E Treatment Treatment Within Treatment 2.069155 + ## 14 D F Treatment Treatment Within Treatment 4.373791 + ## 15 E F Treatment Treatment Within Treatment 3.129488 + +You might have your own distance function that you’d like to use, beyond +the options available in `dist()` or `vegan::vegdist()`. For example, +the RMS distance is kind of like the Euclidean distance, but you take +the mean of the squared differences instead of the sum inside the square +root. Let’s define the distance function: + +``` r +rms_distance <- function (r1, r2) { + sqrt(mean((r2- r1) ^ 2)) +} +``` + +Then, we can pass it to `dist_make()` to create a new distance matrix of +RMS distances. + +``` r +dist_make(X, rms_distance) +``` + + ## A B C D E + ## B 1.1642895 + ## C 0.8145653 0.9156050 + ## D 1.5529017 1.6668670 1.3670972 + ## E 1.1950614 1.1865353 1.2231143 0.9253541 + ## F 1.2716160 1.1373449 1.5068729 1.9560190 1.3995495 + +## Centroid functions + +The `usedist` package contains functions for computing the distance to +group centroid positions. This is accomplished without finding the +location of the centroids themselves, though it is assumed that some +high-dimensional Euclidean space exists where the centroids can be +situated. References for the formulas used can be found in the function +documentation. + +To illustrate, let’s create a set of points in 2-dimensional space. Four +points will be centered around the origin, and four around the point (3, +0). + +``` r +pts <- data.frame( + x = c(-1, 0, 0, 1, 2, 3, 3, 4), + y = c(0, 1, -1, 0, 0, 1, -1, 0), + Item = LETTERS[1:8], + Group = rep(c("Control", "Treatment"), each=4)) + +library(ggplot2) +ggplot(pts, aes(x=x, y=y)) + + geom_point(aes(color=Group)) + + geom_text(aes(label=Item), hjust=1.5) + + coord_equal() +``` + +![](tools/readme/centroid_example-1.png) + +Our goal is to figure out distances for the group centroids using only +the distances between points. First, we need to put the data in matrix +format. + +``` r +pts_matrix <- as.matrix(pts[,c("x", "y")]) +rownames(pts_matrix) <- pts$Item +``` + +Now, we’ll compute the point-to-point distances with `dist()`. + +``` r +pts_distances <- dist(pts_matrix) +pts_distances +``` + + ## A B C D E F G + ## B 1.414214 + ## C 1.414214 2.000000 + ## D 2.000000 1.414214 1.414214 + ## E 3.000000 2.236068 2.236068 1.000000 + ## F 4.123106 3.000000 3.605551 2.236068 1.414214 + ## G 4.123106 3.605551 3.000000 2.236068 1.414214 2.000000 + ## H 5.000000 4.123106 4.123106 3.000000 2.000000 1.414214 1.414214 + +The function `dist_between_centroids()` will calculate the distance +between the centroids of the two groups. Here, we expect to get a +distance of 3. + +``` r +dist_between_centroids( + pts_distances, c("A", "B", "C", "D"), c("E", "F", "G", "H")) +``` + + ## [1] 3 + +The function is only using the distance matrix; it doesn’t know where +the individual points are in space. + +We can use another function, `dist_to_centroids()`, to calculate the +distance from each individual point to the group centroids. Again, this +works without knowing the point locations, only the distances between +points. In our example, the distances within the Control group and +within the Treatment group should all be equal to 1. + +``` r +dist_to_centroids(pts_distances, pts$Group) +``` + + ## Item CentroidGroup CentroidDistance + ## 1 A Control 1.000000 + ## 2 B Control 1.000000 + ## 3 C Control 1.000000 + ## 4 D Control 1.000000 + ## 5 E Control 2.000000 + ## 6 F Control 3.162278 + ## 7 G Control 3.162278 + ## 8 H Control 4.000000 + ## 9 A Treatment 4.000000 + ## 10 B Treatment 3.162278 + ## 11 C Treatment 3.162278 + ## 12 D Treatment 2.000000 + ## 13 E Treatment 1.000000 + ## 14 F Treatment 1.000000 + ## 15 G Treatment 1.000000 + ## 16 H Treatment 1.000000 + +You can use the Pythagorean theorem to check that the other distances +are correct. The distance between point “G” and the centroid for the +*Control* group should be sqrt(32 + 12) = sqrt(10) += 3.162278. + +## Long format data + +Many times, the data is not stored as a matrix, but is represented in +“long” format as a data frame. In this case, one column of the data +frame gives the row label for the matrix, another indicates the column +label, and a third provides the value. To get a real data matrix, we +have to “pivot” the data frame and convert to matrix form. Because this +is such a common operation, `usedist` includes a convenience function, +`pivot_to_matrix()`. + +Here is an example of data in long format: + +``` r +data_long <- data.frame( + row_id = c("A", "A", "A", "B", "B", "C", "C"), + column_id = c("x", "y", "z", "x", "y", "y", "z"), + value = rpois(7, 12)) +data_long +``` + + ## row_id column_id value + ## 1 A x 11 + ## 2 A y 10 + ## 3 A z 10 + ## 4 B x 9 + ## 5 B y 7 + ## 6 C y 15 + ## 7 C z 10 + +The data table has no value for row “B” and column “z”. By default, a +value of 0 is filled in for missing combinations when we convert to +matrix format. Here is how we convert: + +``` r +data_matrix <- pivot_to_matrix(data_long, row_id, column_id, value) +data_matrix +``` + + ## x y z + ## A 11 10 10 + ## B 9 7 0 + ## C 0 15 10 + +Note that we provide bare column names in the call to +`pivot_to_matrix()`. This function requires some extra packages to be +installed. They are listed as suggestions for `usedist`. If the +additional packages are not installed on your system, you’ll get an +error message with the missing packages listed. + +The matrix format is what we need for distance calculations. If you want +to convert from long format and use a custom distance function, you can +combine `pivot_to_matrix()` with `dist_make()`: + +``` r +dist_make(data_matrix, rms_distance) +``` + + ## A B + ## B 6.137318 + ## C 6.976150 9.036961 + +## Parallelization + +Distance calculations can get computationally expensive with large +sample sizes. With the installation of future.apply package, you cna +compute the distances in parallel to save time. + +``` r +library(future.apply) +future::plan(future::multisession) +dist_make(data_matrix, rms_distance) +``` + + ## A B + ## B 6.137318 + ## C 6.976150 9.036961 diff --git a/man/pivot_to_matrix.Rd b/man/pivot_to_matrix.Rd old mode 100644 new mode 100755 index ef94d78..7e51f83 --- a/man/pivot_to_matrix.Rd +++ b/man/pivot_to_matrix.Rd @@ -37,10 +37,10 @@ generate an error, with a message to install the appropriate packages. } \section{Functions}{ \itemize{ -\item \code{pivot_to_numeric_matrix}: Specialized version for numeric values. +\item \code{pivot_to_numeric_matrix()}: Specialized version for numeric values. Deprecated; use \code{pivot_to_matrix} instead. -}} +}} \examples{ longdata <- data.frame( sample_id = paste0("Sample", c(1, 1, 1, 2, 2, 3, 3)), diff --git a/tools/readme/centroid_example-1.png b/tools/readme/centroid_example-1.png old mode 100644 new mode 100755 index 0d83f4be9a6f0856d08a3c9625a930e0ce394523..ceb5ae9c8941f60e948d0fa75232e859132597c7 GIT binary patch literal 4630 zcmeGgeLRzE`w96fR0@%gq)?YUE=G^ZH693`Jqt$#JOBO4<@} zN@5%~^|EG_mSMGJzIA-G_u1*3-+Rsh715p-YCj~BZvbQ@*rT~x$mGG2m$~E0H6SX0sssENe7@9 z07?WafkYyu zkVqsaEi8*c4Jk$J_KN}l)fVZqq%FJ(4FD@jJzQN5#+7~`ElYkrycT{$u{`e8{pu?y zn`bd6=#URKg1$$B?C52|8LJhecJCaUm1x}EO5+>Mj|VwZ7i{fE;9K#?HHVxC8|}Sb zcLT61-8X%KqZ)F`&sLc1-bFT0NlCT-8~BTKkhOHK@#=}$4b=K|>qrVRZnU8{0>>35 zwdM-&krY_oDSJDpn%sqxosW9rJyoKU8wyTa5g1{nMJ44g_2m2q1-Z_{`hkR=1j$Gb zNr~>p_%t)bn>Q}s$=*;N9eKyzho-G0TX$O~sVBigZK{_|mru~YIEuJ$zUKlQJ#~L- ze)0B3YK7xUES-BQ9+JPwPGS+ZW^LgykcqNJlZm3+}o6rI6WP*#V*U3wMvGxM*}-U%V)c^<#wJyUqWv3)CUFkSgba{gNjkF z%1!6W_<5@4bjzj~aw8+6^%Pfg3e(CbVdQq&f7P*pCpzGdk}f2963?#uKQHo&4}8`B zv{HP$jYZREhU!SNGCRT|`{hF{9a#OVYt?YQqrRMo_7Cl#Q(V!-LXy~&@KN?JV;<_I zT_PAoX+`Gz9Lhr9m3wuchCv9P+`GCKUiO-reSI_YBp=V`uSJHuUTGbYxb9{xkI2pQ z*1&EPZ^x>pT=TLIPLc2PTfi!FPKHR{2cR&}lF0YUA=}tOGRNQnTlh*#EIeehF493d z`osQw;&0vit3>^3{W;lM^`b?HemUAqpQq?*~2 zm8;psJPRJC`ufxY?%=nRWBtxhod#pN18HxIj<1+6R;3bp-Blhx zE-4px(W&n&CqpU|?$*)h<(sukejvhI{TA$4QL) zaAS|Q5_YCO^K^dt%R)y;dfJ4sN9zsGpDeuzZ>KVIFBDxZ6cxboMu?hWVI1|3>!dQa z_8WTg3;Dmih8zO4z8FJ3Z|F^Omzh_P?B{KRDn0k|@RAh6toE)o{~4(EAgehP`O|RkMS}CmmWuq&a~*9P@g!r}5BP>-;G+BGp{GbSu*`t*Y1ri6_jsh{gTrKw zN&k&S^QuroD#P?)sE{!GIGLWMi)j|0sq@tqOlB(e#T&A4>_b(qWN$V{(R$!QrRLB8|ebTx~lf)y2wp=1xXhazx5a6YrlVLbAi@ zk?!nB{nbnnE9TZ>&GknH$+yc3OLHb_r{6OFVb6}olNAQygSQeefw!{d2~ z?Fn&12XLzcDO>kAO(Uc^rATSYpAd7HaQUS{v%LMXsZIwDgv8P8pdRuyTWL}oU?W3{ zSzfI0z~1pl$EQ&?a>mvLxzs-~8IYiqB#kGel_t z_vSyx@yCEUM!_F3F#Wz{{@1Z*av@d@tWcRJyrw30r$P^dfnJgiCZJGdPWgT@+uiGO zhz*zmCYD!P?O-W9e63L6-R{%&?PeX-dR_aN467@%TGp8`W>sk|EG^!+pbI}kZ5j&x z{;p}|tqA$3{hJ#8o}^J@-bpk1`x|A%@(>X2YL_lMT{q6N{iNiW_-zn*h! z&1ZXc)7(dfJA5{Pbt8sc?uWg3B3Qv(?m=@BuTpfZ-afY!2XW;uewduS=A8VjzVT-7 zZY^@g8bMa%g~9&1+B6NBhb7bYW_2IO8VAI1^LX7MW52TAJ^{YLoazo*Hb`#n&C}Rc zMW}s+@8eT>lcJWA>5|1J3M&TNy+Irvcxym}Ka$bJNy&*F@BQxFsePLg;--;{7(MKJ z@r};$VIS~QaI>!>W>}@J*vDoAuK+UgjwK*m0an}BUwIplHR1I_iLsi?I zbupc5Xq>z|wDdy|t7i!lONm(SKZ3p7MrIP`ov6%4nsnY!gcL_m{hN-#1t8 zQ18ne5xZG{yN|S5fmQb-bHFfl%d;t_ZP3SGho$o>Mt#(;n%b3?HlKq zVGR`2c`JpKZf1FRp}C8h9MsuCcV?_+7da~G_P7{vmPK37if0TbIJ|gcP#3Hijlyqr r(kp$5F&?n|`-b`#-B7=%_yzskC*&8a9>c?>|93pxkgm6O9!>ohw+bN> literal 21212 zcmeIaWmuG394|VQAO;{H8<0>0q?PWFZloKLZlycK1OaJjqQG8sb(iBE9qqKw+e-CH+5*zWt(*~iaSfA(2bTFv*-MiR-l z=9)m-2m1umqLm!6`1DgFv=2lpkAl_KLv!E>LR<+pb!BrF+rEgDom8LgTB{wL_1Xg6 zxJVL=+leJ?WSmPbZzJA`uzY;=Hl6IsXPqZhEJpcCn7pdjJgN-lCEI*VCP!z_(~4_k zs7PsF&^)VH&Aaw40%JE-d`5Y#=p9V`iI$AF0(ET zzr)eeuzeR#7}4|nNai}~QR$dTy0R_2V0|UyS(4jqdfD?G)AKZkR_oUYhp)N)32CCQ zc1#UXEIjmUi%-5O7DT1VoZ>CnIrvrIS(uJXiX;-d@^NM^5t>7^BV|}Lja$%^mQ`F8 za@T+D6)%EyiS1oK-&&VO>qKEe^!Mae_ir^k`g?mq*e&nt3=)5}JEI&lQuo$YHh2=d zeY~VjPa3@B4*ZxfrzY`CY{@iD=q--6!-e3p?oSX(+Tn^7RdBDbYjkq=;x z%Did7@zcm5XT4wXIvL*xjpA+iNsxN_g;YhTb0AdjO5u|aVX9MXwHxr`C*S+XeA9ei z^A-jLG~KAUiu2&p!?%t2%DPxY<`2ZL%-^1{&JVtMg_ES1ggnILvn3hjT8r?WU}ZvY zYg1dyc}Isy@bQc%WN^PRr~EuCULnHD2hn%RiOqLTCKT+YjfMHf&Kf>&splj~&k9u5r$4>mVpcNL7%yVIS^b8t zn+b(u)Oh;Cq|6Y;%J-yti9DrFr|;P07=Frmz;bGvwGgnZ^Q*?WUB9QkTM4cpcB|@so+5xMDc!b{a>j2-Qt9n$=vCVH3MH#(Pv5CQ|uaOG^{> zJcfx%+ew+fNrL#peM$2VF}#|GUvv_&44KXOVd@PFI>K-Cl6oi(ML(l}@-oLIZmhq{`|C}IQnL%5QsQLUP?mK6Ju)@ zH<@JRxaHRxe*5g8kCOt?-8_=Jw_Eb$uu0V3(c{THqbCZ(P~;|hguxSX_t^(*v%61P zj|lQ6;O!51S;eKY@F~2-7B&J>$*d-Y4i*mFE#2oD=eyUw%uY^v%pIbd_)tAUJ3TqE z{PJ>mw1E)Je_pWff21C8pSb?}zvvHS0wFQsf4)FVmCU7`WJFmZ*Hmjq$VIFB$7=X{{cD=X66NN)4?aH?#RW`Ff!{i+UU zYBa7`x!M4}us3A~{q9?<7|+gWKT-M#T1JS0-NyED1mb$QJ-at)mDbtGm}XcqRBq92 z-kZoX)ac_ec+7C$2ivkY(PDMC(qURVBs8?w@GG?GM8O-mem6q&ByV|nc|wC+Aq0Xc zPo!IE9Vv_)kQ{0bxR9)d;#S?qCM?yUwZX_UW<}P|pY<6&S2u=u&-)#AM$s!~Nxr?g z*MddRfrJfbyTQlY6?6vGBkl;O`~M0xk+N;}J$B~8ETE$Z#w7@6y&Q}#Z&gw}OF>tJd<>(ex^xVP_`82>X z%jYnqS@6xK-d$n4QcU4Jm%9G4ujh#K%fZ>^L|XbVig^79MgZl+S%b9tM1P*rcW231 zzHw>Ch@9<<-Q_r2(6vea?wY1+%qmGci7782JmM8pGvgd;=(+fMWy~aa#BsJRVL?M@ zNMjdF30qau+*sCtz>-@4Qh>=2v)+9xYNrXM(sXuUaJ&<6!J?6)MAX>$qa_H7;P%K> zLAyBt`JsbNqRD!lWFMP^ac8S@t>hmfU#KGp#nZd+H|Lv8G6MDx2)f&-6u*w=j7O~Q zahuGnoUOS{N8L~zg-2=Eh#0A92^oU2J_NdNO^yE5aJ-tdqHjHLOGHo~oLzGJ>+wb* zwCl@@pZC0!QUtZ$+@zxL+TITatHCQJd*_(GqQ_gTU+m(zf0eN=LooI+!GjsiCG2FM z7i%5Dd|mx7y1 zN_h4-uC;!dPsE7(_Gm*z4+ebW>{z|#g_)HmHxrFpc9%gUV`p`V9C>$1CRxnt$wMRN z6U?#l>18ilyjKTMPCB};GfEQ{V^8|Z$%NIvVB-3E4$9B zP^TfcN>gU+nS~JJsNdww%;9@;z8tOE!;Mcxl*0dF2}9ZHKAq7}a)3n?mwJ<8usV&S zv!SKpRWapiHw)$kf9v}d5t`LA)v!murn{)#_&1IVu&7s9Djy#M!S@1Yjia&7_jh8I`&e!0L(L%goTMvzPL$=4aL<(` z6O-JOiN6Bhd^UmeZLLEBZ}P_|JI;@GlX5uDs98GA+((5FD1Gm8J~;?0Hjs;Sd9$+6 zyU?{tn7U?K8m=Q5Y)N&;dHU|$A&sNoS<>C?YVkZwHIb01FE27HZH8sM_c!Pd=L4F3 z5xsm2yUAh!P05~1@AGFVw6fd%Y*dM9ORe-E zh27F54tj<2qsp40I61h$s@gGLY0b4OPv*vQWnCGEv$Cg|IVfjWuf1Z1j2f^jCvwZL z_9ih3?yT)~bHtaKwTH~Xm4k;u+3OxG;zvuAVtPYzD8DRV(J1oRGC=SWKR!3n){a|mPYm<>R?S)H z_Tohu$?eE|(3~0&*>-^777#_4u^^42}HBfR}(@INAm;0gIFZfm$ zZ`*^`Psz(2bPn3Qa18%PEqC|oHLobsFWCGzG$jY=+?SIn-UVhi$Cp?`;u6r&8e(Fn zFWM3)1YhEStOh{(y7(TEO8u<}fR{FUaR|l#35PLfD14K|=M)Fpz3x=C0}Gg)5WNpo zcH><&^Zt}x#S+0+hCWNa?Qr=U^!ziZ6O12YAg+TM($3QkHJyPNSmvMra9Z~KMb2wK zao3y6sxJo2jofUVTdy5>p@WD>y~nwnk3b_#`K#i&%1I2-jH+=z!-Z62!tXGdwTH)o zQ^{*B><~THdOek2^m@DR#%P(Lv$C(B^afRKJ7Cj`CquF(X0u&A+_1 zM}uug3*-UD94XN1E@Z{(WCY8<4yOuYPqWAt=7infPGFQgVxsurOcu*go zmhQVeIh^$0lmx&nI3L?UK2dxDWBo@MWwuNhg(dFdSe;6WAj8Cmo8me}x_42BJMjc! zr)tj5&gQ2_zno#&@rjY@hkgLCrSznU&HFJM)w}a#7DVDuiwvYBoLdXnE7ZA^mi5@$ zo7FmiXRP&98a{t-tW~d~xq9AvzR9mu^JVdUW!hbcW`s&;NQi=#uB>!LMTH_8e5SU7 ze6~mjnyNPk;BW$uopwG1IijJ;lbN|qS>rIxK(XMR!Yw{Fn89XH?U_;qJ+pmd`~PRCMYU3aq(;`=Y} zzgbay^y|-rma3rdLiU!=;Wis$wUdlQKM z0gp;Zq>~wekoO^<{}imB zy;fYBj{2P@7Vu<5gQeOT;g6(4xGfqxg`VCdeTt1#BaZv6sEXrX1N2J^jX&sQK75o8 zy)ED9y}!IaZW14@A}H_J@aqRjiI-==4?Jtnsh3f|e-+q|d9ptzcK$0O+oC&euTN-4 z4q3z5m&B*qCvs>=&T0NG*>|r~G(OALsvpmK)G@#|;_tf`kG*v?-h57P;`{5h@X2~+ z_DfZ1{jakDYjTA%`jM!Im$~^nm`2OpTEK zJwidu?7L|fr?TN<=bkn*UtTD^N8f`biKfXZ0PV|7sO@M2$J+1Mg%pO4)78p#9y{zD z4ZAO|DujlHDiaqCiL_D>H@Uv5rSDQXw%)9dN!S&~$Mxv|_{|MbYL$A@`u9x$2}rvOxvUX^SF zb%vNz<9-}+XTJHKQ;k0rrW!oyBmU(KOXIjMyPTcuXKSZ%S%t4hFo-D60p=|S(Qov2 z<9fEyv>2@_8p+ki{}nrJn=^&W`*7AHsSq|@Ww*vIy6Keib|Z&agGGij6sBOxd|K*W=wcC+-RY<0g^Y z=2_|n%#BYvO=a9G%|auI*Wq3$Kp{O$=LOsZ+}SeXT%4Tu`S|#T(v3wrdV081`Xd?| znazrPPxszcu|ebVtIqE}+lsY5KDNFXc;$91W89cAc4GRfN+)Aw?dKJI-FT@ZX}Qjh`O4 z+AvV>;pBqTV2sQR`aI_j`(zd#nXMS#HB;--Ye#16X8QWCK*nn56f z{8&pXuR3lwWFQPXRcVB-gZ^i}x7_)5`EbH4B0=L=@Rm2cxF+OEH4t1x%o{b0KB{Bn zo!`9!>&WK3kq{ki7H)moQO~)OzwUFom$!bvQjv<12y?SQF8M_0!yjOc{i@Znbl`i? zR>$f6t5LStVfJuxOzgC>udgp0y7T&dynR-BC7Q6A6DvErn4heYlN zY#cza>Bo-en|U65^^V)J#<&N^gP}th=`V_fjT_rwAlU*=hSPes5;G7Zl#6z>g=bo_ zDro`dUKNf=!TeMHUYIzOOK(S9j9Zf@#Vp`;Hw+UcNU*drltzeZK zqJ_@Suq&BsMBQZnft^9(IIbo6b2Ci^j56Yd-UR9Amn}-X1wdAe;0l3D44Ym6P_6bp zwfNcZ%YXz0;^2v6U9occ^BXO)g?m{Ab69y||MV6uxt4yPZySw8@}$I1u`8dc6uth_ zbrvEa7jC(g3^P5<(o5KQW; zgn{1J^8W+hA|=insxkFmde3MM>dYS86+L|eDi7c2w)o{cEb4$Eufq&hw;sQ1V`=$iX0n&xCOqJ^UW>*XGBFu zYb^SL(k_Iq-^(*D{l&XJ4aSdbp2W8`BZm~r!CAhJ$(sZmuf!`$Qn`QUHyI(=G;nIv zvg}J{0(Dh8dB@TXyNz!5=x0PXqF>DEJnkGdf6?_4=bF;I=&x_YD=8jxS5+l6vhH6T zwTqdbpC11(I8+@Jvl}atkJMKWE4Y`Q_L|GG=V9<83wVS3mOf$%5V7v2n3H~w?U+g( zlgLUk6T-PiQj(g=^{bE zyEV2{sKFx9K#Gc)ryDjHx_h--)zjuJffzY15#(GyxF1hW&3erHm790of(?8QUVqU4 zTQjuKDI7sqrTJ-xi=J2uEWRRM9+aht9aa%ya~JJ<*jl>hz7Q~L*T z3>r^6KHhZ(;)HBSdI0pqBb$y4!0oybaH)HVUZqwMBpgPpHhp3?(<;fIhlZet%~DJ( z+?z4ePR(bf@Gl&WY!5!_j0LO-vqg#mKwjK){gB?D?rCw$n^Rw}UrBCJ+J^W1Dk|1SV7Gz9h(|>s9G0x+<97_=DHU+pv; zXxNTv6_5#sO&U=RjJP?Y?Rc+L#0hTGcr6pb0L&LY zm@sqRsv4K+`n0-fXg7C;%T$uO#mW+?KLv%R(S8r3lpR|P;=^P)*p5cTvqw|(otJ;> zwQt-6WG=xaf@d#`lh1l`ADor`A7Y%*l8Ww`m^#fe*412Whw`^8w$Q^u>(d#=&s z9xAt(m`A=M+uZB_ZL;>Yj zDn}{do{o`QtUzIv?I<(XgWUP*EDRVL+K?})ik*!l(>nsRy}mFNZ{lw=hZ;^h*Yrnb z(tPjfw`c9V(Bg-r05C$%H1bq=E!O7y5B^r`cy^ z{<0ARJEhuT*ek&UuZJtAemTp^t}=<`n@5oRckQiH^cRtF<7e{zHcKoL<;C5OCLZcj z&2CmCUcIn>3g5)ai-mcWzM{AEJL7yIE^<^R+0m8mS)%e2O^yb#cU@y@?S}+(w}V`p zTowvib8|$Vpdw#%>vT2OxySY3Zx(i~rhQoJtW?Qr;*Ap`w{WjH2AEpWy_!0PT1gOB zGrHtwNiOaa3SZlycE>vR@629%6Dr#p^#Qp;`y#qk!ETe)Jf9)R%m`7=U}Ib^;oXiK zqKEueIqbE|w~m`HP-3w8=H_)0JP{*kuxnyv7;>bx`n5{g5bXV0T0qn3*-k)w+S#Vn zs#3|7bc$;*VYeIVZX}frr#J0Iye@Q(hF1Y;Bt07(Ubr|vY?kI&fDLB8Jzc;Rdl-~= zC5wG&v7q9`PnrBDnKD^i@VXPxlxsv(fg9M9>!+?#vW$wdj49uQ9B00CWaMk)b`5-% zr0ehYdA^q7ktnkTD4d=&2gN5IySN#*aihpZTR&5e4+QUAt?pMwBJ|%7S5AH8{!;yD zQg_=tR!{Zo`fJ&;YN!xX3fV@Jc7w~?W!PJ#k)Wi52v^3x9FmjhVEdH1gGMta_cKd3 ziih2ixnSH^;*H!B^cr|MRBeUe+peP5Glt{Z?sd8+Q)-KO{$M=WFBMVAal3c>m*L0A z4)+FMC(a^|%?ZoC}901ca929DaqTsE-R2*Ut)!Ql@ANW}a;e4tBGk`1_Ah z^7=!eUR@&`y;pB$TNb|x$;wRhK2!9PKOFTYT=>BtxT4lC+6>C_Dw}s}ln-sX2A;>( z6e$-=*r=;WN}92;zv&n8sy2)~1I4qWhYJ4ICjq|!1}%*@9UVkyLL6%OG@Z)Z5YZjU^R+j6g!sBx4E%Bjoii}va_21 zoSAMB9sXu@XO4dQ@J&FUAt%aBuO>uv{M>l+%k-BWMVW-zNM9;?_-wrglVysS|F7?4 zp;z)Z3Ydl`-)?ds$V<8rDXj@|^W!*I8TbYi9wsR?HXh|>(wqH2s_(E-a7)Q{E3kzg zzP~vIb8Rb^-w{vG?dlgsC=t`xT)5>5*rp8}$wsG^ctrcF#PP3|MNMnhi4C1;c^!ZE zcKX=8gfH%Pq$4Zp#-7tBJ@Ec0B%e%t-Ns>-NylTA7lvo%l@Bn6%JrIdrMOWy9@u*B zV`)Ya^&_-t5bBx&>=@}(IF=u^$Wp!z0`_X=EBzVW@Pdr0;zmRg@?O|m8vhByzgF10 zg{?lB_xp8`hlugtZoY_Vih{0)wzTukbCaF7|tTk!R21!C5@2=%H5D zgT&K~c#~;T$ou@GW`L=GesH0m6~S$;REg&BQ8q=k;gs6c$U&@K-Oj49$>%!UcJ~5~ z4+6?BPy``GmE5(TwOr&~R2FIcRY_Ua7UJnEOmwQNDY)TpOQvte_A;&?(sxg#{Kz@0 zYY3DXe<8wyjgrs2vezG$M0Zw~I1J6LwxiN(zb3RVr_Bs={XjCeXCNQopUQPh{{cN? z3M;))z{lq=T=6rFu`&^Ypox2PP3(bAl_Z2Y8*p)!F`VN{tGrYGxvY2TduVUOl@jlV zYR1ga9cx1^V}I`_>g5`*vDCu)_G4;{M|Kfu7lweAYcbPm(l(|9zNBJ@!B|APhTc~L zG2vm&y@3VQfXZMvlXr?-A%eP60%^-}gS$!s>G~`~9?P;=adOy7Rb)SAebk!5yvAvM z&OT@N)ppc)UfGDJSJ2hRh`Pg%ZPbmo*f4iwX$dh>m@{i|hU(>|axtx}^f#4yJN!_} z@-CQ}4aSTLrsG&H<(cC@Uo2soCGy6F3ry}E#W_SR`pQKV>BydpDch6K|5A^{!hC<# zaJsy980uoaUA>iot7VoPSn09#O+}Sa`yS2QmXgqeCfA}|@^K&-Dw&SdM_c=5%5(e} ztAK)9T?fR>p1g`V26fy~{mA4d<1`3J{KjwCa{h@D<92}u%J1A0l;a3 zK(63skwm_+02IUCOP;b@Br&uq-1szWWjtv+)um4$TYJ#daSM;k$aU~;jm}_gUfpcD zi<+VZVH+dG2%t;y7kE{w>+fomnG(kQPLPqQB0h=Nu=e})Q2D8m1IO>Y3%Ry-%A(@* zR9#G|6mw541jQ7c?i5F#o6*m<$B6jlp&mL=0%{dK2{oJ4Hc#+BU&qUf&kbfNpOqSy zkE^;_XAMC>1I`p(O>U*Q4#=yv!mV72X!C$#xLPf?Ej+Y-$Z(=7yFc3`CSOt(alTJ_g zUBj>xxj&5%$!=i}hGo-a2y$-Zv=qr7K7aHTW zzkW~>Bsl52`=!zmGGLtT4$vrj#}o&Ip4>3m9*s6-wquTwCJ5hW)J6**%L+%j&;m-{M%a8&`Od!|O7>>pj^6Z(r60sBYq!lDUwHI~NCDo24|vl%dF_9hlPNy_O?66$ym508&l% z-pGGs66&!y*XTpf$jCTX^S!Z&Z2B%FEGdQWQiNi@X3y~;W{uAw?lkHd`9OM5Wk*S> znj28J&4J6n5P`NM(6F9PHk`X35T^vSx$PvLj;*Px!Xn+WHlyWkCup@KWJ3}0k?Hq6 zO*|+y4QRT4zRbnM$`tu*&1$2dYA=H34$CZjq{P*&yrN>ZcwYr{HgW0f26=9ojJ&E!61-LQ57p zt^&7yUJtSr0%+Azzsj~8xLcFW0nOYFQmajAy9Bzja*T-q3-ujIQ6VQ#?7?vh#Or+#i!O0?`eI=(N@0 z&yu)$)%FQqTVH;LN2i@GGiF#E0gCzHXM2R0xzOvYp9lkG5(9C6FCDc}99Qge`yu-I z&1n<&vfO?$4r6=Fc67c42%}~%y8%D=_TlM9N8XiK{c0U!FKJ_+EuI8Er&Q|&L6=n} zV$V6h!`Z0F$h&B*GxHP3am3cTtR)Vg36}pIo**;vy^E%o!l%1mb2bI_1$_?e(ZXT_ zqEBcqY5r8^d@sIPY9kw=;a1mz*13ViAl|bjr5%ot0NI}LMGG;8dD|u%8yibGA(i>8 zl9H7%Q_rMZW35$j`gY}I5fj>$hcR&uE7dP(fPi@SL)8W9mfHcm;!$K|m6V5YVN!HowtgBh=xny#cL{@~AyPedZWT!|LU$Zti9zI4Mm%F$@bgZ=vW z_}H>>{Q|gblmc#0FT02N@~V?g6QvR(M?jX~$vv8VqAB74mDp-LoUbx+IsxMFBUxfWkG7mFp{LVbhBb~x@~Wye zB{RJGb_x4mzgCMr%X9;>OK%^g-Nj+LAfpP$(L|G9{RJiB`8D1HyOi9#wE9aQN6?B> z4|Yfu5my*T&3|JnZZRaTH(PK28@rSgzc_dNq#@I{Ruja>>} zF{J#({vrthkNIFa&=RCOyYPuw`+>93G#1$cww6tbPd;6wh&S?~i^M+i1fZ@5@NpY58nZ5lR;Rr0oRm3Y)$eS7JgIPgmSw?B1-2Na z()he46QPvEZR7l9wA+uZms~}}wc4S!@GS~Qq_N6rUbO!j42X)EHgsxx<7h+V{uuKNG$uwmZ=PB#{BN%Zyp-zF+8kocHh5L| z0_kX|+Q9$@OlPIG{}tZ%JD=eRpi&S8QfGX_{Buq0=+GRx9&|=CcE4z+Vg&sx1|)2{ z(j|iLv#=D#%ECNsj`uilW#(?o{uT7uVrjk zoK3tXAtS@pyv7^Nigx7D?F}eDN*@=s|EIifgKGOnt@)_aZ^WuGAbwS^DuOBwFsC3n zj;yhI{UI_kZu$qO8CB*ZKi#Y|$>_h6C3a~!Ca8h;X6E^4`oVtFCwM6MP*{H2Mb^me zT@i>;m2Ir>;f#ySlZYtSnDV^5yqJRU5KAkoc>Nr@gnm&}5zr$0S)Jya8ha9X3N<0# zL-LHB-5kw#yymMNrhDGmm0N(Uiy@sLZN3=XVpxT*2NZM>3DV* zKKh`LDyx}K&u506>2qk0_J7)%u0hS0QaiTFz=0X5AR!U7IZ^gi@;#QGjZ%8BR@MmV z?isI=pfD68`5c4{6kXR!0g<gOzXKJ_Sf(Hg~mVTR6jGOprj6L++Vw6A7{EAEN+{ z?Moo;Hn&P2Z}`mr4FQd#%>T1fuA004z_!`{tY`g`TNepPL;XG)?Y#vQ(@{rU2DPEz z{wSwu8@5Y`k!sBbTLL$mUbz$qQkcIS26~b{_tVDHy>4mMK{yNrd~!NwGU|JO%;|{- zAox+KoW#2a1fm21!EF&yQm9pbYGmXSKT&U_H93!M%;sFHg@(n5TWUhtSB$KnPzg1I8WcH+|*=o^Y=51aD_@J!EUmvcA z1&;9$m%C~@e(DdY`Akd;xnUHoM*C-$t3bKesQ1gKy6gyqf*6bRR2;27tfZx#NE*#C6MkuamKz6$PL?qhcIuvrc>Cc(Ys1u59FpxuBxZ` z?#U;J_*7rDeOBCG{`@MMWvm82mVKBJO>I8|Z1jxG`9Y@ML>5P72!XCf5lY2O$ zZ^(+Mqoyn7Hwt)1dW{e-SX>=*prJvyb|ub0D#CLA0iX5|ZIwZdV`S}1Afv{k#ru(h zQ0`p<2# z7~@spZZD)$`QZ-(sR9eCJ^jDP?tc{|!rdB9v}dU)wy3U}kvNP1t+MCND}dekEa!z0 zJzUZO8N+iK&xn0D%KcOpU^#Zp#`OKJTdd^shi)ETx&gQjTd#wQ)545LQiN|2%PBxW z3{OX61F+kQnV|p>70H&7wG+xU<9QmUf)1g4pyXpOw?+ppT-OIWT%a`EHbdXwtbY9* zLPFDr4l%QV#?Li%bsqrC|1soK?B*zR{=WkX59UnSh?0-iGBR&S1nA5h zhXQl#8OMp#Qy4dL&>RQAhnd{eGf)#r+4|7|*KV$NA3)JIW!T0oD_(aPDbgDj9wmqZ zF;$sn5>y$nr=?(S@68YG(#1=LuBDj+$Kt;}D`!?Ob$&J1V^sHKs|YuiN|Tnm-<&>Y zmsxeLygqWP@#Bx2nCo^`Q?)L+MS2xKuroIHVYz9E8apStV%Lx_fAGzDN}>%&qx&{0 z3$+{@h&EQ&^J(~fZCj|v16m}!zejS!;kFsJfI~X zdw<(lY}hTmBf{U|r}guURIRNuaIn__i^8wnc2oHdOsXfq5Tpa4wkV)4nKT`* zQWAUirHPrS2><$qzw%X~lWc(CT{}&HAv=f!xm;tYkpez){QT)E-%y2>sz>Cpk*BZ6 znMu5H znJ(oPebrP^i$%|$1bYuvj0IIZDjTyR$Hcp@({c#m25~61c!13AqrKH!V|0L#T*$pB zVJ%YvCQW8UqUv6hBn37jFeah)VXf^P0mq75G_H$nulG9WBHkIZ@}fgc3b{HBR)z>g zmqm$4t|jHipn5Xx0!6>O)by+Qaz`XOwp(#Az-3rR;I-R|n{6|klcDo1UO9IppKXnL zwC_U~td;#L2t>0i%&_+gaF8rzwT7hHGWn3)DeUkEataf6Lf`#^ax#hU%#xqbqF8^t zKpAOq=A43}Q}*7R&rT=kFTK;QQj!FmUEzUsUK=o<6jY`LyS)O~dh#I_NX@NgM&yl0 zbD3n_>m_3m0vo?#?>AKeg?WE-9CSg4rmHr9*pF zV6NZpT-3pN_LAuA!`xy|#l9en$&=@9c8PeNjJVk&_#O4LyOWIRxr67(T(N|gc8~(% zvHi!~=$z*B`#~VCBYi+4$fJmL>wrd9{ihrDCJ~+=#Fl`bUL-ltKb9KYIbq%h@t;4! zSP;B(ARTk}w%mqC>)VOf8Gir^-qC%+6CmD=80=W(A*35s1TW;E13M%6(c^lk#I!ll zE&>9vYsCneebE4K`XV&RKdXMEjQ0JMN1kQwgP3BEuR+n-j;;5_BI{Zx8?sFchyrK& zPVcj8B5Y^ts@6s}o z-_5Lau|Qp$e3?yvkM;UZKc+W*jWIWu*^ldKHBtxOH6)F)uE_Y9^reVEW5OI0UQUpy z4B!Mkkl8?E)zkMmc)I9}ulLm)Su(SHj7v#QjL^z3oihwEql)KT-8I)*P^C&aHG zx_6lLj5E2|v3ftY!1;uwuO~3~=e)FPOc*MyIw4_v{|>UIHn+6N{A(p?4fSTr~R$sKwc zKO_+fCPn+xbBEHYp_QbHzd?6D9fGkg5lScO4KPX1;V1av_2>rD0aVdOv8}7s0br`{ z=1xM6-_KFR_lBJtYJWqmU4WR+^2=-N4{@~+YyI&6i+vL4o8wBl1NOL?`vTRI)n&CO zGRd*7GnXy0m}M8x@@|KkSq9M0-N=4X2H{_OPNa^Lhy1aAB3liR8GcbOb1aePMt*LV z$Icvjy9>g-@#9nQJP+lnrgfRO``n}x4-mZ`-p*k|%WXYCdpz0iZ|OXC_AXqIfLeIs zCtKvK9kh`F811R+uRjRNJ?|4(>qz9WE2_WS$i|j`XgB7AH(qQo9(PZp@8i{nnEF9i z1_5Y)jMhhS#w=sI@O*)W1-2+2b{Opc8p>9r>-Snnft*bWb(1ji9qd#IJa5M>;A_JJ-_#Q&WbL`Xea$RDi}#(9}YQy+jB+5KP^Z z*r(py%6qEfs#-GheUda#gPOkP8`_J>u>Mw_CV%>@B;$ZHx8dA?)1Nf;p+{T6Eg->a z@C061QrG2+iXlkK+tb^q!K}6&Rq#I8)LUVLxf(-Gf^ z{7a-pTa_5@aXp$)4kyS9v6@iZ3pQ(V%sXwl1kd@Swwytns~hUS&w^I)lHI3lM4ifP zEt)}AMX&T}(E8f_t3guejUG@f zYy+3AMlX`*x`G1)I@7g2x{zPDLVa)Xg9v)>ld5yHtXFHRyY^QkVX5B2YcgZPvlfut z2g1#yVeCNs*;fKF^e&5zZ|{kvo!XJ;LJOi6k|vJR_<5~4=U;4q9pKzAa`^oD%FEsQ zHb9g~bDl|y;q>wUyzW9U0N}G4Ac?)jb^x-itub^`V$99aeUnXuv&+H*yyglKrL?Cnti4n`n!cl5B&dznn#T9Lq*~PM_ zf?ZG|Zfu+iOgjQ|^j=Y8K-PAvZ_xjtJ9~SJ?_p{ZgRUaShS3zlRW%hoMfHV+`$4sB(Id7g0$pNrYr$$ zC_A+FNG`O-^y?J?pnb2Dce4{yk<2I|N1AFrpP=E^?AFy5W{ils!BKDQ#a4wezMm-3gZrO4u z{MkvN%K^^pN30G;_MhXlumc6Sekye6k9;)-?&iNs`EOM4{#z^H5&7>4!lV7Ku7I-q pU$vzT{I3UrzWKj>XKuN`BH(r(P1=9U3<3Y-rB$R#9zT8ge*p=40g?a!