ammodcomp.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A generic procedure for the computation of an algebraic minimal model
/// of a filtered cell complex stored in a file.
/// This is a complete procedure good for simplicial, cubical,
/// or other complexes.
/// It computes a chain contraction onto a cellular complex
/// with the boundary operator containing only the torsion part.
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2016 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: August 21, 2013.


#ifndef _AMMODCOMP_H_
#define _AMMODCOMP_H_


// include relevant local header files
#include "chaincon/cellnames.h"
#include "chaincon/filtcomplex.h"
#include "chaincon/readfcompl.h"
#include "chaincon/linmap.h"
#include "chaincon/boundary.h"
#include "chaincon/ammodel.h"
#include "chaincon/homcohom.h"

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

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


// --------------------------------------------------
// ---------------- 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.
template <class CellT, class CoefT, class CellArray1, class CellArray2,
        class CellArray3, class CoefArray, class CellNames>
inline void computeAlgMinModel (const tFilteredComplex<CellT> &K,
        bool relativeComplex,
        CellArray1 &H, CellArray2 &A, CellArray3 &B, CoefArray &Q,
        tLinMap<CellT,CellT,CoefT> &pi,
        tLinMap<CellT,CellT,CoefT> &incl,
        tLinMap<CellT,CellT,CoefT> &phi,
        CellNames &homCellNames, CellNames &cohomCellNames,
        bool displayPi, bool displayIncl, bool displayPhi,
        bool displayD, bool verify)
{
        using chomp::homology::sout;

//      typedef typename LinMap::CellDomType CellT;
        typedef tLinMap<CellT,CellT,CoefT> LinMap;

        // compute the minimal model together with the chain contraction
        algMinModel (K, H, A, B, Q, pi, incl, phi);

        // show the computed homology groups
        std::vector<std::vector<CellT> > representants;
        std::vector<int> bettiNumbers;
        std::vector<std::vector<CoefT> > torsion;
        getHomCohom (H, A, B, Q, &representants, &bettiNumbers, &torsion,
                false);
        sout << "Homology groups:\n";
        sout << homcohom2text (bettiNumbers, torsion, "H_", false);

        // show the computed cohomology groups
        std::vector<std::vector<CellT> > cohomRepresentants;
        std::vector<int> cohomBettiNumbers;
        std::vector<std::vector<CoefT> > cohomTorsion;
        getHomCohom (H, A, B, Q, &cohomRepresentants, &cohomBettiNumbers,
                &cohomTorsion, true);
        sout << "Cohomology groups:\n";
        sout << homcohom2text (cohomBettiNumbers, cohomTorsion, "H^", false);

        // prepare correct cell names for homology representants
        int homReprSize = representants. size ();
        for (int dim = 0; dim < homReprSize; ++ dim)
        {
                int size = representants [dim]. size ();
                for (int i = 0; i < size; ++ i)
                        homCellNames (representants [dim] [i]);
        }

        // prepare correct cell names for cohomology representants
        int cohomReprSize = cohomRepresentants. size ();
        for (int dim = 0; dim < cohomReprSize; ++ dim)
        {
                int size = cohomRepresentants [dim]. size ();
                for (int i = 0; i < size; ++ i)
                        cohomCellNames (cohomRepresentants [dim] [i]);
        }

        // show the computed maps if requested to
        tCellOwnNames<CellT> cellOwnNames;
        if (displayPi)
        {
                sout << "The projection map pi:\n" <<
                        cells2names (pi, cellOwnNames, homCellNames);
        }
        if (displayIncl)
        {
                sout << "The inclusion map incl:\n" <<
                        cells2names (incl, homCellNames, cellOwnNames);
        }
        if (displayPhi)
        {
                sout << "The chain contraction phi:\n" << phi;
        }
        if (displayD)
        {
                sout << "The boundary map in the minimal model:\n";
                for (int_t i = 0; i < A. size (); ++ i)
                {
                        sout << homCellNames (A [i]) << " -> " << Q [i] <<
                                " * " << homCellNames (B [i]) << "\n";
                }
        }

        // prepare the boundary map, to be used for verifications
        LinMap boundaryMap;

        // do several verifications to make sure the computed maps are good
        if (verify)
        {
                if (relativeComplex)
                        computeBoundaryMap (K, boundaryMap, K);
                else
                        computeBoundaryMap (K, boundaryMap, Unrestricted ());
                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 bd = bd', "
                                "which is OK if there is torsion.\n";
                }

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

                LinMap composition (boundaryMap * phi);
                composition += phi * boundaryMap;
                composition. negate ();
                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";

                LinMap idHAB;
                addIdentity (H, idHAB);
                addIdentity (A, idHAB);
                addIdentity (B, idHAB);
                if (pi * incl == idHAB)
                        sout << "Verified: 'pi i = id_H+A+B'.\n";
                else
                        sout << "Failed to verify that 'pi i = id_H+A+B'.\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;
} /* computeAlgMinModel */


#endif // _AMMODCOMP_H_