chaincon/homgvf2.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Homology gradient vector field computation: Algorithm 2, new version,
/// with computing the transpose of the projection (faster).
///
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2009-2011 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_HOMGVF2_H_
#define _CHAINCON_HOMGVF2_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/timeused.h"
#include "chomp/system/textfile.h"
#include "chomp/struct/hashsets.h"

// include relevant local header files
#include "chaincon/combchain.h"
#include "chaincon/comblinmap.h"


// --------------------------------------------------
// ------------------ Homology GVF ------------------
// --------------------------------------------------

/// Computes the homology gradient vector field 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>
void homologyGVF2 (const CellArray1 &K, CellArray2 &H,
        tCombLinMap<CellT,CellT> &pi,
        tCombLinMap<CellT,CellT> &incl,
        tCombLinMap<CellT,CellT> &phi)
{
        using chomp::homology::sout;
        using chomp::homology::scon;

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

        // show a message indicating what is being computed
        sout << "Computing a homology gradient vector field...\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 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];

                // compute the modified cell
                tCombChain<CellT> ch (c);
                ch. add (phi (boundary (c)));

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

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

                        // 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 (c) <<
                                ".\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 [0];

                // 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 [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. add (d);

                        // update the transpose of the projection
                        for (int_t k = 0; k < dsize; ++ k)
                        {
                                // if this cell is in the new pi (cj)
                                if (picj. contains (d [k]))
                                {
                                        // if it wasn't there before
                                //      if (!prev. contains (d [k]))
                                        // add it to the transposed matrix
                                        piT. add (d [k], cj);
                                }
                                // if this cell disapeared from pi (cj)
                                else
                                {
                                        // remove it from the transpose of pi
                                        piT. add (d [k], 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 << "Homology gvf computed in " << compTime << ".\n";

        return;
} /* homologyGVF2 */


#endif // _CHAINCON_HOMGVF2_H_