psetjoin.cpp

Go to the documentation of this file.

/// @addtogroup homology
/// @{
 
/////////////////////////////////////////////////////////////////////////////
///
/// @file psetjoin.cpp
///
/// A program for joining topologically disjoint sets of points (or cubes).
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2020 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 November 20, 2006. Last revision: August 23, 2010.
 
 
#include "chomp/system/config.h"
#include "chomp/system/textfile.h"
#include "chomp/cubes/pointset.h"
#include "chomp/struct/hashsets.h"
#include "chomp/system/timeused.h"
#include "chomp/system/arg.h"
 
#include <exception>
#include <new>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <sstream>
#include <vector>
#include <algorithm>
 
using namespace chomp::homology;
 
 
// --------------------------------------------------
// -------------------- OVERTURE --------------------
// --------------------------------------------------
 
const char *title = "\
Join PointSets, ver. 0.02. Copyright (C) 1997-2020 by Pawel Pilarczyk.\n\
This is free software. No warranty. Consult 'license.txt' for details.";
 
const char *helpinfo ="\
Call with: infile1.cub ... infileN.cub prefix.\n\
The program reads the sets of points in the given files,\n\
and joins these sets which are topologically disjoint, that is,\n\
can be separted back by taking their connected components.\n\
The resulting sets are written to prefix001.cub, prefix002.cub, etc.\n\
Switches and additional arguments:\n\
-d N - the minimal distance between the sets (default: 1),\n\
-s N - don't join the largest N input sets,\n\
-h - display this brief help information only and exit.\n\
For more information ask the author at http://www.PawelPilarczyk.com/.";
 
const int maxnames = 16000;
 
class pointsetptr
{
public:
        const coordinate *operator [] (int i) {return (*p) [i];}
 
        pointset *p;
 
}; /* class pointsetptr */
 
inline bool operator < (const pointsetptr &x, const pointsetptr &y)
{
        // ! REVERSE ORDER !
        return (x. p -> size () > y. p -> size ());
} /* operator < */
 
 
// --------------------------------------------------
// -------------------- PSETJOIN --------------------
// --------------------------------------------------
 
int psetbold (pointset &p, int howmany, bool debug = false)
{
        int_t cur = 0;
        while (howmany --)
        {
                int_t maxpoint = p. size ();
                if (debug)
                        sout << "\n" << howmany << " - " << maxpoint;
                if (debug)
                        std::cout << "\n==========\n" << p << "\n==========\n";
                while (cur < maxpoint)
                {
                        // add the neighbors of this point
                        neighbors n (p [cur ++], p. dimension ());
                        const coordinate *c;
                        while ((c = n. get ()) != 0)
                        {
                                if (debug)
                                        sout << "#";
                                if (debug)
                                {
                                        for (int i = 0; i < p. dimension (); ++ i)
                                                sout << (i ? " " : "(") << c [i];
                                        sout << ")";
                                }
                                if (p. check (c))
                                {
                                        if (debug)
                                                sout << "&";
                                        continue;
                                }
                                if (debug)
                                        sout << "+";
                                if (countneighbors (p, c, INSIDE, 1))
                                        p. add (c);
                                if (debug)
                                        sout << "=";
                        }
                }
        }
        return 0;
} /* psetbold */
 
bool disjoint (const pointset &p, const pointset &q)
{
        if (p. size () > q. size ())
                return disjoint (q, p);
        for (int_t i = 0; i < p. size (); ++ i)
        {
                if (q. check (p [i]))
                        return false;
        }
        return true;
} /* disjoint */
 
/// Joins the sets of cubes.
/// Returns: 0 = Ok, -1 = Error.
int joinsets (char *innames [], int nnames, char *prefix, int distance,
        int firstset)
{
        // make sure that the value of "firstset" is correct
        if (firstset < 0)
                firstset = 0;
        if (firstset >= nnames)
                firstset = nnames - 1;
 
        // read the sets of points
        sout << "Reading " << nnames << " input files... ";
        std::vector <pointsetptr> p (nnames);
        int count = 0;
        for (int i = 0; i < nnames; ++ i)
        {
                std::ifstream in (innames [i]);
                if (!in)
                        fileerror (innames [i]);
                p [i]. p = new pointset;
                in >> *(p [i]. p);
                count += p [i]. p -> size ();
        }
        sout << count << " points total read.\n";
 
        // sort the sets according to the number of cubes
        std::sort (p. begin (), p. end ());
 
        // create neighborhoods of the sets of points
        sout << "Enhancing the sets... ";
        pointset **bold = new pointset * [nnames];
        for (int i = firstset; i < nnames; ++ i)
        {
                bold [i] = new pointset;
                bold [i] -> add (*(p [i]. p));
                psetbold (*(bold [i]), distance);
        }
        sout << "Done.\n";
        
        // join the sets which are disjoint
        sout << "Joining the sets... ";
        int nsets = nnames;
        for (int i = firstset; i < nsets; ++ i)
        {
                for (int j = firstset; j < nsets; ++ j)
                {
                        if ((i == j) || !disjoint (*(bold [i]), *(p [j]. p)))
                                continue;
                        bold [i] -> add (*(bold [j]));
                        delete bold [j];
                        p [i]. p -> add (*(p [j]. p));
                        delete p [j]. p;
                        -- nsets;
                        for (int k = j; k < nsets; ++ k)
                        {
                                bold [k] = bold [k + 1];
                                p [k]. p = p [k + 1]. p;
                        }
                }
        }
        sout << "Done.\n";
 
        // write the components to the output file
        sout << "Writing " << nsets << " files... ";
        for (int n = 0; n < nsets; ++ n)
        {
                // prepare a file name for this set
                std::ostringstream s;
                s << prefix;
                for (int digit = 10000; digit > 1; digit /= 10)
                        if (((n + 1) < digit) && (nsets >= digit))
                                s << '0';
                s << (n + 1) << ".cub";
 
                // open a file to write the components to
                std::ofstream out (s. str (). c_str ());
                if (!out)
                        fileerror (s. str (). c_str (), "create");
                out << *(p [n]. p);
                delete (p [n]. p);
                if (n >= firstset)
                        delete bold [n];
        }
        sout << "Done.\n";
 
        return 0;
} /* join sets */
 
 
// --------------------------------------------------
// ---------------------- MAIN ----------------------
// --------------------------------------------------
 
/// The main function of the program.
/// Returns: 0 = Ok, -1 = Error, 1 = Help displayed, 2 = Wrong arguments.
int main (int argc, char *argv [])
{
        char *filenames [maxnames];
        int nnames = 0;
        int distance = 1;
        int firstset = 0;
 
        arguments a;
        arg (a, "d", distance);
        arg (a, "s", firstset);
        arg (a, NULL, filenames, nnames, maxnames);
        arghelp (a);
 
        argstreamprepare (a);
        int argresult = a. analyze (argc, argv);
        argstreamset ();
 
        // show the program's title
        if (argresult >= 0)
                sout << title << '\n';
 
        // if not enough file names, help should be displayed
        if (nnames < 2)
                argresult = 1;
 
        // 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 filenames present, show help information
        if (argresult > 0)
        {
                sout << helpinfo << '\n';
                return 1;
        }
 
        // try running the main function and catch an error message if thrown
        try
        {
                joinsets (filenames, nnames - 1, filenames [nnames - 1],
                        distance, firstset);
                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 */
 
/// @}