chaincon/homgvf0.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Homology gradient vector field computation: Algorithm 0, very slow.
///
/////////////////////////////////////////////////////////////////////////////

// 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_HOMGVF0_H_
#define _CHAINCON_HOMGVF0_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 a given
/// filtered finite cell complex "K".
/// This is an old version of the algorithm, very slow.
/// 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 homologyGVF0 (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;

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

                        // 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 (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];

                // prepare a new phi
                tCombLinMap<CellT,CellT> phiAdd;

                // prepare a new projection pi
                tCombLinMap<CellT,CellT> piAdd;

                // 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 + phi bd cj + bd phi cj)
                        tCombChain<CellT> cjh;
                        cjh. add (cj);
                        cjh. add (phi (boundary (cj)));
                        cjh. add (boundary (phi (cj)));
                        tCombChain<CellT> xi (pi (cjh));

                        // if u does not appear in this chain
                        // then skip the remaining computations for cj
                        if (!xi. contains (u))
                                continue;

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

                        // compute the new image of cj by the projection
                        tCombChain<CellT> combination;
                        combination. add (phi (boundary (cj)));
                        combination. add (boundary (phi (cj)));
                        combination. add (phiAdd (boundary (cj)));
                        combination. add (boundary (phiAdd (cj)));
                        piAdd. add (cj, pi (combination));
                }

                // === DEBUG VERIFICATION BEGIN ===
                // prepare a new projection pi - to test if this is the same
                tCombLinMap<CellT,CellT> piAddTest;

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

                        // compute xi = pi (cj + phi bd cj + bd phi cj)
                        tCombChain<CellT> cjh (boundary (phi (cj)));
                        cjh. add (phi (boundary (cj)));
                        cjh. add (cj);
                        tCombChain<CellT> xi (pi (cjh));

                        // if u does not appear in this chain
                        // then skip the remaining computations for cj
                        if (!xi. contains (u))
                                continue;

                        // compute the new image of cj by the projection
                        tCombChain<CellT> combination;
                        combination. add (phi (boundary (cj)));
                        combination. add (boundary (phi (cj)));
                        combination. add (phiAdd (boundary (cj)));
                        combination. add (boundary (phiAdd (cj)));
                        piAddTest. add (cj, pi (combination));
                }

                // check if both versions of the new pi are the same
                if (!(piAdd == piAddTest))
                {
                        sout << "ERROR: The computed two versions "
                                "of the map pi are not the same!\n";
                        sout << "i = " << i << ", c_i = " << c <<
                                ", u = " << u << ".\n";
                        sout << "The wrong map piAdd:\n" << piAdd;
                        sout << "The correct map piAdd:\n" << piAddTest;
                        throw "Wrong map pi detected.\n";
                }
                // === DEBUG VERIFICATION END ===

                // replace the old phi with the new one
                phi. add (phiAdd);

                // replace the old pi with the new one
                pi. add (piAdd);

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

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

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

        return;
} /* homologyGVF0 */


#endif // _CHAINCON_HOMGVF0_H_