songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Remove assumption in QP_solver that BigInt(double) is implicit. 

Many complicated functors all over the place, I'll be lucky if I didn't
break at least one...
This commit is contained in:
Marc Glisse 2017-11-29 22:32:52 +01:00
parent f69e266263
commit ab04858562
5 changed files with 88 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
Bounding_volumes/include/CGAL/Min_annulus_d.h
View File

@ -466,7 +466,8 @@ public:
{ CGAL_optimisation_precondition(
is_empty() || tco.access_dimension_d_object()( p) == d);
ET sqr_d = sqr_dist( p);
ET h_p_sqr = da_coord(p)[d] * da_coord(p)[d];
ET h_p_sqr(da_coord(p)[d]);
h_p_sqr *= h_p_sqr;
return ( ( sqr_d < h_p_sqr * sqr_i_rad_numer) ||
( h_p_sqr * sqr_o_rad_numer < sqr_d)); }
@ -623,9 +624,10 @@ private:
inner_indices.push_back( 0);
outer_indices.push_back( 0);
center_coords.resize( d+1);
std::copy( da_coord( points[ 0]),
da_coord( points[ 0])+d+1,
center_coords.begin());
std::transform( da_coord( points[ 0]),
da_coord( points[ 0])+d+1,
center_coords.begin(),
typename Coercion_traits<ET,RT>::Cast());
sqr_i_rad_numer = ET( 0);
sqr_o_rad_numer = ET( 0);
sqr_rad_denom = ET( 1);
@ -748,7 +750,8 @@ is_valid( bool verbose, int level) const
// all inner support points on inner boundary?
Inner_support_point_iterator i_pt_it = inner_support_points_begin();
for ( ; i_pt_it != inner_support_points_end(); ++i_pt_it) {
ET h_p_sqr = da_coord (*i_pt_it)[d] * da_coord (*i_pt_it)[d];
ET h_p_sqr(da_coord (*i_pt_it)[d]);
h_p_sqr *= h_p_sqr;
if ( sqr_dist( *i_pt_it) != h_p_sqr * sqr_i_rad_numer)
return CGAL::_optimisation_is_valid_fail( verr,
"annulus does not have all inner support points on its inner boundary");
@ -757,7 +760,8 @@ is_valid( bool verbose, int level) const
// all outer support points on outer boundary?
Outer_support_point_iterator o_pt_it = outer_support_points_begin();
for ( ; o_pt_it != outer_support_points_end(); ++o_pt_it) {
ET h_p_sqr = da_coord (*o_pt_it)[d] * da_coord (*o_pt_it)[d];
ET h_p_sqr(da_coord (*o_pt_it)[d]);
h_p_sqr *= h_p_sqr;
if ( sqr_dist( *o_pt_it) != h_p_sqr * sqr_o_rad_numer)
return CGAL::_optimisation_is_valid_fail( verr,
"annulus does not have all outer support points on its outer boundary");

11
QP_solver/include/CGAL/QP_solver/Initialization.h
View File

@ -549,7 +549,8 @@ void QP_solver<Q, ET, Tags>::
init_solution__b_C(Tag_true)
{
b_C.reserve(qp_m);
std::copy(qp_b, qp_b+qp_m, std::back_inserter(b_C));
std::transform(qp_b, qp_b+qp_m, std::back_inserter(b_C),
typename Coercion_traits<ET,B_entry>::Cast());
}
template < typename Q, typename ET, typename Tags > inline // has ineq.
@ -559,9 +560,11 @@ init_solution__b_C(Tag_false)
b_C.insert(b_C.end(), l, et0);
B_by_index_accessor b_accessor(qp_b); // todo kf: is there some boost
// replacement for this accessor?
std::copy(B_by_index_iterator(C.begin(), b_accessor),
B_by_index_iterator(C.end (), b_accessor),
b_C.begin());
typedef typename std::iterator_traits<B_by_index_iterator>::value_type RT;
std::transform(B_by_index_iterator(C.begin(), b_accessor),
B_by_index_iterator(C.end (), b_accessor),
b_C.begin(),
typename Coercion_traits<ET,RT>::Cast());
}
// initial solution

36
QP_solver/include/CGAL/QP_solver/QP_solver.h
View File

@ -1599,14 +1599,18 @@ ratio_test_1__q_x_S( Tag_false)
// ( A_S_BxB_O * q_x_O) - A_S_Bxj
if ( j < qp_n) {
typedef typename std::iterator_traits<A_by_index_iterator>::value_type RT;
std::transform( q_x_S.begin(),
q_x_S.begin()+S_B.size(),
A_by_index_iterator( S_B.begin(),
A_by_index_accessor( *(qp_A + j))),
q_x_S.begin(),
compose2_2( std::minus<ET>(),
Identity<ET>(),
boost::bind1st( std::multiplies<ET>(), d)));
boost::bind(std::minus<ET>(),
boost::placeholders::_1,
boost::bind(std::multiplies<ET>(), d,
boost::bind(
typename Coercion_traits<ET,RT>::Cast(),
boost::placeholders::_2))));
}
// q_x_S = -+ ( A_S_BxB_O * q_x_O - A_S_Bxj)
@ -1888,9 +1892,11 @@ basis_matrix_stays_regular()
new_row = slack_A[ i-qp_n].first;
A_row_by_index_accessor a_accessor =
boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin());
typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin(),
typename Coercion_traits<ET,RT>::Cast());
inv_M_B.multiply( tmp_x.begin(), // dummy (not used)
tmp_x.begin(), tmp_l_2.begin(), tmp_x_2.begin(),
Tag_false(), // QP
@ -1963,7 +1969,18 @@ compute__x_B_S( Tag_false /*has_equalities_only_and_full_rank*/,
}
namespace QP_solver_impl {
// Writing it with 5 boost::bind was becoming unreadable.
template<class ET, class RT>
struct submul {
ET const& d;
submul(ET const&d):d(d) {}
ET operator()(RT const&x, ET const&y) const {
typename Coercion_traits<ET, RT>::Cast cast;
return cast(x) * d - y;
}
};
}
template < typename Q, typename ET, typename Tags > inline // has inequalities, standard form
void QP_solver<Q, ET, Tags>::
@ -1975,13 +1992,12 @@ compute__x_B_S( Tag_false /*has_equalities_only_and_full_rank*/,
// b_S_B - ( A_S_BxB_O * x_B_O)
B_by_index_accessor b_accessor( qp_b);
typedef typename std::iterator_traits<B_by_index_iterator>::value_type RT;
std::transform( B_by_index_iterator( S_B.begin(), b_accessor),
B_by_index_iterator( S_B.end (), b_accessor),
x_B_S.begin(),
x_B_S.begin(),
compose2_2( std::minus<ET>(),
boost::bind1st( std::multiplies<ET>(), d),
Identity<ET>()));
QP_solver_impl::submul<ET,RT>(d));
// x_B_S = +- ( b_S_B - A_S_BxB_O * x_B_O)
Value_iterator x_it = x_B_S.begin();

49
QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
View File

@ -23,6 +23,7 @@
// Kaspar Fischer
#include <CGAL/QP_solver/Initialization.h>
#include <boost/bind.hpp>
namespace CGAL {
@ -67,9 +68,13 @@ transition( )
// initialize exact version of `-qp_c' (implicit conversion to ET)
C_by_index_accessor c_accessor( qp_c);
typedef typename std::iterator_traits<C_by_index_iterator>::value_type RT;
std::transform( C_by_index_iterator( B_O.begin(), c_accessor),
C_by_index_iterator( B_O.end (), c_accessor),
minus_c_B.begin(), std::negate<ET>());
minus_c_B.begin(),
boost::bind(
typename Coercion_traits<ET,RT>::Cast(),
boost::bind(std::negate<RT>(), boost::placeholders::_1)));
// compute initial solution of phase II
compute_solution(Is_nonnegative());
@ -419,7 +424,8 @@ ratio_test_init__A_Cj( Value_iterator A_Cj_it, int j_, Tag_true)
// store exact version of `A_Cj' (implicit conversion)
if ( j_ < qp_n) { // original variable
CGAL::cpp11::copy_n( *(qp_A + j_), qp_m, A_Cj_it);
CGAL::transform_n( *(qp_A + j_), qp_m, A_Cj_it,
typename Coercion_traits<ET,A_entry>::Cast());
} else { // artificial variable
@ -437,9 +443,11 @@ ratio_test_init__A_Cj( Value_iterator A_Cj_it, int j_, Tag_false)
// store exact version of `A_Cj' (implicit conversion)
if ( j_ < qp_n) { // original variable
A_by_index_accessor a_accessor( *(qp_A + j_));
std::copy( A_by_index_iterator( C.begin(), a_accessor),
A_by_index_iterator( C.end (), a_accessor),
A_Cj_it);
typedef typename std::iterator_traits<A_by_index_iterator>::value_type RT;
std::transform(A_by_index_iterator( C.begin(), a_accessor),
A_by_index_iterator( C.end (), a_accessor),
A_Cj_it,
typename Coercion_traits<ET,RT>::Cast());
} else {
unsigned int k = j_;
@ -461,9 +469,11 @@ ratio_test_init__A_Cj( Value_iterator A_Cj_it, int j_, Tag_false)
} else { // special art.
S_by_index_accessor s_accessor( art_s.begin());
std::copy( S_by_index_iterator( C.begin(), s_accessor),
S_by_index_iterator( C.end (), s_accessor),
A_Cj_it);
typedef typename std::iterator_traits<S_by_index_iterator>::value_type RT;
std::transform(S_by_index_iterator( C.begin(), s_accessor),
S_by_index_iterator( C.end (), s_accessor),
A_Cj_it,
typename Coercion_traits<ET,RT>::Cast());
}
}
}
@ -1347,7 +1357,7 @@ replace_variable_original_original( )
minus_c_B[ k] =
( is_phaseI ?
( j < qp_n ? et0 : -aux_c[j-qp_n-slack_A.size()]) : -ET( *(qp_c+ j)));
( j < qp_n ? et0 : ET(-aux_c[j-qp_n-slack_A.size()])) : -ET( *(qp_c+ j)));
if ( is_phaseI) {
if ( j >= qp_n) ++art_basic;
@ -1437,9 +1447,11 @@ replace_variable_slack_slack( )
// update basis inverse
A_row_by_index_accessor a_accessor =
boost::bind( A_accessor( qp_A, 0, qp_n), _1, new_row);
std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin());
typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin(),
typename Coercion_traits<ET,RT>::Cast());
if ( art_s_i > 0) { // special artificial
tmp_x[ in_B[ art_s_i]] = ET( art_s[ new_row]);
}
@ -1557,7 +1569,7 @@ replace_variable_original_slack( )
minus_c_B[ B_O.size()]
= ( is_phaseI ?
( j < qp_n ? et0 : -aux_c[j-qp_n-slack_A.size()])
( j < qp_n ? et0 : ET(-aux_c[j-qp_n-slack_A.size()]))
: -ET( *(qp_c+ j)));
@ -1588,9 +1600,11 @@ replace_variable_original_slack( )
// update basis inverse
A_row_by_index_accessor a_accessor =
boost::bind (A_accessor( qp_A, 0, qp_n), _1, new_row);
std::copy( A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin());
typedef typename std::iterator_traits<A_row_by_index_iterator>::value_type RT;
std::transform(A_row_by_index_iterator( B_O.begin(), a_accessor),
A_row_by_index_iterator( B_O.end (), a_accessor),
tmp_x.begin(),
typename Coercion_traits<ET,RT>::Cast());
if ( art_s_i > 0) { // special art.
tmp_x[ in_B[ art_s_i]] = ET( art_s[ new_row]);
}
@ -2618,7 +2632,8 @@ multiply__A_S_BxB_O(Value_iterator in, Value_iterator out) const
// foreach row of A in S_B
for ( row_it = S_B.begin(); row_it != S_B.end(); ++row_it,
++out_it) {
*out_it += ET( a_col[ *row_it]) * in_value;
A_entry x = a_col[*row_it];
*out_it += ET(x) * in_value;
}
} else {
if ( *col_it == art_s_i) { // special artificial

13
STL_Extension/include/CGAL/algorithm.h
View File

@ -455,6 +455,19 @@ void nth_element(RandomAccessIterator left,
} // end while(true)
}
// Not a standard function, but close enough
template <class InputIterator, class Size, class OutputIterator, class UnaryFun>
OutputIterator transform_n( InputIterator first, Size n, OutputIterator result, UnaryFun fun)
{
// copies the first `n' transformed items from `first' to `result'.
// Returns the value of `result' after inserting the `n' items.
while( n--) {
*result = fun(*first);
first++;
result++;
}
return result;
}
} //namespace CGAL
#endif // CGAL_ALGORITHM_H