---

mapper.cpp

//$Header: /home/ben/Mapper/c++/RCS/mapper.cpp,v 6.14 2002/07/09 23:03:15 ben Exp $
// Copyright Benedict Adamson 2002.
// This file is part of Mapper.

// Mapper 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 2 of the License, or
// (at your option) any later version.

// Mapper 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 Mapper; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <cassert>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <unistd.h>
#include <vector>

#include <framework/LocalFileInputSource.hpp>
#include <framework/StdInInputSource.hpp>
#include <util/PlatformUtils.hpp>
#include <util/XMLString.hpp>

#include "MapDataDocument.h"
#include "MapperErrorHandler.h"
#include "assoc.h"
#include "curve.h"
#include "minimise.h"
#include "total_energy_function.h"



namespace status {
   enum code {
      OK = 0,
      error = 1
   };
}


namespace {

   void solve(
      total_energy_function &e,
      vector<double> &p, //input, output
      const double ftol,
      unsigned max_iter,
      unsigned &iter, //output
      double &fmin //output
      )
      {
         {//preconditions:
            assert( p.size() );
            assert( 0.0 < ftol && ftol < 1 );
            assert( 0 <  max_iter );
         };
         unsigned iter1;
         if( ftol < 1E-3 && 4 <= max_iter ){
            //find an approximate solution first,
            //before solving to the fine tolerance.
            const double ftol1 = 10.0 * ftol;
            const unsigned max_iter1 = max_iter/2;
            e.set_tol( ftol1 );
            solve(
               e,
               p,
               ftol1, max_iter1,
               iter1, fmin
               );
         }else
            iter1 = 0;
         ;
         unsigned max_iter2 = max_iter2 = max_iter - iter1;
         unsigned iter2;
         e.set_tol( ftol );
         minimise::Fletcher_Reeves_Polak_Ribere(
            e,
            p,
            ftol, 2U, max_iter2,
            iter2, fmin
            );
         iter = iter1 + iter2;
         {//postconditions:
            assert( iter <= max_iter );
            //assert( fmin == f(p) || iter == max_iter); within rounding error.
         };
      }
   
}



int main(
   int argc,
   char **argv
   )
{
   const char *my_name = argv[0];
   {//initialisation
      try{ //Initialise Xerces
         XMLPlatformUtils::Initialize();
      }catch( const XMLException &toCatch ){
         cerr << my_name << ":error: during initialization\n";
         return status::error;
      };
   };
   MapperErrorHandler error_handler( my_name );
   //Initialise command-line options
   const char *input_file = 0; //0 indicates stdin.
   const char *output_file = 0; //0 indicates stdout.
   double ftol = 1E-8;
   unsigned max_iter = 0U;
   {//process the command line arguments
      int c;
      while( (c = getopt( argc, argv, "f:i:o:" )) != -1 ){
         switch( c ){
            case 'o':
               output_file = optarg;
               break;
            case 'f':{
               char *end;
               ftol = strtod( optarg, &end );
               if( *end || ftol <= 0.0 || 1.0 <= ftol ){
                  error_handler.badOptionArgument( 'f', optarg );
                  ftol = 1E-8;
               };
               break;}
            case 'i':{
               char *end;
               unsigned long i;
               i = strtoul( optarg, &end, 10 );
               if( *end ){
                  error_handler.badOptionArgument( 'i', optarg );
                  i = 0;
               };
               max_iter = i;
               break;}
            case '?':
               error_handler.error( "unrecognised command-line option" );
               break; //Give the compiler a break
#ifndef NDEBUG
            default: //unexpected
               abort();
               break; //Give the compiler a break
#endif
         };
      };
      if( optind == argc-1 ) //file given
         input_file = argv[optind];
      else if( optind == argc ) //no file given
            ; //OK
      else{ //several files given
         error_handler.error( "too many files" );
      };
      if( error_handler.errors ){//usage error(s)
         XMLPlatformUtils::Terminate();
         return status::error;
      };
   };

   MapDataDocument doc;

   {//read input file
     InputSource *source;
     if( input_file )
       source = new LocalFileInputSource(
          XMLString::transcode(input_file)
          );
     else
        source = new StdInInputSource;
     ;
     doc.read( *source, error_handler );
     delete source;
     if( error_handler.errors ){//error(s) during parsing
        XMLPlatformUtils::Terminate();
        return status::error;
     };
   };

   unsigned n_constraints;
   {//report on and check input data
      const unsigned n_points = doc.named_points().size();
      const unsigned n_curves = doc.curves().size();
      const unsigned n_fixed_points = doc.fixed_points().size();
      n_constraints = doc.constraints().size();
      assert( n_fixed_points <= n_points );
      cerr << n_points << " points, "
           << n_fixed_points << " fixed points, "
           << n_curves << " curves, "
           << n_constraints << " constraints\n";
      if( n_fixed_points < 2 )
         error_handler.warning( "fewer than 2 fixed points" );
      ;
      if( n_constraints == 0 )
         error_handler.warning( "no constraints" );
      ;
      if( n_constraints < (n_points-n_fixed_points)  )
         error_handler.warning(
            "fewer constraints than named moveable points"
            );
      ;
   };

   if( 0 < max_iter ){//calculation
      total_energy_function e( ftol );
      const set< vector2 * > ps( doc.all_points() );
      const vector< vector2 * > pv( ps.begin(), ps.end() );
      assoc( e, pv, doc.fixed_points(), doc.constraints() );
      try{//calculate point positions
         const unsigned n = e.n_vars();
         if( n && n_constraints ){
            vector<double> p(n);
            p = e.current_var(); //use any hints provided
            unsigned iter;
            double fmin;
            solve(
               e,
               p,
               ftol, max_iter,
               iter, fmin
               );
            e.move_points( p );
            cerr << iter << " iterations, "
                 << fmin << " total energy\n";
            if( max_iter <= iter )
               error_handler.warning( "iteration limit reached" );
            ;
         };
      }catch( ... ){
         error_handler.error( "during minimisation" );
      };
      if( error_handler.errors ){//error(s) during minimisation
         XMLPlatformUtils::Terminate();
         return status::error;
      };
   };

   try{//output
      if( output_file ){
         ofstream file(output_file);
         doc.write( file, ftol );
      }else
         doc.write( cout, ftol );
      ;
   }catch( ... ){
      error_handler.unexpectedError();
#ifndef NDEBUG
      abort();
#endif
   };

   XMLPlatformUtils::Terminate();
   if( error_handler.errors ){
      return status::error;
   }else
      return status::OK;
   ;
}

---