procrecur.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file procrecur.h
///
/// Post-processing of Morse decompositions by quantifying the recurrence.
/// For each Morse set, the minimum recurrence time of each element
/// is computed, unless there is no directed path in the map definition.
/// Please, note that this is experimental code at a developmental stage.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2022 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 <https://www.gnu.org/licenses/>.
 
// Started on October 8, 2008. Last revision: April 23, 2022.
 
 
#ifndef _CMGRAPHS_PROCRECUR_H_
#define _CMGRAPHS_PROCRECUR_H_
 
// include some standard C++ header files
#include <sstream>
#include <fstream>
#include <string>
#include <utility>
 
// 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/digraph.h"
#include "chomp/struct/flatmatr.h"
#include "chomp/struct/multitab.h"
 
// include local header files
#include "config.h"
#include "typedefs.h"
#include "conindex.h"
#include "morsedec.h"
#include "typedyns.h"
#include "utils.h"
 
 
namespace custom {
 
/// Post-processing of Morse decompositions by recurrence quantification.
namespace procRecurrence {
 
 
// --------------------------------------------------
// ------- post-process a Morse decomposition -------
// --------------------------------------------------
 
/// Writes a directed graph in the "dot" format.
template <class GraphType>
inline void graph2dot (const GraphType &g, std::ostream &out)
{
        int_t nVertices = g. countVertices ();
        out << "digraph G {";
 
        // output the vertices
        for (int_t vertex = 0; vertex < nVertices; ++ vertex)
        {
                out << (vertex ? " " : "") << "v" << vertex <<
                        " [label=\"" << vertex << "\"]";
        }
 
        // output the edges
        for (int_t vertex = 0; vertex < nVertices; ++ vertex)
        {
                int_t nEdges = g. countEdges (vertex);
                for (int_t edge = 0; edge < nEdges; ++ edge)
                {
                        int_t target = g. getEdge (vertex, edge);
                        out << " v" << vertex << "->v" << target;
                }
        }
 
        out << "}\n";
        return;
} /* graph2dot */
 
/// Computes minimum recurrence of vertices in the graph.
/// Fills in the array using operator [].
/// The size of the array must be equal to the number of nodes of the graph.
/// If no recurrence can be determined for some vertices (which may happen
/// if the graph is not strongly connected) then sets the value to 0.
template <class ArrayType, class GraphType>
inline void computeRecurrence (ArrayType &recurrence, const GraphType &g)
{
        using chomp::homology::sbug;
 
        // skip the analysis if the size of the set is too big
        const int maxVertices = 30000;
        if (g. countVertices () > maxVertices)
        {
                sbug << "The number of cubes exceeds " << maxVertices <<
                        ". Skipping the analysis.\n";
                return;
        }
 
        // copy the input graph to a weighted graph data structure
        // and make a note of self-edges
        sbug << "Copying the graph...\n";
        typedef chomp::homology::diGraph<int_t> WGraphType;
        WGraphType wg;
        int_t nVertices = g. countVertices ();
        std::vector<bool> selfEdges (nVertices, false);
        for (int_t vertex = 0; vertex < nVertices; ++ vertex)
        {
                wg. addVertex ();
                int_t nEdges = g. countEdges (vertex);
                for (int_t edge = 0; edge < nEdges; ++ edge)
                {
                        int_t target = g. getEdge (vertex, edge);
                        wg. addEdge (target, 1);
                        if (vertex == target)
                                selfEdges [vertex] = true;
                }
        }
 
        // prepare the shortest paths matrix:
        // D[i][j] is the shortest path from vertex i to vertex j
        sbug << "Preparing a " << nVertices << "x" << nVertices <<
                " matrix (" << (nVertices * nVertices * sizeof(int_t) /
                1024 / 1024) << " MB)...\n";
        std::vector<std::vector<signed char> > D (nVertices,
                std::vector<signed char> (nVertices));
 
        sbug << "Running Dijkstra " << nVertices << " times...\n";
        std::vector<int_t> tempResults (nVertices);
        for (int_t source = 0; source < nVertices; ++ source)
        {
                if (source && !(source % 17))
                        showProgress (source, ' ');
                //sbug << "Dijkstra " << source << "\n";
                wg. Dijkstra (source, tempResults);
                //sbug << "Done 1.\n";
                for (int_t i = 0; i < nVertices; ++ i)
                {
                        if (tempResults [i] < -128)
                                D [source] [i] = -128;
                        else if (tempResults [i] > 127)
                                D [source] [i] = 127;
                        else
                                D [source] [i] = tempResults [i];
                }
                //sbug << "Done 2.\n";
        }
 
        sbug << "Computing minimum loop lengths...\n";
        for (int_t source = 0; source < nVertices; ++ source)
        {
                // if there is an edge from source to itself
                if (selfEdges [source])
                {
                        recurrence [source] = 1;
                        continue;
                }
 
                // otherwise find a minimal loop through another vertex
                int_t minLoop (0);
                for (int_t other = 0; other < nVertices; ++ other)
                {
                        if (other == source)
                                continue;
                        // note: negative length == no path
                        int_t len1 = D [source] [other];
                        int_t len2 = D [other] [source];
                        if ((len1 <= 0) || (len2 <= 0))
                                continue;
                        if ((minLoop == 0) || (len1 + len2 < minLoop))
                                minLoop = len1 + len2;
                }
                recurrence [source] = minLoop;
        }
 
        return;
} /* computeRecurrence */
 
/// Computes Vitali variation of the values on the cubical set.
/// The current implementation is restricted to dimension 2.
template <class CubSetType, class Vector>
inline double computeVitaliVariation (const CubSetType &morseCubes,
        const Vector &recurrence, long &counter)
{
        // if the set is empty or the dimension is not supported then return 0
        if (morseCubes.empty() || (morseCubes[0].dim() != 2))
                return 0;
 
        // determine the range of coordinates
        typedef typename CubSetType::value_type cubeT;
        typedef typename cubeT::CoordType coordT;
        coordT coord [2];
        morseCubes[0].coord (coord);
        coordT xMin(coord[0]), xMax(coord[0]+1), yMin(coord[1]), yMax(coord[1]+1);
        for (int_t n = 1; n < morseCubes.size(); ++ n)
        {
                morseCubes[n].coord (coord);
                if (xMin > coord[0])
                        xMin = coord[0];
                if (xMax <= coord[0])
                        xMax = coord[0] + 1;
                if (yMin > coord[1])
                        yMin = coord[1];
                if (yMax <= coord[1])
                        yMax = coord[1] + 1;
        }
 
        // create an array for the recurrence values (-1 means undefined)
        typedef typename Vector::value_type recT;
        std::vector<recT> row (yMax - yMin, -1);
        std::vector<std::vector<recT> > R (xMax - xMin, row);
        for (int_t n = 0; n < morseCubes.size(); ++ n)
        {
                morseCubes[n].coord (coord);
                (R[coord[0]-xMin])[coord[1]-yMin] = recurrence [n];
        }
 
        // compute the variation of the recurrence in the array
        double vvar = 0;
        counter = 0;
        for (int_t x = 0; x < xMax - xMin - 1; ++ x)
        {
                for (int_t y = 0; y < yMax - yMin - 1; ++ y)
                {
                        if ((R[x])[y] <= 0)
                                continue;
                        if (((R[x+1])[y] > 0) && ((R[x])[y+1] > 0) &&
                                ((R[x+1])[y+1] > 0))
                        {
                                double v = (R[x])[y] + (R[x+1])[y+1] -
                                        (R[x+1])[y] - (R[x])[y+1];
                                vvar += (v > 0) ? v : -v;
                                ++ counter;
                        }
                }
        }
 
        return vvar;
} /* computeVitaliVariation */
 
/// Processes a single Morse set by applying the cycle decomposition method.
template <class CubSetType, class GraphType, class CubMapType>
inline void processMorseSet (const CubSetType &morseCubes,
        const CubSetType &allCubes, const GraphType &g,
        const CubMapType &theCubMap, const std::string &baseSetFileName)
{
        using namespace chomp::homology;
        using chomp::homology::timeused;
 
        // measure working time
        timeused workingTime;
        workingTime = 0;
 
        // compute the map graph restricted to the Morse set
        GraphType gm;
        computeRestrictedGraph (gm, g, allCubes, morseCubes);
 
        // compute the period of the graph
        int_t graphPeriod = computePeriod (gm);
        sbug << "Graph period of this Morse set: " << graphPeriod << ".\n";
 
        // compute the recurrence value for each vertex in the graph
        std::vector<int_t> recurrence (morseCubes.size (), 0);
        computeRecurrence (recurrence, gm);
 
        // write summary information
        int_t recMin (0), recMax (0);
        size_t countZero (0);
        for (std::vector<int_t>::const_iterator it = recurrence.begin();
                it != recurrence.end(); ++ it)
        {
                if (*it == 0)
                        ++ countZero;
                if ((*it > 0) && ((recMin == 0) || (*it < recMin)))
                        recMin = *it;
                if (*it > recMax)
                        recMax = *it;
        }
        sbug << "Recurrence between " << recMin << " and " << recMax <<
                " found.\n" << countZero <<
                " non-recurrent cubes were encountered.\n";
 
        // compute Vitali variation of the recurrence
        long counter (0);
        double vvar = computeVitaliVariation (morseCubes, recurrence, counter);
        sbug << "Vitali variation of the recurrence computed on " << counter <<
                " quadruples: " << vvar << ".\n";
 
        // save all the cubes with their recurrence times
        std::string recFileName (baseSetFileName + ".rec");
        sbug << "Saving '" << recFileName << "'...\n";
        std::ofstream f (recFileName. c_str ());
        for (size_t i = 0; i < recurrence.size (); ++ i)
        {
                f << morseCubes [i] << " " << recurrence [i] << "\n";
        }
        f. close ();
 
        // show working time
        sbug << "Recurrence processing completed in " << workingTime << ".\n";
 
        return;
} /* processMorseSet */
 
/// Post-processes a Morse decomposition by computing a decomposition
/// of each Morse set into cycle sets and analyzing this decomposition.
template <class CubSetType, class MorseDecType, class GraphType,
        class CubMapType>
inline void processMorseDec (const MorseDecType &morseDec,
        const CubSetType &allCubes, const GraphType &g,
        const CubMapType &theCubMap, const CubMapType &,
        const std::string &cubesFilePrefix,
        const std::string &procFilePrefix)
{
        using namespace chomp::homology;
 
        if (procFilePrefix.empty() && cubesFilePrefix.empty())
        {
                sbug << "WARNING: Skipping recurrence analysis, "
                        "because its results would not be saved.\n";
                return;
        }
 
        std::string prefix = (procFilePrefix.size() ?
                procFilePrefix : cubesFilePrefix) + "r";
 
        // show a message on what you are doing
        sbug << "\n=================================\n";
        sbug << "====== RECURRENCE ANALYSIS ======\n";
        sbug << "=================================\n";
        sbug << "Analyzing the Morse sets and saving results to '" <<
                prefix << "*'...\n";
        sbug << "Note: cubesFilePrefix = '" << cubesFilePrefix <<
                "', procFilePrefix = '" << procFilePrefix << "'.\n";
 
        // analyze each Morse set
        int nSets = morseDec. count ();
        for (int n = 0; n < nSets; ++ n)
        {
                // restrict the analysis to the given Morse set
                sbug << "\nMorse set no. " << n << " (" <<
                        morseDec [n]. size () << " cubes).\n";
 
                // prepare a base file name for this Morse set
                std::ostringstream s;
                s << prefix << n;
                const std::string baseSetFileName = s. str ();
 
                // process the Morse set
                processMorseSet (morseDec [n], allCubes, g,
                        theCubMap, baseSetFileName);
        }
        sbug << "============== END ==============\n\n";
 
        return;
} /* processMorseDec */
 
 
} // namespace procRecurrence
} // namespace custom
 
 
#endif // _CMGRAPHS_PROCRECUR_H_