Bug fixed.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/periodic_cubical_complex_fix@2821 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: f5816f9d7ca51a30c5f64dc7cd58f40ffc157917

5 files changed, 523 insertions, 73 deletions

diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
index 67735ba1..6c02eff1 100644
--- a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex.cpp
@@ -54,7 +54,7 @@ int main(int argc, char** argv) {
 
   // Compute the persistence diagram of the complex
   Persistent_cohomology pcoh(b);
-  int p = 2;
+  int p = 11;
   double min_persistence = 0;
 
   pcoh.init_coefficients(p);  // initializes the coefficient field for homology
diff --git a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
index 122160a2..32935205 100644
--- a/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
+++ b/src/Bitmap_cubical_complex/example/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
@@ -57,7 +57,7 @@ int main(int argc, char** argv) {
   // Compute the persistence diagram of the complex
   Persistent_cohomology pcoh(b, true);
 
-  int p = 2;
+  int p = 11;
   double min_persistence = 0;
   pcoh.init_coefficients(p);  // initializes the coefficient field for homology
   pcoh.compute_persistent_cohomology(min_persistence);
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index 0442ac34..58e5bbbd 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -1,4 +1,4 @@
-/*    This file is part of the Gudhi Library. The Gudhi library
+/*    This file is part of the Gudhi Library. The Gudhi 
  *    (Geometric Understanding in Higher Dimensions) is a generic C++
  *    library for computational topology.
  *
@@ -100,6 +100,8 @@ class Bitmap_cubical_complex_base {
    * non-negative integer, indicating a position of a cube in the data structure.
    * In the case of functions that compute (co)boundary, the output is a vector if non-negative integers pointing to
    * the positions of (co)boundary element of the input cell.
+   * The boundary elements are guaranteed to be returned so that the 
+   * incidence coeficcients are alternating.
    */
   virtual inline std::vector< size_t > get_boundary_of_a_cell(size_t cell)const;
   /**
@@ -112,12 +114,66 @@ class Bitmap_cubical_complex_base {
    * In the case of functions that compute (co)boundary, the output is a vector if
    * non-negative integers pointing to the
    * positions of (co)boundary element of the input cell.
+   * Note that unlike in the case of boundary, over here the elements are 
+   * not guaranteed to be returned with alternating incidence numbers.
+   * 
    **/
   virtual inline std::vector< size_t > get_coboundary_of_a_cell(size_t cell)const;
   /**
    * In the case of get_dimension_of_a_cell function, the output is a non-negative integer
    * indicating the dimension of a cell.
+   * Note that unlike in the case of boundary, over here the elements are 
+   * not guaranteed to be returned with alternating incidence numbers.
+   * To compute incidence between cells use compute_incidence_between_cells
+   * procedure
    **/
+   
+   /**
+   * This procedure compute incidence numbers between cells. 
+   * Note that first parameter is the address of a cell of dimension n,
+   * and the second parameter is the address of adjusted cell in dimension
+   * n-1. 
+  **/ 
+  virtual int compute_incidence_between_cells( size_t coBoundary , size_t boundary )
+  {	  	  
+	  
+	  //first get the counters for coBoundary and boundary:
+	  std::vector<unsigned> cbd_counter = this->compute_counter_for_given_cell( coBoundary );
+	  std::vector<unsigned> bd_counter = this->compute_counter_for_given_cell( boundary );
+	  
+	  //cbd_counter and bd_counter should agree at all positions except from one:
+	  int number_of_position_in_which_counters_do_not_agree = -1;
+	  size_t number_of_full_faces_that_comes_before = 0;
+	  for ( size_t i = 0 ; i != cbd_counter.size() ; ++i )
+	  {
+		  if ( (cbd_counter[i]%2 == 1)&&(number_of_position_in_which_counters_do_not_agree==-1) )
+		  {
+			  ++number_of_full_faces_that_comes_before;
+		  }
+		  if ( cbd_counter[i] != bd_counter[i] )
+		  {
+			  if ( number_of_position_in_which_counters_do_not_agree != -1 )
+			  {
+				  std::cout <<  "Cells given to compute_incidence_between_cells procedure do not form a pair of coboundary-boundary.\n";
+				  throw "Cells given to compute_incidence_between_cells procedure do not form a pair of coboundary-boundary.";
+			  }
+			  number_of_position_in_which_counters_do_not_agree = i;
+		  }
+	  }
+	  
+	  int incidence = 1;
+	  if ( number_of_full_faces_that_comes_before%2 )incidence = -1;	 	  
+	  //if the boundary cell is on the right from coboundary cell:
+	  if ( cbd_counter[number_of_position_in_which_counters_do_not_agree]+1 == 
+	       bd_counter[number_of_position_in_which_counters_do_not_agree]
+	     )
+	     {
+			 incidence *= -1;
+		 }		
+	 	 	  
+	  return incidence;
+  }
+   
   inline unsigned get_dimension_of_a_cell(size_t cell)const;
   /**
    * In the case of get_cell_data, the output parameter is a reference to the value of a cube in a given position.
@@ -289,7 +345,7 @@ class Bitmap_cubical_complex_base {
    * boundary_simplex_range creates an object of a Boundary_simplex_range class
    * that provides ranges for the Boundary_simplex_iterator.
    **/
-  Boundary_range boundary_range(size_t sh) {
+  Boundary_range boundary_range(size_t sh) {	  	  
     return this->get_boundary_of_a_cell(sh);
   }
 
@@ -668,21 +724,32 @@ Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const char* perseus_
 }
 
 template <typename T>
-std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const {
+std::vector< size_t > Bitmap_cubical_complex_base<T>::get_boundary_of_a_cell(size_t cell)const {	
   std::vector< size_t > boundary_elements;
 
   // Speed traded of for memory. Check if it is better in practice.
   boundary_elements.reserve(this->dimension()*2);
 
-  size_t cell1 = cell;
+  size_t sum_of_dimensions = 0;
+  size_t cell1 = cell;  
   for (size_t i = this->multipliers.size(); i != 0; --i) {
     unsigned position = cell1 / this->multipliers[i - 1];
     if (position % 2 == 1) {
-      boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
-      boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+	  if ( sum_of_dimensions%2 )
+	  {
+		  boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+		  boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+	  }
+	  else
+	  {			  
+		  boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+		  boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+	  }      
+      ++sum_of_dimensions;
     }
     cell1 = cell1 % this->multipliers[i - 1];
-  }
+  }    
+  
   return boundary_elements;
 }
 
@@ -690,23 +757,24 @@ template <typename T>
 std::vector< size_t > Bitmap_cubical_complex_base<T>::get_coboundary_of_a_cell(size_t cell)const {
   std::vector<unsigned> counter = this->compute_counter_for_given_cell(cell);
   std::vector< size_t > coboundary_elements;
-  size_t cell1 = cell;
+  size_t cell1 = cell;  
   for (size_t i = this->multipliers.size(); i != 0; --i) {
     unsigned position = cell1 / this->multipliers[i - 1];
     if (position % 2 == 0) {
-      if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
-        coboundary_elements.push_back(cell - this->multipliers[i - 1]);
-      }
-      if (
-          (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
-        coboundary_elements.push_back(cell + this->multipliers[i - 1]);
-      }
+	  if ((cell > this->multipliers[i - 1]) && (counter[i - 1] != 0)) {
+		coboundary_elements.push_back(cell - this->multipliers[i - 1]);
+	  }
+	  if (
+		  (cell + this->multipliers[i - 1] < this->data.size()) && (counter[i - 1] != 2 * this->sizes[i - 1])) {
+		coboundary_elements.push_back(cell + this->multipliers[i - 1]);
+	  }	     
     }
     cell1 = cell1 % this->multipliers[i - 1];
-  }
+  }     
   return coboundary_elements;
 }
 
+
 template <typename T>
 unsigned Bitmap_cubical_complex_base<T>::get_dimension_of_a_cell(size_t cell)const {
   bool dbg = false;
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
index c3cc93dd..9411f4b2 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h
@@ -88,14 +88,73 @@ class Bitmap_cubical_complex_periodic_boundary_conditions_base : public Bitmap_c
   /**
    * A version of a function that return boundary of a given cell for an object of
    * Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   * The boundary elements are guaranteed to be returned so that the 
+   * incidence coeficcients are alternating.
    */
   virtual std::vector< size_t > get_boundary_of_a_cell(size_t cell) const;
 
   /**
    * A version of a function that return coboundary of a given cell for an object of
    * Bitmap_cubical_complex_periodic_boundary_conditions_base class.
+   * Note that unlike in the case of boundary, over here the elements are 
+   * not guaranteed to be returned with alternating incidence numbers.
+   * To compute incidence between cells use compute_incidence_between_cells
+   * procedure
    */
   virtual std::vector< size_t > get_coboundary_of_a_cell(size_t cell) const;
+  
+  
+  /**
+   * This procedure compute incidence numbers between cells. 
+   * Note that first parameter is the address of a cell of dimension n,
+   * and the second parameter is the address of adjusted cell in dimension
+   * n-1. 
+  **/ 
+  virtual int compute_incidence_between_cells( size_t coBoundary , size_t boundary )
+  {	  	  
+	  //first get the counters for coBoundary and boundary:
+	  std::vector<unsigned> cbd_counter = this->compute_counter_for_given_cell( coBoundary );
+	  std::vector<unsigned> bd_counter = this->compute_counter_for_given_cell( boundary );
+	  
+	  //cbd_counter and bd_counter should agree at all positions except from one:
+	  int number_of_position_in_which_counters_do_not_agree = -1;
+	  size_t number_of_full_faces_that_comes_before = 0;
+	  for ( size_t i = 0 ; i != cbd_counter.size() ; ++i )
+	  {
+		  if ( (cbd_counter[i]%2 == 1)&&(number_of_position_in_which_counters_do_not_agree==-1) )
+		  {
+			  ++number_of_full_faces_that_comes_before;
+		  }
+		  if ( cbd_counter[i] != bd_counter[i] )
+		  {
+			  if ( number_of_position_in_which_counters_do_not_agree != -1 )
+			  {
+				  std::cout <<  "Cells given to compute_incidence_between_cells procedure do not form a pair of coboundary-boundary.\n";
+				  throw "Cells given to compute_incidence_between_cells procedure do not form a pair of coboundary-boundary.";
+			  }
+			  number_of_position_in_which_counters_do_not_agree = i;
+		  }
+	  }
+	  
+	  int incidence = 1;
+	  if ( number_of_full_faces_that_comes_before%2 )incidence = -1;	 	  
+	  //if the boundary cell is on the right from coboundary cell:
+	  if ( (cbd_counter[number_of_position_in_which_counters_do_not_agree]+1 == 
+	       bd_counter[number_of_position_in_which_counters_do_not_agree])
+	       ||
+	       (  
+	       (cbd_counter[number_of_position_in_which_counters_do_not_agree] != 1) 
+	       &&
+	       (bd_counter[number_of_position_in_which_counters_do_not_agree]==0)
+	       )
+	     )
+	     {
+			 incidence *= -1;
+		 }		
+	 	 	  
+	  return incidence;
+  }
+  
 
  protected:
   std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed;
@@ -222,45 +281,72 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T
   bool dbg = false;
   if (dbg) {
     std::cerr << "Computations of boundary of a cell : " << cell << std::endl;
-  }
+  }  
 
   std::vector< size_t > boundary_elements;
+  boundary_elements.reserve(this->dimension()*2);  
   size_t cell1 = cell;
+  size_t sum_of_dimensions = 0;
+  
   for (size_t i = this->multipliers.size(); i != 0; --i) {
     unsigned position = cell1 / this->multipliers[i - 1];
     // this cell have a nonzero length in this direction, therefore we can compute its boundary in this direction.
-
     if (position % 2 == 1) {
       // if there are no periodic boundary conditions in this direction, we do not have to do anything.
-      if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {
-        // std::cerr << "A\n";
-        boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
-        boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+      if (!directions_in_which_periodic_b_cond_are_to_be_imposed[i - 1]) {		  
+		//std::cerr << "A\n";
+        if ( sum_of_dimensions%2 )
+        {
+			boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+			boundary_elements.push_back(cell + this->multipliers[ i - 1 ]); 			
+		}
+		else
+		{
+			boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);        
+			boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+	    }
         if (dbg) {
           std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " ";
         }
       } else {
         // in this direction we have to do boundary conditions. Therefore, we need to check if we are not at the end.
         if (position != 2 * this->sizes[ i - 1 ] - 1) {
-          // std::cerr << "B\n";
-          boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
-          boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+		   //std::cerr << "B\n";			
+          if ( sum_of_dimensions%2 )
+          {
+			  boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+			  boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);			  			  
+		  }
+		  else
+		  {
+			  boundary_elements.push_back(cell + this->multipliers[ i - 1 ]);
+			  boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+	      }
           if (dbg) {
             std::cerr << cell - this->multipliers[ i - 1 ] << " " << cell + this->multipliers[ i - 1 ] << " ";
           }
         } else {
-          // std::cerr << "C\n";
-          boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
-          boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]);
+		  //std::cerr << "C\n";
+          if ( sum_of_dimensions%2 )
+          {
+			 boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);
+			 boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]);	          			 
+		  }
+		  else
+		  {
+			boundary_elements.push_back(cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ]);	          
+			boundary_elements.push_back(cell - this->multipliers[ i - 1 ]);			
+		  }
           if (dbg) {
             std::cerr << cell - this->multipliers[ i - 1 ] << " " <<
                 cell - (2 * this->sizes[ i - 1 ] - 1) * this->multipliers[ i - 1 ] << " ";
           }
         }
       }
