atmodel1.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Algebraic topological model computation: Algorithm 1, new version,
/// but still without computing the transpose of the projection;
/// 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_ATMODEL1_H_
#define _CHAINCON_ATMODEL1_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 algTopModel1 (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;

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

        // show a message indicating what is being computed
        sbug << "[AT Model Algorithm 1.]\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;

        // 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 (phi (boundary (cChain, restr)));
                ch. add (c);

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

                        // skip the rest
                        continue;
                }

                // project the boundary onto the homology generators
                tCombChain<CellT> d (pi (bd_ch));

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

                // go through all the previously processed cells
                for (int_t j = 0; j < i; ++ j)
                {
                        // retrieve the j-th cell
                        const CellT &cj = K [j];

                        // compute the chain pi (cj)
                        tCombChain<CellT> cjh (pi (cj));

                        // if u does not appear in this chain
                        // then skip the remaining computations for cj
                        int_t pos = cjh. position (u);
                        if (pos < 0)
                                continue;

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

                        // compute the new image of cj by the projection
                        pi. add (cj, d);
                }

                // remove the cell from the set of homology representants
                homGener. remove (u);
                incl. 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;
} /* algTopModel1 */


#endif // _CHAINCON_ATMODEL1_H_