atmodel2.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);
/// combinatorial version - for coefficients in Z_2.
///
/////////////////////////////////////////////////////////////////////////////

// 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_ATMODEL2_H_
#define _CHAINCON_ATMODEL2_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/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 CellArray1, class CellArray2, class CellRestrT>
void algTopModel2 (const CellArray1 &K, CellArray2 &H,
        tCombLinMap<CellT,CellT> &pi,
        tCombLinMap<CellT,CellT> &incl,
        tCombLinMap<CellT,CellT> &phi,
        const CellRestrT &restr)
{
        using chomp::homology::sout;
        using chomp::homology::scon;
        using chomp::homology::sbug;
//      using chomp::homology::sseq;

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

        // show a message indicating what is being computed
        sbug << "[AT Model Algorithm 2.]\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 the transpose of the projection pi
        tCombLinMap<CellT,CellT> piT;

        // process all the 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];

                // compute the modified cell
                tCombChain<CellT> cChain;
                cChain. add (c);
                tCombChain<CellT> ch;
                ch. add (c);
                ch += phi (boundary (cChain, restr));

                // compute the boundary of the modified cell
                tCombChain<CellT> bd_ch (boundary (ch, restr));

                // 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);
                        piT. add (c, c);

                        // add the information to the 'seq' stream
                        // (to produce a demo on how it works)
                //      if (sseq. show || sseq. log)
                //              sseq << "add " << c << "\n";

                        // skip the rest
                        continue;
                }

                // project the boundary onto the homology generators
                tCombChain<CellT> d (pi (bd_ch));
                int_t dsize = d. size ();

                // 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.";
                }

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

                // add the information to the 'seq' stream
                // (to produce a demo on how it works)
        //      if (sseq. show || sseq. log)
        //              sseq << "reduce " << c << " and " << u << "\n";

                // 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);

                        // compute new phi (cj) by adding ch to it
                        phi. add (cj, ch);

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

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

                        // compute the new image of cj by the projection
                        picj += d;

                        // 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));
                        //      int_t pos = picj. position (dk);
                        //      if (pos >= 0)
                        //      {
                                        // if it wasn't there before
                                //      if (prev. position (dk) < 0)
                                        // add it to the transposed matrix
                                        piT. add (dk, cj);
                        //      }
                                // if this cell disapeared from pi (cj)
                        //      else
                        //      {
                                        // remove it from the transpose of pi
                        //              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_ATMODEL2_H_