atmodel3r.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Algebraic topological model computation: Algorithm 3, using the SNF;
/// version for an arbitrary commutative ring of coefficients.
///
/////////////////////////////////////////////////////////////////////////////

// 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: June 18, 2012.


#ifndef _CHAINCON_ATMODEL3R_H_
#define _CHAINCON_ATMODEL3R_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/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 relevant local header files
#include "chaincon/chain.h"
#include "chaincon/linmap.h"


// --------------------------------------------------
// -------------------- AT model --------------------
// --------------------------------------------------

#ifndef _product_is_identity_
#define _product_is_identity_
/// Verifies if the product of the given two matrices is an identity matrix.
/// For each column C of A, computes G(C) and checks this chain.
inline bool product_is_identity
        (const chomp::homology::mmatrix<chomp::homology::integer> &A,
        const chomp::homology::mmatrix<chomp::homology::integer> &G)
{
        if (A. getncols () != G. getnrows ())
                return false;
        if (G. getncols () != A. getnrows ())
                return false;

        using chomp::homology::integer;
        using chomp::homology::chain;
        integer one;
        one = 1;

        int_t nCols = A. getncols ();
        for (int_t nCol = 0; nCol < nCols; ++ nCol)
        {
                const chain<integer> &col (A. getcol (nCol));
                chain<integer> img;
                int_t size = col. size ();
                for (int_t pos = 0; pos < size; ++ pos)
                {
                        img. add (G. getcol (col. num (pos)),
                                col. coef (pos));
                }
                if ((img. size () != 1) ||
                        (img. num (0) != nCol) ||
                        (img. coef (0) != one))
                {
                        return false;
                }
        }
        return true;
} /* product_is_identity */
#endif // _product_is_identity_

/// A friend of the 'integer' class (from the CHomP package)
/// that converts 'integer' objects to numbers of type 'int'.
class integerFriend: public chomp::homology::integer
{
public:
        /// Constructor of the object from an 'integer' object.
        integerFriend (const chomp::homology::integer &e):
                chomp::homology::integer (e) {}

        /// Converts to the ordinary integer type.
        operator int () const {return static_cast<int> (num);}
}; /* integerFriend */

/// Translation of CHomP's 'integer' to an 'int' number.
/// This is an ugly function that uses an object of a friend class
/// in order to peek into the internals of the 'integer' object.
inline int integer2int (const chomp::homology::integer &e)
{
        return static_cast<int> (integerFriend (e));
} /* integer2int */

/// Computes an algebraic topological model for the given
/// filtered finite cell complex "K".
/// Cells that represent homology generators are appended to the vector "H".
/// The projection map "pi", the inclusion from "H" to the complex "K",
/// and the homology gradient vector field "phi"
/// are assumed to be initially zero and are constructed.
/// Note: This procedure is not well tested in case of Z_p with p != 2.
template <class CellT, class CoefT, class CellArray1, class CellArray2,
        class CellRestrT>
