paths work

tfields/core.py

@@ -331,9 +331,10 @@ class AbstractNdarray(np.ndarray):
 
             Simply give the file name to save
             >>> p.save(out_file.name)
-            >>> _ = out_file.seek(0)
+            >>> _ = out_file.seek(0)  # this is only necessary in the test
             >>> p1 = tfields.Points3D.load(out_file.name)
             >>> assert p.equal(p1)
+            >>> assert p.coord_sys == p1.coord_sys
 
             The fully nested structure of a TensorMaps object is reconstructed
             >>> out_file_maps = NamedTemporaryFile(suffix='.npz')
@@ -2160,7 +2161,138 @@ class TensorMaps(TensorFields):
 
         """
         maps_list = tfields.lib.sets.disjoint_group_indices(self.maps[mp_idx])
-        return (0, maps_list)
+        return (mp_idx, maps_list)
+
+    def paths(self, mp_idx):
+        """
+        Find the minimal amount of graphs building the original graph with
+        maximum of two links per node i.e. 
+
+            o-----o                       o-----o         
+             \   /                         \   /          
+              \ /                           \ /           
+            o--o--o            o--o          8--o         
+               |                  |
+               |        =         |    +          +       
+               o                  o                    o  
+              / \                /                      \ 
+             /   \              /                        \
+            o     o            o                          o
+
+        where 8 is a duplicated node (one has two links and one has only one.)
+        Examples:
+            >>> import tfields
+            >>> import numpy as np
+            
+            Ascii figure above:
+            >>> a = tfields.TensorMaps([[1, 0], [3, 0], [2, 2], [0, 4], [2, 4],
+            ...                         [4, 4], [1, 6], [3, 6], [2, 2]],
+            ...                        maps=[[[0, 2], [2, 4], [3, 4], [5, 4],
+            ...                               [1, 8], [6, 4], [6, 7], [7, 4]]])
+
+            >>> paths = a.paths(0)
+            >>> assert paths[0].equal([[ 1.,  0.],
+            ...                        [ 2.,  2.],
+            ...                        [ 2.,  4.],
+            ...                        [ 0.,  4.]])
+            >>> assert paths[0].maps[0].equal([[ 0.,  1.,  2.,  3.]])
+            >>> assert paths[1].equal([[ 4.,  4.],
+            ...                        [ 2.,  4.],
+            ...                        [ 1.,  6.],
+            ...                        [ 3.,  6.],
+            ...                        [ 2.,  4.]])
+            >>> assert paths[2].equal([[ 3.,  0.],
+            ...                        [ 2.,  2.]])
+
+        Note:
+            The Longest path problem is a NP-hard problem.
+
+        """
+        obj = self.cleaned()
+        
+        flat_map = np.array(obj.maps[mp_idx].flat)
+        values, counts = np.unique(flat_map, return_counts=True)
+        counts = {v: n for v, n in zip(values, counts)}
+
+        # last is a helper
+        last = np.full(max(flat_map) + 1, -3, dtype=int)
+        duplicat_indices = []
+        d_index = len(obj)
+        for i, val in enumerate(flat_map.copy()):
+            if counts[val] > 2:
+                # The first two occurences are uncritical
+                if last[val] < -1:
+                    last[val] += 1
+                    continue
+
+                # Now we talk about nodes with more than two edges
+                if last[val] == -1:
+                    # append a point and re-link
+                    duplicat_indices.append(val)
+                    flat_map[i] = d_index
+                    last[val] = d_index
+                    d_index += 1
+                else:
+                    # last occurence of val was a duplicate, so we use the same
+                    # value again.
+                    flat_map[i] = last[val]
+                    last[val] = -1
+
+        if duplicat_indices:
+            duplicates = obj[duplicat_indices]
+            obj = type(obj).merged(obj, duplicates)
+        obj.maps = [tfields.Tensors(flat_map.reshape(-1, *obj.maps[mp_idx].shape[1:]))]
+        paths = obj.parts(obj.disjoint_map(0))
+        # paths = [paths[2]]  # this path did not work previously - debugging # with vessel
+
+        # remove duplicate map entries and sort 
+        sorted_paths = []
+        for path in paths:
+            # find start index
+            values, counts = np.unique(path.maps[0].flat, return_counts=True)
+
+            first_node = None
+            for v, c in zip(values, counts):
+                if c == 1:
+                    first_node = v
+                    break
+            edges = [list(edge) for edge in path.maps[0]]
+            if first_node is None:
+                first_node = 0 # edges[0][0]
+                path = path[list(range(len(path))) + [0]]
+                found_first_node = False
+                for edge in edges:
+                    if first_node in edge:
+                        if found_first_node:
+                            edge[edge.index(first_node)] = len(path) - 1
+                            break
+                        else:
+                            found_first_node = True
+
+            # follow the edges until you hit the end
+            chain = [first_node]
+            visited = set()
+            n_edges = len(edges)
+            node = first_node
+            while len(visited) < n_edges:
+                for i, edge in enumerate(edges):
+                    if i in visited:
+                        continue
+                    if node not in edge:
+                        continue
+
+                    # found edge
+                    visited.add(i)
+                    if edge.index(node) != 0:
+                        edge = list(reversed(edge))
+                    chain.extend(edge[1:])
+                    node = edge[-1]
+            path = path[chain]
+            path.maps[0] = tfields.TensorFields([sorted(chain)])
+            sorted_paths.append(path)
+
+        paths = sorted_paths
+        return paths
 
 
 class Container(AbstractNdarray):
@@ -2199,6 +2331,21 @@ class Container(AbstractNdarray):
 
 
 if __name__ == '__main__':  # pragma: no cover
+    # a = tfields.TensorMaps([[1, 0], [3, 0], [2, 2], [0, 4], [2, 4],
+    #                          [4, 4], [1, 6], [3, 6], [2, 2]],
+    #                         maps=[[[0, 2], [2, 4], [3, 4], [5, 4],
+    #                                [1, 8], [6, 4], [6, 7], [7, 4]]])
+    # a = tfields.TensorMaps([[0, 0], [1, 1], [2, 2], [0, 0], [4, 4], [5, 5], [0, 0]],
+    #                         maps=[[[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]])
+    # a = tfields.TensorMaps([[0, 0], [1, 1], [2, 2], [0, 0], [4, 4], [5, 5], [0, 0],
+    #                         [3, 3], [8, 8], [5, 5], [4, 4], [5,5],[8,8]],
+    #                         maps=[[[1, 2], [2, 3], [3, 4], [0, 1], [4, 5], [7, 8], [9, 10], [11, 12]]])
+    # paths = a.paths(0)
+    # print([p for p in paths])
+    # print([p.maps for p in paths])
+
+
+    # quit()
     import doctest
     doctest.testmod()
     # doctest.run_docstring_examples(TensorFields.__getitem__, globals())