Date   
Re: representatives of bars, homology classes. . . perhaps with better formatting

Dmitriy Morozov
 

Hi Ryan,

Your understanding is correct. I just want to add a couple of small things.

This persistent homology of our square has four H_0 bars that merge into one at length-parameter 2.  And it has a single H_1 bar that starts at length parameter 2, and dies as 2\sqrt{2}.  Here is how I merge this understanding with the Dionysus output, line-by-line:
 
<0> 0: 
<1> 0: 
<2> 0: 
<3> 0: 
 
These lines say the four vertices labelled <0>, <1>, <2>, <3> appear at length parameter 0.  The fact that there's nothing to the right of the colon indicates bars are created at length-parameter zero. 

That's right. If there is nothing to the right of the colon (i.e., if the respective column c (=hfp[i]) is empty), then you have a birth in the filtration, when you add that simplex. It's a little tricky to get the representative cycle that's born (it's much easier to get it when it dies). But in the 0-dimensional case, it's just the vertex.
 
The lines:
 
<0,1> 2: 1 * <0> 0 + 1 * <1> 0
<0,3> 2: 1 * <0> 0 + 1 * <3> 0
<1,2> 2: 1 * <1> 0 + 1 * <2> 0
 
say at length-parameter 2, the three edges <0,1>, <0,3>, <1,2> are created and they merge the four 0-cycles into one. 

Right. The precise meaning of the stuff after the colon is that it's the cycle that was non-zero before, but becomes zero after the addition of the simplex before the colon. Moreover, the cycle is reduced, i.e., it's the oldest such cycle. You can get its birth by taking the maximum birth time, or by looking at hfp.pair(i), which will give you the index of the simplex that gave birth to the cycle.
 
The line:
 
<2,3> 2: 
 
indicates the last edge <2,3> is creating a 1-cycle, at side-length 2.  So this helps me.  It's telling me that generally one has to combine these relations to discover the cycle representatives. 

If you want to find representative of the essential cycles, then you need to do extra work. But if you are only looking for representatives for finite bars, then it's best to just check at the time of death.
 
The remaining lines are relatively degenerate:
 
<0,2> 2.82843: 
<1,3> 2.82843: 
<0,1,2> 2.82843: 1 * <0,1> 2 + 1 * <1,2> 2 + 1 * <0,2> 2.82843
<0,1,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,3> 2.82843
<0,2,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,2> 2 + 1 * <2,3> 2
<1,2,3> 2.82843: 
<0,1,2,3> 2.82843: 1 * <0,1,2> 2.82843 + 1 * <0,1,3> 2.82843 + 1 * <0,2,3> 2.82843 + 1 * <1,2,3> 2.82843
 
This says at length 2\sqrt{2} two new 1-cycles are created and simultaneously destroyed, but then also a 2-cycle is created and simultaneously destroyed, i.e. bars of length 0. 

Right, but the line for <0,2,3> also tells you that you have a bar [2,2.82843], and specifically it gives you a cycle that was born when <2,3> entered the filtration and that got killed when <0,2,3> entered the filtration.
 
I hope this helps.
Dmitriy

representatives of bars, homology classes. . . perhaps with better formatting

Ryan Budney
 

Greetings, 
 
Say I've computed persistent homology of some RIPS complex in Dionysus.  How does one go the extra step to find a cycle representative in a homology class?  i.e. I would like to identify the bar, and determine which linear combination of simplices represent the bar. 
 
In the process of writing this question I likely have answered it.  But it would be nice to get confirmation, and perhaps this will be useful to others.
 
To make this concrete and fairly simple, let me use the vertices of a square of side-length two as my data set.
 
import dionysus as di
 
## generate square as a list
 
A = []
for i in [-1,1]:
    for j in [-1, 1]:
        cv = [i, i]
        cv[1] *= j
        A.append(cv)
 
## convert to numpy array
A = np.array(A, float)
 
## generate the associated RIPS complex
X = di.fill_rips(A, 3, 3.0)
 
## compute persistent homology
hfp = di.homology_persistence(X)
 
## here is where I get a bit fuzzy.  My understanding is hfp is a reduced
## matrix that describes the persistent homology.  But I'm not certain how 
## to read it.  Here is some code of Dmitriy's that I lifted from a 2018 thread.
 
for i,c in enumerate(hfp):
  print(X[i],end=': ')
  print(" + ".join(["%d * %s" % (x.element, X[x.index]) for x in c]))
 
If I execute the above code, I get this as output.
 
<0> 0: 
<1> 0: 
<2> 0: 
<3> 0: 
<0,1> 2: 1 * <0> 0 + 1 * <1> 0
<0,3> 2: 1 * <0> 0 + 1 * <3> 0
<1,2> 2: 1 * <1> 0 + 1 * <2> 0
<2,3> 2: 
<0,2> 2.82843: 
<1,3> 2.82843: 
<0,1,2> 2.82843: 1 * <0,1> 2 + 1 * <1,2> 2 + 1 * <0,2> 2.82843
<0,1,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,3> 2.82843
<0,2,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,2> 2 + 1 * <2,3> 2
<1,2,3> 2.82843: 
<0,1,2,3> 2.82843: 1 * <0,1,2> 2.82843 + 1 * <0,1,3> 2.82843 + 1 * <0,2,3> 2.82843 + 1 * <1,2,3> 2.82843
 
What I'd like to do now is see the bars from this output.  My initial impression is this "hfp" object represents the entire chain complex in reduced form. 
 
This persistent homology of our square has four H_0 bars that merge into one at length-parameter 2.  And it has a single H_1 bar that starts at length parameter 2, and dies as 2\sqrt{2}.  Here is how I merge this understanding with the Dionysus output, line-by-line:
 
<0> 0: 
<1> 0: 
<2> 0: 
<3> 0: 
 
These lines say the four vertices labelled <0>, <1>, <2>, <3> appear at length parameter 0.  The fact that there's nothing to the right of the colon indicates bars are created at length-parameter zero. 
 
The lines:
 
<0,1> 2: 1 * <0> 0 + 1 * <1> 0
<0,3> 2: 1 * <0> 0 + 1 * <3> 0
<1,2> 2: 1 * <1> 0 + 1 * <2> 0
 
say at length-parameter 2, the three edges <0,1>, <0,3>, <1,2> are created and they merge the four 0-cycles into one. 
 
The line:
 
<2,3> 2: 
 
indicates the last edge <2,3> is creating a 1-cycle, at side-length 2.  So this helps me.  It's telling me that generally one has to combine these relations to discover the cycle representatives. 
 
The remaining lines are relatively degenerate:
 
<0,2> 2.82843: 
<1,3> 2.82843: 
<0,1,2> 2.82843: 1 * <0,1> 2 + 1 * <1,2> 2 + 1 * <0,2> 2.82843
<0,1,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,3> 2.82843
<0,2,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,2> 2 + 1 * <2,3> 2
<1,2,3> 2.82843: 
<0,1,2,3> 2.82843: 1 * <0,1,2> 2.82843 + 1 * <0,1,3> 2.82843 + 1 * <0,2,3> 2.82843 + 1 * <1,2,3> 2.82843
 
This says at length 2\sqrt{2} two new 1-cycles are created and simultaneously destroyed, but then also a 2-cycle is created and simultaneously destroyed, i.e. bars of length 0. 
 
Thanks, 

representatives of bars, homology classes

Ryan Budney
 

Greetings, 

Say I've computed persistent homology of some RIPS complex in Dionysus.  How does one go the extra step to find a cycle representative in a homology class?  i.e. I would like to identify the bar, and determine which linear combination of simplices represent the bar. 

In the process of writing this question I likely have answered it.  But it would be nice to get confirmation, and perhaps this will be useful to others.

To make this concrete and fairly simple, let me use the vertices of a square of side-length two as my data set.

import dionysus as di

## generate square as a list

A = []
for i in [-1,1]:
    for j in [-1, 1]:
        cv = [i, i]
        cv[1] *= j
        A.append(cv)

## convert to numpy array
A = np.array(A, float)

## generate the associated RIPS complex
X = di.fill_rips(A, 3, 3.0)

## compute persistent homology
hfp = di.homology_persistence(X)

