Parallelization strategies for extraction of (co)homological information of digital objects.

Programa de Doctorado “Matem´aticas” PhD Dissertation Parallelization strategies for extraction of (co)homological information of digital objects Author Ra´ul Reina Molina Supervisor Prof. Dr. Pedro Real Jurado Co-supervisor Dr. Daniel D´ıaz Pernil July 5, 2017 Acknowledgments The work and dedication that requires a work like the one presented in this dissertation cannot be faced alone. From these lines I want to thank all the people who have helped me in some way or another. Not only their technical or theoretical help is important, although I have received it from my supervisors Pedro Real and Daniel D´ıaz, from other fellow doctoral students such as Aldo Gonz´alez Lorenzo or from other doctors like Pawe l Pilarczyk without whose help it would have been impossible to adapt the algorithm suitable for the calculation of Smith’s Normal Form; but it is equally inestimable the emotional help, understanding and effort made by all those who have understood my absences and dedication to this work. In particular, I want to thank the understanding received by my wife without whose help I would not have been able to face this thesis. Last but not least, I can not miss the opportunity to apologize to all those who have received less attention from the accustomed and especially my son born in the course of the development of this dissertation. To all of them, thank you very much. Contents Contents i List of Figures iii 1 Introduction 1 2 Overview of algebraic topology 3 2.1 GeneralTopology........................ 3 2.2 Algebraic Topology . . . . . . . . . . . . . . . . . . . . . . . 5 3 Parallelization strategy for homological calculation based on membrane computing 25 3.1 Membrane computing . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Encoding images as cubical complexes and cubes as tuples . 30 3.3 Membrane computing implementation of homological calculation............................... 32 3.4 Conclusions of the chapter . . . . . . . . . . . . . . . . . . . 45 4 Parallel calculation of an AM-model for a nD digital object 47 4.1 Parallel calculation of GVF in cubical complexes . . . . . . . 47 4.2 Parallel calculation of an AM-model for a cubical complex . 49 4.3 Conclusions ........................... 55 5 Advanced (co)homological information extraction from digital objects 57 5.1 Introduction........................... 57 5.2 (Co)homology information from AM-model . . . . . . . . . . 58 5.3 Cohomology operations from AM-model . . . . . . . . . . . 60 5.4 Conclusions to the chapter . . . . . . . . . . . . . . . . . . . 61 6 Implementation and experimentation 63 6.1 Implementation......................... 63 i 6.2 Experimentation . . . . . . . . . . . . . . . . . . . . . . . . 68 6.3 Empirical complexity analysis . . . . . . . . . . . . . . . . . 76 6.4 Conclusion............................ 78 Bibliography 79 ii List of Figures 2.1 Example of cubical complex . . . . . . . . . . . . . . . . . . . . 13 2.2 Discrete vector field . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1 Example of homology group calculation in a cubical complex. . 39 6.1 Simplicial structure on the Klein bottle . . . . . . . . . . . . . . 70 6.2 Acyclic vector field on the Klein bottle . . . . . . . . . . . . . . 70 6.3 Acyclic vector field on the Klein bottle from an AM-model . . . 71 6.4 Cubical complex decomposition of Bing’s house . . . . . . . . . 74 6.5 Cubical decomposition of torus. . . . . . . . . . . . . . . . . . . 75 6.6 Homology generators at dimension 1 . . . . . . . . . . . . . . . 75 6.7 Empirical complexity graph for 2, 3 and 4 dimensional complexes 77 iii Abstract Digital objects are finite subsets of n-xels within a n-dimensional digital image. The study of the connectivity of these objects, interpreted from a discrete, subdivided or continuous way, has been a priority issue from the very beginning of Digital Imagery. The topological tools that have been exhaustively used in this setting are the notions of connected component, simple point and Euler characteristic. Others topological invariants with a recent increasing popularity are (co)homology groups of digital objects treated as cell complexes. In this thesis, we propose parallelization strategies based on extraction methods of (co)homological information, like Discrete Morse Theory, Effective Homology or AT-model. The first approach is related with Natural Computing, which is a fruitful research area that provides interesting approaches to computational problems inspiring in the way that Nature “computes”. Concretely, we use Membrane Computing, which summarizes with computational rules, the manner living cells work. This area has provided interesting results in theoretical and applied works. The application of these ideas to the process of calculating (co)homology groups of digital objects allows us to develop better algorithms as it brings a natural parallelization of the algorithms implemented in Membrane Computing. Nowadays, there is no current device capable of executing Membrane Computing algorithms, hence the previous processes need to be adapted to be executed by ordinary computing devices. Therefore, in this thesis we present a set of algorithms that provides a compact representation, optimal in some way, of a digital object along with a bidirectional transformation that allows us to compute not only the (co)homology groups but compute some algebraic invariants or operation involving (co)homology classes which can be used as intrinsic information of the digital object. The work presented in this thesis focuses in two main contributions. The first of all is related with Natural Computing. We present a Membrane Computing framework used to make easier the development of Membrane Computing algorithms in Computational Algebraic Topology. This framework is strongly connected with Discrete Morse Theory. The second main contribution is the application of the framework mentioned above for developing a parallel algorithm used to compute a reduction from a cubical cell complex to a CW complex with a v 2.2 Algebraic Topology To find out whether two topological spaces are homeomorphic is sometimes an impossible task. It is preferable to find some properties that remain unchanged under suitable maps and such that they allow to differentiate both spaces. In other words, if both spaces have different values in some property then it can be said that those spaces cannot be the same. Algebraic Topology studies how an algebraic object (called algebraic invariant) can be assigned to any topological space such that homotopy equivalent spaces have isomorphic objects. Therefore, this algebraic property of spaces allows to differentiate them. One of very first algebraic invariant is the fundamental group of a topological space. The fundamental group The fundamental group of a space Xwill be defined so that its elements are loops in Xstarting and ending at a fixed basepoint x0∈X, but two such loops are regarded as being the same if one loop can be continuously deformed to the other within the space X. (All loops that occur during deformations must also start and end at x0.) Definition 2.6 (Path). A path in a space Xis a continuous map γ:I= [0,1] ⊂R→X. Two paths will be considered the same if they are homotopic as continuous maps between topological spaces. As show in Proposition 1.2 in [16], the relation of homotopy on paths with fixed endpoints in any space is an equivalence relation. The equivalence class of a path γunder the equivalence relation of homotopy will be denoted [γ] and called the homotopy class of γ. Given two paths γ1, γ2:I→Xsuch that γ1(1) = γ2(0), there is a composition or product path γ1·γ2that traverses first γ1and then γ2, defined by the formula (γ1·γ2)(s) = γ1(2s),0≤s≤1/2 γ2(2s−1) ,1/2≤s≤1(2.1) This product operation respects homotopy classes (see [16], for example). 5 In particular, we restrict attention to paths γ:I→Xwith the same starting and ending point γ(0) = γ(1) = x0∈X. Such paths are called loops, and the common starting and ending point x0is referred to as the basepoint. The set of all homotopy classes [γ] of loops γ:I→Xat the basepoint x0is denoted π1(X, x0). As is proved in Proposition 1.3 in [16], π1(X, x0) is a group with respect to the product [γ1]·[γ2] = [γ1·γ2]. This group is called the fundamental group of Xat the basepoint x0. π1(X, x0) is the first in a sequence of groups πn(X, x0), called homotopy groups, which are defined in an entirely analogous fashion using the n dimensional cube Inin place of I. Calculation of fundamental group and, in general, higher homotopy groups is a very difficult task even in simple spaces. For example, calculation of higher homotopy groups of spheres involves a non-negligent number of sophisticated algebraic tools. This difficulty made mathematicians opt for other algebraic invariants with simpler computation. Homology groups Homology groups are algebraic invariants with a purely geometric origin. They are defined to capture the geometric idea of hole. A hole is easy to imagine but has a non negligent definition. It can be though as a void in the space, but this idea requires the definition of an ambient space. In order to associate homology groups to a topological space, some decomposition has to be done. This decomposition is known as combinatorialization of the space. We first introduce the required algebraic machinery. Definition 2.7 (Chain complex). Let Rbe a commutative ring with unit. A chain complex Cis a sequence · · · dn+2 −−−→ Cn+1 dn+1 −−−→ Cn dn −→ Cn−1 dn−1 −−−→ · · · where Cnare free R-modules, and dnare R-module homomorphisms, which are called differentials, satisfying the identity dn◦dn+1 = 0, for all n. The subgroup of cycles is defined as Zn(C;R) := ker(dn) and Bn(C;R) := Im(dn+1) is the subgroup of boundaries. The property dn◦dn+1 = 0 implies that Zn(C;R)⊇Bn(C;R), for all n, making the following definition possible. Definition 2.8 (nth homology group). For a chain complex C, the nth homology groups is defined by Hn(C;R) := Zn(C;R)/Bn(C;R) 6 . Note that any object in an homology group represent a class of cycles which are not boundaries so it represents a “hole” in some dimension. Definition 2.9 (Chain map). Let C= (C•,d•)and C0= (C0 •,d0 •)be two chain complexes. A collection of R-module homomorphisms f={f•:C•→C0 •, is called a chain map, written f:C → C0, if for all n: fn−1◦dn= d0 n◦fn(2.2) Proposition 2.10 (3.30 in [21]). Let Cand C0be arbitrary chain complexes. Then any chain map f:C → C0 induces homomorphisms on the homology groups f∗:Hn(C;R)→Hn(C0;R), for all n. In the proof of previous proposition, f∗([c]) := [f(c)] An important part of the theory of chain complexes is the observation that the chain maps can be deformed algebraically. Definition 2.11. Let C= (C•,d•)and C0= (C•,d0 •)be two chain complexes, and let f={f•}and g={g•}be two chain maps f, g :C → C0. A sequence of homomorphisms h={h•}, where hn:Cn→Cn+1 is a chain homotopy between fand gif for all nwe have hn−1◦dn+ dn+1 ◦hn=fn−gn(2.3) Clearly, the existence of a chain homotopy between two maps is an equivalence relation: just replace hwith its negative to show the symmetry, and add two chain homotopy maps to show the transitivity. Without a doubt, the most important use of chain homotopies is to show that two chain maps induce the same homology homomorphisms. Proposition 2.12 (3.32 in [21]). Let C1= (C1 ∗, d1 ∗)and C2= (C2 ∗, d2 ∗)be two chain complexes, and let f, g :C1→ C2be two chain maps such that there exists a chain homotopy hbetween fand g. Then the induced maps f∗, g∗:H•(C1)→H•(C2)are equal. 7 Once the algebraic construction of homology groups has been introduced, it is mandatory to assign a chain complex to every topological space. In order to do that, the space is decomposed into smaller pices called cells in general. Depending on the building blocks used to make that partition we can speak of different types of combinatorialization of topological spaces. The following subsection is devoted to review the most important ways of decomposing a space into smaller pieces. Combinatorialization of topological spaces Each of the decompositions shown below will play a crucial role in the calculation of homology groups of a topological space. For each cellular structure associated to the topological space there is one homology theory. We review here cubical homology, simplicial homology, homology of CW complexes (CW homology) and, finally, singular homology. The later does not become from a combinatorialization of the topological space and, hence, does not require any additional structure. However, is much more difficult to be calculated. Simplicial complexes Definition 2.13 (Simplex). Given p+ 1 points {v0, . . . , vp} ⊂ Rnsuch that the vectors {−−→ v0v1,...,−−→ v0vp} are linearly independent, the p-simplex hv0, . . . , vpiis the convex hull of the set {v0, . . . , vp}, i.e.: hv0, . . . , vpi=(p X i=0 λivi p X i=1 λi= 1 ∧λi≥0) The faces of the simplex hv0, . . . , vpiare the simplices hvi0, . . . , viqiwhere {i0, . . . , iq}⊆{0, . . . , p}and 0≤q≤p. Definition 2.14 (Simplicial complex). A simplicial complex Kin Rnis a set of simplices such that: 1. The faces of all the simplices in Kare also in K. 2. The intersection of any two simplices in Kis empty or a simplex in K. The set of p-simplices in Kis denoted by Kp. 8 A simplicial complex is a collection of sets in Rn, therefore the set |K|=Sσ∈Kσ⊂Rncan be endowed with the relative topology and is, of course, a topological space. Hence the simplicial complex Kcan be viewed as simplicial structure of |K|. A space Xis said to be triangulable if there is a simplicial complex Ksuch that there exists a homemorphism f:|K| → X. Hence, a simplicial structure in Xis given by the set {f(σ)|σ∈K}. For any triangulable space X, the corresponding simplicial structure K enables the calculation of the simplicial homology groups of Xas follows. Definition 2.15 (Simplicial chain complex). Let Kbe a simplicial complex and Ra commutative ring with unit. The R-module of p-chains, Cp(K;R)is the free R-module spanned by Kp. The differential map dpis defined by linear extension of the differential of a simplex, given by the following formula: dp(hv0, . . . , vpi) = p X i=0 (−1)ihv0,...,ˆvi, . . . , vpi(2.4) where the hat over the vertex vimeans the elimination of the vertex vi. The chain complex C(K) = (C•(K;R),d) is the chain complex associated to the simplicial complex K. Is straightforward to prove that the chain complex defined below is, in fact, a chain complex, i.e, the composition of any two consecutive differentials is 0. The p-th homology group Hp(K) can be, then, defined as the p-th homology group of the chain complex C(K). Recall that we have introduced a way of calculating the homology of a simplicial structure Kon a topological space X. For considering the homology groups of X as topological invariant, it must be assured that they do not depend on the triangulation of the space. In the following paragraphs singular homology is presented. This homology theory provides the same results as simplicial homology when boths can be applied. Moreover, singular homology does not depend on any combinatorial structure in the space to be calculated. On the contrary it is not suitable to be used for real calculations. 9 Singular complexes The standard n-simplex, for a non negative integer n, is the set ∆n=((λ0, . . . , λn)∈Rn+1  n+1 X i=0 λi= 1 ∧λi≥0) It is clear that the standard n-simplex is the simplex spanned by the n+ 1 points in Rn+1 (1,0,...,0),...,(0,0,...,1). There are some interesting maps from the standard (n−1)-simplex to the standard n-simplex, called the i-th degeneracy operator and defined by the formula above: i n: ∆n−1→∆n, i n(x0, . . . , xn−1)=(x0, . . . , xi−1,0, xi+1, . . . , xn−1) (2.5) Definition 2.16 (Singular simplex). Let Xbe a topological space. A singular p-simplex is a continuous map σ: ∆p→X. The faces of σare given by the singular p-simplices σ◦i p: ∆p−1→Xfor 0≤i≤p. Once defined the singular simplices, we can define the singular chain complex similarly as the simplicial chain complex. Definition 2.17 (Singular chain complex). Let Xbe a topological space and Ra commutative ring with unit. The R- module of singular p-chains is the free R-module spanned by the singular p-simplices in X. It is denoted as Cp(X;R). Again, the differential map dp is defined as the linear extension of the differential of a singular p-simplex, given in the formula above: dp(σ) = p X i=0 (−1)iσ◦i p(2.6) It is, again, straightforward to check that dp◦dp+1 = 0 and the singular chain complex C(X)=(C∗(X;R), d∗) is, in fact, a chain complex. The singular homology groups of Xare then defined as the homology groups of the singular chain complex C(X). Note that the singular chain complex does not depend on any choice and any two homotopy equivalents spaces have isomorphic homology groups (see [16, Corollary 2.11]). Moreover, if Kis a simplicial structure of X, then both 10 simplicial and singular homology groups are isomorphic (see [16, Theorem 2.27]). Recall that, even though simplicial homology does not provide a practical tool for calculating homology groups of spaces, when a space has some combinatorial structure, e.g. it is a simplicial complex, both calculation agrees as the resulting homology groups are isomorphic. In this sense, singular homology is the reference for all the other homology theories on topological spaces. Cubical complexes We mainly follow the process introduced in [19] with some minor changes. Specifically, Kaczy´nski et al. uses closed cubes as main combinatorial objects while we use open cubes. Acubical cell σis a finite product of intervals: σ=I1× · · · × Ik⊂Rk where Ijis an interval (mj, mj+1) with integer extremes and length one or the singleton {mj}, denoted as (mj), for each j∈ {1, . . . , k}. The interval Ijis referred to as the jth component of σand denoted by Ij(σ). The set of all cubical cells in Rkis denoted by Kk. The set of all cubical cells is K= ∞ [ k=1 Kk(2.7) We usually require the cubical cells to be bounded. Hence, we define Kk ∗,n ={σ∈ Kk: 0 ≤inf Ip(σ),sup Ip(σ)< n, 1≤p≤k}(2.8) where kand nare nonnegative integer numbers. Given a cubical cell σin Rk, its embedding number kis denoted by emb(σ). The dimension of σis defined as the number of unitary intervals in its expression as product of intervals and is denoted by dim(σ). The set of all elementary cubes with dimension pis denoted by Kp. The set of all elementary cubes in Rkwith dimension pis denoted by Kk p. Whenever the dimension of a cubical cell require to be made explicit, it is denoted as a superscript between parenthesis. Therefore, dim σ(p)=p. We also explicitly indicate the dimension of a cell σ(p)writing that σis a p-cell. 11 The closure1of a cubical cell can be decomposed into the union of lowerdimensional cubical cells. If δand σare two cubical cells in Rkof any dimension and δ⊂σ, then δis a face of σand is denoted as δ≤σ. If δis a face of σand δ6=σ, then δis a proper face of σ, denoted as δ < σ. If δis a face of σand dim(δ) = dim(σ)−1 then δis a primary face of σ, denoted by δ∈∂(σ). Therefore, the set of primary faces of a cell σwill be denoted as ∂(σ) As an example, let σbe the open square (0,1) ×(3,4), a 2-dimensional cubical cell in R2. Its closure is the square σ= [0,1] ×[3,4] and is decomposed as follows [0,1] ×[3,4] =(0,1) ×(3,4)∪ (0,1) ×(3) ∪(0,1) ×(4)∪ (0) ×(3,4) ∪(1) ×(3,4)∪ (0) ×(3) ∪(0) ×(4)∪ (1) ×(3) ∪(1) ×(4) Namely, the interior of the square, its four edges and its four vertices. The set of primary faces is given by ∂((0,1) ×(3,4)) = {(0,1) ×(3),(0,1) ×(4),(0) ×(3,4),(1) ×(3,4)} Acubical complex is a collection Kof cubical cells with the same embedding number and such that, for every cubical cell σ∈K, all of its primary faces are in the complex. We denote the set of p-cells in Kas Kp. In figure 2.1 a cubical complex in R2is showed. As usual, using cubical cells as algebraic objects allows us to define the corresponding chain complex whose homology will be called cubical homology of the given cubical complex. Again, as expected, cubical homology agrees with singular homology (see [19] for more details). However, the definition of the differential of a cubical cell is not as simple as the simplicial or singular cases. We will proceed in a similar manner as in the singular case, following the work in [31, 38]. Let σ=I1×· · ·×Inbe a p-cubical cell. Denote Iq= (a0 q, a1 q) with ai qtwo integers such that 0 ≤a1 q−a0 q≤1. Let k1, . . . , kpdenote the indexes of nondegenerated intervals. For any set J⊂ {1, . . . , p}define k(J) = {ki|i∈J} 1Given a set A⊂Rk, its closure is A={p∈Rk:∀ > 0∃q∈A|d(p, q)> }. Informally, the closure of a set is the set itself together with the points that are “infinitely” near to the set. 12 012345678910 0 1 2 3 4 5 6 7 8 9 10 Figure 2.1: Example of cubical complex K0in R2. and for i∈ {0,1}the cubical cell λi J(σ) = I0 1× · · · , I0 q, where I0 j={ai j} if j∈k(J) or I0 j=Ijotherwise. Roughly speaking, for any set of non degenerated intervals of σgiven by indices in J,λi Jrepresents the face given by degenerating those given intervals to the ”lower” value (given for i= 0) or the ”upper” one (given by i= 1). We denote λi {j}as λi j. With the above notation in mind, the differential of a p-cubical cell σis given by: dp(σ) = p X j=1 (−1)j(λ1 j(σ)−λ0 j(σ)) (2.9) It is a simple exercise to check that, in fact, dp◦dp+1 = 0 and, hence, the chain complex given by cubical chains and the differential defined by linear extension of the differential of a cubical cell can be used to define the homology groups of a cubical complex. If we do proceed in a similar manner as we done in the simplicial case, and we call cubical spaces to those homeomorphic to the support of any cubical complex (recall that the support of a cubical complex Kis |K|=Sσ∈Kσ⊂Remb(K)). Again as expected, cubical homology of Kand singular homology of |K|are isomorphic (see [19]). 13 CW Complexes As shown before, computation of homology groups of a space Xrequires some partition of the space into smaller pieces (simplices or cubical cells). CW complexes are a generalization of this procedure that moves away from the simple geometry of cubes or simplices to more abstract shapes. An m-cell is any topological space homeomorphic to an m-dimensional closed unit ball Bm. An open m-cell is a topological space homeomorphic to the interior of a ball Bm. To build a CW complex from its constituent cells, we glue some cells with other cells. The following definition concretes this gluing procedure. Definition 2.18 (Attaching spaces along maps). Let Xand Ybe topological spaces, let A⊆Xbe a closed subspace, and let f:A→Ybe a continuous map. Y∪fXdenotes the quotient space X`Y/ ∼, where the equivalence relation ∼is generated by a∼f(a), for all a∈A. We say that the space Y∪fXis obtained from Yby attaching Xalong f. Observe that the space Y∪fXis equipped with the quotient topology, which means that S⊆Y∪fXis open if and only if q−1(S) is open in X`Y, where q:X`Y→Y∪fXis the quotient map. Definition 2.19 (Constructive definition of CW complex). A CW complex X is obtained by the following inductive construction of the skeletons: (1) The 0-skeleton X(0) is a discrete set. (2) Construct the n-skeleton X(n)by the simultaneous attachment of the n-cells to X(n−1) along their boundaries. In particular, X(n)gets the quotient topology as described above. (3) Equip the space X=S∞ n=0 X(n)with the weak topology: A⊂Xis open if and only if A∩X(n)is open for any n. The definition above asserts that a set A⊆Xis open if and only if f−1 α(A) is open for any cell α, where fα:Bn α→Xis the attachment map (also called the characteristic map). Definition 2.20 (Subcomplex). Let Xbe a CW complex; A⊆Xis called a subcomplex if Ais a union of open cells such that if e⊆A, then the closure e⊆A. 14 (c) f◦h= 0 (d) h◦g= 0 (e) h◦h= 0 This is a classical notion in homological algebra and algebraic topology; see e.g. [10, §12] and comments on the terminology and applications in [14, p. 86]. Note that because of the condition (a), if there exists a chain contraction from Cto C0then their homology modules are isomorphic. Definition 2.33 (AM model). An algebraic minimal model (introduced in [14]), or an AM model for short, of a cell complex Kis a chain contraction from C(K)to a chain complex M= (M∗,dM)such that each Mpis a free R-module and all the non-zero elements in the SNF of each (dM)pare non-invertible in R. An AM model exists for every cell complex, and any two AM models for the same complex are isomorphic. Note that we are not distinguishing here the cellular structure of K, as it can be any of simplicial, cubical or CW complex structure. Cohomology Cohomology is the result of a process of dualization in homology. However, only the chain complex is dualized, instead of homology groups. This dualization operation enriches the homological information of a topological space not only with cohomology groups but with a larger structure called cohomology ring. In this section we review the construction of cohomology groups of any cell complex and how to define some ”products” of chains to induce a ”product” operation in cohomology suitable to define the mentioned ring structure. As usual, given a R-module C, Hom(C,R) denotes the homomorphisms from Cto R. Definition 2.34 (Dual chain complex). Let C= (C•,d•)be a chain complex. Its dual complex C∗= (C•, δ•)is defined as follows: 1. The module of cochains is Cp= Hom(Cp,R) 21 2. The codifferential δ•:C•→C•+1 is defined by δp(ϕ(p)) = ϕ(p)◦dp+1 where ϕ(p)is a p-cochain. Basically, a p-cochain is an homomorphism that assigns an element in Rto every p-cell. As R-modules, Cpis isomorphic to Cp, hence the real difference between a chain complex and a cochain complex is in the connection morphism. In a chain complex, the differential map decreases the dimension while the codifferential increases the dimension. From a matricial point of view, the matrix of δpis the transpose of the matrix of dp+1. This easily proves that δp+1 ◦δp= 0. If we define the submodule of cocycles as Zp= ker δpand the submodule of coboundaries as Bp= Im δp−1, the p-th cohomology group of the cochain complex is defined as Hp=Zp/Bp. Again, this definition does not depend on the choice of the cellularization of the given topological space. If p-th homology generators can be understood as p-dimensional holes (subspaces homotopy equivalent to Sp), p-th cohomology generators can be viewed as ”paths” to cut the complex along them in order to remove a hole. Most of the time, homoloy and cohomology groups are not enough to distinguish two topological spaces. In some of these cases, establishing some relations between (co)homology generators open the possibility of detecting differences between the spaces. The following section is devoted to one of this operations: the cup product. Cup product Sometimes, (co)homology groups are not enough to distinguish two topological spaces. For example, the torus T2=S1×S1and the wedge S2∨S1∨S1have the same (co)homology groups but are not homotopy equivalents (both spaces have non isomorphic fundamental groups, so they cannot be homotopy equivalents as fundamental group is homotopy invariant). Therefore, (co)homology groups need to be enriched with more elaborated algebraic structure. Is in this place where cup product arises. Cup product is an operation on cohomology classes. Given two cohomology classes [ϕ(p)] and [ψ(q)], its cup product is a (p+q) cohomology class. Of course, cup product is closely associated with the relationship between the (co)homology of the Cartesian product of two spaces and the (co)homology of the factors. 22 Cup product on singular or simplicial cohomology is defined as follows: (ϕ(p)`ψ(q))(σ) = ϕ(p)(σp 0)·ψ(q)(σp+q p) (2.13) where σ=hv0, . . . , vp+qiand σj i=hvi, vi+1, . . . , vji. This simple definition is not suitable to be directly translated to cubical or cellular cohomology. Instead, it is required the intervention of the chain map diag induced by the diagonal map diag : x∈X7→ (x, x)∈X×X. Concretely, the cup product is given by the following composition `:C∗(X)× C∗(X)→ C∗(X×X)→ C∗(X) (2.14) where the first map is the K¨uneth map and the second is the map diag∗. Recall that cup product is defined in cohomology as in homology the map diag goes in the opposite direction of the one required. In the previous example, if H1(T2) = Z[a]⊕Z[b] and H1(S2∨S1∨S1) = Z[a0]⊕Z[b0], then a`b6=a0`b0, hence both spaces cannot be homotopy equivalent. In order to compute the cup product, it is helpful the calculation of the approximation of the diagonal map. This approximation map is computed in many different ways depending on the cellularization of the space. The worst case is the CW decomposition, where the diagonal has to be constructed inductively in the dimension, as states the Cellular Approximation Theorem. However, we will not require the application of that result. The approximation of the diagonal map presents a very simple aspect in the case of simplicial homology. In this case, its definition is given below. diagSimp(hv0, . . . , vni) = n X p=0 hv0, . . . , vpi⊗hvp, . . . , vni(2.15) However, the cubical case is not as simple as the simplicial one. For example see [20] where an important proportion of the paper is devoted to the construction of a cubical approximation to the diagonal map. We adopt the ideas in [31, 38], giving as a result the following definition for diagCub diagCub σ(n)=X J⊂{1,...,n} ρJ,J0λ0 J0(σ)⊗λ1 J(σ) (2.16) where ρJ,J0= (−1)ν,ν= # {(i, j)∈J×J0|j < i}and J0={1, . . . , n} \ J. 23 CHAPTER 3 Parallelization strategy for homological calculation based on membrane computing Most applications of topological data analysis (TDA for short) are characterized by using a very large data set from which the internal structure is revealed through TDA. This structure is extracted in the form of connected components, holes and voids of higher dimensions. This characterization is formally described by the notion of homology. Calculation of homology requires a combinatorial representation of space as a cellular complex (usually simplicial or cubical). These applications require calculating real-time homology over a growing data set. This fact makes fundamental the development of parallel algorithms for calculating homology. One parallelization strategy based on membrane computing is shown below. Natural Computing studies new computational paradigms inspired by Nature. It abstracts the way in which Nature computes, giving rise to new models of computation. There are several research fields currently well established in Natural Computing. Among them, Genetic Algorithms introduced by Holland [18] that are inspired by evolution and natural selection in order to find an optimal solution among a large set of feasible solutions; Neural Networks, introduced by McCulloch and Pitts [28], that are based on the interconnection of neurons in brain; or DNA-based molecular computing, initiated by Head [17]. 25 Membrane Computing is a research area in Natural Computing based on the interpretation of the processes that take place inside the cells of living beings as computations. Computational devices in Membrane Computing are called P systems. In general, a P system consists of a membrane structure within which there are multisets of objects, which evolve according to given rules, in an asynchronous non-deterministic maximally parallel manner. From the foundational article [32], different models of P systems have been studied. Based on their architecture, these models can be classified into two broad categories: cell-like P systems and tissue-like P systems. In cell-like P systems, membranes are organized hierarchically into a tree-like structure. The inspiration of this architecture is the set of organs inside the cell. Each and every one of them carries out its biological processes in parallel. Tissue-like P systems were introduced in [26, 27]. These types of P systems are inspired by intercellular communication and cooperation between neurons. The mathematical model of these devices is a network of processors that manipulate symbols and communicate them through preestablished channels. Communication between cells is based on symport / antiport rules. Symport rules move objects across a membrane in a single direction, while antiport rules move objects in opposite directions. In tissue-like P systems the membrane structure is a general undirected graph. The edges of this graph are not given explicitly, but are deduced from the set of rules. Since the initial definition of the tissue-like P systems, several lines of research have been developed and other variants have appeared (see, for example, [1, 2, 13, 23, 30]). Catalyst P systems were introduced by P˘aun in [32]. The fundamental characteristic of these P systems is the presence of objects in the membranes that are not consumed by the application of the rules, but their presence is the membrane is necessary for the rule to be chosen as applicable. Catalysts have been studied extensively in Membrane Computing (see, e.g., [12, 22, 24]). This is not the first relation of Membrane Computing and Computational Algebraic Topology. For instance, in [34] Reina-Molina et al present a Membrane Computing algorithm to solve the problem of thinning 2D and 3D images using cell complexes. In [4, 33], tissue-like P systems are used to segment an image, in [5, 6, 7, 8, 9] Membrane Computing is used for calculating homology groups of 2D digital images. We present here an implementation of a Membrane Computing 26 framework to adapt the ideas behind Discrete Morse Theory in order to develop algorithms in order to solve the problem of computing homology groups of 2D and 3D digital binary images. This work is based in [35]. 3.1 Membrane computing In this section we present various definitions concerned with Membrane Computing. We will use a variant of tissue-like P systems where the application of the rules are regulated by promoters and inhibitors. These promoters have a clear biological inspiration. The rule is applied if the reactants are present, but it is also necessary the presence of all the promoters and none of the inhibitors in the corresponding cell. The promoters are not consumed nor produced by the application of the rule, but if they are not in the cell, the rule cannot be applied. In one step, each reactant in a membrane can be used only for one rule, but if several rules need the presence of the same promoter, then the presence of one unique copy of the promoter suffices for the application of all the rules. In the general case, if there are several possibilities, the rule is non-de- terministically chosen, but sometimes a priority relation between rules is considered, so the concept of priority in our P systems is required. Next, the formal definition of these P systems is recalled. Definition 3.1 (Tissue-like P system with catalysts and priorities). Atissue-like P system with promoters, inhibitors, priorities and input of degree q≥1is a tuple of the form Π = (Γ,Σ,E, w1, . . . , wq,R, Pri, iin, iout) where qis the number of cells (or membranes) of the P system and 1. Γis a finite alphabet, whose symbols are called objects. These objects can be placed in the cells or in the surrounding space (called the environment). 2. Σ⊆Γis the input alphabet. The input of the computation performed by the P system is encoded by using this alphabet. 3. E ⊆ Γis a finite alphabet representing the set of the objects in the environment. Following a biological inspiration, the objects in the environment are available in an arbitrary large amount of copies; 27 4. w1, . . . , wqare strings over Γrepresenting the multisets of objects placed inside the cells at the beginning of the computation; 5. Ris a finite set of rules of the following form: (pro ¬inh |i, u/v, j),for 0≤i6=j≤q, pro, inh, u, v ∈Γ∗ 6. Pri is a finite set of relations Ri> Rj, where Riand Rjare rules from R. It means that if Riand Rjcan be applied, then the application of Rihas priority on the application of Rj. 7. iin ∈ {1,2, . . . , q}denotes the input cell, i.e., the cell where the input of the computation will be placed. 8. iout ∈ {1,2, . . . , q}denotes the output cell, i.e., the cell where the output of the computation will be placed. Informally, a tissue-like P system with promoters, inhibitors and priorities of degree q≥1 can be seen as a set of qcells labeled by 1,2, . . . , q. The cells are the nodes of a virtual graph, where the edges connecting the cells are determined by the communication rules of the P system, i.e., as usual in tissue-like P systems, the edges linking cells are not provided explicitly: if a rule (pro ¬inh |i, u/v, j) is given, then cells iand jare considered linked. The application of a rule (pro ¬inh |i, u/v, j) consists of trading the multiset u(initially in the cell i) against the multiset v(initially in j). After the application of the rule, the multiset udisappears from the cell iand it appears in the cell j. Analogously, the multiset vdisappears from the cell jand it appears in the cell i. The trade can also be between one cell and the environment, labeled by 0. The rule is applied if in the cell with label ithe objects of pro are present in the cell i(promoters), while any of the objects in inh do not appear in the cell (inhibitors). The promoters or the inhibitors are not modified by the application of the rule. If the promoters and inhibitors are empty, we write (i, u/v, j) instead of (∅ ¬∅| i, u/v, j). Finally, we write (pro |i, u/v, j) or (¬inh |i, u/v, j) when only promoters or inhibitors appear, respectively. As usual, some objects not belonging to Ecan arrive to the environment during a computation. So, in a configuration (not initial) it could be found two types of objects in the environment: firstly, those which belong to the environment and appear in an arbitrary large number of copies. Secondly, those which do not belong to the environment but have been sent to the environment by the application of a rule. 28 Rules are used as usual in the framework of membrane computing, that is, in a maximally parallel way (a universal clock is considered). Aconfiguration is an instantaneous description of the P system and it is represented as a tuple (w0, w1, . . . , wq), where w0is the multiset of objects from Γ − E placed in the environment (initially, w0=∅). Given a configuration, we can perform a computation step and obtain a new configuration by applying the rules in a parallel manner as it is shown above. A configuration is halting when no rules can be applied to it. A computation is a sequence of computation steps such that either it is infinite or it is finite and the last step yields a halting configuration (i.e., no rules can be applied to it). Then, a computation halts when the P system reaches a halting configuration. The output of a computation is collected from its halting configuration by reading the objects contained in the output cell. Definition 3.2 (Distributed scheme). AdP scheme (of degree s≥1) is a construct ∆ = (Γ,Σ,Π1,...,Πs, R, iin, iout) where: 1. Γis an alphabet of (communicated) objects. 2. Σ⊂Γis the input alphabet, used as input for the module Πiin . 3. Π1,...,Πsare tissue like P systems with promoters, inhibitors and priorities whose alphabets contain Γ. Each Πiis called a component or module of ∆. 4. Ris a set of distribution rules of the form (pro¬inh|i, u →v, j)where pro,inh,uand vare multisets on Γ,0≤i6=j≤sand i6= 0. The objects in pro are called promoters of the rule and the objects in inh are called inhibitors, as usual in tissue like P systems with promoters and inhibitors. 5. iin and iout represents the modules used as input and output, respectively. Informally, the evolution of a dP scheme can be summarized as follows. A universal common clock is used for all the modules in the scheme. Initially all the modules start evolving together, each of them with their respective initial configuration except the input module, which also uses the input alphabet Σ as input. On each computation step, every time a 29 module reaches its halting configuration, a maximal multiset of distribution rules is selected. Then, these rules are applied and some objects are sent from the output membrane of some modules to the input membrane of another. Then, the evolution starts again until all the modules have reached their halting configuration and no distribution rule can be selected. We use an special destination module, indexed by 0. Every time a rule (pro¬inh|i, u →v, 0) is executed, the multiset uis removed from the output membrane of the ith module. This is used as some kind of garbage collection. 3.2 Encoding images as cubical complexes and cubes as tuples In this section, we work with a cubical cell version of a digital image. A common problem appearing in the development and implementation of algorithms with cubical complexes is the requirement of a huge amount of memory resources. Below, an efficient encoding of cubical cells is described. Given two non-negative integers kand n≥2, the following auxiliary sets are defined Tk n:= {0,1,2, . . . , n −1}k,(3.1) Tk n:= Tk 2n−1,(3.2) In,k := Tnk={0,1,2, . . . , nk−1},(3.3) In,k := T(2n−1)k={0,1,2,...,(2n−1)k−1}.(3.4) The set of points1for the source images is the set Tk n⊂Zkequipped with a cubic neighborhood function, described as follows: two points i= (i1, . . . , ik) and j= (j1, . . . , jk) are said (2k)-adjacent if Pk p=1 |ip−jp|= 1. The neighborhood function is given by N(i1, . . . , ik) = ((j1, . . . , jk)∈Tk n: k X p=1 |ip−jp|= 1) This function restricted to k= 2 defines 4-adjacency and to k= 3 defines 6-adjacency. 1The reader is supposed to be familiar with concepts of Image Algebra. For a detailed text, see [36]. 30 •The input alphabet is given by Σ ={σi:i∈In,k(T(K))}∪{δ± ij :i, j ∈In,k(T(K)) ∧ hσ(i),dσ(j)i=±1}∪ ∪ {V± ij , Vi, Vj:i, j ∈In,k(T(K)) ∧V(σ(i)) = ±σ(j)}∪ ∪ {Ci:i∈In,k(T(K)) ∧V(σ(i)) = 0 ∧V−1(σ(i)) = ∅}∪ ∪ {Eij :i, j ∈In,k(T(K)) ∧ hσ(j),dV(σ(i))i 6= 0 ∧i6=j} •The environment alphabet is given by E= Γ0. •The set of rules Ris given by –R1≡{σiσjσlδs1 il δs2 jl V−s1 il } ¬{φiφj}|1, λ/φ−s1s2 ij φiφj,0for i, j, l ∈ In,k,i6=j6=l6=iand s1, s2∈ {+,−} –R2≡{σiCi} ¬{φ+ ii }|1, λ/φ+ ii ,0for i∈In,k –R3≡{φs1 i1i2Vs2 i2j1φs3 j1j2}|1, Ei2j1/φ−s1s2s3 ij ,0for i, j, l ∈In,k, s1, s2, s3∈ {+,−} and i6=j6=l6=i. –R4≡{δs1 li CiCj} ¬{˜ δs1s2 ji }|1, φs2 lj /˜ δs1s2 ji ˜ Ci˜ Cj˜σi˜σj,0for i, j, l ∈ In,k with i6=j6=l6=i. –R5≡1,˜ δ+ ij ˜ δ− ij /λ, 0for i, j ∈In,k,i6=j. •Pri =∅ •iin =iout = 1. Theorem 3.9. Let K⊂ Kk ∗,n be a cubical complex and Va gradient vector field in K. The module ΠFlow computes the Morse complex for Kassociated to Vin, at most, 4 + blog2νccomputation steps, where ν≤nkis the length of the greatest V-path on K. Remark. The initial configuration for ΠFlow consists in all the elements that codify the cubical complex K(objects σifor cubical cells and Cimarking critical cells), its boundary map (objects δ± ij ) and the gradient vector field (objects V± ij ). There also are objects Eij depicting that σiand σjare distinct maximal faces in V(σi). In this situation, only rules in R1and R2can be selected, as the other rules require the existence of objects φ± ij introduced in the first membrane by these rules. The application of rules R1creates one object φ± ij for each pair of cubical cells σi, σj∈Ksuch that hϕ(σi), σji=±1. Recall that ϕ:C → C given by ϕ= id −d◦Vis the reduced (discrete time) flow of the gradient vector field ([11]). 37 The application of rules R2previously selected introduce objects φ+ ii for each critical cubical cell σi, representing that, for each critical cell µis ϕ(µ) = µ. In the current situation, the second configuration, only rules in R3,R4or R5are available to be selected. Recall that objects φ± ij encodes the presence of a V-path from σito σj with multiplicity equals to ±1. Rules in R3concatenate two V-paths γ and γ0if there is a cell σsuch that the ending cell of γ0, named σ(i), the starting cell of γlay on its boundary and V(σ(i)) = σ. Hence, rules in R3are applied until a maximal path is constructed. This process needs, at most, log2νsteps, where νis the length of the longest maximal V-path. Rules in R4get some objects from the environment, encoding the cells of the Morse complex. Indeed, for each critical path from σito σjwith multiplicity ±1, there is an object φ± ij, hence for each critical cubical cell σlwith hσi,dσli 6= 0 there will be as many objects φ± ij as the number of critical V-paths from σito σjwith multiplicity ±1, respectively. Finally, rules in R5cancel pairs of objects ˜ δ+ ij ,˜ δ− ij . The halting configuration of ΠFlow has as many objects ˜ δ± ij as the incidence of σiin the boundary of σj. Notice that the entire process needs, at most, 4 + blog2νcsteps, where νis the length of the longest critical V-path. We deal here with a membrane algorithm for computing an optimal gradient vector field for a 3D digital object I. Working with a cubical “cellularization” K(I) of I, we show that the critical cells of dimension p (p= 0,1,2,3) of this particular gradient vector field coincide with the pth Betti number of K(I). Discrete Morse Theory establishes a tool to simplify a cell complex4 while the homology groups are kept. This can be used to calculate the homology groups of appropriate cell complexes. Concretely, the technique explained below require the complex to be torsion free. An acyclic vector field over a cubical complex Kis built, removing as many cubical cells as possible. In case that the resulting Morse complex has non-null differential, another acyclic vector field is constructed to accomplish the so called critical cell cancellation (see section 11 in [11]). In this way, a sequence of Morse reductions is built until a complex with null-differential is “reached”. To ensure the finiteness of this process, the 38 (a) Example of cubical complex. (b) Acyclic vector field. (c) Critical cells where two kinds of critical cells are shown: ordinary critical cells representing a hole at any dimension and other critical cells showing a “step” configuration which will be canceled later (d) Critical complex. (e) Final homological complex with null differential. Figure 3.1: Example of homology group calculation in a cubical complex. 39 cubical complex has to be torsion free. Figure 3.1 presents an example of the calculation of homology groups for a 2D cubical complex. In this example the strategy presented above is shown. First, an acyclic vector field is built such that almost every cubical p-cell σis paired to one of its co-faces µ(cubical (p+ 1)-cell with σ∈dµ). The algorithm used to determine the acyclic vector field returns two kind of critical cubical cells: on one hand, cubical p-cells that represent a pdimensional hole and, on the other hand, 1-cells that represents an “step” configuration (as hP11, P12ior hP15, P16iin figure 3.1). Therefore, another acyclic vector field is required to cancel these spurious critical cells. Notice that the first Morse complex is not a cubical complex but a CW complex. However, it is usually much more simpler than the original cubical complex. For binary 3D images, it suffices to build up to two acyclic vector fields to calculate an homotopically equivalent cell complex with null differential, which makes homology group calculation trivial. In the paragraphs below, two P modules are defined to perform this task: ΠVector and ΠVector2. Definition 3.10 (P system to create an acyclic vector field). Let I:T3 n→ {0,1}be a k-D binary image. The module ΠVector is defined as ΠVector ={Γ,Σ,E, w1,R, Pri, iin, iout} where •Γ=Γ0 •Σ = {σi:i∈In,k(T(K(I)))} • E = Γ0 •w1=∅ • R is the following set of rules: –R1≡¬{ViVj}|1, σiσjδ± ij dipU+ ij,dCiCj/σiσjδ± ij U+ ij,dV∓ ij ViVj,0for i, j ∈In,3,i6=j,0≤p≤3and 1≤d≤3 These rules take a critical cubical p-cell (σi), a critical cubical (p+1)-cell (σj) with σi∈dσjand such that the barycentric vector ~vσiσjis 1 2~ud. Then, objects V± ij are brought from the environment (this represents that V(σi) = ±σj) as long as objects Viand Vj (used to prevent that a cubical cell is used in two different vectors in the vector field). 4CW complex in [11] and cubical complex in this paper. 40 •The priorities are given by Pri ={R1[i, j, p, d]> R1[i, j0, p, d0] : i, j, j0∈In,3∧0≤p≤3∧1≤d<d0≤3} ∪ {R1[i, j, p, d]> R1[i0, j, p, d0] : i, i0, j ∈In,3∧0≤p≤3∧1≤d < d0≤3} ∪ {R1[i, j, p, d]> R1[j, j0, p + 1, d0] : i, j, j0∈In,3∧0≤p≤2∧1≤d≤3} where R1[i, j, p, d]represents the corresponding instance of rule R1for the given values of the parameters. The priorities in the first set are used to assign to a cubical cell σia cubical cell σjwith the corresponding barycentric vector parallel to the vector ~udwith the lowest das possible. The priorities of the second kind are used to choose between two cubical cells σiand σi0when another cubical cell σjcan be assigned to be the image by V. The criteria is again to select the cubical cell with the lowest barycentric vector. The last priority type is used to choose whether a cubical cell is chosen as the source cell of a vector or the end cell of another. In this case it is used as the end cell. •iin =iout = 1 Informally speaking, the vector field tries to “go” in the direction of ~u1. When it is not possible to go in that direction, it is used ~u2and, if it cannot be performed, ~u3is used. Hence, the flow of Vcan be thought as a waterfall going downside until it goes from left to right and, finally, rear to back. Definition 3.11 (P system to reduce a vector field). Let I:T3 n→ {0,1}be a k-D binary image. The module ΠVector2is defined as ΠVector2={Γ,Σ,E, w1,R, Pri, iin, iout} where •Γ=Γ0 •Σ = {σi:i∈In,k(T(M))}where Mis the Morse complex built from the acyclic vector field created by ΠVector. • E = Γ0 •w1=∅ 41 • R is the following set of rules: –R1≡¬{ViVj}|1, σiσjδ± ij CiCj/σiσjδ± ij V∓ ij ViVjV, 0for i, j ∈In,3, i6=j. These rules take two cells that can be used to build a vector and create that vector. Notice that, once a cell is used in a vector, it will not be used in another. •Pri =∅. •iin =iout = 1 This P module cancels pairs of critical cells of ΠVector due to the previously mentioned “step” configuration. The dP scheme ∆Bin3DHom = (Γ0,Σ,ΠCub,ΠVector,ΠFlow,ΠVector2,R, iin, iout) with four modules defined below, computes the homology groups for the cubical complex generated from I. •Σ = {Bi:i∈ I−1(1)} •The modules are given by ΠCub This module builds a cubical complex from the 3D binary image I. ΠVector This module creates a maximal acyclic vector field. ΠFlow This module determines the Morse complex associated to the vector field created above. ΠVector2This module creates a maximal acyclic vector field in the Morse complex. •The distribution rules are given by – R1≡(Cub, σiCi→σiCi,Vector) for i∈In,3 – R2≡Cub, δ± ij →δ± ij ,Vectorfor i, j ∈In,3 – R3≡Cub, U+ ij,d →U+ ij,d,Vectorfor i, j ∈In,3and 1 ≤d≤3 – R4≡(Cub,dip →dip,Vector) for i∈In,3and 0 ≤p≤3 – R5≡(Vector, σi→σi,Flow) for i∈In,3 42 – R6≡(Vector,Ci→Ci,Flow) for i∈In,3 – R7≡Vector, δ± ij →δ± ij ,Flowfor i, j ∈In,3, i 6=j – R8≡(Vector,dip →dip,Flow) for i∈In,3and 0 ≤p≤3 – R9≡Vector, V ± ij →V± ij ,Flowfor i, j ∈In,3, i 6=j – R10 ≡Flow,˜σi˜ Ci→σiCi,Vector2for i∈In,3 – R11 ≡Flow,˜ δ± ij →δ± ij ,Vector2for i, j ∈In,3, i 6=j – R12 ≡({V}|Vector2, σi→σi,Flow) for i∈In,3 – R13 ≡({V}|Vector2,Ci→Ci,Flow) for i∈In,3 – R14 ≡{V}|Vector2, δ± ij →δ± ij ,Flowfor i, j ∈In,3, i 6=j – R15 ≡{V}|Vector2, V ± ij →V± ij ,Flowfor i, j ∈In,3, i 6=j •iin = Cub, iout = Vector2 The modules ΠV ector and ΠVector2are the only modules that do not rely on the framework. In other words, the user of the framework only has to deal, most of the times, with writing modules for calculating acyclic vector fields, as the tasks for building the cubical complex associated to the binary image and the building of the Morse complex associated to it, are left to the framework. Theorem 3.12. The dP scheme ∆Bin3DHom defined above calculates the homology groups for a given binary 3D image I:T3 n→ {0,1}in O(log2n)computation steps. Remark. The result is proved by following the behavior of the dP scheme ∆Bin3DHom for an input image I. First of all, the input alphabet for ΠCub is set up from Ias ΠCub is the input module. No distribution rules can be selected until one module reaches its halting stage. By theorem 3.6, ΠCub halts in two steps. Then, the only distribution rules that can be selected are R1,R2,R3and R4. After it application, the module ΠVector starts calculating a vector field. Rules R1in ΠVector build a vector from a cell σito a co-face σj. The priorities ensures that all the p-cells are mapped by V, the vector field, to ap+ 1-cell, prioritizing that the vector joining its barycenter has as lowest index as possible. Also, the rules ensure that the vector field is acyclic, as all the vectors are “parallel” to the vectors in the canonical base in R3 with positive direction. This work is done in only one computation step. 43 This way of building a vector field creates two kind of critical cells. First of all, there are critical cells due to homology generators. Second, there are critical cells due to “steps” in the cubical complex. In figure 3.1a, the cells hP11, P16iand hP11, P12iare in such position. It looks somewhat an “step” in a stair in a bottom view. This configuration of critical cell is due to the priority relation among rules in ΠVector. Hence, we have two kinds of critical cells, those representing an homology generator and those cells that must have to be canceled in later computation steps. Once ΠVector halts, the only distribution rules that can be selected are R5to R9, sending objects to ΠFlow. Concretely, those rules send cells (σi), boundary map (δ± ij ), critical cells (Ci) and vector field V± ij . Module ΠFlow calculates the Morse complex associated to the previous vector field in, at most, 4+3 log2(n−1). To prove this note that the longest V-path in a 3D cubical complex with, at most, (n−1)31-cells have length, at most, (n−1)3and, then, use theorem 3.9. Once the Morse complex is built, it is sent to the module ΠVector2by distribution rules R10 and R11. This module builds an acyclic vector field that “removes” all the cells that do not represent an homology generator. This computation is completed in only one step. If some vector field is built, which means that there are cells that can be cancelled, at least one object Vcomes from the environment that ensures the selection and application of distribution rules R12 to R15, sending the Morse complex and the gradient vector field to the module ΠFlow, where another Morse complex is built and sent to ΠVector2. Once the second Morse complex has came to ΠVector2the dP scheme ∆Bin3DHom reaches the halting state, as no other distribution rule can be selected. This takes place as the only critical cells that are not homology generators are 1-cells that are in a “step” configuration, and those cells are removed by the second vector field. Formally, we have reduced the initial cubical complex to another complex with null differential, which means that all the cells are generators of the homology groups. Recall that the sum of the computation steps required by each module is, in the worst case, at most 12 + 6log2(n−1), which is O(log2n). The full calculation process of homology groups can be followed in Figure 3.1. The module ΠCub builds the cubical complex in Figure 3.1a in two computation steps. Then, the module ΠV ector builds the acyclic vector field in Figure 3.1b in one computation step. Following, the module ΠFlow constructs the (first) Morse complex in 7 steps (because the longest V-path 44 in the complex has length 6). Figure 3.1c shows the two kinds of critical cells. On one hand, P19 is a critical 0-cell due to an homology generator. On the other hand, P12 and P16 are 0-cells due to “step” configuration. hP5, P6i is a critical 1-cell due to “hole” while hP11, P12iand hP15, P16iare critical 1- cells originated by a “step” configuration. Next, the module ΠVector2specify, in one computation step, another acyclic vector field in the complex in Figure 3.1d which leads to the computation in 5 computation steps of the homological complex in Figure 3.1e, that has null differential, so that all the cells are homology generators. 3.4 Conclusions of the chapter In this chapter we have presented a brief introduction to Membrane Computing and how it can be used to efficiently compute the homology generators of 3D digital objects. The presented algorithms are executed sequentially but each of them is massively parallel in the way that all the rules that can be applied in parallel. However, for the best of our knowledge, there is no computing device capable of execute “membrane code” as if it was a living cell and some adaptation have to be performed in those algorithms. On the other hand, Membrane Computing algorithms do not count with algebraic operations in its operations, but precisely Algebra is one of our most useful tool as we are working with Computational Algebraic Topology. Hence we must restrict our research in the development of Membrane computing algorithms in Computational Algebraic Topology to those situations where the algebra can be eluded. This is the case of homology groups of 3D digital objects, which are torsion free. In the following chapters the ideas presented above are used to develop some massively parallel algorithms that provides an efficient tool to extract (co)homological information of digital objects. 45 The algorithm 2 details the construction of an acyclic vector field in the reduced complex. Algorithm 2 Acyclic vector field for a generic CW complex. 1: function AcyclicCWVectorField(K) 2: V←(0 : R[K•]→ R[K•+1]) 3: G←DirectedGraph() 4: for σ∈Kdo 5: G.add node(σ) 6: if σis critical then 7: for τ∈δ(σ)do 8: if V(τ) = 0 ∧τ /∈Im(V)∧ hσ, d(τ)i ∈ R∗then 9: G.add node(τ) 10: G.add edge(σ, τ) 11: for µ∈∂(τ)\ {σ}do 12: if not G.has node (µ)then 13: G.add node(µ) 14: end if 15: G.add edge(τ, µ) 16: end for 17: if G.has cycles() then 18: G.remove node(τ) 19: else 20: V(σ)← hσ, d(τ)iτ 21: break 22: end if 23: end if 24: end for 25: end if 26: end for 27: return V 28: end function Recall that a cell σis critical if it is not either the source or the target of any vector in V. In other words,a critical cell σverifies that V(σ) = 0 and σis not in the support of any chain c∈Im(V). This can be easily checked in the matrices of V. Let us denote by [Vp] th matrix of Vp:Cp→Cp+1 in the standard basis of chains. Then σ(p)is critical if both the p-th column in [Vp] and the p-th row in [Vp−1] are zero. 53 Theorem 4.1. The function Vreturned by algorithm 2 is an acyclic vector field. Proof. To prove that the map Vreturned by algorithm 2 is an acyclic vector field it is mandatory to prove (a) # supp(V(σ)) ≤1 where (c) denotes the support of a chain c, defined by supp Pi∈Iciσi={σi|i∈I}; (b) V2= 0 and (c) there is no closed non-trivial V-paths. (a) This statement is equivalent to prove that any (p+ 1)-cell τis the target of, at most, one vector in V. This fact is guaranteed by line 8 in algorithm 2. (b) To prove that V2= 0 is enough to show that for any cell σ. First of all, suppose that V(σ) = hσ, d(τ)iτ. Let us prove that V(τ) = 0. In fact, line 6 in algorithm 2 only consider assigns a value to V(τ) if τ /∈Im(V). Hence, if τ∈supp(Im(V)) then V(τ) = 0, so V2(σ) = 0. (c) The acyclicity of Vis a direct consequence of the acyclicity of the graph Gin algorithm 2. The SNF reduction Pilarczyk implemented in [31] a SNF algorithm such that the change of bases matrices of two consecutive dimension are mutually inverses, more concretely: [dp] = Ap−1·Dp·Gp(4.22) where [dp] is the matrix of the differential map at dimension p,Dpis in SNF and Ap=G−1 p. This algorithm allows us to construct a reduction of a chain complex to another in which the differential matrices are in SNF, where the projection is given by matrices Apand inclusion is given by matrices Gp. Finally, the composition of all the reductions computed above gives a reduction from the original cubical complex to a reduced complex in Smith Normal Form. 54 4.3 Conclusions In this chapter we have seen a theoretical development to calculate an AM- model for a cell complex. This model is based on iterative calculation of Morse reductions until the differential of the final complex does not have invertible inputs in the ground ring. A Morse reduction is given univocally by an acyclic vector field. It is in this step where massive parallelism is applied. Once the cubical cell complex is drastically reduced, the resulting chain complex is reduced again in order to remove all the pairs σ(p), τ(p+1) such that hσ, d(τ)iis invertible in the ground ring. This cancellations are selected sequentially and performed in parallel. Recall that matricial operations involved in composition of reductions are implemented using the massive parallelism provided by current GPU. In chapter 6 some implementation details are to be discussed. 55 CHAPTER 5 Advanced (co)homological information extraction from digital objects 5.1 Introduction An AM-model not only provides information on the (co)homology groups of a cell complex, but also supplies a relationship between the corresponding (co)chain complexes. This relation, given by the projection and inclusion chain maps, allows an analysis of several characteristics in (co)homology. For example, the inclusion of a reduced complex generator is always a (co)cycle. In addition, the homotopy operator of an AM-model can be interpreted as a ”geometric” homotopy that allows to reduce the original complex into the reduced complex. In this chapter we will how to obtain a large amount of (co)homology information from an AM-model. We also show how to compute cohomology operations (cup and cap products) in the reduced complex so we can get more information in order to distinguish two cell complexes. From here to the end of the chapter, we assume that Kis a cell complex, C={Cp,d}is its associated chain complex, and (f, g, h) : C ⇒⇒ M ={Mp,dM} 57 is an AM-model. For simplifying the presentation, we will assume that the ground ring is Z. 5.2 (Co)homology information from AM-model The first information of relevance at homology level is, of course, the homology and cohomology groups. The AM-model directly provides the (co)homology groups using the two following results. Proposition 5.1. A generator σin Mpis a free homology generator if and only if dM(σ) = 0 and d∗ M(σ)=0. Moreover, if dM(σ)=0and d∗ M(σ) = kτ then σis a homology generator of the torsion subgroup with coefficient k. From now on and depending on the context, a chain σcan also means the corresponding cochain whose value is 1 only in σ. Proof. First of all, if dM(σ) = 0 it is clear that σ∈ker(dMp) for some p. If d∗ M(σ) = 0 then σis not in Im(dMp+1). So σis an homology generator of the free part of the p-th homology group. On the other hand, if dM(σ) = 0 and d∗ M(σ) = kτ, then σcontributes with a factor Zto ker(dMp). Also, as d∗ M(σ) = kτ, there is a factor kZin Im(dMp+1). Therefore, there is a factor Z/kZ=Zkin the quotient that defines homology. By dualization of the proposition below, we can state the following result concerning cohomology groups. Proposition 5.2. A generator σin Mpis a free cohomology generator if and only if d∗ M(σ) = 0 and dM(σ)=0. Moreover, if d∗ M(σ)=0and dM(σ) = kτ then σis a cohomology generator of the torsion subgroup with coefficient k. Let us bear in mind that Mp∼ =M∗ pas any chain can be uniquely viewed as a cochain taking the value 1 on itself. By direct combination of previous results, we can state the following theorem. Theorem 5.3. Let (f, g, h):(Cp,d) ⇒⇒ (Mp,dM)be an AM-model. Then: 58 (a) Hp= M σ∈Mp∧dM(σ)=0∧d∗ M(σ)=0 Z ⊕ M σ∈Mp∧dM(σ)=0∧d∗ M(σ)=kτ Zk  (b) Hp= M σ∈Mp∧d∗ M(σ)=0∧dM(σ)=0 Z ⊕ M σ∈Mp∧d∗ M(σ)=0∧dM(σ)=kτ Zk  Note that the second formula (b) is another way of expressing a consequence of the Universal Coefficient Theorem on Cohomology in the finite dimensional case. In other words, Hp∼ =Zβp⊕Torp−1where βpis the rank of Hp, the p-th Betti number, and Torp−1is the (p−1)-th torsion subgroup. Another homological property is the (co)boundary relationship. Let us consider two generators σ∈Mpand τ∈Mp+1. We say that σis a face of τ in Mif there is a face µof τin Ksuch that f(µ) = σ. We can see here how the boundary relationship in Kis propagated to M. Forman proved in [11, Corollary 3.5] that, given a Morse function1f,Kis homotopy equivalent to a CW complex with exactly one cell per critical cell. This CW complex is, concretely, the complex M. The boundary relationship defined above is another step in the determination of the minimal CW structure for a given cell complex. The invariant factor decomposition for finitely generated Abelian groups states that any finitely generated Abelian group Gis isomorphic to the direct sum Zn⊕Zk1⊕ · · · ⊕ Zkrwhere ki|ki+1. This result, combined with the last theorem, allows us to proof the following theorem. Theorem 5.4. The chain complex Mis minimal in the number of cells. Proof. The construction algorithm for an AM-model (see section 4.2) ensures that the matrices of dMare in SNF and has no ±1 entries. Each non invertible coefficient represents a torsion coefficient. Of course, the matrix representation of the invariant factor decomposition is diagonal diag(1,...,1, k1, . . . , kr,0,...,0) where the 1 is repeated ntimes. Now, as the differentials dMhas no invertible entry, the matrices of dMhas the form diag(k1, . . . , kr,0,...,0). And those kiare precisely the torsion coefficients calculated as in Proposition 5.1. Hence, the chain complex Mis built bearing in mind that any two cells with incidence invertible in the ground 1It is equivalent to give a Morse function or to give an acyclic vector field (see [11]). 59 ring are cancelled. So the remaining cells are those essential for homology. The using of SNF ensures that the differentials are either zero or a non invertible multiple of a cell. 5.3 Cohomology operations from AM-model Cup product (see Section 2.2) defines a ring structure in the graded cohomology group. This operation allows to topologically discriminate, for example, a torus from a wedge of an sphere with two rings. In this case, (co)homology groups cannot distinguish these two spaces. Cup product is easily defined in simplicial homology.In cubical homology it is required a cubical approximation of the diagonal map. We mainly follows the approach in [31]. The definition of cup product on the cohomology ring H∗(K) of a digital object2Kis easy starting from an AM-model (f, g, h) and a cubical approximation of the diagonal map diag :C(K)→ C(K)⊕ C(K) (see equation 2.16). Concretely, given two cochains ϕ(p)and ψ(q), (ϕ(p)`ψ(q))(σ(p+q)) = (µ◦(ϕ⊗ψ)◦∆) (σ) (5.1) where ∆ is the approximation of the diagonal map on the reduced complex Mgiven by ∆ = (f⊗f)◦diag ◦g(5.2) Cap product is another operation involving cohomology classes. Its motivation and definition can be found in, for example, [16, Chapter 3]. It can be thought as a map from C(K)⊗ C∗to C(K) such that assigns a (p−q)-chain to the tensor product of a p-chain and a q-cochain. The following equation express this relation in terms of the K¨unneth formula _:C ⊗ C∗∆⊗id −−−→ C ⊗ C ⊗ C∗id ⊗ −−→ C (5.3) where :C ⊗ C∗→ R is the evaluation map defined below σ(p)⊗ϕ(q)=ϕ(σ) if p=q 0 otherwise (5.4) and the first id is the identity map on cochains while the second one is the identity map on chains. 2Note that digital objects and cubical complexes can be exchanged without problems. 60 If we define the projections onto the factors as πi:C1⊗ C2→ Cifor i= 1,2 with πi(σ1⊗σ2) = σi, then the formula for cap product is given by σ _ ϕ =ϕ(π2(∆(σ))) π1(∆(σ)) (5.5) Note how the equation above can be easily computed in the reduced complex as the aproximation of the diagonal map has been previously calculated. 5.4 Conclusions to the chapter In this chapter, we show the versatility of an AM-model representation in order to get (co)homological information of digital objects. This information extraction begins with the construction of the cubical complex associated to the digital object. This construction is computed in parallel using a quasi literal implementation of the P system ΠCub in definition 3.5. This adaptation basically check all the available cubical cells within the cubes that cover all the points in the digital object and the cubical cells in the complex are those such that all their vertices are in the object. Once the cubical complex associated to the digital object has been calculated, we need to create an AM-model and this is easily performed following the steps in section 4.2. Recall that the construction of an acyclic vector field in the cubical complex is performed in parallel. The design patterns used in the parallelization of the construction of the vector field in the cubical complex guarantee that the vector field is acyclic, however some critical cells that, in fact, do not represent any homology class of the object (called false critical cells), remain. This spurious cells are removed in the next steps. The spanning tree strategy used in the construction of an acyclic vector field can be parallelized using the framework in [40]. We have not implemented this parallelization in the software presented in chapter6 and it is delayed to a future work. The iterative application of the algorithm 2 results in a cell complex, namely a CW complex, whose differential map is either zero or has all the coefficients in the units set of the ground ring. This complex is transformed in another with the same cells but with an equivalent differential in Smith Normal Form. This last step allows to easily calculate all the (co)homology groups of the digital object as the are the same of the last reduced complex. Note that the condition on the differential ensures the minimality of the 61 AM-model. As for each unit coefficient in the differential a reduction can be performed so a pair of cells would disappear. The AM-model brings in two chain maps, projection and inclusion, that allow us to map cells from the reduced complex to (co)cycles in the original complex and cells in the original complex to homology generators in the reduced complex. The approximation of the diagonal given in 5.2 is crucial in the computation of cup and cap products. Both represents operations in (co)homology that leads to homotopy invariants that provides more interesting information of a digital object and relates holes ((co)homology generator) at different dimensions. For instance, the cup product of the two 1-dimensional holes of the torus is precisely its 2-dimensional homology generator. 62 and the differential d1is given by d1 0= 0 (6.5) d1 1= 0 (6.6) d1 2(h0,1,5i) = −h0,8i−h3,5i(6.7) d1 2(h2,5,6i) = h5,9i−h6,9i+h0,8i(6.8) d2 2(h0,4,7i) = h5,9i+h6,9i−h7,9i+h3,5i(6.9) d1 2(h1,7,8i) = h7,9i(6.10) As there are invertible coefficients in the matrices of d1, another acyclic vector field V2can be calculated on M1. The reduced complex is given in the following equations. M2 0=Z[h9i]; d2 0= 0 (6.11) M2 1=Z[h6,9i]⊕Z[h0,8i]; d2 1= 0 (6.12) M2 2=Z[h2,5,6i]; d2 2(h2,5,6i) = −2h6,9i+ 2h0,8i(6.13) Here, we can see that this reduced complex cannot be reduced using a vector field, so let us compute the SNF reduction.This gives as result the final reduced complex Mshown below M0=Z[h9i]; dM0= 0 (6.14) M1=Z[h6,9i]⊕Z[h0,8i]; dM1= 0 (6.15) M2=Z[h2,5,6i]; dM2(h2,5,6i) = 2h6,9i(6.16) This final complex is modified in order to have all the coefficients non negative. The composition of all the reductions calculated above results in the reduction ρAM = (f, g, h) : C ⇒⇒ M defining the AM-model for the Klein bottle. In figure 6.3 the vector field generated by ρAM is shown. Please note that the reduction generated by this vector field is not ρAM as this reduction has been calculated using the SNF reduction. From ρAM several homological information can be extracted. First, the homology groups of the Klein bottle are given by H0(K) = Z[h9i] (6.17) H1(K) = Z[h0,8i]⊕Z2[h6,9i] (6.18) 69 Figure 6.1: Simplicial structure on the Klein bottle. The barycenters of each cell are also drawn. Figure 6.2: Acyclic vector field on the Klein bottle. Critical cells have its barycenters drawn in green. 70 Figure 6.3: Final acyclic vector field on the Klein bottle from the AM- model. Critical cells have its barycenters drawn in green. Note that both homology generators in dimension 1 are not really cycles. Properly speaking, the homology generators are given below g(h0,8i) = −h0,1i−h1,2i+h0,8i−h2,8i(6.19) g(h6,9i) = h8,9i−h0,6i−h6,9i+h0,8i(6.20) However, for the sake of simplicity and bearing in mind that gtakes homology generators in Mto cycles generating homology classes in C, we use cells5in Mas homology classes. From the AM-model, the cohomology groups can also be easily calculated. This gives as result H0=Z[h9i∗] (6.21) H1=Z[h0,8i∗] (6.22) H2=Z2[h2,5,6i∗] (6.23) From the contravariance of cohomology functor, the dual reduction ρ∗ AM is given by ρ∗ AM = (g∗, f∗, h∗) (6.24) 5As every chain complex is a a sequence of free modules, we call cells to the elements in the corresponding basis. 71 Please note that the dual of the inclusion is the projection and viceversa. Hence, cohomology generators in the original complex are given below f∗(h9i∗) = h5i∗+h0i∗+h6i∗+h1i∗+h7i∗+h2i∗+ +h8i∗+h3i∗+h9i∗+h4i∗(6.25) f∗(h0,8i∗) = h6,9i∗− h2,6i∗+h6,8i∗+h4,8i∗− h5,6i∗+h0,3i∗+ +h0,8i∗− h3,4i∗+h0,2i∗− h3,5i∗(6.26) f∗(h2,5,6i∗) = h2,5,6i∗(6.27) Cup product can also be easily calculated using the aproximation of the diagonal given by equation 5.2. On simplicial Klein bottle, the cup product in integer cohomology is given by h9i∗`h9i∗=h9i∗(6.28) h9i∗`h0,8i∗=h0,8i∗(6.29) h9i∗`h2,5,6i∗=h2,5,6i∗(6.30) h0,8i∗`h9i∗=h0,8i∗(6.31) h0,8i∗`h0,8i∗= 0 (6.32) h0,8i∗`h2,5,6i∗= 0 (6.33) h2,5,6i∗`h9i∗=h2,5,6i∗(6.34) h2,5,6i∗`h0,8i∗= 0 (6.35) h2,5,6i∗`h2,5,6i∗= 0 (6.36) Cubical Klein bottle Now we will repeat the same calculation above using a cubical version of Klein bottle found in the source code associated to [31]. This cubical complex has 42 0-cells, 84 1-cells and 42 2-cells in R6. Homology and cohomology groups are given by: H0=Z[(0) ×(1) ×(0) ×(1) ×(1) ×(1)] (6.37) H1=Z[(1) ×(1) ×(1) ×(0) ×(0) ×(0,1)] ⊕ ⊕Z2[(1) ×(1) ×(0) ×(0) ×(0,1) ×(1)] (6.38) H0=Z[(0) ×(1) ×(0) ×(1) ×(1) ×(1)] (6.39) H1=Z[(1) ×(1) ×(1) ×(0) ×(0) ×(0,1)] (6.40) H2=Z2[(0) ×(0,1) ×(0,1) ×(0) ×(0) ×(0)] (6.41) 72 `ν(0) σ(1) τ(2) ν(0) ν σ τ σ(1) σ0 0 τ(2) τ0 0 Table 6.1: Product table on the cohomology ring of the Klein bottle. as expected by previous calculations. However, the (co)homology generators inside the original cubical complex are not presented here as they are too long and do not contribute with any interesting information. Cup product is given by table 6.1, whose results are congruent with the simplicial case. Note that ν,σand τare the generator of cohomology at dimensions 0, 1 and 2, respectively. Bing’s house with two rooms Bing’s house, also known as house with two rooms, is a contractible6 2-complex that is not collapsible (see [3]). In figure 6.4 a cubical decomposition of Bing’s house is shown. The (co)homology groups are given below H0=Z[(5) ×(5) ×(3)] (6.42) H0=Z[((5) ×(5) ×(3))∗] (6.43) The inclusion of the cohomology generator is a chain supported by 90 0- cells. Both results are consistent with the fact that Bing’s house is null homotopic. Of course, the cup product is trivial. 6A cell complex is contractible if it is homotopy equivalent to a point. 73 Figure 6.4: Cubical complex decomposition of Bing’s house. Faces (cubical 2-cells) are coloured using its normal vector. Torus and two handled sphere Consider the torus T=S1×S1with cubical decomposition given in figure 6.5. The (co)homology groups are: H0=Z[(3) ×(3) ×(1)] (6.44) H1=Z[(2) ×(2,3) ×(0),(1,2) ×(1) ×(1)] (6.45) H2=Z[(0,1) ×(0,1) ×(0)] (6.46) H0=Z[(3) ×(3) ×(1)] (6.47) H1=Z[(2) ×(2,3) ×(0),(1,2) ×(1) ×(1)] (6.48) H2=Z[(0,1) ×(0,1) ×(0)] (6.49) 74 Figure 6.5: Cubical decomposition of torus. Figure 6.6: Homology generators at dimension 1 75 The two 1-cycles generating homology group at dimension 1, are c1 1=−(2) ×(2,3) ×(1) + (2) ×(2,3) ×(0) −(2) ×(2) ×(0,1) + +(2) ×(3) ×(0,1) (6.50) c1 2= (2) ×(2,3) ×(1) + (1,2) ×(1) ×(1) −(1) ×(1,2) ×(1) + +(2) ×(1,2) ×(1) −(1) ×(2,3) ×(1) −(1,2) ×(3) ×(1) (6.51) The non trivial cohomology cup product is the one given by `: H1×H1→H2where the non zero products are given by (2) ×(2,3) ×(0) `(1,2) ×(1) ×(1) = (0,1) ×(0,1) ×(0) (6.52) Two handled sphere is the space S1∨S2∨S1. Its homology is trivially isomorphic to the one of the torus. Below are the generators: H0=Z[(0) ×(3) ×(1)] (6.53) H1=Z[(0) ×(0,1) ×(0),(0) ×(2,3) ×(0)] (6.54) H2=Z[(0,1) ×(1,2) ×(0)] (6.55) H0=Z[(0) ×(3) ×(1)] (6.56) H1=Z[(0) ×(0,1) ×(0),(0) ×(2,3) ×(0)] (6.57) H2=Z[(0,1) ×(1,2) ×(0)] (6.58) Therefore, torus and two handled sphere are homologically indistinguishable. However, as previously shown, both spaces are not homotopy equivalent. The cup product is, in two handled sphere, trivial, as can be computed used the equation 5.1. 6.3 Empirical complexity analysis We approach the complexity analysis of the proposed algorithms from an empirical point of view. We choose this approach as the algorithm complexity finally depends on the number of spurious critical cells, that are those critical cells in one step that are removed in the next stage. We perform this analysis as follows. First, a battery of random cubical complexes are created. Then, an AM-model is constructed for each complex. Finally, the point cloud formed by the variables size of the complexes, given by its number of cells, and time used in its computation in the development computer, is approximated using the best fitting curve. This curve represents, precisely, the empirical complexity of the proposed 76 0 20 40 60 80 100 0.0 0.5 1.0 1.5 2.0 2.5 3.0 (a) 2D cubical complexes 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (b) 3D cubical complexes 0 20 40 60 80 100 0 5 10 15 20 25 30 (c) 4D cubical complexes Figure 6.7: Empirical complexity graph for 2, 3 and 4 dimensional complexes. Blue dots represents pairs of (size, time) data. Blue and red lines represents the best fitting polynomial with degree 2 and 3 respectively. algorithm. The curve election has been performed by experimentation with different types of curves, resulting the better results with polynomials of degree two and three. In figure 6.7 it can be observed how the empirical complexity is quadratic (blue line), with practically no difference with the cubical case (red line). 77 6.4 Conclusion Along this chapter a software prototype is presented. This prototype is developed using Python as it is easy to read a write, is multiparadigm and has a lot of implemented modules that help to the developer in many different tasks. However, it is an interpreted language that experience sometimes less speed than other compiled language. Notwithstanding the foregoing, this implementation can be used as an “assembly map” to build another compiled version much faster. Another improvement that can be done to the implemented prototype is related with the parallelization of the vector field simplification stage. As mentioned above, a parallel implementation of the spanning tree algorithm in [40] leads to an important speed up. It also has to be improved one of the most used matrix operation, which is the product of matrices. In cases where the dimension of the matrices are greater than two thousand rows or columns, the multiplication performs so slow. Nonetheless, this can be improved using CUDA libraries. As can be seen from previous paragraphs, many improvements can be done to the software in order to be used in real applications. However, this improvements are left to future work. However, the starting point determined by the work presented here is promising enough, as it shows a good complexity as can be observed in figure 6.7. As final conclusion to this dissertation, two different approaches to the parallelization of (co)homological information extraction are presented. Both approaches lives in two different universes, however they are closely linked to each other. It can be thought as one being the concretion of other to the real world, as it can be easily implemented in current computers. At this point it is beyond doubt the utility of Computational Algebraic Topology in many practical applications. Hence, the work presented here provides a tool to extract many interesting (co)homological information from digital objects regardless the dimension of the object or the ground ring. In many of the literature consulted, the tools and research are relegated to be only in 3D or the ground ring is supposed to be a field. In both cases the homology groups are torsion free. In this dissertation we focus, precisely, in ensuring that the torsion coefficients are calculated. This presents the main obstacle to a full Membrane Computing implementation of the information extraction procedure, as those coefficients have a purely algebraic nature, not combinatorial, which is better to be boarded from Membrane Computing. 78