ssquaressim.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// Computation of the Steenrod squares for simplicial 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: December 10, 2015.


#ifndef _CHAINCON_SSQUARESSIM_H_
#define _CHAINCON_SSQUARESSIM_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/prodcell.h"
#include "chaincon/ez_aw.h"
//#include "chaincon/ez_eml.h"
#include "chaincon/ez_shi.h"


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

/// Swaps the components of all of the product cells in the chain.
/// Adds the result to the resulting chain (which is assumed to be empty).
/// Corresponds to the operator "t" in [CombSteenSq], pg. 94, for p=2.
template <class ProdCellT, class CoefT>
inline void swapProduct (const tChain<ProdCellT,CoefT> &input,
        tChain<ProdCellT,CoefT> &result)
{
        int_t length (input. size ());
        for (int_t i = 0; i < length; ++ i)
        {
                result. add (ProdCellT (input. getCell (i). getRight (),
                        input. getCell (i). getLeft ()),
                        input. getCoef (i));
        }
        return;
} /* swapProduct */

/// Computes the diagonal of a chain.
/// This corresponds to "(x,x)" in [CombSteenSq], pg. 95, formula (3).
template <class CellT, class ProdCellT, class CoefT>
inline void computeDiagonal (const tChain<CellT,CoefT> &x,
        tChain <ProdCellT,CoefT> &result)
{
        int_t xSize = x. size ();
        for (int_t i = 0; i < xSize; ++ i)
        {
                const CellT &xCell (x. getCell (i));
                const CoefT &xCoef (x. getCoef (i));
                result. add (ProdCellT (xCell, xCell), xCoef);
        }
        return;
} /* computeDiagonal */

/// Converts a chain of cells encoded as trivial product cells
/// to the chain of the actual cells.
template <class CellT, class ProdCellT, class CoefT>
inline void flattenProduct (const tChain<ProdCellT,CoefT> &chain,
        tChain<CellT,CoefT> &result)
{
        int_t size (chain. size ());
        for (int_t i = 0; i < size; ++ i)
        {
                result. add (chain. getCell (i). getCell (),
                        chain. getCoef (i));
        }
        return;
} /* flattenProduct */

/// Converts a chain of cells encoded as trivial product cells
/// to the chain of the actual cells.
template <class Cell1T, class Cell2T, class ProdCell1T, class ProdCell2T,
        class CoefT>
inline void flattenProduct
        (const tTensor<ProdCell1T,ProdCell2T,CoefT> &tensor,
        tTensor<Cell1T,Cell2T,CoefT> &result)
{
        int_t size (tensor. size ());
        for (int_t i = 0; i < size; ++ i)
        {
                result. add (tensor. left (i). getCell (),
                        tensor. right (i). getCell (),
                        tensor. coef (i));
        }
        return;
} /* flattenProduct */


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