+      ++sum_of_dimensions;
     }
     cell1 = cell1 % this->multipliers[i - 1];
-  }
+  }         
   return boundary_elements;
 }
 
@@ -295,7 +381,7 @@ std::vector< size_t > Bitmap_cubical_complex_periodic_boundary_conditions_base<T
     }
 
     cell1 = cell1 % this->multipliers[i - 1];
-  }
+  }  
   return coboundary_elements;
 }
 
diff --git a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
index db90eb94..f815a70c 100644
--- a/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
+++ b/src/Bitmap_cubical_complex/test/Bitmap_test.cpp
@@ -42,6 +42,8 @@ Bitmap_cubical_complex_periodic_boundary_conditions_base;
 typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_periodic_boundary_conditions_base>
 Bitmap_cubical_complex_periodic_boundary_conditions;
 
+
+
 BOOST_AUTO_TEST_CASE(check_dimension) {
   std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
 
@@ -100,24 +102,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
   std::vector<double> boundary8;
   boundary8.push_back(1);
   boundary8.push_back(15);
-  boundary8.push_back(7);
   boundary8.push_back(9);
+  boundary8.push_back(7);
   std::vector<double> boundary9;
   boundary9.push_back(2);
   boundary9.push_back(16);
   std::vector<double> boundary10;
   boundary10.push_back(3);
   boundary10.push_back(17);
-  boundary10.push_back(9);
   boundary10.push_back(11);
+  boundary10.push_back(9);
   std::vector<double> boundary11;
   boundary11.push_back(4);
   boundary11.push_back(18);
   std::vector<double> boundary12;
   boundary12.push_back(5);
   boundary12.push_back(19);
-  boundary12.push_back(11);
   boundary12.push_back(13);
+  boundary12.push_back(11);
   std::vector<double> boundary13;
   boundary13.push_back(6);
   boundary13.push_back(20);
@@ -140,24 +142,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
   std::vector<double> boundary22;
   boundary22.push_back(15);
   boundary22.push_back(29);
-  boundary22.push_back(21);
   boundary22.push_back(23);
+  boundary22.push_back(21);
   std::vector<double> boundary23;
   boundary23.push_back(16);
   boundary23.push_back(30);
   std::vector<double> boundary24;
   boundary24.push_back(17);
   boundary24.push_back(31);
-  boundary24.push_back(23);
   boundary24.push_back(25);
+  boundary24.push_back(23);
   std::vector<double> boundary25;
   boundary25.push_back(18);
   boundary25.push_back(32);
   std::vector<double> boundary26;
   boundary26.push_back(19);
   boundary26.push_back(33);
-  boundary26.push_back(25);
   boundary26.push_back(27);
+  boundary26.push_back(25);
   std::vector<double> boundary27;
   boundary27.push_back(20);
   boundary27.push_back(34);
@@ -180,24 +182,24 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
   std::vector<double> boundary36;
   boundary36.push_back(29);
   boundary36.push_back(43);
-  boundary36.push_back(35);
   boundary36.push_back(37);
+  boundary36.push_back(35);
   std::vector<double> boundary37;
   boundary37.push_back(30);
   boundary37.push_back(44);
   std::vector<double> boundary38;
   boundary38.push_back(31);
   boundary38.push_back(45);
-  boundary38.push_back(37);
   boundary38.push_back(39);
+  boundary38.push_back(37);
   std::vector<double> boundary39;
   boundary39.push_back(32);
   boundary39.push_back(46);
   std::vector<double> boundary40;
   boundary40.push_back(33);
   boundary40.push_back(47);
-  boundary40.push_back(39);
   boundary40.push_back(41);
+  boundary40.push_back(39);
   std::vector<double> boundary41;
   boundary41.push_back(34);
   boundary41.push_back(48);
@@ -275,11 +277,13 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_1) {
   for (size_t i = 0; i != increasing.size(); ++i) {
     std::vector< size_t > bd = increasing.get_boundary_of_a_cell(i);
     for (size_t j = 0; j != bd.size(); ++j) {
-      BOOST_CHECK(boundaries[i][j] == bd[j]);
-    }
+      BOOST_CHECK(boundaries[i][j] == bd[j]);            
+    }   
   }
 }
 
+
+
 BOOST_AUTO_TEST_CASE(compute_boundary_test_2) {
   std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
 
@@ -373,16 +377,19 @@ BOOST_AUTO_TEST_CASE(compute_boundary_test_2) {
   coboundaryElements.push_back(40);
   coboundaryElements.push_back(41);
   coboundaryElements.push_back(47);
+   
   size_t number = 0;
   for (size_t i = 0; i != increasing.size(); ++i) {
     std::vector< size_t > bd = increasing.get_coboundary_of_a_cell(i);
     for (size_t j = 0; j != bd.size(); ++j) {
       BOOST_CHECK(coboundaryElements[number] == bd[j]);
-      ++number;
-    }
+      ++number;      
+    }   
   }
 }
 
+
+
 BOOST_AUTO_TEST_CASE(compute_boundary_test_3) {
   std::vector< double > increasingFiltrationOfTopDimensionalCells({1, 2, 3, 4, 5, 6, 7, 8, 9});
 
@@ -578,50 +585,50 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) {
 
   std::vector<double> boundary0;
   std::vector<double> boundary1;
-  boundary1.push_back(0);
   boundary1.push_back(2);
+  boundary1.push_back(0);
   std::vector<double> boundary2;
   std::vector<double> boundary3;
-  boundary3.push_back(2);
   boundary3.push_back(0);
+  boundary3.push_back(2);
   std::vector<double> boundary4;
-  boundary4.push_back(0);
   boundary4.push_back(8);
+  boundary4.push_back(0);
   std::vector<double> boundary5;
-  boundary5.push_back(1);
   boundary5.push_back(9);
+  boundary5.push_back(1);
   boundary5.push_back(4);
   boundary5.push_back(6);
   std::vector<double> boundary6;
-  boundary6.push_back(2);
   boundary6.push_back(10);
+  boundary6.push_back(2);
   std::vector<double> boundary7;
-  boundary7.push_back(3);
   boundary7.push_back(11);
+  boundary7.push_back(3);
   boundary7.push_back(6);
   boundary7.push_back(4);
   std::vector<double> boundary8;
   std::vector<double> boundary9;
-  boundary9.push_back(8);
   boundary9.push_back(10);
+  boundary9.push_back(8);
   std::vector<double> boundary10;
   std::vector<double> boundary11;
-  boundary11.push_back(10);
   boundary11.push_back(8);
+  boundary11.push_back(10);
   std::vector<double> boundary12;
-  boundary12.push_back(8);
   boundary12.push_back(0);
+  boundary12.push_back(8);
   std::vector<double> boundary13;
-  boundary13.push_back(9);
   boundary13.push_back(1);
+  boundary13.push_back(9);
   boundary13.push_back(12);
   boundary13.push_back(14);
   std::vector<double> boundary14;
-  boundary14.push_back(10);
   boundary14.push_back(2);
+  boundary14.push_back(10);
   std::vector<double> boundary15;
-  boundary15.push_back(11);
   boundary15.push_back(3);
+  boundary15.push_back(11);
   boundary15.push_back(14);
   boundary15.push_back(12);
 
@@ -647,7 +654,7 @@ BOOST_AUTO_TEST_CASE(boudary_operator_2d_bitmap_with_periodic_bcond) {
     std::vector< size_t > bd = cmplx.get_boundary_of_a_cell(i);
     for (size_t j = 0; j != bd.size(); ++j) {
       BOOST_CHECK(boundaries[i][j] == bd[j]);
-    }
+    }    
   }
 }
 
