CAPD RedHom Library
C++ API examples

Compute Homology in C++

Cubical complexes

CubicalBettiNumbers.cpp

/////////////////////////////////////////////////////////////////////////////
/// @file CubicalBettiNumbers
///
/// @author Mateusz Juda <mateusz.juda@{ii.uj.edu.pl,gmail.com}>
///
/// @date 2015-04-26
/////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2000-2015 by the CAPD Group.
//
// This file constitutes a part of the CAPD library (capdRedHom),
// distributed under the terms of the GNU General Public License.
// Consult http://capd.ii.uj.edu.pl and  http://redhom.ii.edu.pl/ for details.
/////////////////////////////////////////////////////////////////////////////

#include "capd/apiRedHom/ComplexHelper.h"
#include "capd/apiRedHom/Algorithms.h"

#include <string>
#include <iostream>

using namespace capd::apiRedHom;

std::vector<CubicalHelper::FullCube> kleinCubes = {
   {12, 6, 6, 6}, {11, 6, 6, 6}, {11, 6, 6, 7}, {10, 6, 6, 7}, {9, 6, 6, 7},  {8, 6, 6, 7},  {8, 6, 6, 6},
   {8, 6, 6, 5},  {9, 6, 6, 5},  {10, 6, 6, 5}, {11, 6, 6, 5}, {11, 7, 6, 6}, {11, 7, 6, 7}, {11, 7, 6, 5},
   {11, 8, 6, 6}, {10, 8, 6, 7}, {10, 7, 6, 7}, {9, 7, 6, 7},  {8, 7, 6, 7},  {7, 6, 6, 6},  {8, 7, 6, 5},
   {9, 7, 6, 5},  {10, 7, 6, 5}, {10, 8, 6, 5}, {10, 9, 6, 6}, {10, 9, 6, 7}, {9, 8, 6, 7},  {7, 7, 6, 7},
   {7, 7, 6, 6},  {7, 7, 6, 5},  {9, 8, 6, 5},  {10, 9, 6, 5}, {9, 10, 6, 6}, {9, 9, 6, 7},  {8, 8, 6, 7},
   {8, 8, 6, 5},  {9, 9, 6, 5},  {8, 10, 6, 6}, {8, 10, 6, 7}, {8, 9, 6, 7},  {7, 9, 6, 7},  {7, 8, 6, 7},
   {7, 8, 6, 5},  {7, 9, 6, 5},  {8, 9, 6, 5},  {8, 10, 6, 5}, {7, 11, 7, 6}, {7, 10, 6, 7}, {6, 9, 6, 7},
   {6, 8, 6, 7},  {6, 8, 5, 6},  {6, 8, 6, 5},  {6, 9, 6, 5},  {7, 10, 6, 5}, {6, 11, 7, 6}, {6, 10, 6, 7},
   {6, 8, 5, 7},  {6, 8, 5, 5},  {6, 10, 6, 5}, {5, 10, 7, 6}, {5, 10, 7, 7}, {5, 10, 6, 7}, {5, 9, 6, 7},
   {5, 9, 6, 5},  {5, 10, 6, 5}, {5, 10, 7, 5}, {4, 10, 7, 6}, {4, 10, 7, 7}, {4, 9, 6, 7},  {5, 8, 5, 7},
   {5, 8, 5, 6},  {5, 8, 5, 5},  {4, 9, 6, 5},  {4, 10, 7, 5}, {3, 9, 7, 6},  {3, 9, 7, 7},  {4, 8, 6, 7},
   {4, 8, 5, 7},  {4, 8, 5, 6},  {4, 8, 5, 5},  {4, 8, 6, 5},  {3, 9, 7, 5},  {3, 8, 7, 6},  {3, 8, 7, 7},
   {3, 8, 6, 7},  {4, 7, 5, 6},  {3, 8, 6, 5},  {3, 8, 7, 5},  {2, 7, 7, 6},  {3, 7, 7, 7},  {3, 7, 6, 7},
   {3, 7, 5, 7},  {3, 7, 5, 6},  {3, 7, 5, 5},  {3, 7, 6, 5},  {3, 7, 7, 5},  {2, 6, 7, 6},  {2, 6, 7, 7},
   {3, 6, 6, 7},  {3, 6, 5, 7},  {3, 6, 5, 6},  {3, 6, 5, 5},  {3, 6, 6, 5},  {2, 6, 7, 5},  {2, 6, 8, 6},
   {2, 6, 6, 7},  {2, 6, 6, 5},  {3, 6, 7, 6},  {3, 6, 7, 7},  {2, 6, 5, 7},  {2, 6, 5, 6},  {2, 6, 5, 5},
   {3, 6, 7, 5},  {3, 5, 7, 6},  {3, 5, 7, 7},  {3, 5, 6, 7},  {3, 5, 5, 7},  {2, 5, 5, 6},  {3, 5, 5, 5},
   {3, 5, 6, 5},  {3, 5, 7, 5},  {4, 5, 7, 6},  {4, 4, 7, 6},  {4, 4, 7, 7},  {3, 4, 6, 7},  {3, 4, 5, 7},
   {3, 4, 5, 6},  {3, 4, 5, 5},  {3, 4, 6, 5},  {4, 4, 7, 5},  {4, 4, 6, 7},  {4, 3, 6, 7},  {3, 3, 5, 7},
   {3, 3, 5, 6},  {3, 3, 5, 5},  {4, 3, 6, 5},  {4, 4, 6, 5},  {5, 4, 7, 6},  {5, 4, 7, 7},  {5, 3, 6, 7},
   {4, 2, 5, 7},  {4, 2, 5, 6},  {4, 2, 5, 5},  {5, 3, 6, 5},  {5, 4, 7, 5},  {6, 4, 7, 6},  {6, 3, 7, 7},
   {5, 2, 6, 7},  {5, 2, 5, 7},  {5, 2, 5, 6},  {5, 2, 5, 5},  {5, 2, 6, 5},  {6, 3, 7, 5},  {6, 3, 6, 7},
   {6, 2, 6, 7},  {6, 1, 5, 7},  {6, 1, 5, 6},  {6, 1, 5, 5},  {6, 2, 6, 5},  {6, 3, 6, 5},  {6, 4, 6, 7},
   {6, 4, 6, 5},  {7, 3, 6, 7},  {7, 2, 6, 7},  {7, 1, 6, 7},  {7, 1, 5, 6},  {7, 1, 6, 5},  {7, 2, 6, 5},
   {7, 3, 6, 5},  {7, 5, 6, 6},  {7, 4, 6, 6},  {7, 4, 6, 7},  {8, 3, 6, 7},  {8, 2, 6, 7},  {8, 2, 6, 6},
   {8, 2, 6, 5},  {8, 3, 6, 5},  {7, 4, 6, 5},  {8, 4, 6, 7},  {9, 3, 6, 7},  {9, 2, 6, 7},  {9, 2, 6, 6},
   {9, 2, 6, 5},  {9, 3, 6, 5},  {8, 4, 6, 5},  {8, 5, 6, 7},  {9, 4, 6, 7},  {10, 4, 6, 7}, {10, 3, 6, 7},
   {10, 3, 6, 6}, {10, 3, 6, 5}, {10, 4, 6, 5}, {9, 4, 6, 5},  {8, 5, 6, 5},  {9, 5, 6, 7},  {11, 4, 6, 7},
   {11, 4, 6, 6}, {11, 4, 6, 5}, {9, 5, 6, 5},  {10, 5, 6, 7}, {11, 5, 6, 7}, {11, 5, 6, 6}, {11, 5, 6, 5},
   {10, 5, 6, 5}};