/// Computes the Steenrod squares, given a minimal model
/// and cohomology representants. For cells of simplicial type.
/// Please, make sure to include the appropriate header file with the
/// definition of the Alexander-Whitney diagonal for the cells in use
/// before including the header file containing this procedure.
/// Follows Equation (3) in [CombSteenSq], pg. 95.
/// Warning: This procedure is very inefficient.
/// Please, use "ssquares2" instead.
template <class SimCellT, class CoefT>
inline void ssquares1
        (const std::vector<std::vector<SimCellT> > &cohomRepresentants,
        const tLinMap<SimCellT,SimCellT,CoefT> &pi,
        const tLinMap<SimCellT,SimCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<SimCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<SimCellT,CoefT> > &squareValues)
{
        using chomp::homology::sbug;

        typedef tPair<SimCellT,int> CellSquareT;
        typedef tProdCell<SimCellT> ProdCellT;

        // go through all the cohomology levels
        int maxDim = cohomRepresentants. size ();
        for (int dim = 0; dim < maxDim; ++ dim)
        {
                // show a progress indicator
        //      sbug << "DEBUG: ====== Considering dim " << dim << "...\n";

                // consider all the cohomology representants of the given dimension
                int size = cohomRepresentants [dim]. size ();
                for (int n = 0; n < size; ++ n)
                {
                        // take a cell that represents the cohomology generator
                        const SimCellT &xCell (cohomRepresentants [dim] [n]);
                //      sbug << "DEBUG: xCell " << xCell << ".\n";

                        // show a progress indicator
                //      sbug << "DEBUG: ==== Considering " << xCell <<
                //              "... ====\n";

                        // compute (x, x), the diagonal Delta on the cycle x
                        // corresponding to this cohomology representant
                        tChain<SimCellT,CoefT> xChain (incl (xCell));
                        tChain<ProdCellT,CoefT> xx;
                        computeDiagonal (xChain, xx);
                //      sbug << "DEBUG: xChain = " << xChain << "\n";
                //      sbug << "DEBUG: Delta (xChain) = " << xx << "\n";

                        // compute the components of all the i-th Steenrod squares
                        for (int k = 0; k <= dim; ++ k)
                        {
                                // apply the composition of t with EZ_SHI
                                // for the next round
                                if (k)
                                {
                                //      sbug << "DEBUG: xx = " << xx << ".\n";
                                        tChain<ProdCellT,CoefT> shi_xx;
                                        EZ_SHI (xx, shi_xx);
                                //      sbug << "DEBUG: SHI (xx) = " << shi_xx << ".\n";
                                        tChain<ProdCellT,CoefT> empty;
                                        xx. swap (empty);
                                        swapProduct (shi_xx, xx);
                                //      sbug << "DEBUG: swapped = " << xx << ".\n";
                                }

                                // if the parity of "k" is different from the parity
                                // of "dim" then skip this combination
                                if ((k + dim) & 1)
                                        continue;

                                // calculate the corresponding i and j
                                int squareOrder = (dim - k) / 2;
                                int dim1 = (dim + k) / 2;

                                // show a progress indicator
                        //      sbug << "DEBUG: === Computing for k = " <<
                        //              k << " (i=" << squareOrder <<
                        //              ", j=" << dim1 <<
                        //              ")... ===\n";

                                // apply the AW operator
                                tTensor<ProdCellT,ProdCellT,CoefT> awdiagProd;
                                EZ_AW (xx, awdiagProd);
                                tTensor<SimCellT,SimCellT,CoefT> awdiagFull;
                                flattenProduct (awdiagProd, awdiagFull);
                        //      sbug << "DEBUG: awdiagFull = " << awdiagFull << ".\n";

                                // project the computed tensor to the model
                                tTensor<SimCellT,SimCellT,CoefT> awdiag (pi (awdiagFull));
                                int_t awdiagSize = awdiag. size ();
                        //      sbug << "DEBUG: awdiag = " << awdiag << ".\n";

                                // consider all the cohomology representants
                                // that may appear in the computed formula
                                int size1 = cohomRepresentants [dim1]. size ();
                                for (int i1 = 0; i1 < size1; ++ i1)
                                {
                                        const SimCellT &cell1 =
                                                cohomRepresentants [dim1] [i1];
                                        for (int_t j = 0; j < awdiagSize; ++ j)
                                        {
                                                if (!(awdiag. left (j) == cell1))
                                                        continue;
                                                if (!(awdiag. right (j) == cell1))
                                                        continue;
                                                CoefT coef (awdiag. coef (j));
                                                if (coef == 0)
                                                        continue;
                                                int_t n = cellsSquares. add
                                                        (CellSquareT (cell1, squareOrder));
                                        //      sbug << "DEBUG: Sq^" << squareOrder <<
                                        //              " (" << cell1 << ") += " <<
                                        //              coef << " * " << xCell << ".\n";
                                        //      sbug << "DEBUG: pi (awdiag) = " << awdiag << ".\n";
                                        //      sbug << "DEBUG: cChain = " << incl (cell1) << ".\n";
                                        //      sbug << "DEBUG: xChain = " << xChain << ".\n";
                                        //      sbug << "DEBUG: awdiagFull = " << awdiagFull << ".\n";
                                                squareValues [n]. add (xCell, coef);
                                        }
                                }
                        }
                }
        }

        return;
} /* ssquares1 */

/// Computes the k-th value of the index S in the sum ranges
/// defined in Corollary 3.2 in [CombSteenSq], pg. 97.
inline int ssquare_S (int k, const std::vector<int> &ii, int n, int m)
{
        using chomp::homology::sbug;

        int S (0);
        bool add = true;
        for (int a = k + 1; a <= n; ++ a)
        {
                if (add)
                        S += ii [a];
                else
                        S -= ii [a];
                add = !add;
        }
        if (!((k + n) & 1))
                S += (m + 1) / 2;
        else
                S -= (m + 1) / 2;
        S += k / 2;
        if (S < 0)
                throw "A negative value of S detected.";
        return S;
} /* ssquare_S */

/// Computes the faces of a simplicial cell. The faces appear
/// in the formulas in Corollary 3.2 in [CombSteenSq], pg. 97.
template <class SimCellT>
inline SimCellT ssquare_face (const SimCellT &xCell,
        const std::vector<int> &ii, int n, int m, bool right_hand_side)
{
        using chomp::homology::sbug;

        SimCellT face (xCell);

        // if i+j is even (and so is i-j)
        if (!(m & 1))
        {
                // the first cell
                if (!right_hand_side)
                {
                        for (int a = n; a >= 0; a -= 2)
                        {
                                for (int b = ((a == n) ? m : (ii [a + 1] - 1));
                                        b >= (ii [a] + 1); -- b)
                                {
                                        face = SimCellT (face, b);
                                }
                        }
                }
                // the second cell
                else
                {
                        for (int a = n; a >= 0; a -= 2)
                        {
                                for (int b = ii [a] - 1; b >=
                                        ((a == 0) ? 0 : (ii [a - 1] + 1));
                                        -- b)
                                {
                                        face = SimCellT (face, b);
                                }
                        }
                }
        }
        // if i+j is odd (and so is i-j)
        else
        {
                // the third cell
                if (!right_hand_side)
                {
                        for (int a = n; a >= 1; a -= 2)
                        {
                                for (int b = (ii [a] - 1);
                                        b >= (ii [a - 1] + 1); -- b)
                                {
                                        face = SimCellT (face, b);
                                }
                        }
                }
                // the fourth cell
                else
                {
                        for (int a = n; a >= -1; a -= 2)
                        {
                                for (int b = ((a == n) ? m : (ii [a + 1] - 1));
                                        b >= ((a == -1) ? 0 : (ii [a] + 1));
                                        -- b)
                                {
                                        face = SimCellT (face, b);
                                }
                        }
                }
        }

        return face;
} /* ssquare_face */