@@ -771,6 +778,8 @@ BOOST_AUTO_TEST_CASE(bitmap_2d_with_periodic_bcond_filtration) {
     BOOST_CHECK(filtration[i] == cmplx.get_cell_data(i));
   }
 }
+
+
 BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical_complex_base_check)
 {
     std::vector< double > expected_filtration;
@@ -836,14 +845,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(10);
     expected_boundary.push_back(1);
     expected_boundary.push_back(11);
-    expected_boundary.push_back(5);
     expected_boundary.push_back(7);
+    expected_boundary.push_back(5);
     expected_boundary.push_back(2);
     expected_boundary.push_back(12);
     expected_boundary.push_back(3);
     expected_boundary.push_back(13);
-    expected_boundary.push_back(7);
     expected_boundary.push_back(9);
+    expected_boundary.push_back(7);
     expected_boundary.push_back(4);
     expected_boundary.push_back(14);
     expected_boundary.push_back(10);
@@ -854,14 +863,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(20);
     expected_boundary.push_back(11);
     expected_boundary.push_back(21);
-    expected_boundary.push_back(15);
     expected_boundary.push_back(17);
+    expected_boundary.push_back(15);
     expected_boundary.push_back(12);
     expected_boundary.push_back(22);
     expected_boundary.push_back(13);
     expected_boundary.push_back(23);
-    expected_boundary.push_back(17);
     expected_boundary.push_back(19);
+    expected_boundary.push_back(17);
     expected_boundary.push_back(14);
     expected_boundary.push_back(24);
     expected_boundary.push_back(20);
@@ -936,7 +945,7 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
         Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
         for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
         {
-            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );                                    
             ++bd_it;
         }
 
