experiment1.cpp

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// \file
///
/// A program for experimentation with the formula for the A-W diagonal.
///
/////////////////////////////////////////////////////////////////////////////

// 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: October 21, 2013.


// the choice of cubes or simplices
#define CUBES_NOT_SIMPLICES
#ifdef CUBES_NOT_SIMPLICES
#define SPACE_WRAPPING
#endif

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

// include the string hashing definitions (before including "hashsets.h")
#include "chaincon/stringhash.h"

// 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/system/arg.h"
#include "chomp/struct/hashsets.h"

// include relevant local header files
#include "chaincon/cubcell.h"
#include "chaincon/awdiagcub.h"
#include "chaincon/simplex.h"
#include "chaincon/awdiagsim.h"
#include "chaincon/emptycell.h"
#include "chaincon/chain.h"
#include "chaincon/linmap.h"
#include "chaincon/comblinmap.h"
#include "chaincon/ringzp.h"
#include "chaincon/cellnames.h"
#include "chaincon/awdiag.h"
#include "chaincon/boundary.h"
#ifdef SPACE_WRAPPING
#include "chaincon/wrapping.h"
#endif
#include "atmodcomp.h"


// --------------------------------------------------
// -------------------- OVERTURE --------------------
// --------------------------------------------------

/// The title of the program and licensing information.
const char *title = "\
Experimentation with the formula for the A-W diagonal.\n\
Version 0.00 (Sep 24, 2015). Copyright (C) 1997-2016 by Pawel Pilarczyk.\n\
This is free software. No warranty. Consult 'license.txt' for details.";

/// Brief help information on the program's usage.
const char *helpinfo = "\
This program does some computations to experiment with the A-W diagonal.\n\
The ring of coefficients is Z or Z_p (the integers modulo prime number p).\n\
Call with:\n\
filename - the name of a file that contains a list of cubical cells,\n\
Switches and additional arguments:\n\
filename2 - the name of an additional file for relative (co)homology;\n\
-pN - selection of the coefficients: p=0 for Z, p=1 unsupported in this\n\
\tprogram, p > 1 for Z_p (default: p=2); please, use p <= 181,\n\
-r - compute reduced homology (with the empty set as a cell of dim -1),\n\
-b - don't add boundary cells (use -bx and -ba for X and A selectively),\n\
-w N1,..,Nk - set space wrapping, a.k.a. periodic boundary conditions;\n\
\trepeat the n-tuple for each dimension n; use 0 for no wrapping,\n\
-dpi, -dincl, -dphi - display the computed maps: pi, incl, phi,\n\
-drepr - display homology representants (meaningful for -aN with N < 3),\n\
-aN - homology algorithm: 0 = old (very slow), 1 = new without additional\n\
optimization (relatively fast), 2 = new (default), 3,4 = using the SNF,\n\
--verify - do additional verification of the computed maps,\n\
--log filename - save the output to a file (without progress indicators),\n\
--quiet - suppress data output to the screen (whcih can be still logged),\n\
--help - display this brief help information only and exit.\n\
For more information please consult the accompanying documentation\n\
or ask the program's author at http://www.PawelPilarczyk.com/.";


// --------------------------------------------------
// ------------------- experiment -------------------
// --------------------------------------------------

/// Runs the experimental calculations.
template <class SetT, class CellT, class CoefT, class CellNames>
inline void experiment (const SetT &K,
        const chomp::homology::hashedset<CellT> &H,
        const tLinMap<CellT,CellT,CoefT> &pi,
        const tLinMap<CellT,CellT,CoefT> &incl,
        const tLinMap<CellT,CellT,CoefT> &phi,
        CellNames &cellNames)
{
        using chomp::homology::sout;
        typedef tChain<CellT,CoefT> ChainT;
        typedef tTensor<CellT,CellT,CoefT> TensorT;

        sout << "Checking the experimental formula on " <<
                K. size () << " cells...\n";

        // go through all the cells in the cellular complex
        int_t Ksize (K. size ());
        for (int_t n = 0; n < Ksize; ++ n)
        {
                // retrieve the n-th element of the cellular complex
                const CellT &c (K [n]);

                // compute the projection of c onto the homology module
                ChainT hom_c (pi (c));
                if (hom_c. empty ())
                        continue;
        //      else if (hom_c. size () > 1)
        //      {
        //              sout << "Size of f(" << c << ") is " <<
        //                      hom_c. size ();
        //              continue;
        //      }

                // check the experimental formula for the AW diag in homology
                ChainT bd_phi_c (boundary (phi (c), K));
                ChainT neg_bd_phi_c (bd_phi_c);
                neg_bd_phi_c. negate ();
                CoefT one (1);
                ChainT one_etc;
                one_etc. add (c, one);
                one_etc += neg_bd_phi_c;
                TensorT aw1;
                AWdiagonal (one_etc, aw1);
                TensorT result1 (pi (aw1));

                // compute the result of the AW diag on the hom generators
                TensorT aw2;
                AWdiagonal (incl (hom_c), aw2);
                TensorT result2 (pi (aw2));

                // compare the results and show the information
                bool good = (result1 == result2);
                sout << "Cell " << c << ", size " << hom_c. size () <<
                        (good ? ", formula OK" : ", WRONG RESULT") << ".\n";
                if (!good)
                {
                        sout << " result 1: " << cells2names (result1,
                                cellNames, cellNames) << ".\n";
                        sout << " result 2: " << cells2names (result2,
                                cellNames, cellNames) << ".\n";
                }
        }
        return;
} /* experiment */