## here is where I get a bit fuzzy.  My understanding is hfp is a reduced
## matrix that describes the persistent homology.  But I'm not certain how 
## to read it.  Here is some code of Dmitriy's that I lifted from a 2018 thread.

for i,c in enumerate(hfp):
  print(X[i],end=': ')
  print(" + ".join(["%d * %s" % (x.element, X[x.index]) for x in c]))

If I execute the above code, I get this as output.

<0> 0:
<1> 0:
<2> 0:
<3> 0:
<0,1> 2: 1 * <0> 0 + 1 * <1> 0
<0,3> 2: 1 * <0> 0 + 1 * <3> 0
<1,2> 2: 1 * <1> 0 + 1 * <2> 0
<2,3> 2:
<0,2> 2.82843:
<1,3> 2.82843:
<0,1,2> 2.82843: 1 * <0,1> 2 + 1 * <1,2> 2 + 1 * <0,2> 2.82843
<0,1,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,3> 2.82843
<0,2,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,2> 2 + 1 * <2,3> 2
<1,2,3> 2.82843: 
<0,1,2,3> 2.82843: 1 * <0,1,2> 2.82843 + 1 * <0,1,3> 2.82843 + 1 * <0,2,3> 2.82843 + 1 * <1,2,3> 2.82843

What I'd like to do now is see the bars from this output.  My initial impression is this "hfp" object represents the entire chain complex in reduced form. 

This persistent homology of our square has four H_0 bars that merge into one at length-parameter 2.  And it has a single H_1 bar that starts at length parameter 2, and dies as 2\sqrt{2}.  Here is how I merge this understanding with the Dionysus output, line-by-line:

<0> 0:
<1> 0:
<2> 0:
<3> 0: 

These lines say the four vertices labelled <0>, <1>, <2>, <3> appear at length parameter 0.  The fact that there's nothing to the right of the colon indicates bars are created at length-parameter zero. 

The lines:

<0,1> 2: 1 * <0> 0 + 1 * <1> 0
<0,3> 2: 1 * <0> 0 + 1 * <3> 0
<1,2> 2: 1 * <1> 0 + 1 * <2> 0

say at length-parameter 2, the three edges <0,1>, <0,3>, <1,2> are created and they merge the four 0-cycles into one. 

The line:

<2,3> 2: 

indicates the last edge <2,3> is creating a 1-cycle, at side-length 2.  So this helps me.  It's telling me that generally one has to combine these relations to discover the cycle representatives. 

The remaining lines are relatively degenerate:

<0,2> 2.82843:
<1,3> 2.82843:
<0,1,2> 2.82843: 1 * <0,1> 2 + 1 * <1,2> 2 + 1 * <0,2> 2.82843
<0,1,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,3> 2.82843
<0,2,3> 2.82843: 1 * <0,1> 2 + 1 * <0,3> 2 + 1 * <1,2> 2 + 1 * <2,3> 2
<1,2,3> 2.82843:
<0,1,2,3> 2.82843: 1 * <0,1,2> 2.82843 + 1 * <0,1,3> 2.82843 + 1 * <0,2,3> 2.82843 + 1 * <1,2,3> 2.82843

This says at length 2\sqrt{2} two new 1-cycles are created and simultaneously destroyed, but then also a 2-cycle is created and simultaneously destroyed, i.e. bars of length 0. 

Thanks, 

-ryan

Re: Wasserstein or Bottleneck Mapping

Dmitriy Morozov
 

In the case of the bottleneck distance, only the longest edge is uniquely defined -- there are many matchings that satisfy the constraints. There is a pull request [1] that provides a way to get this longest edge, but I haven't merged it. You could use that fork directly, if you like. But I do need to merge it as well.


In the case of the Wasserstein distance, Hera (the underlying library) doesn't expose this functionality. It's on our to-do list to fix it, but currently it's not available.

Dmitriy

On Fri, Jun 14, 2019 at 7:01 AM <K.A.Garside2@...> wrote:
Hi Dmitriy,

Is it possible to obtain, using the library, the optimal mapping of points for either the bottleneck or wasserstein distance?

Thanks!

Wasserstein or Bottleneck Mapping

K.A.Garside2@...
 

Hi Dmitriy,