@@ -1022,14 +1031,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(10);
     expected_boundary.push_back(1);
     expected_boundary.push_back(11);
-    expected_boundary.push_back(5);
     expected_boundary.push_back(7);
+    expected_boundary.push_back(5);
     expected_boundary.push_back(2);
     expected_boundary.push_back(12);
     expected_boundary.push_back(3);
     expected_boundary.push_back(13);
-    expected_boundary.push_back(7);
     expected_boundary.push_back(9);
+    expected_boundary.push_back(7);
     expected_boundary.push_back(4);
     expected_boundary.push_back(14);
     expected_boundary.push_back(10);
@@ -1040,14 +1049,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(20);
     expected_boundary.push_back(11);
     expected_boundary.push_back(21);
-    expected_boundary.push_back(15);
     expected_boundary.push_back(17);
+    expected_boundary.push_back(15);
     expected_boundary.push_back(12);
     expected_boundary.push_back(22);
     expected_boundary.push_back(13);
     expected_boundary.push_back(23);
-    expected_boundary.push_back(17);
     expected_boundary.push_back(19);
+    expected_boundary.push_back(17);
     expected_boundary.push_back(14);
     expected_boundary.push_back(24);
     expected_boundary.push_back(20);
@@ -1124,7 +1133,7 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
         Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
         for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
         {
-            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );            
             ++bd_it;
         }
 
