1 // 2 // Created by Hassan on 19/11/2015. 3 // 4 5 #ifndef GRAVITY_CONSTANT_H 6 #define GRAVITY_CONSTANT_H 7 #include <iostream> 8 #include <iomanip> 9 #include <vector> 10 #include <forward_list> 11 #include <assert.h> 12 #include <string> 13 #include <map> 14 #include <complex> 15 #include <memory> 16 #include <typeinfo> 17 #include <typeindex> 18 #include <limits> 19 //#include <gravity/types.h> 20 #include <gravity/utils.h> 21 22 23 using namespace std; 24 25 namespace gravity { 26 27 class param_; 28 class func_; 29 template<typename T=double> class func; 30 /** 31 Transform a scalar to a string with user-specified precision. 32 @param[in] a_value number to be transformed. 33 @param[in] n number of decimals in transformation. 34 @return a string with the specified precision. 35 */ 36 template<class T, class = typename enable_if<is_arithmetic<T>::value>::type> to_string_with_precision(const T a_value,const int n)37 string to_string_with_precision(const T a_value, const int n) 38 { 39 std::ostringstream out; 40 if(std::is_same<T,bool>::value && a_value==numeric_limits<T>::lowest()){ 41 return "0"; 42 } 43 if(std::is_same<T,bool>::value && a_value==numeric_limits<T>::max()){ 44 return "1"; 45 } 46 if(std::numeric_limits<T>::is_specialized && a_value==numeric_limits<T>::max()){ 47 return "+∞"; 48 } 49 if(std::numeric_limits<T>::is_specialized && a_value==numeric_limits<T>::lowest()){ 50 return "−∞"; 51 } 52 if(std::numeric_limits<T>::is_specialized && a_value==numeric_limits<T>::max()){ 53 return "+∞"; 54 } 55 out << std::setprecision(n) << a_value; 56 return out.str(); 57 } 58 /** 59 Transform a complex number to a string with user-specified precision. 60 @param[in] a_value complex number to be transformed. 61 @param[in] n number of decimals in transformation. 62 @return a string with the specified precision. 63 */ 64 string to_string_with_precision(const Cpx& a_value, const int n); 65 66 /** Backbone class for constant */ 67 class constant_{ 68 public: 69 CType _type; /**< Constant type: { binary_c, short_c, integer_c, float_c, double_c, long_c, complex_c, par_c, uexp_c, bexp_c, var_c, func_c}*/ set_type(CType type)70 void set_type(CType type){ _type = type;} 71 72 bool _is_transposed = false; /**< True if the constant is transposed */ 73 bool _is_vector = false; /**< True if the constant is a vector or matrix */ 74 size_t _dim[2] = {1,1}; /*< dimension of current object */ 75 76 bool _polar = false; /**< True in case this is a complex number with a polar representation, rectangular representation if false */ 77 ~constant_()78 virtual ~constant_(){}; get_type()79 CType get_type() const { return _type;} 80 81 82 /** Querries */ is_binary()83 bool is_binary() const{ 84 return (_type==binary_c); 85 }; 86 is_short()87 bool is_short() const{ 88 return (_type==short_c); 89 } 90 is_integer()91 bool is_integer() const{ 92 return (_type==integer_c); 93 }; 94 is_float()95 bool is_float() const{ 96 return (_type==float_c); 97 }; 98 is_double()99 bool is_double() const{ 100 return (_type==double_c); 101 }; 102 is_long()103 bool is_long() const{ 104 return (_type==long_c); 105 }; 106 is_complex()107 bool is_complex() const{ 108 return (_type==complex_c); 109 }; 110 eval_all()111 virtual void eval_all() {}; 112 func_is_number()113 virtual bool func_is_number() const{return is_number();}; 114 is_number()115 virtual bool is_number() const{ 116 return (_type!=par_c && _type!=uexp_c && _type!=bexp_c && _type!=var_c && _type!=func_c); 117 } 118 is_param()119 bool is_param() const{ 120 return (_type==par_c); 121 }; 122 123 /** 124 @return true if current object is a unitary expression. 125 */ is_uexpr()126 bool is_uexpr() const{ 127 return (_type==uexp_c); 128 }; 129 130 /** 131 @return true if current object is a binary expression. 132 */ is_bexpr()133 bool is_bexpr() const{ 134 return (_type==bexp_c); 135 }; 136 137 /** 138 @return true if current object is a unitary or binary expression. 139 */ is_expr()140 bool is_expr() const{ 141 return (_type==uexp_c || _type==bexp_c); 142 }; 143 144 is_var()145 bool is_var() const{ 146 return (_type==var_c); 147 }; 148 is_matrix()149 bool is_matrix() const{ 150 return (_dim[0]>1 && _dim[1]>1); 151 } 152 is_function()153 bool is_function() const{ 154 return (_type==func_c); 155 }; 156 update_double_index()157 virtual void update_double_index() {}; 158 get_all_sign()159 virtual Sign get_all_sign() const {return unknown_;}; 160 virtual Sign get_sign(size_t idx=0) const{return unknown_;}; 161 /** Memory allocation */ allocate_mem()162 virtual void allocate_mem(){};/*<< allocates memory for current and all sub-functions */ 163 164 /** Dimension propagation */ propagate_dim(size_t)165 virtual void propagate_dim(size_t){};/*<< Set dimensions to current and all sub-functions */ is_evaluated()166 virtual bool is_evaluated() const{return false;}; evaluated(bool val)167 virtual void evaluated(bool val){}; reverse_sign()168 virtual void reverse_sign(){};/*<< reverses the sign of current object */ get_dim(size_t i)169 virtual size_t get_dim(size_t i) const { 170 if (i>1) { 171 return _dim[0]; 172 throw invalid_argument("In function: size_t constant_::get_dim(size_t i) const, i is out of range!\n"); 173 } 174 return _dim[i]; 175 } 176 177 /** 178 Returns a copy of the current object, detecting the right class, i.e., param, var, func... 179 @return a shared pointer with a copy of the current object 180 */ copy()181 virtual shared_ptr<constant_> copy() const{return nullptr;}; 182 relax(const map<size_t,shared_ptr<param_>> & vars)183 virtual void relax(const map<size_t, shared_ptr<param_>>& vars){}; print()184 virtual void print(){}; uneval()185 virtual void uneval(){}; to_str()186 virtual string to_str() {return string();}; to_str(int prec)187 virtual string to_str(int prec) {return string();}; to_str(size_t idx,int prec)188 virtual string to_str(size_t idx, int prec) {return string();}; to_str(size_t idx1,size_t idx2,int prec)189 virtual string to_str(size_t idx1, size_t idx2, int prec) {return string();}; 190 191 get_dim()192 virtual size_t get_dim() const { 193 size_t dim = _dim[0]; 194 // if(_dim[1]>0){ 195 dim*= _dim[1]; 196 // } 197 return dim; 198 } 199 vectorize()200 virtual void vectorize(){ 201 _is_vector = true; 202 } 203 vec()204 void vec(){ 205 _is_vector = true; 206 } 207 transpose()208 virtual void transpose(){ 209 _is_transposed = !_is_transposed; 210 _is_vector = true; 211 auto temp = _dim[0]; 212 _dim[0] = _dim[1]; 213 _dim[1] = temp; 214 if(get_dim()==1){ 215 _is_vector = false; 216 } 217 } 218 219 220 221 /** 222 Update the dimensions of current object after it is multiplied with c2. 223 @param[in] c2 object multiplying this. 224 */ update_dot_dim(const constant_ & c2)225 void update_dot_dim(const constant_& c2){ 226 update_dot_dim(*this, c2); 227 } 228 is_row_vector()229 bool is_row_vector() const{ 230 return _dim[0]==1 && _dim[1]>1; 231 } 232 is_column_vector()233 bool is_column_vector() const{ 234 return _dim[1]==1 && _dim[0]>1; 235 } 236 is_scalar()237 bool is_scalar() const{ 238 return !_is_vector && _dim[0]==1 && _dim[1]==1; 239 } 240 241 /** 242 Update the dimensions of current object to correspond to c1.c2. 243 @param[in] c1 first element in product. 244 @param[in] c2 second element in product. 245 */ update_dot_dim(const constant_ & c1,const constant_ & c2)246 void update_dot_dim(const constant_& c1, const constant_& c2){ 247 /* If c2 is a scalar or if multiplying a matrix with a row vector */ 248 if(c2.is_scalar() || (c1.is_matrix() && c2.is_row_vector())){ 249 _dim[0] = c1._dim[0]; 250 _dim[1] = c1._dim[1]; 251 _is_vector = c1._is_vector; 252 _is_transposed = c1._is_transposed; 253 return; 254 } 255 /* If c1 is a scalar or if multiplying a row vector with a matrix */ 256 else if(c1.is_scalar() || (c2.is_matrix() && c1.is_column_vector())){ 257 _dim[0] = c2._dim[0]; 258 _dim[1] = c2._dim[1]; 259 _is_vector = c2._is_vector; 260 _is_transposed = c2._is_transposed; 261 return; 262 } 263 /* Both c1 and c2 are non-scalars */ 264 /* If it is a dot product */ 265 if(c1.is_row_vector() && c2.is_column_vector()){ /* c1^T.c2 */ 266 if(!c1.is_matrix_indexed() && !c2.is_matrix_indexed() && c1._dim[1]!=c2._dim[0]){ 267 throw invalid_argument("Dot product with mismatching dimensions"); 268 } 269 _is_transposed = false;/* The result of a dot product is not transposed */ 270 } 271 else if(c1.is_row_vector() && c2.is_row_vector()){ 272 _is_transposed = true;/* this is a term-wise product of transposed vectors */ 273 } 274 _dim[0] = c1._dim[0]; 275 _dim[1] = c2._dim[1]; 276 if(get_dim()==1){ 277 _is_vector = false; 278 } 279 } 280 281 /** 282 Sets the object dimension to the maximum dimension among all arguments. 283 @param[in] p1 first element in list. 284 @param[in] ps remaining elements in list. 285 */ 286 template<typename... Args> set_max_dim(const constant_ & p1,Args &&...ps)287 void set_max_dim(const constant_& p1, Args&&... ps){ 288 _dim[0] = max(_dim[0], p1._dim[0]); 289 list<constant_*> list = {forward<constant_*>((constant_*)&ps)...}; 290 for(auto &p: list){ 291 _dim[0] = max(_dim[0], p->_dim[0]); 292 } 293 } is_matrix_indexed()294 virtual bool is_matrix_indexed() const{return false;}; is_constant()295 virtual bool is_constant() const{return false;}; is_zero()296 virtual bool is_zero() const{return false;}; /**< Returns true if constant equals 0 */ is_unit()297 virtual bool is_unit() const{return false;}; /**< Returns true if constant equals 1 */ is_neg_unit()298 virtual bool is_neg_unit() const{return false;}; /**< Returns true if constant equals -1 */ is_positive()299 virtual bool is_positive() const{return false;}; /**< Returns true if constant is positive */ is_negative()300 virtual bool is_negative() const{return false;}; /**< Returns true if constant is negative */ is_non_positive()301 virtual bool is_non_positive() const{return false;}; /**< Returns true if constant is non positive */ is_non_negative()302 virtual bool is_non_negative() const{return false;}; /**< Returns true if constant is non negative */ 303 }; 304 305 306 /** Polymorphic class constant, can store an arithmetic or a complex number.*/ 307 template<typename type = double> 308 class constant: public constant_{ 309 public: 310 type _val;/**< value of current constant */ 311 312 template<class T, typename enable_if<is_arithmetic<T>::value>::type* = nullptr> zero()313 T zero(){ 314 return (T)(0); 315 } 316 317 template<class T, typename enable_if<is_same<T, Cpx>::value>::type* = nullptr> zero()318 T zero(){ 319 return Cpx(0,0); 320 } 321 322 is_constant()323 bool is_constant() const{ return true;}; 324 /** Constructors */ update_type()325 void update_type(){ 326 if(typeid(type)==typeid(bool)){ 327 set_type(binary_c); 328 return; 329 } 330 if(typeid(type)==typeid(short)) { 331 set_type(short_c); 332 return; 333 } 334 if(typeid(type)==typeid(int)) { 335 set_type(integer_c); 336 return; 337 } 338 if(typeid(type)==typeid(float)) { 339 set_type(float_c); 340 return; 341 } 342 if(typeid(type)==typeid(double)) { 343 set_type(double_c); 344 return; 345 } 346 if(typeid(type)==typeid(long double)) { 347 set_type(long_c); 348 return; 349 } 350 if(typeid(type)==typeid(complex<double>)) { 351 set_type(complex_c); 352 return; 353 } 354 throw invalid_argument("Unknown constant type."); 355 } 356 constant()357 constant(){ 358 _val = zero<type>(); 359 update_type(); 360 } 361 362 ~constant()363 ~constant(){}; 364 365 366 constant& operator=(const constant& c) { 367 _type = c._type; 368 _is_transposed = c._is_transposed; 369 _is_vector = c._is_vector; 370 _val = c._val; 371 return *this; 372 } 373 374 template<class T2, typename std::enable_if<is_convertible<T2, type>::value && sizeof(T2) < sizeof(type)>::type* = nullptr> 375 constant& operator=(const constant<T2>& c) { 376 update_type(); 377 _is_transposed = c._is_transposed; 378 _is_vector = c._is_vector; 379 _val = c._val; 380 return *this; 381 } 382 383 384 template<class T2, typename std::enable_if<is_convertible<T2, type>::value && sizeof(T2) < sizeof(type)>::type* = nullptr> 385 constant(const constant<T2>& c){ /**< Copy constructor */ 386 *this = c; 387 }; 388 constant(const constant & c)389 constant(const constant& c){ /**< Copy constructor */ 390 *this = c; 391 }; 392 copy()393 shared_ptr<constant_> copy() const{return make_shared<constant>(*this);}; 394 constant(const type & val)395 constant(const type& val):constant(){ 396 _val = val; 397 }; 398 range()399 shared_ptr<pair<type,type>> range() const{ 400 return make_shared<pair<type,type>>(_val,_val); 401 } 402 tr()403 constant tr() const{ 404 auto newc(*this); 405 newc.transpose(); 406 return newc; 407 }; 408 409 eval()410 type eval() const { return _val;} 411 set_val(type val)412 void set_val(type val) { 413 _val = val; 414 } 415 reverse_sign()416 void reverse_sign(){ 417 _val *= -1; 418 } 419 get_all_sign()420 Sign get_all_sign() const {return get_sign();}; 421 422 Sign get_sign(size_t idx = 0) const{ 423 return get_sign_(); 424 } 425 get_sign_()426 template<class T=type, class = typename enable_if<is_same<T, Cpx>::value>::type> Sign get_sign_() const{ 427 if (_val == Cpx(0,0)) { 428 return zero_; 429 } 430 if ((_val.real() < 0 && _val.imag() < 0)) { 431 return neg_; 432 } 433 if ((_val.real() > 0 && _val.imag() > 0)) { 434 return pos_; 435 } 436 if (_val.real() >= 0 && _val.imag() >= 0) { 437 return non_neg_; 438 } 439 if (_val.real() <= 0 && _val.imag() <= 0) { 440 return non_pos_; 441 } 442 return unknown_; 443 } 444 445 template<typename T=type, get_sign_()446 typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> Sign get_sign_() const{ 447 if (_val==0) { 448 return zero_; 449 } 450 if (_val > 0) { 451 return pos_; 452 } 453 if (_val < 0) { 454 return neg_; 455 } 456 return unknown_; 457 } 458 459 460 461 /** Operators */ 462 is_zero()463 bool is_zero() const { return zero_val();}; 464 zero_val()465 template<class T=type, class = typename enable_if<is_same<T, Cpx>::value>::type> bool zero_val() const{ 466 return (_val == Cpx(0,0)); 467 } 468 469 template<typename T=type, zero_val()470 typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> bool zero_val() const{ 471 return (_val == 0); 472 } 473 474 is_unit()475 bool is_unit() const{ 476 return unit_val(); 477 } 478 unit_val()479 template<class T=type, class = typename enable_if<is_same<T, Cpx>::value>::type> bool unit_val() const{ 480 return (!_is_vector && !_is_transposed && _val == Cpx(1,0)); 481 } 482 483 template<typename T=type, unit_val()484 typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> bool unit_val() const{ 485 return (!_is_vector && !_is_transposed && _val == 1); 486 } 487 is_negative()488 bool is_negative() const { 489 return get_sign()==neg_; 490 } 491 is_non_negative()492 bool is_non_negative() const { 493 return (get_sign()==zero_||get_sign()==pos_||get_sign()==non_neg_); 494 } 495 is_non_positive()496 bool is_non_positive() const { 497 return (get_sign()==zero_||get_sign()==neg_||get_sign()==non_pos_); 498 } 499 is_positive()500 bool is_positive() const { 501 return get_sign()==pos_; 502 } 503 504 bool operator==(const constant& c) const { 505 return (_type==c._type && _val==c._val); 506 } 507 508 bool operator==(const type& v) const{ 509 return _val==v; 510 } 511 512 constant& operator=(const type& val){ 513 _val = val; 514 return *this; 515 } 516 517 constant& operator+=(const type& v){ 518 _val += v; 519 return *this; 520 } 521 522 constant& operator-=(const type& v){ 523 _val -= v; 524 return *this; 525 } 526 527 constant& operator*=(const type& v){ 528 _val *= v; 529 return *this; 530 } 531 532 constant& operator/=(const type& v){ 533 _val /= v; 534 return *this; 535 } 536 537 friend constant operator+(const constant& c1, const constant& c2){ 538 if(c1._is_vector) 539 return constant(c1) += c2._val; 540 return constant(c2) += c1._val; 541 } 542 543 friend constant operator-(const constant& c1, const constant& c2){ 544 if(c2._is_vector){ 545 return constant(-1*c2) += c1._val; 546 } 547 return constant(c1) -= c2._val; 548 } 549 550 friend constant operator/(const constant& c1, const constant& c2){ 551 if(c2._is_vector){ 552 return constant(1/c2) *= c1._val; 553 } 554 return constant(c1) /= c2._val; 555 } 556 557 friend constant operator*(const constant& c1, const constant& c2){ 558 if(c2._is_vector){ 559 return constant(c2) *= c1._val; 560 } 561 return constant(c1) *= c2._val; 562 } 563 564 friend constant operator^(const constant& c1, const constant& c2){ 565 return constant(pow(c1._val,c2._val)); 566 } 567 568 569 friend constant operator+(const constant& c, type cst){ 570 return constant(c) += cst; 571 } 572 573 friend constant operator-(const constant& c, type cst){ 574 return constant(c) -= cst; 575 } 576 577 friend constant operator*(const constant& c, type cst){ 578 return constant(c) *= cst; 579 } 580 581 582 friend constant operator/(const constant& c, type cst){ 583 return constant(c) /= cst; 584 } 585 586 friend constant operator+(type cst, const constant& c){ 587 return constant(c) += cst; 588 } 589 590 friend constant operator-(type cst, const constant& c){ 591 return constant(-1*c) += cst; 592 } 593 594 friend constant operator*(type cst, const constant& c){ 595 return constant(c) *= cst; 596 } 597 598 599 friend constant operator/(type cst, const constant& c){ 600 return constant(c) *= cst/std::pow(c._val,2); 601 } 602 cos(const constant & c)603 friend constant cos(const constant& c){ 604 return constant(cos(c._val)); 605 } 606 sin(const constant & c)607 friend constant sin(const constant& c){ 608 return constant(sin(c._val)); 609 } 610 sqrt(const constant & c)611 friend constant sqrt(const constant& c){ 612 return constant(sqrt(c._val)); 613 } 614 expo(const constant & c)615 friend constant expo(const constant& c){ 616 return constant(exp(c._val)); 617 } 618 log(const constant & c)619 friend constant log(const constant& c){ 620 return constant(log(c._val)); 621 } 622 623 624 /** Output */ print()625 void print(){ 626 cout << to_str(10); 627 } 628 print(int prec)629 void print(int prec) { 630 cout << to_str(prec); 631 } 632 633 void println(int prec = 10) { 634 cout << to_str(prec) << endl; 635 } 636 to_str()637 string to_str() { 638 return to_string_with_precision(_val,5); 639 } 640 641 string to_str(int prec = 10) { 642 return to_string_with_precision(_val,prec); 643 } 644 645 string to_str(size_t index, int prec = 10) { 646 return to_string_with_precision(_val,prec); 647 } 648 649 650 }; 651 /** 652 Returns the conjugate of cst. 653 @param[in] cst complex number. 654 @return the conjugate of cst. 655 */ 656 constant<Cpx> conj(const constant<Cpx>& cst); 657 constant<double> real(const constant<Cpx>& cst); 658 constant<double> imag(const constant<Cpx>& cst); 659 /** 660 Returns the square magnitude of cst. 661 @param[in] cst complex number. 662 @return the square magnitude of cst. 663 */ 664 constant<double> sqrmag(const constant<Cpx>& cst); 665 constant<double> angle(const constant<Cpx>& cst); 666 667 template<class T, typename enable_if<is_arithmetic<T>::value>::type* = nullptr> unit()668 constant<T> unit(){ 669 return constant<T>(1); 670 } 671 672 template<class T, typename enable_if<is_same<T, Cpx>::value>::type* = nullptr> unit()673 constant<T> unit(){ 674 return constant<T>(Cpx(1,0)); 675 } 676 677 template<class T, typename enable_if<is_arithmetic<T>::value>::type* = nullptr> zero()678 constant<T> zero(){ 679 return constant<T>(0); 680 } 681 682 template<class T, typename enable_if<is_same<T, Cpx>::value>::type* = nullptr> zero()683 constant<T> zero(){ 684 return constant<T>(Cpx(0,0)); 685 } 686 687 } 688 689 690 691 #endif //GRAVITY_CONSTANT_H 692