Is it possible to obtain, using the library, the optimal mapping of points for either the bottleneck or wasserstein distance?

Thanks!

Re: Build Under Windows

Dmitriy Morozov
 

Ah, I see. So it's an issue with the distribution of MinGW. That makes sense, and not something that we can fix.

Thanks again for figuring all this out.
Dmitriy

On Thu, Jun 13, 2019 at 8:26 AM <micheal.shoemaker@...> wrote:

I went back and looked more into the MinGW version issue 8.2.0 issue.  Right now, I'm getting:

CMake Error at bindings/python/pybind11/tools/FindPythonLibsNew.cmake:122 (message):
  Python config failure: Python is 64-bit, chosen compiler is 32-bit

Even though gcc -v says it's x86_64.  This is 8.2 installed from MinGW's main page, but this version of gcc uses thread model win32, which causes all kinds of errors even with gcc 8.1.
If I could find a gcc 8.2 with x86_64 and posix thread model, then I believe it will work fine.

Re: Persistent-based segmentation

Dmitriy Morozov
 

I don't really understand your question. The linked messages give a way to get actual cycles (for example, sequences of edges, not just births/deaths). Do you want to do something else?

But I'm really not sure what you are asking. What is a "border of a cycle"? What are the pixels of the images belonging to a cycle? (Other than getting the cycles as described in the messages I linked.) It's also worth noting that there is no uniquely defined cycle that persistence finds. I think you need to try to formalize mathematically what you want to find, then it would be possible to figure out if it's something Dionysus computes.

I also don't understand the strategy you describe: if you want all pixels between value 5 and 22, you don't need Dionysus for that; NumPy can give you these pretty quickly.

In short, I don't think there is an optimized function in Dionysus to do what you want, besides extracting the representative cycles, as described in the messages I linked originally.

Dmitriy

On Thu, Jun 13, 2019 at 5:20 AM Jack <giacomo@...> wrote:
Hi Dmitriy,

did you have any time to think about my question?

Thanks

Re: Build Under Windows

micheal.shoemaker@...
 

I went back and looked more into the MinGW version issue 8.2.0 issue.  Right now, I'm getting:

CMake Error at bindings/python/pybind11/tools/FindPythonLibsNew.cmake:122 (message):
  Python config failure: Python is 64-bit, chosen compiler is 32-bit

Even though gcc -v says it's x86_64.  This is 8.2 installed from MinGW's main page, but this version of gcc uses thread model win32, which causes all kinds of errors even with gcc 8.1.
If I could find a gcc 8.2 with x86_64 and posix thread model, then I believe it will work fine.

Re: Persistent-based segmentation

Jack <giacomo@...>
 

Hi Dmitriy,

did you have any time to think about my question?

Thanks

Re: Build Under Windows

Dmitriy Morozov
 

This is fantastic. Thank you for these instructions. Building on Windows seems to be the most common question, and now we have this write-up to point people to.

Now if I could figure out how to get the conda script to follow these instructions automatically, it would probably eliminate the questions entirely.

I'm a little confused by the MinGW version issues. Can you say a little more about them? You mention version 8.1.0 and 8.2.0. Do you have a sense what the problem is with MinGW 8.2.0? That's probably worth addressing for future sanity.

Thanks again.
Dmitriy

On Wed, Jun 12, 2019 at 1:46 PM <micheal.shoemaker@...> wrote:
I've spent some time trying to get Dionysus 2 to build under Windows.  Here are my notes on getting it to build.  The versions listed are the ones that I used and I have some notes on versions that didn't work.  I hope this can help some people to get this working.

You'll need these dependencies:

MinGW 8.1.0 (The only other version I tried was 8.2.0, which will error on build.)  Download from: https://sourceforge.net/projects/mingw-w64/
- Install with these settings: Version: 8.1.0, Architecture: x86_64, Threads: posix, Exception: seh, Build version: 0
After installing, Copy C:\Program Files\mingw-w64\x86_64-8.1.0-posix-seh-rt_v6-rev0\mingw64\libexec\gcc\x86_64-w64-mingw32\8.1.0\liblto_plugin-0.dll to C:\Program Files\mingw-w64\x86_64-8.1.0-posix-seh-rt_v6-rev0\mingw64\lib\bfd-plugins