@@ -1208,14 +1217,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(10);
     expected_boundary.push_back(1);
     expected_boundary.push_back(11);
-    expected_boundary.push_back(5);
     expected_boundary.push_back(7);
+    expected_boundary.push_back(5);
     expected_boundary.push_back(2);
     expected_boundary.push_back(12);
     expected_boundary.push_back(3);
     expected_boundary.push_back(13);
-    expected_boundary.push_back(7);
     expected_boundary.push_back(9);
+    expected_boundary.push_back(7);
     expected_boundary.push_back(4);
     expected_boundary.push_back(14);
     expected_boundary.push_back(10);
@@ -1226,14 +1235,14 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
     expected_boundary.push_back(20);
     expected_boundary.push_back(11);
     expected_boundary.push_back(21);
-    expected_boundary.push_back(15);
     expected_boundary.push_back(17);
+    expected_boundary.push_back(15);
     expected_boundary.push_back(12);
     expected_boundary.push_back(22);
     expected_boundary.push_back(13);
     expected_boundary.push_back(23);
-    expected_boundary.push_back(17);
     expected_boundary.push_back(19);
+    expected_boundary.push_back(17);
     expected_boundary.push_back(14);
     expected_boundary.push_back(24);
     expected_boundary.push_back(20);
@@ -1310,7 +1319,7 @@ BOOST_AUTO_TEST_CASE(all_cells_iterator_and_boundary_iterators_in_Bitmap_cubical
         Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
         for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
         {
-            BOOST_CHECK( expected_boundary[bd_it] == *bd );
+            BOOST_CHECK( expected_boundary[bd_it] == *bd );                        
             ++bd_it;
         }
 
@@ -1376,3 +1385,290 @@ BOOST_AUTO_TEST_CASE(Top_dimensional_cells_iterator_range_check)
     }
 }
 
