comphom.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A generic procedure for the computation of chain contraction
/// of a filtered cell complex stored in a file.
/// This is a complete procedure good for simplicial, cubical,
/// or other complexes.
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2011 by Pawel Pilarczyk.
//
// This file is part of my research software package. This is free software:
// you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// This software is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this software; see the file "license.txt". If not,
// please, see <http://www.gnu.org/licenses/>.

// Started on March 24, 2009. Last revision: April 3, 2011.


#ifndef _COMPHOM_H_
#define _COMPHOM_H_


// include some standard C++ header files
#include <istream>
#include <ostream>

// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/system/timeused.h"
#include "chomp/system/textfile.h"
#include "chomp/struct/hashsets.h"

// include relevant local header files
#include "chaincon/comblinmap.h"
#include "chaincon/filtcomplex.h"
#include "chaincon/homgvf.h"


// --------------------------------------------------
// ---------------- compute homology ----------------
// --------------------------------------------------

/// Computes the homology of a cellular complex.
/// The cellular complex is read from the given file.
/// Only top-dimensional cells should be listed in the file,
/// or otherwise the filtration may not be generated correctly.
/// If requested, the homology is only computed, without AW diagonal.
template <class CellT>
void computeHomology (const char *filename,
        chomp::homology::hashedset<CellT> &H,
        tCombLinMap<CellT, CellT> &pi,
        tCombLinMap<CellT, CellT> &incl,
        tCombLinMap<CellT, CellT> &phi,
        bool displayPi, bool displayIncl, bool displayPhi,
        bool verify, int algorithmVersion)
{
        using chomp::homology::sout;

        // read the filtered cell complex
        tFilteredComplex<CellT> K;
        {
                sout << "Reading '" << filename << "'... ";
                std::ifstream in (filename);
                if (!in)
                        chomp::homology::fileerror (filename);
                in >> K;
                in. close ();
                sout << K. size () << " cells read.\n";
        }
//      sout << "Cells read from the input:\n" << K << "\n";

        // make sure that the filtered cell complex is complete
        sout << "Adding boundaries... ";
        addBoundaries (K);
        sout << K. size () << " cells total.\n";
//      sout << "Filtered complex:\n" << K << "==========\n";

        // compute the homology together with the chain contraction
        homologyGVF (K, H, pi, incl, phi, algorithmVersion);

        // show the computed maps
        if (displayPi)
                sout << "The projection map pi:\n" << pi;
        if (displayIncl)
                sout << "The inclusion map incl:\n" << incl;
        if (displayPhi)
                sout << "The homology gradient vector field phi:\n" << phi;

        // show the computed homology representants, ordered by dimension
        sout << "Homology representants (ordered by dimension):\n";
        int_t Hsize = H. size ();
        int_t displayed = 0;
        for (int dim = 0; displayed < Hsize; ++ dim)
        {
                for (int_t i = 0; i < Hsize; ++ i)
                {
                        if (H [i]. dim () != dim)
                                continue;
                        sout << "gen_" << (displayed ++) << ": " <<
                                H [i] << "\n";
                }
        }

        // prepare the boundary map, to be used for verifications
        tCombLinMap<CellT, CellT> boundaryMap;

        // do several verifications to make sure the computed maps are good
        if (verify)
        {
                computeBoundaryMap (K, boundaryMap);
                if (phi * boundaryMap * phi == phi)
                        sout << "Verified: 'phi bd phi = phi'.\n";
                else
                        sout << "Failed to verify that "
                                "'phi pd phi = phi'.\n";

                if (boundaryMap * phi * boundaryMap == boundaryMap)
                        sout << "Verified: 'bd phi bd = bd'.\n";
                else
                        sout << "Failed to verify that 'bd phi pd = bd'.\n";

                tCombLinMap<CellT, CellT> zeroMap;
                if (phi * phi == zeroMap)
                        sout << "Verified: 'phi phi = 0'.\n";
                else
                        sout << "Failed to verify that 'phi phi = 0'.\n";

                tCombLinMap<CellT, CellT> composition (boundaryMap * phi);
                composition. add (phi * boundaryMap);
                addIdentity (K, composition);
                if (incl * pi == composition)
                        sout << "Verified: 'i pi = id - bd phi - phi bd'.\n";
                else
                        sout << "Failed to verify that "
                                "'i pi = id - bd phi - phi bd'.\n";

                tCombLinMap<CellT, CellT> idH;
                addIdentity (H, idH);
                if (pi * incl == idH)
                        sout << "Verified: 'pi i = id_H'.\n";
                else
                        sout << "Failed to verify that 'pi i = id_H'.\n";

                if (pi * phi == zeroMap)
                        sout << "Verified: 'pi phi = 0'.\n";
                else
                        sout << "Failed to verify that 'pi phi = 0'.\n";

                if (phi * incl == zeroMap)
                        sout << "Verified: 'phi i = 0'.\n";
                else
                        sout << "Failed to verify that 'phi i = 0'.\n";
        }

        return;
} /* computeHomology */


#endif // _COMPHOM_H_