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