Boost 1.70 (Boost 1.69 requires a change to the CMakeLists.txt in dionysus.) Download from: https://www.boost.org/

Download and extract to C:\boost\
Add Variable 'BOOST_INCLUDE_DIR' with value 'C:\boost' to the environment variables.
This version of boost will give a warning on build:

BOOST_HEADER_DEPRECATED( "the facilities in <boost/timer/timer.hpp>" )

cmake 3.14.5 Download from: https://cmake.org/download/
Just install, no special notes.

Python 3.7.3 Download from: https://www.python.org/
Be sure to check the box to add to PATH

Build
Then, just clone the git (You may need to install git):

git clone https://github.com/mrzv/dionysus.git

Then run from the dionysus folder (from Powershell):

mkdir build
cd build
cmake -G “MinGW Makefiles” ..
 
mingw32-make.exe 

Then, you can copy the folder C:\path\to\dionysus\build\bindings\python\dionysus into your site-packages folder or project folder.

 

 

 

 

Build Under Windows

micheal.shoemaker@...
 

I've spent some time trying to get Dionysus 2 to build under Windows.  Here are my notes on getting it to build.  The versions listed are the ones that I used and I have some notes on versions that didn't work.  I hope this can help some people to get this working.

You'll need these dependencies:

MinGW 8.1.0 (The only other version I tried was 8.2.0, which will error on build.)  Download from: https://sourceforge.net/projects/mingw-w64/
- Install with these settings: Version: 8.1.0, Architecture: x86_64, Threads: posix, Exception: seh, Build version: 0
After installing, Copy C:\Program Files\mingw-w64\x86_64-8.1.0-posix-seh-rt_v6-rev0\mingw64\libexec\gcc\x86_64-w64-mingw32\8.1.0\liblto_plugin-0.dll to C:\Program Files\mingw-w64\x86_64-8.1.0-posix-seh-rt_v6-rev0\mingw64\lib\bfd-plugins

Boost 1.70 (Boost 1.69 requires a change to the CMakeLists.txt in dionysus.) Download from: https://www.boost.org/

Download and extract to C:\boost\
Add Variable 'BOOST_INCLUDE_DIR' with value 'C:\boost' to the environment variables.
This version of boost will give a warning on build:

BOOST_HEADER_DEPRECATED( "the facilities in <boost/timer/timer.hpp>" )

cmake 3.14.5 Download from: https://cmake.org/download/
Just install, no special notes.

Python 3.7.3 Download from: https://www.python.org/
Be sure to check the box to add to PATH

Build
Then, just clone the git (You may need to install git):

git clone https://github.com/mrzv/dionysus.git

Then run from the dionysus folder (from Powershell):

mkdir build
cd build
cmake -G “MinGW Makefiles” ..
 
mingw32-make.exe 

Then, you can copy the folder C:\path\to\dionysus\build\bindings\python\dionysus into your site-packages folder or project folder.

 

 

 

 

Re: Persistent-based segmentation

Dmitriy Morozov
 

I think you'll have to explain better what you are trying to do before I can help you.

You can get cycles (which are not unique) from Dionysus. That question has come up before on this list. For example, see


I'm not sure what you mean by pixels having persistence larger than a given threshold.

Dmitriy

On Wed, Jun 5, 2019 at 9:57 AM Giacomo Nardi <giacomo@...> wrote:

Hi,

I am using the dionysus package to compute persistent homology for images.  I was wondering if there is any function to perform segmentation based on persistence.

How can I select cycles/pixels having persistence larger than a given threshold ?

Thank you
Giacomo

--
Giacomo Nardi | Research Scientist TDI
Mauna Kea Technologies | 9 rue d Enghien – 75010 Paris 
Tel : + 33 1 48 24 03 45 | Fax : + 33 1 48 24 12 18 | www.maunakeatech.com
Cellvizio iPhone App | Follow us on Twitter | Follow us on LinkedIn 

logomaunakeatech.jpg

This message and any attached documents are CONFIDENTIAL and intended only for the use of the intended recipient.
Ce message et tous les documents attachés sont CONFIDENTIELS et à la seule attention du destinataire pour lequel ils ont été écrits.

