snfchomp.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// An interface to the SNF computation using the CHomP software package.
/// If you have a more efficient software library that can be used
/// to compute the SNF then you are welcome to use this interface
/// as a template for another header file in which you prepare
/// an interface to use the better algorithm. Please, leave this file intact,
/// create another one, and substitute the header file inclusion in "snf.h".
///
/////////////////////////////////////////////////////////////////////////////

// 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 March 24, 2009. Last revision: December 28, 2015.


#ifndef _CHAINCON_SNFCHOMP_H_
#define _CHAINCON_SNFCHOMP_H_


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

// 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/hashsets.h"
#include "chomp/struct/multitab.h"
#include "chomp/struct/integer.h"
#include "chomp/homology/chains.h"

// include relevant local header files
#include "chaincon/chain.h"
#include "chaincon/algcell.h"
#include "chaincon/euclwrap.h"


// --------------------------------------------------
// ------------- SNF computation class --------------
// --------------------------------------------------

/// An interface for the computation of the Smith Normal Form
/// of a series of a few rectangular matrices with coefficients
/// in a Euclidean Domain. A Euclidean Domain is a commutative ring
/// in which division with remainder is well defined; for example,
/// the ring of integers, the ring of polynomials in one variable, or some
/// field (like the rationals or the integers modulo a prime number p).
/// In fact, what is to be computed, are the changes of bases in the complex
/// such that each matrix in the series is in the SNF.
/// The columns and rows of each matrix in the set
/// are indexed by integer numbers starting at 0.
/// In addition to the SNF, the change-of-basis matrices A and G
/// are also computed, such that A is a change of basis from the original one
/// to the new one, and G is a change of basis from the new basis
/// back to the original one. In other words, if {M_q} is the original
/// series of matrices and {S_q} is a series of its SNF
/// then M_q = G_{q-1} S_q A_q, and G_q A_q = A_q G_q = identity.
/// Please, note that it is not possible to make all the boundary matrices
/// in the SNF with the non-zero entries starting at the 'diagonal',
/// because the boundary of a boundary must be zero.
/// Therefore, what is computed here, is a 'simple form' of the boundary matrices
/// in which every column and every row has at most one nonzero entry.
/// Please, note that the division condition required in the SNF is not satisfied
/// either. As a consequence, the collection of coefficients that define
/// a homology group may not be the same for isomorphic groups.
template <class CoefT>
class tMatrixSNF
{
public:
        /// The type of coefficients in the chain.
        typedef CoefT CoefType;

        /// The type of an abstract algebraic cell that is used
        /// to carry the numbers of columns or rows in the matrix.
        /// Please note that only the identifier of this cell is used,
        /// not its dimension, nor its boundary.
        typedef tAlgCell<int_t,unsigned char,CoefT> CellType;

        /// The type of the corresponding image chain.
        typedef tChain<CellType,CoefT> ChainType;

        /// The default constructor of a matrix.
        tMatrixSNF ();

        /// The destructor of a matrix.
        ~tMatrixSNF ();

        /// Sets the number of rows and columns in a given sub-matrix.
        void setSize (int q, int_t numRows, int_t numColumns);

        /// Adds an element to the matrix whose SNF will be computed.
        void add (int q, int_t row, int_t column, const CoefT &element);

        /// Returns the number of rows of the matrix.
        int_t getNumRows (int q) const;

        /// Returns the number of columns of the matrix.
        int_t getNumCols (int q) const;

        /// Computes the almost Smith Normal Form of all the matrices
        /// whose number is provided (indexing stats from 0 and ends at
        /// the given mumber minus 1).
        /// After all the matrices have been defined,
        /// this procedure must be called once to compute the SNF
        /// and all the auxiliary data,
        /// including the change of basis matrices A and G.
        /// Then the results of the computations can be read
        /// with the other methods available in this class.
        void computeSNF (int numMatr);

        /// Returns the coefficient in the given column of the given matrix
        /// in the SNF.
        CoefT getColCoef (int q, int_t n) const;

        /// Returns the cell number in the given column of the given matrix
        /// in the SNF, or -1 if the column is zero.
        int_t getColCell (int q, int_t n) const;

        /// Returns the coefficient in the given row of the given matrix
        /// in the SNF.
        CoefT getRowCoef (int q, int_t n) const;

        /// Returns the cell number in the given row of the given matrix
        /// in the SNF, or -1 if the row is zero.
        int_t getRowCell (int q, int_t n) const;

