1 //
2 // func.cpp
3 // Gravity
4 //
5 // Created by Hijazi, Hassan on 24/10/16.
6 //
7 //
8 #include <cmath>
9 #include <gravity/func.h>
10 //
11 //
12 using namespace std;
13 namespace gravity{
14
get_mag_ang(const func<Cpx> & f)15 pair<func<double>,func<double>> get_mag_ang(const func<Cpx>& f){
16 func<double> f_mag = get_mag(f._cst.get()), f_ang = get_ang(f._cst.get());
17 for (auto &qt: *f._qterms) {
18 auto coef = qt.second._coef;
19 func<double> p1_mag, p2_mag, p1_ang, p2_ang;
20 func<double> c_mag, c_ang;
21 if(!qt.second._p->first->_polar){
22 auto p1_r = get_real(qt.second._p->first.get());
23 auto p1_i = get_imag(qt.second._p->first.get());
24 p1_mag = sqrt(p1_r*p1_r + p1_i*p1_i);
25 p1_ang = atan2(p1_r,p1_i);// TODO fix atan2
26 }
27 else {
28 p1_mag = get_mag(qt.second._p->first.get());
29 p1_ang = get_ang(qt.second._p->first.get());
30 }
31 if(!qt.second._p->second->_polar){
32 auto p2_r = get_real(qt.second._p->second.get());
33 auto p2_i = get_imag(qt.second._p->second.get());
34 p2_mag = sqrt(p2_r*p2_r + p2_i*p2_i);
35 p2_ang = atan2(p2_r,p2_i);// TODO fix atan2
36 }
37 else {
38 p2_mag = get_mag(qt.second._p->second.get());
39 p2_ang = get_ang(qt.second._p->second.get());
40 }
41 if(!coef->_polar){
42 auto c_r = get_real(coef.get());
43 auto c_i = get_imag(coef.get());
44 c_mag = sqrt(c_r*c_r + c_i*c_i);
45 c_ang = atan2(c_r,c_i);// TODO fix atan2
46 }
47 else {
48 c_mag = get_mag(coef.get());
49 c_ang = get_ang(coef.get());
50 }
51 if(coef->_is_transposed){
52 c_mag.transpose();
53 c_ang.transpose();
54 p1_mag._is_vector = true;
55 p1_ang._is_vector = true;
56 p2_mag._is_vector = true;
57 p2_ang._is_vector = true;
58 }
59 if(qt.second._sign){
60 f_mag += c_mag*p1_mag*p2_mag;
61 f_ang += c_ang + p1_ang + p2_ang;
62 }
63 else {
64 f_mag -= c_mag*p1_mag*p2_mag;
65 f_ang -= c_ang + p1_ang + p2_ang;
66 }
67 }
68 for (auto <: *f._lterms) {
69 auto coef = lt.second._coef;
70 func<double> p1_mag, p1_ang;
71 func<double> c_mag, c_ang;
72 if(!lt.second._p->_polar){
73 auto p1_r = get_real(lt.second._p.get());
74 auto p1_i = get_imag(lt.second._p.get());
75 p1_mag = sqrt(p1_r*p1_r + p1_i*p1_i);
76 p1_ang = atan2(p1_r,p1_i);// TODO fix atan2
77 }
78 else {
79 p1_mag = get_mag(lt.second._p.get());
80 p1_ang = get_ang(lt.second._p.get());
81 }
82 if(!coef->_polar){
83 auto c_r = get_real(coef.get());
84 auto c_i = get_imag(coef.get());
85 c_mag = sqrt(c_r*c_r + c_i*c_i);
86 c_ang = atan2(c_r,c_i);// TODO fix atan2
87 }
88 else {
89 c_mag = get_mag(coef.get());
90 c_ang = get_ang(coef.get());
91 }
92 if(coef->_is_transposed){
93 c_mag.transpose();
94 c_ang.transpose();
95 p1_mag._is_vector = true;
96 p1_ang._is_vector = true;
97 }
98 if(lt.second._sign){
99 f_mag += c_mag * p1_mag;
100 f_ang += c_ang + p1_ang;
101 }
102 else {
103 f_mag -= c_mag * p1_mag;
104 f_ang -= c_ang + p1_ang;
105 }
106 }
107 return {f_mag,f_ang};
108 }
109
110 /** WARNING, only call if the variables appearing in the function are complex or double */
get_real_imag(const func<Cpx> & f)111 pair<func<double>,func<double>> get_real_imag(const func<Cpx>& f){
112 func<double> c_r, c_i, c_mag, c_ang;
113 func<double> p1_r, p2_r, p1_i, p2_i;
114 func<double> p1_mag, p2_mag, p1_ang, p2_ang;
115 func<double> f_r = get_real(f._cst.get()), f_i = get_imag(f._cst.get());
116 for (auto &pt: *f._pterms) {
117 func<> fp_r = 1.;
118 func<> fp_i = 0.;
119 auto coef = pt.second._coef;
120 if(coef->_polar){
121 c_mag = get_mag(coef.get());
122 c_ang = get_ang(coef.get());
123 c_r = c_mag*cos(c_ang);
124 c_i = c_mag*sin(c_ang);
125 }
126 else {
127 c_r = get_real(coef.get());
128 c_i = get_imag(coef.get());
129 }
130 if(coef->_is_transposed){
131 c_r.transpose();
132 c_i.transpose();
133 p1_r._is_vector = true;
134 p1_i._is_vector = true;
135 p2_r._is_vector = true;
136 p2_i._is_vector = true;
137 }
138 auto it = pt.second._l->begin();
139 while(it!=pt.second._l->end()){
140 auto vv = it->first;
141 auto expo = it->second;
142 if(vv->_polar){
143 p1_mag = get_mag(vv.get());
144 p1_ang = get_ang(vv.get());
145 p1_r = p1_mag*cos(p1_ang);
146 p1_i = p1_mag*sin(p1_ang);
147 }
148 else {
149 p1_r = get_real(vv.get());
150 p1_i = get_imag(vv.get());
151 }
152 if(coef->_polar || vv->_polar){
153 throw invalid_argument("unsupported yet");
154 }
155 else {
156 for (auto i = 0; i<expo; i++) {
157 auto fp_r_new = fp_r*p1_r - fp_i*p1_i;
158 auto fp_i_new = fp_r*p1_i + fp_i*p1_r;
159 fp_r = fp_r_new;
160 fp_i = fp_i_new;
161 }
162 }
163 it++;
164 }
165 auto fp_r_new = c_r*fp_r - c_i*fp_i;
166 auto fp_i_new = c_r*fp_i + c_i*fp_r;
167 fp_r = fp_r_new;
168 fp_i = fp_i_new;
169 if(!pt.second._sign){
170 fp_r.reverse_sign();
171 fp_i.reverse_sign();
172 }
173 f_r += fp_r;
174 f_i += fp_i;
175 }
176 for (auto &qt: *f._qterms) {
177 auto coef = qt.second._coef;
178 func<double> p1_r, p2_r, p1_i, p2_i;
179 func<double> p1_mag, p2_mag, p1_ang, p2_ang;
180 func<double> c_r, c_i, c_mag, c_ang;
181 func<> fq_r, fq_i;
182 if(qt.second._p->first->_polar){
183 p1_mag = get_mag(qt.second._p->first.get());
184 p1_ang = get_ang(qt.second._p->first.get());
185 p1_r = p1_mag*cos(p1_ang);
186 p1_i = p1_mag*sin(p1_ang);
187 }
188 else {
189 p1_r = get_real(qt.second._p->first.get());
190 p1_i = get_imag(qt.second._p->first.get());
191 }
192 if(qt.second._p->second->_polar){
193 p2_mag = get_mag(qt.second._p->second.get());
194 p2_ang = get_ang(qt.second._p->second.get());
195 p2_r = p2_mag*cos(p2_ang);
196 p2_i = p2_mag*sin(p2_ang);
197 }
198 else {
199 p2_r = get_real(qt.second._p->second.get());
200 p2_i = get_imag(qt.second._p->second.get());
201 }
202 if(qt.second._p->first->_is_transposed){
203 p1_r.transpose();
204 p1_i.transpose();
205 p2_r._is_vector = true;
206 p2_i._is_vector = true;
207 }
208 if(coef->_polar){
209 c_mag = get_mag(coef.get());
210 c_ang = get_ang(coef.get());
211 c_r = c_mag*cos(c_ang);
212 c_i = c_mag*sin(c_ang);
213 }
214 else {
215 c_r = get_real(coef.get());
216 c_i = get_imag(coef.get());
217 }
218 if(coef->_is_transposed){
219 c_r.transpose();
220 c_i.transpose();
221 p1_r._is_vector = true;
222 p1_i._is_vector = true;
223 p2_r._is_vector = true;
224 p2_i._is_vector = true;
225 }
226 if(coef->_polar && qt.second._p->first->_polar && qt.second._p->second->_polar){
227 fq_r = c_mag*p1_mag*p2_mag*cos(c_ang+p1_ang+p2_ang);
228 fq_i = c_mag*p1_mag*p2_mag*sin(c_ang+p1_ang+p2_ang);
229 }
230 else if(!coef->_polar && qt.second._p->first->_polar && qt.second._p->second->_polar){
231 fq_r = c_r*p1_mag*p2_mag*cos(p1_ang+p2_ang) - c_i*p1_mag*p2_mag*sin(p1_ang+p2_ang);
232 fq_i = c_r*p1_mag*p2_mag*sin(p1_ang+p2_ang) + c_i*p1_mag*p2_mag*cos(p1_ang+p2_ang);
233 }
234 else {
235 fq_r = c_r*(p1_r*p2_r - p1_i*p2_i) - c_i*(p1_r*p2_i + p1_i*p2_r);
236 fq_i = c_r*(p1_r*p2_i + p1_i*p2_r) + c_i*(p1_r*p2_r - p1_i*p2_i);
237 }
238 if(!qt.second._sign){
239 fq_r.reverse_sign();
240 fq_i.reverse_sign();
241 }
242 f_r += fq_r;
243 f_i += fq_i;
244 }
245 for (auto <: *f._lterms) {
246 auto coef = lt.second._coef;
247 func<double> p1_r, p1_i, p1_mag, p1_ang;
248 func<double> c_r, c_i, c_mag, c_ang;
249 func<> fl_r, fl_i;
250 if(lt.second._p->_polar){
251 p1_mag = get_mag(lt.second._p.get());
252 p1_ang = get_ang(lt.second._p.get());
253 p1_r = p1_mag*cos(p1_ang);
254 p1_i = p1_mag*sin(p1_ang);
255 }
256 else {
257 p1_r = get_real(lt.second._p.get());
258 p1_i = get_imag(lt.second._p.get());
259 }
260 if(coef->_polar){
261 c_mag = get_mag(coef.get());
262 c_ang = get_ang(coef.get());
263 c_r = c_mag*cos(c_ang);
264 c_i = c_mag*sin(c_ang);
265 }
266 else {
267 c_r = get_real(coef.get());
268 c_i = get_imag(coef.get());
269 }
270 if(coef->_is_transposed){
271 c_r.transpose();
272 c_i.transpose();
273 p1_r._is_vector = true;
274 p1_i._is_vector = true;
275 }
276 if(coef->_polar && lt.second._p->_polar){
277 fl_r = c_mag*p1_mag*cos(c_ang+p1_ang);
278 fl_i = c_mag*p1_mag*sin(c_ang+p1_ang);
279 }
280 else {
281 fl_r = c_r*p1_r - c_i*p1_i;
282 fl_i = c_r*p1_i + c_i*p1_r;
283 }
284 if(!lt.second._sign){
285 fl_r.reverse_sign();
286 fl_i.reverse_sign();
287 }
288 f_r += fl_r;
289 f_i += fl_i;
290 }
291 return {f_r,f_i};
292 }
293 //
294 // bool is_indexed(const constant_* c){
295 // if (c->is_var() || c->is_param()) {
296 // return (((param_*)c)->is_indexed());
297 // }
298 // return true;
299 // }
300 //
301 // size_t get_id(const constant_* c){
302 // if (c->is_var() || c->is_param()) {
303 // return ((param_*)c)->get_vec_id();
304 // }
305 // return 0;
306 // }
307 //
308 // size_t get_id_inst(const constant_* c, unsigned inst){
309 // if (c->is_var() || c->is_param()) {
310 // return ((param_*)c)->get_id_inst(inst);
311 // }
312 // return 0;
313 // }
314 //
315 //// void func_::eval_vector() {
316 //// // try {
317 //// //// if (_val->size()<_dim[0]) {
318 //// //// _val->resize(_dim[0]);
319 //// //// }
320 //// // if (is_constant() && !_expr && _params->size()==1) { // This is a parameter
321 //// // auto p_c = (*_params->begin()).second.first;
322 //// // for (auto i = 0; i < p_c->get_dim(); i++) {
323 //// // set_val(i,t_eval(p_c.get(),i));
324 //// // }
325 //// // _evaluated = true;
326 //// // return;
327 //// // }
328 //// // }
329 //// // catch (const std::out_of_range& e) {
330 //// // cout << "Out of Range before eval.";
331 //// // }
332 //// // if(_expr && _expr->is_bexpr()){
333 //// // auto be = (bexpr*)_expr.get();
334 //// // if (be->_otype==product_ && _lterms->empty() && _pterms->empty() && _qterms->empty()) {//Pure product
335 //// // if (be->_lson->is_matrix() && be->_rson->_is_vector) {//matrix*vect
336 //// // for (size_t i = 0; i < _dim[0]; i++) {
337 //// // double res = 0;
338 //// // for (size_t j = 0; j < be->_lson->_dim[1]; j++) {
339 //// // res += be->_lson->get_val(i,j) * be->_rson->get_val(j);
340 //// // }
341 //// // set_val(i,be->_coef*res);
342 //// // }
343 //// // return;
344 //// // }
345 //// // }
346 //// // }
347 //// if (_is_transposed) {
348 //// for (auto inst = 0; inst < _dim[1]; inst++) {
349 //// set_val(inst,eval(inst));
350 //// }
351 //// }
352 //// else {
353 //// for (auto inst = 0; inst < _dim[0]; inst++) {
354 //// set_val(inst,eval(inst));
355 //// }
356 //// }
357 //// // }
358 //// // catch (const std::out_of_range& e) {
359 //// // cout << "Out of Range error at inst = " << inst << endl;
360 //// // }
361 ////
362 //// // }
363 //// }
364 //
365 //
366 //// void func_::eval_matrix() {
367 //// // if (is_constant() && !_expr && _params->size()==1) { // This is a parameter
368 //// // auto p_c = (*_params->begin()).second.first;
369 //// // for (size_t i = 0; i < _dim[0]; i++) {
370 //// // for (size_t j = 0; j < _dim[1]; j++) {
371 //// // set_val(i,j,t_eval(p_c.get(), i, j));
372 //// // }
373 //// // }
374 //// // _evaluated = true;
375 //// // return;
376 //// // }
377 //// // double res = 0;
378 //// if(!_expr){
379 //// for (size_t i = 0; i < _dim[0]; i++) {
380 //// for (size_t j = 0; j < _dim[1]; j++) {
381 //// set_val(i,j,eval(i,j));
382 //// }
383 //// }
384 //// }
385 //// if(is_constant()){
386 //// _evaluated = true;
387 //// }
388 //// if(_expr){
389 //// auto be = (bexpr*)_expr.get();
390 //// if (be->_otype==product_ && _lterms->empty() && _pterms->empty() && _qterms->empty()) {//Pure product
391 //// if (be->_lson->is_matrix() && be->_rson->is_matrix()) {
392 //// //matrix product
393 //// if (!_is_hessian) {
394 //// for (size_t i = 0; i < _dim[0]; i++) {
395 //// for (size_t j = 0; j < _dim[1]; j++) {
396 //// auto res = 0.;
397 //// for (size_t col = 0; col<be->_lson->_dim[1]; col++) {
398 //// res += be->_lson->get_val(i,col) * be->_rson->get_val(col,j);
399 //// }
400 //// set_val(i,j, be->_coef*res);
401 //// }
402 //// }
403 //// }
404 //// else {
405 //// for (size_t i = 0; i < _dim[0]; i++) {
406 //// for (size_t j = i; j < _dim[1]; j++) {
407 //// auto res = 0.;
408 //// for (size_t col = 0; col<be->_lson->_dim[1]; col++) {
409 //// res += be->_lson->get_val(i,col) * be->_rson->get_val(col,j);
410 //// }
411 //// set_val(i,j, be->_coef*res);
412 //// }
413 //// }
414 //// }
415 //// return;
416 //// }
417 //// if (be->_lson->is_matrix() && !be->_rson->is_matrix() && be->_rson->_is_transposed) {//matrix * transposed vect
418 //// for (size_t i = 0; i < _dim[0]; i++) {
419 //// for (size_t j = 0; j < _dim[1]; j++) {
420 //// set_val(i,j, be->_coef*be->_lson->get_val(i,j) * be->_rson->get_val(j));
421 //// }
422 //// }
423 //// return;
424 //// }
425 //// if (!be->_lson->is_matrix() && !be->_lson->_is_transposed && be->_rson->is_matrix() ) {//vect * matrix
426 //// for (size_t i = 0; i < _dim[0]; i++) {
427 //// for (size_t j = 0; j < _dim[1]; j++) {
428 //// set_val(i,j, be->_coef*(be->_lson->get_val(i) * be->_rson->get_val(i,j)));
429 //// }
430 //// }
431 //// return;
432 //// }
433 //// // if (be->_lson->is_matrix() && be->_rson->_is_vector) {//matrix*vect
434 //// // for (size_t i = 0; i < _dim[0]; i++) {
435 //// // for (size_t j = 0; j < _dim[1]; j++) {
436 //// // set_val(i,j, be->_coef*(be->_lson->get_val(i,j) * be->_rson->get_val(i)));
437 //// // }
438 //// // }
439 //// // return;
440 //// // }
441 //// }
442 //// if (!_is_hessian) {
443 //// for (size_t i = 0; i < _dim[0]; i++) {
444 //// for (size_t j = 0; j < _dim[1]; j++) {
445 //// auto res = 0.;
446 //// for (size_t col = 0; col<be->_lson->_dim[1]; col++) {
447 //// res += be->_lson->get_val(i,col) * be->_rson->get_val(col,j);
448 //// }
449 //// set_val(i,j,res);
450 //// }
451 //// }
452 //// }
453 //// else {
454 //// for (size_t i = 0; i < _dim[0]; i++) {
455 //// for (size_t j = i; j < _dim[1]; j++) {
456 //// auto res = 0.;
457 //// for (size_t col = 0; col<be->_lson->_dim[1]; col++) {
458 //// res += be->_lson->get_val(i,col) * be->_rson->get_val(col,j);
459 //// }
460 //// set_val(i,j,res);
461 //// }
462 //// }
463 //// }
464 //// }
465 //// }
466 //
467 //// double t_eval(const constant_* c, size_t i){
468 //// if (!c) {
469 //// throw invalid_argument("Cannot evaluate nullptr!");
470 //// }
471 //// switch (c->get_type()) {
472 //// case binary_c: {
473 //// return ((constant<bool>*)(c))->eval();
474 //// break;
475 //// }
476 //// case short_c: {
477 //// return ((constant<short>*)(c))->eval();
478 //// break;
479 //// }
480 //// case integer_c: {
481 //// return ((constant<int>*)(c))->eval();
482 //// break;
483 //// }
484 //// case float_c: {
485 //// return ((constant<float>*)(c))->eval();
486 //// break;
487 //// }
488 //// case double_c: {
489 //// return ((constant<double>*)(c))->eval();
490 //// break;
491 //// }
492 //// case long_c: {
493 //// return (double)((constant<long double>*)(c))->eval();
494 //// break;
495 //// }
496 //// case par_c:{
497 //// auto p_c2 = (param_*)(c);
498 //// switch (p_c2->get_intype()) {
499 //// case binary_:
500 //// return ((param<bool>*)p_c2)->eval(i);
501 //// break;
502 //// case short_:
503 //// return ((param<short>*)p_c2)->eval(i);
504 //// break;
505 //// case integer_:
506 //// return ((param<int>*)p_c2)->eval(i);
507 //// break;
508 //// case float_:
509 //// return ((param<float>*)p_c2)->eval(i);
510 //// break;
511 //// case double_:
512 //// return ((param<double>*)p_c2)->eval(i);
513 //// break;
514 //// case long_:
515 //// return (double)((param<long double>*)p_c2)->eval(i);
516 //// break;
517 //// default:
518 //// break;
519 //// }
520 //// break;
521 //// }
522 //// case var_c:{
523 //// auto p_c2 = (param_*)(c);
524 //// switch (p_c2->get_intype()) {
525 //// case binary_:
526 //// return ((var<bool>*)p_c2)->eval(i);
527 //// break;
528 //// case short_:
529 //// return ((var<short>*)p_c2)->eval(i);
530 //// break;
531 //// case integer_:
532 //// return ((var<int>*)p_c2)->eval(i);
533 //// break;
534 //// case float_:
535 //// return ((var<float>*)p_c2)->eval(i);
536 //// break;
537 //// case double_:
538 //// return ((var<double>*)p_c2)->eval(i);
539 //// break;
540 //// case long_:
541 //// return (double)((var<long double>*)p_c2)->eval(i);
542 //// break;
543 //// default:
544 //// break;
545 //// }
546 //// break;
547 //// }
548 //// case uexp_c: {
549 //// return ((uexpr*)c)->eval(i);
550 //// break;
551 //// }
552 //// case bexp_c: {
553 //// return ((bexpr*)c)->eval(i);
554 //// break;
555 //// }
556 //// case func_c: {
557 //// auto f = (func_*)c;
558 //// if (f->is_constant() && f->_evaluated) {
559 //// if (f->is_number()) {
560 //// return f->_val->at(0);
561 //// }
562 //// return f->get_val(i);
563 //// }
564 //// return ((func_*)c)->eval(i);
565 //// break;
566 //// }
567 //// default:
568 //// break;
569 //// }
570 //// return 0;
571 //// }
572 //
573 //// double t_eval(const constant_* c, size_t i, size_t j){
574 //// if (!c) {
575 //// throw invalid_argument("Cannot evaluate nullptr!");
576 //// }
577 //// // if (!c->is_matrix()) {
578 //// // return eval(c, j);
579 //// // }
580 //// switch (c->get_type()) {
581 //// case binary_c: {
582 //// return ((constant<bool>*)(c))->eval();
583 //// break;
584 //// }
585 //// case short_c: {
586 //// return ((constant<short>*)(c))->eval();
587 //// break;
588 //// }
589 //// case integer_c: {
590 //// return ((constant<int>*)(c))->eval();
591 //// break;
592 //// }
593 //// case float_c: {
594 //// return ((constant<float>*)(c))->eval();
595 //// break;
596 //// }
597 //// case double_c: {
598 //// return ((constant<double>*)(c))->eval();
599 //// break;
600 //// }
601 //// case long_c: {
602 //// return (double)((constant<long double>*)(c))->eval();
603 //// break;
604 //// }
605 //// case par_c:{
606 //// auto p_c2 = (param_*)(c);
607 //// switch (p_c2->get_intype()) {
608 //// case binary_:
609 //// return ((param<bool>*)p_c2)->eval(i,j);
610 //// break;
611 //// case short_:
612 //// return ((param<short>*)p_c2)->eval(i,j);
613 //// break;
614 //// case integer_:
615 //// return ((param<int>*)p_c2)->eval(i,j);
616 //// break;
617 //// case float_:
618 //// return ((param<float>*)p_c2)->eval(i,j);
619 //// break;
620 //// case double_:
621 //// return ((param<double>*)p_c2)->eval(i,j);
622 //// break;
623 //// case long_:
624 //// return (double)((param<long double>*)p_c2)->eval(i,j);
625 //// break;
626 //// default:
627 //// break;
628 //// }
629 //// break;
630 //// }
631 //// case var_c:{
632 //// auto p_c2 = (param_*)(c);
633 //// switch (p_c2->get_intype()) {
634 //// case binary_:
635 //// return ((var<bool>*)p_c2)->eval(i,j);
636 //// break;
637 //// case short_:
638 //// return ((var<short>*)p_c2)->eval(i,j);
639 //// break;
640 //// case integer_:
641 //// return ((var<int>*)p_c2)->eval(i,j);
642 //// break;
643 //// case float_:
644 //// return ((var<float>*)p_c2)->eval(i,j);
645 //// break;
646 //// case double_:
647 //// return ((var<double>*)p_c2)->eval(i,j);
648 //// break;
649 //// case long_:
650 //// return (double)((var<long double>*)p_c2)->eval(i,j);
651 //// break;
652 //// default:
653 //// break;
654 //// }
655 //// break;
656 //// }
657 //// case uexp_c: {
658 //// return ((uexpr*)c)->eval(i,j);
659 //// break;
660 //// }
661 //// case bexp_c: {
662 //// return ((bexpr*)c)->eval(i,j);
663 //// break;
664 //// }
665 //// case func_c: {
666 //// auto f = (func_*)c;
667 //// if (f->is_constant() && f->_evaluated) {
668 //// if (f->is_number()) {
669 //// return f->_val->at(0);
670 //// }
671 //// if (f->is_matrix()) {
672 //// return f->get_val(i, j);
673 //// }
674 //// return f->get_val(j);
675 //// }
676 //// return ((func_*)c)->eval(i,j);
677 //// break;
678 //// }
679 //// default:
680 //// break;
681 //// }
682 //// return 0;
683 //// }
684 //
685 //
686 // func_ get_derivative(constant_* c, const param_ &v){
687 // if (!c) {
688 // throw invalid_argument("Cannot evaluate nullptr!");
689 // }
690 // if (c->is_number() || c->is_param()) {
691 // return func_();
692 // }
693 // if (c->is_var()) {
694 // if ((*(param_*)c)==v) {
695 // func_() += 1;
696 // }
697 // return func_();
698 // }
699 //
700 // if(c->is_uexpr()){
701 // return ((uexpr*)c)->get_derivative(v);
702 // }
703 // if(c->is_bexpr()){
704 // return ((bexpr*)c)->get_derivative(v);
705 // }
706 // if(c->is_function()){
707 // return ((func_*)c)->get_derivative(v);
708 // }
709 // return func_();
710 // }
711 //
712 //
713 // Sign constant_::get_all_sign() const{
714 // switch (_type) {
715 // case binary_c: {
716 // return ((constant<bool>*)this)->get_sign();
717 // break;
718 // }
719 // case short_c: {
720 // return ((constant<short>*)this)->get_sign();
721 // break;
722 // }
723 // case integer_c: {
724 // return ((constant<int>*)this)->get_sign();
725 // break;
726 // }
727 // case float_c: {
728 // return ((constant<float>*)this)->get_sign();
729 // break;
730 // }
731 // case double_c: {
732 // return ((constant<double>*)this)->get_sign();
733 // break;
734 // }
735 // case long_c: {
736 // return ((constant<long double>*)this)->get_sign();
737 // break;
738 // }
739 // case par_c:{
740 // return ((param_*)this)->get_all_sign();
741 // break;
742 // }
743 // case uexp_c: {
744 // return ((uexpr*)this)->get_all_sign(); // TO UPDATE
745 // break;
746 // }
747 // case bexp_c: {
748 // return ((bexpr*)this)->get_all_sign(); // TO UPDATE
749 // break;
750 // }
751 // case var_c:{
752 // return ((param_*)this)->get_all_sign();
753 // break;
754 // }
755 // case func_c: {
756 // return ((func_*)this)->get_all_sign();
757 // break;
758 // }
759 // default:
760 // break;
761 // }
762 // return unknown_;
763 // }
764 //
765 //
766 //
767 // Sign constant_::get_sign(size_t idx) const{
768 // switch (_type) {
769 // case binary_c: {
770 // return ((constant<bool>*)this)->get_sign();
771 // break;
772 // }
773 // case short_c: {
774 // return ((constant<short>*)this)->get_sign();
775 // break;
776 // }
777 // case integer_c: {
778 // return ((constant<int>*)this)->get_sign();
779 // break;
780 // }
781 // case float_c: {
782 // return ((constant<float>*)this)->get_sign();
783 // break;
784 // }
785 // case double_c: {
786 // return ((constant<double>*)this)->get_sign();
787 // break;
788 // }
789 // case long_c: {
790 // return ((constant<long double>*)this)->get_sign();
791 // break;
792 // }
793 // case par_c:{
794 // return ((param_*)this)->get_sign(idx);
795 // break;
796 // }
797 // case uexp_c: {
798 // return ((uexpr*)this)->get_sign(idx);
799 // break;
800 // }
801 // case bexp_c: {
802 // return ((bexpr*)this)->get_sign(idx);
803 // break;
804 // }
805 // case var_c:{
806 // return ((param_*)this)->get_sign(idx);
807 // break;
808 // }
809 // case func_c: {
810 // return ((func_*)this)->get_sign(idx);
811 // break;
812 // }
813 // default:
814 // break;
815 // }
816 // return unknown_;
817 // }
818 //
819 // Sign param_::get_sign(size_t idx) const{
820 // switch (_intype) {
821 // case binary_:
822 // if (is_param()) {
823 // return ((param<bool>*)this)->get_sign(idx);
824 // }
825 // return ((var<bool>*)this)->get_sign(idx);
826 // break;
827 // case short_:
828 // if (is_param()) {
829 // return ((param<short>*)this)->get_sign(idx);
830 // }
831 // return ((var<short>*)this)->get_sign(idx);
832 // break;
833 // case integer_:
834 // if (is_param()) {
835 // return ((param<int>*)this)->get_sign(idx);
836 // }
837 // return ((var<int>*)this)->get_sign(idx);
838 // break;
839 // case float_:
840 // if (is_param()) {
841 // return ((param<float>*)this)->get_sign(idx);
842 // }
843 // return ((var<float>*)this)->get_sign(idx);
844 // break;
845 // case double_:
846 // if (is_param()) {
847 // return ((param<double>*)this)->get_sign(idx);
848 // }
849 // return ((var<double>*)this)->get_sign(idx);
850 // break;
851 // case long_:
852 // if (is_param()) {
853 // return ((param<long double>*)this)->get_sign(idx);
854 // }
855 // return ((var<long double>*)this)->get_sign(idx);
856 // break;
857 // default:
858 // throw invalid_argument("Unsupported type");
859 // break;
860 // }
861 // }
862 //
863 // Sign param_::get_all_sign() const{
864 // switch (_intype) {
865 // case binary_:
866 // return pos_;
867 // break;
868 // case short_:
869 // if (is_param()) {
870 // return ((param<short>*)this)->get_all_sign<short>();
871 // }
872 // return ((var<short>*)this)->get_all_sign<short>();
873 // break;
874 // case integer_:
875 // if (is_param()) {
876 // return ((param<int>*)this)->get_all_sign<int>();
877 // }
878 // return ((var<int>*)this)->get_all_sign<int>();
879 // break;
880 // case float_:
881 // if (is_param()) {
882 // return ((param<float>*)this)->get_all_sign<float>();
883 // }
884 // return ((var<float>*)this)->get_all_sign<float>();
885 // break;
886 // case double_:
887 // if (is_param()) {
888 // return ((param<double>*)this)->get_all_sign<double>();
889 // }
890 // return ((var<double>*)this)->get_all_sign<double>();
891 // break;
892 // case long_:
893 // if (is_param()) {
894 // return ((param<long double>*)this)->get_all_sign<long double>();
895 // }
896 // return ((var<long double>*)this)->get_all_sign<long double>();
897 // break;
898 // case complex_:
899 // if (is_param()) {
900 // return ((param<Cpx>*)this)->get_all_sign_cpx();
901 // }
902 // return ((var<Cpx>*)this)->get_all_sign_cpx();
903 // break;
904 // default:
905 // return unknown_;
906 // // throw invalid_argument("Unsupported type");
907 // break;
908 // }
909 //
910 // }
911 //
912 //
913 //
914 //
915 //
916 //
917 //
918 // bool func_::operator==(const func_& f) const{
919 // if (_ftype!=f._ftype || _all_sign!=f._all_sign || _vars->size()!=f._vars->size() || _pterms->size() != f._pterms->size()|| _qterms->size() != f._qterms->size() || _lterms->size() != f._lterms->size()) {
920 // return false;
921 // }
922 // if (this->to_str()!=f.to_str()) {
923 // return false;
924 // }
925 // return true;
926 // }
927 //
928 // func_& func_::operator=(const func_& f){
929 // // _evaluated = f._evaluated;
930 // _to_str = f._to_str;
931 //
932 // delete _all_range;
933 // if (_vars) {
934 // _vars->clear();
935 // }
936 // if (_params){
937 // _params->clear();
938 // }
939 // _dfdx->clear();
940 // delete _range;
941 // delete _convexity;
942 // delete _sign;
943 // delete _lterms;
944 // delete _qterms;
945 // delete _pterms;
946 // delete _DAG;
947 // delete _queue;
948 // delete _cst;
949 //
950 // set_type(func_c);
951 // _lterms = new map<string, lterm>();
952 // _qterms = new map<string, qterm>();
953 // _pterms = new map<string, pterm>();
954 // _DAG = new map<string, expr*>();
955 // _queue = new deque<shared_ptr<expr>>();
956 // _ftype = f._ftype;
957 // _return_type = f._return_type;
958 // _all_convexity = f._all_convexity;
959 // _all_sign = f._all_sign;
960 // _is_transposed = f._is_transposed;
961 // _is_vector = f._is_vector;
962 // // _is_matrix = f._is_matrix;
963 // _is_constraint = f._is_constraint;
964 // _is_hessian = f._is_hessian;
965 // _dim[0] = f._dim[0];
966 // _dim[1] = f._dim[1];
967 // _embedded = f._embedded;
968 // _return_type = f._return_type;
969 // _cst = copy(*f._cst);
970 // if (_cst->is_function()) {
971 // embed(*(func_*)_cst);
972 // }
973 //
974 // _all_range = new pair<double,double>(f._all_range->first, f._all_range->second);
975 // _sign = nullptr;
976 // _convexity = nullptr;
977 // _range = nullptr;
978 //
979 // if(f._sign){
980 // _sign = new vector<Sign>();
981 // *_sign = *f._sign;
982 // }
983 // if (f._convexity) {
984 // _convexity = new vector<Convexity>();
985 // *_convexity = *f._convexity;
986 // }
987 // if (f._range) {
988 // _range= new vector<pair<double,double>>();
989 // *_range = *f._range;// Make sure this creates new pointers inside each pair, otherwise use below.
990 // // for (auto &elem: *f._range) {
991 // // _range->push_back(make_pair<>(copy(elem.first), copy(elem.second)));
992 // // }
993 // }
994 //
995 // for (auto &pair:*f._lterms) {
996 // insert(pair.second);
997 // }
998 // for (auto &pair:*f._qterms) {
999 // insert(pair.second);
1000 // }
1001 // for (auto &pair:*f._pterms) {
1002 // insert(pair.second);
1003 // }
1004 // if (f._val) {
1005 // _val = make_shared<vector<double>>(*f._val);
1006 // }
1007 // if (f._expr) {
1008 // if (f._expr->is_uexpr()) {
1009 // _expr = make_shared<uexpr>(*(uexpr*)(f._expr.get()));
1010 // }
1011 // else {
1012 // _expr = make_shared<bexpr>(*(bexpr*)(f._expr.get()));
1013 // }
1014 // embed(_expr);
1015 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1016 // _queue->push_back(_expr);
1017 // }
1018 // for (auto &df:*f._dfdx) {
1019 // (*_dfdx)[df.first] = make_shared<func_>(*df.second);
1020 // }
1021 // _nnz_j = f._nnz_j;
1022 // _nnz_h = f._nnz_h;
1023 // _hess_link = f._hess_link;
1024 // // _nb_instances = f._nb_instances;
1025 // _nb_vars = f._nb_vars;
1026 // _indices = f._indices;
1027 // if (is_constant()) {
1028 // _evaluated = f._evaluated;
1029 // }
1030 //
1031 // return *this;
1032 // }
1033 //
1034 //
1035 // func_& func_::operator=(func_&& f){
1036 // // _evaluated = f._evaluated;
1037 // _to_str = f._to_str;
1038 //
1039 // delete _all_range;
1040 // if (_vars) {
1041 // _vars->clear();
1042 // delete _vars;
1043 // }
1044 // if (_params){
1045 // _params->clear();
1046 // delete _params;
1047 // }
1048 //
1049 // delete _range;
1050 // delete _convexity;
1051 // delete _sign;
1052 // delete _lterms;
1053 // delete _qterms;
1054 // delete _pterms;
1055 // delete _DAG;
1056 // delete _queue;
1057 // delete _cst;
1058 // set_type(func_c);
1059 // _ftype = f._ftype;
1060 // _return_type = f._return_type;
1061 // _all_convexity = f._all_convexity;
1062 // _all_sign = f._all_sign;
1063 // _all_range = f._all_range;
1064 // f._all_range = nullptr;
1065 // _lterms = f._lterms;
1066 // f._lterms = nullptr;
1067 // _qterms = f._qterms;
1068 // f._qterms = nullptr;
1069 // _pterms = f._pterms;
1070 // f._pterms = nullptr;
1071 // _expr = move(f._expr);
1072 // _DAG = f._DAG;
1073 // f._DAG = nullptr;
1074 // _queue = f._queue;
1075 // f._queue = nullptr;
1076 // _vars = f._vars;
1077 // f._vars = nullptr;
1078 // _params = f._params;
1079 // f._params = nullptr;
1080 // _cst = f._cst;
1081 // f._cst = nullptr;
1082 // _range = f._range;
1083 // f._range = nullptr;
1084 // _convexity = f._convexity;
1085 // f._convexity = nullptr;
1086 // _sign = f._sign;
1087 // f._sign = nullptr;
1088 // _is_transposed = f._is_transposed;
1089 // _is_vector = f._is_vector;
1090 // // _is_matrix = f._is_matrix;
1091 // _is_constraint = f._is_constraint;
1092 // _is_hessian = f._is_hessian;
1093 // _dim[0] = f._dim[0];
1094 // _dim[1] = f._dim[1];
1095 // _embedded = f._embedded;
1096 // _nnz_j = f._nnz_j;
1097 // _nnz_h = f._nnz_h;
1098 // _hess_link = f._hess_link;
1099 // // _nb_instances = f._nb_instances;
1100 // _nb_vars = f._nb_vars;
1101 // _dfdx = move(f._dfdx);
1102 // _val = move(f._val);
1103 // _indices = move(f._indices);
1104 // if (is_constant()) {
1105 // _evaluated = f._evaluated;
1106 // }
1107 // return *this;
1108 // }
1109 //
1110 //
1111 //
1112 // func_::~func_(){
1113 // delete _all_range;
1114 // // for (auto &pt: *_lterms) {
1115 // // if (pt.second._coef->is_function()) {
1116 // // delete pt.second._coef;
1117 // // }
1118 // // }
1119 // // for (auto &pt: *_qterms) {
1120 // // if (pt.second._coef->is_function()) {
1121 // // delete pt.second._coef;
1122 // // }
1123 // // }
1124 // // for (auto &pt: *_pterms) {
1125 // // if (pt.second._coef->is_function()) {
1126 // // delete pt.second._coef;
1127 // // }
1128 // // }
1129 // if (_vars) {
1130 // _vars->clear();
1131 // delete _vars;
1132 // }
1133 // if (_params) {
1134 // _params->clear();
1135 // delete _params;
1136 // }
1137 //
1138 // delete _range;
1139 // delete _convexity;
1140 // delete _sign;
1141 // delete _lterms;
1142 // delete _qterms;
1143 // delete _pterms;
1144 // delete _DAG;
1145 // delete _queue;
1146 // delete _cst;
1147 // };
1148 //
1149 // bool all_zeros(const string& s){
1150 // auto it = s.begin();
1151 // while (it != s.end()) {
1152 // if ((*it)!='-' && (*it)!='0' && (*it)!='.') {
1153 // return false;
1154 // }
1155 // it++;
1156 // }
1157 // return true;
1158 // }
1159 //
1160 // bool constant_::is_zero() const{ /**< Returns true if constant equals 0 */
1161 // // auto a = poly_to_str(this);
1162 // if (is_number() && all_zeros(poly_to_str(this))){
1163 // return true;
1164 // }
1165 // if (is_param() || is_var()) {
1166 // auto p_c = (param_*)this;
1167 // return p_c->get_all_sign()==zero_;
1168 // }
1169 // // if (is_uexpr() || is_bexpr()) {
1170 // // auto e_p = (expr*)this;
1171 // // return e_p->get_all_sign()==zero_;
1172 // // }
1173 // if (is_function()) {
1174 // return ((func_*)(this))->is_zero();
1175 // }
1176 // return false;
1177 // }
1178 //
1179 // bool constant_::is_unit() const{ /**< Returns true if constant equals 1 */
1180 // if(is_number() && t_eval(this)==1 && !_is_vector && !_is_transposed){
1181 // return true;
1182 // }
1183 // // if (is_param()) {
1184 // // auto p_c = (param_*)this;
1185 // // switch (p_c->get_intype()) {
1186 // // case binary_:
1187 // // return ((param<bool>*)p_c)->is_unit();
1188 // // break;
1189 // // case short_:
1190 // // return ((param<short>*)p_c)->is_unit();
1191 // // break;
1192 // // case integer_:
1193 // // return ((param<int>*)p_c)->is_unit();
1194 // // break;
1195 // // case float_:
1196 // // return ((param<float>*)p_c)->is_unit();
1197 // // break;
1198 // // case double_:
1199 // // return ((param<double>*)p_c)->is_unit();
1200 // // break;
1201 // // case long_:
1202 // // return ((param<long double>*)p_c)->is_unit();
1203 // // break;
1204 // // default:
1205 // // break;
1206 // // }
1207 // // }
1208 // if (is_function()) {
1209 // return ((func_*)(this))->is_unit();
1210 // }
1211 // return false;
1212 // }
1213 //
1214 //
1215 // bool constant_::is_positive() const{
1216 // if (get_all_sign()==pos_) {
1217 // return true;
1218 // }
1219 // return false;
1220 // }
1221 //
1222 //
1223 // bool constant_::is_non_positive() const{
1224 // if (get_all_sign()==non_pos_) {
1225 // return true;
1226 // }
1227 // return false;
1228 // }
1229 //
1230 // bool constant_::is_non_negative() const{
1231 // if (get_all_sign()==non_neg_) {
1232 // return true;
1233 // }
1234 // return false;
1235 // }
1236 //
1237 // bool constant_::is_negative() const{
1238 // if (get_all_sign()==neg_) {
1239 // return true;
1240 // }
1241 // return false;
1242 // }
1243 //
1244 // void func_::update_sign(){
1245 // //TODO
1246 // }
1247 //
1248 //
1249 // func_& func_::operator+=(const constant_& c){
1250 // if (c.is_zero()) {
1251 // return *this;
1252 // }
1253 // _evaluated = false;
1254 // if (c.is_number() || (!is_constant() && c.is_param())) {
1255 // _cst = add(_cst, c);
1256 // if (_cst->is_function()) {
1257 // embed(*(func_*)_cst);
1258 // }
1259 // update_sign(*_cst);
1260 // return *this;
1261 // }
1262 // if (c.is_param() || c.is_var()) {
1263 // this->insert(true, constant<double>(1), *(param_*)&c);
1264 // }
1265 // if (c.is_function()) {
1266 // func_* f = (func_*)&c;
1267 // if (is_constant() && !f->is_constant()) {//check _expr
1268 // return *this = (*f + *this);
1269 // }
1270 // if (!is_constant() && f->is_constant()) {//check _expr
1271 // _cst = add(_cst, c);
1272 // if (_cst->is_function()) {
1273 // embed(*(func_*)_cst);
1274 // }
1275 // return *this;
1276 // }
1277 // if (!f->_cst->is_zero()) {
1278 // _cst = add(_cst, *f->_cst);
1279 // if (_cst->is_function()) {
1280 // embed(*(func_*)_cst);
1281 // }
1282 // }
1283 // for (auto &pair:*f->_lterms) {
1284 // this->insert(pair.second);
1285 // }
1286 // for (auto &pair:*f->_qterms) {
1287 // this->insert(pair.second);
1288 // }
1289 // for (auto &pair:*f->_pterms) {
1290 // this->insert(pair.second);
1291 // }
1292 //
1293 // if (_expr && f->_expr) {
1294 // _expr = make_shared<bexpr>(bexpr(plus_, make_shared<func_>((*_expr)), make_shared<func_>((*f->_expr))));
1295 // embed(_expr);
1296 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1297 // _queue->push_back(_expr);
1298 // }
1299 // else if (!_expr && f->_expr) {
1300 // if (f->_expr->is_uexpr()) {
1301 // _expr = make_shared<uexpr>(*(uexpr*)(f->_expr.get()));
1302 // }
1303 // else {
1304 // _expr = make_shared<bexpr>(*(bexpr*)(f->_expr.get()));
1305 // }
1306 // embed(_expr);
1307 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1308 // _queue->push_back(_expr);
1309 // if (!_vars->empty()) {
1310 // _ftype = nlin_;
1311 // }
1312 // }
1313 // update_sign(*f);
1314 // if (_all_convexity==linear_ && f->_all_convexity==convex_) {
1315 // _all_convexity = convex_;
1316 // }
1317 // else if (_all_convexity==linear_ && f->_all_convexity==concave_) {
1318 // _all_convexity = concave_;
1319 // }
1320 // else if (_all_convexity==convex_ && f->_all_convexity==convex_) {
1321 // _all_convexity = convex_;
1322 // }
1323 // else if (_all_convexity==concave_ && f->_all_convexity==concave_) {
1324 // _all_convexity = concave_;
1325 // }
1326 // else if (_all_convexity==convex_ && (f->_all_convexity==convex_ || f->_all_convexity==undet_)) {
1327 // _all_convexity = undet_;
1328 // }
1329 // else if (_all_convexity==concave_ && (f->_all_convexity==concave_ || f->_all_convexity==undet_)) {
1330 // _all_convexity = undet_;
1331 // }
1332 // update_convexity();
1333 // return *this;
1334 // }
1335 // return *this;
1336 // }
1337 //
1338 // func_& func_::operator-=(const constant_& c){
1339 // if (c.is_zero()) {
1340 // return *this;
1341 // }
1342 // _evaluated = false;
1343 // if (c.is_number() || (!is_constant() && c.is_param())) {
1344 // _cst = substract(_cst, c);
1345 // if (_cst->is_function()) {
1346 // embed(*(func_*)_cst);
1347 // }
1348 // update_sign(*_cst);
1349 // return *this;
1350 // }
1351 // if (c.is_param() || c.is_var()) {
1352 // this->insert(false, constant<double>(1), *(param_*)&c);
1353 // }
1354 // if (c.is_function()) {
1355 // func_* f = (func_*)&c;
1356 // if (!is_constant() && f->is_constant()) {
1357 // _cst = substract(_cst, c);
1358 // if (_cst->is_function()) {
1359 // embed(*(func_*)_cst);
1360 // }
1361 // return *this;
1362 // }
1363 // if (is_constant() && !f->is_constant()) {//check _expr
1364 // return *this = (-1*(*f) + *this);
1365 // }
1366 // _cst = substract(_cst, *f->_cst);
1367 // if (_cst->is_function()) {
1368 // embed(*(func_*)_cst);
1369 // }
1370 // for (auto &pair:*f->_lterms) {
1371 // this->insert(!pair.second._sign, *pair.second._coef, *pair.second._p);
1372 // }
1373 // for (auto &pair:*f->_qterms) {
1374 // this->insert(!pair.second._sign, *pair.second._coef, *pair.second._p->first, *pair.second._p->second);
1375 // }
1376 // for (auto &pair:*f->_pterms) {
1377 // this->insert(!pair.second._sign, *pair.second._coef, *pair.second._l);
1378 // }
1379 // if (_expr && f->_expr) {
1380 // _expr = make_shared<bexpr>(bexpr(minus_, make_shared<func_>((*_expr)), make_shared<func_>((*f->_expr))));
1381 // embed(_expr);
1382 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1383 // _queue->push_back(_expr);
1384 // }
1385 // if (!_expr && f->_expr) {
1386 // if (f->_expr->is_uexpr()) {
1387 // _expr = make_shared<uexpr>(*(uexpr*)(f->_expr.get()));
1388 // }
1389 // else {
1390 // _expr = make_shared<bexpr>(*(bexpr*)(f->_expr.get()));
1391 // }
1392 // _expr->_coef *= -1;
1393 // embed(_expr);
1394 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1395 // _queue->push_back(_expr);
1396 // if (!_vars->empty()) {
1397 // _ftype = nlin_;
1398 // }
1399 // }
1400 // update_sign(*f);
1401 // if (_all_convexity==linear_ && f->_all_convexity==concave_) {
1402 // _all_convexity = convex_;
1403 // }
1404 // else if (_all_convexity==linear_ && f->_all_convexity==convex_) {
1405 // _all_convexity = concave_;
1406 // }
1407 // else if (_all_convexity==convex_ && f->_all_convexity==concave_) {
1408 // _all_convexity = convex_;
1409 // }
1410 // else if (_all_convexity==concave_ && f->_all_convexity==convex_) {
1411 // _all_convexity = concave_;
1412 // }
1413 // else if (_all_convexity==convex_ && (f->_all_convexity==concave_ || f->_all_convexity==undet_)) {
1414 // _all_convexity = undet_;
1415 // }
1416 // else if (_all_convexity==concave_ && (f->_all_convexity==convex_ || f->_all_convexity==undet_)) {
1417 // _all_convexity = undet_;
1418 // }
1419 // update_convexity();
1420 // return *this;
1421 // }
1422 // return *this;
1423 // }
1424 //
1425 //
1426 // func_& func_::operator*=(const constant_& c){
1427 // if (is_zero()) {
1428 // return *this;
1429 // }
1430 // if (c.is_zero()) {
1431 // reset();
1432 // return *this;
1433 // }
1434 // if (is_unit()) {
1435 // *this = func_(c);
1436 // return *this;
1437 // }
1438 // if (c.is_unit()) {
1439 // return *this;
1440 // }
1441 //
1442 // /* Case where c is a number */
1443 // if (c.is_number()){
1444 // if (!_cst->is_zero()) {
1445 // _cst = multiply(_cst, c);
1446 // if (_cst->is_function()) {
1447 // embed(*(func_*)_cst);
1448 // }
1449 // }
1450 // for (auto &pair:*_lterms) {
1451 // pair.second._coef = multiply(pair.second._coef, c);
1452 // }
1453 // for (auto &pair:*_qterms) {
1454 // pair.second._coef = multiply(pair.second._coef, c);
1455 // }
1456 // for (auto &pair:*_pterms) {
1457 // pair.second._coef = multiply(pair.second._coef, c);
1458 // }
1459 // if (c.is_negative()) {
1460 // reverse_sign();
1461 // }
1462 // auto val = t_eval(&c);
1463 // if (_expr) {
1464 // _expr->_coef *= val;
1465 // }
1466 // if (_evaluated) {
1467 // for (unsigned i = 0; i<_val->size(); i++) {
1468 // _val->at(i) *= val;
1469 // }
1470 // }
1471 // // this->update_dot_dim(c);
1472 // return *this;
1473 // }
1474 //// if (c.is_matrix() || is_matrix()) {
1475 //// *this = func_(bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c)));
1476 //// _evaluated = false;
1477 //// return *this;
1478 //// }
1479 // /* Case where the current function is not constant and the other operand is */
1480 // if(!is_constant() && (c.is_param() || (c.is_function() && ((func_*)&c)->is_constant()))) {
1481 // bool transp = false;
1482 // auto fc = func_(c);
1483 // if(is_linear() && _is_transposed){// Situation where var^T * c is transformed into (c^T*var)^T
1484 // fc.transpose();
1485 // this->transpose();
1486 // transp = true;
1487 // }
1488 // if (!_cst->is_zero()) {
1489 // _cst = multiply(_cst, fc);
1490 // if (_cst->is_function()) {
1491 // embed(*(func_*)_cst);
1492 // }
1493 // }
1494 // for (auto &pair:*_lterms) {
1495 // pair.second._coef = multiply(pair.second._coef, fc);
1496 // if (pair.second._coef->is_function()) {
1497 // embed(*(func_*)pair.second._coef);
1498 // }
1499 // // update_nb_instances(pair.second);
1500 //
1501 // }
1502 // for (auto &pair:*_qterms) {
1503 // pair.second._coef = multiply(pair.second._coef, fc);
1504 // if (pair.second._coef->is_function()) {
1505 // embed(*(func_*)pair.second._coef);
1506 // }
1507 // // update_nb_instances(pair.second);
1508 // }
1509 // for (auto &pair:*_pterms) {
1510 // pair.second._coef = multiply(pair.second._coef, fc);
1511 // if (pair.second._coef->is_function()) {
1512 // embed(*(func_*)pair.second._coef);
1513 // }
1514 // // update_nb_instances(pair.second); // TODO
1515 // }
1516 // if (_expr) {
1517 // if (c.is_function() && ((func_*)&c)->is_number()) {
1518 // _expr->_coef *= t_eval(&c);
1519 // }
1520 // else {
1521 // _expr = make_shared<bexpr>(bexpr(product_, make_shared<func_>((*_expr)), make_shared<func_>((c))));
1522 // // _dim[0] = _expr->_dim[0];
1523 // // _dim[1] = _expr->_dim[1];
1524 // embed(_expr);
1525 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1526 // _queue->push_back(_expr);
1527 // }
1528 // }
1529 // if (c.is_negative()) {
1530 // reverse_sign();
1531 // }
1532 // if (c.get_all_sign()==unknown_) {
1533 // _all_sign = unknown_;
1534 // if (!_qterms->empty()) {
1535 // _all_convexity = undet_;
1536 // }
1537 // }
1538 // _evaluated = false;
1539 // if(transp){
1540 // this->transpose();
1541 // }
1542 // else {
1543 // this->update_dot(c);
1544 // }
1545 // return *this;
1546 // }
1547 // /* Case where the current function is constant and the other operand is not. */
1548 // if (is_constant() && (c.is_var() || (c.is_function() && !((func_*)&c)->is_constant()))) {
1549 // func_ f(c);
1550 // if (!f._cst->is_zero()) {
1551 // if (f._cst->is_function()) {
1552 // auto fc = (func_*)f._cst;
1553 // *fc = (*this)* (*fc);
1554 // f.embed(*fc);
1555 // }
1556 // else {
1557 // f._cst = multiply(f._cst, *this);
1558 // }
1559 // }
1560 // for (auto &pair:*f._lterms) {
1561 // if (pair.second._coef->is_function()) {
1562 // auto fc = (func_*)pair.second._coef;
1563 // *fc = (*this)* (*fc);
1564 // f.embed(*fc);
1565 // }
1566 // else {
1567 // pair.second._coef = multiply(pair.second._coef, *this);
1568 // }
1569 // if (pair.second._coef->_is_transposed) {
1570 // pair.second._p->_is_vector = true;
1571 // if (!pair.second._coef->is_number() && pair.second._coef->_dim[1]!=pair.second._p->_dim[0]) {
1572 // DebugOn("vector dot product with mismatching dimensions, check your param/var dimensions");
1573 // }
1574 // }
1575 // // f.update_nb_instances(pair.second);
1576 //
1577 // }
1578 // for (auto &pair:*f._qterms) {
1579 // if (pair.second._coef->is_function()) {
1580 // auto fc = (func_*)pair.second._coef;
1581 // *fc = (*this)* (*fc);
1582 // f.embed(*fc);
1583 // }
1584 // else {
1585 // pair.second._coef = multiply(pair.second._coef, *this);
1586 // }
1587 // if (pair.second._coef->_is_transposed) {
1588 // pair.second._p->first->_is_vector = true;
1589 // pair.second._p->second->_is_vector = true;
1590 // if (!pair.second._coef->is_number() && pair.second._coef->_dim[1]!=pair.second._p->first->_dim[0]) {
1591 // DebugOn("vector dot product with mismatching dimensions, check your param/var dimensions");
1592 // }
1593 //
1594 // }
1595 // // update_nb_instances(pair.second);
1596 // }
1597 // for (auto &pair:*f._pterms) {
1598 // if (pair.second._coef->is_function()) {
1599 // auto fc = (func_*)pair.second._coef;
1600 // *fc = (*this)* (*fc);
1601 // f.embed(*fc);
1602 // }
1603 // else {
1604 // pair.second._coef = multiply(pair.second._coef, *this);
1605 // }
1606 // // update_nb_instances(pair.second); // TODO
1607 // }
1608 // if (f._expr) {
1609 // if (this->is_number()) {
1610 // f._expr->_coef *= t_eval(this);
1611 // }
1612 // else {
1613 // f._expr = make_shared<bexpr>(bexpr(product_, make_shared<func_>(*this), make_shared<func_>((*f._expr))));
1614 // f.embed(f._expr);
1615 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1616 // f._queue->push_back(f._expr);
1617 // }
1618 // }
1619 // if (this->is_negative()) {
1620 // f.reverse_sign();
1621 // }
1622 // if (this->get_all_sign()==unknown_) {
1623 // f._all_sign = unknown_;
1624 // if (!f._qterms->empty()) {
1625 // f._all_convexity = undet_;
1626 // }
1627 // }
1628 // f._evaluated = false;
1629 // if(update_dot(c)){
1630 // f._dim[0] = _dim[0];
1631 // f._dim[1] = _dim[1];
1632 // }
1633 // return *this = move(f);
1634 // }
1635 // if (c.is_param() || c.is_var()) {
1636 //// if (c.is_matrix() || is_matrix()) {
1637 //// *this = func_(bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c)));
1638 //// }
1639 //// else {
1640 // if(_is_transposed && !c._is_vector){
1641 // auto new_c = copy(c);
1642 // auto new_p = (param_*)new_c;
1643 // new_p->_is_vector = true;
1644 // new_p->_name = "["+new_p->_name+"]";
1645 // func_ f(*new_p);
1646 // *this *= f;
1647 // delete new_c;
1648 // }
1649 // else {
1650 // func_ f(c);
1651 // *this *= f;
1652 // }
1653 //// }
1654 // _evaluated = false;
1655 // return *this;
1656 // }
1657 // if (_expr || (c.is_function() && ((func_*)&c)->_expr)) {
1658 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1659 // *this = func_(be);
1660 // _evaluated = false;
1661 // return *this;
1662 // }
1663 // /* Case where the multiplication invlolves multiplying variables/parameters together, i.e., they are both parametric or both include variables */
1664 // if (c.is_function()) {
1665 // func_* f = (func_*)&c;
1666 // constant_* coef;
1667 // vector<bool>* is_sum = nullptr;
1668 // func_ res;
1669 // for (auto& t1: *_pterms) {
1670 // if (t1.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial function), we cannot factor the coefficients. Just create a binary expression and return it.
1671 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1672 // *this = func_(be);
1673 // _evaluated = false;
1674 // return *this;
1675 // }
1676 // for (auto& t2: *f->_pterms) {
1677 // is_sum = nullptr;
1678 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial function), see comment above.
1679 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1680 // *this = func_(be);
1681 // _evaluated = false;
1682 // return *this;
1683 // }
1684 // auto newl(*t1.second._l);
1685 // for (auto& it: *t2.second._l) {
1686 // newl.push_back(make_pair<>(it.first, it.second));
1687 // }
1688 // coef = copy(*t1.second._coef);
1689 // coef = multiply(coef, *t2.second._coef);
1690 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1691 // delete coef;
1692 // }
1693 // for (auto& t2: *f->_qterms) {
1694 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial function), see comment above.
1695 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1696 // *this = func_(be);
1697 // _evaluated = false;
1698 // return *this;
1699 // }
1700 // auto newl(*t1.second._l);
1701 // newl.push_back(make_pair<>((t2.second._p->first), 1));
1702 // newl.push_back(make_pair<>((t2.second._p->second), 1));
1703 // coef = copy(*t1.second._coef);
1704 // coef = multiply(coef, *t2.second._coef);
1705 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1706 // delete coef;
1707 // }
1708 // for (auto& t2: *f->_lterms) {
1709 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial function)
1710 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1711 // *this = func_(be);
1712 // _evaluated = false;
1713 // return *this;
1714 // }
1715 // auto newl(*t1.second._l);
1716 // newl.push_back(make_pair<>((t2.second._p), 1));
1717 // coef = copy(*t1.second._coef);
1718 // coef = multiply(coef, *t2.second._coef);
1719 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1720 // delete coef;
1721 // }
1722 // if (!f->_cst->is_zero()) {
1723 // auto newl(*t1.second._l);
1724 // coef = copy(*f->_cst);
1725 // coef = multiply(coef, *t1.second._coef);
1726 // res.insert(t1.second._sign, *coef, newl);
1727 // delete coef;
1728 // }
1729 // }
1730 //
1731 // for (auto& t1: *_qterms) {
1732 // if (t1.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(Quadratic term)
1733 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1734 // *this = func_(be);
1735 // _evaluated = false;
1736 // return *this;
1737 // }
1738 // for (auto& t2: *f->_pterms) {
1739 // is_sum = nullptr;
1740 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1741 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1742 // *this = func_(be);
1743 // _evaluated = false;
1744 // return *this;
1745 // }
1746 // auto newl(*t2.second._l);
1747 // newl.push_front(make_pair<>(t1.second._p->first, 1));
1748 // newl.push_front(make_pair<>(t1.second._p->second, 1));
1749 // coef = copy(*t1.second._coef);
1750 // coef = multiply(coef, *t2.second._coef);
1751 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1752 // delete coef;
1753 // }
1754 // for (auto& t2: *f->_qterms) {
1755 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1756 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1757 // *this = func_(be);
1758 // _evaluated = false;
1759 // return *this;
1760 // }
1761 // coef = copy(*t1.second._coef);
1762 // coef = multiply(coef, *t2.second._coef);
1763 // list<pair<param_*, int>> newl;
1764 // newl.push_back(make_pair<>(t1.second._p->first, 1));
1765 // newl.push_back(make_pair<>(t1.second._p->second, 1));
1766 // newl.push_back(make_pair<>(t2.second._p->first, 1));
1767 // newl.push_back(make_pair<>(t2.second._p->second, 1));
1768 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1769 // delete coef;
1770 // }
1771 // for (auto& t2: *f->_lterms) {
1772 // is_sum = nullptr;
1773 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1774 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1775 // *this = func_(be);
1776 // _evaluated = false;
1777 // return *this;
1778 // }
1779 // coef = copy(*t1.second._coef);
1780 // coef = multiply(coef, *t2.second._coef);
1781 // list<pair<param_*, int>> newl;
1782 // newl.push_back(make_pair<>(t1.second._p->first, 1));
1783 // newl.push_back(make_pair<>(t1.second._p->second, 1));
1784 // newl.push_back(make_pair<>(t2.second._p, 1));
1785 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1786 // delete coef;
1787 // }
1788 // if (!f->_cst->is_zero()) {
1789 // coef = copy(*t1.second._coef);
1790 // coef = multiply(coef, *f->_cst);
1791 // res.insert(t1.second._sign, *coef, *t1.second._p->first, *t1.second._p->second);
1792 // delete coef;
1793 // }
1794 //
1795 // }
1796 // for (auto& t1: *_lterms) {
1797 //// if (t1.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(Quadratic term)
1798 //// auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1799 //// *this = func_(be);
1800 //// _evaluated = false;
1801 //// return *this;
1802 //// }
1803 // for (auto& t2: *f->_pterms) {
1804 // is_sum = nullptr;
1805 // if (t1.second._coef->_is_transposed || t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1806 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1807 // *this = func_(be);
1808 // _evaluated = false;
1809 // return *this;
1810 // }
1811 // auto newl(*t2.second._l);
1812 // newl.push_front(make_pair<>((t1.second._p), 1));
1813 // coef = copy(*t1.second._coef);
1814 // coef = multiply(coef, *t2.second._coef);
1815 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1816 // delete coef;
1817 // }
1818 // for (auto& t2: *f->_qterms) {
1819 // if (t1.second._coef->_is_transposed || t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1820 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1821 // *this = func_(be);
1822 // _evaluated = false;
1823 // return *this;
1824 // }
1825 // coef = copy(*t1.second._coef);
1826 // coef = multiply(coef, *t2.second._coef);
1827 // list<pair<param_*, int>> newl;
1828 // newl.push_back(make_pair<>(t1.second._p, 1));
1829 // newl.push_back(make_pair<>(t2.second._p->first, 1));
1830 // newl.push_back(make_pair<>(t2.second._p->second, 1));
1831 // res.insert(!(t1.second._sign^t2.second._sign), *coef, newl);
1832 // delete coef;
1833 // }
1834 // for (auto& t2: *f->_lterms) {
1835 //// if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1836 //// auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1837 //// *this = func_(be);
1838 //// _evaluated = false;
1839 //// return *this;
1840 //// }
1841 // coef = copy(*t1.second._coef);
1842 // coef = multiply(coef, *t2.second._coef);
1843 // res.insert(!(t1.second._sign^t2.second._sign), *coef, *t1.second._p, *t2.second._p, _is_transposed);
1844 // delete coef;
1845 // }
1846 // if (!f->_cst->is_zero()) {
1847 // coef = copy(*t1.second._coef);
1848 // coef = multiply(coef, *f->_cst);
1849 // res.insert(t1.second._sign, *coef, *t1.second._p);
1850 // delete coef;
1851 // }
1852 // }
1853 // if (!_cst->is_zero()) {
1854 // for (auto& t2: *f->_pterms) {
1855 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1856 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1857 // *this = func_(be);
1858 // _evaluated = false;
1859 // return *this;
1860 // }
1861 // coef = copy(*_cst);
1862 // coef = multiply(coef, *t2.second._coef);
1863 // res.insert(t2.second._sign, *coef, *t2.second._l);
1864 // delete coef;
1865 // }
1866 // for (auto& t2: *f->_qterms) {
1867 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1868 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1869 // *this = func_(be);
1870 // _evaluated = false;
1871 // return *this;
1872 // }
1873 // coef = copy(*_cst);
1874 // coef = multiply(coef, *t2.second._coef);
1875 // res.insert(t2.second._sign, *coef, *t2.second._p->first, *t2.second._p->second);
1876 // delete coef;
1877 // }
1878 // for (auto& t2: *f->_lterms) {
1879 // if (t2.second._coef->_is_transposed) {// If the coefficient in front is transposed: a^T.(polynomial term)
1880 // auto be = bexpr(product_, make_shared<func_>(*this), make_shared<func_>(c));
1881 // *this = func_(be);
1882 // _evaluated = false;
1883 // return *this;
1884 // }
1885 // coef = copy(*_cst);
1886 // coef = multiply(coef, *t2.second._coef);
1887 // res.insert(t2.second._sign, *coef, *t2.second._p);
1888 // delete coef;
1889 // }
1890 // if (!f->_cst->is_zero()) {
1891 // coef = copy(*_cst);
1892 // coef = multiply(coef, *f->_cst);
1893 // delete res._cst;
1894 // res._cst = coef;
1895 // if (_cst->is_function()) {
1896 // embed(*(func_*)_cst);
1897 // }
1898 // }
1899 // }
1900 // res.update_dot_dim(*this, c);
1901 // *this = move(res);
1902 // }
1903 // _evaluated = false;
1904 // return *this;
1905 // }
1906 //
1907 // func_& func_::operator/=(const constant_& c){
1908 // if (is_zero()) {
1909 // return *this;
1910 // }
1911 // if (c.is_unit()) {
1912 // return *this;
1913 // }
1914 // if (c.is_zero()) {
1915 // throw invalid_argument("dividing by zero!\n");
1916 // }
1917 // _evaluated = false;
1918 // /* Case where c is a number */
1919 // if (c.is_number()){
1920 // if (!_cst->is_zero()) {
1921 // _cst = divide(_cst, c);
1922 // if (_cst->is_function()) {
1923 // embed(*(func_*)_cst);
1924 // }
1925 // }
1926 // for (auto &pair:*_lterms) {
1927 // pair.second._coef = divide(pair.second._coef, c);
1928 // }
1929 // for (auto &pair:*_qterms) {
1930 // pair.second._coef = divide(pair.second._coef, c);
1931 // }
1932 // for (auto &pair:*_pterms) {
1933 // pair.second._coef = divide(pair.second._coef, c);
1934 // }
1935 // if (_expr) {
1936 // _expr = make_shared<bexpr>(bexpr(div_, make_shared<func_>((*_expr)), make_shared<func_>((c))));
1937 // embed(_expr);
1938 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1939 // _queue->push_back(_expr);
1940 // }
1941 // if (c.is_negative()) {
1942 // reverse_convexity();
1943 // reverse_sign();
1944 // }
1945 // }
1946 // /* Case where the current function is not constant and the other operand is */
1947 // if(!is_constant() && (c.is_param() || (c.is_function() && ((func_*)&c)->is_constant()))) {
1948 // if (!_cst->is_zero()) {
1949 // _cst = divide(_cst, c);
1950 // if (_cst->is_function()) {
1951 // embed(*(func_*)_cst);
1952 // }
1953 // }
1954 // for (auto &pair:*_lterms) {
1955 // pair.second._coef = divide(pair.second._coef, c);
1956 // }
1957 // for (auto &pair:*_qterms) {
1958 // pair.second._coef = divide(pair.second._coef, c);
1959 // }
1960 // for (auto &pair:*_pterms) {
1961 // pair.second._coef = divide(pair.second._coef, c);
1962 // }
1963 // if (c.is_negative()) {
1964 // reverse_convexity();
1965 // reverse_sign();
1966 // }
1967 // if (c.get_all_sign()==unknown_) {
1968 // _all_sign = unknown_;
1969 // if (!_qterms->empty()) {
1970 // _all_convexity = undet_;
1971 // }
1972 // }
1973 // if (_expr) {
1974 // _expr = make_shared<bexpr>(bexpr(div_, make_shared<func_>((*_expr)), make_shared<func_>((c))));
1975 // embed(_expr);
1976 // // _DAG->insert(make_pair<>(_expr->get_str(), _expr));
1977 // _queue->push_back(_expr);
1978 // }
1979 // return *this;
1980 // }
1981 // auto be = bexpr(div_, make_shared<func_>(*this), make_shared<func_>((c)));
1982 // *this = func_(be);
1983 // return *this;
1984 // }
1985 //
1986 // // void func_::update_nb_ind(){
1987 // // _nb_instances = 0;
1988 // // for (auto &p: *_vars) {
1989 // // if (!p.second.first->_is_vector) {
1990 // // _nb_instances = max(_nb_instances, p.second.first->get_dim());
1991 // // }
1992 // // }
1993 // // }
1994 //
1995 //
1996 //
1997
1998
1999
2000
2001
2002
2003 /**
2004 Reverse the convexity property of the current function
2005 */
reverse_convexity()2006 void func_::reverse_convexity(){
2007 if (_all_convexity==convex_) {
2008 _all_convexity=concave_;
2009 }
2010 else if (_all_convexity==concave_) {
2011 _all_convexity=convex_;
2012 }
2013 }
2014
2015
2016
2017 //
update_sign_add(const constant_ & c)2018 void func_::update_sign_add(const constant_& c){
2019 _all_sign = sign_add(_all_sign, c.get_all_sign());
2020 }
2021
update_sign_multiply(const constant_ & c)2022 void func_::update_sign_multiply(const constant_& c){
2023 _all_sign = sign_product(_all_sign, c.get_all_sign());
2024 }
2025 //
2026 // void func_::set_max_dim(const lterm& l){
2027 // // assert(_dim[0] <= l._p->_dim[0] && _dim[1] <= l._p->_dim[1]);
2028 // if(l._p->is_indexed() && l._p->_indices->_ids->size()>1){//This is a multi-instance product
2029 // return;
2030 // }
2031 // if(!l._coef->_is_transposed && !l._p->_is_vector){
2032 // _dim[0] = max(_dim[0], l._p->_dim[0]);
2033 // }
2034 // if(l._coef->_is_vector){
2035 // _dim[0] = l._coef->_dim[0];
2036 // }
2037 // if(l._p->_is_vector){
2038 // _dim[1] = l._p->_dim[1];
2039 // }
2040 // /* Instructions above work for matrix and vector dot products, below we check if it's a component-wise vector,matrix product */
2041 // if(!l._coef->is_matrix() && l._p->is_matrix()){
2042 // _dim[0] = l._p->_dim[0];
2043 // }
2044 // if(l._coef->is_matrix() && !l._p->is_matrix() && l._p->_is_transposed){
2045 // _dim[1] = l._coef->_dim[1];
2046 // }
2047 // }
2048 //
2049 // void func_::set_max_dim(const qterm& q){
2050 // if(q._p->first->is_indexed() && q._p->first->_indices->_ids->size()>1){//This is a multi-instance product
2051 // return;
2052 // }
2053 // if(q._c_p1_transposed){ // q = (coef*p1)^T*p2
2054 // _dim[0] = max(_dim[0],q._p->first->_dim[1]);
2055 // _dim[1] = max(_dim[1],q._p->second->_dim[1]);
2056 // return;
2057 // }
2058 // if(!q._coef->_is_transposed && !q._p->first->_is_vector){
2059 // _dim[0] = max(_dim[0], q._p->first->_dim[0]);
2060 // }
2061 // if (q._p->first->_is_vector) {
2062 // _dim[0] = q._p->first->_dim[0];
2063 // }
2064 // if (!q._coef->_is_transposed && q._p->second->_is_vector) {
2065 // _dim[1] = q._p->second->_dim[1];
2066 // }
2067 // /* Instructions above work for matrix and vector dot products, below we check if it's a component-wise vector,matrix product */
2068 // if(!q._p->first->is_matrix() && q._p->second->is_matrix()){
2069 // _dim[0] = q._p->second->_dim[0];
2070 // }
2071 // if(q._p->first->is_matrix() && !q._p->second->is_matrix() && q._p->second->_is_transposed){
2072 // _dim[1] = q._p->first->_dim[1];
2073 // }
2074 // }
2075 //
2076 //
2077 // void func_::update_sign(const lterm& l){
2078 // Sign sign = get_all_sign(l);
2079 // if (sign==unknown_ || ((_all_sign==non_neg_ || _all_sign==pos_) && sign!=non_neg_ && sign!=pos_)) {
2080 // _all_sign = unknown_;
2081 // }
2082 // else if((_all_sign==non_pos_ || _all_sign==neg_) && sign!=non_pos_ && sign!=neg_){
2083 // _all_sign = unknown_;
2084 // }
2085 // else if(_all_sign==zero_ || _all_sign==pos_ || _all_sign==neg_){// take weaker sign
2086 // _all_sign = sign;
2087 // }
2088 // }
2089 //
2090 // void func_::update_sign(const qterm& q){
2091 // Sign sign = get_all_sign(q);
2092 // if (sign==unknown_ || ((_all_sign==non_neg_ || _all_sign==pos_) && sign!=non_neg_ && sign!=pos_)) {
2093 // _all_sign = unknown_;
2094 // }
2095 // else if((_all_sign==non_pos_ || _all_sign==neg_) && sign!=non_pos_ && sign!=neg_){
2096 // _all_sign = unknown_;
2097 // }
2098 // else if(_all_sign==zero_ || _all_sign==pos_ || _all_sign==neg_){// take weaker sign
2099 // _all_sign = sign;
2100 // }
2101 // }
2102 //
2103 // void func_::update_sign(const pterm& p){
2104 // Sign sign = get_all_sign(p);
2105 // if (sign==unknown_ || ((_all_sign==non_neg_ || _all_sign==pos_) && sign!=non_neg_ && sign!=pos_)) {
2106 // _all_sign = unknown_;
2107 // }
2108 // else if((_all_sign==non_pos_ || _all_sign==neg_) && sign!=non_pos_ && sign!=neg_){
2109 // _all_sign = unknown_;
2110 // }
2111 // else if(_all_sign==zero_ || _all_sign==pos_ || _all_sign==neg_){// take weaker sign
2112 // _all_sign = sign;
2113 // }
2114 // }
2115 //
2116 // void func_::update_convexity(const qterm& q){
2117 // Convexity conv = get_convexity(q);
2118 // if (_all_convexity==undet_ || conv ==undet_ || (_all_convexity==convex_ && conv==concave_) || (_all_convexity==concave_ && conv==convex_)) {
2119 // _all_convexity = undet_;
2120 // }
2121 // else {
2122 // _all_convexity = conv;
2123 // }
2124 // }
2125 //
2126 // bool func_::insert(bool sign, const constant_& coef, const param_& p);
2127
2128
2129 //
2130 // bool func_::insert(bool sign, const constant_& coef, const param_& p1, const param_& p2, bool c_p1_transposed);
2131
2132 // bool func_::insert(const qterm& term)
2133
2134 // bool func_::insert(bool sign, const constant_& coef, const param_& p, int exp)
2135
2136 // bool func_::insert(bool sign, const constant_& coef, const list<pair<shared_ptr<param_>, int>>& l)
2137
2138
2139 //
2140 // // void func_::insert(expr& e){
2141 // // // insert(term._sign, *term._coef, *term._l);
2142 // // auto name = e.get_str();
2143 // //// auto pair_it = _DAG->find(name);
2144 // //
2145 // // _ftype = nlin_;
2146 // //
2147 // // // if (pair_it == _DAG->end()) {
2148 // // // update_sign(e);
2149 // // // update_convexity(e);
2150 // // // _DAG->insert(make_pair<>(name, e));
2151 // // // }
2152 // // // else {
2153 // // // if (pair_it->second._sign == e.sign) {
2154 // // //// pair_it->second._coef = add(pair_it->second._coef, coef);
2155 // // //// if (!pair_it->second._sign) { // both negative
2156 // // //// pair_it->second._sign = true;
2157 // // //// }
2158 // // // }
2159 // // // else{
2160 // // //// pair_it->second._coef = substract(pair_it->second._coef, coef);
2161 // // // }
2162 // // //// if (pair_it->second._coef->is_zero()) {
2163 // // //// if (p1.is_var()) {
2164 // // //// decr_occ_var(s1);
2165 // // //// }
2166 // // //// else {
2167 // // //// decr_occ_param(s1);
2168 // // //// }
2169 // // //// if (p2.is_var()) {
2170 // // //// decr_occ_var(s2);
2171 // // //// }
2172 // // //// else {
2173 // // //// decr_occ_param(s2);
2174 // // //// }
2175 // // //// _qterms->erase(pair_it);
2176 // // //// update_sign();
2177 // // //// update_convexity();
2178 // // //// }
2179 // // //// else {
2180 // // //// update_sign(pair_it->second);
2181 // // //// update_convexity(pair_it->second);
2182 // // //// }
2183 // // //// return false;
2184 // // // }
2185 // // //// _DAG->insert(<#const value_type &__v#>)
2186 // // }
2187 //
2188 //
2189 //
2190 //
2191 // func_ cos(const constant_& c){
2192 // func_ res;
2193 // if (c.is_number()) {
2194 // res._val->resize(1, std::cos(t_eval(&c)));
2195 // }
2196 // else {
2197 // res._evaluated = false;
2198 // }
2199 // res._expr = make_shared<uexpr>(uexpr(cos_, make_shared<func_>((c))));
2200 // res.embed(res._expr);
2201 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2202 // res._queue->push_back(res._expr);
2203 // if (!res._vars->empty()) {
2204 // res._ftype = nlin_;
2205 // }
2206 // res._dim[0] = c._dim[0];
2207 // res._dim[1] = c._dim[1];
2208 // return res;
2209 // };
2210 //
2211 // func_ cos(constant_&& c){
2212 // func_ res;
2213 // if (c.is_number()) {
2214 // res._val->resize(1, std::cos(t_eval(&c)));
2215 // }
2216 // else {
2217 // res._evaluated = false;
2218 // }
2219 // res._expr = make_shared<uexpr>(uexpr(cos_, make_shared<func_>((move(c)))));
2220 // res.embed(res._expr);
2221 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2222 // res._queue->push_back(res._expr);
2223 // if (!res._vars->empty()) {
2224 // res._ftype = nlin_;
2225 // }
2226 // res._dim[0] = c._dim[0];
2227 // res._dim[1] = c._dim[1];
2228 // return res;
2229 // };
2230 //
2231 //
2232 // func_ sin(const constant_& c){
2233 // func_ res;
2234 // if (c.is_number()) {
2235 // res._val->resize(1, std::sin(t_eval(&c)));
2236 // }
2237 // else {
2238 // res._evaluated = false;
2239 // }
2240 // res._expr = make_shared<uexpr>(uexpr(sin_, make_shared<func_>((c))));
2241 // res.embed(res._expr);
2242 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2243 // res._queue->push_back(res._expr);
2244 // if (!res._vars->empty()) {
2245 // res._ftype = nlin_;
2246 // }
2247 // res._dim[0] = c._dim[0];
2248 // res._dim[1] = c._dim[1];
2249 // return res;
2250 // };
2251 //
2252 // func_ sin(constant_&& c){
2253 // func_ res;
2254 // if (c.is_number()) {
2255 // res._val->resize(1, std::sin(t_eval(&c)));
2256 // }
2257 // else {
2258 // res._evaluated = false;
2259 // }
2260 // res._expr = make_shared<uexpr>(uexpr(sin_, make_shared<func_>((move(c)))));
2261 // res.embed(res._expr);
2262 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2263 // res._queue->push_back(res._expr);
2264 // if (!res._vars->empty()) {
2265 // res._ftype = nlin_;
2266 // }
2267 // res._dim[0] = c._dim[0];
2268 // res._dim[1] = c._dim[1];
2269 // return res;
2270 // };
2271 //
2272 //
2273 // func_ sqrt(const constant_& c){
2274 // func_ res;
2275 // if (c.is_number()) {
2276 // res._val->resize(1, std::sqrt(t_eval(&c)));
2277 // }
2278 // else {
2279 // res._evaluated = false;
2280 // }
2281 // auto exp = make_shared<uexpr>(uexpr(sqrt_,make_shared<func_>((c))));
2282 // res._expr = exp;
2283 // res.embed(res._expr);
2284 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2285 // res._queue->push_back(res._expr);
2286 // if (!res._vars->empty()) {
2287 // res._ftype = nlin_;
2288 // if (exp->_son->is_concave()) {
2289 // res._all_convexity = concave_;
2290 // }
2291 // else {
2292 // res._all_convexity = undet_;
2293 // }
2294 // }
2295 // res._dim[0] = c._dim[0];
2296 // res._dim[1] = c._dim[1];
2297 // return res;
2298 // };
2299 //
2300 // func_ sqrt(constant_&& c){
2301 // func_ res;
2302 // if (c.is_number()) {
2303 // res._val->resize(1, std::sqrt(t_eval(&c)));
2304 // }
2305 // else {
2306 // res._evaluated = false;
2307 // }
2308 // auto exp = make_shared<uexpr>(uexpr(sqrt_, make_shared<func_>((move(c)))));
2309 // res._expr = exp;
2310 // res.embed(res._expr);
2311 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2312 // res._queue->push_back(res._expr);
2313 // if (!res._vars->empty()) {
2314 // res._ftype = nlin_;
2315 // if (exp->_son->is_concave()) {
2316 // res._all_convexity = concave_;
2317 // }
2318 // else {
2319 // res._all_convexity = undet_;
2320 // }
2321 // }
2322 // res._dim[0] = c._dim[0];
2323 // res._dim[1] = c._dim[1];
2324 // return res;
2325 // };
2326 //
2327 //
2328 // func_ expo(const constant_& c){
2329 // func_ res;
2330 // if (c.is_number()) {
2331 // res._val->resize(1, std::exp(t_eval(&c)));
2332 // }
2333 // else {
2334 // res._evaluated = false;
2335 // }
2336 // res._is_vector = c._is_vector;
2337 // // res._is_matrix = c._is_matrix;
2338 // res._is_transposed = c._is_transposed;
2339 // res._dim[0] = c._dim[0];
2340 // res._dim[1] = c._dim[1];
2341 // auto exp = make_shared<uexpr>(uexpr(exp_, make_shared<func_>((c))));
2342 // res._expr = exp;
2343 // res.embed(res._expr);
2344 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2345 // res._queue->push_back(res._expr);
2346 // if (!res._vars->empty()) {
2347 // res._ftype = nlin_;
2348 // if (exp->_son->is_convex()) {
2349 // res._all_convexity = convex_;
2350 // }
2351 // else {
2352 // res._all_convexity = undet_;
2353 // }
2354 // }
2355 // return res;
2356 // };
2357 //
2358 //
2359 // func_ expo(constant_&& c){
2360 // func_ res;
2361 // if (c.is_number()) {
2362 // res._val->resize(1, std::exp(t_eval(&c)));
2363 // }
2364 // else {
2365 // res._evaluated = false;
2366 // }
2367 // res._is_vector = c._is_vector;
2368 // // res._is_matrix = c._is_matrix;
2369 // res._is_transposed = c._is_transposed;
2370 // res._dim[0] = c._dim[0];
2371 // res._dim[1] = c._dim[1];
2372 // auto exp = make_shared<uexpr>(uexpr(exp_, make_shared<func_>((move(c)))));
2373 // res._expr = exp;
2374 // res.embed(res._expr);
2375 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2376 // res._queue->push_back(res._expr);
2377 // if (!res._vars->empty()) {
2378 // res._ftype = nlin_;
2379 // if (exp->_son->is_convex()) {
2380 // res._all_convexity = convex_;
2381 // }
2382 // else {
2383 // res._all_convexity = undet_;
2384 // }
2385 // }
2386 // return res;
2387 // };
2388 //
2389 //
2390 //
2391 // func_ log(const constant_& c){
2392 // func_ res;
2393 // if (c.is_number()) {
2394 // res._val->resize(1, std::log(t_eval(&c)));
2395 // }
2396 // else {
2397 // res._evaluated = false;
2398 // }
2399 //
2400 //
2401 // auto exp = make_shared<uexpr>(uexpr(log_, make_shared<func_>((c))));
2402 // res._expr = exp;
2403 // res.embed(res._expr);
2404 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2405 // res._queue->push_back(res._expr);
2406 //
2407 // if (!res._vars->empty()) {
2408 // res._ftype = nlin_;
2409 // if (exp->_son->is_concave()) {
2410 // res._all_convexity = concave_;
2411 // }
2412 // else {
2413 // res._all_convexity = undet_;
2414 // }
2415 // }
2416 // res._dim[0] = c._dim[0];
2417 // res._dim[1] = c._dim[1];
2418 // return res;
2419 // };
2420 //
2421 // func_ log(constant_&& c){
2422 // func_ res;
2423 // if (c.is_number()) {
2424 // res._val->resize(1, std::log(t_eval(&c)));
2425 // }
2426 // else {
2427 // res._evaluated = false;
2428 // }
2429 // auto exp = make_shared<uexpr>(uexpr(log_, make_shared<func_>((move(c)))));
2430 // res._expr = exp;
2431 // res.embed(res._expr);
2432 // // res._DAG->insert(make_pair<>(res._expr->get_str(), res._expr));
2433 // res._queue->push_back(res._expr);
2434 // if (!res._vars->empty()) {
2435 // res._ftype = nlin_;
2436 // if (exp->_son->is_concave()) {
2437 // res._all_convexity = concave_;
2438 // }
2439 // else {
2440 // res._all_convexity = undet_;
2441 // }
2442 // }
2443 // res._dim[0] = c._dim[0];
2444 // res._dim[1] = c._dim[1];
2445 // return res;
2446 // };
2447 //
2448 // constant_* copy(const constant_& c2){/**< Return a copy c2 detecting the right class, i.e., constant<>, param<>, uexpr , bexpr or func. WARNING: allocates memory! */
2449 //
2450 // switch (c2.get_type()) {
2451 // case binary_c: {
2452 // return new constant<bool>(*(constant<bool>*)(&c2));
2453 // break;
2454 // }
2455 // case short_c: {
2456 // return new constant<short>(*(constant<short>*)(&c2));
2457 // break;
2458 // }
2459 // case integer_c: {
2460 // return new constant<int>(*(constant<int>*)(&c2));
2461 // break;
2462 // }
2463 // case float_c: {
2464 // return new constant<float>(*(constant<float>*)(&c2));
2465 // break;
2466 // }
2467 // case double_c: {
2468 // return new constant<double>(*(constant<double>*)(&c2));
2469 // break;
2470 // }
2471 // case long_c: {
2472 // return new constant<long double>(*(constant<long double>*)(&c2));
2473 // break;
2474 // }
2475 // case complex_c: {
2476 // return new constant<Cpx>(*(constant<Cpx>*)(&c2));
2477 // break;
2478 // }
2479 // case par_c:{
2480 // auto p_c2 = (param_*)(&c2);
2481 // switch (p_c2->get_intype()) {
2482 // case binary_:
2483 // return new param<bool>(*(param<bool>*)p_c2);
2484 // break;
2485 // case short_:
2486 // return new param<short>(*(param<short>*)p_c2);
2487 // break;
2488 // case integer_:
2489 // return new param<int>(*(param<int>*)p_c2);
2490 // break;
2491 // case float_:
2492 // return new param<float>(*(param<float>*)p_c2);
2493 // break;
2494 // case double_:
2495 // return new param<double>(*(param<double>*)p_c2);
2496 // break;
2497 // case long_:
2498 // return new param<long double>(*(param<long double>*)p_c2);
2499 // break;
2500 // case complex_:
2501 // return new param<Cpx>(*(param<Cpx>*)p_c2);
2502 // break;
2503 // default:
2504 // break;
2505 // }
2506 // break;
2507 // }
2508 // case uexp_c: {
2509 // return new uexpr(*(uexpr*)&c2);
2510 // break;
2511 // }
2512 // case bexp_c: {
2513 // return new bexpr(*(bexpr*)&c2);
2514 // break;
2515 // }
2516 // case var_c:{
2517 // auto p_c2 = (param_*)(&c2);
2518 // switch (p_c2->get_intype()) {
2519 // case binary_:
2520 // return new var<bool>(*(var<bool>*)p_c2);
2521 // break;
2522 // case short_:
2523 // return new var<short>(*(var<short>*)p_c2);
2524 // break;
2525 // case integer_:
2526 // return new var<int>(*(var<int>*)p_c2);
2527 // break;
2528 // case float_:
2529 // return new var<float>(*(var<float>*)p_c2);
2530 // break;
2531 // case double_:
2532 // return new var<double>(*(var<double>*)p_c2);
2533 // break;
2534 // case long_:
2535 // return new var<long double>(*(var<long double>*)p_c2);
2536 // break;
2537 // case complex_:
2538 // return new var<Cpx>(*(var<Cpx>*)p_c2);
2539 // break;
2540 // default:
2541 // break;
2542 // }
2543 // break;
2544 // }
2545 // case func_c: {
2546 // return new func_(c2);
2547 // break;
2548 // }
2549 //
2550 // default:
2551 // break;
2552 // }
2553 // return nullptr;
2554 // }
2555 //
2556 // bool equals(const constant_* c1, const constant_* c2){/**< Checks if c2 equals c1 detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2557 // if ((!c1 && !c2) || (c1==c2)) {
2558 // return true;
2559 // }
2560 // if ((c1 && !c2) || (!c1 && c2)) {
2561 // return false;
2562 // }
2563 //
2564 // switch (c2->get_type()) {
2565 // case binary_c: {
2566 // return (c1->is_binary() && *(constant<bool>*)c1 == *(constant<bool>*)c2);
2567 // break;
2568 // }
2569 // case short_c: {
2570 // return (c1->is_short() && *(constant<short>*)c1 == *(constant<short>*)c2);
2571 // break;
2572 // }
2573 // case integer_c: {
2574 // return (c1->is_integer() && *(constant<int>*)c1 == *(constant<int>*)c2);
2575 // break;
2576 // }
2577 // case float_c: {
2578 // return (c1->is_float() && *(constant<float>*)c1 == *(constant<float>*)c2);
2579 // break;
2580 // }
2581 // case double_c: {
2582 // return (c1->is_double() && *(constant<double>*)c1 == *(constant<double>*)c2);
2583 // break;
2584 // }
2585 // case long_c: {
2586 // return (c1->is_long() && *(constant<long double>*)c1 == *(constant<long double>*)c2);
2587 // break;
2588 // }
2589 // case complex_c: {
2590 // return (c1->is_complex() && *(constant<Cpx>*)c1 == *(constant<Cpx>*)c2);
2591 // break;
2592 // }
2593 // case par_c:{
2594 // return (c1->is_param() && *(param_ *)c1 == *(param_ *)c2);
2595 // break;
2596 // }
2597 // case uexp_c: {
2598 // return (c1->is_uexpr() && *(uexpr *)c1 == *(uexpr *)c2);
2599 // break;
2600 // }
2601 // case bexp_c: {
2602 // return (c1->is_bexpr() && *(bexpr *)c1 == *(bexpr *)c2);
2603 // break;
2604 // }
2605 // case var_c:{
2606 // return (c1->is_var() && *(param_ *)c1 == *(param_ *)c2);
2607 // break;
2608 // }
2609 // case func_c:{
2610 // if (c1->is_function()){
2611 // auto f1 = (func_ *)c1;
2612 // auto f2 = (func_ *)c2;
2613 // return *f1==*f2;
2614 // }
2615 // break;
2616 // }
2617 // default:
2618 // break;
2619 // }
2620 // return false;
2621 // }
2622 //
2623 // void set_val(param_* p, double val, size_t i){ /**< Polymorphic set_val */
2624 // switch (p->get_intype()) {
2625 // case binary_:
2626 // ((param<bool>*)p)->set_val(i, val);
2627 // break;
2628 // case short_:
2629 // ((param<short>*)p)->set_val(i, val);
2630 // break;
2631 // case integer_:
2632 // ((param<int>*)p)->set_val(i, val);
2633 // break;
2634 // case float_:
2635 // ((param<float>*)p)->set_val(i, val);
2636 // break;
2637 // case double_:
2638 // ((param<double>*)p)->set_val(i, val);
2639 // break;
2640 // case long_:
2641 // ((param<long double>*)p)->set_val(i, val);
2642 // break;
2643 // default:
2644 // throw invalid_argument("Polymorphic set_val only works for arithmetic types");
2645 // break;
2646 // }
2647 // }
2648 //
2649 //
2650 // /* Polymorphic functions */
2651 //
2652 // //void reverse_sign(constant_* c){ /**< Reverses the sign of the constant. */
2653 // // switch (c->get_type()) {
2654 // // case binary_c: {
2655 // // ((constant<bool>*)c)->set_val(!((constant<bool>*)c)->eval());
2656 // // break;
2657 // // }
2658 // // case short_c: {
2659 // // ((constant<short>*)c)->set_val(-1*((constant<short>*)c)->eval());
2660 // // break;
2661 // // }
2662 // // case integer_c: {
2663 // // ((constant<int>*)c)->set_val(-1*((constant<int>*)c)->eval());
2664 // // break;
2665 // // }
2666 // // case float_c: {
2667 // // ((constant<float>*)c)->set_val(-1*((constant<float>*)c)->eval());
2668 // // break;
2669 // // }
2670 // // case double_c: {
2671 // // ((constant<double>*)c)->set_val(-1*((constant<double>*)c)->eval());
2672 // // break;
2673 // // }
2674 // // case long_c: {
2675 // // ((constant<long double>*)c)->set_val(-1*((constant<long double>*)c)->eval());
2676 // // break;
2677 // // }
2678 // // default:
2679 // // throw invalid_argument("Cannot reverse sign of non-numeric constant");
2680 // // break;
2681 // // }
2682 // //
2683 // //}
2684 //
2685 // void poly_print(const constant_* c){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2686 // cout << poly_to_str(c);
2687 // }
2688 //
2689 // string matrix_to_str(const constant_* c){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2690 //
2691 // if (!c) {
2692 // return "null";
2693 // }
2694 // if(!c->is_matrix()){
2695 // return poly_to_str(c);
2696 // }
2697 // switch (c->get_type()) {
2698 // case par_c:{
2699 // return ((param_*)(c))->get_name(false,false);
2700 // break;
2701 // }
2702 // case uexp_c: {
2703 // return ((uexpr*)c)->get_str();
2704 // break;
2705 // }
2706 // case bexp_c: {
2707 // return ((bexpr*)c)->get_str();
2708 // break;
2709 // }
2710 // case var_c: {
2711 // return ((param_*)(c))->get_name(false,false);
2712 // break;
2713 // }
2714 // case func_c: {
2715 // return ((func_*)c)->to_str();
2716 // break;
2717 // }
2718 // default:
2719 // break;
2720 // }
2721 // return "null";
2722 // }
2723 //
2724 // string poly_to_str(const constant_* c){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2725 //
2726 // if (!c) {
2727 // return "null";
2728 // }
2729 // switch (c->get_type()) {
2730 // case binary_c: {
2731 // return ((constant<bool>*)(c))->to_str();
2732 // break;
2733 // }
2734 // case short_c: {
2735 // return ((constant<short>*)(c))->to_str();
2736 // break;
2737 // }
2738 // case integer_c: {
2739 // return ((constant<int>*)(c))->to_str();
2740 // break;
2741 // }
2742 // case float_c: {
2743 // return ((constant<float>*)(c))->to_str();
2744 // break;
2745 // }
2746 // case double_c: {
2747 // return ((constant<double>*)(c))->to_str();
2748 // break;
2749 // }
2750 // case long_c: {
2751 // return ((constant<long double>*)(c))->to_str();
2752 // break;
2753 // }
2754 // case complex_c: {
2755 // return ((constant<Cpx>*)(c))->to_str();
2756 // break;
2757 // }
2758 // case par_c:{
2759 // return ((param_*)(c))->get_name(false,false);
2760 // break;
2761 // }
2762 // case uexp_c: {
2763 // return ((uexpr*)c)->get_str();
2764 // break;
2765 // }
2766 // case bexp_c: {
2767 // return ((bexpr*)c)->get_str();
2768 // break;
2769 // }
2770 // case var_c: {
2771 // return ((param_*)(c))->get_name(false,false);
2772 // break;
2773 // }
2774 // case func_c: {
2775 // return ((func_*)c)->to_str();
2776 // break;
2777 // }
2778 // default:
2779 // break;
2780 // }
2781 // return "null";
2782 // }
2783 //
2784 // // string poly_to_str(const constant_* c, size_t inst1, size_t inst2){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2785 // //
2786 // // if (!c) {
2787 // // return "null";
2788 // // }
2789 // // switch (c->get_type()) {
2790 // // case binary_c: {
2791 // // return ((constant<bool>*)(c))->to_str();
2792 // // break;
2793 // // }
2794 // // case short_c: {
2795 // // return ((constant<short>*)(c))->to_str();
2796 // // break;
2797 // // }
2798 // // case integer_c: {
2799 // // return ((constant<int>*)(c))->to_str();
2800 // // break;
2801 // // }
2802 // // case float_c: {
2803 // // return ((constant<float>*)(c))->to_str();
2804 // // break;
2805 // // }
2806 // // case double_c: {
2807 // // return ((constant<double>*)(c))->to_str();
2808 // // break;
2809 // // }
2810 // // case long_c: {
2811 // // return ((constant<long double>*)(c))->to_str();
2812 // // break;
2813 // // }
2814 // // case complex_c: {
2815 // // return ((constant<Cpx>*)(c))->to_str();
2816 // // break;
2817 // // }
2818 // // case par_c:{
2819 // // auto p_c = (param_*)(c);
2820 // // switch (p_c->get_intype()) {
2821 // // case binary_:
2822 // // return ((param<bool>*)p_c)->to_str(inst1,inst2);
2823 // // break;
2824 // // case short_:
2825 // // return ((param<short>*)p_c)->to_str(inst1,inst2);
2826 // // break;
2827 // // case integer_:
2828 // // return ((param<int>*)p_c)->to_str(inst1,inst2);
2829 // // break;
2830 // // case float_:
2831 // // return ((param<float>*)p_c)->to_str(inst1,inst2);
2832 // // break;
2833 // // case double_:
2834 // // return ((param<double>*)p_c)->to_str(inst1,inst2);
2835 // // break;
2836 // // case long_:
2837 // // return ((param<long double>*)p_c)->to_str(inst1,inst2);
2838 // // break;
2839 // // case complex_:
2840 // // return ((param<Cpx>*)p_c)->to_str(inst1,inst2);
2841 // // break;
2842 // // default:
2843 // // break;
2844 // // }
2845 // // break;
2846 // // }
2847 // // case uexp_c: {
2848 // // return ((uexpr*)c)->to_str(inst2);
2849 // // break;
2850 // // }
2851 // // case bexp_c: {
2852 // // return ((bexpr*)c)->to_str(inst2);
2853 // // break;
2854 // // }
2855 // // case var_c: {
2856 // // auto p_c = (param_*)(c);
2857 // // switch (p_c->get_intype()) {
2858 // // case binary_:
2859 // // return ((var<bool>*)p_c)->get_name(inst1,inst2);
2860 // // break;
2861 // // case short_:
2862 // // return ((var<short>*)p_c)->get_name(inst1,inst2);
2863 // // break;
2864 // // case integer_:
2865 // // return ((var<int>*)p_c)->get_name(inst1,inst2);
2866 // // break;
2867 // // case float_:
2868 // // return ((var<float>*)p_c)->get_name(inst1,inst2);
2869 // // break;
2870 // // case double_:
2871 // // return ((var<double>*)p_c)->get_name(inst1,inst2);
2872 // // break;
2873 // // case long_:
2874 // // return ((var<long double>*)p_c)->get_name(inst1,inst2);
2875 // // break;
2876 // // case complex_:
2877 // // return ((var<Cpx>*)p_c)->get_name(inst1,inst2);
2878 // // break;
2879 // // default:
2880 // // break;
2881 // // }
2882 // // break;
2883 // // }
2884 // //
2885 // // case func_c: {
2886 // // return ((func_*)c)->to_str(inst2);
2887 // // break;
2888 // // }
2889 // // default:
2890 // // break;
2891 // // }
2892 // // return "null";
2893 // // }
2894 // //
2895 // //
2896 // // string poly_to_str(const constant_* c, size_t inst){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
2897 // //
2898 // // if (!c) {
2899 // // return "null";
2900 // // }
2901 // // switch (c->get_type()) {
2902 // // case binary_c: {
2903 // // return ((constant<bool>*)(c))->to_str();
2904 // // break;
2905 // // }
2906 // // case short_c: {
2907 // // return ((constant<short>*)(c))->to_str();
2908 // // break;
2909 // // }
2910 // // case integer_c: {
2911 // // return ((constant<int>*)(c))->to_str();
2912 // // break;
2913 // // }
2914 // // case float_c: {
2915 // // return ((constant<float>*)(c))->to_str();
2916 // // break;
2917 // // }
2918 // // case double_c: {
2919 // // return ((constant<double>*)(c))->to_str();
2920 // // break;
2921 // // }
2922 // // case long_c: {
2923 // // return ((constant<long double>*)(c))->to_str();
2924 // // break;
2925 // // }
2926 // // case complex_c: {
2927 // // return ((constant<Cpx>*)(c))->to_str();
2928 // // break;
2929 // // }
2930 // // case par_c:{
2931 // // auto p_c = (param_*)(c);
2932 // // switch (p_c->get_intype()) {
2933 // // case binary_:
2934 // // return ((param<bool>*)p_c)->to_str(inst);
2935 // // break;
2936 // // case short_:
2937 // // return ((param<short>*)p_c)->to_str(inst);
2938 // // break;
2939 // // case integer_:
2940 // // return ((param<int>*)p_c)->to_str(inst);
2941 // // break;
2942 // // case float_:
2943 // // return ((param<float>*)p_c)->to_str(inst);
2944 // // break;
2945 // // case double_:
2946 // // return ((param<double>*)p_c)->to_str(inst);
2947 // // break;
2948 // // case long_:
2949 // // return ((param<long double>*)p_c)->to_str(inst);
2950 // // break;
2951 // // case complex_:
2952 // // return ((param<Cpx>*)p_c)->to_str(inst);
2953 // // break;
2954 // // default:
2955 // // break;
2956 // // }
2957 // // break;
2958 // // }
2959 // // case uexp_c: {
2960 // // return ((uexpr*)c)->to_str(inst);
2961 // // break;
2962 // // }
2963 // // case bexp_c: {
2964 // // return ((bexpr*)c)->to_str(inst);
2965 // // break;
2966 // // }
2967 // // case var_c: {
2968 // // auto p_c = (param_*)(c);
2969 // // switch (p_c->get_intype()) {
2970 // // case binary_:
2971 // // return ((var<bool>*)p_c)->get_name(inst);
2972 // // break;
2973 // // case short_:
2974 // // return ((var<short>*)p_c)->get_name(inst);
2975 // // break;
2976 // // case integer_:
2977 // // return ((var<int>*)p_c)->get_name(inst);
2978 // // break;
2979 // // case float_:
2980 // // return ((var<float>*)p_c)->get_name(inst);
2981 // // break;
2982 // // case double_:
2983 // // return ((var<double>*)p_c)->get_name(inst);
2984 // // break;
2985 // // case long_:
2986 // // return ((var<long double>*)p_c)->get_name(inst);
2987 // // break;
2988 // // case complex_:
2989 // // return ((var<Cpx>*)p_c)->get_name(inst);
2990 // // break;
2991 // // default:
2992 // // break;
2993 // // }
2994 // // break;
2995 // // }
2996 // // case func_c: {
2997 // // return ((func_*)c)->to_str(inst);
2998 // // break;
2999 // // }
3000 // // default:
3001 // // break;
3002 // // }
3003 // // return "null";
3004 // // }
3005 //
3006 //
3007 //
3008 //
3009 // func_ operator+(const constant_& c1, const constant_& c2){
3010 // return func_(c1) += c2;
3011 // }
3012 //
3013 // // func_ operator+(func_&& f, const constant_& c){
3014 // // return f += c;
3015 // // }
3016 // //
3017 // // func_ operator+(const constant_& c, func_&& f){
3018 // // return f += c;
3019 // // }
3020 //
3021 //
3022 // func_ operator-(const constant_& c1, const constant_& c2){
3023 // return func_(c1) -= c2;
3024 // }
3025 //
3026 // // func_ operator-(func_&& f, const constant_& c){
3027 // // return f -= c;
3028 // // }
3029 // //
3030 // // func_ operator-(const constant_& c, func_&& f){
3031 // // return (f *= -1) += c;
3032 // // }
3033 // //
3034 //
3035 // func_ operator*(const constant_& c1, const constant_& c2){// Rewrite this to change res after the multiplication is done, make sure both vars are now vecs.
3036 // if(c1.is_number()) {
3037 // if (c1._is_transposed && !c1.is_matrix()) {//TODO check when matrix
3038 // auto new_c2 = copy(c2);
3039 // new_c2->_is_vector = true;
3040 // auto res = func_(c1);
3041 // res._dim[1] = c2._dim[0];
3042 // res *= move(*new_c2);
3043 // res._is_vector = false;
3044 // res._is_transposed = false;
3045 // delete new_c2;
3046 // return res;
3047 // }
3048 // else {
3049 // return func_(c1) *= c2;
3050 // }
3051 // }
3052 // else {
3053 //// if (c1._is_transposed) {//TODO check when matrix
3054 //// auto new_c2 = copy(c2);
3055 //// new_c2->_is_vector = true;
3056 //// return func_(c1) *= *new_c2;
3057 //// }
3058 //// if(c1.is_number()) {
3059 //// if (c1._is_transposed && !c1.is_matrix()) {//TODO check when matrix
3060 //// auto new_c2 = copy(c2);
3061 //// new_c2->_is_vector = true;
3062 //// auto res = func_(c1);
3063 //// res._dim[1] = c2._dim[0];
3064 //// res *= move(*new_c2);
3065 //// res._is_vector = false;
3066 //// res._is_transposed = false;
3067 //// delete new_c2;
3068 //// return res;
3069 //// }
3070 //// else {
3071 //// return func_(c1) *= c2;
3072 //// }
3073 //// }
3074 //// else {
3075 // // if (c1._is_transposed) {
3076 // //// auto new_c2 = copy(c2);
3077 // //// new_c2->_is_vector = true;
3078 // //// auto res = func_(c1) *= move(*new_c2);
3079 // //// auto res = ;
3080 // //// delete new_c2;
3081 // // return func_(c1) *= c2;
3082 // // }
3083 // // if (c2.is_var()) {
3084 // // return func_(c2) *= c1;
3085 // // }
3086 // // if(c1.is_function()){
3087 // // return (*(func_*)&c1) *= c2;
3088 // // }
3089 // // else {
3090 // return func_(c1) *= c2;
3091 // // }
3092 // }
3093 // }
3094 //
3095 // // func_ operator*(func_&& f, const constant_& c){
3096 // // return f *= c;
3097 // // }
3098 // //
3099 // // func_ operator*(const constant_& c, func_&& f){
3100 // // return f *= c;
3101 // // }
3102 //
3103 // func_ operator/(const constant_& c1, const constant_& c2){
3104 // return func_(c1) /= c2;
3105 // }
3106 //
3107 // // func_ operator/(func_&& f, const constant_& c){
3108 // // return f /= c;
3109 // // }
3110 //
3111 // //
3112 // //func_ operator+(const func_& f1, const func_& f2){
3113 // // return func_(f1) += f2;
3114 // //}
3115 // //
3116 // //func_ operator+(func_&& f1, const func_& f2){
3117 // // return f1 += f2;
3118 // //}
3119 // //
3120 // //func_ operator+(const func_& f1, func_&& f2){
3121 // // return f2 += f1;
3122 // //}
3123 // //
3124 // //
3125 // //func_ operator-(const func_& f1, const func_& f2){
3126 // // return func_(f1) -= f2;
3127 // //}
3128 // //
3129 // //func_ operator-(func_&& f1, const func_& f2){
3130 // // return f1 -= f2;
3131 // //}
3132 // //
3133 // //func_ operator-(const func_& f1, func_&& f2){
3134 // // return (f2 *= -1) += f1;
3135 // //}
3136 // //
3137 // constant_* add(constant_* c1, const func_& f){
3138 // switch (c1->get_type()) {
3139 // case binary_c: {
3140 // auto res = new func_(f);
3141 // *res += ((constant<bool>*)c1)->eval();
3142 // delete c1;
3143 // return c1 = res;
3144 // break;
3145 // }
3146 // case short_c: {
3147 // auto res = new func_(f);
3148 // *res += ((constant<short>*)c1)->eval();
3149 // delete c1;
3150 // return c1 = res;
3151 // break;
3152 // }
3153 // case integer_c: {
3154 // auto res = new func_(f);
3155 // *res += ((constant<int>*)c1)->eval();
3156 // delete c1;
3157 // return c1 = res;
3158 // break;
3159 // }
3160 // case float_c: {
3161 // auto res = new func_(f);
3162 // *res += ((constant<float>*)c1)->eval();
3163 // delete c1;
3164 // return c1 = res;
3165 // break;
3166 // }
3167 // case double_c: {
3168 // auto res = new func_(f);
3169 // *res += ((constant<double>*)c1)->eval();
3170 // delete c1;
3171 // return c1 = res;
3172 // break;
3173 // }
3174 // case long_c: {
3175 // auto res = new func_(f);
3176 // *res += ((constant<long double>*)c1)->eval();
3177 // delete c1;
3178 // return c1 = res;
3179 // break;
3180 // }
3181 // case complex_c:{
3182 // auto res = new func_(f);
3183 // *res += *((constant<Cpx>*)c1);
3184 // delete c1;
3185 // return c1 = res;
3186 // break;
3187 // }
3188 // case par_c:{
3189 // auto res = new func_(f);
3190 // if (f.is_constant()) {
3191 // res->insert(true, constant<double>(1), *(param_*)c1);
3192 // }
3193 // else{
3194 // auto cst = res->get_cst();
3195 // cst = add(cst, *c1);
3196 // }
3197 // delete c1;
3198 // return c1 = res;
3199 // break;
3200 // }
3201 // case var_c:{
3202 // auto res = new func_(f);
3203 // delete c1;
3204 // res->insert(true, constant<double>(1), *(param_*)c1);
3205 // return c1 = res;
3206 // break;
3207 // }
3208 // case uexp_c: {
3209 // auto res = new bexpr(plus_, make_shared<func_>((*c1)), make_shared<func_>(f));
3210 // delete c1;
3211 // c1 = (constant_*)res;
3212 // return c1;
3213 // break;
3214 // }
3215 // case bexp_c: {
3216 // auto res = new bexpr(plus_, make_shared<func_>((*c1)), make_shared<func_>(f));
3217 // delete c1;
3218 // c1 = (constant_*)res;
3219 // return c1;
3220 // break;
3221 // }
3222 // case func_c: {
3223 // auto res = new func_(f);
3224 // *res += *((func_*)c1);
3225 // delete c1;
3226 // return c1 = res;
3227 // break; }
3228 //
3229 // default:
3230 // break;
3231 // }
3232 // return c1;
3233 // }
3234 //
3235 // constant_* add(constant_* c1, const constant_& c2){ /**< adds c2 to c1, updates its type and returns the result **/
3236 // switch (c2.get_type()) {
3237 // case binary_c: {
3238 // return add(c1, *(constant<bool>*)&c2);
3239 // break;
3240 // }
3241 // case short_c: {
3242 // return add(c1, *(constant<short>*)&c2);
3243 // break;
3244 // }
3245 // case integer_c: {
3246 // return add(c1, *(constant<int>*)&c2);
3247 // break;
3248 // }
3249 // case float_c: {
3250 // return add(c1, *(constant<float>*)&c2);
3251 // break;
3252 // }
3253 // case double_c: {
3254 // return add(c1, *(constant<double>*)&c2);
3255 // break;
3256 // }
3257 // case long_c: {
3258 // return add(c1, *(constant<long double>*)&c2);
3259 // break;
3260 // }
3261 // case complex_c: {
3262 // return add(c1, *(constant<Cpx>*)&c2);
3263 // break;
3264 // }
3265 // case par_c:{
3266 // auto pc2 = (param_*)(&c2);
3267 // switch (pc2->get_intype()) {
3268 // case binary_:
3269 // return add(c1, *(param<bool>*)pc2);
3270 // break;
3271 // case short_:
3272 // return add(c1, *(param<short>*)pc2);
3273 // break;
3274 // case integer_:
3275 // return add(c1, *(param<int>*)pc2);
3276 // break;
3277 // case float_:
3278 // return add(c1, *(param<float>*)pc2);
3279 // break;
3280 // case double_:
3281 // return add(c1, *(param<double>*)pc2);
3282 // break;
3283 // case long_:
3284 // return add(c1, *(param<long double>*)pc2);
3285 // break;
3286 // case complex_:
3287 // return add(c1, *(param<Cpx>*)pc2);
3288 // break;
3289 // default:
3290 // break;
3291 // }
3292 // break;
3293 // }
3294 // case var_c:{
3295 // auto pc2 = (param_*)(&c2);
3296 // switch (pc2->get_intype()) {
3297 // case binary_:
3298 // return add(c1, *(var<bool>*)pc2);
3299 // break;
3300 // case short_:
3301 // return add(c1, *(var<short>*)pc2);
3302 // break;
3303 // case integer_:
3304 // return add(c1, *(var<int>*)pc2);
3305 // break;
3306 // case float_:
3307 // return add(c1, *(var<float>*)pc2);
3308 // break;
3309 // case double_:
3310 // return add(c1, *(var<double>*)pc2);
3311 // break;
3312 // case long_:
3313 // return add(c1, *(var<long double>*)pc2);
3314 // break;
3315 // case complex_:
3316 // return add(c1, *(var<Cpx>*)pc2);
3317 // break;
3318 // default:
3319 // break;
3320 // }
3321 // break;
3322 // }
3323 // case uexp_c: {
3324 // auto res = new bexpr(plus_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3325 // delete c1;
3326 // c1 = (constant_*)res;
3327 // return c1;
3328 // break;
3329 // }
3330 // case bexp_c: {
3331 // auto res = new bexpr(plus_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3332 // delete c1;
3333 // c1 = (constant_*)res;
3334 // return c1;
3335 // break;
3336 // }
3337 // case func_c: {
3338 // return add(c1, *(func_*)&c2);
3339 // break;
3340 // }
3341 // default:
3342 // break;
3343 // }
3344 // return nullptr;
3345 // }
3346 //
3347 // constant_* substract(constant_* c1, const constant_& c2){ /**< substracts c2 from c1, updates its type and returns the result **/
3348 // switch (c2.get_type()) {
3349 // case binary_c: {
3350 // return substract(c1, *(constant<bool>*)&c2);
3351 // break;
3352 // }
3353 // case short_c: {
3354 // return substract(c1, *(constant<short>*)&c2);
3355 // break;
3356 // }
3357 // case integer_c: {
3358 // return substract(c1, *(constant<int>*)&c2);
3359 // break;
3360 // }
3361 // case float_c: {
3362 // return substract(c1, *(constant<float>*)&c2);
3363 // break;
3364 // }
3365 // case double_c: {
3366 // return substract(c1, *(constant<double>*)&c2);
3367 // break;
3368 // }
3369 // case long_c: {
3370 // return substract(c1, *(constant<long double>*)&c2);
3371 // break;
3372 // }
3373 // case complex_c: {
3374 // return substract(c1, *(constant<Cpx>*)&c2);
3375 // break;
3376 // }
3377 // case par_c:{
3378 // auto pc2 = (param_*)(&c2);
3379 // switch (pc2->get_intype()) {
3380 // case binary_:
3381 // return substract(c1, *(param<bool>*)pc2);
3382 // break;
3383 // case short_:
3384 // return substract(c1, *(param<short>*)pc2);
3385 // break;
3386 // case integer_:
3387 // return substract(c1, *(param<int>*)pc2);
3388 // break;
3389 // case float_:
3390 // return substract(c1, *(param<float>*)pc2);
3391 // break;
3392 // case double_:
3393 // return substract(c1, *(param<double>*)pc2);
3394 // break;
3395 // case long_:
3396 // return substract(c1, *(param<long double>*)pc2);
3397 // break;
3398 // case complex_:
3399 // return substract(c1, *(param<Cpx>*)pc2);
3400 // break;
3401 // default:
3402 // break;
3403 // }
3404 // break;
3405 // }
3406 // case var_c:{
3407 // auto pc2 = (param_*)(&c2);
3408 // switch (pc2->get_intype()) {
3409 // case binary_:
3410 // return substract(c1, *(var<bool>*)pc2);
3411 // break;
3412 // case short_:
3413 // return substract(c1, *(var<short>*)pc2);
3414 // break;
3415 // case integer_:
3416 // return substract(c1, *(var<int>*)pc2);
3417 // break;
3418 // case float_:
3419 // return substract(c1, *(var<float>*)pc2);
3420 // break;
3421 // case double_:
3422 // return substract(c1, *(var<double>*)pc2);
3423 // break;
3424 // case long_:
3425 // return substract(c1, *(var<long double>*)pc2);
3426 // break;
3427 // case complex_:
3428 // return substract(c1, *(var<Cpx>*)pc2);
3429 // break;
3430 // default:
3431 // break;
3432 // }
3433 // break;
3434 // }
3435 // case uexp_c: {
3436 // auto res = new bexpr(minus_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3437 // delete c1;
3438 // c1 = (constant_*)res;
3439 // return c1;
3440 // break;
3441 // }
3442 // case bexp_c: {
3443 // auto res = new bexpr(minus_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3444 // delete c1;
3445 // c1 = (constant_*)res;
3446 // return c1;
3447 // break;
3448 // }
3449 // case func_c: {
3450 // auto f2 = (func_*)&c2;
3451 // if (f2->is_number()) {
3452 // return substract(c1, *f2->get_cst());
3453 // }
3454 // auto f = new func_(*c1);
3455 // delete c1;
3456 // *f -= *f2;
3457 // return c1 = (constant_*)f;
3458 // break;
3459 // break;
3460 // }
3461 // }
3462 // return nullptr;
3463 // }
3464 //
3465 // constant_* multiply(constant_* c1, const constant_& c2){ /**< adds c2 to c1, updates its type and returns the result **/
3466 // switch (c2.get_type()) {
3467 // case binary_c: {
3468 // return multiply(c1, *(constant<bool>*)&c2);
3469 // break;
3470 // }
3471 // case short_c: {
3472 // return multiply(c1, *(constant<short>*)&c2);
3473 // break;
3474 // }
3475 // case integer_c: {
3476 // return multiply(c1, *(constant<int>*)&c2);
3477 // break;
3478 // }
3479 // case float_c: {
3480 // return multiply(c1, *(constant<float>*)&c2);
3481 // break;
3482 // }
3483 // case double_c: {
3484 // return multiply(c1, *(constant<double>*)&c2);
3485 // break;
3486 // }
3487 // case long_c: {
3488 // return multiply(c1, *(constant<long double>*)&c2);
3489 // break;
3490 // }
3491 // case complex_c: {
3492 // return multiply(c1, *(constant<Cpx>*)&c2);
3493 // break;
3494 // }
3495 // case uexp_c: {
3496 // auto res = new bexpr(product_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3497 // delete c1;
3498 // c1 = (constant_*)res;
3499 // return c1;
3500 // break;
3501 // }
3502 // case bexp_c: {
3503 // auto res = new bexpr(product_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3504 // delete c1;
3505 // c1 = (constant_*)res;
3506 // return c1;
3507 // break;
3508 // }
3509 // case par_c:{
3510 // auto pc2 = (param_*)(&c2);
3511 // switch (pc2->get_intype()) {
3512 // case binary_:
3513 // return multiply(c1, *(param<bool>*)pc2);
3514 // break;
3515 // case short_:
3516 // return multiply(c1, *(param<short>*)pc2);
3517 // break;
3518 // case integer_:
3519 // return multiply(c1, *(param<int>*)pc2);
3520 // break;
3521 // case float_:
3522 // return multiply(c1, *(param<float>*)pc2);
3523 // break;
3524 // case double_:
3525 // return multiply(c1, *(param<double>*)pc2);
3526 // break;
3527 // case long_:
3528 // return multiply(c1, *(param<long double>*)pc2);
3529 // break;
3530 // case complex_:
3531 // return multiply(c1, *(param<Cpx>*)pc2);
3532 // break;
3533 // default:
3534 // break;
3535 // }
3536 // break;
3537 // }
3538 // case var_c:{
3539 // auto pc2 = (param_*)(&c2);
3540 // switch (pc2->get_intype()) {
3541 // case binary_:
3542 // return multiply(c1, *(var<bool>*)pc2);
3543 // break;
3544 // case short_:
3545 // return multiply(c1, *(var<short>*)pc2);
3546 // break;
3547 // case integer_:
3548 // return multiply(c1, *(var<int>*)pc2);
3549 // break;
3550 // case float_:
3551 // return multiply(c1, *(var<float>*)pc2);
3552 // break;
3553 // case double_:
3554 // return multiply(c1, *(var<double>*)pc2);
3555 // break;
3556 // case long_:
3557 // return multiply(c1, *(var<long double>*)pc2);
3558 // break;
3559 // case complex_:
3560 // return multiply(c1, *(var<Cpx>*)pc2);
3561 // break;
3562 // default:
3563 // break;
3564 // }
3565 // break;
3566 // }
3567 // case func_c: {
3568 // auto f = new func_(c2);
3569 // *f *= *c1;
3570 // delete c1;
3571 // return c1 = (constant_*)f;
3572 // break;
3573 // }
3574 // default:
3575 // break;
3576 // }
3577 // return nullptr;
3578 // }
3579 //
3580 //
3581 // constant_* divide(constant_* c1, const constant_& c2){
3582 // switch (c2.get_type()) {
3583 // case binary_c: {
3584 // return divide(c1, *(constant<bool>*)&c2);
3585 // break;
3586 // }
3587 // case short_c: {
3588 // return divide(c1, *(constant<short>*)&c2);
3589 // break;
3590 // }
3591 // case integer_c: {
3592 // return divide(c1, *(constant<int>*)&c2);
3593 // break;
3594 // }
3595 // case float_c: {
3596 // return divide(c1, *(constant<float>*)&c2);
3597 // break;
3598 // }
3599 // case double_c: {
3600 // return divide(c1, *(constant<double>*)&c2);
3601 // break;
3602 // }
3603 // case long_c: {
3604 // return divide(c1, *(constant<long double>*)&c2);
3605 // break;
3606 // }
3607 // case complex_c: {
3608 // return divide(c1, *(constant<Cpx>*)&c2);
3609 // break;
3610 // }
3611 // case uexp_c: {
3612 // auto res = new bexpr(product_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3613 // delete c1;
3614 // c1 = (constant_*)res;
3615 // return c1;
3616 // break;
3617 // }
3618 // case bexp_c: {
3619 // auto res = new bexpr(product_, make_shared<func_>((*c1)), make_shared<func_>((c2)));
3620 // delete c1;
3621 // c1 = (constant_*)res;
3622 // return c1;
3623 // break;
3624 // }
3625 // default:{
3626 // auto f = new func_(*c1);
3627 // *f /= func_(c2);
3628 // delete c1;
3629 // return c1 = (constant_*)f;
3630 // break;
3631 // }
3632 // }
3633 // return nullptr;
3634 // }
3635 //
3636 //
3637 // shared_ptr<func_> func_::get_stored_derivative(const string& vid) const{
3638 // auto it = _dfdx->find(vid);
3639 // if (it!=_dfdx->end()) {
3640 // return it->second;
3641 // }
3642 // else {
3643 // throw invalid_argument("No derivatives stored!\n");
3644 // }
3645 // }
3646 //
3647 // func_ func_::get_dfdx(const param_& v){
3648 // _dfdx->clear();
3649 // compute_derivatives();
3650 // return *get_stored_derivative(v._name).get();
3651 // }
3652 //
3653 // func_* func_::compute_derivative(const param_ &v){
3654 // auto vid = v._name;
3655 // if(_dfdx->count(vid)!=0){
3656 // return _dfdx->at(vid).get();
3657 // }
3658 //
3659 // auto df = new func_(get_derivative(v));
3660 // df->allocate_mem();
3661 // (*_dfdx)[vid] = shared_ptr<func_>(df);
3662 // DebugOff( "First derivative with respect to " << v.get_name(false,false) << " = " << df->to_str() << endl);
3663 // return df;
3664 // }
3665 //
3666 // void func_::compute_derivatives(){ /**< Computes and stores the derivative of f with respect to all variables. */
3667 // size_t vid = 0, vjd = 0;
3668 // param_* vi;
3669 // param_* vj;
3670 // DebugOff( "Computing derivatives for " << to_str() << endl);
3671 // for (auto &vp: *_vars) {
3672 // vi = vp.second.first.get();
3673 // vid = vi->get_id();
3674 // auto vi_name = vp.first;
3675 // auto df = compute_derivative(*vi);
3676 // // if (is_nonlinear()) {
3677 // DebugOff( "First derivative with respect to " << vp.first << " = " << df->to_str() << endl);
3678 // // }
3679 // for (auto &vp2: *df->_vars) {
3680 // vj = vp2.second.first.get();
3681 // vjd = vj->get_id();
3682 // auto vj_name = vp2.first;
3683 // if (vi_name.compare(vj_name) <= 0) { //only store lower left part of hessian matrix since it is symmetric.
3684 // auto d2f = df->compute_derivative(*vj);
3685 // DebugOff( "Second derivative with respect to " << vp2.first << " and " << vp.first << " = " << d2f->to_str() << endl);
3686 // // d2f->print();
3687 // }
3688 // }
3689 //
3690 // }
3691 // };
3692 //
3693 /* Build Compositions of positive integer n given k positive integers. All credits go to Dr. Martin von Gagern
3694 @param[in] k: the sum
3695 @param[in] n: the number of integers
3696 */
build_compositions(int k,int n)3697 vector<vector<int>> build_compositions(int k, int n) {
3698 int pos = 0, last = n - 1;
3699 vector<vector<int>> res;
3700 vector<int> a;
3701 a.resize(n);
3702 a[0] = k;
3703 res.push_back(a);
3704 while (true) {
3705 if (pos != last) {
3706 --a[pos];
3707 ++pos;
3708 a[pos] = 1;
3709 }
3710 else {
3711 if (a[last] == k)
3712 return res;
3713 for (--pos; a[pos] == 0; --pos);
3714 --a[pos];
3715 int tmp = 1 + a[last];
3716 ++pos;
3717 a[last] = 0;
3718 a[pos] = tmp;
3719 }
3720 res.push_back(a);
3721 }
3722 }
3723
3724 /* Build all monomials of degree d in dimension n (size of vars) */
get_monomials(unsigned d)3725 vector<pterm> func_::get_monomials(unsigned d){
3726 vector<pterm> res;
3727 int n = _vars->size();
3728 auto comp = build_compositions(d, n);
3729 auto m = comp.size();
3730 res.resize(m);
3731 for (auto i = 0; i<m; i++) {
3732 auto l = make_shared<list<pair<shared_ptr<param_>, int>>>();
3733 auto row = comp[i];
3734 auto v_it = _vars->begin();
3735 for (auto j = 0; j<n; j++) {
3736 auto exp = row[j];
3737 auto v = v_it->second.first;
3738 if (exp>0) {
3739 l->push_back(make_pair<>(v, exp));
3740 }
3741 v_it++;
3742 }
3743 res[i] = pterm(true, make_shared<constant<>>(1), l);
3744 cout << res[i].to_str(0, 1) << " | ";
3745 }
3746 cout << endl;
3747 return res;
3748 }
3749 //
3750 //
3751 //
3752 // void func_::reset_val(){
3753 // _evaluated = false;
3754 // if (_expr) {
3755 // _expr->reset_val();
3756 // }
3757 // }
3758 //
3759 // void func_::reindex(param_* v){
3760 //
3761 // }
3762 //
3763 // void func_::replace(param_* v, func_& f){
3764 // auto vin = get_var(*v->_vec_id);
3765 // f.reindex(v);
3766 // //index f here!
3767 // for (auto &l_p:*_lterms) {
3768 // auto < = l_p.second;
3769 // if (lt._p==vin) {
3770 // //replace here.
3771 // }
3772 // }
3773 //
3774 // if (_expr) {
3775 // if (_expr->is_uexpr()) {
3776 // auto ue = (uexpr*)_expr.get();
3777 // ue->_son->replace(v,f);
3778 // }
3779 // else {
3780 // auto be = (bexpr*)_expr.get();
3781 // be->_lson->replace(v,f);
3782 // be->_rson->replace(v,f);
3783 // }
3784 // }
3785 // _vars->erase(v->_name);
3786 // DebugOn(to_str());
3787 // // print();
3788 // }
3789 //
get_var(size_t vid) const3790 shared_ptr<param_> func_::get_var(size_t vid) const{
3791 for (auto &v_p:*_vars) {
3792 if (*v_p.second.first->_vec_id == vid) {
3793 return v_p.second.first;
3794 }
3795 }
3796 return nullptr;
3797 }
3798
3799
has_sym_var(const param_ & v) const3800 bool func_::has_sym_var(const param_& v) const{
3801 return get_var(*v._vec_id)!=nullptr;
3802 }
3803 //
has_var(const param_ & v) const3804 bool func_::has_var(const param_& v) const{
3805 return get_var(v.get_name(true, false))!=nullptr;
3806 }
3807
3808 //
3809 //
3810 // bool func_::has_var(const string& name) const{
3811 // return _vars->count(name)>0;
3812 // }
3813 //
3814 // func_ func_::get_derivative(const param_ &v) const{
3815 // func_ res;
3816 // if(!has_var(v)){
3817 // return res;
3818 // }
3819 // for (auto <: *_lterms) {
3820 // if (*lt.second._p == v) {
3821 // auto coef = copy(*lt.second._coef);
3822 // if ((coef->_is_vector && coef->_is_transposed)) {
3823 // coef->transpose();
3824 // }
3825 // res = move(*coef);
3826 // delete coef;
3827 // if(!lt.second._sign){
3828 // res *= -1;
3829 // }
3830 // }
3831 // }
3832 // for (auto <: *_qterms) {
3833 // if (*lt.second._p->first == v) {
3834 // auto coef = copy(*lt.second._coef);
3835 // if (coef->_is_vector && coef->_is_transposed) {
3836 // coef->transpose();
3837 // }
3838 // if(lt.second._sign) {
3839 //// if(lt.second._c_p1_transposed){
3840 //// res += (*coef*(*lt.second._p->second)).tr();
3841 //// }
3842 //// else {
3843 // res += *coef*(*lt.second._p->second);
3844 //// }
3845 // }
3846 // else {
3847 //// if(lt.second._c_p1_transposed){
3848 //// res -= (*coef*(*lt.second._p->second)).tr();
3849 //// }
3850 //// else {
3851 // res -= *coef*(*lt.second._p->second);
3852 //// }
3853 // }
3854 // delete coef;
3855 // if (lt.second._coef->_is_vector || lt.second._coef->is_matrix()) {
3856 // res._is_vector = true;
3857 // // res._nb_instances = lt.second._coef->get_dim();
3858 // }
3859 // }
3860 // if (*lt.second._p->second == v) {
3861 // auto coef = copy(*lt.second._coef);
3862 // if (coef->_is_vector && coef->_is_transposed) {
3863 // coef->transpose();
3864 // }
3865 // if(lt.second._sign) {
3866 // res += *coef*(*lt.second._p->first);
3867 // }
3868 // else {
3869 // res -= *coef*(*lt.second._p->first);
3870 // }
3871 // delete coef;
3872 // if (lt.second._coef->_is_vector || lt.second._coef->is_matrix()) {
3873 // res._is_vector = true;
3874 // }
3875 // }
3876 //
3877 // //CHECK THIS
3878 // // if (lt.second._coef->_is_vector || lt.second._coef->_is_matrix) {
3879 // // res._is_vector = true;
3880 // // res._nb_instances = lt.second._coef->get_dim();
3881 // // res._dim[0] = lt.second._coef->_dim[0]; res._dim[1] = lt.second._coef->_dim[1];
3882 // // }
3883 // }
3884 // for (auto <: *_pterms) {
3885 // for (auto &p: *lt.second._l) {
3886 // if (*p.first == v) {
3887 // func_ pterm = constant<double>(1);
3888 // if (!lt.second._sign) {
3889 // pterm = constant<double>(-1);
3890 // }
3891 // auto expo = p.second;
3892 // if (expo > 1) {
3893 // pterm *= expo;
3894 // pterm *= *lt.second._coef;
3895 // pterm *= (*p.first);
3896 // for (int i = 1; i<expo-1; i++) {
3897 // pterm *= *p.first;
3898 // }
3899 // }
3900 // for (auto &p2: *lt.second._l) {
3901 // if (p2!=p) {
3902 // func_ pterm2(*p2.first);
3903 // for (int i = 1; i<p2.second; i++) {
3904 // pterm2 *= *p2.first;
3905 // }
3906 // pterm *= pterm2;
3907 // }
3908 // }
3909 // res += *lt.second._coef*pterm;
3910 // }
3911 // }
3912 // }
3913 // if (!_expr) {
3914 // // res.untranspose();
3915 // return res;
3916 // }
3917 // else { // f is a composition of functions
3918 // res += _expr->get_derivative(v);
3919 // // res.untranspose();
3920 // return res;
3921 // }
3922 // }
3923 //
3924 // double func_::get_val(size_t i) const{
3925 // // if (is_constant()) {
3926 // // if(eval(i)!=_val->at(i)){
3927 // // throw invalid_argument("error");
3928 // // }
3929 // // }
3930 // // if (is_number()) {
3931 // // return _val->at(0);
3932 // // }
3933 // // else {
3934 // // if (i>=_val->size()) {
3935 // // throw invalid_argument("error");
3936 // // }
3937 // if (_val->size()<=i){
3938 // throw invalid_argument("Func get_val out of range");
3939 // }
3940 // return _val->at(i);
3941 // // }
3942 // }
3943 //
3944 // double func_::get_val(size_t i, size_t j) const{
3945 //
3946 // // if () {
3947 // // return _val->at(0);
3948 // // }
3949 // // if (!_is_matrix) {
3950 // // // if (_is_transposed) {
3951 // // // return get_val(i);
3952 // // // }
3953 // // throw invalid_argument("get_val i,j");
3954 // // return get_val(j);
3955 // // }
3956 // // if (is_number()) {
3957 // // return _val->at(0);
3958 // // }
3959 // // else {
3960 // if (_is_transposed) {
3961 // // if (j*_dim[0]+i>=_val->size()) {
3962 // // throw invalid_argument("error");
3963 // // }
3964 // return _val->at(j*_dim[0]+i);
3965 // }
3966 // // if (i*_dim[1]+j>=_val->size()) {
3967 // // throw invalid_argument("error");
3968 // // }
3969 // return _val->at(i*_dim[1]+j);
3970 // // }
3971 // }
3972 //
3973 // // double func_::force_eval(size_t i){
3974 // //
3975 // // double res = 0;
3976 // // for (auto &pair:*_pterms) {
3977 // // res += pair.second.eval(i);
3978 // // }
3979 // // for (auto &pair:*_qterms) {
3980 // // res += pair.second.eval(i);
3981 // // }
3982 // // for (auto &pair:*_lterms) {
3983 // // res += pair.second.eval(i);
3984 // // }
3985 // // res += eval(_cst,i);
3986 // // if (_expr) {
3987 // // res += _expr->eval(i);
3988 // // }
3989 // // return res;
3990 // // }
3991 //
3992 // double func_::eval(size_t i, size_t j){
3993 // // if (_val->size()<_nb_instances) {
3994 // // _val->resize(_nb_instances);
3995 // // }
3996 // if (!is_matrix() && (!_indices || !_indices->_ids || _indices->_ids->size()==1)) {
3997 // // if (!is_matrix()) {
3998 // // if (_is_transposed) {
3999 // // return eval(i);
4000 // // }
4001 // return eval(j);
4002 // }
4003 //
4004 // if (is_constant() && _evaluated) {
4005 // if (is_number()) {
4006 // return _val->at(0);
4007 // }
4008 // else {
4009 // if (_is_transposed) {
4010 // return _val->at(j*_dim[0]+i);
4011 // }
4012 // return _val->at(i*_dim[1]+j);
4013 // }
4014 // }
4015 // double res = 0;
4016 // for (auto &pair:*_pterms) {
4017 // res += pair.second.eval(i,j);
4018 // }
4019 // for (auto &pair:*_qterms) {
4020 // res += pair.second.eval(i,j);
4021 // }
4022 // for (auto &pair:*_lterms) {
4023 // res += pair.second.eval(i,j);
4024 // }
4025 // res += t_eval(_cst,i,j);
4026 // if (_expr) {
4027 // res += _expr->eval(i,j);
4028 // }
4029 // if (is_number()) {
4030 // _val->at(0) = res;
4031 // }
4032 // else {
4033 // if (is_constant()) {
4034 // if (_is_transposed) {
4035 // if(_dim[0]*j + i ==_val->size()-1){
4036 // _evaluated = true;
4037 // }
4038 //
4039 // }
4040 // else {
4041 // if(_dim[1]*i + j ==_val->size()-1){
4042 // _evaluated = true;
4043 // }
4044 // }
4045 //
4046 // }
4047 // set_val(i,j,res);
4048 // }
4049 // return res;
4050 // }
4051 //
4052 // double func_::eval(size_t i){
4053 // if (_val->size()<=i) {
4054 // allocate_mem();
4055 // }
4056 // // if (!_val) {
4057 // // throw invalid_argument("_val not defined for function.\n");
4058 // // }
4059 // if (is_constant() && _evaluated) {
4060 // if (is_number()) {
4061 // return _val->at(0);
4062 // }
4063 // // if (i>=_val->size()) {
4064 // // throw invalid_argument("error");
4065 // // }
4066 // // if (_val->at(i) != force_eval(i)) {
4067 // // throw invalid_argument("error");
4068 // // }
4069 // if (_val->size()<=i){
4070 // throw invalid_argument("Func eval out of range");
4071 // }
4072 // return _val->at(i);
4073 // }
4074 // double res = 0;
4075 // for (auto &pair:*_pterms) {
4076 // res += pair.second.eval(i);
4077 // }
4078 // for (auto &pair:*_qterms) {
4079 // res += pair.second.eval(i);
4080 // }
4081 // for (auto &pair:*_lterms) {
4082 // res += pair.second.eval(i);
4083 // }
4084 // res += t_eval(_cst,i);
4085 // if (_expr) {
4086 // // if (_expr->is_uexpr()) {
4087 // // auto ue = (uexpr*)_expr.get();
4088 // //// if (ue->_son->is_constant()) {
4089 // //// _nb_instances = max(_nb_instances, ue->_son->_nb_instances);
4090 // //// _val->resize(max(_val->size(),ue->_son->_nb_instances));//TODO is this necessary?
4091 // //// }
4092 // //
4093 // // }
4094 // // else {
4095 // // auto be = (bexpr*)_expr.get();
4096 // // if(!be->is_inner_product()) {
4097 // // _nb_instances = max(_nb_instances, max(be->_lson->_nb_instances,be->_rson->_nb_instances));
4098 // // }
4099 // //// if (be->_lson->is_constant() && be->_rson->is_constant()) {
4100 // //// _val->resize(max(_val->size(),max(be->_lson->_nb_instances,be->_rson->_nb_instances)));
4101 // //// }
4102 // //
4103 // // }
4104 // res += _expr->eval(i);
4105 // }
4106 // if (is_number()) {
4107 // _val->at(0) = res;
4108 // _evaluated = true;//TODO fix this
4109 // }
4110 // else {
4111 // // if (i>=_val->size()) {
4112 // // throw invalid_argument("error");
4113 // // }
4114 // if (is_constant() && i==get_dim()-1) {
4115 // _evaluated = true;
4116 // }
4117 // if (_val->size()<=i){
4118 // throw invalid_argument("Param eval out of range");
4119 // }
4120 // _val->at(i) = res;
4121 // }
4122 // return res;
4123 // }
4124 //
4125 // string func_::to_str() const{ return string();};
4126 // string str;
4127 // int ind = 0;
4128 // string sign = " + ";
4129 // // for (int inst = 0; inst < _nb_instances ; inst++) {
4130 // ind = 0;
4131 //
4132 // for (auto &pair:*_pterms) {
4133 // str += pair.second.to_str(ind++);
4134 // }
4135 // if (!_pterms->empty() && (!_qterms->empty() || !_lterms->empty())) {
4136 // str += " + ";
4137 // }
4138 // ind = 0;
4139 // for (auto &pair:*_qterms) {
4140 // str += pair.second.to_str(ind++);
4141 // }
4142 // if (!_qterms->empty() && !_lterms->empty()) {
4143 // str += " + ";
4144 // }
4145 // ind = 0;
4146 // for (auto &pair:*_lterms) {
4147 // str += pair.second.to_str(ind++);
4148 // }
4149 // if (_cst->is_number()) {
4150 // auto val = poly_to_str(_cst);
4151 // if (val.front()=='-') {
4152 // str += " - " + val.substr(1);
4153 // }
4154 // else if (val != "0"){
4155 // if (!_pterms->empty() || !_qterms->empty() || !_lterms->empty()) {
4156 // str += " + ";
4157 // }
4158 // str += val;
4159 // }
4160 // }
4161 // else {
4162 // if (!_pterms->empty() || !_qterms->empty() || !_lterms->empty()) {
4163 // str += " + ";
4164 // }
4165 // str += "(";
4166 // str += poly_to_str(_cst);
4167 // str += ")";
4168 // }
4169 // if (_expr && (!_pterms->empty() || !_qterms->empty() || !_lterms->empty() || !_cst->is_zero())) {
4170 // str += " + ";
4171 // }
4172 // if (_expr) {
4173 // str += _expr->get_str();
4174 // }
4175 // if (_is_vector && (is_number() || _vars->size()>1 || _params->size()>1)) {
4176 // str = "[" + str +"]";
4177 // }
4178 // if (_is_transposed) {
4179 // str += "\u1D40";
4180 // }
4181 // if(str.size()==0){
4182 // str = "0";
4183 // }
4184 // return str;
4185 // }
4186 //
4187 //
4188
4189 /** Return the number of linear terms in this function */
nb_linear_terms() const4190 unsigned func_::nb_linear_terms() const{
4191 return _lterms->size();
4192 }
4193
4194
4195
4196
4197 template<typename type>
get_scale_factor(double unit)4198 double func<type>::get_scale_factor(double unit){
4199 double absmax = 0;
4200 for (auto <:this->get_lterms()) {
4201 if (lt.second._coef->is_function()) {
4202 auto coef = static_pointer_cast<func<type>>(lt.second._coef);
4203 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4204 }
4205 else if(lt.second._coef->is_param()) {
4206 auto coef = static_pointer_cast<param<type>>(lt.second._coef);
4207 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4208 }
4209 else if(lt.second._coef->is_number()) {
4210 auto coef = static_pointer_cast<constant<type>>(lt.second._coef);
4211 absmax = std::max(absmax,std::abs(coef->eval()));
4212 }
4213 }
4214 for (auto &qt:this->get_qterms()) {
4215 if (qt.second._coef->is_function()) {
4216 auto coef = static_pointer_cast<func<type>>(qt.second._coef);
4217 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4218 }
4219 else if(qt.second._coef->is_param()) {
4220 auto coef = static_pointer_cast<param<type>>(qt.second._coef);
4221 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4222 }
4223 else if(qt.second._coef->is_number()) {
4224 auto coef = static_pointer_cast<constant<type>>(qt.second._coef);
4225 absmax = std::max(absmax,std::abs(coef->eval()));
4226 }
4227 }
4228 for (auto &pt:this->get_pterms()) {
4229 if (pt.second._coef->is_function()) {
4230 auto coef = static_pointer_cast<func<type>>(pt.second._coef);
4231 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4232 }
4233 else if(pt.second._coef->is_param()) {
4234 auto coef = static_pointer_cast<param<type>>(pt.second._coef);
4235 absmax = std::max(absmax,std::max(std::abs(coef->_range->first),std::abs(coef->_range->second)));
4236 }
4237 else if(pt.second._coef->is_number()) {
4238 auto coef = static_pointer_cast<constant<type>>(pt.second._coef);
4239 absmax = std::max(absmax,std::abs(coef->eval()));
4240 }
4241 }
4242 if(_expr){
4243 /* TODO: NL part*/
4244 }
4245 if(absmax>unit)
4246 return unit/absmax;
4247 return 1;
4248 }
4249
4250 /* Make sure all coefficients are in [-unit,unit] */
4251 template<typename type>
scale_coefs(double unit)4252 void func<type>::scale_coefs(double unit){
4253 double factor = get_scale_factor(unit);
4254 *this *= factor;
4255 }
4256
4257 template<typename type>
4258 template<typename T>
replace(const var<T> & v,const func<T> & f) const4259 func<type> func<type>::replace(const var<T>& v, const func<T>& f) const{/**< Replace v with function f everywhere it appears. Assumes v has only unique keys! */
4260 if(!has_sym_var(v)) return *this;
4261 func<type> new_this = *this;
4262 vector<bool> keep_ids, diff_ids;
4263 bool updated_ids = false;
4264 for (auto <:this->get_lterms()) {
4265 auto vv = lt.second._p;
4266 if(v.get_vec_id()==vv->get_vec_id()){
4267 auto f_cpy = f;
4268 if(vv->_indices->is_subset(*v._indices)){/* Case where vv index set is a subset of v, we can safely remove the linear term including vv */
4269 if(vv->is_matrix_indexed() && !v.is_matrix_indexed()){
4270 f_cpy.update_var_indices(*vv->_indices, *v._indices);
4271 }
4272 else {
4273 keep_ids = v._indices->get_common_refs(*vv->_indices); /* indices to keep */
4274 if(!updated_ids && f_cpy._indices && !v.is_matrix_indexed() && vv->_indices->size() != f_cpy._indices->size()){
4275 f_cpy.update_rows(keep_ids);
4276 f_cpy.keep_unique_keys();
4277 updated_ids = true;
4278 }
4279 }
4280 auto vcast = *static_pointer_cast<var<T>>(vv);
4281 new_this.insert(!lt.second._sign,*lt.second._coef,*lt.second._p);/* Remove coef*vv from new_f */
4282 if (lt.second._coef->is_function()) {
4283 auto coef = *static_pointer_cast<func<T>>(lt.second._coef);
4284 if(vv->is_matrix_indexed() && coef._is_transposed){
4285 new_this += std::pow(-1,1-lt.second._sign)*coef.tr()*f_cpy;
4286 }
4287 else {
4288 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4289 }
4290 }
4291 else if(lt.second._coef->is_param()) {
4292 auto coef = *static_pointer_cast<param<T>>(lt.second._coef);
4293 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4294 }
4295 else if(lt.second._coef->is_number()) {
4296 auto coef = *static_pointer_cast<constant<T>>(lt.second._coef);
4297 if(vv->is_matrix_indexed() && coef._is_transposed){
4298 new_this += std::pow(-1,1-lt.second._sign)*coef.tr()*f_cpy;
4299 }
4300 else {
4301 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4302 }
4303 }
4304 }
4305 else if(vv->_indices->is_superset(*v._indices)){ /* vv has more indices */
4306 keep_ids = vv->_indices->get_common_refs(*v._indices); /* indices to keep */
4307 diff_ids = vv->_indices->get_diff_refs(*v._indices); /* indices not in v */
4308 auto diff_set = indices(*vv->_indices);
4309 diff_set.filter_refs(diff_ids);
4310 if(!updated_ids && new_this._indices && !vv->is_matrix_indexed()){
4311 new_this.update_rows(keep_ids);
4312 new_this.keep_unique_keys();
4313 updated_ids = true;
4314 }
4315 auto vcast = *static_pointer_cast<var<T>>(vv);
4316 new_this.insert(!lt.second._sign,*lt.second._coef,*lt.second._p);/* Remove coef*v from new_f */
4317 if (lt.second._coef->is_function()) {
4318 auto coef = *static_pointer_cast<func<T>>(lt.second._coef);
4319 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4320 if(vcast._is_vector){
4321 new_this += coef*vcast.in(diff_set);
4322 }
4323 }
4324 else if(lt.second._coef->is_param()) {
4325 auto coef = *static_pointer_cast<param<T>>(lt.second._coef);
4326 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4327 if(vcast._is_vector){
4328 new_this += coef*vcast.in(diff_set);
4329 }
4330 }
4331 else if(lt.second._coef->is_number()) {
4332 auto coef = *static_pointer_cast<constant<T>>(lt.second._coef);
4333 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4334 if(vcast._is_vector){
4335 new_this += coef*vcast.in(diff_set);
4336 }
4337 }
4338 }
4339 else {
4340 if(vv->is_matrix_indexed() && !v.is_matrix_indexed()){
4341 f_cpy.update_var_indices(*vv->_indices, *v._indices);
4342 }
4343 new_this.insert(!lt.second._sign,*lt.second._coef,*lt.second._p);/* Remove coef*vv from new_f */
4344 if (lt.second._coef->is_function()) {
4345 auto coef = *static_pointer_cast<func<T>>(lt.second._coef);
4346 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4347 }
4348 else if(lt.second._coef->is_param()) {
4349 auto coef = *static_pointer_cast<param<T>>(lt.second._coef);
4350 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4351 }
4352 else if(lt.second._coef->is_number()) {
4353 auto coef = *static_pointer_cast<constant<T>>(lt.second._coef);
4354 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy;
4355 }
4356 }
4357 }
4358 }
4359 for (auto <:this->get_qterms()) {
4360 auto vv1 = lt.second._p->first;
4361 auto vv2 = lt.second._p->second;
4362 if(v.get_vec_id()==vv1->get_vec_id() || v.get_vec_id()==vv2->get_vec_id()){ /* Quadratic term with v in it*/
4363 auto f_cpy = f;
4364 auto v_cpy = v;
4365 shared_ptr<param_> vv = nullptr;
4366 if(v.get_vec_id()==vv1->get_vec_id()){
4367 vv = vv1;
4368 }
4369 else{
4370 vv = vv2;
4371 }
4372 if(vv->_indices->is_subset(*v._indices)){/* Case where vv index set is a subset of v, wecan safely remove the lienear term including vv */
4373 keep_ids = v._indices->get_common_refs(*vv->_indices); /* indices to keep */
4374 if(!updated_ids && f_cpy._indices && !v.is_matrix_indexed()){
4375 f_cpy.update_rows(keep_ids);
4376 f_cpy.keep_unique_keys();
4377 updated_ids = true;
4378 }
4379 auto vcast = *static_pointer_cast<var<T>>(vv);
4380 new_this.insert(!lt.second._sign,*lt.second._coef,*vv1,*vv2);/* Remove coef*v*v from new_f */
4381 if (lt.second._coef->is_function()) {
4382 auto coef = *static_pointer_cast<func<T>>(lt.second._coef);
4383 if(vv->is_matrix_indexed() && coef._is_transposed){//?
4384 coef.transpose();
4385 }
4386 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4387 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4388 }
4389 else if(v.get_vec_id()==vv1->get_vec_id()){
4390 auto v2 = *static_pointer_cast<var<T>>(vv2);
4391 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4392 }
4393 else {
4394 auto v1 = *static_pointer_cast<var<T>>(vv1);
4395 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4396 }
4397 }
4398 else if(lt.second._coef->is_param()) {
4399 auto coef = *static_pointer_cast<param<T>>(lt.second._coef);
4400 if(vv->is_matrix_indexed() && coef._is_transposed){//?
4401 coef.transpose();
4402 }
4403 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4404 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4405 }
4406 else if(v.get_vec_id()==vv1->get_vec_id()){
4407 auto v2 = *static_pointer_cast<var<T>>(vv2);
4408 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4409 }
4410 else {
4411 auto v1 = *static_pointer_cast<var<T>>(vv1);
4412 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4413 }
4414 }
4415 else if(lt.second._coef->is_number()) {
4416 auto coef = *static_pointer_cast<constant<T>>(lt.second._coef);
4417 if(vv->is_matrix_indexed() && coef._is_transposed){//?
4418 coef.transpose();
4419 }
4420 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4421 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4422 }
4423 else if(v.get_vec_id()==vv1->get_vec_id()){
4424 auto v2 = *static_pointer_cast<var<T>>(vv2);
4425 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4426 }
4427 else {
4428 auto v1 = *static_pointer_cast<var<T>>(vv1);
4429 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4430 }
4431 }
4432 }
4433 else if(v._indices->is_subset(*vv->_indices)){ /* vv has more indices */
4434 if(f_cpy.has_matrix_indexed_vars()){/* cannot do this yet */
4435 return *this;
4436 }
4437 keep_ids = vv->_indices->get_common_refs(*v._indices); /* indices to keep */
4438 diff_ids = vv->_indices->get_diff_refs(*v._indices); /* indices not in v */
4439 auto diff_set = indices(*vv->_indices);
4440 diff_set.filter_refs(diff_ids);
4441 if(!updated_ids && new_this._indices){
4442 if(!vv->is_matrix_indexed()){
4443 new_this.update_rows(keep_ids);
4444 new_this.keep_unique_keys();
4445 updated_ids = true;
4446 }
4447 else {
4448 new_this.update_indices(keep_ids);
4449 new_this.keep_unique_keys();
4450 new_this.print();
4451 updated_ids = true;
4452 }
4453 }
4454 if(vv->is_matrix_indexed() && f_cpy.has_matrix_indexed_vars()){/* cannot do this yet */
4455 return *this;
4456 // f_cpy.update_var_indices(*new_this.get_var(vv->get_name(false,false))->_indices, *v._indices);
4457 // f_cpy.print();
4458 }
4459 auto v1 = *static_pointer_cast<var<T>>(new_this.get_var(vv1->get_name(false,false)));
4460 auto v2 = *static_pointer_cast<var<T>>(new_this.get_var(vv2->get_name(false,false)));
4461 new_this.insert(!lt.second._sign,*lt.second._coef,*vv1,*vv2);/* Remove coef*v*v from new_f */
4462 if (lt.second._coef->is_function()) {
4463 auto coef = *static_pointer_cast<func<T>>(lt.second._coef);
4464 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4465 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4466 }
4467 else if(v.get_vec_id()==vv1->get_vec_id()){
4468 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4469 }
4470 else {
4471 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4472 }
4473 if(v1._is_vector || v2._is_vector){
4474 new_this += coef*v1.in(diff_set)*v2.in(diff_set);
4475 }
4476 }
4477 else if(lt.second._coef->is_param()) {
4478 auto coef = *static_pointer_cast<param<T>>(lt.second._coef);
4479 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4480 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4481 }
4482 else if(v.get_vec_id()==vv1->get_vec_id()){
4483 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4484 }
4485 else {
4486 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4487 }
4488 if(v1._is_vector || v2._is_vector){
4489 new_this += coef*v1.in(diff_set)*v2.in(diff_set);
4490 }
4491 }
4492 else if(lt.second._coef->is_number()) {
4493 auto coef = *static_pointer_cast<constant<T>>(lt.second._coef);
4494 if(v.get_vec_id()==vv1->get_vec_id() && v.get_vec_id()==vv2->get_vec_id()){
4495 new_this += std::pow(-1,1-lt.second._sign)*coef*pow(f_cpy,2);
4496 }
4497 else if(v.get_vec_id()==vv1->get_vec_id()){
4498 new_this += std::pow(-1,1-lt.second._sign)*coef*f_cpy*v2;
4499 }
4500 else {
4501 new_this += std::pow(-1,1-lt.second._sign)*coef*v1*f_cpy;
4502 }
4503 if(v1._is_vector || v2._is_vector){
4504 new_this += coef*v1.in(diff_set)*v2.in(diff_set);
4505 }
4506 }
4507 }
4508 }
4509 }
4510 for (auto &term:this->get_pterms()) {
4511 for (auto &v_p:*term.second._l) {
4512 auto vv = v_p.first;
4513 if(v.get_vec_id()==vv->get_vec_id()){
4514 auto f_cpy = f;
4515 /* Get polynomial term excluding v */
4516 auto new_term = func<type>(term.second.exclude(vv));
4517 if(vv->_indices->is_subset(*v._indices)){/* Case where vv index set is a subset of v, we can safely remove the linear term including vv */
4518 if(vv->is_matrix_indexed() && !v.is_matrix_indexed()){
4519 f_cpy.update_var_indices(*vv->_indices, *v._indices);
4520 }
4521 else {
4522 keep_ids = v._indices->get_common_refs(*vv->_indices); /* indices to keep */
4523 if(!updated_ids && f_cpy._indices && !v.is_matrix_indexed() && vv->_indices->size() != f_cpy._indices->size()){
4524 f_cpy.update_rows(keep_ids);
4525 f_cpy.keep_unique_keys();
4526 updated_ids = true;
4527 }
4528 }
4529 auto vcast = *static_pointer_cast<var<T>>(vv);
4530 new_this.insert(!term.second._sign,*term.second._coef,*term.second._l);/* Remove polynomial term from new_f */
4531 if (term.second._coef->is_function()) {
4532 auto coef = *static_pointer_cast<func<T>>(term.second._coef);
4533 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4534 }
4535 else if(term.second._coef->is_param()) {
4536 auto coef = *static_pointer_cast<param<T>>(term.second._coef);
4537 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4538 }
4539 else if(term.second._coef->is_number()) {
4540 auto coef = *static_pointer_cast<constant<T>>(term.second._coef);
4541 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4542 }
4543 }
4544 else if(vv->_indices->is_superset(*v._indices)){ /* vv has more indices */
4545 keep_ids = vv->_indices->get_common_refs(*v._indices); /* indices to keep */
4546 diff_ids = vv->_indices->get_diff_refs(*v._indices); /* indices not in v */
4547 auto diff_set = indices(*vv->_indices);
4548 diff_set.filter_refs(diff_ids);
4549 if(!updated_ids && new_this._indices && !vv->is_matrix_indexed()){
4550 new_this.update_rows(keep_ids);
4551 new_this.keep_unique_keys();
4552 updated_ids = true;
4553 }
4554 auto vcast = *static_pointer_cast<var<T>>(vv);
4555 new_this.insert(!term.second._sign,*term.second._coef,*term.second._l);/* Remove coef*v from new_f */
4556 if (term.second._coef->is_function()) {
4557 auto coef = *static_pointer_cast<func<T>>(term.second._coef);
4558 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4559 if(vcast._is_vector){
4560 auto new_term = func<type>(term.second);
4561 new_term.update_var_indices(diff_set);
4562 new_this += coef*new_term.in(diff_set);
4563 }
4564 }
4565 else if(term.second._coef->is_param()) {
4566 auto coef = *static_pointer_cast<param<T>>(term.second._coef);
4567 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4568 if(vcast._is_vector){
4569 auto new_term = func<type>(term.second);
4570 new_term.update_var_indices(diff_set);
4571 new_this += coef*new_term.in(diff_set);
4572 }
4573 }
4574 else if(term.second._coef->is_number()) {
4575 auto coef = *static_pointer_cast<constant<T>>(term.second._coef);
4576 new_this += std::pow(-1,1-term.second._sign)*coef*new_term*pow(f_cpy,v_p.second);
4577 if(vcast._is_vector){
4578 auto new_term = func<type>(term.second);
4579 new_term.update_var_indices(diff_set);
4580 new_this += coef*new_term.in(diff_set);
4581 }
4582 }
4583 }
4584 }
4585 }
4586 }
4587
4588 /* TODO: nonlinear parts */
4589 // if (this->_expr) {
4590 // if (this->_expr->is_uexpr()) {
4591 // auto ue = static_pointer_cast<uexpr<type>>(this->_expr);
4592 // }
4593 // else {
4594 // auto be = static_pointer_cast<bexpr<type>>(this->_expr);
4595 // }
4596 // }
4597 return new_this;
4598 }
4599
print_symbolic(bool endline,bool display_input)4600 void func_::print_symbolic(bool endline, bool display_input){
4601 string str;
4602 if (display_input) {
4603 if (is_constant()) {
4604 cout << " (Constant";
4605 }
4606 else if (is_linear()) {
4607 cout << " (Linear";
4608 }
4609 else if (is_convex()) {
4610 cout << " (Convex";
4611 }
4612 else if (is_concave()){
4613 cout << " (Concave";
4614 }
4615 else {
4616 cout << " (Unknown";
4617 }
4618 if (is_complex()) {
4619 cout << " Complex) : ";
4620 }
4621 else {
4622 cout << ") : ";
4623 }
4624 if (!_embedded && !is_constant()) {
4625 str += "f(";
4626 for (auto pair_it = _vars->begin(); pair_it != _vars->end();) {
4627 str += pair_it->second.first->get_name(false,false);
4628 if (next(pair_it) != _vars->end()) {
4629 str += ",";
4630 }
4631 pair_it++;
4632 }
4633 str += ") = ";
4634 }
4635 }
4636 str += to_str();
4637 _to_str = str;
4638 cout << this->_to_str;
4639 if (endline)
4640 cout << endl;
4641 }
4642
4643
4644
4645
4646
4647 //
4648
4649 //
4650 // void func_::print(size_t index) {
4651 // cout << to_str(index);
4652 // }
4653 //
4654 // void func_::print(){
4655 // string str;
4656 // if (is_constant()) {
4657 // str += " (Constant) : ";
4658 // }
4659 // else if (is_complex()) {
4660 // cout << " (Complex) : ";
4661 // }
4662 // else if (is_linear()) {
4663 // str += " (Linear) : ";
4664 // }
4665 // else if (is_convex()) {
4666 // str += " (Convex) : ";
4667 // }
4668 // else if (is_concave()){
4669 // str += " (Concave) : ";
4670 // }
4671 // else {
4672 // str += " (Unknown) : ";
4673 // }
4674 // if (!_embedded && !is_constant()) {
4675 // str += "f(";
4676 // for (auto pair_it = _vars->begin(); pair_it != _vars->end();) {
4677 // // if (!pair_it->second.first->_is_vector) {
4678 // str += pair_it->second.first->get_name(false,true);
4679 // // str += pair_it->second.first->get_name(false,false)+"[";
4680 // // str += to_string(pair_it->second.first->get_id_inst())+"]";
4681 // if (next(pair_it) != _vars->end()) {
4682 // str += ",";
4683 // }
4684 // // }
4685 // pair_it++;
4686 // }
4687 // str += ") = ";
4688 // }
4689 // auto space_size = str.size();
4690 // auto nb_inst = _dim[0];
4691 // allocate_mem();
4692 // for (size_t inst = 0; inst<nb_inst; inst++) {
4693 // eval(inst);
4694 // if (inst>0) {
4695 // str.insert(str.end(), space_size, ' ');
4696 // }
4697 // str += to_str(inst);
4698 // str += "\n";
4699 // }
4700 // str += "\n";
4701 // cout << str;
4702 // }
4703 //
4704 //
4705 // void func_::update_to_str(bool input){
4706 // _to_str = to_str();
4707 // }
4708
4709
4710 /**
4711 Returns a pointer to the constant part of the function.
4712 @return a pointer to the constant part of the function.
4713 */
get_cst() const4714 shared_ptr<constant_> func_::get_cst() const{
4715 return _cst;
4716 }
4717
4718 /**
4719 Returns a pointer to the variable matching the name provided.
4720 @param[in] name variable name.
4721 @return a pointer to the variable matching the name provided.
4722 */
get_var(const string & name) const4723 shared_ptr<param_> func_::get_var(const string& name) const{
4724 if (_vars->empty()) {
4725 return nullptr;
4726 }
4727 auto pair_it = _vars->find(name);
4728 if (pair_it==_vars->end()) {
4729 return nullptr;
4730 }
4731 else {
4732 return get<1>(*pair_it).first;
4733 }
4734 }
4735
4736 /**
4737 Returns a pointer to the parameter matching the name provided.
4738 @param[in] name variable name.
4739 @return a pointer to the parameter matching the name provided.
4740 */
get_param(string name) const4741 shared_ptr<param_> func_::get_param(string name) const{
4742 auto pair_it = _params->find(name);
4743 if (pair_it==_params->end()) {
4744 return nullptr;
4745 }
4746 else {
4747 return get<1>(*pair_it).first;
4748 }
4749 }
4750
add_var(shared_ptr<param_> v,int nb)4751 void func_::add_var(shared_ptr<param_> v, int nb){/**< Inserts the variable in this function input list. nb represents the number of occurences v has. WARNING: Assumes that v has not been added previousely!*/
4752 if (_vars->count(v->get_name(false,false))!=0) {
4753 throw invalid_argument("In function add_var(v,nb): Variable already contained in function");
4754 }
4755 // if(_vars->empty()){
4756 // _indices = v->_indices;
4757 // }
4758 _vars->insert(make_pair<>(v->get_name(false,false), make_pair<>(v, nb)));
4759 if (v->_is_vector) {// i.e., it appears in a sum
4760 if (v->is_matrix()) {
4761 if (v->_is_transposed) {
4762 _nb_vars += v->get_dim(0);
4763 }
4764 _nb_vars += v->get_dim(1);
4765 }
4766 else {
4767 _nb_vars += v->get_dim();
4768 }
4769 }
4770 else {
4771 _nb_vars++;
4772 }
4773 }
4774
get_mag(constant_ * c)4775 func<double> get_mag(constant_* c){
4776 if(c->is_double()){
4777 func<> res = *((constant<double>*)c);
4778 return res;
4779 }
4780 if(c->is_complex()){
4781 func<> res = std::abs(((constant<Cpx>*)c)->eval());
4782 if(c->_is_transposed){
4783 res.transpose();
4784 }
4785 return res;
4786 }
4787 if(c->is_param() || c->is_var()){
4788 auto p = (param_*)c;
4789 if(p->_is_imag || p->_is_angle || p->_is_real || p->_is_sqrmag){
4790 throw invalid_argument("unsupported");
4791 }
4792 if(p->get_intype()==complex_){
4793 if(p->_mag){
4794 if(p->is_var()){
4795 return *((var<>*)p->_mag.get());
4796 }
4797 else {
4798 return *((param<>*)p->_mag.get());
4799 }
4800 }
4801 else {
4802 return 0;
4803 }
4804 }
4805 else if(p->get_intype()==double_){
4806 if(p->is_var()){
4807 return *((var<>*)p);
4808 }
4809 else {
4810 return *((param<>*)p);
4811 }
4812 }
4813 }
4814 if(c->is_function()){
4815 auto f = (func_*)c;
4816 if(f->get_return_type()==double_){
4817 return *(func<>*)f;
4818 }
4819 if(f->get_return_type()==complex_){
4820 return get_mag_ang(*(func<Cpx>*)f).first;
4821 }
4822 }
4823 throw invalid_argument("unsupported");
4824 }
4825
get_real(constant_ * c)4826 func<double> get_real(constant_* c){
4827 if(c->is_double()){
4828 func<> res = ((constant<double>*)c)->eval();
4829 if(c->_is_transposed){
4830 res.transpose();
4831 }
4832 return res;
4833 }
4834 if(c->is_complex()){
4835 func<> res = ((constant<Cpx>*)c)->eval().real();
4836 if(c->_is_transposed){
4837 res.transpose();
4838 }
4839 return res;
4840 }
4841 if(c->is_param() || c->is_var()){
4842 auto p = (param_*)c;
4843 if(p->_is_imag || p->_is_angle || p->_is_real || p->_is_sqrmag){
4844 throw invalid_argument("unsupported");
4845 }
4846 if(p->get_intype()==complex_){
4847 if(p->_real){
4848 if(p->is_var()){
4849 return *((var<>*)p->_real.get());
4850 }
4851 else {
4852 return *((param<>*)p->_real.get());
4853 }
4854 }
4855 else {
4856 return 0;
4857 }
4858 }
4859 else if(p->get_intype()==double_){
4860 if(p->is_var()){
4861 return *((var<>*)p);
4862 }
4863 else {
4864 return *((param<>*)p);
4865 }
4866 }
4867 }
4868 if(c->is_function()){
4869 auto f = (func_*)c;
4870 if(f->get_return_type()==double_){
4871 return *(func<>*)f;
4872 }
4873 if(f->get_return_type()==complex_){
4874 return get_real_imag(*(func<Cpx>*)f).first;
4875 }
4876 }
4877 throw invalid_argument("unsupported");
4878 }
4879
get_ang(constant_ * c)4880 func<double> get_ang(constant_* c){
4881 if(c->is_double()){
4882 return std::atan2(0,((constant<double>*)c)->eval());
4883 }
4884 if(c->is_complex()){
4885 return arg(((constant<Cpx>*)c)->eval());
4886 }
4887 if(c->is_param() || c->is_var()){
4888 auto p = (param_*)c;
4889 if(p->_is_angle){
4890 if(p->is_var()){
4891 return *((var<>*)p);
4892 }
4893 else {
4894 return *((param<>*)p);
4895 }
4896 }
4897 if(p->_is_imag || p->_is_real || p->_is_sqrmag){
4898 throw invalid_argument("unsupported");
4899 }
4900 if(p->get_intype()==complex_){
4901 if(p->_ang){
4902 if(p->is_var()){
4903 if(p->_is_conjugate){
4904 func<> res = (*((var<>*)p->_ang.get()));
4905 res.reverse_sign();
4906 return res;
4907 }
4908 else {
4909 return *((var<>*)p->_ang.get());
4910 }
4911 }
4912 else {
4913 if(p->_is_conjugate){
4914 func<> res = (*((param<>*)p->_ang.get()));
4915 res.reverse_sign();
4916 return res;
4917 }
4918 else {
4919 return *((param<>*)p->_ang.get());
4920 }
4921 }
4922 }
4923 else {
4924 return 0;
4925 }
4926 }
4927 else {
4928 return 0;
4929 }
4930 }
4931 if(c->is_function()){
4932 auto f = (func_*)c;
4933 if(f->get_return_type()==complex_){
4934 return get_mag_ang(*(func<Cpx>*)f).second;
4935 }
4936 return 0;
4937 }
4938 throw invalid_argument("unsupported");
4939 }
4940
get_imag(constant_ * c)4941 func<double> get_imag(constant_* c){
4942 if(c->is_double()){
4943 return 0;
4944 }
4945 if(c->is_complex()){
4946 return ((constant<Cpx>*)c)->eval().imag();
4947 }
4948 if(c->is_param() || c->is_var()){
4949 auto p = (param_*)c;
4950 if(p->_is_imag || p->_is_angle || p->_is_real || p->_is_sqrmag){
4951 throw invalid_argument("unsupported");
4952 }
4953 if(p->get_intype()==complex_){
4954 if(p->_imag){
4955 if(p->is_var()){
4956 if(p->_is_conjugate){
4957 func<> res = (*((var<>*)p->_imag.get()));
4958 res.reverse_sign();
4959 return res;
4960 }
4961 else {
4962 return *((var<>*)p->_imag.get());
4963 }
4964 }
4965 else {
4966 if(p->_is_conjugate){
4967 func<> res = (*((param<>*)p->_imag.get()));
4968 res.reverse_sign();
4969 return res;
4970 }
4971 else {
4972 return *((param<>*)p->_imag.get());
4973 }
4974 }
4975 }
4976 else {
4977 return 0;
4978 }
4979 }
4980 else {
4981 return 0;
4982 }
4983 }
4984 if(c->is_function()){
4985 auto f = (func_*)c;
4986 if(f->get_return_type()==complex_){
4987 return get_real_imag(*(func<Cpx>*)f).second;
4988 }
4989 return 0;
4990 }
4991 throw invalid_argument("unsupported");
4992 }
4993
4994
add_param(shared_ptr<param_> p,int nb)4995 void func_::add_param(shared_ptr<param_> p, int nb){/**< Inserts the parameter in this function input list. WARNING: Assumes that p has not been added previousely!*/
4996 if (_params->count(p->get_name(false,false))!=0) {
4997 throw invalid_argument("In function add_param(v,nb): parameter already contained in function");
4998 }
4999 // set_max_dim(*p);
5000 _params->insert(make_pair<>(p->get_name(false,false), make_pair<>(p, nb)));
5001 }
5002
5003 template<>
get_OA_symbolic(const vector<param<double>> & c,const param<double> & c0,const indices & Pert)5004 func<double> func<double>:: get_OA_symbolic(const vector<param<double>>& c, const param<double>& c0, const indices& Pert){ /**< Returns an outer-approximation of the function using the current value of the variables **/
5005 func<double> res; // res = gradf(x*)*(x-x*) + f(x*)
5006 int j=0;
5007 res=0;
5008 for(auto &it: *_vars)
5009 {
5010 auto v = it.second.first;
5011 auto v_ids = indices(Pert,*v->_indices);
5012 switch (v->get_intype()) {
5013 case binary_:
5014 res += c[j]*(static_pointer_cast<var<bool>>(v)->from_ith(2,v_ids));
5015 break;
5016 case short_:
5017 res += c[j]*(static_pointer_cast<var<short>>(v)->from_ith(2,v_ids));
5018 break;
5019 case integer_:
5020 res += c[j]*(static_pointer_cast<var<int>>(v)->from_ith(2,v_ids));
5021 break;
5022 case float_:
5023 res += c[j]*(static_pointer_cast<var<float>>(v)->from_ith(2,v_ids));
5024 break;
5025 break;
5026 case double_:
5027 res += c[j]*(static_pointer_cast<var<double>>(v)->from_ith(2,v_ids));
5028 break;
5029 default:
5030 break;
5031 }
5032 j++;
5033 }
5034 res += c0;
5035 res.in(indices(Pert,*this->_indices));
5036 merge_vars(res);
5037 return res;
5038 }
5039
5040
5041 // template<>
5042 // double func<double>::eval_double(size_t i) {
5043 //// if(is_zero()){
5044 //// return _range->first;
5045 //// }
5046 // if (_vars->empty() && _evaluated) {
5047 // if (_params->empty()){
5048 // return _val->at(0);
5049 // }
5050 // return _val->at(i);
5051 // }
5052 // double res = 0;
5053 //// if(!_cst->is_zero())
5054 // res += eval_cst(i);
5055 //// if(!_lterms->empty()){
5056 // for (auto &pair:*_lterms) {
5057 // if ((pair.second._coef->_is_transposed || pair.second._coef->is_matrix() || (pair.second._p->is_indexed() && pair.second._p->_indices->_ids->size()>1)) && !pair.second._p->is_matrix()) {
5058 // auto dim = pair.second._p->get_dim(i);
5059 // if (pair.second._sign) {
5060 // for (size_t j = 0; j<dim; j++) {
5061 // res += eval_coef(pair.second._coef,i,j) * ((param<double>*)(pair.second._p.get()))->eval(i,j);
5062 // }
5063 // }
5064 // else {
5065 // for (size_t j = 0; j<dim; j++) {
5066 // res -= eval_coef(pair.second._coef,i,j) * ((param<double>*)(pair.second._p.get()))->eval(i,j);
5067 // }
5068 // }
5069 // }
5070 // else {
5071 // if (pair.second._sign) {
5072 // res += eval_coef(pair.second._coef,i) * ((param<double>*)(pair.second._p.get()))->eval(i);
5073 // }
5074 // else {
5075 // res -= eval_coef(pair.second._coef,i) * ((param<double>*)(pair.second._p.get()))->eval(i);
5076 // }
5077 // }
5078 // }
5079 //// }
5080 // // res += eval_lterms(i);
5081 //// if(!_qterms->empty()){
5082 // for (auto &pair:*_qterms) {
5083 // double qval = 0;
5084 // if (pair.second._coef_p1_tr) { // qterm = (coef*p1)^T*p2
5085 // assert(pair.second._p->first->_dim[1]==1 && pair.second._coef->_dim[0]==pair.second._p->second->_dim[0]);
5086 // for (auto i = 0; i<pair.second._p->first->_dim[0]; i++) {
5087 // for (auto j = 0; j<pair.second._p->first->_dim[0]; j++) {
5088 // qval += eval_coef(pair.second._coef,i,j) * ((param<double>*)(pair.second._p->first.get()))->eval(i) * ((param<double>*)(pair.second._p->second.get()))->eval(j);
5089 // }
5090 // }
5091 // }
5092 // else if (pair.second._p->first->is_matrix() && !pair.second._p->second->is_matrix() && !pair.second._p->second->_is_transposed) {//matrix * vect
5093 // for (size_t j = 0; j<pair.second._p->second->_dim[0]; j++) {
5094 // qval += ((param<double>*)(pair.second._p->first.get()))->eval(i,j) * ((param<double>*)(pair.second._p->second.get()))->eval(j);
5095 // }
5096 // qval *= eval_coef(pair.second._coef,i);
5097 // }
5098 // else if (!pair.second._p->first->is_matrix() && pair.second._p->first->_is_transposed && pair.second._p->second->is_matrix() ) {//transposed vect * matrix
5099 // for (size_t j = 0; j<pair.second._p->first->_dim[0]; j++) {
5100 // qval += ((param<double>*)(pair.second._p->first.get()))->eval(j) * ((param<double>*)(pair.second._p->second.get()))->eval(j,i);
5101 // }
5102 // qval *= eval_coef(pair.second._coef,i);
5103 // }
5104 // else if (!pair.second._p->first->is_matrix() && pair.second._p->first->_is_transposed && !pair.second._p->second->is_matrix() && i==0) {//transposed vect * vec, a dot product of two vectors
5105 // for (size_t j = 0; j<pair.second._p->first->_dim[1]; j++) {
5106 // qval += ((param<double>*)(pair.second._p->first.get()))->eval(j) * ((param<double>*)(pair.second._p->second.get()))->eval(j);
5107 // }
5108 // qval *= eval_coef(pair.second._coef,i);
5109 // }
5110 // else if (!pair.second._coef->is_matrix() && pair.second._coef->_is_transposed && !pair.second._p->first->is_matrix()) {//transposed vect * vec, a dot product of two vectors
5111 // for (size_t j = 0; j<pair.second._p->first->_dim[0]; j++) {
5112 // qval += eval_coef(pair.second._coef,j) * ((param<double>*)(pair.second._p->first.get()))->eval(j) * ((param<double>*)(pair.second._p->second.get()))->eval(j);
5113 // }
5114 // }
5115 // else {
5116 // qval += eval_coef(pair.second._coef,i) * ((param<double>*)(pair.second._p->first.get()))->eval(i) * ((param<double>*)(pair.second._p->second.get()))->eval(i);
5117 // }
5118 // if (!pair.second._sign) {
5119 // qval *= -1;
5120 // }
5121 // res += qval;
5122 // }
5123 //// }
5124 // // res += eval_qterms(i);
5125 //// if(!_pterms->empty()){
5126 // for (auto &pair:*_pterms) {
5127 // double pval = 1;
5128 // for (auto &vpair: *pair.second._l) {
5129 // pval *= std::pow(eval(vpair.first, i), vpair.second);
5130 // }
5131 // pval *= eval_coef(pair.second._coef,i);
5132 // if (!pair.second._sign) {
5133 // pval *= -1;
5134 // }
5135 // res += pval;
5136 // }
5137 // // res += eval_pterms(i);
5138 //// }
5139 // if(_expr)
5140 // res += eval_expr(_expr,i);
5141 // if (_vars->empty() && _params->empty()) {
5142 // _val->at(0) = res;
5143 // _evaluated = true;
5144 // }
5145 // else {
5146 // if (_vars->empty() && i==_val->size()-1) {
5147 // _evaluated = true;
5148 // }
5149 // _val->at(i) = res;
5150 // }
5151 // return res;
5152 // }
5153 //
5154 //
5155 //
5156 // void func_::delete_var(const string& vid){
5157 // auto vit = _vars->find(vid);
5158 // if (vit==_vars->end()) {
5159 // return;
5160 // }
5161 // _vars->erase(vit);
5162 // }
5163 //
5164 // void func_::delete_param(const string& vid){
5165 // auto vit = _params->find(vid);
5166 // if (vit==_params->end()) {
5167 // return;
5168 // }
5169 // _params->erase(vit);
5170 // }
5171 //
5172 //
5173 //
5174 //
nb_occ_var(string name) const5175 unsigned func_::nb_occ_var(string name) const{/**< Returns the number of occurences the variable has in this function. */
5176 auto pair_it = _vars->find(name);
5177 if (pair_it==_vars->end()) {
5178 return 0;
5179 }
5180 else {
5181 return get<1>(*pair_it).second;
5182 }
5183 }
5184
nb_occ_param(string name) const5185 unsigned func_::nb_occ_param(string name) const{/**< Returns the number of occurences the parameter has in this function. */
5186 auto pair_it = _params->find(name);
5187 if (pair_it==_params->end()) {
5188 return 0;
5189 }
5190 else {
5191 return get<1>(*pair_it).second;
5192 }
5193 }
5194
incr_occ_var(string str)5195 void func_::incr_occ_var(string str){/**< Increases the number of occurences the variable has in this function. */
5196 auto pair_it = _vars->find(str);
5197 if (pair_it==_vars->end()) {
5198 throw invalid_argument("Non-existing variable in function!\n");
5199 }
5200 else {
5201 get<1>(*pair_it).second++;
5202 }
5203 }
5204
incr_occ_param(string str)5205 void func_::incr_occ_param(string str){/**< Increases the number of occurences the parameter has in this function. */
5206 auto pair_it = _params->find(str);
5207 if (pair_it==_params->end()) {
5208 throw invalid_argument("Non-existing variable in function!\n");
5209 }
5210 else {
5211 get<1>(*pair_it).second++;
5212 }
5213 }
5214
decr_occ_var(string str,int nb)5215 void func_::decr_occ_var(string str, int nb){/**< Decreases the number of occurences the variable has in this function by nb. */
5216 auto pair_it = _vars->find(str);
5217 if (pair_it==_vars->end()) {
5218 return;
5219 }
5220 else {
5221 get<1>(*pair_it).second-=nb;
5222 if (get<1>(*pair_it).second==0) {
5223 _vars->erase(pair_it);
5224 }
5225 }
5226 }
5227
decr_occ_param(string str,int nb)5228 void func_::decr_occ_param(string str, int nb){/**< Decreases the number of occurences the parameter has in this function by nb. */
5229 auto pair_it = _params->find(str);
5230 if (pair_it==_params->end()) {
5231 return;
5232 }
5233 else {
5234 get<1>(*pair_it).second -= nb;
5235 if (get<1>(*pair_it).second==0) {
5236 _params->erase(pair_it);
5237 }
5238 }
5239 }
5240
5241
5242
5243
5244
is_convex(size_t idx) const5245 bool func_::is_convex(size_t idx) const{
5246 return (_convexity->at(idx)==convex_ || _convexity->at(idx)==linear_);
5247 };
5248
is_concave(size_t idx) const5249 bool func_::is_concave(size_t idx) const{
5250 return (_convexity->at(idx)==concave_ || _convexity->at(idx)==linear_);
5251 };
5252
5253
5254 // bool func_::is_number() const{
5255 // return (_vars->empty() && _params->empty());
5256 // }
5257
5258
5259
is_linear() const5260 bool func_::is_linear() const{
5261 return (_ftype==lin_);
5262 };
5263
is_quadratic() const5264 bool func_::is_quadratic() const{
5265 return (_ftype==quad_);
5266 };
5267
is_polynomial() const5268 bool func_::is_polynomial() const{
5269 return (_ftype==pol_);
5270 };
5271
is_nonlinear() const5272 bool func_::is_nonlinear() const{
5273 return (_ftype==nlin_);
5274 };
5275
5276 // bool func_::is_complex() const{
5277 // for(auto &it: *_vars){
5278 // if (it.second.first->is_complex()) {
5279 // return true;
5280 // }
5281 // }
5282 // for(auto &it: *_params){
5283 // if (it.second.first->is_complex()) {
5284 // return true;
5285 // }
5286 // }
5287 // return false;
5288 // };
5289 //
5290 //
5291 // bool func_::is_zero() const{/*<< A function is zero if it is constant and equals zero or if it is a sum of zero valued parameters */
5292 // if(is_number() && !_is_vector && !_is_transposed && t_eval(this)==0){
5293 // return true;
5294 // }
5295 // if (_ftype==const_ && _cst->is_zero()){
5296 // for (auto& it:*_params) {
5297 // if (!it.second.first->is_zero()) {
5298 // return false;
5299 // }
5300 // }
5301 // return true;
5302 // }
5303 // return false;
5304 // }
5305 //
5306 // bool func_::is_transposed() const {
5307 // return _is_transposed;
5308 // }
5309 //
get_ftype() const5310 FType func_::get_ftype() const { return _ftype;}
5311
5312 //
5313 //
get_square(shared_ptr<param_> p)5314 qterm* func_::get_square(shared_ptr<param_> p){ /**< Returns the quadratic term containing a square of p or nullptr if none exists. **/
5315 for (auto pair_it = _qterms->begin(); pair_it != _qterms->end(); pair_it++) {
5316 if (pair_it->second._p->first==p && pair_it->second._p->second==p) {
5317 return &pair_it->second;
5318 }
5319 }
5320 return nullptr;
5321 }
5322 //
5323 // func_ func_::get_outer_app(){ /**< Returns an outer-approximation of the function using the current value of the variables **/
5324 // func_ res; // res = gradf(x*)*(x-x*) + f(x*)
5325 // param_* v;
5326 // for(auto &it: *_vars){
5327 // v = it.second.first.get();
5328 // res += (get_stored_derivative(v->_name)->eval())*((*v) - t_eval(v));
5329 // }
5330 // res += eval();
5331 // return res;
5332 // }
5333 //
5334 //
5335 // pair<double,double>* func_::get_all_range() const{
5336 // return _all_range;
5337 // }
5338 //
5339 // Sign func_::get_all_sign() const{
5340 // return _all_sign;
5341 // }
5342 //
5343 // Sign func_::get_sign(size_t idx) const{
5344 // return _sign->at(idx);
5345 // }
5346 //
5347 // Sign func_::get_all_sign(const lterm& l) {
5348 // if (l._coef->is_zero()) {
5349 // return zero_;
5350 // }
5351 // if (l._coef->get_all_sign()==unknown_ || l._p->get_all_sign()==unknown_) {
5352 // return unknown_;
5353 // }
5354 // auto s = l._coef->get_all_sign() * l._p->get_all_sign();
5355 // if(s == 1 || s == 2) {
5356 // if (l._sign) {
5357 // return non_neg_;
5358 // }
5359 // else {
5360 // return non_pos_;
5361 // }
5362 // }
5363 // if(s == 4) {
5364 // if (l._sign) {
5365 // return pos_;
5366 // }
5367 // else {
5368 // return neg_;
5369 // }
5370 // }
5371 // if(s == -1 || s == -2) {
5372 // if (l._sign) {
5373 // return non_pos_;
5374 // }
5375 // else{
5376 // return non_neg_;
5377 // }
5378 // }
5379 // if(s == -4) {
5380 // if (l._sign) {
5381 // return neg_;
5382 // }
5383 // else {
5384 // return pos_;
5385 // }
5386 // }
5387 // return unknown_;
5388 // }
5389 //
5390 // Sign func_::get_all_sign(const qterm& l) {
5391 // if (l._coef->is_zero()) {
5392 // return zero_;
5393 // }
5394 // if (l._coef->get_all_sign()==unknown_ || l._p->first->get_all_sign()==unknown_ || l._p->second->get_all_sign()==unknown_) {
5395 // return unknown_;
5396 // }
5397 // auto s = l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign();
5398 // if(s == 1 || s == 2 || s == 4) {
5399 // if (l._sign) {
5400 // return non_neg_;
5401 // }
5402 // else {
5403 // return non_pos_;
5404 // }
5405 // }
5406 // if(s == 8) {
5407 // if (l._sign) {
5408 // return pos_;
5409 // }
5410 // else {
5411 // return neg_;
5412 // }
5413 // }
5414 // if(s == -1 || s == -2 || s == -4) {
5415 // if (l._sign) {
5416 // return non_pos_;
5417 // }
5418 // else{
5419 // return non_neg_;
5420 // }
5421 // }
5422 // if(s == -8) {
5423 // if (l._sign) {
5424 // return neg_;
5425 // }
5426 // else {
5427 // return pos_;
5428 // }
5429 // }
5430 // return unknown_;
5431 // }
5432 //
5433 // Sign func_::get_all_sign(const pterm& l) {
5434 // if (l._coef->is_zero()) {
5435 // return zero_;
5436 // }
5437 // // if (l._coef->get_all_sign()==unknown_ || l._p->first->get_all_sign()==unknown_ || l._p->second->get_all_sign()==unknown_) {
5438 // // return unknown_;
5439 // // }
5440 // // if (l._sign) {
5441 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == 2) {
5442 // // return non_neg_;
5443 // // }
5444 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == 4) {
5445 // // return pos_;
5446 // // }
5447 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == -2) {
5448 // // return non_pos_;
5449 // // }
5450 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == -4) {
5451 // // return neg_;
5452 // // }
5453 // // }
5454 // // else {
5455 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == 2) {
5456 // // return non_pos_;
5457 // // }
5458 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == 4) {
5459 // // return neg_;
5460 // // }
5461 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == -2) {
5462 // // return non_neg_;
5463 // // }
5464 // // if(l._coef->get_all_sign() * l._p->first->get_all_sign() * l._p->second->get_all_sign() == -4) {
5465 // // return pos_;
5466 // // }
5467 // // }
5468 // return unknown_;
5469 // }
5470 //
5471 //
5472 //
5473
5474
5475
5476 /**
5477 Index the function and its variables/parameters using nodes of a graph
5478 @param[in] vec vector of nodes
5479 @return current function
5480 */
5481 // func_& func_::in(const vector<Node*>& vec){
5482 // _nb_vars = 0;
5483 // string key;
5484 // auto iter = _vars->begin();
5485 // while (iter!=_vars->end()) {
5486 // auto pair = (*iter++);
5487 // auto v = pair.second.first;
5488 // if(!v->_indices){
5489 // v->index_in(vec);
5490 // }
5491 // else if(v->_indices->_type==in_arcs_){
5492 // v->index_in_arcs(vec);
5493 // }
5494 // else if(v->_indices->_type==out_arcs_){
5495 // v->index_out_arcs(vec);
5496 // }
5497 // else if(v->_indices->_type==in_gens_){
5498 // v->index_in_aux(vec,"gens");
5499 // }
5500 // if (!v->_is_vector) {// i.e., it is not transposed
5501 // _nb_vars++;
5502 // }
5503 // else {
5504 // _nb_vars += v->get_dim();
5505 // }
5506 // }
5507 // iter = _params->begin();
5508 // while (iter!=_params->end()) {
5509 // auto pair = (*iter++);
5510 // auto p = pair.second.first;
5511 // if(!p->_indices){
5512 // p->index_in(vec);
5513 // }
5514 // else if(p->_indices->_type==in_arcs_){
5515 // p->index_in_arcs(vec);
5516 // }
5517 // else if(p->_indices->_type==out_arcs_){
5518 // p->index_out_arcs(vec);
5519 // }
5520 // else if(p->_indices->_type==in_gens_){
5521 // p->index_in_aux(vec,"gens");
5522 // }
5523 // }
5524 // return *this;
5525 // }
5526
5527 /**
5528 Index the function and its variables/parameters using the indices in ids
5529 @param[in] ids indices
5530 @return current function
5531 */
5532 // func_& func_::in(const indices& ids) {
5533 // _nb_vars = 0;
5534 // string key;
5535 // auto iter = _vars->begin();
5536 // while (iter!=_vars->end()) {
5537 // auto pair = (*iter++);
5538 // auto v = pair.second.first;
5539 // if(!v->is_indexed()){
5540 // v->index_in(ids);
5541 // }
5542 // if (!v->_is_vector) {// i.e., it is not transposed
5543 // _nb_vars++;
5544 // }
5545 // else {
5546 // _nb_vars += v->get_dim();
5547 // }
5548 // }
5549 // iter = _params->begin();
5550 // while (iter!=_params->end()) {
5551 // auto pair = (*iter++);
5552 // auto p = pair.second.first;
5553 // if(!p->is_indexed()){
5554 // p->index_in(ids);
5555 // }
5556 // }
5557 // _indices = make_shared<indices>(ids);
5558 // return *this;
5559 // }
5560 /** The function iterates over key references in _ids and keeps only the unique entries */
keep_unique_keys()5561 void func_::keep_unique_keys() {
5562 if(_indices){
5563 _indices->keep_unique_keys();
5564 }
5565 }
5566
5567 /**
5568 Update the variables/parameters indexing using the set ids and the references in old_ids.
5569 @param[in] ids indext set
5570 */
5571 template<typename type>
update_var_indices(const indices & ids,const indices & old_ids)5572 void func<type>::update_var_indices(const indices& ids, const indices& old_ids){
5573 _nb_vars = 0;
5574 string key;
5575 auto iter = _vars->begin();
5576 while (iter!=_vars->end()) {
5577 auto pair = (*iter++);
5578 auto p = pair.second.first;
5579 p->match_matrix_ids(ids,old_ids);
5580 }
5581 iter = _params->begin();
5582 while (iter!=_params->end()) {
5583 auto pair = (*iter++);
5584 auto p = pair.second.first;
5585 p->match_matrix_ids(ids,old_ids);
5586 }
5587 for (auto &pair:*_lterms) {
5588 auto coef = pair.second._coef;
5589 if (coef->is_function()) {
5590 auto c = static_pointer_cast<func<type>>(coef);
5591 c->update_var_indices(ids,old_ids);
5592 }
5593 }
5594 for (auto &pair:*_qterms) {
5595 auto coef = pair.second._coef;
5596 if (coef->is_function()) {
5597 auto c = static_pointer_cast<func<type>>(coef);
5598 c->update_var_indices(ids,old_ids);
5599 }
5600 }
5601 for (auto &pair:*_pterms) {
5602 auto coef = pair.second._coef;
5603 if (coef->is_function()) {
5604 auto c = static_pointer_cast<func<type>>(coef);
5605 c->update_var_indices(ids,old_ids);
5606 }
5607 }
5608 this->index_in(range(1,ids.size()));
5609 /* TODO: nonlinear part*/
5610 }
5611
5612 /**
5613 Update the variables/parameters indexing using the set ids.
5614 @param[in] ids indext set
5615 */
5616 template<typename type>
update_var_indices(const indices & ids)5617 void func<type>::update_var_indices(const indices& ids){
5618 _nb_vars = 0;
5619 string key;
5620 auto iter = _vars->begin();
5621 while (iter!=_vars->end()) {
5622 auto pair = (*iter++);
5623 auto p = pair.second.first;
5624 switch (p->get_intype()) {
5625 case binary_:
5626 static_pointer_cast<var<bool>>(p)->index_in(ids);
5627 break;
5628 case short_:
5629 static_pointer_cast<var<short>>(p)->index_in(ids);
5630 break;
5631 case integer_:
5632 static_pointer_cast<var<int>>(p)->index_in(ids);
5633 break;
5634 case float_:
5635 static_pointer_cast<var<float>>(p)->index_in(ids);
5636 break;
5637 case double_:
5638 static_pointer_cast<var<double>>(p)->index_in(ids);
5639 break;
5640 case long_:
5641 static_pointer_cast<var<long double>>(p)->index_in(ids);
5642 break;
5643 case complex_:
5644 static_pointer_cast<var<Cpx>>(p)->index_in(ids);
5645 break;
5646 default:
5647 break;
5648 }
5649 }
5650 iter = _params->begin();
5651 while (iter!=_params->end()) {
5652 auto pair = (*iter++);
5653 auto p = pair.second.first;
5654 switch (p->get_intype()) {
5655 case binary_:
5656 static_pointer_cast<param<bool>>(p)->index_in(ids);
5657 break;
5658 case short_:
5659 static_pointer_cast<param<short>>(p)->index_in(ids);
5660 break;
5661 case integer_:
5662 static_pointer_cast<param<int>>(p)->index_in(ids);
5663 break;
5664 case float_:
5665 static_pointer_cast<param<float>>(p)->index_in(ids);
5666 break;
5667 case double_:
5668 static_pointer_cast<param<double>>(p)->index_in(ids);
5669 break;
5670 case long_:
5671 static_pointer_cast<param<long double>>(p)->index_in(ids);
5672 break;
5673 case complex_:
5674 static_pointer_cast<param<Cpx>>(p)->index_in(ids);
5675 break;
5676 default:
5677 break;
5678 }
5679 }
5680 for (auto &pair:*_lterms) {
5681 auto coef = pair.second._coef;
5682 if (coef->is_function()) {
5683 auto c = static_pointer_cast<func<type>>(coef);
5684 c->index_in(ids);
5685 }
5686 }
5687 for (auto &pair:*_qterms) {
5688 auto coef = pair.second._coef;
5689 if (coef->is_function()) {
5690 auto c = static_pointer_cast<func<type>>(coef);
5691 c->index_in(ids);
5692 }
5693 }
5694 for (auto &pair:*_pterms) {
5695 auto coef = pair.second._coef;
5696 if (coef->is_function()) {
5697 auto c = static_pointer_cast<func<type>>(coef);
5698 c->index_in(ids);
5699 }
5700 }
5701 /* TODO: nonlinear part*/
5702 }
5703 /**
5704 Update the function indexing and its variables/parameters using the keep_ids vector of bool, only keep an index if it corresponding entry in keep_id is true.
5705 @param[in] keep_ids vector of booleans, specifying which ids to keep
5706 @return current function
5707 */
5708 template<typename type>
update_indices(const vector<bool> & keep_ids)5709 void func<type>::update_indices(const vector<bool>& keep_ids) {
5710 _nb_vars = 0;
5711 string key;
5712 bool has_matrix_ids = has_matrix_indexed_vars();
5713 auto iter = _vars->begin();
5714 while (iter!=_vars->end()) {
5715 auto pair = (*iter++);
5716 auto v = pair.second.first;
5717 if(has_matrix_ids && !v->is_matrix_indexed()){
5718 // v->_indices->to_matrix();
5719 v->_indices->filter_matrix_rows(keep_ids);
5720 }
5721 else {
5722 v->_indices->filter_refs(keep_ids);
5723 }
5724 if (!v->_is_vector) {// TODO: matrix indexed case
5725 _nb_vars++;
5726 }
5727 else {
5728 _nb_vars += v->get_dim();
5729 }
5730 }
5731 iter = _params->begin();
5732 while (iter!=_params->end()) {
5733 auto pair = (*iter++);
5734 auto p = pair.second.first;
5735 if(has_matrix_ids && !p->is_matrix_indexed()){
5736 // p->_indices->to_matrix();
5737 p->_indices->filter_matrix_rows(keep_ids);
5738 }
5739 else {
5740 p->_indices->filter_refs(keep_ids);
5741 }
5742 }
5743 if(_indices) {
5744 if(has_matrix_ids && !_indices->is_matrix_indexed()){
5745 _indices->filter_matrix_rows(keep_ids);
5746 }
5747 else {
5748 _indices->filter_refs(keep_ids);
5749 }
5750 _dim[0] = _indices->size();
5751 }
5752 auto f_cpy = *this;
5753 bool update = false;
5754 for (auto &pair:*_lterms) {
5755 if(pair.second._p->is_matrix_indexed() && pair.second._p->_indices->nb_keys()==0){
5756 f_cpy.insert(!pair.second._sign, *pair.second._coef, *pair.second._p);
5757 update = true;
5758 }
5759 }
5760 for (auto &pair:*_qterms) {
5761 if((pair.second._p->first->is_matrix_indexed() && pair.second._p->first->_indices->nb_keys()==0) || (pair.second._p->first->is_matrix_indexed() && pair.second._p->first->_indices->nb_keys()==0)){
5762 f_cpy.insert(!pair.second._sign, *pair.second._coef, *pair.second._p->first, *pair.second._p->second);
5763 update = true;
5764 }
5765 }
5766 // for (auto &pair:*_pterms) {
5767 // insert(pair.second);
5768 // }
5769 if(update){
5770 *this = f_cpy;
5771 }
5772 }
5773
5774
5775 /**
5776 Update the function indexing and its variables/parameters using the keep_ids vector of bool, only keep an index if it corresponding entry in keep_id is true.
5777 @param[in] keep_ids vector of booleans, specifying which ids to keep
5778 @return current function
5779 */
5780 template<typename type>
update_rows(const vector<bool> & keep_ids)5781 void func<type>::update_rows(const vector<bool>& keep_ids) {
5782 _nb_vars = 0;
5783 string key;
5784 auto iter = _vars->begin();
5785 while (iter!=_vars->end()) {
5786 auto pair = (*iter++);
5787 auto v = pair.second.first;
5788 if(v->get_nb_inst()>1)
5789 v->_indices->filter_rows(keep_ids);/* TODO: update name */
5790 if (!v->_is_vector) {// i.e., it is not transposed
5791 _nb_vars++;
5792 }
5793 else {
5794 _nb_vars += v->get_dim();
5795 }
5796 }
5797 iter = _params->begin();
5798 while (iter!=_params->end()) {
5799 auto pair = (*iter++);
5800 auto p = pair.second.first;
5801 if(p->get_nb_inst()>1)
5802 p->_indices->filter_rows(keep_ids);
5803 }
5804 if(_indices && get_nb_inst()>1) {
5805 _indices->filter_rows(keep_ids);
5806 _dim[0] = _indices->size();
5807 }
5808 auto f_cpy = *this;
5809 bool update = false;
5810 for (auto &pair:*_lterms) {
5811 if(pair.second._p->is_matrix_indexed() && pair.second._p->_indices->nb_keys()==0){
5812 f_cpy.insert(!pair.second._sign, *pair.second._coef, *pair.second._p);
5813 update = true;
5814 }
5815 }
5816 for (auto &pair:*_qterms) {
5817 if((pair.second._p->first->is_matrix_indexed() && pair.second._p->first->_indices->nb_keys()==0) || (pair.second._p->first->is_matrix_indexed() && pair.second._p->first->_indices->nb_keys()==0)){
5818 f_cpy.insert(!pair.second._sign, *pair.second._coef, *pair.second._p->first, *pair.second._p->second);
5819 update = true;
5820 }
5821 }
5822 for (auto &pair:*_pterms) {
5823 bool has_zero_var = false;
5824 for (auto &p: *pair.second._l) {
5825 if (p.first->is_matrix_indexed() && p.first->_indices->nb_keys()==0) {
5826 has_zero_var = true;
5827 break;
5828 }
5829 }
5830 if(has_zero_var){
5831 f_cpy.insert(!pair.second._sign, *pair.second._coef, *pair.second._l);
5832 update = true;
5833 }
5834 }
5835 if(update){
5836 *this = f_cpy;
5837 }
5838 }
5839
5840
5841
5842
5843
5844 //
5845 //
5846 // string poly_to_str(const constant_* c, size_t inst){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
5847 //
5848 // if (!c) {
5849 // return "null";
5850 // }
5851 // switch (c->get_type()) {
5852 // case binary_c: {
5853 // return ((constant<bool>*)(c))->to_str();
5854 // break;
5855 // }
5856 // case short_c: {
5857 // return ((constant<short>*)(c))->to_str();
5858 // break;
5859 // }
5860 // case integer_c: {
5861 // return ((constant<int>*)(c))->to_str();
5862 // break;
5863 // }
5864 // case float_c: {
5865 // return ((constant<float>*)(c))->to_str();
5866 // break;
5867 // }
5868 // case double_c: {
5869 // return ((constant<double>*)(c))->to_str();
5870 // break;
5871 // }
5872 // case long_c: {
5873 // return ((constant<long double>*)(c))->to_str();
5874 // break;
5875 // }
5876 // case complex_c: {
5877 // return ((constant<Cpx>*)(c))->to_str();
5878 // break;
5879 // }
5880 // case par_c:{
5881 // auto p_c = (param_*)(c);
5882 // switch (p_c->get_intype()) {
5883 // case binary_:
5884 // return ((param<bool>*)p_c)->to_str(inst);
5885 // break;
5886 // case short_:
5887 // return ((param<short>*)p_c)->to_str(inst);
5888 // break;
5889 // case integer_:
5890 // return ((param<int>*)p_c)->to_str(inst);
5891 // break;
5892 // case float_:
5893 // return ((param<float>*)p_c)->to_str(inst);
5894 // break;
5895 // case double_:
5896 // return ((param<double>*)p_c)->to_str(inst);
5897 // break;
5898 // case long_:
5899 // return ((param<long double>*)p_c)->to_str(inst);
5900 // break;
5901 // case complex_:
5902 // return ((param<Cpx>*)p_c)->to_str(inst);
5903 // break;
5904 // default:
5905 // break;
5906 // }
5907 // break;
5908 // }
5909 // case uexp_c: {
5910 // return ((uexpr*)c)->to_str(inst);
5911 // break;
5912 // }
5913 // case bexp_c: {
5914 // return ((bexpr*)c)->to_str(inst);
5915 // break;
5916 // }
5917 // case var_c: {
5918 // auto p_c = (param_*)(c);
5919 // switch (p_c->get_intype()) {
5920 // case binary_:
5921 // return ((var<bool>*)p_c)->get_name(inst);
5922 // break;
5923 // case short_:
5924 // return ((var<short>*)p_c)->get_name(inst);
5925 // break;
5926 // case integer_:
5927 // return ((var<int>*)p_c)->get_name(inst);
5928 // break;
5929 // case float_:
5930 // return ((var<float>*)p_c)->get_name(inst);
5931 // break;
5932 // case double_:
5933 // return ((var<double>*)p_c)->get_name(inst);
5934 // break;
5935 // case long_:
5936 // return ((var<long double>*)p_c)->get_name(inst);
5937 // break;
5938 // case complex_:
5939 // return ((var<Cpx>*)p_c)->get_name(inst);
5940 // break;
5941 // default:
5942 // break;
5943 // }
5944 // break;
5945 // }
5946 // case func_c: {
5947 // return ((func_*)c)->to_str(inst);
5948 // break;
5949 // }
5950 // default:
5951 // break;
5952 // }
5953 // return "null";
5954 // }
5955 //
5956 // string poly_to_str(const constant_* c, size_t inst1, size_t inst2){/**< printing c, detecting the right class, i.e., constant<>, param<>, uexpr or bexpr. */
5957 //
5958 // if (!c) {
5959 // return "null";
5960 // }
5961 // switch (c->get_type()) {
5962 // case binary_c: {
5963 // return ((constant<bool>*)(c))->to_str();
5964 // break;
5965 // }
5966 // case short_c: {
5967 // return ((constant<short>*)(c))->to_str();
5968 // break;
5969 // }
5970 // case integer_c: {
5971 // return ((constant<int>*)(c))->to_str();
5972 // break;
5973 // }
5974 // case float_c: {
5975 // return ((constant<float>*)(c))->to_str();
5976 // break;
5977 // }
5978 // case double_c: {
5979 // return ((constant<double>*)(c))->to_str();
5980 // break;
5981 // }
5982 // case long_c: {
5983 // return ((constant<long double>*)(c))->to_str();
5984 // break;
5985 // }
5986 // case complex_c: {
5987 // return ((constant<Cpx>*)(c))->to_str();
5988 // break;
5989 // }
5990 // case par_c:{
5991 // auto p_c = (param_*)(c);
5992 // switch (p_c->get_intype()) {
5993 // case binary_:
5994 // return ((param<bool>*)p_c)->to_str(inst1,inst2);
5995 // break;
5996 // case short_:
5997 // return ((param<short>*)p_c)->to_str(inst1,inst2);
5998 // break;
5999 // case integer_:
6000 // return ((param<int>*)p_c)->to_str(inst1,inst2);
6001 // break;
6002 // case float_:
6003 // return ((param<float>*)p_c)->to_str(inst1,inst2);
6004 // break;
6005 // case double_:
6006 // return ((param<double>*)p_c)->to_str(inst1,inst2);
6007 // break;
6008 // case long_:
6009 // return ((param<long double>*)p_c)->to_str(inst1,inst2);
6010 // break;
6011 // case complex_:
6012 // return ((param<Cpx>*)p_c)->to_str(inst1,inst2);
6013 // break;
6014 // default:
6015 // break;
6016 // }
6017 // break;
6018 // }
6019 // case uexp_c: {
6020 // return ((uexpr*)c)->to_str(inst2);
6021 // break;
6022 // }
6023 // case bexp_c: {
6024 // return ((bexpr*)c)->to_str(inst2);
6025 // break;
6026 // }
6027 // case var_c: {
6028 // auto p_c = (param_*)(c);
6029 // switch (p_c->get_intype()) {
6030 // case binary_:
6031 // return ((var<bool>*)p_c)->get_name(inst1,inst2);
6032 // break;
6033 // case short_:
6034 // return ((var<short>*)p_c)->get_name(inst1,inst2);
6035 // break;
6036 // case integer_:
6037 // return ((var<int>*)p_c)->get_name(inst1,inst2);
6038 // break;
6039 // case float_:
6040 // return ((var<float>*)p_c)->get_name(inst1,inst2);
6041 // break;
6042 // case double_:
6043 // return ((var<double>*)p_c)->get_name(inst1,inst2);
6044 // break;
6045 // case long_:
6046 // return ((var<long double>*)p_c)->get_name(inst1,inst2);
6047 // break;
6048 // case complex_:
6049 // return ((var<Cpx>*)p_c)->get_name(inst1,inst2);
6050 // break;
6051 // default:
6052 // break;
6053 // }
6054 // break;
6055 // }
6056 //
6057 // case func_c: {
6058 // if(c->is_matrix()){
6059 // return to_string(((func_*)c)->get_val(inst1,inst2));
6060 // }
6061 // return ((func_*)c)->to_str(inst2);
6062 // break;
6063 // }
6064 // default:
6065 // break;
6066 // }
6067 // return "null";
6068 // }
6069 //
6070 //
6071 // func_ conj(const func_& f){
6072 // func_ newf(f);
6073 // newf._is_conjugate = !newf._is_conjugate;
6074 // return newf;
6075 // }
6076 //
6077 // func_ ang(const func_& f){
6078 // func_ newf(f);
6079 // newf._is_angle = true;
6080 // return newf;
6081 // }
6082 //
6083 // func_ sqrmag(const func_& f){
6084 // func_ newf(f);
6085 // newf._is_sqrmag = true;
6086 // return newf;
6087 // }
6088 //
6089 // func_ real(const func_& f){
6090 // func_ newf(f);
6091 // newf._is_real = true;
6092 // return newf;
6093 // }
6094 //
6095 // func_ imag(const func_& f){
6096 // func_ newf(f);
6097 // newf._is_imag = true;
6098 // return newf;
6099 // }
6100 //
6101
6102 //
6103 // template func<double> constant<double>::get_real() const;
6104 // template func<double> constant<double>::get_imag() const;
6105 // template func<double> constant<Cpx>::get_real() const;
6106 // template func<double> constant<Cpx>::get_imag() const;
6107 // template func<double> var<double>::get_real() const;
6108 // template func<double> var<double>::get_imag() const;
6109 // template func<double> var<Cpx>::get_real() const;
6110 // template func<double> var<Cpx>::get_imag() const;
6111 // template <>
6112 // func<double> var<double>::get_real() const;
6113 // template <>
6114 // func<double> var<double>::get_imag() const;
6115 // template <>
6116 // func<double> param<double>::get_real() const;
6117 // template <>
6118 // func<double> param<double>::get_imag() const;
6119 //
6120 // template <>
6121 // func<double> constant<Cpx>::get_real() const;
6122 // template <>
6123 // func<double> constant<Cpx>::get_imag() const;
6124 // template <>
6125 // func<double> var<Cpx>::get_real() const;
6126 // template <>
6127 // func<double> var<Cpx>::get_imag() const;
6128 // template <>
6129 // func<double> param<Cpx>::get_real() const;
6130 // template <>
6131 // func<double> param<Cpx>::get_imag() const;
6132 template func<double> func<double>::replace<double>(const var<double>&, const func<double>&) const;
6133 template void func<double>::scale_coefs(double);
6134 template double func<double>::get_scale_factor(double);
6135 template void func<double>::update_rows(vector<bool> const&);
6136 template void func<double>::update_indices(vector<bool> const&);
6137 // template func<double> min(const func<double>& p1, const func<double>& p2);
6138 // template func<double> max(const func<double>& p1, const func<double>& p2);
6139
6140 }
6141