ssquarescub.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Computation of the Steenrod squares for cubical cells.
///
/////////////////////////////////////////////////////////////////////////////

// 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 July 4, 2013. Last revision: January 17, 2016.


#ifndef _CHAINCON_SSQUARESCUB_H_
#define _CHAINCON_SSQUARESCUB_H_


// include some standard C++ header files
#include <istream>
#include <ostream>
#include <vector>
//#include <algorithm>

// include selected header files from the CHomP library
#include "chomp/system/config.h"
#include "chomp/struct/hashsets.h"
#include "chomp/struct/multitab.h"

// include relevant local header files
#include "chaincon/chain.h"
#include "chaincon/tensor.h"
#include "chaincon/linmap.h"
#include "chaincon/pair.h"
#include "chaincon/shuffles.h"


// --------------------------------------------------
// -------------- auxiliary procedures --------------
// --------------------------------------------------

/// An auxiliary class for iterating all the possible disjoint subsets
/// A and B of integers in [0, ..., range - 1], both of the given length.
class STuples
{
public:
        /// The only constructor allowed.
        STuples (int length, int range);

        /// Gets a pair of sets A and B.
        /// Returns true if more pairs are available,
        /// or false if this was the last one (and rewinds the iterator).
        template <class ArrayA, class ArrayB>
        bool get (ArrayA &A, ArrayB &B);

private:
        /// The length of each tuple: A and B.
        int _length;

        /// The range of numbers that may appear in the tuples.
        /// The actual numbers are between 0 to _range - 1, inclusive.
        int _range;

        /// The shuffle iterator that defines the current choice of A.
        Shuffles sA;

        /// Remembers if the next choice of A is available.
        bool nextA;

        /// The current choice of A.
        std::vector<int> curA;

        /// The current choice of the complement of A.
        std::vector<int> curAcompl;

        /// The shuffle iterator that defines the current choice of B.
        Shuffles sB;

        /// Remembers if the next choice of B is available.
        bool nextB;

        /// The current choice of B.
        std::vector<int> curB;

        /// The current choice of the complement of B.
        std::vector<int> curBcompl;

}; /* class STuples */

// --------------------------------------------------

inline STuples::STuples (int length, int range):
        _length (length), _range (range), sA (length, range - length),
        curA (length), curAcompl (range - length),
        sB (length, range - 2 * length), curB (length),
        curBcompl (range - 2 * length)
{
        // define the initial choice of the sets A and B
        int junk (0);
        nextA = sA. get (curA, curAcompl, junk);
        nextB = sB. get (curB, curBcompl, junk);

        return;
} /* STuples::STuples */

template <class ArrayA, class ArrayB>
inline bool STuples::get (ArrayA &A, ArrayB &B)
{
        int length (A. size ());

        // extract the sets A and B from the current choice
        for (int i = 0; i < length; ++ i)
        {
                A [i] = curA [i];
                B [i] = curAcompl [curB [i]];
        }

        // if this was the next choice of A and B then inform about this
        if (!nextA && !nextB)
                return false;

        // prepare the next choice of the sets A and B
        int junk (0);
        nextB = sB. get (curB, curBcompl, junk);
        if (!nextB)
                nextA = sA. get (curA, curAcompl, junk);

        // inform about the existence of the next choice
        return true;
} /* STuples::get */

// --------------------------------------------------

/// Computes the sign that appears in the formula for the face.
/// Note that the indices in the code count from 0, not from 1.
/// Returns true if the sign is positive or false if it is negative.
inline bool ssquare_sign (int x, const std::vector<int> &A,
        const std::vector<int> &B)
{
        // count the number of elements in A and B smaller or equal x
        int length (A. size ());
        int countA (0);
        while ((countA < length) && (A [countA] <= x))
                ++ countA;
        int countB (0);
        while ((countB < length) && (B [countB] <= x))
                ++ countB;
        // if the number of elements in {0,...,x} \ (A u B) is even
        // then return true, otherwise return false
        return ((x + countA + countB) & 1);
} /* ssquare_sign */

