atmodel2r.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Algebraic topological model computation: Algorithm 2, new version,
/// with computing the transpose of the projection (faster);
/// 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: March 24, 2011.


#ifndef _CHAINCON_ATMODEL2R_H_
#define _CHAINCON_ATMODEL2R_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 relevant local header files
#include "chaincon/chain.h"
#include "chaincon/linmap.h"
#include "chaincon/combchain.h"
#include "chaincon/comblinmap.h"


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

/// 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.
template <class CellT, class CoefT, class CellArray1, class CellArray2,
        class CellRestrT>
void algTopModel2 (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::sout;
        using chomp::homology::scon;
        using chomp::homology::sbug;

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

        // show a message indicating what is being computed
        sbug << "[AT Model Algorithm 2r.]\n";
        sout << "Computing an algebraic topological model...\n";

        // start measuring computation time
        chomp::homology::timeused compTime;
        compTime = 0;

        // prepare a set of cells that represent homology generators
        chomp::homology::hashedset<CellT> homGener;

        // prepare a combinatorial version of the transpose
        // of the projection pi, which is a matrix that stores
        // all the cells that have a given cell in their images by pi
        tCombLinMap<CellT,CellT> piT;

        // process all the remaining cells in the filter
        int_t KSize = K. size ();
        for (int_t i = 0; i < KSize; ++ i)
        {
                // show progress indicator
                if (!(i % 17))
                        scon << i << " out of " << KSize << ", " <<
                                (i * 100 / KSize) << "% done.\r";

                // retrieve the i-th cell in the filter
                const CellT &c = K [i];
        //      sout << "\nDEBUG: i = " << i << ", c = " << c << ".\n";

                // compute the modified cell
                tChain<CellT,CoefT> cChain;
                cChain. add (c, CoefT (1));
                tChain<CellT,CoefT> ch;
                ch. add (c, CoefT (1));
                ch. negate ();
                ch += phi (boundary (cChain, restr));
        //      sout << "DEBUG: bd_c = " << boundary (cChain, restr) << ".\n";
                ch. negate ();
        //      sout << "DEBUG: ch = c - phi (bd (c)) = " << ch << ".\n";

                // compute the boundary of the modified cell
                tChain<CellT,CoefT> bd_ch (boundary (ch, restr));
        //      sout << "DEBUG: bd_ch = " << bd_ch << ".\n";

                // if the boundary is zero
                if (bd_ch. empty ())
                {
                        // add this cell to the set of homology representants
                        homGener. add (c);

                        // define the inclusion for this cell
                        incl. add (c, ch);

                        // define the projection on this cell
                        pi. add (c, c, CoefT (1));

                        // add the cell to the list in the transpose
                        piT. add (c, c);

                        // skip the rest
                        continue;
                }

                // project the boundary onto the homology generators
                tChain<CellT,CoefT> d (pi (bd_ch));
                int_t dsize = d. size ();
        //      sout << "DEBUG: d = pi (bd_ch) = " << d << ".\n";

                // make sure that this projection is nontrivial
                if (d. empty ())
                {
                        sout << "Debug: c = " << c << ".\n";
                        sout << "Debug: bd_c = " <<
                                boundary (cChain, restr) << ".\n";
                        sout << "Debug: ch = " << ch << ".\n";
                        sout << "Debug: bd_ch = " << bd_ch << ".\n";
                        throw "Empty boundary in the homology algorithm 2r.";
                }

                // choose any element of this boundary
                const CellT &u = d. getCell (0);

                // compute the inverse of its coefficient
                CoefT lambdaInv (d. getCoef (0));
                lambdaInv. invert ();

                // go through all the cells which have this element
                // in their image by pi
                tCombChain<CellT> piTu (piT (u));
                int_t size = piTu. size ();
                for (int_t j = 0; j < size; ++ j)
                {
                        // retrieve the j-th cell
                        const CellT &cj = piTu. getCell (j);

                        // retrieve the image of the cell cj for editing
                        tChain<CellT,CoefT> &picj = pi. getImage (cj);

                        // determine the coefficient to use
                        int_t pos = picj. position (u);
                        if (pos < 0)
                                throw "ATmodel: u not found in pi(cj).";
                        CoefT etaLambdaInv (picj. getCoef (pos));
                        etaLambdaInv *= lambdaInv;

                        // compute new phi (cj)
                        tChain<CellT,CoefT> elich (ch);
                        elich *= etaLambdaInv;
                        phi. add (cj, elich);
                //      sout << "DEBUG: phi (" << cj << ") := " <<
                //              phi (cj) << ".\n";

                        // remember the previous image of cj by pi
                        tChain<CellT,CoefT> prev (picj);

                        // compute the new image of cj by the projection
                        tChain<CellT,CoefT> melid (d);
                        etaLambdaInv. negate ();
                        melid *= etaLambdaInv;
                        picj += melid;
                //      sout << "DEBUG: pi(" << cj << "): " << prev <<
                //              " -> " << picj << ".\n";

                        // update the transpose of the projection
                        for (int_t k = 0; k < dsize; ++ k)
                        {
                                // if this cell is in the new pi (cj)
                                const CellT &dk (d. getCell (k));
                                bool newContains =
                                        (picj. position (dk) >= 0);
                                bool prevContains =
                                        (prev. position (dk) >= 0);
                                if (newContains != prevContains)
                                        piT. add (dk, cj);
                        }
                }

                // remove the cell from the set of homology representants
                homGener. remove (u);
                incl. remove (u);
                piT. remove (u);
        }

        // copy the computed homology generators
        int_t genSize = homGener. size ();
        for (int_t i = 0; i < genSize; ++ i)
                H. push_back (homGener [i]);

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

        return;
} /* algTopModel2 */


#endif // _CHAINCON_ATMODEL2R_H_