Commit 66b436b4 authored by Martin Marinov's avatar Martin Marinov

Merge commit '46e532f8' into marinom/merge-from-ReForm

parents 67240c66 46e532f8
Pipeline #15402 passed with stages
in 4 minutes and 56 seconds
......@@ -22,7 +22,6 @@
#include <cmath>
#include <climits>
#include <CoMISo/Utils/VSToolsT.hh>
#include <CoMISo/Utils/gmm.hh>
#include <CoMISo/Config/CoMISoDefines.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_;}
......
......@@ -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(
......
......@@ -40,9 +40,8 @@
//== INCLUDES =================================================================
#include "Eigen_Tools.hh"
#include <queue>
#include <CoMISo/Utils/VSToolsT.hh>
#include <gmm/gmm.h>
#include <queue>
//== NAMESPACES ===============================================================
......
......@@ -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 ===============================================================
......
......@@ -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);
......
......@@ -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
......
......@@ -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);
}
//! 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(double_round(_x) - _x);
}
//! get if _x is rounded within _tol
inline double rounded(const double _x, const double _tol)
{
return round_residue(_x) <= _tol;
}
//! compare two double values within _tol
inline bool 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
//=============================================================================
/*===========================================================================*\
* *
* 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
//=============================================================================
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment