00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BonQuadRow.hpp"
00011 #include <cfloat>
00012
00013 namespace Bonmin{
00014
00015 QuadRow::QuadRow():
00016 c_(0),
00017 a_(),
00018 Q_(),
00019 grad_evaled_(false)
00020 {
00021 }
00022
00023 QuadRow::QuadRow(const QuadRow &other):
00024 c_(other.c_),
00025 a_(other.a_),
00026 Q_(other.Q_),
00027 g_(),
00028 a_grad_idx_(),
00029 Q_row_grad_idx_(),
00030 Q_col_grad_idx_(),
00031 Q_hessian_idx_(),
00032 grad_evaled_(false)
00033 {
00034 initialize();
00035 }
00036
00037 QuadRow & QuadRow::operator=(const QuadRow &rhs){
00038 if(this != &rhs){
00039 c_ = rhs.c_;
00040 a_ = rhs.a_;
00041 Q_ = rhs.Q_;
00042 Q_hessian_idx_.clear();
00043 g_.clear();
00044 a_grad_idx_.clear();
00045 Q_row_grad_idx_.clear();
00046 Q_col_grad_idx_.clear();
00047 initialize();
00048
00049 grad_evaled_ = false;
00050 }
00051 return (*this);
00052 }
00053
00054 QuadRow::QuadRow(const QuadCut &cut):
00055 c_(0),
00056 a_(cut.row()),
00057 Q_(cut.Q(), cut.type())
00058 {
00059 initialize();
00060 }
00061
00062 QuadRow&
00063 QuadRow::operator=(const QuadCut &cut){
00064 c_ = cut.c();
00065 a_ = cut.row();
00066 Q_ = cut.Q();
00067 Q_.make_upper_triangular(cut.type());
00068 g_.clear();
00069 a_grad_idx_.clear();
00070 Q_row_grad_idx_.clear();
00071 Q_col_grad_idx_.clear();
00072
00073
00074 initialize();
00075 return (*this);
00076 }
00077
00078
00079 QuadRow::QuadRow(const OsiRowCut &cut):
00080 c_(0),
00081 a_(cut.row()),
00082 Q_()
00083 {
00084 initialize();
00085 }
00086
00087 QuadRow&
00088 QuadRow::operator=(const OsiRowCut &cut){
00089 c_ = 0;
00090 a_ = cut.row();
00091 Q_ = TMat();
00092 g_.empty();
00093 a_grad_idx_.empty();
00094 Q_row_grad_idx_.clear();
00095 Q_col_grad_idx_.clear();
00096
00097
00098 initialize();
00099 return (*this);
00100 }
00101
00102 void
00103 QuadRow::initialize(){
00104
00105 for(int i = 0 ; i < Q_.nnz_ ; i++){
00106 assert(Q_.jCol_[i] >= Q_.iRow_[i]);}
00107 grad_evaled_ = false;
00108
00109 int n;
00110
00111 n = a_.getNumElements();
00112 a_grad_idx_.reserve(n);
00113
00114
00115 const int * indices = a_.getIndices();
00116 const double * elems = a_.getElements();
00117
00118 for(int i = 0 ; i < n ; i++){
00119 std::pair<gStore::iterator, bool> res = g_.insert(std::make_pair(indices[i], std::make_pair(elems[i],0.)));
00120 a_grad_idx_.push_back(res.first);
00121 }
00122
00123 n = Q_.numNonEmptyRows();
00124 const TMat::RowS& nonEmptyRows = Q_.nonEmptyRows();
00125 Q_row_grad_idx_.reserve(n);
00126
00127 for(TMat::RowS::const_iterator i = nonEmptyRows.begin() ; i != nonEmptyRows.end() ; i++){
00128 std::pair<gStore::iterator, bool> res = g_.insert(std::make_pair(i->first, std::make_pair(0.,0.)));
00129 Q_row_grad_idx_.push_back(res.first);
00130 }
00131
00132
00133 n = Q_.numNonEmptyCols();
00134 const TMat::RowS& nonEmptyCols = Q_.nonEmptyCols();
00135 Q_col_grad_idx_.reserve(n);
00136
00137 for(TMat::RowS::const_iterator i = nonEmptyCols.begin() ; i != nonEmptyCols.end() ; i++){
00138 std::pair<gStore::iterator, bool> res = g_.insert(std::make_pair(i->first, std::make_pair(0.,0.)));
00139 Q_col_grad_idx_.push_back(res.first);
00140 }
00141
00142 }
00143
00145 void
00146 QuadRow::print(){
00147 std::cout<<"constant term "<<c_<<std::endl;
00148 int * a_ind = a_.getIndices();
00149 const double * a_el = a_.getElements();
00150 int n = a_.getNumElements();
00151 std::cout<<"Linear term (size "<<n<<"): ";
00152 for(int i = 0 ; i < n ; i++)
00153 {
00154 std::cout<<a_el[i]<<" * x["<<a_ind[i]<<"]\t";
00155 if(i && ( (i % 5) == 0 ) ) std::cout<<std::endl<<"\t\t";
00156 }
00157 }
00159 double
00160 QuadRow::eval_f(const double *x, bool new_x){
00161
00162 internal_eval_grad(x);
00163 double value = c_;
00164
00165
00166 int * a_ind = a_.getIndices();
00167 const double * a_el = a_.getElements();
00168 int n = a_.getNumElements();
00169 for(int i = 0 ; i < n ; i++)
00170 {
00171 value += a_el[i] * x[a_ind[i]];
00172
00173 }
00174
00175
00176 for(gStore::iterator i = g_.begin() ; i != g_.end() ; i++){
00177 value += i->second.second * x[i->first];
00178 }
00179 return value;
00180 }
00181
00183 int
00184 QuadRow::nnz_grad(){
00185 return static_cast<int>(g_.size());}
00187 void
00188 QuadRow::gradiant_struct(const int nnz, int * indices, bool offset){
00189 int n = 0;
00190 for(gStore::iterator i = g_.begin() ; i != g_.end() ; i++){
00191 indices[n++] = i->first + offset;
00192 }
00193 assert(n == nnz);
00194 assert(nnz == (int) g_.size());
00195 }
00196
00198 void
00199 QuadRow::eval_grad(const int nnz, const double * x, bool new_x, double * values){
00200
00201 #ifdef DEBUG
00202
00203 for(gStore::iterator i = g_.begin() ; i != g_.end() ; i++){
00204 printf("x[%i] = %g, ",i->first, x[i->first]);
00205 }
00206 #endif
00207
00208 internal_eval_grad(x);
00209 int n = 0;
00210 #ifdef DEBUG
00211 std::cout<<"Computing gradient"<<std::endl;
00212 #endif
00213 for(gStore::iterator i = g_.begin() ; i != g_.end() ; i++){
00214 #ifdef DEBUG
00215 printf("%i: %g, %g\n", i->first, i->second.second, i->second.first);
00216 #endif
00217 values[n++] = 2*i->second.second + i->second.first;
00218 }
00219 assert (nnz == (int) g_.size());
00220 }
00221
00222 void
00223 QuadRow::internal_eval_grad(const double *x){
00224
00225 for(gStore::iterator i = g_.begin() ; i != g_.end() ; i++){
00226 i->second.second = 0;
00227 }
00228
00229
00230 const TMat::RowS & nonEmptyRows = Q_.nonEmptyRows();
00231 int k = 0;
00232 for(TMat::RowS::const_iterator ii = nonEmptyRows.begin() ; ii != nonEmptyRows.end();
00233 ii++, k++){
00234 double value = 0;
00235 assert(ii->first == Q_.iRow_[Q_.rowOrdering_[ii->second]]);
00236 for(int i = ii->second ; i < Q_.nnz_ && ii->first == Q_.iRow_[Q_.rowOrdering_[i]] ; i++)
00237 {
00238 value += x[Q_.jCol_[Q_.rowOrdering_[i]]] * Q_.value_[Q_.rowOrdering_[i]];
00239 }
00240 Q_row_grad_idx_[k]->second.second += value;
00241 assert(Q_row_grad_idx_[k]->first == ii->first);
00242 }
00243 const TMat::RowS & nonEmptyCols = Q_.nonEmptyCols();
00244 k = 0;
00245 for(TMat::RowS::const_iterator ii = nonEmptyCols.begin() ; ii != nonEmptyCols.end();
00246 ii++, k++){
00247 double value = 0;
00248 assert(ii->first == Q_.jCol_[Q_.columnOrdering_[ii->second]]);
00249 for(int i = ii->second ; i < Q_.nnz_ && ii->first == Q_.jCol_[Q_.columnOrdering_[i]] ; i++)
00250 {
00251 if(Q_.iRow_[Q_.columnOrdering_[i]] != Q_.jCol_[Q_.columnOrdering_[i]])
00252 value += x[Q_.iRow_[Q_.columnOrdering_[i]]] * Q_.value_[Q_.columnOrdering_[i]];
00253 }
00254 Q_col_grad_idx_[k]->second.second += value;
00255 assert(Q_col_grad_idx_[k]->first == ii->first);
00256 }
00257
00258 grad_evaled_ = true;
00259 }
00260
00261 void
00262 QuadRow::add_to_hessian(AdjustableMat &H, bool offset){
00263 assert(Q_hessian_idx_.empty());
00264 for(int i = 0 ; i < Q_.nnz_ ; i++){
00265 std::pair<int, int> e;
00266 e = std::make_pair(Q_.jCol_[i] + offset, Q_.iRow_[i] + offset);
00267 AdjustableMat::iterator pos = H.find(e);
00268 if(pos != H.end()){
00269 if(pos->second.second != -1)
00270 pos->second.second++;
00271 Q_hessian_idx_.push_back(pos);
00272 }
00273 else {
00274 std::pair<AdjustableMat::iterator, bool> res =
00275 H.insert(std::make_pair(e, std::make_pair(H.size(), 1)));
00276 assert(res.second == true);
00277 Q_hessian_idx_.push_back(res.first);
00278 }
00279 }
00280 }
00281
00282 void
00283 QuadRow::remove_from_hessian(AdjustableMat &H){
00284 for(int i = 0 ; i < Q_.nnz_ ; i++){
00285 if(Q_hessian_idx_[i]->second.second != -1)
00286 Q_hessian_idx_[i]->second.second--;
00287 if(Q_hessian_idx_[i]->second.second == 0){
00288 H.erase(Q_hessian_idx_[i]);
00289 }
00290 }
00291 Q_hessian_idx_.clear();
00292 }
00293
00295 void
00296 QuadRow::eval_hessian(double lambda, double * values){
00297 for(int i = 0 ; i < Q_.nnz_ ; i++){
00298 #ifdef DEBUG
00299 printf("iRow %i, jCol %i, value %g , nnz %i\n",
00300 Q_hessian_idx_[i]->first.second,
00301 Q_hessian_idx_[i]->first.first,
00302 Q_.value_[i],
00303 Q_hessian_idx_[i]->second.first);
00304 #endif
00305 values[Q_hessian_idx_[i]->second.first] += (lambda * 2 * Q_.value_[i]);
00306 }
00307 }
00308
00309 }
00310