/// Computes a given face of a simplicial cell.
/// See the formula in the paper for details.
template <class CubCellT>
inline CubCellT ssquare_face (const CubCellT &xCell,
        const std::vector<int> &A, const std::vector<int> &B,
        int range, bool right_hand_side)
{
        using chomp::homology::sbug;

        int length (A. size ());

        CubCellT face (xCell);
        if (right_hand_side)
        {
                for (int i = length - 1; i >= 0; -- i)
                {
                        bool sign (ssquare_sign (B [i], A, B));
                        face = CubCellT (face, sign ? B [i] :
                            (B [i] + face. dim ()));
                }
        }
        else
        {
                for (int i = length - 1; i >= 0; -- i)
                {
                        bool sign (!ssquare_sign (A [i], A, B));
                        face = CubCellT (face, sign ? A [i] :
                            (A [i] + face. dim ()));
                }
        }

        return face;
} /* ssquare_face */


// --------------------------------------------------
// ---------------- Steenrod squares ----------------
// --------------------------------------------------

/// Computes a single Steenrod square.
/// Follows the formula provided by Marek Krcal.
template <class CubCellT, class CoefT>
inline CoefT ssquare (int i, int j,
        const CubCellT &cCell, const tChain<CubCellT,CoefT> &xChain,
        const tLinMap<CubCellT,CubCellT,CoefT> &pi)
{
        using chomp::homology::sbug;

        // prepare the zero coefficient to return if necessary
        CoefT coef (0);
        if ((i < 0) || (j < 0) || (i > j))
                return coef;

        // define an iterator for the sets A and B
//      sbug << "DEBUG: Creating sTuples (" << i << ", " << (i + j) << ").\n";
        STuples sTuples (i, i + j);

        // for all the index combinations in the big sum
        bool moreAvailable (true);
        while (moreAvailable)
        {
                // get the current combination of A and B
                std::vector<int> A (i);
                std::vector<int> B (i);
                moreAvailable = sTuples. get (A, B);

                // verify the tuples if they are appropriate
                bool wrongTuple (false);
                if (true)
                {
                        for (int k = 0; k < i; ++ k)
                        {
                                if ((A [k] < i + j) && (B [k] < i + j))
                                        continue;
                                sbug << "WRONG TUPLE: k = " << k <<
                                        " and one of these: " <<
                                        A [k] << " or " << B [k] <<
                                        " is not < " << (i + j) << ".\n";
                                wrongTuple = true;
                                break;
                        }
                }

                if ((false || wrongTuple) && (sbug. show || sbug. log))
                {
                        sbug << "DEBUG: A = (";
                        for (size_t c = 0; c < A. size (); ++ c)
                                sbug << (c ? "," : "") << A [c];
                        sbug << "), B = (";
                        for (size_t c = 0; c < B. size (); ++ c)
                                sbug << (c ? "," : "") << B [c];
                        sbug << ").\n";
                }

                // for all the elements of the chain "x"
                int_t size = xChain. size ();
                for (int_t s = 0; s < size; ++ s)
                {
                        // take the cell at the given position, and the coef
                        const CubCellT &sCell (xChain. getCell (s));
                        CoefT xCoef (xChain. getCoef (s));

                        // check the dimension of the cell
                        if (sCell. dim () != (i + j))
                        {
                                sbug << "WRONG CELL DIMENSION: s = " << s <<
                                        ", sCell = " << sCell << ".\n";
                                sbug << "The dimension is " << sCell. dim () <<
                                        ", and should be " << i << ".\n";
                        }
                //      sbug << "DEBUG: s = " << s << ", sCell = " << sCell << ".\n";

                        // square the coef, because it appears twice in the product
                        xCoef *= xCoef;

                        // apply the face operators for the first cell
                        CubCellT xFace1 (ssquare_face (sCell, A, B, j, false));

                        // multiply by the corresponding coeficient
                        tChain<CubCellT,CoefT> xFace1proj (pi (xFace1));
                        int_t n1 (xFace1proj. position (cCell));
                //      sbug << "DEBUG: xFace1 = " << xFace1 << ", proj = " <<
                //              xFace1proj << ", pos = " << n1 <<
                //              ((n1 >= 0) ? "!" : ".") << "\n";
                        if (n1 < 0)
                                continue;
                        xCoef *= xFace1proj. getCoef (n1);

                        // apply the face operators for the second cell
                        CubCellT xFace2 (ssquare_face (sCell, A, B, j, true));

                        // multiply by the corresponding coeficient
                        tChain<CubCellT,CoefT> xFace2proj (pi (xFace2));
                        int_t n2 (xFace2proj. position (cCell));
                //      sbug << "DEBUG: xFace2 = " << xFace2 << ", proj = " <<
                //              xFace2proj << ", pos = " << n2 <<
                //              ((n2 >= 0) ? "!!" : ".") << "\n";
                        if (n2 < 0)
                                continue;
                        xCoef *= xFace2proj. getCoef (n2);

                        sbug << "DEBUG: xFace1 = " << xFace1 << ", proj = " <<
                                xFace1proj << ", pos = " << n1 <<
                                ((n1 >= 0) ? "!" : ".") << "\n";
                        sbug << "DEBUG: xFace2 = " << xFace2 << ", proj = " <<
                                xFace2proj << ", pos = " << n2 <<
                                ((n2 >= 0) ? "!!" : ".") << "\n";

                        // add the product to the result
                        coef += xCoef;
                }
        }

        // return the computed result
        return coef;
} /* ssquare */