        /// For a cell in the original basis, returns the
        /// corresponding combination of elements in the new basis
        /// (the column of the change-of-basis matrix in this direction).
        ChainType getNewChain (int q, int_t n) const;

        /// For a cell in the original basis, returns the given row
        /// of the change-of-basis matrix towards the new basis.
        ChainType getNewChainRow (int q, int_t n) const;

        /// For a cell in the new basis, returns a combination
        /// of the original cells that correspond to this cell
        /// (the column of the change-of-basis matrix in this direction).
        ChainType getOldChain (int q, int_t n) const;

        /// For a cell in the new basis, returns the given row
        /// of the change-of-basis matrix towards the original basis.
        ChainType getOldChainRow (int q, int_t n) const;

private:
        /// Copy constructor is not allowed.
        tMatrixSNF (const tMatrixSNF<CoefT> &other);

        /// Assignment operator is not allowed.
        tMatrixSNF<CoefT> &operator = (const tMatrixSNF<CoefT> &other);

        /// The matrices whose Smith Normal Form will be computed.
        chomp::homology::multitable
                <chomp::homology::mmatrix<tEuclWrap<CoefT> > > matrices;

        /// Change of basis matrices from the original one to the SNF one.
        chomp::homology::multitable
                <chomp::homology::mmatrix<tEuclWrap<CoefT> > > chgBasisA;

        /// Change of basis matrices from the SNF one to the original one.
        chomp::homology::multitable
                <chomp::homology::mmatrix<tEuclWrap<CoefT> > > chgBasisG;

        /// Coverts a column of a matrix to an abstract algebraic chain.
        static ChainType abstractChain
                (const chomp::homology::chain<tEuclWrap<CoefT> > &ch);

        /// Remembers if the SNF has already been computed.
        bool computedSNF;

}; /* class tMatrixSNF */

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

template <class CoefT>
inline tMatrixSNF<CoefT>::tMatrixSNF (): computedSNF (false)
{
        return;
} /* tMatrixSNF::tMatrixSNF */

template <class CoefT>
inline tMatrixSNF<CoefT>::~tMatrixSNF ()
{
        return;
} /* tMatrixSNF::~tMatrixSNF */

template <class CoefT>
inline tMatrixSNF<CoefT>::tMatrixSNF (const tMatrixSNF<CoefT> &other)
{
        throw "Trying to use the undefined copy constructor for tMatrixSNF.";
        return;
} /* tMatrixSNF::tMatrixSNF */

template <class CoefT>
inline tMatrixSNF<CoefT> &tMatrixSNF<CoefT>::operator =
        (const tMatrixSNF<CoefT> &other)
{
        throw "Trying to use the undefined operator = for tMatrixSNF.";
        return *this;
} /* tMatrixSNF::operator = */

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

template <class CoefT>
inline void tMatrixSNF<CoefT>::setSize (int q,
        int_t numRows, int_t numColumns)
{
        // define the size of the matrix
        matrices [q]. define (numRows, numColumns);

        return;
} /* tMatrixSNF::setSize */

template <class CoefT>
inline void tMatrixSNF<CoefT>::add (int q, int_t row, int_t column,
        const CoefT &element)
{
        matrices [q]. add (row, column, tEuclWrap<CoefT> (element));
        return;
} /* tMatrixSNF::add */

template <class CoefT>
inline int_t tMatrixSNF<CoefT>::getNumRows (int q) const
{
        return matrices [q]. getnrows ();
} /* tMatrixSNF::getNumRows */

template <class CoefT>
inline int_t tMatrixSNF<CoefT>::getNumCols (int q) const
{
        return matrices [q]. getncols ();
} /* tMatrixSNF::getNumCols */