void algTopModel3 (const CellArray1 &K, CellArray2 &H,
        tLinMap<CellT,CellT,CoefT> &pi,
        tLinMap<CellT,CellT,CoefT> &incl,
        tLinMap<CellT,CellT,CoefT> &phi,
        const CellRestrT &restr)
{
        using chomp::homology::chaincomplex;
        using chomp::homology::chain;
        using chomp::homology::mmatrix;
        using chomp::homology::timeused;
        using chomp::homology::sout;
//      using chomp::homology::scon;
        using chomp::homology::sbug;
        using chomp::homology::integer;
        using chomp::homology::hashedset;
        using chomp::homology::multitable;

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

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

        // construct a chain complex
        sbug << "[AT Model Algorithm 3r.]\n";
        sout << "Constructing the chain complex of K...\n";
        timeused constrTime;
        constrTime = 0;

        // set up the appropriate ring of coefficients
        integer::initialize (CoefT::getp ());
        integer zero;
        zero = 0;
        integer one;
        one = 1;

        // create the chain complex
        const int_t KSize = K. size ();
        const int dim = K [KSize - 1]. dim ();
        const int trace_generators = 1;
        const int trace_bases = 1;
        chaincomplex<integer> cx (dim, trace_generators, trace_bases);

        // go through all the cells in the filter and compute boundaries
        multitable<hashedset<CellT> > cells;
        for (int_t i = 0; i < KSize; ++ i)
        {
                // retrieve the i-th cell in the filter
                const CellT &c = K [i];

                // determine the dimension of this cell
                int d = c. dim ();
                if (d < 0)
                {
                        throw "Reduced homology not supported "
                                "in Algorithm 3.\n"
                                "Please, use Algorithm 4 instead.";
                }

                // add the cell to the right array of cells
                int_t cNum = cells [d]. size ();
                cells [d]. add (c);

                // 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 number of the boundary cell
                        int_t bNum = cells [d - 1]. getnumber (bCell);

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

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

                        // add the boundary relation to the chain complex
                        cx. add (d, bNum, cNum, bCoef);
                }
        }

        // make sure that the size of the chain complex is enough
        // to contain all the cells
        for (int d = 0; d <= dim; ++ d)
        {
                cx. def_gen (d, cells [d]. size ());
        }

        // show timing information
        sout << "The chain complex created in " << constrTime << ".\n";

        // compute the SNF
        sout << "Computing the SNF of the boundary operator...\n";
        timeused snfTime;
        snfTime = 0;

        int *level = 0;
        bool quiet = false;
        cx. simple_form (level, quiet);

//      for (int q = dim; q >= 0; -- q)
//      {
//              sbug << "DEBUG: Simple form [" << q << "]:\n" <<
//                      cx. getboundary (q);
//      }

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

        // DEBUG: Check if the change of basis matrices are correct.
        if (true)
        {
                timeused verifTime;
                verifTime = 0;
                sout << "Change of basis check: ";
                for (int q = 0; q <= dim; ++ q)
                {
                        const mmatrix<integer> &G (cx. gethomgen (q));
                        const mmatrix<integer> &A (cx. getchgbasis (q));
                        if (!product_is_identity (G, A) ||
                                !product_is_identity (A, G))
                        {
                                sout << "Failed at level " << q << ".\n";
                                throw "Wrong change of basis detected.";
                        }
                }
                sout << "Completed successfully in " << verifTime << ".\n";
        }

        // compute the g.v.f. together with the projection and the inclusion
        sout << "Composing an algebraic topological model from the SNF...\n";
        timeused composingTime;
        composingTime = 0;

        // prepare the lists of chain complex generators
        // which are also homology generators
