diff --git a/README.md b/README.md index 93af01f..f665f0b 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,4 @@ # MuRiT -Copyright © 2022 Michael Bleher ## Description `MuRiT` is a standalone program that provides Vietoris-Rips transformations for multifiltered flag complexes. @@ -36,8 +35,7 @@ go run main.go [--options] For help in building your own executables, [see here](https://go.dev/doc/tutorial/compile-install). ### Contribute -We envision a broad interaction with the research community. -If you have a special use-case for `MuRiT` that isn't covered by the current implementation or have an idea for additional features, get in touch or submit your pull requests. +If you have a special use-case for `MuRiT` that isn't covered by the current implementation or have an idea for additional features, feel free to get in touch or submit your pull requests. ## Using `MuRiT` @@ -45,7 +43,7 @@ If you have a special use-case for `MuRiT` that isn't covered by the current imp - a comma-separated lower-triangular distance matrix, - a pointwise annotation file, -- and a one-dimensional subfiltration $(VR_0, i_0, j_0, ...)-- ... --(VR_n, i_n, j_n, ...)$. +- and a one-dimensional subfiltration $[VR_0, i_0, j_0, ...]-- ... --[VR_n, i_n, j_n, ...]$. It's output is @@ -55,10 +53,10 @@ It's output is A call to `MuRiT` looks as follows ``` -murit --dist [filename] --pt_fltr [filename] --sub_fltr (VR_0, i_0, j_0, ...)-- ... --(VR_n, i_n, j_n, ...) [--options] +murit --dist --minima --path [VR_0, i_0, j_0, ...]-- ... --[VR_n, i_n, j_n, ...] [--options] ``` -*Note: MuRiT relies on the standard partial order on $\mathbb{R}^n$ for some $n$. Make sure you have translated* +*Note: MuRiT uses the standard partial order on $\mathbb{R}^n$ for some $n$.* ### List of Command Line Arguments @@ -67,11 +65,11 @@ murit --dist [filename] --pt_fltr [filename] --sub_fltr (VR_0, i_0, j_0, ...)-- file name of lower-triangular distance matrix. -##### `--pt_fltr` +##### `--minima` -file name of pointwise filtration annotation. +file name of pointwise minima annotation. - file content: on row 'i' a comma-separated list of minimal filtration values for data point 'i'. + file content: on row `i` a comma-separated list of minimal filtration values for data point $x_i$. example file: @@ -79,12 +77,12 @@ file name of pointwise filtration annotation. (1,1,1) // minimum of point 2 ... -##### `--sub_fltr` +##### `--path` - command line input of sub-filtration along which to compute 1d persistence. + command line input of the path along which to compute 1d persistence. example: - (VR_0, i_0, j_0, k_0,...)-- ... --(VR_n, i_n, j_n, k_n,...) + `[VR_0, i_0, j_0, k_0,...]-- ... --[VR_n, i_n, j_n, k_n,...]` ##### `--threads` number of threads (default: runtime.NumCPU()) @@ -94,9 +92,10 @@ Show status messages (default: false) Show help message ##### `--ripser` run ripser on auxiliary distance matrix (default: false). -Requires a local ripser installation in PATH -Further Ripser options +Requires a local ripser installation in PATH. + +Ripser options (handed directly to ripser) - `--dim` compute persistent homology up to dimension k (default: 1). - `--threshold` compute persistent homology up to threshold t (in auxiliary distance matrix, default: enclosing radius). - `--modulus` compute homology with coefficients in the prime field Z/pZ (default: 2). @@ -110,22 +109,26 @@ _The distance matrix and filtration file for the following example is provided i Consider the metric space of a diamond consisting of four points with side length $1$ and diagonal length $2$. Assume that for each point there is additional information about time $t$ at which the point entered the data set (say in seconds), and some threshold height $h$ above which it could be observed. -The following list of filtration tuples $(t, h)$ for each of the four points is provided in `examples/diamonds.fltr`. +The following list of filtration tuples $(t, h)$ for each of the four points is provided in `examples/diamond.minima`. ``` (0,0) (1.2,0) (1.6,1) -(3.1,1) +(2,1) + ``` This means that the first point was observed from the beginning and already at lowest height, while the second point only appears after 1.2 seconds and also at lowest height, et cetera. Running `MuRiT` on this data produces the distance matrix of a Vietoris-Rips transformation ``` -$ murit --dist examples/diamond.dist --pt_fltr examples/diamond.fltr --verbose -Read pointwise annotation file -Build subfiltration -[[0 0 0] [1 0 0] [1 1.2 0] [1 1.6 1] [1 2 1]] -Build auxiliary Distance Matrix +$ murit --dist examples/diamond.dist --minima examples/diamond.minima --verbose +Pointwise minima +[[[0 0]] [[1.2 0]] [[1.6 1]] [[2 1]]] + +Path +[0,0,0]--[1,0,0]--[1,1.2,0]--[1,1.6,1]--[1,2,1] + +Auxiliary Distance Matrix 3 4,5 5,5,5 @@ -134,31 +137,34 @@ _Note: If no subfiltration is specified, `MuRiT` will build some auxiliary initi If we add the `--ripser` flag, the distance matrix is automatically handed over to `Ripser` to calculate persistent homology along the one-dimensional subfiltration. ``` -$ murit --dist examples/diamond.dist --pt_fltr examples/diamond.fltr --verbose --ripser -Read pointwise annotation file -Build subfiltration -[[0 0 0] [1 0 0] [1 1.2 0] [1 1.6 1] [1 2 1]] -Build auxiliary Distance Matrix ---- -Run Ripser +$ murit --dist examples/diamond.dist --minima examples/diamond.minima --verbose --ripser +Pointwise minima +[[[0 0]] [[1.2 0]] [[1.6 1]] [[2 1]]] + +Path +[0,0,0]--[1,0,0]--[1,1.2,0]--[1,1.6,1]--[1,2,1] + +Ripser value range: [3,5] distance matrix with 4 points, using threshold at enclosing radius 5 persistent homology intervals in dim 1: - ``` Note that we do not find non-trivial homology in this example. This makes sense, because the last point enters the filtration only at the last step in the one-dimensional subfiltration. Only from that point on we might expect to see the cycle around the diamond. -Adding two further filtration steps $[2,2,1]--[2.1,2,1]$ at higher Vietoris-Rips parameters to the automatically generated subfiltration, we get: +Adding two further filtration steps `[2,2,1]--[2.1,2,1]` at higher Vietoris-Rips parameters to the automatically generated subfiltration, we get: ``` -$ murit --dist examples/diamond.dist --pt_fltr examples/diamond.fltr --verbose --ripser --sub_fltr [0,0,0]--[1,0,0]--[1,1,0]--[1,2,0]--[1,2,1]--[2,2,1]--[2.1,2,1] -Read pointwise annotation file -Read subfiltration -[[0 0 0] [1 0 0] [1 1 0] [1 2 0] [1 2 1] [2 2 1] [2.1 2 1]] -Build auxiliary Distance Matrix ---- -Run Ripser +$ murit --dist examples/diamond.dist --minima examples/diamond.minima --verbose --ripser + --path [0,0,0]--[1,0,0]--[1,1,0]--[1,2,0]--[1,2,1]--[2,2,1]--[2.1,2,1] + +Pointwise minima +[[[0 0]] [[1.2 0]] [[1.6 1]] [[2 1]]] + +Path +[0,0,0]--[1,0,0]--[1,1,0]--[1,2,0]--[1,2,1]--[2,2,1]--[2.1,2,1] + +Ripser value range: [4,6] distance matrix with 4 points, using threshold at enclosing radius 6 persistent homology intervals in dim 1: @@ -166,11 +172,11 @@ persistent homology intervals in dim 1: {[0,2] (5), [1,3] (5), [2,3] (5), [0,1] (4)} {[0,2] (5), [1,3] (5), [2,3] (5), [0,1] (4)} ``` -Indeed we find that there is a homology class in dimension one, which is born at filtration value $[1,2,1]$, i.e. once the fourth point has been observed. -The homology class dies at filtration value $[2,2,1]$, which is when the diagonals are added to the Vietoris-Rips complex. +Indeed we find that there is a homology class in dimension one, which is born at filtration value `[1,2,1]`, i.e. once the fourth point has been observed. +The homology class dies at filtration value `[2,2,1]`, which is when the diagonals are added to the Vietoris-Rips complex. -_Note: By using the `--ripser` flag, the result is automatically converted back to the subfiltration values._ +_Note: By using the `--ripser` flag, the result is automatically converted back to filtration values._ ## Citing diff --git a/examples/diamond.fltr b/examples/diamond.minima similarity index 100% rename from examples/diamond.fltr rename to examples/diamond.minima diff --git a/examples/test.fltr b/examples/test.minima similarity index 100% rename from examples/test.fltr rename to examples/test.minima diff --git a/main.go b/main.go index cca47fc..7fbfc4f 100644 --- a/main.go +++ b/main.go @@ -25,8 +25,8 @@ type workload struct { } type args struct{ - sub_fltr [][]float64 - fltr_list [][]float64 + path [][]float64 + minima_list [][][]float64 } @@ -63,6 +63,21 @@ func equal(poset_element1 []float64, poset_element2 []float64) bool { } +// find index for which a data point first enters a totally ordered subfiltration from list of filtration minima +func get_index_of_entry(minima [][]float64, path [][]float64) int{ + // for each filtration step along the subfiltration, check if one of the minima of the data point lies below. + // If so, we found the point of entry + for k, x := range path { + for _, minimum := range minima{ + if leq(minimum, x){ + return k+1 // shift from zero- to one-indexing! + } + } + // if the point is not contained in the path at all, set index to (maximal filtration value + 1) <- equivalent to infty + } + return len(path) // shift from zero- to one-indexing! +} + // // Reader function @@ -115,17 +130,28 @@ func worker(in chan workload, out chan string, b args) { if err != nil { log.Fatalf("Distance conversion error: %v", err) } - // Goal: Determine modified distance for pair of datapoints (x_i,x_j), see article. - // find sub-filtration point in which the given edge is present and set distance to that value. - modified_distance := len(b.sub_fltr) - fltr_point_i := append([]float64{distance}, b.fltr_list[w.i]...) - fltr_point_j := append([]float64{distance}, b.fltr_list[j]...) - for k, x := range b.sub_fltr { - if leq(fltr_point_i, x) && leq(fltr_point_j, x) { - modified_distance = k+1 // adjusting from zero- to one-indexing to get a well-behaved semi distance matrix - break - } - } + // Append distance of x_i and x_j to the start of each filtration value of the given minimum + var minima_i [][]float64 + for _, minimum := range b.minima_list[w.i]{ + minima_i = append(minima_i, append([]float64{distance}, minimum...)) + } + var minima_j [][]float64 + for _, minimum := range b.minima_list[j]{ + minima_j = append(minima_j, append([]float64{distance}, minimum...)) + } + // Determine modified distance for pair of datapoints (x_i,x_j). + // modified distance is the index of entry for the edge (x_i,x_j) + // the edge is present as soon as both x_i and x_j have entered + var modified_distance int + index_of_entry_i := get_index_of_entry(minima_i, b.path) + index_of_entry_j := get_index_of_entry(minima_j, b.path) + // use maximum of the two + if index_of_entry_i >= index_of_entry_j { + modified_distance = index_of_entry_i + } else { + modified_distance = index_of_entry_j + } + // fmt.Println(minima_i, index_of_entry_i, "-", minima_j, index_of_entry_j, "-", modified_distance) // Concatenate modified distance to current matrix line sb.WriteString(strconv.Itoa(modified_distance)) // write modified distance sb.WriteByte(',') // write separator @@ -190,8 +216,8 @@ func hash(s string) string { //------------------------------------------------------------------------- func main() { var dist_file_name string - var fltr_file_name string - var sub_fltr_input string + var minima_file_name string + var path_input string var threads string var verbose bool @@ -213,24 +239,24 @@ func main() { flag.StringVar(&dist_file_name, "dist", "", "file name of lower-triangular distance matrix.") - aux_description=`file name of pointwise filtration annotation. + aux_description=`file name of pointwise minima annotation. file content: on row 'i' a comma-separated list of minimal filtration values for data point 'i'. standard partial order on R^n. example: - (0,0,1) // minimum of point 1 - (1,1,1) // minimum of point 2 + (0,0,1), (1,0,0) // minima of point 1 + (1,1,1) // minima of point 2 ... ` - flag.StringVar(&fltr_file_name, "pt_fltr", "", aux_description) + flag.StringVar(&minima_file_name, "minima", "", aux_description) aux_description=`command line input of sub-filtration along which to compute 1d persistence. example: - (VR_0, i_0, j_0, k_0,...)-- ... --(VR_n, i_n, j_n, k_n,...) + [VR_0, i_0, j_0, k_0,...]-- ... --[VR_n, i_n, j_n, k_n,...] ` - flag.StringVar(&sub_fltr_input, "sub_fltr", "", aux_description) + flag.StringVar(&path_input, "path", "", aux_description) flag.StringVar(&threads, "threads", "", "number of threads (default: runtime.NumCPU())") flag.BoolVar(&verbose, "verbose", false, "Show status messages (default: false)") @@ -246,11 +272,11 @@ func main() { flag.Usage = func() { flagSet := flag.CommandLine aux_description = `Usage: -murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- ... --(VR_n, i_n, j_n, ...) [--options] +murit --dist --minima --path (VR_0, i_0, j_0, ...)-- ... --(VR_n, i_n, j_n, ...) [--options] ` fmt.Printf(aux_description) fmt.Printf("\nCommand Line Arguments\n") - arguments := []string{"dist", "pt_fltr", "sub_fltr", "threads", "verbose", "help", "ripser"} + arguments := []string{"dist", "minima", "path", "threads", "verbose", "help", "ripser"} for _, name := range arguments { flag := flagSet.Lookup(name) fmt.Printf("-%s\n", flag.Name) @@ -277,42 +303,49 @@ murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- if dist_file_name == "" { log.Fatal("dist file name required (--dist_file)") } - if fltr_file_name == "" { - log.Fatal("filtration file name required (--fltr_file)") + if minima_file_name == "" { + log.Fatal("filtration file name required (--minima_file)") } // Read filtration file - if verbose {fmt.Println("Read pointwise annotation file")} - fltr_file, err := os.Open(fltr_file_name) + if verbose {fmt.Println("Pointwise minima")} + minima_file, err := os.Open(minima_file_name) if err != nil { - log.Fatalf("Failed to open file '%s': %v", fltr_file_name, err) + log.Fatalf("Failed to open file '%s': %v", minima_file_name, err) } - var fltr_list [][]float64 - fltr_scanner := bufio.NewScanner(fltr_file) - for fltr_scanner.Scan() { - var row []float64 - for _, v := range strings.Split(fltr_scanner.Text(),",") { - v = strings.Trim(v, " ()") - s, err := strconv.ParseFloat(v,64) + var minima_list [][][]float64 + minima_scanner := bufio.NewScanner(minima_file) + for minima_scanner.Scan() { + var minima [][]float64 + // assume each line is a comma-separated list of minima in the format (a_1, a_2, ...), (b_1, b_2, ...) + line := minima_scanner.Text() + // split lines into separate minima (a_1, a_2, ...) + for _, x := range strings.Split(line, "),("){ + // convert the separated minimum into a list of floats value, by value + var minimum []float64 + for _, value_str := range strings.Split(strings.Trim(x, " ()"), ",") { + value_float, err := strconv.ParseFloat(value_str, 64) if err != nil { log.Fatalf("Filtration parse error: %v", err) } - row = append(row, s) + minimum = append(minimum, value_float) + } + minima = append(minima, minimum) } - fltr_list = append(fltr_list, row) + minima_list = append(minima_list, minima) } + if verbose {fmt.Println(minima_list)} // Close filtration file - fltr_file.Close() + minima_file.Close() // Read sub filtration from command line OR create default sub filtration - var sub_fltr [][]float64 - if sub_fltr_input != "" { + var path [][]float64 + if path_input != "" { // Parse sub filtration from command line input - if verbose {fmt.Println("Read subfiltration")} - for _, p := range strings.Split(sub_fltr_input,"--"){ + for _, p := range strings.Split(path_input,"--"){ var point []float64 for _, q := range strings.Split(strings.Trim(p, " []"),","){ s, err := strconv.ParseFloat(q,64) @@ -321,41 +354,51 @@ murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- } point = append(point, s) } - sub_fltr = append(sub_fltr, point) + path = append(path, point) } // Check if sub filtration is valid (i.e. in lexicographical order) - for i, j := 0, 1; j < len(sub_fltr); i, j = i+1, j+1 { - if !(leq(sub_fltr[i], sub_fltr[j])) { - log.Fatalf("Invalid sub-filtration: sub_fltr[%v] = %v !<= %v = sub_fltr[%v]", i, sub_fltr[i], sub_fltr[j], j) + for i, j := 0, 1; j < len(path); i, j = i+1, j+1 { + if !(leq(path[i], path[j])) { + log.Fatalf("Invalid Path: path[%v] = %v !<= %v = path[%v]", i, path[i], path[j], j) } } } else { // Extract a valid sub filtration from the filtration file (traverse through list once and successively add larger elements) - if verbose {fmt.Println("Build subfiltration")} - sub_fltr = append(sub_fltr, append([]float64{0}, fltr_list[0]...)) - for _, x := range fltr_list { - x = append([]float64{1}, x...) - if leq(sub_fltr[len(sub_fltr)-1], x) && !equal(sub_fltr[len(sub_fltr)-1], x) { - sub_fltr = append(sub_fltr, x) + path = append(path, append([]float64{0}, minima_list[0][0]...)) + for _, x := range minima_list { + if leq(path[len(path)-1], append([]float64{1}, x[0]...)) && !equal(path[len(path)-1], append([]float64{1}, x[0]...)) { + path = append(path, append([]float64{1}, x[0]...)) } } } + if verbose { + fmt.Println("\nPath") + outer_sep := "" + for _, fltr_point := range path { + inner_sep := "" + fmt.Print(outer_sep,"[") + for _, value := range fltr_point{ + fmt.Print(inner_sep, value) + inner_sep = "," + } + fmt.Print("]") + outer_sep = "--" + } + fmt.Print("\n") + } - if verbose {fmt.Println(sub_fltr)} - // // Prepare Input and Communication Channels // - if verbose {fmt.Println("Build auxiliary Distance Matrix")} - // For future development: - // Set filename of auxiliary distance matrix in dependence of sub_fltr - // aux_file_name = filepath.Dir(dist_file_name)+"/"+hash(sub_fltr_input)+".aux" + // In future development: + // Set filename of auxiliary distance matrix in dependence of path + // aux_file_name = filepath.Dir(dist_file_name)+"/"+hash(path_input)+".aux" aux_file_name = filepath.Dir(dist_file_name)+"/"+strconv.FormatInt(time.Now().UTC().UnixNano(), 10)+".aux" // Concatenate background information from above for workers - b := args{sub_fltr, fltr_list} + b := args{path, minima_list} // Open distance matrix file in_file, err := os.Open(dist_file_name) @@ -392,6 +435,7 @@ murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- if ripser { out_writer = bufio.NewWriter(aux_file) } else { + if verbose {fmt.Println("\nAuxiliary Distance Matrix")} out_writer = bufio.NewWriter(os.Stdout) } @@ -420,8 +464,7 @@ murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- // Run ripser on auxiliary distance matrix var ripser_output []byte if ripser { - if verbose {fmt.Println("---")} - if verbose {fmt.Println("Run Ripser")} + if verbose {fmt.Println("\nRipser")} ripser_arguments := []string{"--format", "lower-distance"} if ripser_dim != ""{ @@ -474,17 +517,17 @@ murit --dist --pt_fltr --sub_fltr (VR_0, i_0, j_0, ...)-- if err != nil { log.Fatalf("parsing error: %v", err) } - //adjust from one- to zero-indexing + // shift from one- to zero-indexing birth = birth-1 death, err := strconv.Atoi(x[1]) if err != nil { log.Fatalf("parsing error: %v", err) } - //adjust from one- to zero-indexing + // shift from one- to zero-indexing death = death-1 - fmt.Println(" [", sub_fltr[birth], ",", sub_fltr[death], "):") + fmt.Println(" [", path[birth], ",", path[death], "):") } else { fmt.Println(line) }