/// Computes the Steenrod squares, given a minimal model
/// and cohomology representants. For cells of cubical type,
/// using the direct formula provided by Marek Krcal.
template <class CubCellT, class CoefT>
inline void ssquares
        (const std::vector<std::vector<CubCellT> > &cohomRepresentants,
        const tLinMap<CubCellT,CubCellT,CoefT> &pi,
        const tLinMap<CubCellT,CubCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<CubCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<CubCellT,CoefT> > &squareValues)
{
        using chomp::homology::sbug;

        typedef tPair<CubCellT,int> CellSquareT;

        // iterate all the valid combinations of i, c, x;
        // Sq^i maps from H^j to H^{j+i}, c is in H^j, x is in H_{j+i}
        int maxDim = cohomRepresentants. size ();
        // consider all the valid values of "j"
        for (int j = 0; j < maxDim; ++ j)
        {
                if (cohomRepresentants [j]. empty ())
                        continue;

                // consider all the valid values of "i"
                for (int i = 0; (i <= j) && (i + j < maxDim); ++ i)
                {
                        if (cohomRepresentants [i + j]. empty ())
                        {
                                continue;
                        }

                        // consider all the representants of H^j
                        int size_c = cohomRepresentants [j]. size ();
                        for (int n = 0; n < size_c; ++ n)
                        {
                                // take a representant of the cohomology generator
                                const CubCellT &cCell (cohomRepresentants [j] [n]);
                                sbug << "DEBUG: cCell " << cCell << ".\n";
                                tChain<CubCellT,CoefT> cChain (incl (cCell));

                                // consider all the representants of H^{j+i}
                                int size_x = cohomRepresentants [j + i]. size ();
                                for (int m = 0; m < size_x; ++ m)
                                {
                                        // take a representant of the cohomology generator
                                        const CubCellT &xCell (cohomRepresentants [j + i] [m]);
                                        sbug << "DEBUG: xCell " << xCell << ".\n";
                                        tChain<CubCellT,CoefT> xChain (incl (xCell));

                                        // show a progress indicator
                                        sbug << "DEBUG: i = " << i << ", j = " << j << ", j+i = " << (j + i) <<
                                                ", c = " << cCell << ", x = " << xCell << "...\n";

                                        // compute the Steenrod square Sq^i (c) (x)
                                        // following the formula provided by Marek Krcal
                                        CoefT coef (ssquare (i, j, cCell, xChain, pi));
                                        if (coef == 0)
                                                continue;
                                        int_t num = cellsSquares. add
                                                (CellSquareT (cCell, i));
                                        sbug << "DEBUG: Sq^" << i << " (" << cCell <<
                                                ") += " << coef << " * " << xCell << ".\n";
                                        squareValues [num]. add (xCell, coef);
                                }
                        }
                }
        }

        return;
} /* ssquares */

/// Computes the Steenrod squares using the only procedure in this case.
template <class CubCellT, class CoefT>
inline void ssquares1
        (const std::vector<std::vector<CubCellT> > &cohomRepresentants,
        const tLinMap<CubCellT,CubCellT,CoefT> &pi,
        const tLinMap<CubCellT,CubCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<CubCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<CubCellT,CoefT> > &squareValues)
{
        return ssquares (cohomRepresentants, pi, incl,
                cellsSquares, squareValues);
} /* ssquares1 */

/// Computes the Steenrod squares using the only procedure in this case.
template <class CubCellT, class CoefT>
inline void ssquares2
        (const std::vector<std::vector<CubCellT> > &cohomRepresentants,
        const tLinMap<CubCellT,CubCellT,CoefT> &pi,
        const tLinMap<CubCellT,CubCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<CubCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<CubCellT,CoefT> > &squareValues)
{
        return ssquares (cohomRepresentants, pi, incl,
                cellsSquares, squareValues);
} /* ssquares2 */


#endif // _CHAINCON_SSQUARESCUB_H_