ammodel.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Computation of the algebraic minimal model using the SNF.
/// This algorithm works for an arbitrary commutative ring of coefficients
/// that is an Euclidean domain.
///
/////////////////////////////////////////////////////////////////////////////

// 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: January 13, 2013.


#ifndef _CHAINCON_AMMODEL_H_
#define _CHAINCON_AMMODEL_H_


// include relevant local header files
#include "chaincon/chain.h"
#include "chaincon/linmap.h"
#include "chaincon/comblinmap.h"
#include "chaincon/filtcomplex.h"
#include "chaincon/algcell.h"
#include "chaincon/snf.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 "chomp/struct/multitab.h"
#include "chomp/struct/integer.h"
#include "chomp/homology/chains.h"

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


// --------------------------------------------------
// -------------------- AM Model --------------------
// --------------------------------------------------

/// Computes a minimal model for a given filtered finite cell complex "K".
/// Cells that represent those generators of M whose boundary and coboundary
/// are both zero are appended to the vector "H"; note that they represent
/// (co)homology generators corresponding to the free part
/// of the (co)homology module.
/// The other generators of "M" are appended to the vectors "A" and "B"
/// with the corresponding coefficients stored in the vector "Q"
/// in such a way that the boundary of "A[i]" equals the coefficient "Q[i]"
/// times "B[i]".
/// The projection map "pi" onto M, the inclusion from "H+A+B" to "K",
/// and the chain contraction "phi" on "K" are assumed to be initially zero
/// and are constructed.
template <class CellT, class CoefT, class CellArray1, class CellArray2,
        class CellArray3, class CoefArray>
