Skip to content

GitLab

CoMISo
CoMISo
Commit 46e532f8 authored by Martin Marinov's avatar Martin Marinov Committed by GitHub Enterprise
1 Browse files
Options

Tidy up the round functionality (#4) 

* Clean up the round functionality
* Add int64_round() and switch double_round() to use it to work correctly for large FP values
* Remove some dead code that is using the round functionality
parent 95c47070
Hide whitespace changes
Inline Side-by-side
Showing with 67 additions and 270 deletions
Solver/ConstrainedSolver.hh
View file @ 46e532f8
......@@ -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
View file @ 46e532f8
......@@ -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/MISolver.cc
View file @ 46e532f8
......@@ -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
View file @ 46e532f8
......@@ -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
View file @ 46e532f8
......@@ -26,23 +26,67 @@
#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 ? _x - 0.5 : _x + 0.5);
}
inline double round(double _x)
//! 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 _x < 0 ? int(_x - 0.5) : int(_x + 0.5);
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; }
//=============================================================================
......
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