template <class CoefT>
inline void tMatrixSNF<CoefT>::computeSNF (int numMatr)
{
        using chomp::homology::sout;
        using chomp::homology::sbug;

        if (computedSNF)
                throw "Trying to compute the SNF twice.\n";

        for (int q = 0; q < numMatr; ++ q)
        {
                // determine the size of the matrix
        //      int_t nRows = matrices [q]. getnrows ();
                int_t nCols = matrices [q]. getncols ();

                // create the identity matrices for tracing changes of bases
                chgBasisA [q]. identity (nCols);
                chgBasisG [q]. identity (nCols);

                // connect the change of basis matrices to the main matrix
                matrices [q]. dom_img. add (chgBasisA [q]);
                matrices [q]. dom_dom. add (chgBasisG [q]);
                if (q > 0)
                {
                        matrices [q]. img_img. add (chgBasisA [q - 1]);
                        matrices [q]. img_dom. add (chgBasisG [q - 1]);
                }

                // connect the matrices between each other
                if (q < numMatr - 1)
                        matrices [q]. dom_img. add (matrices [q + 1]);
                if (q > 0)
                        matrices [q]. img_dom. add (matrices [q - 1]);
        }

//      for (int q = 0; q < numMatr; ++ q)
        for (int q = numMatr - 1; q >= 0; -- q)
        {
                if (matrices [q]. getnrows () == 0)
                        continue;
                if (matrices [q]. getncols () == 0)
                        continue;
        //      sout << "DEBUG: Original matrix:\n" << matrices [q];
                sout << "SNF " << q << ": ";
                matrices [q]. simple_form ();
                sout << "\n";
        //      matrices [q]. simple_form_to_SNF ();
        //      sbug << "DEBUG: nearly SNF:\n" << matrices [q];
        //      sbug << "DEBUG: chgBasisA:\n" << chgBasisA [q];
        //      sbug << "DEBUG: chgBasisG:\n" << chgBasisG [q];
        }

        // remember that the SNF has been computed
        computedSNF = true;

        return;
} /* tMatrixSNF::computeSNF */

template <class CoefT>
inline CoefT tMatrixSNF<CoefT>::getColCoef (int q, int_t n) const
{
        const chomp::homology::chain<tEuclWrap<CoefT> > &column =
                matrices [q]. getcol (n);
        if (column. size () < 1)
                return CoefT (0);
        else
                return tEuclWrap<CoefT> (column. coef (0)). getCoef ();
} /* tMatrixSNF::getColCoef */

template <class CoefT>
inline int_t tMatrixSNF<CoefT>::getColCell (int q, int_t n) const
{
        const chomp::homology::chain<tEuclWrap<CoefT> > &column =
                matrices [q]. getcol (n);
        if (column. size () < 1)
                return -1;
        else
                return column. num (0);
} /* tMatrixSNF::getColCell */

template <class CoefT>
inline CoefT tMatrixSNF<CoefT>::getRowCoef (int q, int_t n) const
{
        const chomp::homology::chain<tEuclWrap<CoefT> > &row =
                matrices [q]. getrow (n);
        if (row. size () < 1)
                return CoefT (0);
        else
                return tEuclWrap<CoefT> (row. coef (0)). getCoef ();
} /* tMatrixSNF::getRowCoef */

template <class CoefT>
inline int_t tMatrixSNF<CoefT>::getRowCell (int q, int_t n) const
{
        const chomp::homology::chain<tEuclWrap<CoefT> > &row =
                matrices [q]. getrow (n);
        if (row. size () < 1)
                return -1;
        else
                return row. num (0);
} /* tMatrixSNF::getRowCell */

template <class CoefT>
inline typename tMatrixSNF<CoefT>::ChainType tMatrixSNF<CoefT>::abstractChain
        (const chomp::homology::chain<tEuclWrap<CoefT> > &ch)
{
        ChainType chResult;
        int_t size = ch. size ();
        for (int_t i = 0; i < size; ++ i)
        {
                chResult. add (CellType (ch. num (i)),
                        tEuclWrap<CoefT> (ch. coef (i)). getCoef ());
        }
        return chResult;
} /* tMatrixSNF::abstractChain */

template <class CoefT>
inline typename tMatrixSNF<CoefT>::ChainType
        tMatrixSNF<CoefT>::getNewChain (int q, int_t n) const
{
        return abstractChain (chgBasisA [q]. getcol (n));
} /* tMatrixSNF::getNewChain */

template <class CoefT>
inline typename tMatrixSNF<CoefT>::ChainType
        tMatrixSNF<CoefT>::getNewChainRow (int q, int_t n) const
{
        return abstractChain (chgBasisA [q]. getrow (n));
} /* tMatrixSNF::getNewChainRow */

template <class CoefT>
inline typename tMatrixSNF<CoefT>::ChainType
        tMatrixSNF<CoefT>::getOldChain (int q, int_t n) const
{
        return abstractChain (chgBasisG [q]. getcol (n));
} /* tMatrixSNF::getOldChain */

template <class CoefT>
inline typename tMatrixSNF<CoefT>::ChainType
        tMatrixSNF<CoefT>::getOldChainRow (int q, int_t n) const
{
        return abstractChain (chgBasisG [q]. getrow (n));
} /* tMatrixSNF::getOldChainRow */


#endif // _CHAINCON_SNFCHOMP_H_