// --------------------------------------------------
// ---------------------- main ----------------------
// --------------------------------------------------

/// The main procedure of the program.
/// Returns: 0 = Ok, -1 = Error, 1 = Help displayed, 2 = Wrong arguments.
int main (int argc, char *argv [])
{
#ifdef CUBES_NOT_SIMPLICES
        typedef tCubCell<int_t,SettableWrapping<int_t>,
                SettableEmptyCell> CellT;
#else
        typedef tSimplex<int_t,SettableEmptyCell> CellT;
#endif
        typedef tZp<short> CoefT;
        typedef tLinMap<CellT,CellT,CoefT> MapT;

        using namespace chomp::homology;

        // turn on a message that will appear if the program does not finish
        program_time = "Aborted after";

        // prepare user-configurable data
        char *Xname = 0;
        char *Aname = 0;
        int p = 2;
        bool reduced = false;
        bool addBoundariesX = true;
        bool addBoundariesA = true;
        bool addNoBoundaries = false;
        bool displayPi = false;
        bool displayIncl = false;
        bool displayPhi = false;
        bool displayRepr = false;
        bool verify = false;
        int algorithmVersion = 2;
#ifdef SPACE_WRAPPING
        const int maxWrapping = 100;
        char *wrapping [maxWrapping];
        int nWrapping = 0;
#endif

        // analyze the command line
        arguments a;
        arg (a, NULL, Xname);
        arg (a, NULL, Aname);
        arg (a, "p", p);
        arg (a, "a", algorithmVersion);
#ifdef SPACE_WRAPPING
        arg (a, "w", wrapping, nWrapping, maxWrapping);
#endif
        argswitch (a, "r", reduced, true);
        argswitch (a, "bx", addBoundariesX, false);
        argswitch (a, "ba", addBoundariesA, false);
        argswitch (a, "b", addNoBoundaries, true);
        argswitch (a, "dpi", displayPi, true);
        argswitch (a, "dincl", displayIncl, true);
        argswitch (a, "dphi", displayPhi, true);
        argswitch (a, "drepr", displayRepr, true);
        argswitch (a, "-verify", verify, true);
        arghelp (a);

        argstreamprepare (a);
        int argresult = a. analyze (argc, argv);
        argstreamset ();

        // show the program's title
        if (argresult >= 0)
                sout << title << '\n';

        // if something was incorrect, show an additional message and exit
        if (argresult < 0)
        {
                sout << "Call with '--help' for help.\n";
                return 2;
        }

        // if help requested or no filename present, show help information
        if ((argresult > 0) || !Xname || (p < 0))
        {
                sout << helpinfo << '\n';
                return 1;
        }

        // try running the main function and catch an error message if thrown
        try
        {
                // set up the ring of coefficients and make a correction to p
                CoefT::setp (p);
                p = CoefT::getp ();

                // set the existence of the empty cell if desired
                if (reduced)
                        CellT::EmptyType::setExistence (true);

#ifdef SPACE_WRAPPING
                // set wrapping as desired
                for (int i = 0; i < nWrapping; ++ i)
                {
                        setWrapping<CellT::WrapType>
                                (std::string (wrapping [i]));
                }
#endif

                // read the filtered cell complexes
                tFilteredComplex<CellT> K;
                tFilteredComplex<CellT> L;
                readFilteredComplexes (Xname, Aname, K, L,
                        addBoundariesX && !addNoBoundaries,
                        addBoundariesA && !addNoBoundaries);
                bool relativeComplex = !L. empty ();

                // prepare the data for storing the result of computation
                chomp::homology::hashedset<CellT> H;
                MapT pi, incl, phi;
                tCellNames<CellT> cellNames;
        //      cellNames. setPrefix ("g");

                // compute an AT model
                sout << "Carrying out computations in the ring " <<
                        CoefT::ringsymbol () << ".\n";
                computeAlgTopModel (K, relativeComplex,
                        H, pi, incl, phi, cellNames,
                        displayPi, displayIncl, displayPhi,
                        displayRepr, verify, algorithmVersion);

                // verify the experimental formulas
                experiment (K, H, pi, incl, phi, cellNames);

                program_time = "Total time used:";
                program_time = 1;
                return 0;
        }
        catch (const char *msg)
        {
                sout << "ERROR: " << msg << '\n';
                return -1;
        }
        catch (const std::exception &e)
        {
                sout << "ERROR: " << e. what () << '\n';
                return -1;
        }
        catch (...)
        {
                sout << "ABORT: An unknown error occurred.\n";
                return -1;
        }

} /* main */