fmatvec 0.0.0

ast.h
1 #ifndef _FMATVEC_AST_H_
2 #define _FMATVEC_AST_H_
3
4 #include "fmatvec/function.h"
5 #include <map>
6 #include <array>
7 #include <memory>
8 #include <boost/uuid/uuid.hpp>
9 #include <boost/uuid/uuid_generators.hpp>
10 #include <boost/container/small_vector.hpp>
11 #include <utility>
12 #include <fmatvec/stream.h>
13
14 namespace fmatvec {
15
16 // forward declarations
17
18 class IndependentVariable;
19
20 namespace AST {
21 class Vertex;
22 class Symbol;
23 class Operation;
24 class NativeFunction;
25 template<class T> class Constant;
26 FMATVEC_EXPORT SymbolicExpression substScalar(const SymbolicExpression &se, const IndependentVariable& a, const SymbolicExpression &b);
27 }
28
29 template<class... Arg> class Eval;
30 class EvalHelper;
31
37 class FMATVEC_EXPORT SymbolicExpression : public std::shared_ptr<const AST::Vertex> {
38 template<class... Arg> friend class Eval;
39 friend class EvalHelper;
40 friend class AST::Operation;
41 friend class AST::Constant<long>;
42 friend class AST::Constant<double>;
43 friend class AST::NativeFunction;
44 friend FMATVEC_EXPORT SymbolicExpression parDer(const SymbolicExpression &dep, const IndependentVariable &indep);
45 friend SymbolicExpression AST::substScalar(const SymbolicExpression &se,
46 const IndependentVariable& a, const SymbolicExpression &b);
47 protected:
48 template<class T> SymbolicExpression(const shared_ptr<T> &x);
49 #ifndef SWIG
50 static const struct ConstructSymbol{} constructSymbol; // just used for tag dispatching
51 #endif
52 SymbolicExpression(ConstructSymbol);
53
54 // do not allow to call these functions inherited from std::shared_ptr
55 using std::shared_ptr<const AST::Vertex>::get;
56 using std::shared_ptr<const AST::Vertex>::operator*;
57 using std::shared_ptr<const AST::Vertex>::operator->;
58 public:
60 SymbolicExpression();
62 SymbolicExpression(const std::string &str);
64 SymbolicExpression(int x);
65 SymbolicExpression(long x);
67 SymbolicExpression(double x);
68
72 static void garbageCollect();
73
74 // default ctors and assignment operators
75 SymbolicExpression(const SymbolicExpression& x) = default;
76 SymbolicExpression(SymbolicExpression&& x) = default;
77 SymbolicExpression& operator=(const SymbolicExpression& x) = default;
78 SymbolicExpression& operator=(SymbolicExpression&& x) = default;
79 ~SymbolicExpression() = default;
80
81 // operators
82 #ifndef SWIG
83 SymbolicExpression operator+(const SymbolicExpression &b) const;
84 SymbolicExpression operator-(const SymbolicExpression &b) const;
85 SymbolicExpression operator*(const SymbolicExpression &b) const;
86 SymbolicExpression operator/(const SymbolicExpression &b) const;
87 #endif
88 SymbolicExpression operator-() const;
89 SymbolicExpression& operator+=(const SymbolicExpression &b);
90 SymbolicExpression& operator-=(const SymbolicExpression &b);
91 SymbolicExpression& operator*=(const SymbolicExpression &b);
92 SymbolicExpression& operator/=(const SymbolicExpression &b);
93
94 // We implement operator<<= to be able to use SymbolicExpression fully equivalent to the vector/matrix classes
95 // which have a assign method which redimension before assigning. Redimensioning is not required for a scalar
96 // type like assign
97 SymbolicExpression& operator<<=(const SymbolicExpression &src);
98
99 #if defined(FMATVEC_DEBUG) && !defined(SWIG)
100 static signed long evalOperationsCount;
101 #endif
102 };
103
107 class FMATVEC_EXPORT IndependentVariable : public SymbolicExpression {
108 friend class AST::Symbol;
109 friend FMATVEC_EXPORT std::istream& operator>>(std::istream& s, IndependentVariable &v);
110 public:
112 IndependentVariable();
114 IndependentVariable(const std::string &str);
115
116 #ifndef SWIG
117 inline IndependentVariable& operator^=(double x);
119 #endif
120
121 private:
122 IndependentVariable(const shared_ptr<const AST::Symbol> &x);
123 };
124
125 // define the operator results of SymbolicExpression and IndependentVariable with other AT types.
126 FMATVEC_OPERATORRESULT1(SymbolicExpression, SymbolicExpression)
127
128 FMATVEC_OPERATORRESULT2(SymbolicExpression, int, SymbolicExpression)
129 FMATVEC_OPERATORRESULT2(SymbolicExpression, long, SymbolicExpression)
130 FMATVEC_OPERATORRESULT2(SymbolicExpression, double, SymbolicExpression)
131
132 FMATVEC_OPERATORRESULT1(IndependentVariable, SymbolicExpression)
133
134 FMATVEC_OPERATORRESULT2(IndependentVariable, int, SymbolicExpression)
135 FMATVEC_OPERATORRESULT2(IndependentVariable, long, SymbolicExpression)
136 FMATVEC_OPERATORRESULT2(IndependentVariable, double, SymbolicExpression)
137 FMATVEC_OPERATORRESULT2(IndependentVariable, SymbolicExpression, SymbolicExpression)
138
139 // Some member function definition of SymbolicExpression/IndependentVariable are moved to the end of this file
140 // since they need the defintion of the other class defined in this file.
141
142 FMATVEC_EXPORT SymbolicExpression operator+(double a, const SymbolicExpression &b);
143 FMATVEC_EXPORT SymbolicExpression operator-(double a, const SymbolicExpression &b);
144 FMATVEC_EXPORT SymbolicExpression operator*(double a, const SymbolicExpression &b);
145 FMATVEC_EXPORT SymbolicExpression operator/(double a, const SymbolicExpression &b);
146 FMATVEC_EXPORT SymbolicExpression operator+(int a, const SymbolicExpression &b);
147 FMATVEC_EXPORT SymbolicExpression operator-(int a, const SymbolicExpression &b);
148 FMATVEC_EXPORT SymbolicExpression operator*(int a, const SymbolicExpression &b);
149 FMATVEC_EXPORT SymbolicExpression operator/(int a, const SymbolicExpression &b);
150 FMATVEC_EXPORT SymbolicExpression operator+(long a, const SymbolicExpression &b);
151 FMATVEC_EXPORT SymbolicExpression operator-(long a, const SymbolicExpression &b);
152 FMATVEC_EXPORT SymbolicExpression operator*(long a, const SymbolicExpression &b);
153 FMATVEC_EXPORT SymbolicExpression operator/(long a, const SymbolicExpression &b);
154
157 FMATVEC_EXPORT SymbolicExpression parDer(const SymbolicExpression &dep, const IndependentVariable &indep);
158
160 FMATVEC_EXPORT std::ostream& operator<<(std::ostream& s, const SymbolicExpression& se);
162 FMATVEC_EXPORT std::istream& operator>>(std::istream& s, SymbolicExpression &se);
164 FMATVEC_EXPORT std::istream& operator>>(std::istream& s, IndependentVariable &v);
165
166 // function operations overloaded for SymbolicExpression
167 FMATVEC_EXPORT SymbolicExpression pow(const SymbolicExpression &a, const SymbolicExpression &b);
168 FMATVEC_EXPORT SymbolicExpression log(const SymbolicExpression &a);
169 FMATVEC_EXPORT SymbolicExpression sqrt(const SymbolicExpression &a);
170 FMATVEC_EXPORT SymbolicExpression sin(const SymbolicExpression &a);
171 FMATVEC_EXPORT SymbolicExpression cos(const SymbolicExpression &a);
172 FMATVEC_EXPORT SymbolicExpression tan(const SymbolicExpression &a);
173 FMATVEC_EXPORT SymbolicExpression sinh(const SymbolicExpression &a);
174 FMATVEC_EXPORT SymbolicExpression cosh(const SymbolicExpression &a);
175 FMATVEC_EXPORT SymbolicExpression tanh(const SymbolicExpression &a);
176 FMATVEC_EXPORT SymbolicExpression asin(const SymbolicExpression &a);
177 FMATVEC_EXPORT SymbolicExpression acos(const SymbolicExpression &a);
178 FMATVEC_EXPORT SymbolicExpression atan(const SymbolicExpression &a);
179 FMATVEC_EXPORT SymbolicExpression atan2(const SymbolicExpression &y, const SymbolicExpression &x);
180 FMATVEC_EXPORT SymbolicExpression asinh(const SymbolicExpression &a);
181 FMATVEC_EXPORT SymbolicExpression acosh(const SymbolicExpression &a);
182 FMATVEC_EXPORT SymbolicExpression atanh(const SymbolicExpression &a);
183 FMATVEC_EXPORT SymbolicExpression exp(const SymbolicExpression &a);
184 FMATVEC_EXPORT SymbolicExpression sign(const SymbolicExpression &a);
185 FMATVEC_EXPORT SymbolicExpression heaviside(const SymbolicExpression &a);
186 FMATVEC_EXPORT SymbolicExpression abs(const SymbolicExpression &a);
187 FMATVEC_EXPORT SymbolicExpression min(const SymbolicExpression &a, const SymbolicExpression &b);
188 FMATVEC_EXPORT SymbolicExpression max(const SymbolicExpression &a, const SymbolicExpression &b);
189 FMATVEC_EXPORT SymbolicExpression condition(const SymbolicExpression &c, const SymbolicExpression &gt, const SymbolicExpression &lt);
190
191 #ifndef SWIG
192 namespace AST { // internal namespace
193
194 /* ***** Struct for a "bytecode" for fast runtime evaluation *****
195 * This class uses very low-level technices like raw-pointers to gain maximal performance
196 * during runtime evaluation.
197 *
198 * This class store a kind of "bytecode" which can be used to evaluate an mathematical expression very fast.
199 * This class is not copyable and not moveable since it uses intern pointers between instances of this class.
200 * For optimal performance this class should be used in continous memory e.g. as std::vector<ByteCode>.
201 * For the same reason, to keep all data near together for best processor memory cache usage, all members
202 * of these class also to not use allocated memory: std::array instead of std::vector is used.
203 * The usage of std::vector<ByteCode> requires that the move-ctor must be provided at compile time
204 * (it cannot be delete'd) but the move-ctor throws at runtime when called to ensure its not used.
205 * Hence std::vector<ByteCode> can only be used without reallocation: std::vector::reserve must be called before use.
206 * To add an element only std::vector::emplace_back() can be used.
207 */
208 struct FMATVEC_EXPORT ByteCode {
209 static constexpr size_t N { 12 };
210 ByteCode(size_t n);
211 ByteCode(const ByteCode &) = delete;
212 ByteCode(ByteCode &&) noexcept;
213 ByteCode &operator=(const ByteCode &) = delete;
214 ByteCode &operator=(ByteCode &&) = delete;
215 ~ByteCode() = default;
216 using Arg = boost::container::small_vector<double*,N>;
217 std::function<void(double*, const Arg&)> func; // the operation this byteCode entry executes
218 double retValue; // storage of the return value of the operation: retPtr may point to this value
219 double* retPtr; // a pointer to which the operation can write its result to
220 Arg argsPtr; // pointers from which the operation reads its arguments
221 };
222
223 template<class Func, class ArgS>
224 struct SymbolicFuncWrapArg1 {
225 static SymbolicExpression call(
226 const std::shared_ptr<Function<Func>> &func,
227 const ArgS &arg);
228 };
229
230 template<>
231 struct SymbolicFuncWrapArg1<double(double), SymbolicExpression> {
232 static SymbolicExpression call(
233 const std::shared_ptr<Function<double(double)>> &func,
234 const SymbolicExpression &arg);
235 };
236
237 template<class Type>
238 struct SymbolicFuncWrapArg1<double(Vector<Type, double>), Vector<Type, SymbolicExpression>> {
239 static SymbolicExpression call(
240 const std::shared_ptr<Function<double(Vector<Type, double>)>> &func,
241 const Vector<Type, SymbolicExpression> &arg);
242 };
243
244 template<class Func, class Arg1S, class Arg2S>
245 struct SymbolicFuncWrapArg2 {
246 static SymbolicExpression call(
247 const std::shared_ptr<Function<Func>> &func,
248 const Arg1S &arg1, const Arg2S &arg2);
249 };
250
251 template<>
252 struct SymbolicFuncWrapArg2<double(double,double), SymbolicExpression, SymbolicExpression> {
253 static SymbolicExpression call(
254 const std::shared_ptr<Function<double(double,double)>> &func,
255 const SymbolicExpression &arg1, const SymbolicExpression &arg2);
256 };
257
258 template<class Type>
259 struct SymbolicFuncWrapArg2<double(Vector<Type,double>,double), Vector<Type, SymbolicExpression>, SymbolicExpression> {
260 static SymbolicExpression call(
261 const std::shared_ptr<Function<double(Vector<Type,double>,double)>> &func,
262 const Vector<Type, SymbolicExpression> &arg1, const SymbolicExpression &arg2);
263 };
264
265 template<class Type>
266 struct SymbolicFuncWrapArg2<double(double,Vector<Type,double>), SymbolicExpression, Vector<Type, SymbolicExpression>> {
267 static SymbolicExpression call(
268 const std::shared_ptr<Function<double(double,Vector<Type,double>)>> &func,
269 const SymbolicExpression &arg1, const Vector<Type, SymbolicExpression> &arg2);
270 };
271
272 template<class Type1, class Type2>
273 struct SymbolicFuncWrapArg2<double(Vector<Type1,double>,Vector<Type2,double>), Vector<Type1, SymbolicExpression>, Vector<Type2, SymbolicExpression>> {
274 static SymbolicExpression call(
275 const std::shared_ptr<Function<double(Vector<Type1,double>,Vector<Type2,double>)>> &func,
276 const Vector<Type1, SymbolicExpression> &arg1, const Vector<Type2, SymbolicExpression> &arg2);
277 };
278
279 // ***** Vertex *****
280
282 class FMATVEC_EXPORT Vertex {
283 friend Operation;
284 friend NativeFunction;
285 public:
286
288 virtual SymbolicExpression parDer(const IndependentVariable &x) const=0;
290 inline virtual bool isConstantInt() const;
293 bool isZero() const;
296 bool isOne() const;
297
298 virtual std::vector<ByteCode>::iterator dumpByteCode(std::vector<ByteCode> &byteCode,
299 std::map<const Vertex*, std::vector<AST::ByteCode>::iterator> &existingVertex) const=0;
300
301 virtual void walkVertex(const std::function<void(const std::shared_ptr<const Vertex>&)> &func) const=0;
302
303 protected:
304
305 // helper function to make it easy to implement new expression optimizations. See ast.cc Operation::create for details.
306 // Returns true if this Vertex the Vertex (SymbolicExpression) b. Every Symbol variables in this Vertex are free. This
307 // means that true is also returned if the Symbols in this Vertex can be replaced by anything such that the expressions
308 // are equal. All these required replacements are stored in m.
309 // NOTE:
310 // We define our own less operator for std::map<IndependentVariable, ...> since we abort in the overloaded function
311 // inline bool operator<(const fmatvec::SymbolicExpression& a, const fmatvec::SymbolicExpression& b) noexcept,
312 // see below. This abort is needed to avoid that e.g. std::min/std::max can be called with SymbolicExpression arguments.
313 // If your code calls this abort function you may call somewhere std::min(SymbolicExpression, SymbolicExpression) instead of
314 // fmatvec::min(SymbolicExpression, SymbolicExpression). The same applies for std::max/fmatvec::max.
315 // I was not able to solve this in a better way!
316 struct LessIV {
317 bool operator()(const IndependentVariable& a, const IndependentVariable& b) const {
318 return a.std::shared_ptr<const AST::Vertex>::get() < b.std::shared_ptr<const AST::Vertex>::get();
319 }
320 };
321 using MapIVSE = std::map<IndependentVariable, SymbolicExpression, LessIV>;
322 virtual bool equal(const SymbolicExpression &b, MapIVSE &m) const=0;
323 };
324
325 inline bool Vertex::isConstantInt() const {
326 return false;
327 }
328
329 // ***** Constant *****
330
332 template<class T>
333 class FMATVEC_EXPORT Constant : public Vertex, public std::enable_shared_from_this<Constant<T>> {
334 friend SymbolicExpression;
335 public:
336
337 static SymbolicExpression create(const T& c_);
338 SymbolicExpression parDer(const IndependentVariable &x) const override;
339 inline bool isConstantInt() const override;
341 inline const T& getValue() const;
342 std::vector<ByteCode>::iterator dumpByteCode(std::vector<ByteCode> &byteCode,
343 std::map<const Vertex*, std::vector<AST::ByteCode>::iterator> &existingVertex) const override;
344
345 void walkVertex(const std::function<void(const std::shared_ptr<const Vertex>&)> &func) const override;
346
347 private:
348
349 Constant(const T& c_);
350 const T c;
351 bool equal(const SymbolicExpression &b, MapIVSE &m) const override;
352 using CacheKey = T;
353 static std::map<CacheKey, std::weak_ptr<const Constant>> cache;
354 };
355
356 template<>
357 inline bool Constant<double>::isConstantInt() const {
358 return false;
359 }
360
361 template<>
362 inline bool Constant<long>::isConstantInt() const {
363 return true;
364 }
365
366 template<class T>
367 const T& Constant<T>::getValue() const {
368 return c;
369 }
370
371 // ***** Symbol *****
372
374 class FMATVEC_EXPORT Symbol : public Vertex, public std::enable_shared_from_this<Symbol> {
375 friend SymbolicExpression;
376 friend IndependentVariable;
377 public:
378
379 static IndependentVariable create(const boost::uuids::uuid& uuid_=boost::uuids::random_generator()());
380 SymbolicExpression parDer(const IndependentVariable &x) const override;
383 inline void setValue(double x_) const;
384
385 std::string getUUIDStr() const;
386
387 std::vector<ByteCode>::iterator dumpByteCode(std::vector<ByteCode> &byteCode,
388 std::map<const Vertex*, std::vector<AST::ByteCode>::iterator> &existingVertex) const override;
389
390 void walkVertex(const std::function<void(const std::shared_ptr<const Vertex>&)> &func) const override;
391
392 private:
393
394 Symbol(const boost::uuids::uuid& uuid_);
395 bool equal(const SymbolicExpression &b, MapIVSE &m) const override;
396 mutable double x = 0.0;
397 boost::uuids::uuid uuid; // each variable has a uuid (this is only used when the AST is serialized and for caching)
398 using CacheKey = boost::uuids::uuid;
399 static std::map<CacheKey, std::weak_ptr<const Symbol>> cache;
400 };
401
402 void Symbol::setValue(double x_) const {
403 x=x_;
404 }
405
406 // ***** ScalarFunctionWrapArg *****
407
408 class ScalarFunctionWrapArg {
409 public:
410 virtual double operator()(const ByteCode::Arg &arg) = 0;
411 virtual double dirDer(const ByteCode::Arg &arg) = 0;
412 virtual double dirDerDirDer(const ByteCode::Arg &arg) = 0;
413 };
414
415 class ScalarFunctionWrapArgS : public ScalarFunctionWrapArg {
416 public:
417 ScalarFunctionWrapArgS(std::shared_ptr<fmatvec::Function<double(double)>> func_) : func(std::move(func_)) {}
418 double operator()(const ByteCode::Arg &arg) override {
419 return (*func)(*arg[0]);
420 }
421 double dirDer(const ByteCode::Arg &arg) override {
422 return func->dirDer(*arg[1], *arg[0]);
423 }
424 double dirDerDirDer(const ByteCode::Arg &arg) override {
425 return func->dirDerDirDer(*arg[1], *arg[2], *arg[0]);
426 }
427 private:
428 std::shared_ptr<fmatvec::Function<double(double)>> func;
429 };
430
431 template<class Type>
432 class ScalarFunctionWrapArgV : public ScalarFunctionWrapArg {
433 public:
434 ScalarFunctionWrapArgV(const std::shared_ptr<fmatvec::Function<double(fmatvec::Vector<Type,double>)>> &func_, int argSize) : func(func_) {
435 argV.resize(argSize);
436 dir1V.resize(argSize);
437 dir2V.resize(argSize);
438 }
439 double operator()(const ByteCode::Arg &arg) override {
440 for(size_t i=0; i<argV.size(); ++i)
441 argV(i) = *arg[i];
442 return (*func)(argV);
443 }
444 double dirDer(const ByteCode::Arg &arg) override {
445 size_t argVsize=argV.size();
446 for(size_t i=0; i<argVsize; ++i)
447 argV(i) = *arg[i];
448 for(size_t i=0; i<argVsize; ++i)
449 dir1V(i) = *arg[i+argVsize];
450 return func->dirDer(dir1V, argV);
451 }
452 double dirDerDirDer(const ByteCode::Arg &arg) override {
453 size_t argVsize=argV.size();
454 for(size_t i=0; i<argVsize; ++i)
455 argV(i) = *arg[i];
456 for(size_t i=0; i<argVsize; ++i)
457 dir1V(i) = *arg[i+argVsize];
458 for(size_t i=0; i<argVsize; ++i)
459 dir2V(i) = *arg[i+2*argVsize];
460 return func->dirDerDirDer(dir1V, dir2V, argV);
461 }
462 private:
463 std::shared_ptr<fmatvec::Function<double(fmatvec::Vector<Type,double>)>> func;
464 fmatvec::Vector<Type,double> argV;
465 fmatvec::Vector<Type,double> dir1V;
466 fmatvec::Vector<Type,double> dir2V;
467 };
468
469 class ScalarFunctionWrapArgSS : public ScalarFunctionWrapArg {
470 public:
471 ScalarFunctionWrapArgSS(std::shared_ptr<fmatvec::Function<double(double,double)>> func_) : func(std::move(func_)) {}
472 double operator()(const ByteCode::Arg &arg) override {
473 return (*func)(*arg[0], *arg[1]);
474 }
475 double dirDer(const ByteCode::Arg &arg) override {
476 return func->dirDer1(*arg[2], *arg[0], *arg[1]) +
477 func->dirDer2(*arg[3], *arg[0], *arg[1]);
478 }
479 double dirDerDirDer(const ByteCode::Arg &arg) override {
480 return func->dirDer1DirDer1(*arg[2], *arg[4], *arg[0], *arg[1]) +
481 func->dirDer2DirDer1(*arg[5], *arg[2], *arg[0], *arg[1]) +
482 func->dirDer2DirDer1(*arg[3], *arg[4], *arg[0], *arg[1]) +
483 func->dirDer2DirDer2(*arg[3], *arg[5], *arg[0], *arg[1]);
484 }
485 private:
486 std::shared_ptr<fmatvec::Function<double(double,double)>> func;
487 };
488
489 template<class Type>
490 class ScalarFunctionWrapArgVS : public ScalarFunctionWrapArg {
491 public:
492 ScalarFunctionWrapArgVS(const std::shared_ptr<fmatvec::Function<double(Vector<Type,double>,double)>> &func_, int argSize) : func(func_) {
493 argV.resize(argSize);
494 dir1V.resize(argSize);
495 dir2V.resize(argSize);
496 }
497 double operator()(const ByteCode::Arg &arg) override {
498 for(size_t i=0; i<argV.size(); ++i)
499 argV(i) = *arg[i];
500 return (*func)(argV, *arg[argV.size()]);
501 }
502 double dirDer(const ByteCode::Arg &arg) override {
503 auto argVsize=argV.size();
504 for(size_t i=0; i<argVsize; ++i)
505 argV(i) = *arg[i];
506 for(size_t i=0; i<argVsize; ++i)
507 dir1V(i) = *arg[i+1+argVsize];
508 return func->dirDer1(dir1V , argV, *arg[argVsize]) +
509 func->dirDer2(*arg[2*argVsize+1], argV, *arg[argVsize]);
510 }
511 double dirDerDirDer(const ByteCode::Arg &arg) override {
512 auto argVsize=argV.size();
513 for(size_t i=0; i<argVsize; ++i)
514 argV(i) = *arg[i];
515 for(size_t i=0; i<argVsize; ++i)
516 dir1V(i) = *arg[i+1+argVsize];
517 for(size_t i=0; i<argVsize; ++i)
518 dir2V(i) = *arg[i+2+2*argVsize];
519 return func->dirDer1DirDer1(dir1V, dir2V, argV, *arg[argVsize]) +
520 func->dirDer2DirDer1(*arg[1+2*argVsize], dir2V , argV, *arg[argVsize]) +
521 func->dirDer2DirDer1(*arg[2+3*argVsize], dir1V , argV, *arg[argVsize]) +
522 func->dirDer2DirDer2(*arg[1+2*argVsize], *arg[2+3*argVsize], argV, *arg[argVsize]);
523 }
524 private:
525 std::shared_ptr<fmatvec::Function<double(Vector<Type,double>,double)>> func;
526 fmatvec::Vector<Type,double> argV;
527 fmatvec::Vector<Type,double> dir1V;
528 fmatvec::Vector<Type,double> dir2V;
529 };
530
531 template<class Type>
532 class ScalarFunctionWrapArgSV : public ScalarFunctionWrapArg {
533 public:
534 ScalarFunctionWrapArgSV(const std::shared_ptr<fmatvec::Function<double(double,Vector<Type,double>)>> &func_, int argSize) : func(func_) {
535 argV.resize(argSize);
536 dir1V.resize(argSize);
537 dir2V.resize(argSize);
538 }
539 double operator()(const ByteCode::Arg &arg) override {
540 for(size_t i=0; i<argV.size(); ++i)
541 argV(i) = *arg[i+1];
542 return (*func)(*arg[0], argV);
543 }
544 double dirDer(const ByteCode::Arg &arg) override {
545 auto argVsize=argV.size();
546 for(size_t i=0; i<argVsize; ++i)
547 argV(i) = *arg[i+1];
548 for(size_t i=0; i<argVsize; ++i)
549 dir1V(i) = *arg[i+2+argVsize];
550 return func->dirDer1(*arg[argVsize+1], *arg[0], argV) +
551 func->dirDer2(dir1V , *arg[0], argV);
552 }
553 double dirDerDirDer(const ByteCode::Arg &arg) override {
554 auto argVsize=argV.size();
555 for(size_t i=0; i<argVsize; ++i)
556 argV(i) = *arg[i+1];
557 for(size_t i=0; i<argVsize; ++i)
558 dir1V(i) = *arg[i+2+argVsize];
559 for(size_t i=0; i<argVsize; ++i)
560 dir2V(i) = *arg[i+3+2*argVsize];
561 return func->dirDer1DirDer1(*arg[argVsize+1], *arg[2*argVsize+2], *arg[0], argV) +
562 func->dirDer2DirDer1(dir2V , *arg[argVsize+1] , *arg[0], argV) +
563 func->dirDer2DirDer1(dir1V , *arg[2*argVsize+2], *arg[0], argV) +
564 func->dirDer2DirDer2(dir1V , dir2V , *arg[0], argV);
565 }
566 private:
567 std::shared_ptr<fmatvec::Function<double(double,Vector<Type,double>)>> func;
568 fmatvec::Vector<Type,double> argV;
569 fmatvec::Vector<Type,double> dir1V;
570 fmatvec::Vector<Type,double> dir2V;
571 };
572
573 template<class Type1, class Type2>
574 class ScalarFunctionWrapArgVV : public ScalarFunctionWrapArg {
575 public:
576 ScalarFunctionWrapArgVV(const std::shared_ptr<fmatvec::Function<double(Vector<Type1,double>,Vector<Type2,double>)>> &func_,
577 int argaSize, int argbSize) : func(func_) {
578 argV1.resize(argaSize);
579 argV2.resize(argbSize);
580 dir1V1.resize(argaSize);
581 dir1V2.resize(argbSize);
582 dir2V1.resize(argaSize);
583 dir2V2.resize(argbSize);
584 }
585 double operator()(const ByteCode::Arg &arg) override {
586 auto argV1size = argV1.size();
587 auto argV2size = argV2.size();
588 for(int i=0; i<argV1size; ++i)
589 argV1(i)=*arg[i];
590 for(int i=0; i<argV2size; ++i)
591 argV2(i)=*arg[argV1size+i];
592 return (*func)(argV1, argV2);
593 }
594 double dirDer(const ByteCode::Arg &arg) override {
595 auto argV1size = argV1.size();
596 auto argV2size = argV2.size();
597 for(int i=0; i<argV1size; ++i)
598 argV1(i)=*arg[i];
599 for(int i=0; i<argV2size; ++i)
600 argV2(i)=*arg[argV1size+i];
601 for(size_t i=0; i<argV1size; ++i)
602 dir1V1(i) = *arg[argV1size+argV2size+i];
603 for(size_t i=0; i<argV2size; ++i)
604 dir1V2(i) = *arg[2*argV1size+argV2size+i];
605 return func->dirDer1(dir1V1, argV1, argV2) +
606 func->dirDer2(dir1V2, argV1, argV2);
607 }
608 double dirDerDirDer(const ByteCode::Arg &arg) override {
609 auto argV1size = argV1.size();
610 auto argV2size = argV2.size();
611 for(int i=0; i<argV1size; ++i)
612 argV1(i)=*arg[i];
613 for(int i=0; i<argV2size; ++i)
614 argV2(i)=*arg[argV1size+i];
615 for(size_t i=0; i<argV1size; ++i)
616 dir1V1(i) = *arg[argV1size+argV2size+i];
617 for(size_t i=0; i<argV2size; ++i)
618 dir1V2(i) = *arg[2*argV1size+argV2size+i];
619 for(size_t i=0; i<argV1size; ++i)
620 dir2V1(i) = *arg[2*argV1size+2*argV2size+i];
621 for(size_t i=0; i<argV2size; ++i)
622 dir2V2(i) = *arg[3*argV1size+2*argV2size+i];
623 return func->dirDer1DirDer1(dir1V1, dir2V1, argV1, argV2) +
624 func->dirDer2DirDer1(dir2V2, dir1V1, argV1, argV2) +
625 func->dirDer2DirDer1(dir1V2, dir2V1, argV1, argV2) +
626 func->dirDer2DirDer2(dir1V2, dir2V2, argV1, argV2);
627 }
628 private:
629 std::shared_ptr<fmatvec::Function<double(Vector<Type1,double>,Vector<Type2,double>)>> func;
630 fmatvec::Vector<Type1,double> argV1;
631 fmatvec::Vector<Type2,double> argV2;
632 fmatvec::Vector<Type1,double> dir1V1;
633 fmatvec::Vector<Type2,double> dir1V2;
634 fmatvec::Vector<Type1,double> dir2V1;
635 fmatvec::Vector<Type2,double> dir2V2;
636 };
637
638 // ***** Function *****
639
641 class FMATVEC_EXPORT NativeFunction : public Vertex, public std::enable_shared_from_this<NativeFunction> {
642 public:
643 static SymbolicExpression create(const std::shared_ptr<ScalarFunctionWrapArg> &funcWrapper,
644 const std::vector<SymbolicExpression> &argS,
645 const std::vector<SymbolicExpression> &dir1S={},
646 const std::vector<SymbolicExpression> &dir2S={});
647 SymbolicExpression parDer(const IndependentVariable &x) const override;
648
649 std::vector<ByteCode>::iterator dumpByteCode(std::vector<ByteCode> &byteCode,
650 std::map<const Vertex*,
651 std::vector<AST::ByteCode>::iterator> &existingVertex) const override;
652
653 void walkVertex(const std::function<void(const std::shared_ptr<const Vertex>&)> &func) const override;
654
655 private:
656 NativeFunction(const std::shared_ptr<ScalarFunctionWrapArg> &func_,
657 const std::vector<SymbolicExpression> &argS,
658 const std::vector<SymbolicExpression> &dir1S,
659 const std::vector<SymbolicExpression> &dir2S);
660 bool equal(const SymbolicExpression &b, MapIVSE &m) const override;
661
662 std::shared_ptr<ScalarFunctionWrapArg> funcWrapper;
663 const std::vector<SymbolicExpression> argS;
664 const std::vector<SymbolicExpression> dir1S;
665 const std::vector<SymbolicExpression> dir2S;
666
667 using CacheKey = std::tuple<std::weak_ptr<ScalarFunctionWrapArg>,
668 std::vector<std::weak_ptr<const Vertex>>,
669 std::vector<std::weak_ptr<const Vertex>>,
670 std::vector<std::weak_ptr<const Vertex>> >;
671 struct CacheKeyComp {
672 bool operator()(const CacheKey& l, const CacheKey& r) const;
673 private:
674 std::owner_less<std::weak_ptr<const Vertex>> ol;
675 std::owner_less<std::weak_ptr<const ScalarFunctionWrapArg>> olf;
676 };
677 static std::map<CacheKey, std::weak_ptr<const NativeFunction>, CacheKeyComp> cache;
678 };
679
680 // ***** Operation *****
681
683 class FMATVEC_EXPORT Operation : public Vertex, public std::enable_shared_from_this<Operation> {
684 friend SymbolicExpression;
685 friend SymbolicExpression fmatvec::AST::substScalar(const SymbolicExpression &se,
686 const IndependentVariable& a, const SymbolicExpression &b);
687 friend boost::spirit::qi::rule<boost::spirit::istream_iterator, SymbolicExpression()>&
688 fmatvec::getBoostSpiritQiRule<SymbolicExpression>();
689 friend boost::spirit::karma::rule<std::ostream_iterator<char>, SymbolicExpression()>&
690 fmatvec::getBoostSpiritKarmaRule<SymbolicExpression>();
691 public:
692
694 enum Operator { Plus, Minus, Mult, Div, Pow, Log, Sqrt, Neg, Sin, Cos, Tan, Sinh, Cosh, Tanh, ASin, ACos, ATan, ATan2, ASinh, ACosh, ATanh, Exp, Sign, Heaviside, Abs, Min, Max, Condition };
695 static SymbolicExpression create(Operator op_, const std::vector<SymbolicExpression> &child_);
696 SymbolicExpression parDer(const IndependentVariable &x) const override;
697
698 Operator getOp() const { return op; }
699 const std::vector<SymbolicExpression>& getChilds() const { return child; }
700 std::vector<ByteCode>::iterator dumpByteCode(std::vector<ByteCode> &byteCode,
701 std::map<const Vertex*, std::vector<AST::ByteCode>::iterator> &existingVertex) const override;
702
703 void walkVertex(const std::function<void(const std::shared_ptr<const Vertex>&)> &func) const override;
704
705 private:
706
707 Operation(Operator op_, const std::vector<SymbolicExpression> &child_);
708 bool equal(const SymbolicExpression &b, MapIVSE &m) const override;
709 Operator op;
710 std::vector<SymbolicExpression> child;
711 using CacheKey = std::pair<Operator, std::vector<std::weak_ptr<const Vertex>>>;
712 struct CacheKeyComp {
713 bool operator()(const CacheKey& l, const CacheKey& r) const;
714 private:
715 std::owner_less<std::weak_ptr<const Vertex>> ol;
716 };
717 static std::map<CacheKey, std::weak_ptr<const Operation>, CacheKeyComp> cache;
718
719 struct OpMap {
720 std::string funcName; // used to dump/read the expression using boost spirit/karma: e.g.
721 // "plus"
722 std::function<void(double*, const ByteCode::Arg&)> func; // used for runtime evaluation: e.g.
723 // [](double* r, const ByteCode::Arg& a){ *r = *a[0] + *a[1]; }
724 };
725 static const std::map<Operator, OpMap> opMap;
726 };
727
728 inline SymbolicExpression SymbolicFuncWrapArg1<double(double), SymbolicExpression>::call(
729 const std::shared_ptr<fmatvec::Function<double(double)>> &func,
730 const SymbolicExpression &arg) {
731 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgS>(func), std::vector<SymbolicExpression>{arg});
732 }
733
734 template<class Type>
735 SymbolicExpression SymbolicFuncWrapArg1<double(Vector<Type, double>), Vector<Type, SymbolicExpression>>::call(
736 const std::shared_ptr<Function<double(Vector<Type, double>)>> &func,
737 const Vector<Type, SymbolicExpression> &arg) {
738 std::vector<SymbolicExpression> argV(arg.size());
739 for(int i=0; i<arg.size(); ++i)
740 argV[i]=arg(i);
741 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgV<Type>>(func, arg.size()), argV);
742 }
743
744 inline SymbolicExpression SymbolicFuncWrapArg2<double(double,double), SymbolicExpression, SymbolicExpression>::call(
745 const std::shared_ptr<fmatvec::Function<double(double,double)>> &func,
746 const SymbolicExpression &arg1, const SymbolicExpression &arg2) {
747 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgSS>(func), std::vector<SymbolicExpression>{arg1,arg2});
748 }
749
750 template<class Type>
751 SymbolicExpression SymbolicFuncWrapArg2<double(Vector<Type,double>,double), Vector<Type, SymbolicExpression>, SymbolicExpression>::call(
752 const std::shared_ptr<Function<double(Vector<Type,double>,double)>> &func,
753 const Vector<Type, SymbolicExpression> &arg1, const SymbolicExpression &arg2) {
754 std::vector<SymbolicExpression> arg(arg1.size()+1);
755 for(int i=0; i<arg1.size(); ++i)
756 arg[i]=arg1(i);
757 arg[arg1.size()]=arg2;
758 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgVS<Type>>(func, arg1.size()), arg);
759 }
760
761 template<class Type>
762 SymbolicExpression SymbolicFuncWrapArg2<double(double,Vector<Type,double>), SymbolicExpression, Vector<Type, SymbolicExpression>>::call(
763 const std::shared_ptr<Function<double(double,Vector<Type,double>)>> &func,
764 const SymbolicExpression &arg1, const Vector<Type, SymbolicExpression> &arg2) {
765 std::vector<SymbolicExpression> arg(1+arg2.size());
766 arg[0]=arg1;
767 for(int i=0; i<arg2.size(); ++i)
768 arg[i+1]=arg2(i);
769 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgSV<Type>>(func, arg2.size()), arg);
770 }
771
772 template<class Type1, class Type2>
773 SymbolicExpression SymbolicFuncWrapArg2<double(Vector<Type1,double>,Vector<Type2,double>), Vector<Type1, SymbolicExpression>, Vector<Type2, SymbolicExpression>>::call(
774 const std::shared_ptr<Function<double(Vector<Type1,double>,Vector<Type2,double>)>> &func,
775 const Vector<Type1, SymbolicExpression> &arg1, const Vector<Type2, SymbolicExpression> &arg2) {
776 std::vector<SymbolicExpression> arg(arg1.size()+arg2.size());
777 for(int i=0; i<arg1.size(); ++i)
778 arg[i]=arg1(i);
779 for(int i=0; i<arg2.size(); ++i)
780 arg[arg1.size()+i]=arg2(i);
781 return AST::NativeFunction::create(std::make_shared<AST::ScalarFunctionWrapArgVV<Type1,Type2>>(func, arg1.size(), arg2.size()), arg);
782 }
783
784 } // end namespace AST
785 #endif
786
787 inline bool operator<(const fmatvec::SymbolicExpression& a, const fmatvec::SymbolicExpression& b) noexcept {
788 // See LessIV why we need to overload this operator< and call abort!
789 std::cerr<<"ABORT: See LessIV: have you called std::min/max instead of fmatvev::min/max with type SymbolicExpression?"<<std::endl;
790 assert(((void)"See LessIV: have you called std::min/max instead of fmatvev::min/max with type SymbolicExpression?", false));
791 abort();
792 }
793
794 IndependentVariable& IndependentVariable::operator^=(double x) {
795 static_cast<const AST::Symbol*>(get())->setValue(x);
796 return *this;
797 }
798
799 template<> boost::spirit::qi::rule<boost::spirit::istream_iterator, IndependentVariable()>& getBoostSpiritQiRule<IndependentVariable>();
800 template<> boost::spirit::qi::rule<boost::spirit::istream_iterator, SymbolicExpression()>& getBoostSpiritQiRule<SymbolicExpression>();
801
802 template<> boost::spirit::karma::rule<std::ostream_iterator<char>, IndependentVariable()>& getBoostSpiritKarmaRule<IndependentVariable>();
803 template<> boost::spirit::karma::rule<std::ostream_iterator<char>, SymbolicExpression()>& getBoostSpiritKarmaRule<SymbolicExpression>();
804
805 // for vectors/matrices of type IndependentVariable/SymbolicExpression throw exceptions instead of assert.
806 // these exceptions are thrown regardless whether NDEBUG is defined or not.
807 // this is done since vectors/matrix operations of this type are usually user input, should not abort/assert on errors,
808 // and it not time critical since it just build the expression once which is than evaluated at runtime.
809 template<> struct AssertUseException<IndependentVariable> { constexpr static bool value = true; };
810 template<> struct AssertUseException<SymbolicExpression> { constexpr static bool value = true; };
811
812 } // end namespace fmatvec
813
814 #endif
fmatvec::AST::Symbol::setValue
void setValue(double x_) const
Definition: ast.h:402
fmatvec::IndependentVariable::operator^=
IndependentVariable & operator^=(double x)
Set the double value of the independent value.
Definition: ast.h:794
fmatvec::AST::SymbolicFuncWrapArg1
Definition: ast.h:224
fmatvec::AST::Constant
A vertex of the AST representing a constant (long or double)
Definition: ast.h:25
fmatvec::AST::ScalarFunctionWrapArgSS
Definition: ast.h:469
fmatvec::AST::Constant::getValue
const T & getValue() const
Get the constant value of the vertex.
Definition: ast.h:367
fmatvec::AST::Constant::isConstantInt
bool isConstantInt() const override
Rreturn true if this Vertex is a constant integer.
fmatvec::AST::NativeFunction::CacheKeyComp
Definition: ast.h:671
fmatvec::AST::Vertex::isConstantInt
virtual bool isConstantInt() const
Rreturn true if this Vertex is a constant integer.
Definition: ast.h:325
fmatvec::AST::ByteCode
Definition: ast.h:208
fmatvec::AST::Operation::CacheKeyComp
Definition: ast.h:712
fmatvec::AST::SymbolicFuncWrapArg2
Definition: ast.h:245
fmatvec::AST::ScalarFunctionWrapArgV
Definition: ast.h:432
fmatvec::Function
Definition: function.h:202
fmatvec::AST::ScalarFunctionWrapArg
Definition: ast.h:408
fmatvec::AST::ScalarFunctionWrapArgVV
Definition: ast.h:574
fmatvec::AST::Operation::OpMap
Definition: ast.h:719
fmatvec::AST::NativeFunction
A vertex of the AST representing an arbitary function.
Definition: ast.h:641
std::shared_ptr
Definition: types.h:81
fmatvec::Vector
Definition: fmatvec.h:48
fmatvec::IndependentVariable
Definition: ast.h:107
fmatvec::AST::Symbol
A vertex of the AST representing a independent variable.
Definition: ast.h:374
fmatvec::Eval
Definition: ast.h:29
boost::spirit::basic_istream_iterator
Definition: stream.h:32
fmatvec::AST::Operation
A vertex of the AST representing an operation.
Definition: ast.h:683
fmatvec::AssertUseException
Definition: types.h:219
fmatvec
Namespace fmatvec.
Definition: _memory.cc:28
fmatvec::AST::ScalarFunctionWrapArgVS
Definition: ast.h:490
fmatvec::AST::ScalarFunctionWrapArgSV
Definition: ast.h:532
fmatvec::AST::Vertex
A abstract class for a Vertex of the AST (abstract syntax tree).
Definition: ast.h:282
fmatvec::SymbolicExpression
Definition: ast.h:37
fmatvec::AST::Operation::Operator
Operator
Defined operations.
Definition: ast.h:694
fmatvec::operator<<
ostream & operator<<(ostream &s, const SymbolicExpression &se)
Write a SymbolicExpression to a stream using serialization.
Definition: ast.cc:190
fmatvec::AST::ScalarFunctionWrapArgS
Definition: ast.h:415
fmatvec::parDer
SymbolicExpression parDer(const SymbolicExpression &dep, const IndependentVariable &indep)
Definition: ast.cc:299
fmatvec::AST::Vertex::LessIV
Definition: ast.h:316