//      multitable <hashedset<int_t> > homGener;

        // prepare the matrices for psi
        multitable <mmatrix<integer> > psi;

        // scan through the entire chain complex in the SNF
        const mmatrix<integer> zeroMatrix;
        for (int q = 0; q <= dim; ++ q)
        {
                // prepare a number for subsequent generator representants
                int_t gen_number = 0;

                // get the matrices computed for the Smith Normal Form:
                // the boundary at level q
                const mmatrix<integer> &Dq = cx. getboundary (q);
                // the boundary at level q + 1
                const mmatrix<integer> &Dq1 = (q < dim) ?
                        cx. getboundary (q + 1) : zeroMatrix;
                // the change of basis from the new one to the original one
                const mmatrix<integer> &Gq = cx. gethomgen (q);
                // the change of basis from the original one to the new one
                const mmatrix<integer> &Aq = cx. getchgbasis (q);

                // set the size of the matrix psi in dimension q - 1
                if (q)
                {
                        psi [q - 1]. define (cells [q]. size (),
                                cells [q - 1]. size ());
                }

                // process the bases of the chain complex in the SNF
                int ncells = cells [q]. size ();
                for (int i = 0; i < ncells; ++ i)
                {
                        // get the column of the matrix Dq
                        const chain<integer> &column = Dq. getcol (i);

                        // if the column in the SNF is nonzero
                        if (!column. empty ())
                        {
                                integer e (column. coef (0));
                                if (e. delta () != 1)
                                {
                                        throw "Non-invertible coefficient "
                                                "encountered in the SNF.";
                                }
                                int b = column. num (0);
                                if (!q)
                                        throw "Negative psi index.";
                                psi [q - 1]. add (i, b, e /* one? */);
                        }

                        // if the column in the SNF is zero
                        else if ((q == dim) || (Dq1. getrow (i). empty ()))
                        {
                                // add the number of this cell
                                // to the list of homology generators
                        //      homGener [q]. add (i);

                                // determine the real homology generator
                                const chain<integer> &genChain =
                                        Gq. getcol (i);

                                // define a representant of this generator
                                const CellT &repr = cells [q] [gen_number];
                                ++ gen_number;
                                H. push_back (repr);

                                // compute the inclusion map on it
                                int_t genSize = genChain. size ();
                                for (int_t j = 0; j < genSize; ++ j)
                                {
                                        int_t n = genChain. num (j);
                                        integer e = genChain. coef (j);
                                        const CellT &genCell = cells [q] [n];
                                        incl. add (repr, genCell,
                                                CoefT (integer2int (e)));
                                }

                                // add this term to the projection
                                const chain<integer> &piRow =
                                        Aq. getrow (i);
                                int_t piRowSize = piRow. size ();
                                for (int_t j = 0; j < piRowSize; ++ j)
                                {
                                        int_t num = piRow. num (j);
                                        integer e = piRow. coef (j);
                                        const CellT &cell = cells [q] [num];
                                        pi. add (cell, repr,
                                                CoefT (integer2int (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 maps to the original basis...\n";
        for (int q = 1; q <= dim; ++ q)
        {
                // consider the three matrices to be multiplied
                const mmatrix<integer> &M_A = cx. getchgbasis (q - 1);
                const mmatrix<integer> &M_psi = psi [q - 1];
                const mmatrix<integer> &M_A1 = cx. gethomgen (q);

                // determine the numbers of columns in the matrices
                int_t sizeA = M_A. getncols ();
                int_t sizePsi = M_psi. getncols ();
                int_t sizeA1 = M_A1. getncols ();

                // for all the columns in the first matrix on the right...
                for (int_t i = 0; i < sizeA; ++ i)
                {
                        // take the column
                        const chain<integer> &column = M_A. getcol (i);

                        // compute its image by psi
                        chain<integer> psiC;
                        int_t sizeC = column. size ();
                        for (int_t j = 0; j < sizeC; ++ j)
                        {
                                int_t n = column. num (j);
                                if (n >= sizePsi)
                                        continue;
                                const integer &e = column. coef (j);
                                psiC. add (M_psi. getcol (n), e /* one? */);
                        }

                        // compute the image of this image by A1
                        chain<integer> phiColumn;
                        int_t sizePsiC = psiC. size ();
                        for (int j = 0; j < sizePsiC; ++ j)
                        {
                                int_t n = psiC. num (j);
                                if (n >= sizeA1)
                                        continue;
                                const integer &e = psiC. coef (j);
                                phiColumn. add (M_A1. getcol (n), e /* one? */);
                        }

                        // decode the resulting chain and put it in phi
                        int_t sizePhiColumn = phiColumn. size ();
                        for (int j = 0; j < sizePhiColumn; ++ j)
                        {
                                int_t n = phiColumn. num (j);
                                integer e = phiColumn. coef (j);
                                phi. add (cells [q - 1] [i], cells [q] [n],
                                        CoefT (integer2int (e)));
                        }
                }
        }

        // show timing information
        sout << "AT model composed from the SNF in " << composingTime <<
                ".\n";

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

        return;
} /* algTopModel3 */


#endif // _CHAINCON_ATMODEL3R_H_