/// Computes a single Steenrod square.
/// Follows the formulas in Corollary 3.2 in [CombSteenSq], pg. 97.
template <class SimCellT, class CoefT>
inline CoefT ssquare (int i, int j,
        const SimCellT &cCell, const tChain<SimCellT,CoefT> &xChain,
        const tLinMap<SimCellT,SimCellT,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;

        // determine the values of the auxiliary variables m and n
        const int n = j - i;
        const int m = i + j;
//      sbug << "\nDEBUG: Computing Sq^" << i << " (c) (x):\nc = " <<
//              cChain << ",\nx = " << xChain << ".\n";

        // prepare the iteration indices
        std::vector<int> ii (n + 1);
        for (int k = n; k >= 0; -- k)
        {
                ii [k] = ssquare_S (k, ii, n, m);
        //      sbug << "DEBUG: ii [" << k << "] = " << ii [k] << ".\n";
        }

        // for all the index combinations in the big sum of sums
        while (1)
        {
        //      sbug << "DEBUG: Using ii = (";
        //      for (int a = 0; a <= n; ++ a)
        //              sbug << (a ? ", " : "") << ii [a];
        //      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 SimCellT &sCell (xChain. getCell (s));
                        CoefT xCoef (xChain. getCoef (s));
                //      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
                        SimCellT xFace1 (ssquare_face (sCell, ii,
                                n, m, false));

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

                        // apply the face operators for the second cell
                        SimCellT xFace2 (ssquare_face (sCell, ii,
                                n, m, true));

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

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

                // take the next combination of indices or quit the loop
                bool valid (true);
                for (int k = n; k >= 1; -- k)
                {
                        ++ (ii [k]);
                        if (ii [k] <= ((k == n) ? m : (ii [k + 1] - 1)))
                        {
                                ii [0] = ssquare_S (0, ii, n, m);
                                break;
                        }
                        ii [k] = ssquare_S (k, ii, n, m);
                        if (k == 1)
                                valid = false;
                }
                if (!n || !valid)
                        break;
        }

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

/// Computes the Steenrod squares, given a minimal model
/// and cohomology representants. For cells of simplicial type.
/// Please, make sure to include the appropriate header file with the
/// definition of the Alexander-Whitney diagonal for the cells in use
/// before including the header file containing this procedure.
/// Follows the formulas in Corollary 3.2 in [CombSteenSq], pg. 97.
template <class SimCellT, class CoefT>
inline void ssquares2
        (const std::vector<std::vector<SimCellT> > &cohomRepresentants,
        const tLinMap<SimCellT,SimCellT,CoefT> &pi,
        const tLinMap<SimCellT,SimCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<SimCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<SimCellT,CoefT> > &squareValues)
{
        using chomp::homology::sbug;

        typedef tPair<SimCellT,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 SimCellT &cCell (cohomRepresentants [j] [n]);
                        //      sbug << "DEBUG: cCell " << cCell << ".\n";
                        //      tChain<SimCellT,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 SimCellT &xCell (cohomRepresentants [j + i] [m]);
                                //      sbug << "DEBUG: xCell " << xCell << ".\n";
                                        tChain<SimCellT,CoefT> xChain (incl (xCell));

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

                                        // compute the Steenrod square Sq^i (c) (x)
                                        // following the formula in Corollary 3.2 of the paper
                                        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;
} /* ssquares2 */

/// Computes the Steenrod squares, given a minimal model
/// and cohomology representants. For cells of simplicial type.
/// Please, make sure to include the appropriate header file with the
/// definition of the Alexander-Whitney diagonal for the cells in use
/// before including the header file containing this procedure.
/// Follows Corollary 3.2 in [CombSteenSq], pg. 97.
template <class SimCellT, class CoefT>
inline void ssquares
        (const std::vector<std::vector<SimCellT> > &cohomRepresentants,
        const tLinMap<SimCellT,SimCellT,CoefT> &pi,
        const tLinMap<SimCellT,SimCellT,CoefT> &incl,
        chomp::homology::hashedset<tPair<SimCellT,int> > &cellsSquares,
        chomp::homology::multitable<tChain<SimCellT,CoefT> > &squareValues)
{
        return ssquares2 (cohomRepresentants, pi, incl,
                cellsSquares, squareValues);
} /* ssquares */


#endif // _CHAINCON_SSQUARESSIM_H_