Commit a0e47ef4 authored by Martin Marinov's avatar Martin Marinov

Minor improvements to code whitespace style and clarity. Add some progress ticks.

parent dcc31ac2
......@@ -4,7 +4,7 @@
* 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 *
......@@ -20,7 +20,7 @@
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
\*===========================================================================*/
//=============================================================================
......@@ -48,24 +48,22 @@ namespace COMISO {
//== IMPLEMENTATION ==========================================================
template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::
solve_const( const RMatrixT& _constraints,
const CMatrixT& _A,
VectorT& _x,
const VectorT& _rhs,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
ConstrainedSolver::solve_const(const RMatrixT& _constraints,
const CMatrixT& _A,
VectorT& _x,
const VectorT& _rhs,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
// copy matrices
RMatrixT constraints( gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
RMatrixT constraints(gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
gmm::copy(_constraints, constraints);
CMatrixT A( gmm::mat_nrows(_A), gmm::mat_ncols(_A));
CMatrixT A(gmm::mat_nrows(_A), gmm::mat_ncols(_A));
gmm::copy(_A, A);
// ... and vectors
......@@ -74,13 +72,13 @@ solve_const( const RMatrixT& _constraints,
// call non-const function
solve(constraints,
A,
_x,
rhs,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
A,
_x,
rhs,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
}
......@@ -89,20 +87,19 @@ solve_const( const RMatrixT& _constraints,
template<class RMatrixT, class VectorT, class VectorIT >
void
ConstrainedSolver::
solve_const( const RMatrixT& _constraints,
const RMatrixT& _B,
VectorT& _x,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
ConstrainedSolver::solve_const(const RMatrixT& _constraints,
const RMatrixT& _B,
VectorT& _x,
const VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
// copy matrices
RMatrixT constraints( gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
RMatrixT constraints(gmm::mat_nrows(_constraints), gmm::mat_ncols(_constraints));
gmm::copy(_constraints, constraints);
RMatrixT B( gmm::mat_nrows(_B), gmm::mat_ncols(_B));
RMatrixT B(gmm::mat_nrows(_B), gmm::mat_ncols(_B));
gmm::copy(_B, B);
// ... and vectors
......@@ -110,12 +107,12 @@ solve_const( const RMatrixT& _constraints,
// call non-const function
solve(constraints,
B,
_x,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
B,
_x,
idx_to_round,
_reg_factor,
_show_miso_settings,
_show_timings);
}
......@@ -123,16 +120,15 @@ 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 )
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;
......@@ -140,59 +136,59 @@ solve(
COMISO_GMM::factored_to_quadratic(_B, A, rhs);
// solve
solve( _constraints, A, _x, rhs,
_idx_to_round, _reg_factor,
_show_miso_settings,
_show_timings);
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);
// 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();
// if( _show_timings) std::cerr << __FUNCTION__ << "\n Initial dimension: " << nrows << " x " << ncols << ", number of constraints: " << ncons << std::endl;
// // c_elim[i] = index of variable which is eliminated in condition i
// // or -1 if condition is invalid
// std::vector<int> c_elim( ncons);
// // StopWatch for Timings
// Base::StopWatch sw, sw2; sw.start(); sw2.start();
// // 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();
// // c_elim[i] = index of variable which is eliminated in condition i
// // or -1 if condition is invalid
// std::vector<int> c_elim( ncons);
// // eliminate conditions and return column matrix Bcol
// gmm::col_matrix< gmm::rsvector< double > > Bcol( nrows, ncols);
// // 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();
// // reindexing vector
// std::vector<int> new_idx;
// // eliminate conditions and return column matrix Bcol
// gmm::col_matrix< gmm::rsvector< double > > Bcol( nrows, ncols);
// eliminate_constraints( _constraints, _B, _idx_to_round, c_elim, new_idx, Bcol);
// double time_eliminate = sw.stop()/1000.0;
// // reindexing vector
// std::vector<int> new_idx;
// if( _show_timings) std::cerr << "Eliminated dimension: " << gmm::mat_nrows(Bcol) << " x " << gmm::mat_ncols(Bcol) << std::endl;
// eliminate_constraints( _constraints, _B, _idx_to_round, c_elim, new_idx, Bcol);
// double time_eliminate = sw.stop()/1000.0;
// // setup and solve system
// double time_setup = setup_and_solve_system( Bcol, _x, _idx_to_round, _reg_factor, _show_miso_settings);
// sw.start();
// if( _show_timings) std::cerr << "Eliminated dimension: " << gmm::mat_nrows(Bcol) << " x " << gmm::mat_ncols(Bcol) << std::endl;
// // 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);
// // 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_resubstitute = sw.stop()/1000.0; sw.start();
// // double time_setup_solve = sw.stop()/1000.0; sw.start();
// // double time_total = sw2.stop()/1000.0;
// // restore eliminated vars to fulfill the given conditions
// restore_eliminated_vars( _constraints, _x, c_elim, new_idx);
// 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;
// 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;
}
......@@ -200,52 +196,51 @@ solve(
template<class RMatrixT, class CMatrixT, class VectorT, class VectorIT>
void
ConstrainedSolver::
solve(
RMatrixT& _constraints,
CMatrixT& _A,
VectorT& _x,
VectorT& _rhs,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings )
void
ConstrainedSolver::solve(
RMatrixT& _constraints,
CMatrixT& _A,
VectorT& _x,
VectorT& _rhs,
VectorIT& _idx_to_round,
double _reg_factor,
bool _show_miso_settings,
bool _show_timings)
{
DEB_enter_func;
// show options dialog
if( _show_miso_settings)
if (_show_miso_settings)
miso_.show_options_dialog();
gmm::size_type nrows = gmm::mat_nrows(_A);
gmm::size_type ncols = gmm::mat_ncols(_A);
gmm::size_type ncons = gmm::mat_nrows(_constraints);
DEB_out_if( _show_timings, 1, "Initital dimension: " << nrows << " x " << ncols
<< ", number of constraints: " << ncons
<< " use reordering: " << use_constraint_reordering() << "\n")
DEB_out_if(_show_timings, 1, "Initital dimension: " << nrows << " x " << ncols
<< ", number of constraints: " << ncons
<< " use reordering: " << use_constraint_reordering() << "\n")
// StopWatch for Timings
Base::StopWatch sw, sw2; sw.start(); sw2.start();
// 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);
std::vector<int> c_elim(ncons);
// apply sparse gauss elimination to make subsequent _conditions independent
if(use_constraint_reordering())
make_constraints_independent_reordering( _constraints, _idx_to_round, c_elim);
if (use_constraint_reordering())
make_constraints_independent_reordering(_constraints, _idx_to_round, c_elim);
else
make_constraints_independent( _constraints, _idx_to_round, c_elim);
make_constraints_independent(_constraints, _idx_to_round, c_elim);
double time_gauss = sw.stop()/1000.0; sw.start();
double time_gauss = sw.stop() / 1000.0; sw.start();
// re-indexing vector
std::vector<int> new_idx;
gmm::csc_matrix< double > Acsc;
eliminate_constraints( _constraints, _A, _x, _rhs, _idx_to_round, c_elim, new_idx, Acsc);
double time_eliminate = sw.stop()/1000.0;
eliminate_constraints(_constraints, _A, _x, _rhs, _idx_to_round, c_elim, new_idx, Acsc);
double time_eliminate = sw.stop() / 1000.0;
/// TODO: temporary disable this since it was causing performance issues
//DEB_out_if( _show_timings, 2,
......@@ -253,21 +248,21 @@ solve(
// << "\n#nonzeros: " << gmm::nnz(Acsc) << "\n");
sw.start();
miso_.solve( Acsc, _x, _rhs, _idx_to_round);
double time_miso = sw.stop()/1000.0; sw.start();
miso_.solve(Acsc, _x, _rhs, _idx_to_round);
double time_miso = sw.stop() / 1000.0; sw.start();
rhs_update_table_.store(_constraints, c_elim, new_idx);
// restore eliminated vars to fulfill the given conditions
restore_eliminated_vars( _constraints, _x, c_elim, new_idx);
restore_eliminated_vars(_constraints, _x, c_elim, new_idx);
double time_resubstitute = sw.stop()/1000.0; sw.start();
double time_resubstitute = sw.stop() / 1000.0; sw.start();
double time_total = time_gauss + time_eliminate + time_miso + time_resubstitute;
DEB_out_if( _show_timings, 1,"Timings: \n\t" <<
"\tGauss Elimination " << time_gauss << " s\n\t" <<
"\tSystem Elimination " << time_eliminate << " s\n\t" <<
"\tMi-Solver " << time_miso << " s\n\t" <<
"\tResubstitution " << time_resubstitute << " s\n\t" <<
"\tTotal " << time_total << "\n\n");
DEB_out_if(_show_timings, 1, "Timings: \n\t" <<
"\tGauss Elimination " << time_gauss << " s\n\t" <<
"\tSystem Elimination " << time_eliminate << " s\n\t" <<
"\tMi-Solver " << time_miso << " s\n\t" <<
"\tResubstitution " << time_resubstitute << " s\n\t" <<
"\tTotal " << time_total << "\n\n");
}
......@@ -276,18 +271,17 @@ solve(
template<class RMatrixT, class VectorT >
void
ConstrainedSolver::
resolve(
const RMatrixT& _B,
VectorT& _x,
VectorT* _constraint_rhs,
bool _show_timings )
ConstrainedSolver::resolve(
const RMatrixT& _B,
VectorT& _x,
VectorT* _constraint_rhs,
bool _show_timings)
{
// extract rhs from quadratic system
VectorT rhs;
// gmm::col_matrix< gmm::rsvector< double > > A;
// COMISO_GMM::factored_to_quadratic(_B, A, rhs);
//TODO only compute rhs, not complete A for efficiency
// gmm::col_matrix< gmm::rsvector< double > > A;
// COMISO_GMM::factored_to_quadratic(_B, A, rhs);
//TODO only compute rhs, not complete A for efficiency
gmm::size_type m = gmm::mat_nrows(_B);
gmm::size_type n = gmm::mat_ncols(_B);
......@@ -297,20 +291,20 @@ resolve(
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
typedef typename gmm::linalg_traits<CRowT>::value_type VecValT;
gmm::resize(rhs, n-1);
gmm::resize(rhs, n - 1);
gmm::clear(rhs);
for(unsigned int i = 0; i < m; ++i)
for (unsigned int i = 0; i < m; ++i)
{
// get current condition row
CRowT row = gmm::mat_const_row( _B, i);
RIter row_it = gmm::vect_const_begin( row);
RIter row_end = gmm::vect_const_end( row);
CRowT row = gmm::mat_const_row(_B, i);
RIter row_it = gmm::vect_const_begin(row);
RIter row_end = gmm::vect_const_end(row);
if(row_end == row_it) continue;
if (row_end == row_it) continue;
--row_end;
if(row_end.index() != n-1) continue;
if (row_end.index() != n - 1) continue;
VecValT n_i = *row_end;
while(row_end != row_it)
while (row_end != row_it)
{
--row_end;
rhs[row_end.index()] -= (*row_end) * n_i;
......@@ -319,7 +313,7 @@ resolve(
// solve
resolve(_x, _constraint_rhs, &rhs,
_show_timings);
_show_timings);
}
......@@ -328,12 +322,11 @@ resolve(
template<class VectorT >
void
ConstrainedSolver::
resolve(
VectorT& _x,
VectorT* _constraint_rhs,
VectorT* _rhs ,
bool _show_timings )
ConstrainedSolver::resolve(
VectorT& _x,
VectorT* _constraint_rhs,
VectorT* _rhs,
bool _show_timings)
{
DEB_enter_func;
// StopWatch for Timings
......@@ -341,7 +334,7 @@ resolve(
sw.start();
// apply stored updates and eliminations to exchanged rhs
if(_constraint_rhs)
if (_constraint_rhs)
{
// apply linear transformation of Gaussian elimination
rhs_update_table_.cur_constraint_rhs_.resize(gmm::mat_nrows(rhs_update_table_.D_));
......@@ -349,10 +342,10 @@ resolve(
// update rhs of stored constraints
gmm::size_type nc = gmm::mat_ncols(rhs_update_table_.constraints_p_);
for(gmm::size_type i=0; i<rhs_update_table_.cur_constraint_rhs_.size(); ++i)
rhs_update_table_.constraints_p_(i,nc-1) = -rhs_update_table_.cur_constraint_rhs_[i];
for (gmm::size_type i = 0; i < rhs_update_table_.cur_constraint_rhs_.size(); ++i)
rhs_update_table_.constraints_p_(i, nc - 1) = -rhs_update_table_.cur_constraint_rhs_[i];
}
if(_rhs)
if (_rhs)
rhs_update_table_.cur_rhs_ = *_rhs;
std::vector<double> rhs_red = rhs_update_table_.cur_rhs_;
......@@ -369,17 +362,17 @@ resolve(
miso_.resolve(_x, rhs_red);
double time_miso = sw.stop()/1000.0; sw.start();
double time_miso = sw.stop() / 1000.0; sw.start();
// restore eliminated vars to fulfill the given conditions
restore_eliminated_vars( rhs_update_table_.constraints_p_, _x, rhs_update_table_.c_elim_, rhs_update_table_.new_idx_);
restore_eliminated_vars(rhs_update_table_.constraints_p_, _x, rhs_update_table_.c_elim_, rhs_update_table_.new_idx_);
double time_resubstitute = sw.stop()/1000.0; sw.start();
double time_resubstitute = sw.stop() / 1000.0; sw.start();
double time_total = time_miso + time_resubstitute;
DEB_out_if( _show_timings, 1, "Timings: \n\t" <<
"\tMi-Solver " << time_miso << " s\n\t" <<
"\tResubstitution " << time_resubstitute << " s\n\t" <<
"\tTotal " << time_total << "\n\n");
DEB_out_if(_show_timings, 1, "Timings: \n\t" <<
"\tMi-Solver " << time_miso << " s\n\t" <<
"\tResubstitution " << time_resubstitute << " s\n\t" <<
"\tTotal " << time_total << "\n\n");
}
......@@ -387,19 +380,18 @@ resolve(
template<class RMatrixT, class VectorIT >
void
ConstrainedSolver::
make_constraints_independent(
RMatrixT& _constraints,
VectorIT& _idx_to_round,
std::vector<int>& _c_elim)
void
ConstrainedSolver::make_constraints_independent(
RMatrixT& _constraints,
VectorIT& _idx_to_round,
std::vector<int>& _c_elim)
{
DEB_enter_func;
// setup linear transformation for rhs, start with identity
gmm::size_type nr = gmm::mat_nrows(_constraints);
gmm::resize(rhs_update_table_.D_, nr, nr);
gmm::clear(rhs_update_table_.D_);
for(gmm::size_type i=0; i<nr; ++i) rhs_update_table_.D_(i,i) = 1.0;
for (gmm::size_type i = 0; i < nr; ++i) rhs_update_table_.D_(i, i) = 1.0;
// Base::StopWatch sw;
// number of variables
......@@ -407,11 +399,11 @@ make_constraints_independent(
// TODO Check: HZ added 14.08.09
_c_elim.clear();
_c_elim.resize( gmm::mat_nrows(_constraints), -1);
_c_elim.resize(gmm::mat_nrows(_constraints), -1);
// build round map
std::vector<bool> roundmap( n_vars, false);
for(unsigned int i=0; i<_idx_to_round.size(); ++i)
std::vector<bool> roundmap(n_vars, false);
for (unsigned int i = 0; i < _idx_to_round.size(); ++i)
roundmap[_idx_to_round[i]] = true;
// copy constraints into column matrix (for faster update via iterators)
......@@ -422,7 +414,7 @@ make_constraints_independent(
gmm::copy(_constraints, constraints_c);
// for all conditions
for(unsigned int i=0; i<gmm::mat_nrows(_constraints); ++i)
for (unsigned int i = 0; i < gmm::mat_nrows(_constraints); ++i)
{
// get elimination variable
int elim_j = -1;
......@@ -438,28 +430,28 @@ make_constraints_independent(
typedef typename gmm::linalg_traits<CRowT>::const_iterator RIter;
// get current condition row
CRowT row = gmm::mat_const_row( _constraints, i);
RIter row_it = gmm::vect_const_begin( row);
RIter row_end = gmm::vect_const_end( row);
CRowT row = gmm::mat_const_row(_constraints, i);
RIter row_it = gmm::vect_const_begin(row);
RIter row_end = gmm::vect_const_end(row);
double elim_val = FLT_MAX;
double max_elim_val = -FLT_MAX;
// new: gcd
std::vector<int> v_gcd;
v_gcd.resize(gmm::nnz(row),-1);
v_gcd.resize(gmm::nnz(row), -1);
int n_ints(0);
bool gcd_update_valid(true);
for(; row_it != row_end; ++row_it)
for (; row_it != row_end; ++row_it)
{
int cur_j = static_cast<int>(row_it.index());
// do not use the constant part
if (cur_j != (int)n_vars - 1 )
if (cur_j != (int)n_vars - 1)
{
// found real valued var? -> finished (UPDATE: no not any more, find biggest real value to avoid x/1e-13)
if( !roundmap[ cur_j ])
if (!roundmap[cur_j])
{
if( fabs(*row_it) > max_elim_val)
if (fabs(*row_it) > max_elim_val)
{
elim_j = (int)cur_j;
max_elim_val = fabs(*row_it);
......@@ -477,7 +469,7 @@ make_constraints_independent(
// the warning below to DEB_line at high verbosity.
if ((double(int(cur_row_val)) - cur_row_val) != 0.0)
{
DEB_line(11, "coefficient of integer variable is NOT integer : " <<
DEB_line(11, "coefficient of integer variable is NOT integer : " <<
cur_row_val);
gcd_update_valid = false;
}
......@@ -486,63 +478,55 @@ make_constraints_independent(
++n_ints;
// store integer closest to 1, must be greater than epsilon_
if( fabs(cur_row_val-1.0) < elim_val && cur_row_val > epsilon_)
{
elim_int_j = cur_j;
elim_val = fabs(cur_row_val-1.0);
if (fabs(cur_row_val - 1.0) < elim_val && cur_row_val > epsilon_)
{
elim_int_j = cur_j;
elim_val = fabs(cur_row_val - 1.0);
}
}
}
}
// first try to eliminate a valid (>epsilon_) real valued variable (safer)
if( max_elim_val > epsilon_)
{}
else // use the best found integer
elim_j = elim_int_j;
if (max_elim_val <= epsilon_)
elim_j = elim_int_j; // use the best found integer
// if no integer or real valued variable greater than epsilon_ existed, then
// elim_j is now -1 and this row is not considered as a valid constraint
// store result
_c_elim[i] = elim_j;
// error check result
if( elim_j == -1)
{
// redundant or incompatible?
if( noisy_ > 0)
DEB_warning_if( (fabs(gmm::mat_const_row(_constraints, i)[n_vars-1]) > epsilon_ ), 1,
"incompatible condition: " << fabs(gmm::mat_const_row(_constraints, i)[n_vars-1]) )
// else
// std::cerr << "Warning: redundant condition:\n";
if (elim_j == -1)
{// redundant or incompatible?
DEB_warning_if((noisy_ > 0) &&
(fabs(gmm::mat_const_row(_constraints, i)[n_vars - 1]) > epsilon_), 1,
"incompatible condition: " <<
fabs(gmm::mat_const_row(_constraints, i)[n_vars - 1]));
}
else
if(roundmap[elim_j] && elim_val > 1e-6)
else if (roundmap[elim_j] && elim_val > 1e-6)
{
if (do_gcd_ && gcd_update_valid)
{
if( do_gcd_ && gcd_update_valid)
{
// perform gcd update
bool gcd_ok = update_constraint_gcd( _constraints, i, elim_j, v_gcd, n_ints);
DEB_warning_if( (noisy_ > 0) && !gcd_ok, 1," GCD update failed! "
<< DEB_os_str(gmm::mat_const_row(_constraints, i)) )
}