Bug in mapper2D$adjacency

[pcamara]

Hi,

There seems to be a bug in mapper2D$adjacency (and I did not check mapper1D$adjacency; potentially it is also there). For instance, consider the following example:

m2 <- mapper2D(
    distance_matrix = dist(data.frame( x=2*cos(1:100), y=sin(1:100) )),
    filter_values = list( 2*cos(1:100), sin(1:100) ),
    num_intervals = c(15,15),
    percent_overlap = 70,
    num_bins_when_clustering = 10)

and let us look at the first 7 vertices:

m2$points_in_vertex[1:7]

[[1]]
[1]  4 48 92

[[2]]
[1] 29 73

[[3]]
[1]  4 48 92

[[4]]
[1] 23 67

[[5]]
[1] 29 73

[[6]]
[1]  4 48 92

[[7]]
[1] 23 67

There should be edges at (1,3), (1,6), (3,6), (2,5), and (4,7). However, m2$adjacency misses the (1,6) edge:

m2$adjacency[1:7,1:7]

      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    1    0    0    0    0
[2,]    0    0    0    0    1    0    0
[3,]    1    0    0    0    0    1    0
[4,]    0    0    0    0    0    0    1
[5,]    0    1    0    0    0    0    0
[6,]    0    0    1    0    0    0    0
[7,]    0    0    0    1    0    0    0

The following code I think would correctly compute the adjacency matrix and can be potentially useful to fix the bug (although I am sure there are more efficient ways to do it):

library(Matrix)
adjacency <- function(m2) {
  l <- length(m2$points_in_vertex)
  simps <- Matrix(0, l, l, sparse=TRUE)
  simps[lower.tri(simps, diag=FALSE)] <- as.numeric(
    lapply(apply(as.matrix(combn(l,2)), 2, 
                 function(x) {c(m2$points_in_vertex[x[1]], m2$points_in_vertex[x[2]])}),
           function(x) {length(intersect(x[[1]], x[[2]])) > 0}))
  return(as.matrix(simps+t(simps)))
}

adjacency(m2)[1:7,1:7]

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    1    0    0    1    0
[2,]    0    0    0    0    1    0    0
[3,]    1    0    0    0    0    1    0
[4,]    0    0    0    0    0    0    1
[5,]    0    1    0    0    0    0    0
[6,]    1    0    1    0    0    0    0
[7,]    0    0    0    1    0    0    0

By the way, thanks for this repository. It is really useful.

Best,

Pablo

[paultpearson]

Hi Pablo,

Are you using the version of TDAmapper on github or CRAN? Given that you're filing a github issue, I think it's likely that you're using the github version, but would like to verify that before looking at fixing the issue. (The CRAN version is likely out of date and may have more bugs, but the github version is the best version available.)

Thanks!

Paul

[pcamara]

Hi Paul,

I am using the Github version.

Thanks,

Pablo

[pcamara]

Sorry for the multiple emails... This is a much more efficient implementation of the function adjacency() using Rcpp. This is the content of the file adjacency.cpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix adjacencyCpp(List x) {
  List xlist(x);
  int n = xlist.size();
  NumericMatrix res(n,n);  

  for(int i=0; i<n-1; i++) {
    NumericVector loc1 = xlist[i];
    for(int j=i+1; j<n; j++) {
      NumericVector loc2 = xlist[j];
      if (intersect(loc1,loc2).length() > 0) {
        res(i,j) = 1;
        res(j,i) = 1;
      };
    };
  };
  return res;
}

Then, in R:

library(Rcpp)
sourceCpp('adjacency.cpp')

adjacencyCpp(m2$points_in_vertex)[1:7,1:7]

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    1    0    0    1    0
[2,]    0    0    0    0    1    0    0
[3,]    1    0    0    0    0    1    0
[4,]    0    0    0    0    0    0    1
[5,]    0    1    0    0    0    0    0
[6,]    1    0    1    0    0    0    0
[7,]    0    0    0    1    0    0    0

[peekxc]

Correct me if I'm wrong on any of this, but I think perhaps this is actually not due to the computation of the adjacency matrix, but in the admissibility criteria of which level sets are considered? It looks like maybe here only adjacent level sets are considered, which works when the overlap is < 50%.

You can kind of see that here: m2$level_of_vertex[1:7] [1] 1 1 2 2 2 3 3

The nodes pairs (1,3), (3,6), (2,5), and (4,7) are all in adjacent level sets: t(sapply(m1$level_of_vertex[1:7], function(lsfi) TDAmapper:::lsmi_from_lsfi(lsfi = lsfi, num_intervals = c(15, 15))))

     [,1] [,2]
[1,]    1    1
[2,]    1    1
[3,]    2    1
[4,]    2    1
[5,]    2    1
[6,]    3    1
[7,]    3    1

but the pair (1, 6) is not.

I have a pull request that (implicitly) takes care of this issue. Using the code from that PR:

m2 <- mapper_new(
  distance_matrix = dist(data.frame( x=2*cos(1:100), y=sin(1:100) )),
    filter_values = list( 2*cos(1:100), sin(1:100) ),
    num_intervals = c(15,15),
    percent_overlap = 70,
    num_bins_when_clustering = 10
)
m2$adjacency[1:7, 1:7]

I get the following. As you can see, (1, 6) are linked.

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    1    0    0    1    0
[2,]    0    0    0    0    1    0    0
[3,]    1    0    0    0    0    1    0
[4,]    0    0    0    0    0    0    1
[5,]    0    1    0    0    0    0    0
[6,]    1    0    1    0    0    0    0
[7,]    0    0    0    1    0    0    0