Commit cd35d6f6 authored by Henrik Zimmer's avatar Henrik Zimmer

Added noisy-parameter to set amout of output

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@27 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 91eff718
...@@ -62,7 +62,7 @@ public: ...@@ -62,7 +62,7 @@ public:
/// default Constructor /// default Constructor
ConstrainedSolver() { epsilon_ = 1e-6; } ConstrainedSolver() { epsilon_ = 1e-6; noisy_ = 1; }
/// Destructor /// Destructor
~ConstrainedSolver() { } ~ConstrainedSolver() { }
...@@ -206,6 +206,9 @@ public: ...@@ -206,6 +206,9 @@ public:
/// Set numerical epsilon for valid constraint coefficient /// Set numerical epsilon for valid constraint coefficient
void set_epsilon( double _epsilon) { epsilon_ = _epsilon;} void set_epsilon( double _epsilon) { epsilon_ = _epsilon;}
/// Set noise-level (how much std output is given) 0 basically none, 1 important stuff (warning/timing, is default), 2+ not so important
void set_noisy( int _noisy) { noisy_ = _noisy;}
/** @name Verify the result. /** @name Verify the result.
* Functions to verify the result of the constrained solver. Are the constraints met, are the correct variables correctly rounded ... * Functions to verify the result of the constrained solver. Are the constraints met, are the correct variables correctly rounded ...
*/ */
...@@ -272,6 +275,7 @@ private: ...@@ -272,6 +275,7 @@ private:
ConstrainedSolver& operator=(const ConstrainedSolver& _rhs); ConstrainedSolver& operator=(const ConstrainedSolver& _rhs);
double epsilon_; double epsilon_;
int noisy_;
}; };
......
...@@ -268,18 +268,21 @@ make_constraints_independent( ...@@ -268,18 +268,21 @@ make_constraints_independent(
// store result // store result
_c_elim[i] = elim_j; _c_elim[i] = elim_j;
// error check result // error check result
if( elim_j == -1) if( noisy_ > 0)
{ {
// redundant or incompatible? if( elim_j == -1)
if( fabs(gmm::mat_const_row(_constraints, i)[n_vars-1]) > epsilon_ ) {
std::cerr << "Warning: incompatible condition:\n"; // redundant or incompatible?
if( fabs(gmm::mat_const_row(_constraints, i)[n_vars-1]) > epsilon_ )
std::cerr << "Warning: incompatible condition:\n";
else
std::cerr << "Warning: redundant condition:\n";
}
else else
std::cerr << "Warning: redundant condition:\n"; if(roundmap[elim_j] && elim_val > 1e-6)
std::cerr << "Warning: eliminate non +-1 integer -> correct rounding cannot be guaranteed:\n"
<< gmm::mat_const_row(_constraints, i) << std::endl;
} }
else
if(roundmap[elim_j] && elim_val > 1e-6)
std::cerr << "Warning: eliminate non +-1 integer -> correct rounding cannot be guaranteed:\n"
<< gmm::mat_const_row(_constraints, i) << std::endl;
// is this condition dependent? // is this condition dependent?
if( elim_j != -1 ) if( elim_j != -1 )
...@@ -322,7 +325,6 @@ eliminate_constraints( ...@@ -322,7 +325,6 @@ eliminate_constraints(
std::vector<int>& _new_idx, std::vector<int>& _new_idx,
gmm::col_matrix<SVector3T>& _Bcol) gmm::col_matrix<SVector3T>& _Bcol)
{ {
std::cerr << __FUNCTION__ << std::endl;
// copy into column matrix // copy into column matrix
gmm::resize(_Bcol, gmm::mat_nrows(_B), gmm::mat_ncols(_B)); gmm::resize(_Bcol, gmm::mat_nrows(_B), gmm::mat_ncols(_B));
gmm::copy( _B, _Bcol); gmm::copy( _B, _Bcol);
...@@ -384,9 +386,13 @@ eliminate_constraints( ...@@ -384,9 +386,13 @@ eliminate_constraints(
std::sort(_idx_to_round.begin(), _idx_to_round.end()); std::sort(_idx_to_round.begin(), _idx_to_round.end());
_idx_to_round.resize( std::unique(_idx_to_round.begin(), _idx_to_round.end()) -_idx_to_round.begin()); _idx_to_round.resize( std::unique(_idx_to_round.begin(), _idx_to_round.end()) -_idx_to_round.begin());
std::cerr << "remaining variables: " << gmm::mat_ncols(_Bcol) << std::endl;
std::cerr << "remaining integer variables: " << _idx_to_round.size() << std::endl; if( noisy_ > 2 )
std::cerr << std::endl; {
std::cerr << __FUNCTION__ << "remaining variables: " << gmm::mat_ncols(_Bcol) << std::endl;
std::cerr << __FUNCTION__ << "remaining integer variables: " << _idx_to_round.size() << std::endl;
std::cerr << __FUNCTION__ << std::endl;
}
} }
...@@ -406,8 +412,6 @@ eliminate_constraints( ...@@ -406,8 +412,6 @@ eliminate_constraints(
std::vector<int>& _new_idx, std::vector<int>& _new_idx,
CSCMatrixT& _Acsc) CSCMatrixT& _Acsc)
{ {
std::cerr << __FUNCTION__ << std::endl;
ACG::StopWatch sw; ACG::StopWatch sw;
sw.start(); sw.start();
// define iterator on matrix A and on constraints C // define iterator on matrix A and on constraints C
...@@ -486,7 +490,8 @@ eliminate_constraints( ...@@ -486,7 +490,8 @@ eliminate_constraints(
constraint[ constraint.size()-1] = constraint_rhs; constraint[ constraint.size()-1] = constraint_rhs;
} }
} }
std::cerr << "Constraints integrated " << sw.stop()/1000.0 << std::endl; if( noisy_ > 2)
std::cerr << __FUNCTION__ << " Constraints integrated " << sw.stop()/1000.0 << std::endl;
// eliminate vars // eliminate vars
_Acsc.init_with_good_format(_A); _Acsc.init_with_good_format(_A);
...@@ -494,7 +499,8 @@ eliminate_constraints( ...@@ -494,7 +499,8 @@ eliminate_constraints(
std::vector< double > elim_varvals( elim_varids.size(), 0); std::vector< double > elim_varvals( elim_varids.size(), 0);
gmm::eliminate_csc_vars2( elim_varids, elim_varvals, _Acsc, _x, _rhs); gmm::eliminate_csc_vars2( elim_varids, elim_varvals, _Acsc, _x, _rhs);
std::cerr << "Constraints eliminated " << sw.stop()/1000.0 << std::endl; if( noisy_ > 2)
std::cerr << __FUNCTION__ << " Constraints eliminated " << sw.stop()/1000.0 << std::endl;
sw.start(); sw.start();
// init _new_idx vector // init _new_idx vector
_new_idx.resize( mat_ncols(_constraints)); _new_idx.resize( mat_ncols(_constraints));
...@@ -520,7 +526,9 @@ eliminate_constraints( ...@@ -520,7 +526,9 @@ eliminate_constraints(
std::sort(_idx_to_round.begin(), _idx_to_round.end()); std::sort(_idx_to_round.begin(), _idx_to_round.end());
_idx_to_round.resize( std::unique(_idx_to_round.begin(), _idx_to_round.end()) -_idx_to_round.begin()); _idx_to_round.resize( std::unique(_idx_to_round.begin(), _idx_to_round.end()) -_idx_to_round.begin());
std::cerr << "Indices reindexed " << sw.stop()/1000.0 << std::endl << std::endl;
if( noisy_ > 2)
std::cerr << __FUNCTION__ << "Indices reindexed " << sw.stop()/1000.0 << std::endl << std::endl;
} }
...@@ -579,7 +587,6 @@ setup_and_solve_system( CMatrixT& _B, ...@@ -579,7 +587,6 @@ setup_and_solve_system( CMatrixT& _B,
double _reg_factor, double _reg_factor,
bool _show_miso_settings) bool _show_miso_settings)
{ {
std::cerr << __FUNCTION__ << std::endl;
ACG::StopWatch s1; ACG::StopWatch s1;
ACG::StopWatch sw; sw.start(); ACG::StopWatch sw; sw.start();
unsigned int m = gmm::mat_nrows(_B); unsigned int m = gmm::mat_nrows(_B);
...@@ -590,30 +597,36 @@ setup_and_solve_system( CMatrixT& _B, ...@@ -590,30 +597,36 @@ setup_and_solve_system( CMatrixT& _B,
CMatrixT Bt; CMatrixT Bt;
gmm::resize( Bt, n, m); gmm::resize( Bt, n, m);
gmm::copy( gmm::transposed( _B), Bt); gmm::copy( gmm::transposed( _B), Bt);
std::cerr << "Bt took " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1 )
std::cerr << __FUNCTION__ << " Bt took " << s1.stop()/1000.0 << std::endl;
s1.start(); s1.start();
// setup BtB // setup BtB
CMatrixT BtB; CMatrixT BtB;
gmm::resize( BtB, n, n); gmm::resize( BtB, n, n);
gmm::mult( Bt, _B, BtB); gmm::mult( Bt, _B, BtB);
std::cerr << "BtB took " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1 )
std::cerr << __FUNCTION__ << " BtB took " << s1.stop()/1000.0 << std::endl;
s1.start(); s1.start();
// extract rhs // extract rhs
std::vector< double > rhs( n); std::vector< double > rhs( n);
gmm::copy( gmm::scaled(gmm::mat_const_col( BtB, n - 1),-1.0), rhs); gmm::copy( gmm::scaled(gmm::mat_const_col( BtB, n - 1),-1.0), rhs);
rhs.resize( n - 1); rhs.resize( n - 1);
std::cerr << "rhs extract resize " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1)
std::cerr << __FUNCTION__ << " rhs extract resize " << s1.stop()/1000.0 << std::endl;
s1.start(); s1.start();
// resize BtB to only contain the actual system matrix (and not the rhs) // resize BtB to only contain the actual system matrix (and not the rhs)
gmm::resize( BtB, n - 1, n - 1); gmm::resize( BtB, n - 1, n - 1);
std::cerr << "BtB resize took " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1)
std::cerr << __FUNCTION__ << " BtB resize took " << s1.stop()/1000.0 << std::endl;
s1.start(); s1.start();
_x.resize( n - 1); _x.resize( n - 1);
std::cerr << "x resize took " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1)
std::cerr << __FUNCTION__ << " x resize took " << s1.stop()/1000.0 << std::endl;
// regularize if necessary // regularize if necessary
if(_reg_factor != 0.0) if(_reg_factor != 0.0)
...@@ -624,7 +637,8 @@ setup_and_solve_system( CMatrixT& _B, ...@@ -624,7 +637,8 @@ setup_and_solve_system( CMatrixT& _B,
CSCMatrix BtBCSC; CSCMatrix BtBCSC;
BtBCSC.init_with_good_format( BtB); BtBCSC.init_with_good_format( BtB);
std::cerr << "CSC init " << s1.stop()/1000.0 << std::endl; if( noisy_ > 1)
std::cerr << __FUNCTION__ << " CSC init " << s1.stop()/1000.0 << std::endl;
double setup_time = sw.stop()/1000.0; double setup_time = sw.stop()/1000.0;
// create solver // create solver
...@@ -637,7 +651,8 @@ setup_and_solve_system( CMatrixT& _B, ...@@ -637,7 +651,8 @@ setup_and_solve_system( CMatrixT& _B,
misw.start(); misw.start();
// miso solve // miso solve
miso.solve( BtBCSC, _x, rhs, _idx_to_round); miso.solve( BtBCSC, _x, rhs, _idx_to_round);
std::cerr << "Miso Time " << misw.stop()/1000.0 << "s." << std::endl << std::endl; if( noisy_ > 1)
std::cerr << __FUNCTION__ << " Miso Time " << misw.stop()/1000.0 << "s." << std::endl << std::endl;
return setup_time; return setup_time;
} }
...@@ -664,7 +679,7 @@ restore_eliminated_vars( RMatrixT& _constraints, ...@@ -664,7 +679,7 @@ restore_eliminated_vars( RMatrixT& _constraints,
if( _new_idx[i] != -1) if( _new_idx[i] != -1)
{ {
// error handling // error handling
if( i < _new_idx[i]) std::cerr << "Warning: UNSAVE Ordering!!!\n"; if( i < _new_idx[i] && noisy_ > 0) std::cerr << "Warning: UNSAFE Ordering!!!\n";
_x[i] = _x[_new_idx[i]]; _x[i] = _x[_new_idx[i]];
} }
......
...@@ -566,8 +566,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar, ...@@ -566,8 +566,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar,
typedef typename gmm::linalg_traits<MatrixT>::const_sub_col_type ColT; typedef typename gmm::linalg_traits<MatrixT>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter; typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
ACG::StopWatch sw;
sw.start();
// unsigned int m = gmm::mat_nrows( _A); // unsigned int m = gmm::mat_nrows( _A);
unsigned int n = gmm::mat_ncols( _A ); unsigned int n = gmm::mat_ncols( _A );
...@@ -613,8 +611,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar, ...@@ -613,8 +611,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar,
} }
} }
std::cerr << " before matrix took " << sw.stop()/1000.0 << std::endl;
sw.start();
// delete last element // delete last element
_rhs.resize( n_new ); _rhs.resize( n_new );
_x.resize( n_new ); _x.resize( n_new );
...@@ -724,7 +720,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar, ...@@ -724,7 +720,6 @@ void eliminate_vars( const std::vector<IntegerT>& _evar,
gmm::copy( Atmp, _A); gmm::copy( Atmp, _A);
*/ */
std::cerr << " after matrix took " << sw.stop()/1000.0 << std::endl;
} }
......
...@@ -40,298 +40,11 @@ ...@@ -40,298 +40,11 @@
namespace ACG { namespace ACG {
//== IMPLEMENTATION ========================================================== //== IMPLEMENTATION ==========================================================
/*
template< class VecT>
void
MISolver::solve(
gmm::col_matrix< VecT >& _B,
VecT& _x,
Veci& _to_round,
bool _fixed_order = false )
{
typedef typename gmm::row_matrix< VecT > RMatrixT;
typedef typename gmm::col_matrix< VecT > CMatrixT;
// some statistics
int n_local = 0;
int n_cg = 0;
int n_full = 0;
Veci to_round(_to_round);
// if the order is not fixed, uniquify the indices
if( !_fixed_order)
{
// copy to round vector and make it unique
std::sort(to_round.begin(), to_round.end());
Veci::iterator last_unique;
last_unique = std::unique(to_round.begin(), to_round.end());
int r = last_unique - to_round.begin();
to_round.resize( r);
}
// setup quadratic csc system Ax=rhs
CSCMatrix A;
Vecd rhs;
setup_quadratic_csc_system(_B, A, _x, rhs );
// initalize old indices
Vecui old_idx(rhs.size());
for(unsigned int i=0; i<old_idx.size(); ++i)
old_idx[i] = i;
// Setup Cholmod solver used for full solution
ACG::CholmodSolver chol;
if( initial_full_solution_ || direct_rounding_)
{
if( noisy_ > 2) std::cerr << "initial full solution" << std::endl;
chol.calc_system_gmm(A);
chol.solve(_x, rhs);
++n_full;
}
if( no_rounding_ || _to_round.empty() )
return;
// preconditioner for CG
gmm::identity_matrix PS, PR;
// neighbors for local optimization
Vecui neigh_i;
// Vector for reduced solution
Vecd xr(_x);
// direct rounding?
if( direct_rounding_ && _fixed_order )
std::cerr << "Rounding order is fixed, direct rounding does not make sense!" << std::endl;
if( direct_rounding_ && !_fixed_order )
{
int n_ints = to_round.size();
int n_syst = _B.ncols();
int rhs_idx= n_syst-1;
// get rounded values and setup constraint matrix
Veci c_elim; //columns to be eliminated
RMatrixT C( n_ints, n_syst);
n_ints = 0;
for( uint j=0; j < to_round.size(); ++j)
{
if( to_round[j] > -1)
{
int cur_idx = to_round[j];
double rnd_x = ROUND(xr[cur_idx]);
C( j, cur_idx ) = -1;
C( j, rhs_idx ) = rnd_x;
c_elim.push_back( cur_idx);
++n_ints;
}
}
if( to_round.size() != n_ints)
gmm::resize( C, n_ints, n_syst);
// eliminate all from _B (only economic when rounding very many)
Veci empty_round_vector;
Veci new_idx;
eliminate_conditions( C, _B, empty_round_vector, c_elim, new_idx);
// form new A
setup_quadratic_csc_system( _B, A, _x, rhs);
// perform final solve
chol.calc_system_gmm(A);
chol.solve(_x, rhs);
++n_full;
// restore eliminated variables
restore_eliminated_vars( C, _x, c_elim, new_idx);
}
if( !direct_rounding_)
// loop until solution computed
for(unsigned int i=0; i<to_round.size(); ++i)
{
if( noisy_ > 0)
{
std::cerr << "Integer DOF's left: " << to_round.size()-(i+1) << " ";
if( noisy_ > 1)
std::cerr << "residuum_norm: " << gmm::residuum_norm( A, xr, rhs) << std::endl;
}
// index to eliminate
unsigned int i_best = 0;
if( _fixed_order ) // if order is fixed, simply get next index from to_round
{
i_best = to_round[i];
}
else // else search for best rounding candidate
{
// find index yielding smallest rounding error
double r_best = FLT_MAX;
for(unsigned int j=0; j<to_round.size(); ++j)
{
if( to_round[j] != -1)
{
int cur_idx = to_round[j];
double rnd_error = fabs( ROUND(xr[cur_idx]) - xr[cur_idx]);
if( rnd_error < r_best)
{
i_best = cur_idx;
r_best = rnd_error;
}
}
}
}
// store rounded value
double rnd_x = ROUND(xr[i_best]);
_x[ old_idx[i_best] ] = rnd_x;
// compute neighbors
neigh_i.clear();
Col col = gmm::mat_const_col(A, i_best);
ColIter it = gmm::vect_const_begin( col);
ColIter ite = gmm::vect_const_end ( col);
for(; it!=ite; ++it)
neigh_i.push_back(it.index());
// eliminate var
gmm::eliminate_var(i_best, rnd_x, A, xr, rhs);
gmm::eliminate_var_idx(i_best, to_round);
// update old_idx
old_idx.erase( old_idx.begin()+i_best );
if( direct_rounding_ && !_fixed_order) continue;
// current error
double cur_error = FLT_MAX;
// compute new solution
unsigned int n_its = max_local_iters_;
if(max_local_iters_ > 0)
{
if( noisy_ > 2)std::cerr << "use local iteration ";
n_its = gmm::gauss_seidel_local(A, xr, rhs, neigh_i, max_local_iters_, max_local_error_);
++n_local;
}
// update error bound
// if gauss_seidel did not reach max iters then error must be less than
// tolerance (max_local_error_)
if( n_its < max_local_iters_) cur_error = max_local_error_;
if( cur_error > max_cg_error_)
{
if( noisy_ > 2) std::cerr << ", cg ";
gmm::iteration iter(max_cg_error_);
iter.set_maxiter(max_cg_iters_);
iter.set_noisy(std::max( int(noisy_)-4, int(0)));
// conjugate gradient
if( max_cg_iters_ > 0)
{
gmm::cg( A, xr, rhs, PS, PR, iter);
cur_error = iter.get_res();
++n_cg;
}
}
if( cur_error > max_full_error_ )
{
if( noisy_ > 2)std::cerr << ", full ";
if( gmm::mat_ncols( A) > 0)
{
chol.calc_system_gmm(A);
chol.solve(xr,rhs);
++n_full;
}
}
if( noisy_ > 2)std::cerr << std::endl;
}
// final full solution?
if( final_full_solution_ || direct_rounding_)
{
if( noisy_ > 2) std::cerr << "final full solution" << std::endl;
if( gmm::mat_ncols( A) > 0)
{
chol.calc_system_gmm(A);
chol.solve( xr, rhs);
++n_full;
}
}
// store solution values to result vector
for(unsigned int i=0; i<old_idx.size(); ++i)
{
_x[ old_idx[i] ] = xr[i];
}
// output statistics
if( stats_)
{
std::cerr << "\t" << __FUNCTION__ << " *** Statistics of MiSo Solver ***\n";
std::cerr << "\t\t Number of CG iterations = " << n_cg << std::endl;
std::cerr << "\t\t Number of LOCAL iterations = " << n_local << std::endl;
std::cerr << "\t\t Number of FULL iterations = " << n_full << std::endl;
std::cerr << "\t\t Number of ROUNDING = " << _to_round.size() << std::endl;
std::cerr << std::endl;
}
}
template< class CMatrixT >
void
MISolver::setup_quadratic_csc_system(
CMatrixT& _B
CSCMatrix& _A,
Vecd& _x,
Vecd& _rhs )
{
// clear _A
_A.do_clear();
// clear rhs
_rhs.clear();
unsigned int m = gmm::mat_nrows(_B);
unsigned int n = gmm::mat_ncols(_B);
// set up B transposed
CMatrixT Bt( n ,m);
gmm::copy( gmm::transposed( _B), Bt);
// setup BtB
CMatrixT BtB( n, n);
gmm::mult( Bt, _B, BtB);
// extract rhs
_rhs.resize( n);
gmm::copy( gmm::scaled(gmm::mat_const_col( BtB, n - 1),-1.0), rhs);
_rhs.resize( n - 1);
// resize BtB to only contain the actual system matrix (and not the rhs)
gmm::resize( BtB, n - 1, n - 1);
_x.resize( n - 1);
// regularize if necessary
//if(_reg_factor != 0.0)
// gmm::regularize_hack(BtB, _reg_factor);
// BtB -> CSC
_A.init_with_good_format( BtB);
}
*/
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
//============================================================================= //=============================================================================
} // namespace ACG } // namespace ACG
//============================================================================= //=============================================================================
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