Re: Difference between R TDA package and python Dionysus2

Dmitriy Morozov
 

The R package uses Dionysus 1, for sure. Dionysus 1 doesn't have any code to construct triangulations on grids: it's left to the user. How the R package implements this step is a question for its authors.


On Fri, May 31, 2019 at 6:23 AM <K.A.Garside2@...> wrote:
The R package states that it uses Dionysus as a library however? I assume that means it is not using Dionysus 2?

" For that, this package provides an R interface for the efficient algorithms of the C++ libraries 'GUDHI' <http://gudhi.gforge.inria.fr/>, 'Dionysus' <http://www.mrzv.org/software/dionysus/>, and 'PHAT' <https://bitbucket.org/phat-code/phat/>. "

Re: Difference between R TDA package and python Dionysus2

K.A.Garside2@...
 

The R package states that it uses Dionysus as a library however? I assume that means it is not using Dionysus 2?

" For that, this package provides an R interface for the efficient algorithms of the C++ libraries 'GUDHI' <http://gudhi.gforge.inria.fr/>, 'Dionysus' <http://www.mrzv.org/software/dionysus/>, and 'PHAT' <https://bitbucket.org/phat-code/phat/>. "

Re: Difference between R TDA package and python Dionysus2

Dmitriy Morozov
 

I don't know what the R package does, but if I had to guess, the two are using different triangulations of the domain. Dionysus 2 uses the Freudenthal triangulation; I don't know what the R package does. 


On Thu, May 30, 2019 at 11:23 AM <K.A.Garside2@...> wrote:
I've just been comparing the output between the Dionysus version in the TDA package in R and the python Dionysus 2 version and there are slightly different results.

I ran it on both sub and super level set filtration on a random field and got slightly different results. Is there a reason for the differences between the R and python implementation?

The 6x6 random field was as follows:

0.424929227922883 0.0128508228081597 -0.294177345297659 -0.49885453849692 -0.442318052990762 -0.626453810742332
0.201060103852361 -0.00595309362438814 -0.286812470518481 -0.845635155782885 -0.328884989662953 -0.539250210004794
0.172361359893672 0.0791446149595664 0.0298949263564899 -0.208043570082747 -0.280335952479328 -0.770728999699832
-0.0263036491145396 -0.158923878454222 0.239839365534883 -0.124288608608425 -0.289240676712963 -0.241020973140624
-0.423771084401562 -0.155098230023893 0.239843914614278 0.0208525497197847 0.14293072464965 -0.127618257296956
-0.371494249200723 -0.0954999907866207 -0.213831055439932 0.268132722688926 0.0128755839751564 -0.374495881349471


R gives the following diagrams:
# Sub level set:
     dimension      Birth       Death
[1,]         0 -0.8456352  0.42492923
[2,]         0 -0.7707290 -0.44231805
[3,]         0 -0.6264538 -0.53925021
[4,]         0 -0.4237711  0.02989493
[5,]         0 -0.3744959 -0.12761826
[6,]         0 -0.2892407 -0.28033595
[7,]         0 -0.2138311 -0.09549999

# Super level set
     dimension       Death      Birth
[1,]         0 -0.84563516  0.4249292
[2,]         0  0.02989493  0.2681327
[3,]         1 -0.28924068 -0.2803360
[4,]         1 -0.84563516 -0.4988545


Whilst the python Dionysus2 version gives:
# Sub level set 
0 -0.845635175704956 Inf
0 -0.7707290053367615 -0.44231805205345154
0 -0.6264538168907166 -0.5392501950263977
0 -0.42377108335494995 0.07914461195468903
0 -0.37449589371681213 -0.12761825323104858
0 -0.28924068808555603 -0.2803359627723694
0 -0.21383105218410492 -0.15509822964668274
0 -0.15892387926578522 -0.15509822964668274
1 0.020852549001574516 0.14293073117733002

# Super level set
0 0.4249292314052582 -Inf
0 0.2681327164173126 0.07914461195468903
0 0.14293073117733002 0.020852549001574516
0 -0.09549999237060547 -0.15509822964668274
1 -0.15509822964668274 -0.15892387926578522
1 -0.2803359627723694 -0.28924068808555603
1 -0.49885454773902893 -0.845635175704956
 