+
+
+
+BOOST_AUTO_TEST_CASE( check_if_boundary_of_boundary_is_zero_non_periodic_case_3_d )
+{	
+    std::vector< unsigned > sizes(3);
+    sizes[0] = 5;
+    sizes[1] = 5;
+    sizes[2] = 5;
+
+    std::vector< double > data(125 , 0);
+    
+    int number_of_all_elements = (2*sizes[0]+1)*(2*sizes[1]+1)*(2*sizes[2]+1);  
+    Bitmap_cubical_complex_base ba( sizes , data );
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );   
+        int i = 1;
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {			
+            Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd);
+            int j = 1;
+            for ( Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += i*j;
+				j *= -1;
+			}
+            i *= -1;
+        }        
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {
+			BOOST_CHECK( elems_in_boundary[i] == 0 );
+		}          
+    }
+}
+
+
+
+
+BOOST_AUTO_TEST_CASE( check_if_boundary_of_boundary_is_zero_non_periodic_case_4_d )
+{	
+    std::vector< unsigned > sizes(4);
+    sizes[0] = 3;
+    sizes[1] = 3;
+    sizes[2] = 3;
+    sizes[3] = 3;
+
+    std::vector< double > data(81		 , 0);
+    
+    int number_of_all_elements = (2*sizes[0]+1)*(2*sizes[1]+1)*(2*sizes[2]+1)*(2*sizes[3]+1);  
+    Bitmap_cubical_complex_base ba( sizes , data );
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {					
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );   
+        int i = 1;
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {			
+            Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd);
+            int j = 1;
+            for ( Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += i*j;
+				j *= -1;
+			}
+            i *= -1;
+        }        
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {
+			BOOST_CHECK( elems_in_boundary[i] == 0 );
+		}          
+    }
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( check_if_boundary_of_boundary_is_zero_periodic_case_2d )
+{	
+    std::vector< unsigned > sizes(2);
+    sizes[0] = 10;
+    sizes[1] = 10;    
+    
+    std::vector<bool> directions_of_periodicity(2,true);
+    std::vector< double > data(100 , 0);
+    
+    int number_of_all_elements = (2*sizes[0])*(2*sizes[1]);// *(2*sizes[2]);  
+    Bitmap_cubical_complex_periodic_boundary_conditions ba( sizes , data , directions_of_periodicity );
+    for ( Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );   
+        int i = 1;       
+        
+        //std::cout << "Element : " << *it << std::endl;         
+        
+        Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+			//std::cout <<  *bd << " ";
+            Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd);
+            int j = 1;
+            for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += i*j;
+				j *= -1;
+			}
+            i *= -1;
+        }   
+        //getchar();        
+             
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {	
+			BOOST_CHECK( elems_in_boundary[i] == 0 );		
+		}          
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE( check_if_boundary_of_boundary_is_zero_periodic_case_3d )
+{	
+    std::vector< unsigned > sizes(3);
+    sizes[0] = 5;
+    sizes[1] = 5;
+    sizes[2] = 5;
+    
+    std::vector<bool> directions_of_periodicity(3,true);
+
+    std::vector< double > data(125 , 0);
+    
+    int number_of_all_elements = (2*sizes[0])*(2*sizes[1])*(2*sizes[2]);  
+    Bitmap_cubical_complex_periodic_boundary_conditions ba( sizes , data , directions_of_periodicity );
+    for ( Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+		//std::cout << "Element : " << *it << std::endl;
+		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );   
+        int i = 1;                
+        
+        Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd);
+            //std::cout << *bd << " ";
+            int j = 1;
+            for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{				
+				elems_in_boundary[ *bd2 ] += i*j;
+				j *= -1;
+			}
+            i *= -1;
+        }                   
+        //getchar();
+             
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {
+			BOOST_CHECK( elems_in_boundary[i] == 0 );		
+		}          
+    }
+}
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( check_if_boundary_of_boundary_is_zero_periodic_case_4d )
+{	
+    std::vector< unsigned > sizes(4);
+    sizes[0] = 4;
+    sizes[1] = 4;
+    sizes[2] = 4;
+    sizes[3] = 4;
+    
+    std::vector<bool> directions_of_periodicity(4,true);
+
+    std::vector< double > data(256 , 0);
+    
+    int number_of_all_elements = (2*sizes[0])*(2*sizes[1])*(2*sizes[2])*(2*sizes[3]);  
+    Bitmap_cubical_complex_periodic_boundary_conditions ba( sizes , data , directions_of_periodicity );
+    for ( Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );   
+        int i = 1;                
+        
+        Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd);
+            int j = 1;
+            for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += i*j;
+				j *= -1;
+			}
+            i *= -1;
+        }           
+             
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {
+			BOOST_CHECK( elems_in_boundary[i] == 0 );		
+		}          
+    }
+}
+
+
+
+
+
+BOOST_AUTO_TEST_CASE( compute_incidence_between_cells_test )
+{	
+    std::vector< unsigned > sizes(3);
+    sizes[0] = 5;
+    sizes[1] = 5;    
+    sizes[2] = 5;
+        
+    std::vector< double > data(125 , 0);
+    
+    int number_of_all_elements = (2*sizes[0]+1)*(2*sizes[1]+1)*(2*sizes[1]+1);
+    Bitmap_cubical_complex_base ba( sizes , data );
+    for ( Bitmap_cubical_complex_base::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );                                         
+        Bitmap_cubical_complex_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_base::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+            Bitmap_cubical_complex_base::Boundary_range second_bdrange = ba.boundary_range(*bd);         
+            for ( Bitmap_cubical_complex_base::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += 
+				ba.compute_incidence_between_cells( *it , *bd ) * 
+				ba.compute_incidence_between_cells( *bd , *bd2 );					
+			}
+        }                        
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {	
+			BOOST_CHECK( elems_in_boundary[i] == 0 );		
+		}          
+    }
+}
+
+
+
+BOOST_AUTO_TEST_CASE( compute_incidence_between_cells_test_periodic_boundary_conditions )
+{	
+    std::vector< unsigned > sizes(3);
+    sizes[0] = 5;
+    sizes[1] = 5;    
+    sizes[2] = 5;    
+    
+    std::vector<bool> directions_of_periodicity(2,true);
+    std::vector< double > data(125 , 0);
+    
+    int number_of_all_elements = (2*sizes[0])*(2*sizes[1])*(2*sizes[2]);  
+    Bitmap_cubical_complex_periodic_boundary_conditions ba( sizes , data , directions_of_periodicity );
+    for ( Bitmap_cubical_complex_periodic_boundary_conditions::All_cells_iterator it = ba.all_cells_iterator_begin() ; it != ba.all_cells_iterator_end() ; ++it )
+    {		
+        std::vector< int > elems_in_boundary( number_of_all_elements,0 );                            	        
+        Bitmap_cubical_complex_periodic_boundary_conditions_base::Boundary_range bdrange = ba.boundary_range(*it);
+        for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd = bdrange.begin() ; bd != bdrange.end() ; ++bd )
+        {
+			//std::cout <<  *bd << " ";
+            Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_range second_bdrange = ba.boundary_range(*bd);            
+            for ( Bitmap_cubical_complex_periodic_boundary_conditions::Boundary_iterator bd2 = second_bdrange.begin() ; bd2 != second_bdrange.end() ; ++bd2 )
+			{
+				elems_in_boundary[ *bd2 ] += 
+				ba.compute_incidence_between_cells( *it , *bd ) * 
+				ba.compute_incidence_between_cells( *bd , *bd2 );				
+			}
+        }                
+        //check if there is anything nonzero in elems_in_boundary
+        for ( size_t i = 0 ; i != elems_in_boundary.size() ; ++i )
+        {	
+			BOOST_CHECK( elems_in_boundary[i] == 0 );		
+		}          
+    }
+}
+