void algMinModel (const tFilteredComplex<CellT> &K, CellArray1 &H,
        CellArray2 &A, CellArray3 &B, CoefArray &Q,
        tLinMap<CellT,CellT,CoefT> &pi,
        tLinMap<CellT,CellT,CoefT> &incl,
        tLinMap<CellT,CellT,CoefT> &phi)
{
        using chomp::homology::sout;
//      using chomp::homology::scon;
        using chomp::homology::sbug;
        using chomp::homology::timeused;
        using chomp::homology::multitable;
        using chomp::homology::hashedset;

        typedef typename tMatrixSNF<CoefT>::CellType AlgCellType;
        typedef typename tMatrixSNF<CoefT>::ChainType AlgChainType;

        // don't do anything for the empty complex
        if (K. empty ())
                return;

        // start measuring total computation time of this routine
        timeused compTime;
        compTime = 0;

        // split K into cells of each dimension
        sout << "Splitting K into sets of cells of each dimension...\n";
        timeused splitTime;
        splitTime = 0;

        // prepare an array of sets of cells
        multitable<hashedset<CellT> > cells;

        // prepare a strict upper bound on the dimension of the cells
        int maxDim = 0;

        // split the set of cells in the filtration into the subsets
        const int_t KSize = K. size ();
        for (int_t k = 0; k < KSize; ++ k)
        {
                const CellT &cell (K [k]);
                int dim (cell. dim ());
                if (dim >= 0)
                        cells [dim]. add (cell);
                if (dim >= maxDim)
                        maxDim = dim + 1;
        }

        // show timing information
        sout << "K split in " << splitTime << ".\n";

        // construct the matrices of the boundary operator
        sout << "Creating the matrices of the boundary operator on K...\n";
        timeused constrTime;
        constrTime = 0;

        // create the collection of matrices for the boundary operator
        tMatrixSNF<CoefT> matrix;

        // set the correct sizes of the matrices
        for (int q = 0; q < maxDim; ++ q)
        {
                int_t nRows = (q > 0) ? cells [q - 1]. size () :
                        (CellT::EmptyType::exists () ? 1 : 0);
                int_t nCols = cells [q]. size ();
                matrix. setSize (q, nRows, nCols);
        //      sbug << "DEBUG: Boundary " << q << " matrix: " <<
        //              nRows << " rows, " << nCols << " columns.\n";
        }

        // go through all the cells in the filter and compute boundaries
        for (int q = 0; q < maxDim; ++ q)
        {
                if ((q == 0) && !CellT::EmptyType::exists ())
                        continue;
                int_t size (cells [q]. size ());
                for (int_t cNum = 0; cNum < size; ++ cNum)
                {
                        // retrieve the i-th cell in the filter
                        const CellT &c = cells [q] [cNum];

                        // compute the boundary of this cell
                        int_t bLen = c. boundaryLength ();
                        for (int_t j = 0; j < bLen; ++ j)
                        {
                                // create the j-th boundary cell
                                CellT bCell (c, j);

                                // determine the coefficient
                                CoefT bCoef (c. boundaryCoef (j));

                                // determine the number of the boundary cell
                                int_t bNum = (q > 0) ?
                                        cells [q - 1]. getnumber (bCell) : 0;

                                // ignore the cell is it is not in the set
                                if (bNum < 0)
                                        continue;

                                // add the boundary coefficient to the matrix
                        //      sbug << "DEBUG: Boundary matrix " << q <<
                        //              " add " << bCoef << " at row " <<
                        //              bNum << " and col " << cNum << ".\n";
                                matrix. add (q, bNum, cNum, bCoef);
                        }
                }
        }

        // show timing information
        sout << "The boundary matrices created in " << constrTime << ".\n";

        // compute the SNF
        sout << "Computing the SNF of the boundary operator...\n";
        timeused snfTime;
        snfTime = 0;
        matrix. computeSNF (maxDim);

        // show timing information
        sout << "The SNF computed in " << snfTime << ".\n";

        // DEBUG: Check if the change of basis matrices are correct.
        if (true)
                verifyChangeOfBasisSNF (matrix, maxDim);

        // compute the g.v.f. together with the projection and the inclusion
        sout << "Composing a chain contraction from the SNF...\n";
        timeused composingTime;
        composingTime = 0;

        // prepare the matrices for psi
        typedef tLinMap<AlgCellType,AlgCellType,CoefT> AlgMapType;
        multitable<AlgMapType> psi;

        // prepare matrices for storing representant cells for generators
        multitable<tCombLinMap<AlgCellType,CellT> > representants;

        // prepare data for storing torsion relations:
        // torsionUp (cell) = higher-dimensional cell with the torsion coef
//      multitable<AlgMapType> torsionUp;

        // scan through the boundary matrices in the SNF
        for (int q = 0; q < maxDim; ++ q)
        {
                // prepare a counter of homology generator representants
                int_t genNumber (0);

                // scan the coefficients in this matrix
                int_t nCols = cells [q]. size ();
                for (int_t cell = 0; cell < nCols; ++ cell)
                {
                        // determine the row with a nonzero coef if any
                        int_t bcell = matrix. getColCell (q, cell);
                        bool torsionCell = false;
                //      CoefT torsionCoef (0);

                        // if the column is non-zero
                        if (bcell >= 0)
                        {
                                const CoefT e (matrix. getColCoef (q, cell));
                                const int delta (e. delta ());
                                if (delta == 1)
                                {
                                        const AlgCellType dom (bcell);
                                        const AlgCellType img (cell);
                                        if (q)
                                                psi [q - 1]. add (dom, img, e);
                                        continue;
                                }
                                else
                                {
                                        const AlgCellType c (bcell);
                                        B. push_back (representants [q - 1].
                                                getImage (c). getCell (0));
                                        Q. push_back (e);
                                        torsionCell = true;
                                //      torsionCoef = e;
                                }
                        }

                        // otherwise check the preimage of the cell
                        // by the boundary operator in the SNF
                        bool torsionImage = false;
                        if (q < maxDim - 1)
                        {
                                CoefT e (matrix. getRowCoef (q + 1, cell));

                                // if the cell is a boundary then skip it
                                if (e. delta () == 1)
                                        continue;

                                // if the cell is an image of torsion part
                                else if (e. delta () > 1)
                                {
                                //      const AlgCellType higher (matrix.
                                //              getRowCell (q + 1, cell));
                                //      torsionUp [q]. add (cell, higher, e);
                                        torsionImage = true;
                                //      torsionCoef = e;
                                }
                        }

                        // define a representant of this generator
                        const CellT &repr = cells [q] [genNumber];
                        representants [q]. add (AlgCellType (cell), repr);
                        ++ genNumber;
                        if (torsionCell)
                                A. push_back (repr);
                        else if (!torsionImage)
                                H. push_back (repr);

                        // determine the real homology generator
                        const AlgChainType genChain
                                (matrix. getOldChain (q, cell));

                        // compute the inclusion map on this generator
                        int_t genSize = genChain. size ();
                        for (int_t j = 0; j < genSize; ++ j)
                        {
                                int_t n = genChain. getCell (j). getId ();
                                CoefT e = genChain. getCoef (j);
                                const CellT &genCell = cells [q] [n];
                                incl. add (repr, genCell, e);
                        }

                        // add this generator to the image of all the cells
                        // by the projection, whenever the cells have it
                        // in their images by the change-of-basis matrix
                        const AlgChainType piRow
                                (matrix. getNewChainRow (q, cell));
                        int_t piRowSize = piRow. size ();
                        for (int_t j = 0; j < piRowSize; ++ j)
                        {
                                int_t n = piRow. getCell (j). getId ();
                                CoefT e = piRow. getCoef (j);
                                const CellT &cell = cells [q] [n];
                                pi. add (cell, repr, e);
                        }
                }
        }

        // compute the matrices for phi_q and compose them into one map:
        // \phi_{q - 1} := A^{-1}_q \circ \psi_{q - 1} \circ A_{q - 1}
        sout << "Transferring the computed chain contraction "
                "to the original basis...\n";
        for (int q = 1; q < maxDim; ++ q)
        {
                // for all the cells of dimension q - 1
                int_t nCells (cells [q - 1]. size ());
                for (int_t cell = 0; cell < nCells; ++ cell)
                {
                        // compute the chain in the new basis
                        AlgChainType ch (matrix. getNewChain (q - 1, cell));

                        // apply psi to this chain
                        AlgChainType psich (psi [q - 1] (ch));

                        // compute the resulting chain in the original basis
                        AlgChainType algImage;
                        int_t size (psich. size ());
                        for (int_t i = 0; i < size; ++ i)
                        {
                                // compute the algebraic chain for this part
                                int_t n (psich. getCell (i). getId ());
                                const CoefT &e (psich. getCoef (i));
                                AlgChainType a (matrix. getOldChain (q, n));
                                a *= e;
                                algImage += a;
                        }

                        // translate the algebraic chain to the real one
                        tChain<CellT,CoefT> image;
                        int_t algSize (algImage. size ());
                        for (int_t i = 0; i < algSize; ++ i)
                        {
                                const CoefT &coef (algImage. getCoef (i));
                                int_t n (algImage. getCell (i). getId ());
                                const CellT &cell (cells [q] [n]);
                                image. add (cell, coef);
                        }

                        // set the image in the map phi
                        phi. getImage (cells [q - 1] [cell]). swap (image);
                }
        }

        // show timing information
        sout << "Chain contraction composed in " << composingTime << ".\n";

        // show information on the total computation time of this routine
        sout << "AM model computed in " << compTime << ".\n";

        return;
} /* algMinModel */


#endif // _CHAINCON_AMMODEL_H_