void KleinBetti()
{
   std::vector<int> betti = CubicalHelper::BettiNumbers(kleinCubes);
   std::vector<int> expectedBetti = {1, 1};
   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

void KleinBettiZ2()
{
   std::vector<int> betti = CubicalHelper::BettiNumbers(kleinCubes, 2);
   std::vector<int> expectedBetti = {1, 2, 1};
   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

void KleinBettiZp()
{
   std::vector<int> betti = CubicalHelper::BettiNumbers(kleinCubes, 1019);
   std::vector<int> expectedBetti = {1, 1};
   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

void BuildComplexBetti()
{
   std::vector<size_t> complexSize = {20, 20, 20, 20};  // size in each dimension
   CubicalComplex complex(complexSize);

   for (auto cube : kleinCubes) {
      complex.insert(cube);  // insert full cube
   }

   complex.fillWithBoundaries();  // generates all cubes in boundary
   std::vector<int> betti = ComputeBettiNumbers(complex);
   std::vector<int> expectedBetti = {1, 1};

   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

void BuildComplexBettiZ2()
{
   std::vector<size_t> complexSize = {20, 20, 20, 20};  // size in each dimension
   CubicalComplex complex(complexSize);

   for (auto cube : kleinCubes) {
      complex.insert(cube);  // insert full cube
   }

   complex.fillWithBoundaries();  // generates all cubes in boundary
   std::vector<int> betti = ComputeBettiNumbers(complex, 2);
   std::vector<int> expectedBetti = {1, 2, 1};

   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

void BuildComplexBettiZp()
{
   std::vector<size_t> complexSize = {20, 20, 20, 20};  // size in each dimension
   CubicalComplex complex(complexSize);

   for (auto cube : kleinCubes) {
      complex.insert(cube);  // insert full cube
   }

   complex.fillWithBoundaries();  // generates all cubes in boundary
   std::vector<int> betti = ComputeBettiNumbers(complex, 1019);
   std::vector<int> expectedBetti = {1, 1};

   assert(expectedBetti == betti);
   std::cout << "OK\n";
}

int main()
{
   try {
      KleinBetti();
      KleinBettiZp();
      KleinBettiZp();
      BuildComplexBetti();
      BuildComplexBettiZ2();
      BuildComplexBettiZp();

      return 0;
   }
   catch (std::exception& ex) {
      std::cerr << "Error: " << ex.what() << std::endl;
      return 1;
   }
   catch (const char* ex) {
      std::cerr << "Error: " << ex << std::endl;
      return 1;
   }
   catch (const std::string& ex) {
      std::cerr << "Error: " << ex << std::endl;
      return 1;
   }
   catch (...) {
      std::cerr << "Unknown exception" << std::endl;
      return 1;
   }

   return 1;
}

Cubical complexes in parallel

CubicalBettiNumbersParallel.cpp

/////////////////////////////////////////////////////////////////////////////
/// @file CubicalBettiNumbersParallel
///
/// @author Mateusz Juda <mateusz.juda@{ii.uj.edu.pl,gmail.com}>
///
/// @date 2016-06-01
/////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2000-2016 by the CAPD Group.
//
// This file constitutes a part of the CAPD library (capdRedHom),
// distributed under the terms of the GNU General Public License.
// Consult http://capd.ii.uj.edu.pl and  http://redhom.ii.edu.pl/ for details.
/////////////////////////////////////////////////////////////////////////////

#include <capd/apiRedHom/Algorithms.h>
#include <capd/apiRedHom/Complex.h>
#include <capd/apiRedHom/TheConfig.h>
#include <iostream>

using namespace capd;

int main(int /*argc*/, char** /*argv*/)
{
   try {
      const size_t size = 50;
      std::cout << "Create empty GridGraphCubicalComplex, size: " << size << std::endl;
      // std::vector<uint64_t> memory((2 * (size + 1) - 1) * (2 * (size + 1) - 1) * (2 * (size + 1) - 1), 1);
      // apiRedHom::CubicalComplex complex({size, size, size},
      //                                   boost::make_iterator_range(&memory.front(), &memory.back() + 1), 41);

      apiRedHom::CubicalComplex complex({size, size, size}, apiRedHom::CubicalComplex::CellCodeRange(), 40);
      std::cout << "Insert some cells" << std::endl;
      for (size_t z = 0; z < 2 * size + 1; ++z) {
         for (size_t y = 0; y < 2 * size + 1; ++y) {
            for (size_t x = 0; x < 2 * size + 1; ++x) {
               if ((x * y * z) % 100 != 0) {
                  complex.insert({{x / 2, x % 2}, {y / 2, y % 2}, {z / 2, z % 2}});
               }
            }
         }
      }

      apiRedHom::TheConfig::instance().setDefaultComputationModel(2);  // Make sure we use TBB == 2, sequential == 1
      // Make sure we allow parallel implementation with connected components, 4-th position responsible for it
      apiRedHom::TheConfig::instance().setBettiNumbersFlags("11111");
      // set threashold for recursive call
      apiRedHom::TheConfig::instance().setULong("ComputeConnectedComponentsOfGridGraph.threashold", 500000ll);

      std::cout << "Start algorithm" << std::endl;
      std::vector<int> betti = apiRedHom::ComputeBettiNumbers(complex);
      std::cout << "Betti numbers:\n";

      for (auto b : betti) {
         std::cout << b << " ";
      }

      return 0;
   }
   catch (std::exception& ex) {
      std::cerr << "Error: " << ex.what() << std::endl;
      return 1;
   }
   catch (const char* ex) {
      std::cerr << "Error: " << ex << std::endl;
      return 1;
   }
   catch (const std::string& ex) {
      std::cerr << "Error: " << ex << std::endl;
      return 1;
   }
   catch (...) {
      std::cerr << "Unknown exception" << std::endl;
      return 1;
   }

   return 1;
}