1 //
2 // Type.h
3 // Gravity
4 //
5 // Created by Hassan on 3 Jan 2016.
6 //
7
8 #ifndef Gravity___Type_h
9 #define Gravity___Type_h
10 #include <memory>
11 #include <list>
12 #define _USE_MATH_DEFINES
13 #include <math.h>
14 #include <map>
15 #include <set>
16 #include <assert.h>
17 #include <string>
18 #include <iostream>
19 #include <algorithm>
20 #include <complex> // std::complex
21
nthOccurrence(const std::string & str,const std::string & findMe,int nth)22 inline int nthOccurrence(const std::string& str, const std::string& findMe, int nth)
23 {
24 size_t pos = 0;
25 int cnt = 0;
26
27 while( cnt != nth )
28 {
29 pos+=1;
30 pos = str.find(findMe, pos);
31 if ( pos == std::string::npos )
32 return -1;
33 cnt++;
34 }
35 return pos;
36 }
37
38 namespace gravity{
39 #define EPS 0.00001
40 #define Cpx complex<double>
41 //#define Real double
42 //#define Integer integer //same name with a boost graph library.
43 #define Binary bool
44 #define Debug(x)
45 #define DebugOn(x) cout << x
46 #define Warning(x) cout << x
47 #define WarningOff(x)
48 #define DebugOff(x)
49
50 typedef unsigned int uint; /* Index type */
51 //typedef std::set<ind> indx; /* Set of indices type */
52
53 typedef enum { linear_, convex_, concave_, undet_} Convexity; /* Convexity Type */
54 typedef enum { neg_ = -2, non_pos_ = -1, zero_ = 0, non_neg_ = 1, pos_ = 2, unknown_ = 3} Sign; /* Sign Type */
55 typedef enum { binary_, short_, integer_, float_, double_, long_, complex_} NType; /* Number Type */
56 typedef enum { binary_c, short_c, integer_c, float_c, double_c, long_c, par_c, uexp_c, bexp_c, var_c, func_c, complex_c} CType; /* Constant type, ancestor to parameter, var and function */
57 typedef enum { Unknown,
58 PrimalAndDualFeasible,
59 PrimalFeasible,
60 DualFeasible,
61 PrimalInfeasible,
62 DualInfeasible,
63 PrimalAndDualInfeasible,
64 IllPosed,
65 PrimalInfeasibleOrUnbounded } Outcome;
66 typedef enum { geq, leq, eq } ConstraintType;
67 typedef enum { const_, lin_, quad_, pol_, nlin_ } FType; /* Function type in constraint: Constant, Linear, Quadratic, Polynomial or Nonlinear function */
68 typedef enum { lin_m, quad_m, pol_m, nlin_m } MType; /* Model type: Linear, Quadratic, Polynomial or Nonlinear function */
69 typedef enum { minimize, maximize } ObjectiveType;
70 typedef enum { id_, plus_, minus_, product_, div_, power_, cos_, sin_, sqrt_, exp_, log_, tan_, acos_, asin_, atan_, atan2_, abs_, df_abs_, relu_, unit_step_, min_, max_} OperatorType; /* Operation type in the expression tree */
71
72 typedef enum { R_, R_p_, C_} SpaceType; /* Real, Positive Reals, Complex */
73
74 typedef enum { ordered_pairs_, unordered_ } SetType;
75 // typedef enum { vec_=0, in_ordered_pairs_=1, from_ordered_pairs_=2, to_ordered_pairs_=3, in_arcs_=4, from_arcs_=5, to_arcs_=6, in_nodes_=7, in_set_=8, mask_=9, in_bags_=10, time_expand_ = 11, in_set_at_} IndexType; /* Index type */
76
77 typedef enum { unindexed_, in_, in_pairs_, out_, from_, to_, prev_, in_at_, in_time_, from_time_, to_time_, in_arcs_, out_arcs_, in_gens_, in_pot_gens_, in_bats_, in_pot_bats_,in_wind_, in_pv_, min_time_, excl_, matrix_} IndexType; /* Index type */
78
79
80
81
82
83 using namespace std;
84
85 static double pi = 4.*atan(1.);
86
87
88 class space{
89 public:
90 SpaceType _type;
91 vector<size_t> _dim;
92 };
93
94 class R: public space{
95 public:
R()96 R(){};
97 template<typename... Args>
R(size_t t1,Args &&...args)98 R(size_t t1, Args&&... args) {
99 _type = R_;
100 list<size_t> dims = {forward<size_t>(args)...};
101 dims.push_front(t1);
102 size_t size = dims.size();
103 _dim.resize(size);
104 auto it = dims.begin();
105 size_t index = 0;
106 while (it!=dims.end()) {
107 _dim[index++] = *it++;
108 }
109 }
110
111 R operator^(size_t n){return R(n);};
112
113 };
114
115 class R_p: public space{
116 public:
R_p()117 R_p(){};
118 template<typename... Args>
R_p(size_t t1,Args &&...args)119 R_p(size_t t1, Args&&... args) {
120 _type = R_p_;
121 list<size_t> dims = {forward<size_t>(args)...};
122 dims.push_front(t1);
123 size_t size = dims.size();
124 _dim.resize(size);
125 auto it = dims.begin();
126 size_t index = 0;
127 while (it!=dims.end()) {
128 _dim[index++] = *it++;
129 }
130 }
131
132 R_p operator^(size_t n){return R_p(n);};
133
134 };
135
136 class C: public space{
137 public:
C()138 C(){};
139 template<typename... Args>
C(size_t t1,Args &&...args)140 C(size_t t1, Args&&... args) {
141 _type = C_;
142 list<size_t> dims = {forward<size_t>(args)...};
143 dims.push_front(t1);
144 size_t size = dims.size();
145 _dim.resize(size);
146 auto it = dims.begin();
147 size_t index = 0;
148 while (it!=dims.end()) {
149 _dim[index++] = *it++;
150 }
151 }
152
153 C operator^(size_t n){return C(n);};
154 /* TODO */
155 };
156
157
158
159 /** Class for manipulating indices */
160 class index_{
161 public:
162 string _name;
163 string _type_name="indices";
164 bool _active = true;
_name(name)165 index_(const string& name, bool active=true):_name(name), _active(active){};
index_(const index_ & idx)166 index_(const index_& idx):_name(idx._name), _active(idx._active){};
167 template<typename... Args>
index_(string t1,Args &&...args)168 index_(string t1, Args&&... args) {
169 list<string> indices;
170 indices = {forward<Args>(args)...};
171 indices.push_front(t1);
172 auto it = indices.begin();
173 for (size_t i= 0; i < indices.size(); i++) {
174 _name += *it;
175 if (i< indices.size()-1) {
176 _name += ",";
177 }
178 it++;
179 }
180 }
181 };
182
183 class index_pair{
184 public:
185 string _name;
186 string _type_name="index_pairs";
187 bool _active = true;
188 index_* _src = nullptr;
189 index_* _dest = nullptr;
190 index_pair(const index_& src, const index_& dest, bool active = true):_name(src._name+","+dest._name), _active(active), _src(new index_(src)), _dest(new index_(dest)){};
191 index_pair(const string& src, const string& dest, bool active = true):_name(src+","+dest), _active(active), _src(new index_(src)), _dest(new index_(dest)){};
index_pair(const index_pair & p)192 index_pair(const index_pair& p):_name(p._name), _active(p._active), _src(new index_(*p._src)), _dest(new index_(*p._dest)){};
~index_pair()193 ~index_pair(){
194 delete _src;
195 _src = nullptr;
196 delete _dest;
197 _dest = nullptr;
198 }
199 };
200
201 class ordered_pairs{
202
203 public:
204 size_t _first;
205 size_t _last;
206 string _type_name="ordered_pairs";
207 std::vector<index_pair*> _keys;
ordered_pairs(size_t p1,size_t p2)208 ordered_pairs(size_t p1 ,size_t p2){
209 _first = p1;
210 _last = p2;
211 auto n = p2 - p1 + 1;
212 assert(n >= 0);
213 _keys.resize(n*(n-1)/2);
214 string key;
215 size_t index = 0;
216 for (int i = p1-1; i < p2; i++){
217 for (int j = i+1; j < p2; j++){
218 _keys[index++] = new index_pair(index_(to_string(i)), index_(to_string(j)));
219 }
220 }
221 }
ordered_pairs(size_t p1,size_t p2,bool rev)222 ordered_pairs(size_t p1 ,size_t p2, bool rev){
223 _first = p1;
224 _last = p2;
225 auto n = p2 - p1 + 1;
226 assert(n >= 0);
227 _keys.resize(n*(n-1)/2);
228 string key;
229 size_t index = 0;
230 for (int i = p1-1; i < p2; i++){
231 for (int j = i+1; j < p2; j++){
232 _keys[index++] = new index_pair(index_(to_string(j)), index_(to_string(i)));
233 }
234 }
235 }
~ordered_pairs()236 ~ordered_pairs(){
237 for (auto p: _keys) { delete p;}
238 }
239 };
240
241 class indices{
242 private:
243 string _name;/**< index set can be given a name */
244
245 public:
246
247
248 IndexType _type = unindexed_;/**< index type */
249 bool _time_extended = false;/*<< indices are time extended */
250 size_t _time_pos = 0;/*<< number of commas before time extension */
251 shared_ptr<vector<size_t>> _dim = nullptr;/*<< A vector storing the dimension of sub-indices */
252 shared_ptr<vector<string>> _keys = nullptr; /*<< A vector storing all the keys */
253
254 shared_ptr<map<string,size_t>> _keys_map = nullptr; /*<< A map storing all the indices, the size_t number indicates the right position in the _keys vector */
255
256 set<size_t> _excluded_keys; /*<< A set storing all indices that should be excluded */
257 shared_ptr<vector<vector<size_t>>> _ids = nullptr;
258
get_name()259 string get_name() const{
260 string name = _name;
261 if(_type==from_){
262 name = "from("+name+")";
263 }
264 else if(_type==to_){
265 name = "to("+name+")";
266 }
267 return name;
268 }
269
indices(string name)270 indices(string name){
271 _name = name;
272 _keys_map = make_shared<map<string,size_t>>();
273 _keys = make_shared<vector<string>>();
274 _dim = make_shared<vector<size_t>>();
275 _dim->resize(2,0);
276 }
277
set_name(const string & name)278 void set_name(const string& name){
279 _name = name;
280 }
281
282
indices(int name)283 indices(int name){
284 _name = to_string(name);
285 _keys_map = make_shared<map<string,size_t>>();
286 _keys = make_shared<vector<string>>();
287 _dim = make_shared<vector<size_t>>();
288 _dim->resize(2,0);
289 }
290
indices()291 indices(){
292 _keys_map = make_shared<map<string,size_t>>();
293 _keys = make_shared<vector<string>>();
294 _dim = make_shared<vector<size_t>>();
295 _dim->resize(2,0);
296 }
297
indices(const ordered_pairs & pairs)298 indices(const ordered_pairs& pairs){
299 auto n = pairs._keys.size();
300 _keys_map = make_shared<map<string,size_t>>();
301 _keys = make_shared<vector<string>>();
302 _dim = make_shared<vector<size_t>>();
303 _keys->resize(n);
304 size_t index = 0;
305 string key;
306 for (int i = 0; i < n; i++){
307 key = pairs._keys.at(index)->_name;
308 (*_keys_map)[key]= index;
309 (*_keys)[index++] = key;
310 }
311 _dim->resize(1);
312 _dim->at(0) = n;
313 _dim->resize(2,0);
314 }
315
316
317 /** Returns the number of comma-separated fields in each key */
get_nb_entries()318 unsigned get_nb_entries() const{
319 if(_keys->empty()){
320 return 0;
321 }
322 return count(_keys->front().begin(), _keys->front().end(), ',') + 1;
323 }
324
325 /** Retain the entries starting at the ith position
326 @param[in] start_position Use the substring starting after the start_position comma separator
327 @param[in] nb_entries Number of comma separated entries to keep
328 */
from_ith(unsigned start_position,unsigned nb_entries)329 indices from_ith(unsigned start_position, unsigned nb_entries) const{
330 indices res(this->get_name()+"from_ith("+to_string(start_position)+")");
331 string key;
332 res._ids = make_shared<vector<vector<size_t>>>();
333 res._ids->resize(1);
334 if(_type == matrix_){/* If ids is matrix indexed */
335 res._type = matrix_;
336 auto nb_rows = this->get_nb_rows();
337 res._ids->resize(nb_rows);
338 for (size_t i = 0; i<nb_rows; i++) {
339 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
340 auto key_ref = _ids->at(i).at(j);
341 key = _keys->at(key_ref);
342 auto pos = nthOccurrence(key, ",", start_position);
343 if(pos>0){
344 key = key.substr(pos+1);
345 }
346 pos = nthOccurrence(key, ",", nb_entries);
347 if(pos>0){
348 key = key.substr(0,pos);
349 }
350 auto it = res._keys_map->find(key);
351 if (it == res._keys_map->end()){
352 res.insert(key);
353 res._ids->at(i).push_back(res._keys->size()-1);
354 }
355 else{
356 res._ids->at(i).push_back(it->second);
357 }
358 }
359 }
360 }
361 else if(is_indexed()){/* If ids has key references, use those */
362 for(auto &key_ref: _ids->at(0)){
363 key = _keys->at(key_ref);
364 auto pos = nthOccurrence(key, ",", start_position);
365 if(pos>0){
366 key = key.substr(pos+1);
367 }
368 pos = nthOccurrence(key, ",", nb_entries);
369 if(pos>0){
370 key = key.substr(0,pos);
371 }
372 auto it = res._keys_map->find(key);
373 if (it == res._keys_map->end()){
374 res.insert(key);
375 res._ids->at(0).push_back(res._keys->size()-1);
376 }
377 else{
378 res._ids->at(0).push_back(it->second);
379 }
380 }
381
382 }
383 else {
384 for(auto key: *_keys){
385 auto pos = nthOccurrence(key, ",", start_position);
386 if(pos>0){
387 key = key.substr(pos+1);
388 }
389 pos = nthOccurrence(key, ",", nb_entries);
390 if(pos>0){
391 key = key.substr(0,pos);
392 }
393 auto it = res._keys_map->find(key);
394 if (it == res._keys_map->end()){
395 res.insert(key);
396 res._ids->at(0).push_back(res._keys->size()-1);
397 }
398 else{
399 res._ids->at(0).push_back(it->second);
400 }
401 }
402 }
403 return res;
404 }
405
406 /** Remove keys starting at the ith position and spanning nb_entries
407 @param[in] start_position
408 */
ignore_ith(unsigned start_position,unsigned nb_entries)409 indices ignore_ith(unsigned start_position, unsigned nb_entries) const{
410 indices res(this->get_name()+"ignore_ith("+to_string(start_position)+","+to_string(nb_entries)+")");
411 string key;
412 res._ids = make_shared<vector<vector<size_t>>>();
413 res._ids->resize(1);
414 string first_part, last_part;
415 if(_type == matrix_){/* If ids is matrix indexed */
416 res._type = matrix_;
417 auto nb_rows = this->get_nb_rows();
418 res._ids->resize(nb_rows);
419 for (size_t i = 0; i<nb_rows; i++) {
420 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
421 auto key_ref = _ids->at(i).at(j);
422 key = _keys->at(key_ref);
423 auto pos = nthOccurrence(key, ",", start_position);
424 first_part = key.substr(0,pos);
425 if(pos>0){
426 key = key.substr(pos+1);
427 }
428 pos = nthOccurrence(key, ",", nb_entries);
429 if(pos>0){
430 last_part = key.substr(pos+1);
431 }
432 if(first_part.size()>0 && last_part.size()>0){ /* stitch them together */
433 key = first_part+","+last_part;
434 }
435 else {
436 key = first_part+last_part;
437 }
438 auto it = res._keys_map->find(key);
439 if (it == res._keys_map->end()){
440 res.insert(key);
441 res._ids->at(i).push_back(res._keys->size()-1);
442 }
443 else{
444 res._ids->at(i).push_back(it->second);
445 }
446 }
447 }
448 }
449 else if(is_indexed()){/* If current index set has key references, use those */
450 for(auto &key_ref: _ids->at(0)){
451 key = _keys->at(key_ref);
452 auto pos = nthOccurrence(key, ",", start_position);
453 first_part = key.substr(0,pos);
454 if(pos>0){
455 key = key.substr(pos+1);
456 }
457 pos = nthOccurrence(key, ",", nb_entries);
458 if(pos>0){
459 last_part = key.substr(pos+1);
460 }
461 if(first_part.size()>0 && last_part.size()>0){ /* stitch them together */
462 key = first_part+","+last_part;
463 }
464 else {
465 key = first_part+last_part;
466 }
467 auto it = res._keys_map->find(key);
468 if (it == res._keys_map->end()){
469 res.insert(key);
470 res._ids->at(0).push_back(res._keys->size()-1);
471 }
472 else{
473 res._ids->at(0).push_back(it->second);
474 }
475 }
476
477 }
478 else {
479 string first_part, last_part;
480 for(auto key: *_keys){
481 auto pos = nthOccurrence(key, ",", start_position);
482 first_part = key.substr(0,pos);
483 if(pos>0){
484 key = key.substr(pos+1);
485 }
486 pos = nthOccurrence(key, ",", nb_entries);
487 if(pos>0){
488 last_part = key.substr(pos+1);
489 }
490 if(first_part.size()>0 && last_part.size()>0){ /* stitch them together */
491 key = first_part+","+last_part;
492 }
493 else {
494 key = first_part+last_part;
495 }
496 auto it = res._keys_map->find(key);
497 if (it == res._keys_map->end()){
498 res.insert(key);
499 res._ids->at(0).push_back(res._keys->size()-1);
500 }
501 else{
502 res._ids->at(0).push_back(it->second);
503 }
504 }
505 }
506 return res;
507 }
508
509
510
511 /** Returns a vector of bools indicating if the ith reference is in ids and in this. The function iterates over key references in _ids. */
get_common_refs(const indices & ids)512 vector<bool> get_common_refs(const indices& ids) const{
513 vector<bool> res;
514 // assert(_ids);
515 set<size_t> unique_ids;
516 if(_type == matrix_){/* If ids is matrix indexed */
517 auto nb_rows = this->get_nb_rows();
518 for (size_t i = 0; i<nb_rows; i++) {
519 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
520 auto idx = _ids->at(i).at(j);
521 if(ids.has_key(_keys->at(idx))){
522 res.push_back(true);
523 }
524 else {
525 res.push_back(false);
526 }
527 }
528 }
529 }
530 else if(is_indexed()){/* If ids has key references, use those */
531 for (auto &idx:_ids->at(0)) {
532 if(ids.has_key(_keys->at(idx))){
533 res.push_back(true);
534 }
535 else {
536 res.push_back(false);
537 }
538 }
539 }
540 else {
541 for (auto &key:*_keys) {
542 if(ids.has_key(key)){
543 res.push_back(true);
544 }
545 else {
546 res.push_back(false);
547 }
548 }
549 }
550 return res;
551 }
552
has_key(const string & key)553 bool has_key(const string& key) const{
554 if(is_indexed()){/* If ids has key references, use those */
555 auto idx0 = _keys_map->at(key);
556 for (auto i= 0; i<_ids->size();i++) {
557 for (auto const idx :_ids->at(i)) {
558 if(idx0==idx)
559 return true;
560 }
561 }
562 return false;
563 }
564 return _keys_map->count(key)!=0;
565 }
566
567 /** Returns a vector of bools indicating if the ith reference is in this but not in ids. The function iterates over key references in _ids. */
get_diff_refs(const indices & ids)568 vector<bool> get_diff_refs(const indices& ids) const{
569 vector<bool> res;
570 // assert(_ids);
571 set<size_t> unique_ids;
572 if(_type == matrix_){/* If ids is matrix indexed */
573 auto nb_rows = this->get_nb_rows();
574 for (size_t i = 0; i<nb_rows; i++) {
575 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
576 auto idx = _ids->at(i).at(j);
577 if(!ids.has_key(_keys->at(idx))){
578 res.push_back(true);
579 }
580 else {
581 res.push_back(false);
582 }
583 }
584 }
585 }
586 else if(is_indexed()){/* If ids has key references, use those */
587 for (auto &idx:_ids->at(0)) {
588 if(!ids.has_key(_keys->at(idx))){
589 res.push_back(true);
590 }
591 else {
592 res.push_back(false);
593 }
594 }
595 }
596 else {
597 for (auto &key:*_keys) {
598 if(!ids.has_key(key)){
599 res.push_back(true);
600 }
601 else {
602 res.push_back(false);
603 }
604 }
605 }
606 return res;
607 }
608
609 /** Returns an index set based on the references stored in _ids. The function iterates over key references in _ids and keeps only the unique entries */
get_unique_keys()610 indices get_unique_keys() const{
611 indices res(_name);
612 res._type = unindexed_;
613 res._keys = _keys;
614 res._keys_map = _keys_map;
615 res._excluded_keys = _excluded_keys;
616 res._dim = make_shared<vector<size_t>>();
617 res._dim->resize(1);
618 res._dim->at(0) = res._keys->size();
619 return res;
620 }
621
622 /* Delete rows where keep[i] is false. */
filter_rows(const vector<bool> & keep)623 void filter_rows(const vector<bool>& keep){
624 string excluded = "\\{";
625 if(_type == matrix_){/* If ids is matrix indexed */
626 if(keep.size()!=get_nb_rows()){
627 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
628 }
629 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
630 auto nb_rows = this->get_nb_rows();
631 for (size_t i = 0; i<nb_rows; i++) {
632 if(keep[i]==true){
633 new_ids->push_back(_ids->at(i));
634 }
635 }
636 _ids = new_ids;
637 }
638 else if(_ids){
639 if(keep.size()!=_ids->at(0).size()){
640 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
641 }
642 vector<vector<size_t>> new_ids;
643 new_ids.resize(1);
644 for (auto idx = 0; idx<keep.size();idx++) {
645 if(keep[idx]){
646 new_ids.at(0).push_back(_ids->at(0).at(idx));
647 }
648 else {
649 excluded += "("+_keys->at(_ids->at(0).at(idx)) +"),";
650 }
651 }
652 *_ids = new_ids;
653 _name += excluded.substr(0,excluded.size()-1) + "}";
654 }
655 else {
656 if(keep.size()!=_keys->size()){
657 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
658 }
659 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
660 new_ids->resize(1);
661 for (auto idx = 0; idx<keep.size();idx++) {
662 if(keep[idx]){
663 new_ids->at(0).push_back(idx);
664 }
665 else {
666 excluded += "("+_keys->at(idx) +"),";
667 }
668 }
669 _ids = new_ids;
670 _name += excluded.substr(0,excluded.size()-1) + "}";
671 }
672 }
673
674 /* Delete row if all correspondong cols in keep[i] are false. */
filter_matrix_rows(const vector<bool> & keep)675 void filter_matrix_rows(const vector<bool>& keep){
676 string excluded = "\\{";
677 if(_type == matrix_){
678 auto nb_rows = _ids->size();
679 auto nb_cols = keep.size()/nb_rows;
680 vector<vector<size_t>> new_ids;
681 new_ids.resize(nb_rows);
682 for (auto row = 0; row<nb_rows;row++) {
683 bool exclude_row = true;/* Excluderow if all cols are false */
684 for (auto col = 0; col<nb_cols;col++) {
685 if(keep[row*nb_cols+col]){
686 exclude_row = false;
687 }
688 }
689 if(exclude_row){
690 excluded += "("+_keys->at(_ids->at(row).at(0)) +"),";
691 }
692 else {
693 new_ids.at(row).push_back(_ids->at(row).at(0));
694 }
695 }
696 *_ids = new_ids;
697 _name += excluded.substr(0,excluded.size()-1) + "}";
698 remove_empty_rows();
699 }
700 else if(_ids){
701 if(keep.size()<=_ids->at(0).size()){
702 throw invalid_argument("in filter_matrix_rows(const vector<bool>& keep): dimension mismatch");
703 }
704 auto nb_rows = _ids->at(0).size();
705 auto nb_cols = keep.size()/nb_rows;
706 vector<vector<size_t>> new_ids;
707 new_ids.resize(1);
708 for (auto row = 0; row<nb_rows;row++) {
709 bool exclude_row = true;/* Excluderow if all cols are false */
710 for (auto col = 0; col<nb_cols;col++) {
711 if(keep[row*nb_cols+col]){
712 exclude_row = false;
713 }
714 }
715 if(exclude_row){
716 excluded += "("+_keys->at(_ids->at(0).at(row)) +"),";
717 }
718 else {
719 new_ids.at(0).push_back(_ids->at(0).at(row));
720 }
721 }
722 *_ids = new_ids;
723 _name += excluded.substr(0,excluded.size()-1) + "}";
724 }
725 else {
726 if(keep.size()<=_keys->size()){
727 throw invalid_argument("in filter_matrix_rows(const vector<bool>& keep): dimension mismatch");
728 }
729 auto nb_rows = _keys->size();
730 auto nb_cols = keep.size()/nb_rows;
731 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
732 new_ids->resize(1);
733 for (auto row = 0; row<nb_rows;row++) {
734 bool exclude_row = true;/* Excluderow if all cols are false */
735 for (auto col = 0; col<nb_cols;col++) {
736 if(keep[row*nb_cols+col]){
737 exclude_row = false;
738 }
739 }
740 if(exclude_row){
741 excluded += "("+_keys->at(row) +"),";
742 }
743 else {
744 new_ids->at(0).push_back(row);
745 }
746 }
747 _ids = new_ids;
748 _name += excluded.substr(0,excluded.size()-1) + "}";
749 }
750 }
751
752 /* Delete indices where keep[i] is false. */
filter_refs(const vector<bool> & keep)753 void filter_refs(const vector<bool>& keep){
754 string excluded = "\\{";
755 if(_type == matrix_){/* If ids is matrix indexed */
756 if(keep.size()!=nb_keys()){
757 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
758 }
759 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
760 auto nb_rows = this->get_nb_rows();
761 new_ids->resize(nb_rows);
762 size_t idx = 0;
763 for (size_t i = 0; i<nb_rows; i++) {
764 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
765 if(keep[idx++]==true){
766 new_ids->at(i).push_back(_ids->at(i).at(j));
767 }
768 }
769 }
770 // shared_ptr<vector<vector<size_t>>> nnz_ids = make_shared<vector<vector<size_t>>>();
771 // for (size_t i = 0; i<nb_rows; i++) {
772 // if(new_ids->at(i).size()>0){
773 // nnz_ids->push_back(new_ids->at(i));
774 // }
775 // }
776 _ids = new_ids;
777 remove_empty_rows();
778 }
779 else if(_ids){
780 if(keep.size()!=_ids->at(0).size()){
781 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
782 }
783 vector<vector<size_t>> new_ids;
784 new_ids.resize(1);
785 for (auto idx = 0; idx<keep.size();idx++) {
786 if(keep[idx]){
787 new_ids.at(0).push_back(_ids->at(0).at(idx));
788 }
789 else {
790 excluded += "("+_keys->at(_ids->at(0).at(idx)) +"),";
791 }
792 }
793 *_ids = new_ids;
794 _name += excluded.substr(0,excluded.size()-1) + "}";
795 }
796 else {
797 if(keep.size()!=_keys->size()){
798 throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
799 }
800 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
801 new_ids->resize(1);
802 for (auto idx = 0; idx<keep.size();idx++) {
803 if(keep[idx]){
804 new_ids->at(0).push_back(idx);
805 }
806 else {
807 excluded += "("+_keys->at(idx) +"),";
808 }
809 }
810 _ids = new_ids;
811 _name += excluded.substr(0,excluded.size()-1) + "}";
812 }
813 }
814
is_matrix_indexed()815 bool is_matrix_indexed() const{
816 return _type == matrix_;
817 }
818
remove_empty_rows()819 void remove_empty_rows() {
820 if(_type!=matrix_){
821 throw invalid_argument("clean_empty_rows() can only be called on a matrix indexed set");
822 }
823 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
824 auto nb_rows = this->get_nb_rows();
825 for (size_t i = 0; i<nb_rows; i++) {
826 if(_ids->at(i).size()>0){
827 new_ids->push_back(move(_ids->at(i)));
828 }
829 }
830 _ids = new_ids;
831 }
832
833 /** Returns a vector of bools indicating if a reference is unique so far. The function iterates over key references in _ids. */
get_unique_refs()834 vector<bool> get_unique_refs() const{
835 vector<bool> res;
836 set<size_t> unique_ids;
837 if(_ids){
838 for (auto &idx:_ids->at(0)) {
839 if(unique_ids.insert(idx).second){
840 res.push_back(true);
841 }
842 else {
843 res.push_back(false);
844 }
845 }
846 }
847 return res;
848 }
849 /* transform a matrix indexed set to a vector indexed one */
flatten()850 void flatten() {
851 if(_type != matrix_){
852 return;
853 }
854 shared_ptr<vector<vector<size_t>>> new_ids = make_shared<vector<vector<size_t>>>();
855 new_ids->resize(1);
856 for (size_t i = 0; i<get_nb_rows(); i++) {
857 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
858 new_ids->at(0).push_back(_ids->at(i).at(j));
859 }
860 }
861 _ids = new_ids;
862 _type = in_;
863 }
864 // /** The function iterates over the ith key references in _ids and deletes the ones where keep[i] is false. */
865 // void filter_refs(const vector<bool>& keep) const{
866 // if(_ids){
867 // if(keep.size()!=_ids->at(0).size()){
868 // throw invalid_argument("in filter_refs(const vector<bool>& keep): keep has a different size than index set");
869 // }
870 // vector<vector<size_t>> new_ids;
871 // new_ids.resize(1);
872 // for (auto idx = 0; idx<keep.size();idx++) {
873 // if(keep[idx]){
874 // new_ids.at(0).push_back(_ids->at(0).at(idx));
875 // }
876 // }
877 // *_ids = new_ids;
878 // }
879 // }
880
881
882
883 /* Returns true if current index set is a subset of ids */
is_subset(const indices & ids)884 bool is_subset(const indices & ids) const{
885 if(_type == matrix_){/* If ids is matrix indexed */
886 auto nb_rows = this->get_nb_rows();
887 for (size_t i = 0; i<nb_rows; i++) {
888 for (size_t j = 0; j<this->_ids->at(i).size(); j++) {
889 auto key_ref = _ids->at(i).at(j);
890 auto key = _keys->at(key_ref);
891 if(!ids.has_key(key)){
892 return false;
893 }
894 }
895 }
896 if(ids._type==matrix_){
897 if(nb_keys()==ids.nb_keys())
898 return false;
899 }
900 else {
901 if(nb_keys()==ids.size())
902 return false;
903 }
904 }
905 else if(_ids){
906 for (auto idx = 0; idx<size();idx++) {
907 if(!ids.has_key(_keys->at(_ids->at(0).at(idx)))){
908 return false;
909 }
910 }
911 }
912 else {
913 for (auto idx = 0; idx<size();idx++) {
914 if(!ids.has_key(_keys->at(idx))){
915 return false;
916 }
917 }
918 }
919 if(ids._type==matrix_){
920 if(size()==ids.nb_keys())
921 return false;
922 }
923 else {
924 if(size()==ids.size())
925 return false;
926 }
927 return true;
928 }
929
930 /* Returns true if current index set is a superset of ids */
is_superset(const indices & ids)931 bool is_superset(const indices & ids) const{
932 return ids.is_subset(*this);
933 }
934
935 /* Transform the current index set into a matrix indexed set. If transformed (it was not a matrix already), the resulting set will have unidimentional columns. */
to_matrix()936 void to_matrix(){
937 if(_type!=matrix_){
938 *this = get_matrix_form();
939 }
940 }
941
942 /* Get a matrix version of the current index set. If transformed (it was not a matrix already), the resulting set will have unidimentional columns. */
get_matrix_form()943 indices get_matrix_form() const {
944 if(_type==matrix_){
945 return *this;
946 }
947 indices res(_name);
948 res._type = matrix_;
949 res._keys = _keys;
950 res._keys_map= _keys_map;
951 auto nb_rows = this->size();
952 for (size_t i = 0; i<nb_rows; i++) {
953 res.add_empty_row();
954 if(_ids){
955 res.add_in_row(i, _keys->at(_ids->at(0).at(i)));
956 }
957 else {
958 res.add_in_row(i, _keys->at(i));
959 }
960 }
961 return res;
962 }
963 /** The function iterates over key references in _ids and keeps only the unique entries */
keep_unique_keys()964 void keep_unique_keys(){
965 indices res(_name);
966 set<size_t> unique_ids;
967 if(_ids){
968 for (auto &idx:_ids->at(0)) {
969 if(unique_ids.insert(idx).second){
970 res.add(_keys->at(idx));
971 }
972 }
973 *this=res;
974 }
975 }
976 // indices(const space& s){
977 // list<indices> l;
978 // _dim->resize(l.size());
979 // for (auto i = 0; i<s._dim.size(); i++) {
980 // l.push_back(indices(0,s._dim[i]-1));
981 // _dim->at(i) = s._dim[i];
982 // }
983 // *this = indices(l);
984 // }
985
986
987
988 template<typename... Args>
init(string idx1,Args &&...args)989 void init(string idx1, Args&&... args) {
990 list<string> indices;
991 indices = {forward<string>(args)...};
992 indices.push_front(idx1);
993 auto n = indices.size();
994 _keys_map = make_shared<map<string,size_t>>();
995 _keys = make_shared<vector<string>>();
996 _keys->resize(n);
997 _dim = make_shared<vector<size_t>>();
998 _dim->resize(1);
999 _dim->at(0) = n;
1000 auto it = indices.begin();
1001 for (size_t i= 0; i< n; i++) {
1002 (*_keys_map)[*it]= i;
1003 (*_keys)[i] = (*it);
1004 it++;
1005 }
1006 }
1007
1008
1009 bool operator==(const indices& cpy) const{
1010 if (_name != cpy._name || _type != cpy._type || _time_extended!=cpy._time_extended || _time_pos != cpy._time_pos || *_dim!=*cpy._dim || _excluded_keys != cpy._excluded_keys || *_keys_map != *cpy._keys_map) return false;
1011 if(_ids==cpy._ids) return true; /* accounts for both being nullptr */
1012 if((_ids && !cpy._ids) || (cpy._ids && !_ids) || (*_ids != *cpy._ids)) return false;
1013 return true;
1014 }
1015
1016 bool operator!=(const indices& cpy) const{
1017 return !(*this==cpy);
1018 }
1019
shallow_copy(shared_ptr<indices> cpy)1020 void shallow_copy(shared_ptr<indices> cpy){
1021 _name = cpy->_name;
1022 _type = cpy->_type;
1023 _keys_map = cpy->_keys_map;
1024 _keys = cpy->_keys;
1025 _dim = cpy->_dim;
1026 _excluded_keys = cpy->_excluded_keys;
1027 if(cpy->_ids){
1028 _ids = make_shared<vector<vector<size_t>>>(*cpy->_ids);
1029 }
1030 _time_extended = cpy->_time_extended;
1031 _time_pos = cpy->_time_pos;
1032 }
1033
time_expand(const indices & T)1034 void time_expand(const indices& T) {//Fix this to expand ids size
1035 _time_extended = true;
1036 auto dim = this->size()*(T.size() - T._excluded_keys.size());
1037 /* update the indices of the old parameter*/
1038 string key;
1039 auto keys = *_keys;
1040 //CLEAR OLD ENTRIES
1041 _keys->clear();
1042 _keys_map->clear();
1043 _keys->resize(dim);
1044 _ids->at(0).clear();
1045 //STORE NEW ENTRIES
1046 unsigned index = 0;
1047 for (auto i = 0; i<keys.size();i++) {
1048 for(unsigned t = 0; t < T.size(); t++ ) {
1049 if(T._excluded_keys.count(t)!=0){
1050 continue;
1051 }
1052 key = keys[i];
1053 key += ",";
1054 key += T._keys->at(t);
1055 _keys_map->insert(make_pair<>(key, index));
1056 _keys->at(index++) = key;
1057 }
1058 }
1059 _name += ".time_expanded";
1060 }
1061
1062 indices& operator=(const indices& cpy){
1063 if(_name.empty())
1064 _name = cpy._name;
1065 _type = cpy._type;
1066 _keys_map = cpy._keys_map;
1067 _keys = cpy._keys;
1068 _dim = cpy._dim;
1069 _excluded_keys = cpy._excluded_keys;
1070 if(cpy._ids){
1071 _ids = make_shared<vector<vector<size_t>>>(*cpy._ids);
1072 }
1073 else{
1074 _ids = nullptr;
1075 }
1076 _time_extended = cpy._time_extended;
1077 _time_pos = cpy._time_pos;
1078 return *this;
1079 }
1080
deep_copy()1081 indices deep_copy() const{
1082 indices cpy;
1083 cpy._name = _name;
1084 cpy._type = _type;
1085 cpy._dim = _dim;
1086 cpy._excluded_keys = _excluded_keys;
1087 if(_ids){
1088 cpy._ids = make_shared<vector<vector<size_t>>>(*_ids);
1089 }
1090 if(_keys){
1091 cpy._keys = make_shared<vector<string>>(*_keys);
1092 }
1093 if(_keys_map){
1094 cpy._keys_map = make_shared<map<string,size_t>>(*_keys_map);
1095 }
1096 cpy._time_extended = _time_extended;
1097 cpy._time_pos = _time_pos;
1098 return cpy;
1099 }
1100
1101 indices& operator=(indices&& cpy){
1102 if(_name.empty())
1103 _name = cpy._name;
1104 _type = cpy._type;
1105 _keys_map = move(cpy._keys_map);
1106 _excluded_keys = move(cpy._excluded_keys);
1107 _keys = move(cpy._keys);
1108 _dim = move(cpy._dim);
1109 _ids = move(cpy._ids);
1110 _time_extended = cpy._time_extended;
1111 _time_pos = cpy._time_pos;
1112 return *this;
1113 }
1114
indices(const indices & cpy)1115 indices(const indices& cpy){
1116 *this=cpy;
1117 }
1118
indices(indices && cpy)1119 indices(indices&& cpy){
1120 *this=move(cpy);
1121 }
1122
1123 // template<typename Tobj>
1124 // indices(const vector<Tobj*>& vec){
1125 // _keys_map = make_shared<map<string,size_t>>();
1126 // _keys = make_shared<vector<string>>();
1127 // size_t i = 0;
1128 // for (auto idx:vec) {
1129 // if(idx->_active){
1130 // _keys->push_back(idx->_name);
1131 // (*_keys_map)[idx->_name]= i;
1132 // i++;
1133 // }
1134 // }
1135 // }
1136 //
1137 // template<typename Tobj>
1138 // indices(const vector<Tobj>& vec){
1139 // _keys_map = make_shared<map<string,size_t>>();
1140 // _keys = make_shared<vector<string>>();
1141 // size_t i = 0;
1142 // for (auto idx:vec) {
1143 // if(idx._active){
1144 // _keys->push_back(idx._name);
1145 // (*_keys_map)[idx._name]= i;
1146 // i++;
1147 // }
1148 // }
1149 // }
1150
1151 template<typename Tobj>
1152 indices(const vector<Tobj*>& vec, bool include_inactive = false){
1153 _keys_map = make_shared<map<string,size_t>>();
1154 _keys = make_shared<vector<string>>();
1155 _dim = make_shared<vector<size_t>>();
1156 _dim->resize(1);
1157 size_t i = 0;
1158 for (auto idx:vec) {
1159 if(include_inactive || idx->_active){
1160 (*_keys_map)[idx->_name]= i;
1161 _keys->push_back(idx->_name);
1162 i++;
1163 }
1164 }
1165 if (_keys->size()>0) {
1166 _name = vec.front()->_type_name;
1167 }
1168 _dim->at(0) = _keys->size();
1169 }
1170
1171 template<typename Tobj>
1172 indices(const vector<Tobj>& vec, bool include_inactive = false){
1173 _keys_map = make_shared<map<string,size_t>>();
1174 _keys = make_shared<vector<string>>();
1175 _dim = make_shared<vector<size_t>>();
1176 _dim->resize(1);
1177 size_t i = 0;
1178 for (auto idx:vec) {
1179 if(include_inactive || idx._active){
1180 (*_keys_map)[idx._name]= i;
1181 _keys->push_back(idx._name);
1182 i++;
1183 }
1184 }
1185 if (_keys->size()>0) {
1186 _name = vec.front()->_type_name;
1187 }
1188 _dim->at(0) = _keys->size();
1189 }
1190
insert(const string & key)1191 void insert(const string& key){
1192 _keys->push_back(key);
1193 _keys_map->insert(make_pair<>(key,_keys->size()-1));
1194 }
1195
extend(const indices & T)1196 void extend(const indices& T) {
1197 if(!_ids){
1198 *this = indices(*this,T);
1199 return;
1200 }
1201 if (_ids->size()>1) {//double indexed
1202 auto dim = _ids->size();
1203 auto new_dim = T.size()*dim;
1204 _ids->resize(new_dim);
1205 for (auto i =dim; i<new_dim; i++) {
1206 _ids->at(i).resize(_ids->at(i%dim).size());
1207 for (auto j =0; j<_ids->at(i).size(); j++) {
1208 _ids->at(i).at(j) = _ids->at(i%dim).at(j);
1209 }
1210 }
1211 }
1212 else {
1213 auto dim = _ids->at(0).size();
1214 auto new_dim = T.size()*dim;
1215 _ids->at(0).resize(new_dim);
1216 for (auto i = dim; i<new_dim; i++) {
1217 _ids->at(0).at(i) = _ids->at(0).at(i%dim);
1218 }
1219 }
1220 }
1221
indices(const list<indices> & vecs)1222 indices(const list<indices>& vecs) {
1223 // if (vecs.size()==2) {
1224 // _type = matrix_;
1225 // }
1226 _keys_map = make_shared<map<string,size_t>>();
1227 _keys = make_shared<vector<string>>();
1228 _dim = make_shared<vector<size_t>>();
1229 _dim->resize(vecs.size());
1230 size_t dim = 1;
1231 size_t time_pos= 0, nb_ids = 0, idx = 0;
1232 vector<size_t> dims;
1233 for(auto &vec: vecs){
1234 if(vec.empty()){
1235 WarningOff("\n\nWARNING: Defining indices with an empty vector!\n\n");
1236 // exit(-1);
1237 }
1238 if(vec._time_extended){
1239 _time_extended = true;
1240 time_pos = vec._time_pos;
1241 }
1242 else{
1243 nb_ids++;
1244 }
1245 dim *= vec.size();
1246 dims.push_back(vec.size());
1247 _name += vec._name+",";
1248 _dim->at(idx++) = vec.size();
1249 }
1250 auto vec1 = vecs.front();
1251 _name = _name.substr(0,_name.size()-1); /* remove last comma */
1252 if(_time_extended && !vec1._time_extended){
1253 if(time_pos==0 && !vec1.empty()) {//TODO CHECK
1254 _time_pos = std::count(vec1._keys->front().begin(), vec1._keys->front().end(), ',')+1;
1255 }
1256 else {
1257 _time_pos = time_pos+nb_ids;
1258 }
1259 }
1260 _keys->resize(dim);
1261 size_t den = 1;
1262 size_t real_idx = 0;
1263 bool excluded = false;
1264 for(size_t idx = 0; idx < dim ; idx++){
1265 string key;
1266 den = dim;
1267 excluded = false;
1268 for(auto it = vecs.begin(); it!= vecs.end(); it++) {
1269 auto vec = &(*it);
1270 den /= vec->size();
1271 real_idx = (idx/den)%vec->size();
1272 if (vec->_excluded_keys.count(real_idx)==1) {
1273 excluded = true;
1274 }
1275 if(vec->is_indexed()){
1276 real_idx = vec->_ids->at(0).at(real_idx);
1277 }
1278 key += vec->_keys->at(real_idx);
1279 if(next(it)!=vecs.end()){
1280 key += ",";
1281 }
1282 }
1283 (*_keys)[idx] = key;
1284 (*_keys_map)[key] = idx;
1285 if (excluded) {
1286 _excluded_keys.insert(idx);
1287 }
1288 }
1289 }
1290
1291
1292
1293
1294
1295
1296 template<typename... Args>
insert(const string & s1,Args &&...args)1297 void insert(const string& s1, Args&&... args) {
1298 add(s1,args...);
1299 }
1300
insert(const vector<string> & all_keys)1301 void insert(const vector<string>& all_keys) {
1302 add(all_keys);
1303 }
1304
add_empty_row()1305 void add_empty_row() {
1306 _type = matrix_;
1307 if (!_ids) {
1308 _ids = make_shared<vector<vector<size_t>>>();
1309 }
1310 _ids->resize(_ids->size()+1);
1311 }
1312
1313 /** Add a key in the specified row (for sparse matrix indexing)*/
add_in_row(size_t row_nb,const string & key)1314 void add_in_row(size_t row_nb, const string& key) {
1315 _type = matrix_;
1316 if (!_ids) {
1317 _ids = make_shared<vector<vector<size_t>>>();
1318 }
1319 _ids->resize(row_nb+1);
1320 auto ref_it = _keys_map->find(key);
1321 if(ref_it==_keys_map->end()){
1322 auto idx = _keys->size();
1323 _keys_map->insert(make_pair<>(key,idx));
1324 _keys->push_back(key);
1325 _ids->at(row_nb).push_back(idx);
1326 }
1327 else{
1328 _ids->at(row_nb).push_back(ref_it->second);
1329 }
1330 }
1331
1332 template<typename... Args>
add(const string & s1,Args &&...args)1333 void add(const string& s1, Args&&... args) {
1334 add(vector<string>({s1,args...}));
1335 }
1336
1337
1338 template<typename... Args>
indices(const indices & vec1,Args &&...args)1339 indices(const indices& vec1, Args&&... args) {
1340 list<indices> vecs;
1341 vecs = {vec1,forward<Args>(args)...};
1342 *this = indices(vecs);
1343 }
1344
1345 /** Adds all new keys found in ids
1346 @return index set of added indices.
1347 */
1348 template<typename... Args>
add(const indices & ids)1349 indices add(const indices& ids) {
1350 indices added("added");
1351 if(ids.is_matrix_indexed()){
1352 if(!is_matrix_indexed()){
1353 throw invalid_argument("calling add(ids) with a matrix indexed set while current set is not matrix indexed.");
1354 }
1355 auto nb_rows = ids.get_nb_rows();
1356 for (size_t i = 0; i<nb_rows; i++) {
1357 for (size_t j = 0; j<ids._ids->at(i).size(); j++) {
1358 auto key = ids._keys->at(ids._ids->at(i).at(j));
1359 auto idx = _keys->size();
1360 auto pp = _keys_map->insert(make_pair<>(key,idx));
1361 if (pp.second) {//new index inserted
1362 _keys->push_back(key);
1363 added.add(key);
1364 }
1365 }
1366 }
1367 }
1368 else if(ids.is_indexed()){
1369 for(auto &key_ref: ids._ids->at(0)){
1370 auto key = ids._keys->at(key_ref);
1371 auto idx = _keys->size();
1372 auto pp = _keys_map->insert(make_pair<>(key,idx));
1373 if (pp.second) {//new index inserted
1374 _keys->push_back(key);
1375 added.add(key);
1376 }
1377 }
1378 }
1379 else {
1380 auto it = ids._keys->begin();
1381 for (size_t i= 0; i < ids.size(); i++) {
1382 auto idx = _keys->size();
1383 auto pp = _keys_map->insert(make_pair<>(*it,idx));
1384 if (pp.second) {//new index inserted
1385 _keys->push_back(*it);
1386 added.add(*it);
1387 }
1388 it++;
1389 }
1390 }
1391 return added;
1392 }
1393
1394
1395 /*
1396 Add refs to all keys found in ids
1397 */
1398 template<typename... Args>
add_refs(const indices & ids)1399 void add_refs(const indices& ids) {
1400 if(ids.is_matrix_indexed()){
1401 if(!is_matrix_indexed()){
1402 throw invalid_argument("calling add_refs(ids) with a matrix indexed set while current set is not matrix indexed.");
1403 }
1404 auto nb_rows = ids.get_nb_rows();
1405 for (size_t i = 0; i<nb_rows; i++) {
1406 _ids->push_back(ids._ids->at(i));
1407 }
1408 }
1409 else if(ids.is_indexed()){
1410 for(auto &key_ref: ids._ids->at(0)){
1411 add_ref(ids._keys->at(key_ref));
1412 }
1413 }
1414 else {
1415 for(auto &key: *ids._keys){
1416 add_ref(key);
1417 }
1418 }
1419 }
1420
1421
1422
1423 /** Adds a reference in _ids
1424 */
add_ref(int key_id)1425 void add_ref(int key_id){
1426 if(!_ids){
1427 _ids = make_shared<vector<vector<size_t>>>();
1428 _ids->resize(1);
1429 }
1430 _ids->at(0).push_back(key_id);
1431 }
1432
1433
1434 /** Adds a reference to the key specified as argument, i.e., adds the corresponding index in _ids
1435 @throw invalid_argument if key is not part of _keys_map
1436 */
add_ref(const string & key)1437 void add_ref(const string& key){
1438 if(!_ids){
1439 _ids = make_shared<vector<vector<size_t>>>();
1440 _ids->resize(1);
1441 }
1442 auto ref_it = _keys_map->find(key);
1443 if(ref_it==_keys_map->end()){
1444 throw invalid_argument("in indices::add_ref(string), unknown key: " + key);
1445 }
1446 _ids->at(0).push_back(ref_it->second);
1447 }
1448
add(const vector<string> & all_keys)1449 void add(const vector<string>& all_keys) {
1450 for (auto &key: all_keys) {
1451 auto idx = _keys->size();
1452 auto pp = _keys_map->insert(make_pair<>(key,idx));
1453 if (pp.second) {//new key inserted
1454 _keys->push_back(key);
1455 }
1456 else{
1457 throw invalid_argument("in indices::add(string...) cannot add same key twice: " + key);
1458 }
1459 }
1460 }
1461
1462
1463 inline size_t get_id_inst(size_t inst = 0) const {
1464 if (_ids) {
1465 if(_ids->at(0).size() <= inst){
1466 throw invalid_argument("indices::get_id_inst(size_t inst) inst is out of range");
1467 }
1468 return _ids->at(0).at(inst);
1469 }
1470 return inst;
1471 };
1472
1473
1474
1475
1476
is_indexed()1477 bool is_indexed() const{
1478 return (_ids!=nullptr);
1479 }
1480
remove_excluded()1481 void remove_excluded(){
1482 _ids = nullptr;
1483 map<string,size_t> new_keys_map;
1484 for(auto &key_id: _excluded_keys){
1485 auto key = _keys->at(key_id);
1486 _keys_map->erase(key);
1487 }
1488 _keys->clear();
1489 _keys->resize(_keys_map->size());
1490 size_t idx = 0;
1491 for(auto &key_id: *_keys_map){
1492 _keys->at(idx) = key_id.first;
1493 new_keys_map[key_id.first] = idx++;
1494 }
1495 *_keys_map = new_keys_map;
1496 _excluded_keys.clear();
1497 _dim->resize(1);
1498 _dim->at(0) = _keys->size();
1499 }
1500
reindex()1501 void reindex(){
1502 _ids->at(0).clear();
1503 for(auto idx = 0; idx<_keys->size();idx++){
1504 if(_excluded_keys.count(idx)==0){
1505 _ids->at(0).push_back(idx);
1506 }
1507 }
1508 }
1509
exclude(string key)1510 indices exclude(string key){
1511 auto res = *this;
1512 res._keys_map = make_shared<map<string,size_t>>(*_keys_map);
1513 res._keys = make_shared<vector<string>>(*_keys);
1514 res._excluded_keys.insert(res._keys_map->at(key));
1515 // if(!is_indexed()){
1516 // res._ids = make_shared<vector<vector<size_t>>>();
1517 // res._ids->resize(1);
1518 // }
1519 // res.reindex();
1520 res._name += "\\{" + key+"}";
1521 res.remove_excluded();
1522 return res;
1523 }
1524
get_max_nb_columns()1525 size_t get_max_nb_columns() const {
1526 assert(is_indexed());
1527 if(_type != matrix_){
1528 throw invalid_argument("cannot call get_nb_cols() on a non-indexed set");
1529 }
1530 auto nb_inst = _ids->size();
1531 size_t max_dim = 0;
1532 for(size_t inst = 0; inst<nb_inst; inst++){
1533 auto nb_idx = _ids->at(inst).size();
1534 max_dim = max(max_dim,nb_idx);
1535 }
1536 return max_dim;
1537 }
1538
get_nb_rows()1539 size_t get_nb_rows() const {
1540 if(_type != matrix_){
1541 throw invalid_argument("cannot call get_nb_rows() on a non-indexed set");
1542 }
1543 return _ids->size();
1544 };
1545
1546
1547 /** @return size of index set, if this is a matrix indexing, it returns the number of rows, if the set is a mask (has _ids) it returns the mask size, otherwize it return the total number of keys. **/
size()1548 size_t size() const {
1549 if(is_indexed()){
1550 if(_type == matrix_){ //Matrix-indexed
1551 return get_nb_rows();
1552 }
1553 return _ids->at(0).size();
1554 }
1555 return _keys->size();
1556 };
1557
nb_keys()1558 size_t nb_keys() const {
1559 if(_type != matrix_){
1560 throw invalid_argument("cannot call nb_keys() on a non-matrix index set");
1561 }
1562 size_t n = 0;
1563 for(auto &vec: *_ids){
1564 n+=vec.size();
1565 }
1566 return n;
1567 };
1568
nb_active_keys()1569 size_t nb_active_keys() const {return _keys->size() - _excluded_keys.size();};
1570
empty()1571 bool empty() const {
1572 return _keys->size() - _excluded_keys.size() == 0;
1573 }
print()1574 void print() const{
1575 cout << endl;
1576 auto i = 0;
1577 for(auto &key:*_keys){
1578 if (_excluded_keys.count(i++)==0) {
1579 cout << key << " ";
1580 }
1581 }
1582 cout << endl;
1583 }
1584
first()1585 string first() const{
1586 return _keys->front();
1587 }
1588
last()1589 string last() const{
1590 return _keys->back();
1591 }
1592 };
1593
1594 /** Adds all new keys found in ids
1595 @return index set of added indices.
1596 **/
1597 template<typename... Args>
union_ids(const indices & ids1,Args &&...args)1598 indices union_ids(const indices& ids1, Args&&... args) {
1599 vector<indices> all_ids;
1600 all_ids = {ids1,forward<Args>(args)...};
1601 indices res("Union(");
1602 auto nb_entries = ids1.get_nb_entries();
1603 if(!ids1.is_indexed()){ //if the index set is not indexed, we assume the rest are not indexed as well
1604 res._keys_map = make_shared<map<string,size_t>>(*ids1._keys_map);
1605 res._keys = make_shared<vector<string>>(*ids1._keys);
1606 for (size_t i= 1; i < all_ids.size(); i++) {
1607 auto ids = all_ids[i];
1608 if(nb_entries!=ids.get_nb_entries()){
1609 throw invalid_argument("union cannot be applied to index sets with different number of entries");
1610 }
1611 res.set_name(res.get_name() + ids.get_name()+",");
1612 auto it = ids._keys->begin();
1613 for (size_t i= 0; i < ids.size(); i++) {
1614 auto kkey = res._keys->size();
1615 auto pp = res._keys_map->insert(make_pair<>(*it,kkey));
1616 if (pp.second) {//new index inserted
1617 res._keys->push_back(*it);
1618 }
1619 it++;
1620 }
1621 }
1622 }
1623 else{ //means the ids1 is indexed so we assume all the rest are indexed as well
1624 // res._ids = make_shared<vector<vector<size_t>>>(*ids1._ids);
1625 for (size_t i= 0; i < all_ids.size(); i++) {
1626 auto ids = all_ids[i];
1627 if(nb_entries!=ids.get_nb_entries()){
1628 throw invalid_argument("union cannot be applied to index sets with different number of entries");
1629 }
1630 res.set_name(res.get_name() + ids.get_name()+",");
1631 auto it_ids = ids._ids->begin()->begin();
1632 for (size_t i= 0; i < ids.size(); i++) {
1633 auto kkey = res._keys->size();
1634 auto pp = res._keys_map->insert(make_pair<>(ids._keys->at(*it_ids),kkey));
1635 if (pp.second) {//new index inserted
1636 res._keys->push_back(ids._keys->at(*it_ids));
1637 }
1638 it_ids++;
1639 }
1640 }
1641 }
1642 auto name = res.get_name();
1643 res.set_name(name.substr(0,name.size()-1) + ")");
1644 return res;
1645 }
1646
1647 indices operator-(const indices& s1, const indices& s2);
1648
1649 class node_pairs{
1650
1651 public:
1652 string _name;
1653 std::vector<gravity::index_pair*> _keys;
node_pairs()1654 node_pairs(){
1655 _keys.resize(0);
1656 };
node_pairs(string name)1657 node_pairs(string name):_name(name){_keys.resize(0);};
~node_pairs()1658 ~node_pairs(){
1659 clear();
1660 }
clear()1661 void clear() {
1662 for (auto p: _keys) { delete p;}
1663 _keys.clear();
1664 }
1665 };
1666
1667 typedef enum { ipopt, gurobi, bonmin, cplex, sdpa, _mosek, clp} SolverType; /* Solver type */
1668
1669 // settings of solvers. used by solvers like sdpa.
1670 typedef enum {unsolved = -1, penalty=0, fast=1, medium=2, stable=3} SolverSettings;
1671
1672 typedef pair<shared_ptr<size_t>,shared_ptr<indices>> unique_id; /* A unique identifier is defined as a pair<variable id, index address> */
1673
1674 template <class T>
type_name(const T & t)1675 std::string type_name(const T& t) {
1676 return t._type_name;
1677 }
1678
1679 indices range(size_t i, size_t j);
1680
1681
1682 /** Combine each ith key of each index set with the other ith keys in the other index sets
1683 \warning all index sets must have the same size
1684 */
1685 template<typename... Args>
combine(const indices & ids1,Args &&...args)1686 indices combine(const indices& ids1, Args&&... args){
1687 bool matrix_indexed = false;
1688 indices res;
1689 string res_name;
1690 list<indices> all_ids = {ids1,forward<Args>(args)...};
1691 for (auto &ids: all_ids) {
1692 if(ids.is_matrix_indexed()){
1693 matrix_indexed = true;
1694 }
1695 res_name += ids.get_name()+",";
1696 }
1697 res_name = res_name.substr(0, res_name.size()-1);/* remove last comma */
1698 res.set_name(res_name);
1699 if(matrix_indexed){
1700 auto nb_rows = ids1.get_nb_rows();
1701 res._ids = make_shared<vector<vector<size_t>>>();
1702 res._ids->resize(nb_rows);
1703 for (auto &ids: all_ids) {
1704 if(!ids.is_matrix_indexed())
1705 throw invalid_argument("In combine(ids..) all or none of the index sets should be matrix indexed");
1706 if(ids.get_nb_rows()!=nb_rows)
1707 throw invalid_argument("In combine(ids..) all indices should have the same number of rows");
1708 }
1709 for (size_t i = 0; i<nb_rows; i++) {
1710 for (size_t j = 0; j< ids1._ids->at(i).size(); j++) {
1711 string combined_key, key = "";
1712 for (auto &ids: all_ids) {
1713 if(ids._ids->at(i).size()!=ids1._ids->at(i).size())
1714 throw invalid_argument("In combine(ids..) all indices should have the same number of entries per row");
1715 auto idx = ids._ids->at(i).at(j);
1716 key = ids._keys->at(idx);
1717 combined_key += key +",";
1718 }
1719 combined_key = combined_key.substr(0, combined_key.size()-1);/* remove last comma */
1720 res.add_in_row(i,combined_key);
1721 }
1722 }
1723 }
1724 else {
1725 auto nb_keys = ids1.size();
1726 for (auto i = 0; i< nb_keys; i++) {
1727 string combined_key, key = "", prev_key = "";
1728 // bool same_key = true;
1729 for (auto &ids: all_ids) {
1730 auto idx = ids.get_id_inst(i);
1731 key = ids._keys->at(idx);
1732 // if(prev_key!="" && key!=prev_key){
1733 // same_key = false;
1734 // }
1735 // prev_key = key;
1736 combined_key += key +",";
1737 }
1738 // if(same_key){
1739 // combined_key = key;
1740 // }
1741 // else {
1742 combined_key = combined_key.substr(0, combined_key.size()-1);/* remove last comma */
1743 // }
1744 res.add(combined_key);
1745 }
1746 }
1747 return res;
1748 }
1749
1750
1751 }
1752
1753
1754
1755 #endif
1756