Difference between R TDA package and python Dionysus2

K.A.Garside2@...
 

I've just been comparing the output between the Dionysus version in the TDA package in R and the python Dionysus 2 version and there are slightly different results.

I ran it on both sub and super level set filtration on a random field and got slightly different results. Is there a reason for the differences between the R and python implementation?

The 6x6 random field was as follows:

0.424929227922883 0.0128508228081597 -0.294177345297659 -0.49885453849692 -0.442318052990762 -0.626453810742332
0.201060103852361 -0.00595309362438814 -0.286812470518481 -0.845635155782885 -0.328884989662953 -0.539250210004794
0.172361359893672 0.0791446149595664 0.0298949263564899 -0.208043570082747 -0.280335952479328 -0.770728999699832
-0.0263036491145396 -0.158923878454222 0.239839365534883 -0.124288608608425 -0.289240676712963 -0.241020973140624
-0.423771084401562 -0.155098230023893 0.239843914614278 0.0208525497197847 0.14293072464965 -0.127618257296956
-0.371494249200723 -0.0954999907866207 -0.213831055439932 0.268132722688926 0.0128755839751564 -0.374495881349471


R gives the following diagrams:
# Sub level set:
     dimension      Birth       Death
[1,]         0 -0.8456352  0.42492923
[2,]         0 -0.7707290 -0.44231805
[3,]         0 -0.6264538 -0.53925021
[4,]         0 -0.4237711  0.02989493
[5,]         0 -0.3744959 -0.12761826
[6,]         0 -0.2892407 -0.28033595
[7,]         0 -0.2138311 -0.09549999

# Super level set
     dimension       Death      Birth
[1,]         0 -0.84563516  0.4249292
[2,]         0  0.02989493  0.2681327
[3,]         1 -0.28924068 -0.2803360
[4,]         1 -0.84563516 -0.4988545


Whilst the python Dionysus2 version gives:
# Sub level set 
0 -0.845635175704956 Inf
0 -0.7707290053367615 -0.44231805205345154
0 -0.6264538168907166 -0.5392501950263977
0 -0.42377108335494995 0.07914461195468903
0 -0.37449589371681213 -0.12761825323104858
0 -0.28924068808555603 -0.2803359627723694
0 -0.21383105218410492 -0.15509822964668274
0 -0.15892387926578522 -0.15509822964668274
1 0.020852549001574516 0.14293073117733002

# Super level set
0 0.4249292314052582 -Inf
0 0.2681327164173126 0.07914461195468903
0 0.14293073117733002 0.020852549001574516
0 -0.09549999237060547 -0.15509822964668274
1 -0.15509822964668274 -0.15892387926578522
1 -0.2803359627723694 -0.28924068808555603
1 -0.49885454773902893 -0.845635175704956
 

Re: Simplex Representation

Dmitriy Morozov
 

Well, that's a very weak upper bound, since that number is astronomical. The actual computer memory is a far stronger limitation.


On Mon, Apr 15, 2019 at 11:19 AM tom <gebhart@...> wrote:
This implies the largest number of 0-dimensional simplices Dionysus can handle is the highest integer representation of an unsigned long int? 

Re: Simplex Representation

tom
 

This implies the largest number of 0-dimensional simplices Dionysus can handle is the highest integer representation of an unsigned long int? 

Re: Simplex Representation

Dmitriy Morozov
 

This would be terrible for efficiency. I made a deliberate choice to disallow this. You can always add an extra layer of indirection: store your preferred type in a list and use that list's indices as the vertices for the simplices. 


On Thu, Mar 28, 2019 at 10:10 PM tom <gebhart@...> wrote:
Hi Dimitry,

Would it be difficult to allow for other immutable objects/primitives to be used as Simplex identifiers in the python bindings? It would be super convenient to be able to index simplices by tuples or strings instead of just integers. 

Thanks 

Simplex Representation

tom
 

Hi Dimitry,

Would it be difficult to allow for other immutable objects/primitives to be used as Simplex identifiers in the python bindings? It would be super convenient to be able to index simplices by tuples or strings instead of just integers. 

Thanks