Merge branch 'marinom/merge-from-ReForm' into 'master'

Clean up the math tools and round functionality

See merge request !57

NSolver/NPDerivativeChecker.hh

@@ -22,7 +22,6 @@
 #include <cmath>
 #include <climits>
 
-#include <CoMISo/Utils/VSToolsT.hh>
 #include <CoMISo/Utils/gmm.hh>
 
 #include <CoMISo/Config/CoMISoDefines.hh>

Solver/ConstrainedSolver.hh

@@ -41,13 +41,10 @@
 #include "MISolver.hh"
 #include <vector>
 
-//== FORWARDDECLARATIONS ======================================================
-
-//== DEFINES ==================================================================
-#define ROUND(x) ((x)<0?int((x)-0.5):int((x)+0.5))
 //== NAMESPACES ===============================================================
 
-namespace COMISO {
+namespace COMISO 
+{
 //== CLASS DEFINITION =========================================================
 
 /** \class ConstrainedSolver ConstrainedSolver.hh <ACG/.../ConstrainedSolver.hh>
@@ -263,37 +260,6 @@ public:
   // Get/Set whether the constraint reordering is used (default true)
   bool& use_constraint_reordering() { return miso_.use_constraint_reordering(); }
 
-/** @name Verify the result.
- * Functions to verify the result of the constrained solver. Are the constraints met, are the correct variables correctly rounded ...
- */
-/*@{*/
-
-
-  template<class RMatrixT, class CMatrixT, class VectorT>
-  double verify_constrained_system( 
-				   const RMatrixT& _conditions,
-				   const CMatrixT& _A,
-				   const VectorT&  _x,
-				   const VectorT&  _rhs,
-				   double          _eps = 1e-9);
-
-  template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT>
-  double verify_constrained_system_round( 
-					 const RMatrixT& _conditions,
-					 const CMatrixT& _A,
-					 const VectorT&  _x,
-					 const VectorT&  _rhs,
-					 const VectorIT& _idx_to_round,
-					 double          _eps = 1e-9);
-
-  template<class RMatrixT, class VectorT, class VectorIT>
-  void verify_mi_factored( const RMatrixT& _conditions,
-			   const RMatrixT& _B, 
-			   const VectorT&  _x,
-			   const VectorIT& _idx_to_round );
-/*@}*/
-
-
   /// Access the MISolver (e.g. to change settings)
   COMISO::MISolver& misolver() { return miso_;}
 

Solver/ConstrainedSolverT.cc

@@ -125,76 +125,19 @@ ConstrainedSolver::solve_const(const RMatrixT& _constraints,
 //-----------------------------------------------------------------------------
 
 
-template<class RMatrixT, class VectorT, class VectorIT >
-void
-ConstrainedSolver::solve(
-  RMatrixT& _constraints,
-  RMatrixT& _B,
-  VectorT&  _x,
-  VectorIT& _idx_to_round,
-  double    _reg_factor,
-  bool      _show_miso_settings,
-  bool      _show_timings)
+template <class RMatrixT, class VectorT, class VectorIT>
+void ConstrainedSolver::solve(RMatrixT& _constraints, RMatrixT& _B, VectorT& _x,
+    VectorIT& _idx_to_round, double _reg_factor, bool _show_miso_settings,
+    bool _show_timings)
 {
   // convert into quadratic system
   VectorT rhs;
-  gmm::col_matrix< gmm::rsvector< double > > A;
+  gmm::col_matrix<gmm::rsvector<double>> A;
   COMISO_GMM::factored_to_quadratic(_B, A, rhs);
 
   // solve
-  solve(_constraints, A, _x, rhs,
-    _idx_to_round, _reg_factor,
-    _show_miso_settings,
-    _show_timings);
-
-  //   int nrows = gmm::mat_nrows(_B);
-  //   int ncols = gmm::mat_ncols(_B);
-  //   int ncons = gmm::mat_nrows(_constraints);
-
-  //   if( _show_timings) std::cerr << __FUNCTION__ << "\n Initial dimension: " << nrows << " x " << ncols << ", number of constraints: " << ncons << std::endl;
-
-  //   // StopWatch for Timings
-  //   Base::StopWatch sw, sw2; sw.start(); sw2.start();
-
-  //   // c_elim[i] = index of variable which is eliminated in condition i
-  //   // or -1 if condition is invalid
-  //   std::vector<int> c_elim( ncons);
-
-  //   // apply sparse gauss elimination to make subsequent _constraints independent
-  //   make_constraints_independent( _constraints, _idx_to_round, c_elim);
-  //   double time_gauss = sw.stop()/1000.0; sw.start();
-
-  //   // eliminate conditions and return column matrix Bcol
-  //   gmm::col_matrix< gmm::rsvector< double > > Bcol( nrows, ncols);
-
-  //   // reindexing vector
-  //   std::vector<int>                          new_idx;
-
-  //   eliminate_constraints( _constraints, _B, _idx_to_round, c_elim, new_idx, Bcol);
-  //   double time_eliminate = sw.stop()/1000.0; 
-
-  //   if( _show_timings) std::cerr << "Eliminated dimension: " << gmm::mat_nrows(Bcol) << " x " << gmm::mat_ncols(Bcol) << std::endl;
-
-  //   // setup and solve system
-  //   double time_setup = setup_and_solve_system( Bcol, _x, _idx_to_round, _reg_factor, _show_miso_settings);
-  //   sw.start();
-
-  //   //  double time_setup_solve = sw.stop()/1000.0; sw.start();
-
-  //   // restore eliminated vars to fulfill the given conditions
-  //   restore_eliminated_vars( _constraints, _x, c_elim, new_idx);
-
-  //   double time_resubstitute = sw.stop()/1000.0; sw.start();
-
-  //   //  double time_total = sw2.stop()/1000.0;
-
-  //   if( _show_timings) std::cerr << "Timings: \n\t" <<
-  //     "Gauss Elimination  " << time_gauss          << " s\n\t" <<
-  //     "System Elimination " << time_eliminate      << " s\n\t" <<
-  //     "Setup              " << time_setup          << " s\n\t" <<
-  //    // "Setup + Mi-Solver  " << time_setup_solve    << " s\n\t" <<
-  //     "Resubstitution     " << time_resubstitute   << " s\n\t" << std::endl << std::endl;
-  //     //"Total              " << time_total          << std::endl;
+  solve(_constraints, A, _x, rhs, _idx_to_round, _reg_factor,
+      _show_miso_settings, _show_timings);
 }
 
 
@@ -1266,149 +1209,6 @@ ConstrainedSolver::restore_eliminated_vars(
 //-----------------------------------------------------------------------------
 
 
-template<class RMatrixT, class VectorT, class VectorIT >
-void
-ConstrainedSolver::verify_mi_factored(
-  const RMatrixT& _conditions,
-  const RMatrixT& _B,
-  const VectorT&  _x,
-  const VectorIT& _idx_to_round)
-{
-  DEB_enter_func;
-  DEB_out(2, "######### Verify Constrained Solver Result ############\n");
-
-  // create extended x vector
-  std::vector<double> x(_x);
-  x.resize(x.size() + 1);
-  x.back() = 1.0;
-
-  // verify conditions
-  std::vector<double> a(gmm::mat_nrows(_conditions));
-
-  gmm::mult(_conditions, x, a);
-
-  int conditions_not_ok = 0;
-  for (unsigned int i = 0; i < a.size(); ++i)
-    if (a[i] > 1e-6)
-    {
-      ++conditions_not_ok;
-    }
-
-  DEB_out_if(conditions_not_ok == 0, 2, "all conditions are ok!\n")
-    DEB_out_if(conditions_not_ok != 0, 1, " conditions are not fullfilled:\n ")
-
-    // verify rounding
-    int roundings_not_ok = 0;
-  for (unsigned int i = 0; i < _idx_to_round.size(); ++i)
-  {
-    double d = _x[_idx_to_round[i]];
-    if (fabs(d - round(d)) > 1e-6)
-      ++roundings_not_ok;
-  }
-
-  DEB_out_if(roundings_not_ok, 1, roundings_not_ok << " Integer variables are not rounded\n")
-    DEB_out_if(!roundings_not_ok, 2, "all Integer roundings are ok\n")
-
-    // evaluate energy
-    VectorT Bx(x);
-  gmm::mult(_B, x, Bx);
-  DEB_out(1, "Total energy: " << gmm::vect_sp(Bx, Bx) << "\n");
-  DEB_out(2, "######### FINISHED ############\n");
-}
-
-
-
-//-----------------------------------------------------------------------------
-
-
-template<class RMatrixT, class CMatrixT, class VectorT>
-double
-ConstrainedSolver::verify_constrained_system(
-  const RMatrixT& _conditions,
-  const CMatrixT& _A,
-  const VectorT&  _x,
-  const VectorT&  _rhs,
-  double          _eps)
-{
-  DEB_enter_func;
-  typedef typename linalg_traits<RMatrixT>::const_sub_row_type RowT;
-  typedef typename linalg_traits<RowT>::const_iterator RIter;
-
-  VectorT Ax(_x.size());
-  gmm::mult(_A, _x, Ax);
-
-  gmm::add(_rhs, gmm::scaled(Ax, -1.0), Ax);
-  double norm = gmm::vect_norm2(Ax);
-  //std::cerr << __FUNCTION__ << ": Error residual: " << norm << " vector : " << Ax << std::endl;
-
-  DEB_out(2, ": Checking constraints...\n");
-
-  unsigned int row_cond = gmm::mat_nrows(_conditions);
-  unsigned int col_cond = gmm::mat_ncols(_conditions);
-  bool all_conditions_ok = true;
-  for (unsigned int r = 0; r < row_cond; ++r)
-  {
-    double cond_value = 0.0;
-    RowT row = gmm::mat_const_row(_conditions, r);
-    RIter row_it = gmm::vect_const_begin(row);
-    RIter row_end = gmm::vect_const_end(row);
-    //std::cerr << "\t checking row : " << row << std::endl;
-
-    for (; row_it != row_end; ++row_it)
-    {
-      if (row_it.index() == col_cond - 1)
-        cond_value += (*row_it);
-      else
-        cond_value += _x[row_it.index()] * (*row_it);
-    }
-    //std::cerr << "\t Value is : " << cond_value << std::endl;
-    //std::cerr << "--- --- --- --- ---\n";
-    if (fabs(cond_value) > _eps)
-    {
-      DEB_out(1, "\t Error on row " << r << " with vector " << row
-        << " and condition value " << cond_value << "\n");
-      all_conditions_ok = false;
-    }
-  }
-  DEB_out(1,
-    (all_conditions_ok ? ": All conditions ok!" : ": Some conditions not ok!") << "\n")
-    return norm;
-}
-
-//-----------------------------------------------------------------------------
-
-
-template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT>
-double
-ConstrainedSolver::verify_constrained_system_round(
-  const RMatrixT& _conditions,
-  const CMatrixT& _A,
-  const VectorT&  _x,
-  const VectorT&  _rhs,
-  const VectorIT& _idx_to_round,
-  double          _eps)
-{
-  DEB_enter_func;
-  // test integer roundings
-  DEB_out(2, ": Testing integer roundings...\n");
-  bool all_roundings_ok = true;
-
-  for (unsigned int i = 0; i < _idx_to_round.size(); ++i)
-    if (fabs(ROUND(_x[_idx_to_round[i]]) - _x[_idx_to_round[i]]) != 0.0)
-    {
-      DEB_out(1, "\t Warning: variable " << _idx_to_round[i] << " was not rounded!"
-        << " Value is = " << _x[_idx_to_round[i]] << "\n")
-        all_roundings_ok = false;
-    }
-  DEB_out(1, (all_roundings_ok ? ": All roundings ok!" : ": Some roundings not ok!") << "\n")
-
-    // also test other stuff
-    return verify_constrained_system(_conditions, _A, _x, _rhs, _eps);
-}
-
-//-----------------------------------------------------------------------------
-
-
 template<class CMatrixT>
 void
 ConstrainedSolver::eliminate_columns(

Solver/Eigen_ToolsT.cc

@@ -40,9 +40,8 @@
 //== INCLUDES =================================================================
 
 #include "Eigen_Tools.hh"
-#include <queue>
-#include <CoMISo/Utils/VSToolsT.hh>
 #include <gmm/gmm.h>
+#include <queue>
 
 //== NAMESPACES ===============================================================
 

Solver/GMM_ToolsT.cc

@@ -37,10 +37,9 @@
 #include "GMM_Tools.hh"
 
 #define GMM_USES_LAPACK
-#include <gmm/gmm_lapack_interface.h>
 
-#include <CoMISo/Utils/VSToolsT.hh>
 #include <Base/Debug/DebOut.hh>
+#include <gmm/gmm_lapack_interface.h>
 #include <queue>
 
 //== NAMESPACES ===============================================================

Solver/MISolver.cc

@@ -47,6 +47,7 @@ ILOSTLBEGIN
 #endif
 
 #include <CoMISo/Utils/gmm.hh>
+#include <CoMISo/Utils/Tools.hh>
 
 #include <Base/Debug/DebTime.hh>
 #include <Base/Utils/StopWatch.hh>
@@ -54,8 +55,6 @@ ILOSTLBEGIN
 #include <queue>
 #include <float.h>
 
-#define ROUND(x) ((x) < 0 ? int((x) - 0.5) : int((x) + 0.5))
-
 namespace COMISO
 {
 namespace
@@ -383,7 +382,7 @@ void MISolver::solve_direct_rounding(
   Vecd elim_v;
   for (unsigned int i = 0; i < to_round.size(); ++i)
   {
-    _x[to_round[i]] = ROUND(_x[to_round[i]]);
+    _x[to_round[i]] = double_round(_x[to_round[i]]);
     elim_i.push_back(to_round[i]);
     elim_v.push_back(_x[to_round[i]]);
     // update old idx
@@ -494,7 +493,7 @@ void MISolver::solve_iterative(
         if (to_round[j] != -1)
         {
           int cur_idx = to_round[j];
-          double rnd_error = fabs(ROUND(xr[cur_idx]) - xr[cur_idx]);
+          const auto rnd_error = round_residue(xr[cur_idx]);
           if (rnd_error < r_best)
           {
             i_best = cur_idx;
@@ -507,7 +506,7 @@ void MISolver::solve_iterative(
     }
 
     // store rounded value
-    double rnd_x = ROUND(xr[i_best]);
+    const auto rnd_x = double_round(xr[i_best]);
     _x[old_idx[i_best]] = rnd_x;
 
     // compute neighbors
@@ -672,11 +671,7 @@ void MISolver::solve_multiple_rounding(
 
     // find index yielding smallest rounding error
     for (unsigned int j = 0; j < to_round.size(); ++j)
-    {
-      const auto xr_j = _x[to_round[j]];
-      const auto rnd_error = fabs(ROUND(xr_j) - xr_j);
-      rndg_queue.add(j, rnd_error);
-    }
+      rndg_queue.add(j, round_residue(_x[to_round[j]]));
     rndg_queue.get_ids(tr_best);
 
     DEB_only(time_search_next_integer += sw.stop());
@@ -693,7 +688,7 @@ void MISolver::solve_multiple_rounding(
     for (const auto tr_indx : tr_best)
     {
       const unsigned i_cur = static_cast<unsigned>(to_round[tr_indx]);
-      const double rnd_x = ROUND(_x[i_cur]); // store rounded value
+      const double rnd_x = double_round(_x[i_cur]); // store rounded value
 
       // compute neighbors
       const Col col = gmm::mat_const_col(_A, i_cur);

Solver/MISolver.hh

@@ -33,7 +33,6 @@
 #ifndef COMISO_MISOLVER_HH
 #define COMISO_MISOLVER_HH
 
-
 //== INCLUDES =================================================================
 #include <CoMISo/Config/CoMISoDefines.hh>
 #include <CoMISo/Config/config.hh>
@@ -51,19 +50,12 @@
 
 #include <vector>
 
-#define ROUND_MI(x) ((x)<0?int((x)-0.5):int((x)+0.5))
-
-
-//== FORWARDDECLARATIONS ======================================================
+//== NAMESPACES ===============================================================
 
+namespace COMISO 
+{
 
-namespace COMISO {
 class MISolverDialog;
-}
-
-//== NAMESPACES ===============================================================
-
-namespace COMISO {
 
 //== CLASS DEFINITION =========================================================
 
@@ -284,13 +276,13 @@ private:
   unsigned int noisy_;
   bool         stats_;
 
-  // time limit for gurobi solver (in seconds)
+  // time limit for Gurobi solver (in seconds)
   double       gurobi_max_time_;
 
   // flag
   bool         cholmod_step_done_;
 
-  // declar direct solver depending on availability
+  // declare direct solver depending on availability
 #if COMISO_SUITESPARSE_AVAILABLE
   COMISO::CholmodSolver   direct_solver_;
 #elif COMISO_EIGEN3_AVAILABLE

Utils/Tools.hh

@@ -26,20 +26,69 @@
 #define TOOLS_HH
 
 #include <math.h>
-#include <float.h>
+#include <stdint.h>
 
+//! get the closest integer to _x, much faster than round(), can overflow!
 inline int int_round(const double _x) 
 { 
-  return int(_x >= 0.0 ? _x + 0.5 : _x - 0.5);
+  return int(_x < 0 ? _x - 0.5 : _x + 0.5);
 }
 
-inline bool same(const double _x, const double _y, const double _tol)
+//! a popular macro to access int_round()
+#define ROUND_MI(x) int_round(x) 
+
+//! get the closest 64-it integer to _x, much faster than round(), no overflow!
+inline int64_t int64_round(const double _x)
+{
+  return int64_t(_x < 0 ? _x - 0.5 : _x + 0.5);
+}
+
+/*! 
+Same as above, but cast to a double, much faster than round()!
+Pass through an int64 cast to avoid the possibility of a 32-bit overflow.
+
+The Release assembly of double_round() with MSVC:
+
+00007FFDA3D7F60C  movaps      xmm0,xmm1  
+00007FFDA3D7F60F  comisd      xmm8,xmm1  
+00007FFDA3D7F614  jbe         07FFDA3D7F61Ch
+00007FFDA3D7F616  subsd       xmm0,xmm7  
+00007FFDA3D7F61A  jmp         07FFDA3D7F620h
+00007FFDA3D7F61C  addsd       xmm0,xmm7  
+00007FFDA3D7F620  cvttsd2si   rax,xmm0   // eax if int_round() is used
+00007FFDA3D7F625  xorps       xmm6,xmm6 
+00007FFDA3D7F628  cvtsi2sd    xmm6,rax   // eax if int_round() is used
+
+The only difference between int_round and int64_round() is eax vs rax.
+
+Measurements of double_round() in MISolver::solve_multiple_rounding() for 
+large CF data on a x64 system: 
+int64_round time: 78.85s
+int_round time:   79.52s
+*/
+inline double double_round(const double _x) { return double(int64_round(_x)); }
+
+//! get the residual after rounding
+inline double round_residue(const double _x) 
 {
-  return fabs(_x - _y) < _tol;
+  return fabs(double_round(_x) - _x);
 }
+
+//! get if _x is rounded within _tol
+inline bool is_rounded(const double _x, const double _tol)
+{
+  return round_residue(_x) <= _tol;
+}
+
+//! compare two double values within _tol
+inline bool are_same(const double _x, const double _y, const double _tol)
+{
+  return fabs(_x - _y) <= _tol;
+}
+
+//! get _a * _a
 template <typename T> inline T sqr(const T& _a) { return _a * _a; }
 
 //=============================================================================
 #endif // TOOLS_HH defined
 //=============================================================================
-

Utils/VSToolsT.hh

-/*===========================================================================*\
- *                                                                           *
- *                               CoMISo                                      *
- *      Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen      *
- *                           www.rwth-graphics.de                            *
- *                                                                           *
- *---------------------------------------------------------------------------* 
- *  This file is part of CoMISo.                                             *
- *                                                                           *
- *  CoMISo 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.                                      *
- *                                                                           *
- *  CoMISo 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 CoMISo.  If not, see <http://www.gnu.org/licenses/>.          *
- *                                                                           *
-\*===========================================================================*/ 
-
-#ifndef VSTOOLS_HH
-#define VSTOOLS_HH
-
-
-//== FORWARDDECLARATIONS ======================================================
-
-//== NAMESPACES ===============================================================
-
-//== DEFINITION =========================================================
-
-/** These functions are required for Visual Studio to work around missing 
-    functions. Basic equivalent functions for double exist in the float 
-    header but are named different. So this wrapper makes them standard compatible.
-    */
-#ifdef WIN32
-#include <math.h>
-#include <float.h>
-
-namespace std {
-
-inline int isnan(double x)
-{
-  return _isnan(x);
-} 
-
-// Which idiot defines isinf as a macro somewhere?
-#ifdef isinf 
-#undef isinf
-#endif
-
-inline int isinf(double x)
-{
-  return !_finite(x);
-} 
-
-inline int isfinite(double x)
-{
-  return _finite(x);
-}
-
-}// namespace std
-
-inline double nearbyint(double x) 
-{
-   if( x >= 0.0 )
-     return int( x + 0.5 );
-   else
-     return int( x - 0.5 );
-}
-
-inline double round ( double _value ) 
-{
-   return nearbyint(_value);
-}
-
-#endif //WIN32
-
-//=============================================================================
-#endif // VSTOOLS_HH defined
-//=============================================================================
-