INMOST
A toolkit for distributed mathematical modeling
inmost_expression.h
1 #ifndef INMOST_AUTODIFF_ETEXPR_H_INCLUDED
2 #define INMOST_AUTODIFF_ETEXPR_H_INCLUDED
3 #include "inmost_common.h"
4 #include "inmost_sparse.h"
5 #include <sstream> //for debug
6 #include <new>
7 
8 //TODO:
9 // 1. Incorporate tables
10 // 2. (ok, test) implement condition
11 // 3. (ok, test) implement stencil
12 // 4. (???) copying of basic_dynamic_variable
13 // 5. Consider optimization by checking zero variation multipliers, check that assembly do not degrade.
14 // 6. floor, ceil, atan, acos, asin, max, min functions
15 // 7. choice of directional derivatives at discontinuities for abs, pow, max, min (see ADOL-C)
16 // 8. replace stencil with foreach for provided iterators
17 // 9. enclose in namespace
18 //10. CheckCurrentAutomatizator -> CheckCurrentRowMerger
19 //10.0 structure/service to handle multiple RowMerger objects in openmp environment
20 //10.1 user should be able to provide RowMerger when Automatizator is not compiled
21 //10.2 Automatizator may provide internal structure for RowMerger
22 
23 
24 #ifdef _MSC_VER
25 #pragma warning(disable : 4503)
26 #endif
27 
28 #define CNT_USE_MERGER 0
29 //#define CNT_USE_MERGER ENUMUNDEF
30 
31 #if defined(USE_AUTODIFF)
32 namespace INMOST
33 {
34 
35 
36  class basic_expression
37  {
38  protected:
39  static thread_private<Sparse::RowMerger> merger;
40  public:
41  basic_expression() {}
42  virtual INMOST_DATA_REAL_TYPE GetValue() const = 0;
43  virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const = 0;
44  virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const = 0;
45  virtual void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const = 0;
46  virtual INMOST_DATA_ENUM_TYPE GetCount() const = 0;
47  virtual ~basic_expression() {}
48  };
49 
50  template<class Derived>
51  class shell_expression : virtual public basic_expression
52  {
53  public:
54  shell_expression() {}
55  __INLINE virtual INMOST_DATA_REAL_TYPE GetValue() const {return static_cast<const Derived *>(this)->GetValue(); }
56  __INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const { if( mult ) return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
57  __INLINE virtual void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const { if( mult ) return static_cast<const Derived *>(this)->GetJacobian(mult,r); }
58  __INLINE virtual void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {return static_cast<const Derived *>(this)->GetHessian(multJ,J,multH,H); }
59  __INLINE virtual INMOST_DATA_ENUM_TYPE GetCount() const { return static_cast<const Derived*>(this)->GetCount(); }
60  operator Derived & () {return *static_cast<Derived *>(this);}
61  operator const Derived & () const {return *static_cast<const Derived *>(this);}
62  ~shell_expression() {}
63  };
64 
65 
66  class const_expression : public shell_expression<const_expression>
67  {
68  INMOST_DATA_REAL_TYPE value;
69  public:
70  const_expression(const const_expression & other) :value(other.value) {}
71  const_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue) {}
72  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
73  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const {(void)mult; (void)r;}
74  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const {(void)mult; (void)r;}
75  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {(void)multJ; (void)J; (void)multH; (void)H;}
76  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return 0; }
77  __INLINE const_expression & operator =(const_expression const & other)
78  {
79  value = other.value;
80  return *this;
81  }
82  };
83 
84 
85  class var_expression : public shell_expression<var_expression>
86  {
87  INMOST_DATA_REAL_TYPE value;
88  INMOST_DATA_ENUM_TYPE index;
89  public:
90  var_expression() : value(0), index(ENUMUNDEF) {}
91  var_expression(const var_expression & other) :value(other.value), index(other.index) {}
92  var_expression(INMOST_DATA_REAL_TYPE pvalue, INMOST_DATA_ENUM_TYPE pindex) : value(pvalue), index(pindex) {}
93  var_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue), index(ENUMUNDEF) {}
94  __INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
95  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
96  __INLINE INMOST_DATA_ENUM_TYPE GetIndex() const { return index; }
97  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const {if( mult && index != ENUMUNDEF ) r[index] += mult;}
98  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const {if( mult && index != ENUMUNDEF ) r[index] += mult;}
99  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const {if( index != ENUMUNDEF ) J.Push(index,multJ); (void)multH; (void)H;}
100  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return index == ENUMUNDEF ? 0 : 1; }
101  __INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE i) const {return index == i? 1.0 : 0.0;}
102  __INLINE var_expression & operator =(var_expression const & other)
103  {
104  value = other.value;
105  index = other.index;
106  return *this;
107  }
108  __INLINE var_expression & operator =(INMOST_DATA_REAL_TYPE other)
109  {
110  value = other;
111  index = ENUMUNDEF;
112  return *this;
113  }
114  bool check_nans() const
115  {
116  return value != value;
117  }
118  bool check_infs() const
119  {
120  return __isinf__(value);
121  }
122  };
123 
124 #if defined(PACK_ARRAY)
125 #pragma pack(push,r1,4)
126 #endif
130  class multivar_expression : public shell_expression<multivar_expression>
131  {
132  INMOST_DATA_REAL_TYPE value; //< Value of the variable.
133  Sparse::Row entries; //< Sparse vector of variations.
134  public:
135  multivar_expression() :value(0) {}
136  multivar_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue) {}
137  multivar_expression(const multivar_expression & other) : value(other.value), entries(other.entries) {}
138  multivar_expression(INMOST_DATA_REAL_TYPE pvalue, Sparse::Row & pentries)
139  : value(pvalue), entries(pentries) {}
140  multivar_expression(INMOST_DATA_REAL_TYPE pvalue, INMOST_DATA_ENUM_TYPE pindex, INMOST_DATA_REAL_TYPE pdmult = 1.0)
141  : value(pvalue)
142  {
143  entries.Push(pindex,pdmult);
144  }
145  multivar_expression(const basic_expression & expr)
146  {
147  value = expr.GetValue();
148  INMOST_DATA_ENUM_TYPE cnt = expr.GetCount();
149  //if (cnt >= CNT_USE_MERGER)
150  {
151  merger->Resize(cnt);
152  expr.GetJacobian(1.0, *merger);
153  merger->RetrieveRow(entries);
154  merger->Clear();
155  }
156  //else expr.GetJacobian(1.0, entries);
157  }
158  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
159  __INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
160  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
161  {
162  r.AddRow(mult, entries);
163  }
164  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
165  {
166  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + r.Size();
167  if (cnt >= CNT_USE_MERGER)
168  {
169  merger->Resize(cnt);
170  merger->AddRow(mult, entries);
171  merger->AddRow(1.0, r);
172  merger->RetrieveRow(r);
173  merger->Clear();
174  }
175  else for (Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
176  r[it->first] += it->second * mult;
177  }
178  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
179  {
180  J = entries;
181  if( !J.isSorted() ) std::sort(J.Begin(),J.End());
182  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= multJ;
183  H.Clear();
184  (void)multH;
185  }
186  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return entries.Size(); }
187  __INLINE multivar_expression & operator = (INMOST_DATA_REAL_TYPE pvalue)
188  {
189  value = pvalue;
190  entries.Clear();
191  return *this;
192  }
193  __INLINE multivar_expression & operator = (multivar_expression const & other)
194  {
195  value = other.value;
196  entries = other.entries;
197  return *this;
198  }
199  /*
200  __INLINE multivar_expression & operator = (basic_expression const & expr)
201  {
202  value = expr.GetValue();
203  if( CheckCurrentAutomatizator() )
204  FromBasicExpression(entries,expr);
205  else
206  {
207  Sparse::Row tmp;
208  expr.GetJacobian(1.0,tmp);
209  entries.Swap(tmp);
210  }
211  return *this;
212  }
213  */
214  __INLINE Sparse::Row & GetRow() {return entries;}
215  __INLINE const Sparse::Row & GetRow() const {return entries;}
216  __INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE index) const {return GetRow().get_safe(index);}
217  __INLINE bool operator < (INMOST_DATA_REAL_TYPE right) const {return value < right;}
218  __INLINE bool operator > (INMOST_DATA_REAL_TYPE right) const {return value > right;}
219  __INLINE bool operator <=(INMOST_DATA_REAL_TYPE right) const {return value <= right;}
220  __INLINE bool operator >=(INMOST_DATA_REAL_TYPE right) const {return value >= right;}
221  __INLINE bool operator ==(INMOST_DATA_REAL_TYPE right) const {return value == right;}
222  __INLINE bool operator !=(INMOST_DATA_REAL_TYPE right) const {return value != right;}
223  __INLINE bool operator < (basic_expression const & expr) const {return value < expr.GetValue();}
224  __INLINE bool operator > (basic_expression const & expr) const {return value > expr.GetValue();}
225  __INLINE bool operator <=(basic_expression const & expr) const {return value <= expr.GetValue();}
226  __INLINE bool operator >=(basic_expression const & expr) const {return value >= expr.GetValue();}
227  __INLINE bool operator ==(basic_expression const & expr) const {return value == expr.GetValue();}
228  __INLINE bool operator !=(basic_expression const & expr) const {return value != expr.GetValue();}
229  __INLINE bool operator < (multivar_expression const & expr) const {return value < expr.GetValue();}
230  __INLINE bool operator > (multivar_expression const & expr) const {return value > expr.GetValue();}
231  __INLINE bool operator <=(multivar_expression const & expr) const {return value <= expr.GetValue();}
232  __INLINE bool operator >=(multivar_expression const & expr) const {return value >= expr.GetValue();}
233  __INLINE bool operator ==(multivar_expression const & expr) const {return value == expr.GetValue();}
234  __INLINE bool operator !=(multivar_expression const & expr) const {return value != expr.GetValue();}
235  __INLINE multivar_expression & operator +=(basic_expression const & expr)
236  {
237  value += expr.GetValue();
238  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + expr.GetCount();
239  //if (cnt >= CNT_USE_MERGER)
240  {
241  merger->Resize(cnt);
242  merger->PushRow(1.0, entries);
243  expr.GetJacobian(1.0, *merger);
244  merger->RetrieveRow(entries);
245  merger->Clear();
246  }
247  //else expr.GetJacobian(1.0, entries);
248  return *this;
249  }
250  __INLINE multivar_expression & operator -=(basic_expression const & expr)
251  {
252  value -= expr.GetValue();
253  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + expr.GetCount();
254  //if (cnt >= CNT_USE_MERGER)
255  {
256  merger->Resize(cnt);
257  merger->PushRow(1.0, entries);
258  expr.GetJacobian(-1.0, *merger);
259  merger->RetrieveRow(entries);
260  merger->Clear();
261  }
262  //else expr.GetJacobian(-1.0, entries);
263  return *this;
264  }
265  __INLINE multivar_expression & operator *=(basic_expression const & expr)
266  {
267  INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
268  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + expr.GetCount();
269  //if (cnt >= CNT_USE_MERGER)
270  {
271  merger->Resize(cnt);
272  merger->PushRow(rval, entries);
273  expr.GetJacobian(lval, *merger);
274  merger->RetrieveRow(entries);
275  merger->Clear();
276  }
277  //else
278  //{
279  // for (Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
280  // it->second *= rval;
281  // expr.GetJacobian(lval, entries);
282  //}
283  value *= rval;
284  return *this;
285  }
286  __INLINE multivar_expression & operator /=(basic_expression const & expr)
287  {
288  INMOST_DATA_REAL_TYPE rval = expr.GetValue();
289  INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
290  value *= reciprocial_rval;
291  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + expr.GetCount();
292  //if (cnt >= CNT_USE_MERGER)
293  {
294  merger->Resize(cnt);
295  merger->PushRow(reciprocial_rval, entries);
296  expr.GetJacobian(-value * reciprocial_rval, *merger);
297  merger->RetrieveRow(entries);
298  merger->Clear();
299  }
300  //else
301  //{
302  // for (Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it)
303  // it->second *= reciprocial_rval;
304  // expr.GetJacobian(-value * reciprocial_rval, entries);
305  //}
306  return *this;
307  }
308  __INLINE multivar_expression & operator +=(INMOST_DATA_REAL_TYPE right)
309  {
310  value += right;
311  return *this;
312  }
313  __INLINE multivar_expression & operator -=(INMOST_DATA_REAL_TYPE right)
314  {
315  value -= right;
316  return *this;
317  }
318  __INLINE multivar_expression & operator *=(INMOST_DATA_REAL_TYPE right)
319  {
320  value *= right;
321  for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second *= right;
322  return *this;
323  }
324  __INLINE multivar_expression & operator /=(INMOST_DATA_REAL_TYPE right)
325  {
326  value /= right;
327  for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second /= right;
328  return *this;
329  }
330  bool check_nans() const
331  {
332  if( value != value ) return true;
333  for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
334  if( it->second != it->second ) return true;
335  return false;
336  }
337  bool check_infs() const
338  {
339  if( __isinf__(value) ) return true;
340  for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
341  if( __isinf__(it->second) ) return true;
342  return false;
343  }
349  INMOST_DATA_ENUM_TYPE Record(Sparse::Row::entry * v) const
350  {
351  INMOST_DATA_ENUM_TYPE k = 0;
352  v[k].first = entries.Size();
353  v[k].second = value;
354  k++;
355  for(INMOST_DATA_ENUM_TYPE r = 0; r < entries.Size(); ++r)
356  {
357  v[k].first = entries.GetIndex(r);
358  v[k].second = entries.GetValue(r);
359  k++;
360  }
361  return k;
362  }
364  INMOST_DATA_ENUM_TYPE RecordSize() const
365  {
366  return 1 + entries.Size();
367  }
372  INMOST_DATA_ENUM_TYPE Retrieve(const Sparse::Row::entry * v)
373  {
374  int k = 0;
375  value = v[k].second;
376  entries.Resize(v[k].first);
377  k++;
378  for(int r = 0; r < (int)entries.Size(); ++r)
379  {
380  entries.GetIndex(r) = v[k].first;
381  entries.GetValue(r) = v[k].second;
382  k++;
383  }
384  return k;
385  }
387  static INMOST_DATA_ENUM_TYPE RetrieveSize(const Sparse::Row::entry * v)
388  {
389  return 1 + v[0].first;
390  }
391  void Print(double eps = -1, std::ostream & sout = std::cout) const
392  {
393  sout << value << " ";// std::endl;
394  entries.Print(eps, sout);
395  }
396  void swap(multivar_expression & b)
397  {
398  std::swap(value,b.value);
399  entries.Swap(b.entries);
400  }
401  friend class multivar_expression_reference;
402  };
403 
407  class hessian_multivar_expression : public shell_expression<hessian_multivar_expression>
408  {
409  INMOST_DATA_REAL_TYPE value;
410  Sparse::Row entries;
411  Sparse::HessianRow hessian_entries;
412  public:
413  void swap(hessian_multivar_expression & b)
414  {
415  std::swap(value,b.value);
416  entries.Swap(b.entries);
417  hessian_entries.Swap(b.hessian_entries);
418  }
420  hessian_multivar_expression() :value(0) {}
422  hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue) : value(pvalue) {}
424  hessian_multivar_expression(const hessian_multivar_expression & other) : value(other.value), entries(other.entries), hessian_entries(other.hessian_entries) {}
426  //hessian_multivar_expression(const multivar_expression & other) : value(other.GetValue()), entries(other.GetRow()) {}
428  hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, Sparse::Row & pentries)
429  : value(pvalue), entries(pentries) {}
431  hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, Sparse::Row & pentries, Sparse::HessianRow & phentries)
432  : value(pvalue), entries(pentries), hessian_entries(phentries) {}
434  hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, INMOST_DATA_ENUM_TYPE pindex, INMOST_DATA_REAL_TYPE pdmult = 1.0)
435  : value(pvalue)
436  {
437  entries.Push(pindex,pdmult);
438  }
440  hessian_multivar_expression(const basic_expression & expr)
441  {
442  value = expr.GetValue();
443  expr.GetHessian(1.0,entries,1.0,hessian_entries);
444  }
445  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
446  __INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
447  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
448  {
449  r.AddRow(mult,entries);
450  //for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
451  // r[it->first] += it->second*mult;
452  }
453  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
454  {
455  INMOST_DATA_ENUM_TYPE cnt = entries.Size() + r.Size();
456  if (cnt >= CNT_USE_MERGER)
457  {
458  merger->Resize(cnt);
459  merger->AddRow(mult, entries);
460  merger->AddRow(1.0, r);
461  merger->RetrieveRow(r);
462  merger->Clear();
463  }
464  else for (Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
465  r[it->first] += it->second * mult;
466  }
467  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
468  {
469  J = entries;
470  if( !J.isSorted() ) std::sort(J.Begin(),J.End());
471  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= multJ;
472  H = hessian_entries;
473  for(Sparse::HessianRow::iterator it = H.Begin(); it != H.End(); ++it) it->second *= multH;
474  }
475  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return entries.Size(); }
476  __INLINE multivar_expression GetVariable(INMOST_DATA_ENUM_TYPE index)
477  {
478  multivar_expression ret(0);
479  for(int k = 0; k < (int)entries.Size(); ++k)
480  {
481  if( entries.GetIndex(k) == index )
482  {
483  ret.SetValue(entries.GetValue(k));
484  Sparse::Row & r = ret.GetRow();
485  for(int q = 0; q < (int)hessian_entries.Size(); ++q)
486  {
487  Sparse::HessianRow::index & i = hessian_entries.GetIndex(q);
488  if( i.first == index )
489  r.Push(i.second,hessian_entries.GetValue(q)*(i.first == i.second ? 1.0 : 0.5));
490  else if( i.second == index )
491  r.Push(i.first,hessian_entries.GetValue(q)*(i.first == i.second ? 1.0 : 0.5));
492  }
493  break;
494  }
495  }
496  return ret;
497  }
498  __INLINE hessian_multivar_expression & operator = (INMOST_DATA_REAL_TYPE pvalue)
499  {
500  value = pvalue;
501  entries.Clear();
502  hessian_entries.Clear();
503  return *this;
504  }
505  /*
506  __INLINE hessian_multivar_expression & operator = (basic_expression const & expr)
507  {
508  value = expr.GetValue();
509  Sparse::Row tmp;
510  Sparse::HessianRow htmp;
511  expr.GetHessian(1.0,tmp,1.0,htmp);
512  entries.Swap(tmp);
513  hessian_entries.Swap(htmp);
514  return *this;
515  }
516  */
517  /*
518  __INLINE hessian_multivar_expression & operator = (multivar_expression const & other)
519  {
520  value = other.GetValue();
521  entries = other.GetRow();
522  hessian_entries.Clear();
523  return *this;
524  }
525  */
526  __INLINE hessian_multivar_expression & operator = (hessian_multivar_expression const & other)
527  {
528  value = other.value;
529  entries = other.entries;
530  hessian_entries = other.hessian_entries;
531  return *this;
532  }
533  __INLINE Sparse::Row & GetRow() {return entries;}
534  __INLINE Sparse::HessianRow & GetHessianRow() {return hessian_entries;}
535  __INLINE const Sparse::Row & GetRow() const {return entries;}
536  __INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE index) const {return GetRow().get_safe(index);}
537  __INLINE const Sparse::HessianRow & GetHessianRow() const {return hessian_entries;}
538  __INLINE hessian_multivar_expression & operator +=(basic_expression const & expr)
539  {
540  value += expr.GetValue();
541  Sparse::Row tmpr, tmp;
542  Sparse::HessianRow htmpr, htmp;
543  expr.GetHessian(1.0,tmpr,1.0,htmpr);
544  Sparse::Row::MergeSortedRows(1.0,entries,1.0,tmpr,tmp);
545  entries.Swap(tmp);
546  Sparse::HessianRow::MergeSortedRows(1.0,hessian_entries,1.0,htmpr,htmp);
547  hessian_entries.Swap(htmp);
548  return *this;
549  }
550  __INLINE hessian_multivar_expression & operator -=(basic_expression const & expr)
551  {
552  value -= expr.GetValue();
553  Sparse::Row tmpr, tmp;
554  Sparse::HessianRow htmpr, htmp;
555  expr.GetHessian(1.0,tmpr,1.0,htmpr);
556  Sparse::Row::MergeSortedRows(1.0,entries,-1.0,tmpr,tmp);
557  entries.Swap(tmp);
558  Sparse::HessianRow::MergeSortedRows(1.0,hessian_entries,-1.0,htmpr,htmp);
559  hessian_entries.Swap(htmp);
560  return *this;
561  }
562  __INLINE hessian_multivar_expression & operator *=(basic_expression const & expr)
563  {
564  throw NotImplemented; //check code below is correct
565  INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
566  Sparse::Row tmp, tmpr;
567  Sparse::HessianRow htmp, htmpr;
568  expr.GetHessian(1.0,tmpr,1.0,htmpr);
569  Sparse::Row::MergeSortedRows(rval,entries,lval,tmpr,tmp);
570  entries.Swap(tmp);
571  Sparse::HessianRow::MergeSortedRows(lval,hessian_entries,rval,htmpr,htmp);
572  hessian_entries.Swap(htmp);
573  return *this;
574  }
575  __INLINE hessian_multivar_expression & operator /=(basic_expression const & expr)
576  {
577  throw NotImplemented; //check code below is correct
578  INMOST_DATA_REAL_TYPE rval = expr.GetValue();
579  INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
580  value *= reciprocial_rval;
581  Sparse::Row tmp, tmpr;
582  Sparse::HessianRow htmp, htmpr;
583  expr.GetHessian(-value*reciprocial_rval,tmpr,1.0,htmpr);
584  Sparse::Row::MergeSortedRows(reciprocial_rval,entries,-value*reciprocial_rval,tmpr,tmp);
585  entries.Swap(tmp);
586  Sparse::HessianRow::MergeSortedRows(reciprocial_rval,hessian_entries,-value*reciprocial_rval,htmpr,htmp);
587  hessian_entries.Swap(htmp);
588  return *this;
589  }
590  __INLINE hessian_multivar_expression & operator +=(INMOST_DATA_REAL_TYPE right)
591  {
592  value += right;
593  return *this;
594  }
595  __INLINE hessian_multivar_expression & operator -=(INMOST_DATA_REAL_TYPE right)
596  {
597  value -= right;
598  return *this;
599  }
600  __INLINE hessian_multivar_expression & operator *=(INMOST_DATA_REAL_TYPE right)
601  {
602  value *= right;
603  for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second *= right;
604  for(Sparse::HessianRow::iterator it = hessian_entries.Begin(); it != hessian_entries.End(); ++it) it->second *= right;
605  return *this;
606  }
607  __INLINE hessian_multivar_expression & operator /=(INMOST_DATA_REAL_TYPE right)
608  {
609  value /= right;
610  for(Sparse::Row::iterator it = entries.Begin(); it != entries.End(); ++it) it->second /= right;
611  for(Sparse::HessianRow::iterator it = hessian_entries.Begin(); it != hessian_entries.End(); ++it) it->second /= right;
612  return *this;
613  }
614  bool check_nans() const
615  {
616  if( value != value ) return true;
617  for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
618  if( it->second != it->second ) return true;
619  for(Sparse::HessianRow::const_iterator it = hessian_entries.Begin(); it != hessian_entries.End(); ++it)
620  if( it->second != it->second ) return true;
621  return false;
622  }
623  bool check_infs() const
624  {
625  if( __isinf__(value)) return true;
626  for(Sparse::Row::const_iterator it = entries.Begin(); it != entries.End(); ++it)
627  if( __isinf__(it->second) ) return true;
628  for(Sparse::HessianRow::const_iterator it = hessian_entries.Begin(); it != hessian_entries.End(); ++it)
629  if( __isinf__(it->second) ) return true;
630  return false;
631  }
632  friend class hessian_multivar_expression_reference;
633  };
634 #if defined(PACK_ARRAY)
635 #pragma pack(pop,r1)
636 #endif
637 
638  static INMOST_DATA_REAL_TYPE stub_multivar_expression_reference_value; //for default constructor in multivar_expression_reference
639 
640 
641  class multivar_expression_reference : public shell_expression<multivar_expression_reference>
642  {
643  INMOST_DATA_REAL_TYPE & value;
644  Sparse::Row * entries;
645  public:
647  multivar_expression_reference() : value(stub_multivar_expression_reference_value), entries(NULL) {}
649  multivar_expression_reference(INMOST_DATA_REAL_TYPE & _value, Sparse::Row * _entries)
650  : value(_value), entries(_entries) {}
652  multivar_expression_reference(const multivar_expression_reference & other)
653  : value(other.value), entries(other.entries) {}
655  multivar_expression_reference(multivar_expression & other)
656  : value(other.value), entries(&other.entries) {}
658  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
660  __INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
662  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
663  {
664  if( entries )
665  {
666  r.AddRow(mult,*entries);
667  //~ for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
668  //~ r[it->first] += it->second*mult;
669  }
670  }
672  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
673  {
674  if( entries )
675  {
676  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + r.Size();
677  if (cnt >= CNT_USE_MERGER)
678  {
679  merger->Resize(cnt);
680  merger->AddRow(mult, *entries);
681  merger->AddRow(1.0, r);
682  merger->RetrieveRow(r);
683  merger->Clear();
684  }
685  else for (Sparse::Row::const_iterator it = entries->Begin(); it != entries->End(); ++it)
686  r[it->first] += it->second * mult;
687  }
688  }
689  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J,INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
690  {
691  if( entries )
692  {
693  J = *entries;
694  if( !J.isSorted() ) std::sort(J.Begin(),J.End());
695  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= multJ;
696  H.Clear();
697  }
698  (void)multH;
699  }
700  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return entries->Size(); }
701  multivar_expression_reference & operator = (INMOST_DATA_REAL_TYPE pvalue)
702  {
703  value = pvalue;
704  entries->Clear();
705  return *this;
706  }
707 
708  __INLINE multivar_expression_reference & operator = (basic_expression const & expr)
709  {
710  value = expr.GetValue();
711  if (entries)
712  {
713  INMOST_DATA_ENUM_TYPE cnt = expr.GetCount();
714  //if (cnt > CNT_USE_MERGER)
715  {
716  merger->Resize(cnt);
717  expr.GetJacobian(1.0, *merger);
718  merger->RetrieveRow(*entries);
719  merger->Clear();
720  }
721  //else
722  //{
723  // entries->Clear();
724  // expr.GetJacobian(1.0, *entries);
725  //}
726  }
727  return *this;
728  }
729 
730  __INLINE multivar_expression_reference & operator = (multivar_expression_reference const & other)
731  {
732  value = other.GetValue();
733  *entries = other.GetRow();
734  return *this;
735  }
736  /*
737  __INLINE multivar_expression_reference & operator = (multivar_expression const & other)
738  {
739  value = other.GetValue();
740  *entries = other.GetRow();
741  return *this;
742  }
743  */
744  __INLINE Sparse::Row & GetRow() {return *entries;}
745  __INLINE const Sparse::Row & GetRow() const {return *entries;}
746  __INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE index) const {return GetRow().get_safe(index);}
747  __INLINE multivar_expression_reference & operator +=(basic_expression const & expr)
748  {
749  value += expr.GetValue();
750  if (entries)
751  {
752  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + expr.GetCount();
753  //if (cnt >= CNT_USE_MERGER)
754  {
755  merger->Resize(cnt);
756  merger->PushRow(1.0, *entries);
757  expr.GetJacobian(1.0, *merger);
758  merger->RetrieveRow(*entries);
759  merger->Clear();
760  }
761  //else expr.GetJacobian(1.0, *entries);
762  }
763  return *this;
764  }
765  __INLINE multivar_expression_reference & operator -=(basic_expression const & expr)
766  {
767  value -= expr.GetValue();
768  if (entries)
769  {
770  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + expr.GetCount();
771  //if (cnt >= CNT_USE_MERGER)
772  {
773  merger->Resize(cnt);
774  merger->PushRow(1.0, *entries);
775  expr.GetJacobian(-1.0, *merger);
776  merger->RetrieveRow(*entries);
777  merger->Clear();
778  }
779  //else expr.GetJacobian(-1.0, *entries);
780  }
781  return *this;
782  }
783  __INLINE multivar_expression_reference & operator *=(basic_expression const & expr)
784  {
785  INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
786  if (entries)
787  {
788  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + expr.GetCount();
789  //if (cnt >= CNT_USE_MERGER)
790  {
791  merger->Resize(cnt);
792  merger->PushRow(rval, *entries);
793  expr.GetJacobian(lval, *merger);
794  merger->RetrieveRow(*entries);
795  merger->Clear();
796  }
797  //else
798  //{
799  // for (Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
800  // it->second *= rval;
801  // expr.GetJacobian(lval, *entries);
802  //}
803  }
804  value *= rval;
805  return *this;
806  }
807  __INLINE multivar_expression_reference & operator /=(basic_expression const & expr)
808  {
809  INMOST_DATA_REAL_TYPE rval = expr.GetValue();
810  INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
811  value *= reciprocial_rval;
812  if (entries)
813  {
814  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + expr.GetCount();
815  //if (cnt >= CNT_USE_MERGER)
816  {
817  merger->Resize(cnt);
818  merger->PushRow(reciprocial_rval, *entries);
819  expr.GetJacobian(-value * reciprocial_rval, *merger);
820  merger->RetrieveRow(*entries);
821  merger->Clear();
822  }
823  //else
824  //{
825  // for (Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
826  // it->second *= reciprocial_rval;
827  // expr.GetJacobian(-value * reciprocial_rval, *entries);
828  //}
829  }
830  return *this;
831  }
832  __INLINE multivar_expression_reference & operator +=(INMOST_DATA_REAL_TYPE right)
833  {
834  value += right;
835  return *this;
836  }
837  __INLINE multivar_expression_reference & operator -=(INMOST_DATA_REAL_TYPE right)
838  {
839  value -= right;
840  return *this;
841  }
842  __INLINE multivar_expression_reference & operator *=(INMOST_DATA_REAL_TYPE right)
843  {
844  value *= right;
845  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second *= right;
846  return *this;
847  }
848  __INLINE multivar_expression_reference & operator /=(INMOST_DATA_REAL_TYPE right)
849  {
850  value /= right;
851  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second /= right;
852  return *this;
853  }
854  bool check_nans() const
855  {
856  if( value != value ) return true;
857  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
858  if( it->second != it->second ) return true;
859  return false;
860  }
861  bool check_infs() const
862  {
863  if( __isinf__(value) ) return true;
864  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
865  if( __isinf__(it->second) ) return true;
866  return false;
867  }
868  };
869 
870  class value_reference
871  {
872  INMOST_DATA_REAL_TYPE& value;
873  public:
875  value_reference() : value(stub_multivar_expression_reference_value) {}
877  value_reference(INMOST_DATA_REAL_TYPE& _value)
878  : value(_value) {}
880  value_reference(const value_reference& other)
881  : value(other.value) {}
882  __INLINE value_reference& operator = (INMOST_DATA_REAL_TYPE pvalue)
883  {
884  value = pvalue;
885  return *this;
886  }
887  __INLINE value_reference& operator = (value_reference const& other)
888  {
889  value = other.value;
890  return *this;
891  }
892  /*
893  __INLINE value_reference& operator = (basic_expression const& expr)
894  {
895  value = expr.GetValue();
896  return *this;
897  }
898  __INLINE value_reference& operator +=(basic_expression const& expr)
899  {
900  value += expr.GetValue();
901  return *this;
902  }
903  __INLINE value_reference& operator -=(basic_expression const& expr)
904  {
905  value -= expr.GetValue();
906  return *this;
907  }
908  __INLINE value_reference& operator *=(basic_expression const& expr)
909  {
910  value *= expr.GetValue();
911  return *this;
912  }
913  __INLINE value_reference& operator /=(basic_expression const& expr)
914  {
915  value /= expr.GetValue();
916  return *this;
917  }
918  */
919  __INLINE value_reference& operator +=(INMOST_DATA_REAL_TYPE right)
920  {
921  value += right;
922  return *this;
923  }
924  __INLINE value_reference& operator -=(INMOST_DATA_REAL_TYPE right)
925  {
926  value -= right;
927  return *this;
928  }
929  __INLINE value_reference& operator *=(INMOST_DATA_REAL_TYPE right)
930  {
931  value *= right;
932  return *this;
933  }
934  __INLINE value_reference& operator /=(INMOST_DATA_REAL_TYPE right)
935  {
936  value /= right;
937  return *this;
938  }
939  bool check_nans() const
940  {
941  if (value != value) return true;
942  return false;
943  }
944  bool check_infs() const
945  {
946  if (__isinf__(value)) return true;
947  return false;
948  }
949  operator INMOST_DATA_REAL_TYPE() const { return value; }
950  operator INMOST_DATA_REAL_TYPE& () { return value; }
951  };
952 
953 
954 
955  class hessian_multivar_expression_reference : public shell_expression<hessian_multivar_expression_reference>
956  {
957  INMOST_DATA_REAL_TYPE & value;
958  Sparse::Row * entries;
959  Sparse::HessianRow * hentries;
960  public:
962  hessian_multivar_expression_reference(INMOST_DATA_REAL_TYPE & _value, Sparse::Row * _entries, Sparse::HessianRow * _hentries)
963  : value(_value), entries(_entries), hentries(_hentries) {}
965  hessian_multivar_expression_reference(const hessian_multivar_expression_reference & other)
966  : value(other.value), entries(other.entries), hentries(other.hentries) {}
968  hessian_multivar_expression_reference(hessian_multivar_expression & other)
969  : value(other.value), entries(&other.entries), hentries(&other.hessian_entries) {}
971  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
973  __INLINE void SetValue(INMOST_DATA_REAL_TYPE val) { value = val; }
975  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
976  {
977  r.AddRow(mult,*entries);
978  //~ for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
979  //~ r[it->first] += it->second*mult;
980  }
982  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
983  {
984  INMOST_DATA_ENUM_TYPE cnt = entries->Size() + r.Size();
985  if (cnt >= CNT_USE_MERGER)
986  {
987  merger->Resize(cnt);
988  merger->AddRow(mult, *entries);
989  merger->AddRow(1.0, r);
990  merger->RetrieveRow(r);
991  merger->Clear();
992  }
993  else for (Sparse::Row::const_iterator it = entries->Begin(); it != entries->End(); ++it)
994  r[it->first] += it->second * mult;
995  }
996  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J,INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
997  {
998  J = *entries;
999  if( !J.isSorted() ) std::sort(J.Begin(),J.End());
1000  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= multJ;
1001  H = *hentries;
1002  for(Sparse::HessianRow::iterator it = H.Begin(); it != H.End(); ++it) it->second *= multH;
1003  }
1004  __INLINE INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_ENUM_TYPE index) const {return GetRow().get_safe(index);}
1005  __INLINE multivar_expression GetVariable(INMOST_DATA_ENUM_TYPE index)
1006  {
1007  multivar_expression ret(0);
1008  for(int k = 0; k < (int)entries->Size(); ++k)
1009  {
1010  if( entries->GetIndex(k) == index )
1011  {
1012  ret.SetValue(entries->GetValue(k));
1013  Sparse::Row & r = ret.GetRow();
1014  for(int q = 0; q < (int)hentries->Size(); ++q)
1015  {
1016  Sparse::HessianRow::index & i = hentries->GetIndex(q);
1017  if( i.first == index )
1018  r.Push(i.second,hentries->GetValue(q)*(i.first == i.second ? 1.0 : 0.5));
1019  else if( i.second == index )
1020  r.Push(i.first,hentries->GetValue(q)*(i.first == i.second ? 1.0 : 0.5));
1021  }
1022  break;
1023  }
1024  }
1025  return ret;
1026  }
1027  __INLINE hessian_multivar_expression_reference & operator = (INMOST_DATA_REAL_TYPE pvalue)
1028  {
1029  value = pvalue;
1030  entries->Clear();
1031  hentries->Clear();
1032  return *this;
1033  }
1034  __INLINE hessian_multivar_expression_reference & operator = (basic_expression const & expr)
1035  {
1036  value = expr.GetValue();
1037  Sparse::Row tmp;
1038  Sparse::HessianRow htmp;
1039  expr.GetHessian(1.0,tmp,1.0,htmp);
1040  entries->Swap(tmp);
1041  hentries->Swap(htmp);
1042  return *this;
1043  }
1044  __INLINE hessian_multivar_expression_reference & operator = (multivar_expression_reference const & other)
1045  {
1046  value = other.GetValue();
1047  *entries = other.GetRow();
1048  hentries->Clear();
1049  return *this;
1050  }
1051  __INLINE hessian_multivar_expression_reference & operator = (multivar_expression const & other)
1052  {
1053  value = other.GetValue();
1054  *entries = other.GetRow();
1055  hentries->Clear();
1056  return *this;
1057  }
1058  __INLINE hessian_multivar_expression_reference & operator = (hessian_multivar_expression_reference const & other)
1059  {
1060  value = other.GetValue();
1061  *entries = other.GetRow();
1062  *hentries = other.GetHessianRow();
1063  return *this;
1064  }
1065  __INLINE hessian_multivar_expression_reference & operator = (hessian_multivar_expression const & other)
1066  {
1067  value = other.GetValue();
1068  *entries = other.GetRow();
1069  *hentries = other.GetHessianRow();
1070  return *this;
1071  }
1072  __INLINE Sparse::Row & GetRow() {return *entries;}
1073  __INLINE Sparse::HessianRow & GetHessianRow() {return *hentries;}
1074  __INLINE const Sparse::Row & GetRow() const {return *entries;}
1075  __INLINE const Sparse::HessianRow & GetHessianRow() const {return *hentries;}
1076  __INLINE hessian_multivar_expression_reference & operator +=(basic_expression const & expr)
1077  {
1078  value += expr.GetValue();
1079  Sparse::Row tmpr, tmp;
1080  Sparse::HessianRow htmpr, htmp;
1081  expr.GetHessian(1.0,tmpr,1.0,htmpr);
1082  Sparse::Row::MergeSortedRows(1.0,*entries,1.0,tmpr,tmp);
1083  entries->Swap(tmp);
1084  Sparse::HessianRow::MergeSortedRows(1.0,*hentries,1.0,htmpr,htmp);
1085  hentries->Swap(htmp);
1086  return *this;
1087  }
1088  __INLINE hessian_multivar_expression_reference & operator -=(basic_expression const & expr)
1089  {
1090  value -= expr.GetValue();
1091  Sparse::Row tmpr, tmp;
1092  Sparse::HessianRow htmpr, htmp;
1093  expr.GetHessian(1.0,tmpr,1.0,htmpr);
1094  Sparse::Row::MergeSortedRows(1.0,*entries,-1.0,tmpr,tmp);
1095  entries->Swap(tmp);
1096  Sparse::HessianRow::MergeSortedRows(1.0,*hentries,-1.0,htmpr,htmp);
1097  hentries->Swap(htmp);
1098  return *this;
1099  }
1100  __INLINE hessian_multivar_expression_reference & operator *=(basic_expression const & expr)
1101  {
1102  throw NotImplemented; //check code below is correct
1103  INMOST_DATA_REAL_TYPE lval = value, rval = expr.GetValue();
1104  Sparse::Row tmp, tmpr;
1105  Sparse::HessianRow htmp, htmpr;
1106  expr.GetHessian(1.0,tmpr,1.0,htmpr);
1107  Sparse::Row::MergeSortedRows(rval,*entries,lval,tmpr,tmp);
1108  entries->Swap(tmp);
1109  Sparse::HessianRow::MergeSortedRows(lval,*hentries,rval,htmpr,htmp);
1110  hentries->Swap(htmp);
1111  value *= rval;
1112  return *this;
1113  }
1114  __INLINE hessian_multivar_expression_reference & operator /=(basic_expression const & expr)
1115  {
1116  throw NotImplemented; //check code below is correct
1117  INMOST_DATA_REAL_TYPE rval = expr.GetValue();
1118  INMOST_DATA_REAL_TYPE reciprocial_rval = 1.0/rval;
1119  value *= reciprocial_rval;
1120  Sparse::Row tmp, tmpr;
1121  Sparse::HessianRow htmp, htmpr;
1122  expr.GetHessian(-value*reciprocial_rval,tmpr,1.0,htmpr);
1123  Sparse::Row::MergeSortedRows(reciprocial_rval,*entries,-value*reciprocial_rval,tmpr,tmp);
1124  entries->Swap(tmp);
1125  Sparse::HessianRow::MergeSortedRows(reciprocial_rval,*hentries,-value*reciprocial_rval,htmpr,htmp);
1126  hentries->Swap(htmp);
1127  return *this;
1128  }
1129  __INLINE hessian_multivar_expression_reference & operator +=(INMOST_DATA_REAL_TYPE right)
1130  {
1131  value += right;
1132  return *this;
1133  }
1134  __INLINE hessian_multivar_expression_reference & operator -=(INMOST_DATA_REAL_TYPE right)
1135  {
1136  value -= right;
1137  return *this;
1138  }
1139  __INLINE hessian_multivar_expression_reference & operator *=(INMOST_DATA_REAL_TYPE right)
1140  {
1141  value *= right;
1142  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second *= right;
1143  for(Sparse::HessianRow::iterator it = hentries->Begin(); it != hentries->End(); ++it)
1144  it->second *= right;
1145  return *this;
1146  }
1147  __INLINE hessian_multivar_expression_reference & operator /=(INMOST_DATA_REAL_TYPE right)
1148  {
1149  value /= right;
1150  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it) it->second /= right;
1151  for(Sparse::HessianRow::iterator it = hentries->Begin(); it != hentries->End(); ++it) it->second /= right;
1152  return *this;
1153  }
1154  bool check_nans() const
1155  {
1156  if( value != value ) return true;
1157  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
1158  if( it->second != it->second ) return true;
1159  for(Sparse::HessianRow::iterator it = hentries->Begin(); it != hentries->End(); ++it)
1160  if( it->second != it->second ) return true;
1161  return false;
1162  }
1163  bool check_infs() const
1164  {
1165  if( __isinf__(value) ) return true;
1166  for(Sparse::Row::iterator it = entries->Begin(); it != entries->End(); ++it)
1167  if( __isinf__(it->second) ) return true;
1168  for(Sparse::HessianRow::iterator it = hentries->Begin(); it != hentries->End(); ++it)
1169  if( __isinf__(it->second) ) return true;
1170  return false;
1171  }
1172  };
1173 
1174 
1175  template<class A>
1176  class const_multiplication_expression : public shell_expression<const_multiplication_expression<A> >
1177  {
1178  const A & arg;
1179  INMOST_DATA_REAL_TYPE value, dmult;
1180  public:
1181  const_multiplication_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE pdmult) : arg(parg), dmult(pdmult)
1182  {
1183  value = arg.GetValue()*dmult;
1184  }
1185  const_multiplication_expression(const const_multiplication_expression & other) : arg(other.arg), value(other.value), dmult(other.dmult) {}
1186  const_multiplication_expression(const const_multiplication_expression & other, const A & parg) : arg(parg), value(other.value), dmult(other.dmult) {}
1187  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1188  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1189  {
1190  arg.GetJacobian(mult*dmult,r);
1191  }
1192  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1193  {
1194  arg.GetJacobian(mult*dmult,r);
1195  }
1196  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1197  {
1198  arg.GetHessian(multJ*dmult,J,multH*dmult,H);
1199  }
1200  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1201  };
1202 
1203  template<class A>
1204  class variation_multiplication_expression : public shell_expression<variation_multiplication_expression<A> >
1205  {
1206  const A & arg;
1207  INMOST_DATA_REAL_TYPE value, dmult;
1208  public:
1209  variation_multiplication_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE pdmult) : arg(parg), dmult(pdmult)
1210  {
1211  value = arg.GetValue();
1212  }
1213  variation_multiplication_expression(const variation_multiplication_expression & other) : arg(other.arg), value(other.value), dmult(other.dmult) {}
1214  variation_multiplication_expression(const variation_multiplication_expression & other, const A & parg) : arg(parg), value(other.value), dmult(other.dmult) {}
1215  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1216  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1217  {
1218  arg.GetJacobian(mult*dmult,r);
1219  }
1220  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1221  {
1222  arg.GetJacobian(mult*dmult,r);
1223  }
1224  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1225  {
1226  arg.GetHessian(multJ*dmult,J,multH*dmult,H);
1227  }
1228  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1229  };
1230 
1231 
1232  template<class A>
1233  class const_division_expression : public shell_expression<const_division_expression<A> >
1234  {
1235  const A & arg;
1236  INMOST_DATA_REAL_TYPE value, dmult;
1237  public:
1238  const_division_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE pdmult) : arg(parg), dmult(pdmult)
1239  {
1240  dmult = 1.0/dmult;
1241  value = arg.GetValue()*dmult;
1242  }
1243  const_division_expression(const const_division_expression & other) : arg(other.arg), value(other.value), dmult(other.dmult) {}
1244  const_division_expression(const const_division_expression & other, const A & parg) : arg(parg), value(other.value), dmult(other.dmult) {}
1245  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1246  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1247  {
1248  arg.GetJacobian(mult*dmult,r);
1249  }
1250  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1251  {
1252  arg.GetJacobian(mult*dmult,r);
1253  }
1254  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1255  {
1256  arg.GetHessian(multJ*dmult,J,multH*dmult,H);
1257  }
1258  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1259  };
1260 
1261  template<class A>
1262  class const_addition_expression : public shell_expression<const_addition_expression<A> >
1263  {
1264  const A & arg;
1265  INMOST_DATA_REAL_TYPE value;
1266  public:
1267  const_addition_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE padd) : arg(parg)
1268  {
1269  value = arg.GetValue()+padd;
1270  }
1271  const_addition_expression(const const_addition_expression & other) : arg(other.arg), value(other.value) {}
1272  const_addition_expression(const const_addition_expression & other, const A & parg) : arg(parg), value(other.value) {}
1273  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1274  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1275  {
1276  arg.GetJacobian(mult,r);
1277  }
1278  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1279  {
1280  arg.GetJacobian(mult,r);
1281  }
1282  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1283  {
1284  arg.GetHessian(multJ,J,multH,H);
1285  }
1286  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1287  };
1288 
1289  template<class A>
1290  class const_subtraction_expression : public shell_expression<const_subtraction_expression<A> >
1291  {
1292  const A & arg;
1293  INMOST_DATA_REAL_TYPE value;
1294  public:
1295  const_subtraction_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE pleft) : arg(parg)
1296  {
1297  value = pleft-arg.GetValue();
1298  }
1299  const_subtraction_expression(const const_subtraction_expression & other) : arg(other.arg), value(other.value) {}
1300  const_subtraction_expression(const const_subtraction_expression & other, const A & parg) : arg(parg), value(other.value) {}
1301  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1302  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1303  {
1304  arg.GetJacobian(-mult,r);
1305  }
1306  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1307  {
1308  arg.GetJacobian(-mult,r);
1309  }
1310  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1311  {
1312  arg.GetHessian(-multJ,J,-multH,H);
1313  }
1314  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1315  };
1316 
1319  template<class A>
1320  class reciprocal_expression : public shell_expression<reciprocal_expression<A> >
1321  {
1322  const A & arg;
1323  INMOST_DATA_REAL_TYPE value, reciprocial_val;
1324  public:
1325  reciprocal_expression(const shell_expression<A> & parg,INMOST_DATA_REAL_TYPE coef) : arg(parg)
1326  {
1327  reciprocial_val = 1.0/arg.GetValue();
1328  value = coef*reciprocial_val;
1329  }
1330  reciprocal_expression(const reciprocal_expression & other)
1331  : arg(other.arg), value(other.value), reciprocial_val(other.reciprocial_val) {}
1332  reciprocal_expression(const reciprocal_expression & other, const A & parg)
1333  : arg(parg), value(other.value), reciprocial_val(other.reciprocial_val) {}
1334  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1335  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1336  {
1337  arg.GetJacobian(-mult*value*reciprocial_val,r);
1338  }
1339  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1340  {
1341  arg.GetJacobian(-mult*value*reciprocial_val,r);
1342  }
1343  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1344  {
1345  Sparse::HessianRow ArgH;
1346  double coefJ = -multJ*value*reciprocial_val;
1347  double signJ = coefJ < 0 ? -1 : 1;
1348  arg.GetHessian(signJ,J,-2*multH*value*reciprocial_val,ArgH);
1349  Sparse::HessianRow::MergeJacobianHessian(2*value*reciprocial_val*reciprocial_val*signJ,J,J,1.0,ArgH,H);
1350  for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
1351  }
1352  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1353  };
1354 
1355  template<class A>
1356  class unary_minus_expression : public shell_expression<unary_minus_expression<A> >
1357  {
1358  const A & arg;
1359  INMOST_DATA_REAL_TYPE value;
1360  public:
1361  unary_minus_expression(const shell_expression<A> & parg) : arg(parg) {value = -arg.GetValue();}
1362  unary_minus_expression(const unary_minus_expression & b) : arg(b.arg) {}
1363  unary_minus_expression(const unary_minus_expression & b, const A & parg) : arg(parg) {}
1364  __INLINE INMOST_DATA_REAL_TYPE GetValue() const {return value;}
1365  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1366  {
1367  arg.GetJacobian(-mult,r);
1368  }
1369  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1370  {
1371  arg.GetJacobian(-mult,r);
1372  }
1373  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1374  {
1375  arg.GetHessian(-multJ,J,-multH,H);
1376  }
1377  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1378  };
1379 
1380  template<class A>
1381  class unary_plus_expression : public shell_expression<unary_plus_expression<A> >
1382  {
1383  const A & arg;
1384  INMOST_DATA_REAL_TYPE value;
1385  public:
1386  unary_plus_expression(const shell_expression<A> & parg) : arg(parg) {value = arg.GetValue();}
1387  unary_plus_expression(const unary_plus_expression & b) : arg(b.arg) {}
1388  unary_plus_expression(const unary_plus_expression & b, const A & parg) : arg(parg) {}
1389  __INLINE INMOST_DATA_REAL_TYPE GetValue() const {return value;}
1390  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1391  {
1392  arg.GetJacobian(mult,r);
1393  }
1394  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1395  {
1396  arg.GetJacobian(mult,r);
1397  }
1398  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1399  {
1400  arg.GetHessian(multJ,J,multH,H);
1401  }
1402  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1403  };
1404 
1405  template<class A>
1406  class abs_expression : public shell_expression<abs_expression<A> >
1407  {
1408  const A & arg;
1409  INMOST_DATA_REAL_TYPE value, dmult;
1410  public:
1411  abs_expression(const shell_expression<A> & parg) : arg(parg)
1412  {
1413  value = arg.GetValue();
1414  dmult = value < 0.0 ? -1.0 : 1.0;
1415  value = ::fabs(value);
1416  }
1417  abs_expression(const abs_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1418  abs_expression(const abs_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1419  __INLINE INMOST_DATA_REAL_TYPE GetValue() const {return value;}
1420  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1421  {
1422  arg.GetJacobian( (value == 0 ? (mult < 0.0 ? -1 : 1) : 1) * mult * dmult, r);
1423  }
1424  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1425  {
1426  arg.GetJacobian( (value == 0 ? (mult < 0.0 ? -1 : 1) : 1) * mult * dmult, r);
1427  }
1428  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ,Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1429  {
1430  double a = (value == 0 ? (multJ < 0.0 ? -1 : 1) : 1);
1431  double b = (value == 0 ? (multH < 0.0 ? -1 : 1) : 1);
1432  arg.GetHessian( a * multJ * dmult,J, b*multH * dmult,H);
1433  }
1434  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1435  };
1436 
1437  // ex
1438  template<class A>
1439  class exp_expression : public shell_expression<exp_expression<A> >
1440  {
1441  const A & arg;
1442  INMOST_DATA_REAL_TYPE value;
1443  public:
1444  exp_expression(const shell_expression<A> & parg) : arg(parg)
1445  {
1446  value = arg.GetValue();
1447  value = ::exp(value);
1448  }
1449  exp_expression(const exp_expression & b) : arg(b.arg), value(b.value) {}
1450  exp_expression(const exp_expression & b, const A & parg) : arg(parg), value(b.value) {}
1451  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1452  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1453  {
1454  arg.GetJacobian( mult * value, r);
1455  }
1456  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1457  {
1458  arg.GetJacobian( mult * value, r);
1459  }
1460  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1461  {
1462  Sparse::HessianRow ArgH;
1463  double coefJ = multJ*value;
1464  double signJ = coefJ < 0.0 ? -1 : 1;
1465  arg.GetHessian(signJ, J, multH*value, ArgH); //check
1466  Sparse::HessianRow::MergeJacobianHessian(coefJ*signJ,J,J,1.0,ArgH,H);
1467  for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
1468  }
1469  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1470  };
1471 
1472  template<class A>
1473  class log_expression : public shell_expression<log_expression<A> >
1474  {
1475  const A & arg;
1476  INMOST_DATA_REAL_TYPE value, dmult;
1477  public:
1478  log_expression(const shell_expression<A> & parg) : arg(parg)
1479  {
1480  value = arg.GetValue();
1481  dmult = 1.0/value;
1482  value = ::log(value);
1483  }
1484  log_expression(const log_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1485  log_expression(const log_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1486  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1487  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1488  {
1489  arg.GetJacobian(mult*dmult,r);
1490  }
1491  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1492  {
1493  arg.GetJacobian(mult*dmult,r);
1494  }
1495  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1496  {
1497  Sparse::HessianRow ArgH;
1498  double coefJ = multJ*dmult;
1499  double signJ = coefJ < 0.0 ? -1 : 1;
1500  arg.GetHessian(signJ, J, 2*multH*dmult, ArgH); //check
1501  Sparse::HessianRow::MergeJacobianHessian(-coefJ*signJ*dmult,J,J,1.0,ArgH,H);
1502  for(INMOST_DATA_ENUM_TYPE k = 0; k < J.Size(); ++k) J.GetValue(k) *= coefJ*signJ;
1503  }
1504  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1505  };
1506 
1507 
1508  template<class A>
1509  class sin_expression : public shell_expression<sin_expression<A> >
1510  {
1511  const A & arg;
1512  INMOST_DATA_REAL_TYPE value, dmult;
1513  public:
1514  sin_expression(const shell_expression<A> & parg) : arg(parg)
1515  {
1516  value = arg.GetValue();
1517  dmult = ::cos(value);
1518  value = ::sin(value);
1519  }
1520  sin_expression(const sin_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1521  sin_expression(const sin_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1522  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1523  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1524  {
1525  arg.GetJacobian(mult*dmult,r);
1526  }
1527  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1528  {
1529  arg.GetJacobian(mult*dmult,r);
1530  }
1531  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1532  {
1533  Sparse::HessianRow htmp;
1534  arg.GetHessian(1,J,1,htmp);
1535  assert(J.isSorted());
1536  assert(htmp.isSorted());
1537  Sparse::HessianRow::MergeJacobianHessian(-value*multH,J,J,dmult*multH,htmp,H);
1538  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult*multJ;
1539  assert(H.isSorted());
1540  }
1541  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1542  };
1543 
1544  template<class A>
1545  class cos_expression : public shell_expression<cos_expression<A> >
1546  {
1547  const A & arg;
1548  INMOST_DATA_REAL_TYPE value, dmult;
1549  public:
1550  cos_expression(const shell_expression<A> & parg) : arg(parg)
1551  {
1552  value = arg.GetValue();
1553  dmult = -(::sin(value));
1554  value = ::cos(value);
1555  }
1556  cos_expression(const cos_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1557  cos_expression(const cos_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1558  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1559  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1560  {
1561  arg.GetJacobian(mult*dmult,r);
1562  }
1563  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1564  {
1565  arg.GetJacobian(mult*dmult,r);
1566  }
1567  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1568  {
1569  //arg.GetHessian(multJ*dmult,J,-multH*value,H);
1570  Sparse::HessianRow htmp;
1571  arg.GetHessian(1,J,1,htmp);
1572  Sparse::HessianRow::MergeJacobianHessian(-value*multH,J,J,dmult*multH,htmp,H);
1573  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second*=dmult*multJ;
1574  }
1575  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1576  };
1577 
1578  template<class A>
1579  class sqrt_expression : public shell_expression<sqrt_expression<A> >
1580  {
1581  const A & arg;
1582  INMOST_DATA_REAL_TYPE value;
1583  public:
1584  sqrt_expression(const shell_expression<A> & parg) : arg(parg)
1585  {
1586  value = ::sqrt(arg.GetValue());
1587  }
1588  sqrt_expression(const sqrt_expression & b) : arg(b.arg), value(b.value) {}
1589  sqrt_expression(const sqrt_expression & b, const A & parg) : arg(parg), value(b.value) {}
1590  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1591  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1592  {
1593  if( value ) arg.GetJacobian(0.5*mult/value,r);
1594  }
1595  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1596  {
1597  if( value ) arg.GetJacobian(0.5*mult/value,r);
1598  }
1599  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1600  {
1601  //general formula:
1602  // (F(G))'' = F'(G) G'' + F''(G) G'.G'
1603  if( value )
1604  {
1605  Sparse::HessianRow htmp;
1606  arg.GetHessian(1,J,1,htmp);
1607  Sparse::HessianRow::MergeJacobianHessian(-0.25/::pow(value,3.0)*multH,J,J,0.5/value*multH,htmp,H);
1608  for(Sparse::Row::iterator it = J.Begin(); it != J.End(); ++it) it->second *= 0.5/value*multJ;
1609  }
1610  //arg.GetHessian(0.5*multJ/value,J,-0.25*multH/::pow(value,3),H);
1611  }
1612  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1613  };
1614 
1615 
1616  template<class A>
1617  class soft_abs_expression : public shell_expression<soft_abs_expression<A> >
1618  {
1619  const A & arg;
1620  INMOST_DATA_REAL_TYPE value, dmult;
1621  public:
1622  soft_abs_expression(const shell_expression<A> & parg, INMOST_DATA_REAL_TYPE tol) : arg(parg)
1623  {
1624  INMOST_DATA_REAL_TYPE lval = arg.GetValue();
1625  value = ::sqrt(lval*lval+tol*tol);
1626  dmult = value ? lval/value : 0.0;
1627  }
1628  soft_abs_expression(const soft_abs_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1629  soft_abs_expression(const soft_abs_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1630  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1631  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1632  {
1633  arg.GetJacobian(mult*dmult,r);
1634  }
1635  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1636  {
1637  arg.GetJacobian(mult*dmult,r);
1638  }
1639  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1640  {
1641  (void)multJ,(void)J,(void)multH,(void)H;
1642  throw NotImplemented;
1643  }
1644  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1645  };
1646 
1647  template<class A>
1648  class soft_sign_expression : public shell_expression<soft_sign_expression<A> >
1649  {
1650  const A & arg;
1651  INMOST_DATA_REAL_TYPE value, dmult;
1652  public:
1653  soft_sign_expression(const shell_expression<A> & parg, INMOST_DATA_REAL_TYPE tol) : arg(parg)
1654  {
1655  INMOST_DATA_REAL_TYPE lval = arg.GetValue(), lval2 = lval*lval;
1656  INMOST_DATA_REAL_TYPE div = lval2+tol*tol;
1657  INMOST_DATA_REAL_TYPE sdiv = ::sqrt(div);
1658  value = sdiv ? lval/sdiv : 0.0;
1659  dmult = sdiv ? (1.0 - (div ? lval2/div : 0.0))/sdiv : 0.0;
1660  }
1661  soft_sign_expression(const soft_sign_expression & b) : arg(b.arg), value(b.value), dmult(b.dmult) {}
1662  soft_sign_expression(const soft_sign_expression & b, const A & parg) : arg(parg), value(b.value), dmult(b.dmult) {}
1663  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; };
1664  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1665  {
1666  arg.GetJacobian(mult*dmult,r);
1667  }
1668  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1669  {
1670  arg.GetJacobian(mult*dmult,r);
1671  }
1672  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1673  {
1674  (void)multJ,(void)J,(void)multH,(void)H;
1675  throw NotImplemented;
1676  }
1677  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
1678  };
1679 
1680  template<class A, class B>
1681  class soft_max_expression : public shell_expression<soft_max_expression<A,B> >
1682  {
1683  const A & left;
1684  const B & right;
1685  INMOST_DATA_REAL_TYPE value, ldmult, rdmult;
1686  public:
1687  soft_max_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright, INMOST_DATA_REAL_TYPE tol) : left(pleft), right(pright)
1688  {
1689  INMOST_DATA_REAL_TYPE lval = left.GetValue(), rval = right.GetValue();
1690  INMOST_DATA_REAL_TYPE diff = lval-rval, root = ::sqrt(diff*diff+tol*tol);
1691  value = 0.5*(lval+rval + root);
1692  if (root)
1693  {
1694  ldmult = 0.5 * (1 + diff / root);
1695  rdmult = 0.5 * (1 - diff / root);
1696  }
1697  else ldmult = rdmult = 0.5;
1698  }
1699  soft_max_expression(const soft_max_expression & other)
1700  : left(other.left), right(other.right), value(other.value), ldmult(other.ldmult), rdmult(other.rdmult) {}
1701  soft_max_expression(const soft_max_expression & other, const A & pleft, const B & pright)
1702  : left(pleft), right(pright), value(other.value), ldmult(other.ldmult), rdmult(other.rdmult) {}
1703 
1704  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1705  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1706  {
1707  left.GetJacobian(mult*ldmult,r);
1708  right.GetJacobian(mult*rdmult,r);
1709  }
1710  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1711  {
1712  left.GetJacobian(mult*ldmult,r);
1713  right.GetJacobian(mult*rdmult,r);
1714  }
1715  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1716  {
1717  (void)multJ,(void)J,(void)multH,(void)H;
1718  throw NotImplemented;
1719  }
1720  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
1721  };
1722 
1723  template<class A>
1724  class soft_max_const_expression : public shell_expression<soft_max_const_expression<A> >
1725  {
1726  const A& left;
1727  INMOST_DATA_REAL_TYPE right;
1728  INMOST_DATA_REAL_TYPE value, ldmult;
1729  public:
1730  soft_max_const_expression(const shell_expression<A>& pleft, INMOST_DATA_REAL_TYPE pright, INMOST_DATA_REAL_TYPE tol) : left(pleft), right(pright)
1731  {
1732  INMOST_DATA_REAL_TYPE lval = left.GetValue(), rval = right;
1733  INMOST_DATA_REAL_TYPE diff = lval - rval, root = ::sqrt(diff * diff + tol * tol);
1734  value = 0.5 * (lval + rval + root);
1735  ldmult = 0.5 * (1 + (root ? diff / root : 0.0));
1736  }
1737  soft_max_const_expression(const soft_max_const_expression& other)
1738  : left(other.left), right(other.right), value(other.value), ldmult(other.ldmult) {}
1739  soft_max_const_expression(const soft_max_const_expression& other, const A& pleft, INMOST_DATA_REAL_TYPE pright)
1740  : left(pleft), right(pright), value(other.value), ldmult(other.ldmult) {}
1741  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1742  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger& r) const {left.GetJacobian(mult * ldmult, r);}
1743  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row& r) const {left.GetJacobian(mult * ldmult, r);}
1744  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row& J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow& H) const
1745  {
1746  (void)multJ, (void)J, (void)multH, (void)H;
1747  throw NotImplemented;
1748  }
1749  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount(); }
1750  };
1751 
1752  template<class A, class B>
1753  class soft_min_expression : public shell_expression<soft_min_expression<A,B> >
1754  {
1755  const A & left;
1756  const B & right;
1757  INMOST_DATA_REAL_TYPE value, ldmult, rdmult;
1758  public:
1759  soft_min_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright, INMOST_DATA_REAL_TYPE tol) : left(pleft), right(pright)
1760  {
1761  INMOST_DATA_REAL_TYPE lval = left.GetValue(), rval = right.GetValue();
1762  INMOST_DATA_REAL_TYPE diff = lval-rval, root = ::sqrt(diff*diff+tol*tol);
1763  value = 0.5*(lval+rval - root);
1764  if (root)
1765  {
1766  ldmult = 0.5 * (1 - diff / root);
1767  rdmult = 0.5 * (1 + diff / root);
1768  }
1769  else ldmult = rdmult = 0.5;
1770  }
1771  soft_min_expression(const soft_min_expression & other)
1772  : left(other.left), right(other.right), value(other.value), ldmult(other.ldmult), rdmult(other.rdmult) {}
1773  soft_min_expression(const soft_min_expression & other, const A & pleft, const B & pright)
1774  : left(pleft), right(pright), value(other.value), ldmult(other.ldmult), rdmult(other.rdmult) {}
1775  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1776  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1777  {
1778  left.GetJacobian(mult*ldmult,r);
1779  right.GetJacobian(mult*rdmult,r);
1780  }
1781  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1782  {
1783  left.GetJacobian(mult*ldmult,r);
1784  right.GetJacobian(mult*rdmult,r);
1785  }
1786  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1787  {
1788  (void)multJ,(void)J,(void)multH,(void)H;
1789  throw NotImplemented;
1790  }
1791  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
1792  };
1793 
1794 
1795  template<class A>
1796  class soft_min_const_expression : public shell_expression<soft_min_const_expression<A> >
1797  {
1798  const A& left;
1799  INMOST_DATA_REAL_TYPE right;
1800  INMOST_DATA_REAL_TYPE value, ldmult;
1801  public:
1802  soft_min_const_expression(const shell_expression<A>& pleft, INMOST_DATA_REAL_TYPE pright, INMOST_DATA_REAL_TYPE tol) : left(pleft), right(pright)
1803  {
1804  INMOST_DATA_REAL_TYPE lval = left.GetValue(), rval = right;
1805  INMOST_DATA_REAL_TYPE diff = lval - rval, root = ::sqrt(diff * diff + tol * tol);
1806  value = 0.5 * (lval + rval - root);
1807  ldmult = 0.5 * (1 - (root ? diff / root : 0.0));
1808  }
1809  soft_min_const_expression(const soft_min_const_expression& other)
1810  : left(other.left), right(other.right), value(other.value), ldmult(other.ldmult) {}
1811  soft_min_const_expression(const soft_min_const_expression& other, const A& pleft, INMOST_DATA_REAL_TYPE pright)
1812  : left(pleft), right(pright), value(other.value), ldmult(other.ldmult) {}
1813  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1814  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger& r) const {left.GetJacobian(mult * ldmult, r);}
1815  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row& r) const { left.GetJacobian(mult * ldmult, r); }
1816  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row& J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow& H) const
1817  {
1818  (void)multJ, (void)J, (void)multH, (void)H;
1819  throw NotImplemented;
1820  }
1821  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount(); }
1822  };
1823 
1824  template<class A, class B>
1825  class multiplication_expression : public shell_expression<multiplication_expression<A,B> >
1826  {
1827  const A & left;
1828  const B & right;
1829  INMOST_DATA_REAL_TYPE value;
1830  public:
1831  multiplication_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
1832  {
1833  value = left.GetValue()*right.GetValue();
1834  }
1835  multiplication_expression(const multiplication_expression & other)
1836  : left(other.left), right(other.right), value(other.value) {}
1837  multiplication_expression(const multiplication_expression & other, const A & pleft, const B & pright)
1838  : left(pleft), right(pright), value(other.value) {}
1839  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1840  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1841  {
1842  left.GetJacobian(mult*right.GetValue(),r);
1843  right.GetJacobian(mult*left.GetValue(),r);
1844  }
1845  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1846  {
1847  left.GetJacobian(mult*right.GetValue(),r);
1848  right.GetJacobian(mult*left.GetValue(),r);
1849  }
1850  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1851  {
1852  // (F*G)'' = (F'G+G'F)' = (F''G + F'G' + G''F + G'F') = (F''G + G''F + 2F'G')
1853  Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
1854  Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
1855  left.GetHessian(1,JL,1,HL); //retrieve jacobian row and hessian matrix of the left expression
1856  assert(JL.isSorted());
1857  assert(HL.isSorted());
1858  right.GetHessian(1,JR,1,HR); //retrieve jacobian row and hessian matrix of the right expression
1859  assert(JR.isSorted());
1860  assert(HR.isSorted());
1861  //assume rows are sorted (this is to be ensured by corresponding GetHessian functions)
1862  //preallocate J to JL.Size+JR.Size
1863  //perform merging of two sorted arrays
1864  //resize to correct size
1865  Sparse::Row::MergeSortedRows(right.GetValue()*multJ,JL,left.GetValue()*multJ,JR,J);
1866  assert(J.isSorted());
1867  //preallocate H to HL.Size+HR.Size+JL.Size*JR.Size
1868  //merge sorted
1869  Sparse::HessianRow::MergeJacobianHessian(2.0*multH,JL,JR,right.GetValue()*multH,HL,left.GetValue()*multH,HR,H);
1870  assert(H.isSorted());
1871  }
1872  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
1873  };
1874 
1875  template<class A, class B>
1876  class division_expression : public shell_expression<division_expression<A,B> >
1877  {
1878  const A & left;
1879  const B & right;
1880  INMOST_DATA_REAL_TYPE value, reciprocal_rval;
1881  public:
1882  division_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
1883  {
1884  INMOST_DATA_REAL_TYPE lval = left.GetValue();
1885  INMOST_DATA_REAL_TYPE rval = right.GetValue();
1886  if( rval )
1887  {
1888  reciprocal_rval = 1.0 / rval;
1889  value = lval * reciprocal_rval;
1890  }
1891  else
1892  {
1893  reciprocal_rval = 0;
1894  value = 0;
1895  }
1896  }
1897  division_expression(const division_expression & other) : left(other.left), right(other.right), value(other.value), reciprocal_rval(other.reciprocal_rval) {}
1898  division_expression(const division_expression & other, const A & pleft, const B & pright) :
1899  left(pleft), right(pright), value(other.value), reciprocal_rval(other.reciprocal_rval) {}
1900 
1901  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1902  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1903  {
1904  left.GetJacobian(mult * reciprocal_rval,r);
1905  right.GetJacobian(- mult * value * reciprocal_rval,r);
1906  }
1907  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1908  {
1909  left.GetJacobian(mult * reciprocal_rval,r);
1910  right.GetJacobian(- mult * value * reciprocal_rval,r);
1911  }
1912  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1913  {
1914  throw NotImplemented;
1915  Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
1916  Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
1917  left.GetHessian(multJ,JL,multH,HL); //retrieve jacobian row and hessian matrix of the left expression
1918  right.GetHessian(multJ,JR,multH,HR); //retrieve jacobian row and hessian matrix of the right expression
1919  //assume rows are sorted (this is to be ensured by corresponding GetHessian functions)
1920  //preallocate J to JL.Size+JR.Size
1921  //perform merging of two sorted arrays
1922  //resize to correct size
1923  Sparse::Row::MergeSortedRows(reciprocal_rval,JL,-value * reciprocal_rval,JR,J);
1924  //preallocate H to HL.Size+HR.Size+JL.Size*JR.Size
1925  //merge sorted
1926  Sparse::HessianRow::MergeJacobianHessian(2.0,JL,JR,reciprocal_rval,HL,2*left.GetValue()*multH,HR,H);
1927  }
1928  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return reciprocal_rval ? left.GetCount() + right.GetCount() : 0; }
1929  };
1930 
1931  template<class A, class B>
1932  class addition_expression : public shell_expression<addition_expression<A,B> >
1933  {
1934  const A & left;
1935  const B & right;
1936  INMOST_DATA_REAL_TYPE value;
1937  public:
1938  addition_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
1939  {
1940  value = left.GetValue() + right.GetValue();
1941  }
1942  addition_expression(const addition_expression & other)
1943  : left(other.left), right(other.right), value(other.value) {}
1944  addition_expression(const addition_expression & other, const A & pleft, const B & pright)
1945  : left(pleft), right(pright), value(other.value) {}
1946  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1947  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1948  {
1949  left.GetJacobian(mult,r);
1950  right.GetJacobian(mult,r);
1951  }
1952  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1953  {
1954  left.GetJacobian(mult,r);
1955  right.GetJacobian(mult,r);
1956  }
1957  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
1958  {
1959  Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
1960  Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
1961  left.GetHessian(1,JL,1,HL); //retrieve jacobian row and hessian matrix of the left expression
1962  assert(JL.isSorted());
1963  assert(HL.isSorted());
1964  right.GetHessian(1,JR,1,HR); //retrieve jacobian row and hessian matrix of the right expression
1965  assert(JR.isSorted());
1966  assert(HR.isSorted());
1967  Sparse::Row::MergeSortedRows(multJ,JL,multJ,JR,J);
1968  assert(J.isSorted());
1969  Sparse::HessianRow::MergeSortedRows(multH,HL,multH,HR,H);
1970  assert(H.isSorted());
1971  }
1972  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
1973  };
1974 
1975  template<class A, class B>
1976  class subtraction_expression : public shell_expression<subtraction_expression<A,B> >
1977  {
1978  const A & left;
1979  const B & right;
1980  INMOST_DATA_REAL_TYPE value;
1981  public:
1982  subtraction_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
1983  {
1984  value = left.GetValue() - right.GetValue();
1985  }
1986  subtraction_expression(const subtraction_expression & other)
1987  : left(other.left), right(other.right),value(other.value) {}
1988  subtraction_expression(const subtraction_expression & other, const A & pleft, const B & pright)
1989  : left(pleft), right(pright),value(other.value) {}
1990  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
1991  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
1992  {
1993  left.GetJacobian(mult,r);
1994  right.GetJacobian(-mult,r);
1995  }
1996  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
1997  {
1998  left.GetJacobian(mult,r);
1999  right.GetJacobian(-mult,r);
2000  }
2001  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2002  {
2003  //(F-G)'' = (F'-G')' = (F''-G'')
2004  Sparse::Row JL, JR; //temporary jacobian rows from left and right expressions
2005  Sparse::HessianRow HL, HR; //temporary hessian rows form left and right expressions
2006  left.GetHessian(1,JL,1,HL); //retrieve jacobian row and hessian matrix of the left expression
2007  right.GetHessian(1,JR,1,HR); //retrieve jacobian row and hessian matrix of the right expression
2008  Sparse::Row::MergeSortedRows(multJ,JL,-multJ,JR,J);
2009  Sparse::HessianRow::MergeSortedRows(multH,HL,-multH,HR,H);
2010  }
2011  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
2012  };
2013 
2014  template<class A, class B>
2015  class pow_expression : public shell_expression<pow_expression<A,B> >
2016  {
2017  const A & left;
2018  const B & right;
2019  INMOST_DATA_REAL_TYPE value, ldmult, rdmult;
2020  public:
2021  pow_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
2022  {
2023  INMOST_DATA_REAL_TYPE lval = left.GetValue();
2024  INMOST_DATA_REAL_TYPE rval = right.GetValue();
2025  value = ::pow(lval,rval);
2026  if( lval != 0 )
2027  {
2028  ldmult = value * rval / lval;
2029  rdmult = value * ::log(lval);
2030  }
2031  else
2032  {
2033  ldmult = 0;
2034  rdmult = 0;
2035  }
2036  }
2037  pow_expression(const pow_expression & other)
2038  :left(other.left), right(other.right), value(other.value),
2039  ldmult(other.ldmult), rdmult(other.rdmult) {}
2040  pow_expression(const pow_expression & other, const A & pleft, const B & pright)
2041  :left(pleft), right(pright), value(other.value),
2042  ldmult(other.ldmult), rdmult(other.rdmult) {}
2043  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2044  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2045  {
2046  left.GetJacobian(mult*ldmult,r);
2047  right.GetJacobian(mult*rdmult,r);
2048  }
2049  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2050  {
2051  left.GetJacobian(mult*ldmult,r);
2052  right.GetJacobian(mult*rdmult,r);
2053  }
2054  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2055  {
2056  throw NotImplemented;
2057  }
2058  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
2059  };
2060 
2061 
2062  template<class A, class B>
2063  class atan2_expression : public shell_expression<atan2_expression<A,B> >
2064  {
2065  const A & left;
2066  const B & right;
2067  INMOST_DATA_REAL_TYPE value, ldmult, rdmult;
2068  public:
2069  atan2_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
2070  {
2071  INMOST_DATA_REAL_TYPE lval = left.GetValue();
2072  INMOST_DATA_REAL_TYPE rval = right.GetValue();
2073  value = ::atan2(lval,rval);
2074  ldmult = rval/(rval*rval+lval*lval);
2075  rdmult = -lval/(rval*rval+lval*lval);
2076  }
2077  atan2_expression(const atan2_expression & other)
2078  :left(other.left), right(other.right), value(other.value),
2079  ldmult(other.ldmult), rdmult(other.rdmult) {}
2080  atan2_expression(const atan2_expression & other, const A & pleft, const B & pright)
2081  :left(pleft), right(pright), value(other.value),
2082  ldmult(other.ldmult), rdmult(other.rdmult) {}
2083  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2084  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2085  {
2086  left.GetJacobian(mult*ldmult,r);
2087  right.GetJacobian(mult*rdmult,r);
2088  }
2089  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2090  {
2091  left.GetJacobian(mult*ldmult,r);
2092  right.GetJacobian(mult*rdmult,r);
2093  }
2094  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2095  {
2096  throw NotImplemented;
2097  }
2098  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount() + right.GetCount(); }
2099  };
2100 
2101  template<class A>
2102  class pow_const_expression : public shell_expression<pow_const_expression<A> >
2103  {
2104  const A & left;
2105  INMOST_DATA_REAL_TYPE value, ldmult, ldmult2;
2106  public:
2107  pow_const_expression(const shell_expression<A> & pleft, INMOST_DATA_REAL_TYPE pright) : left(pleft)
2108  {
2109  INMOST_DATA_REAL_TYPE lval = left.GetValue();
2110  value = ::pow(lval,pright);
2111  if( lval != 0 )
2112  {
2113  ldmult = value * pright / lval;
2114  ldmult2 = ldmult * (pright-1) / lval;
2115  }
2116  else
2117  ldmult = ldmult2 = 0;
2118  }
2119  pow_const_expression(const pow_const_expression & other)
2120  :left(other.left), value(other.value), ldmult(other.ldmult) {}
2121  pow_const_expression(const pow_const_expression & other, const A & pleft)
2122  :left(pleft), value(other.value), ldmult(other.ldmult) {}
2123  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2124  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2125  {
2126  left.GetJacobian(mult*ldmult,r);
2127  }
2128  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2129  {
2130  left.GetJacobian(mult*ldmult,r);
2131  }
2132  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2133  {
2134  //(F(G))'' = (F'(G)G')' = (F''(G)G'+F'(G)G'')
2135  //(g(x)^n)'' = n g(x)^(n - 2) (g(x) g''(x) + (n - 1) g'(x)^2)
2136  Sparse::Row JL;
2137  Sparse::HessianRow HL;
2138  left.GetHessian(1,JL,1,HL); //retrieve jacobian row and hessian matrix of the left expression
2139  Sparse::HessianRow::MergeJacobianHessian(multH*ldmult2,JL,JL,multH*ldmult,HL,H);
2140  for(Sparse::Row::iterator it = JL.Begin(); it != JL.End(); ++it) it->second *= ldmult*multJ;
2141  (void)J;
2142  }
2143  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return left.GetCount(); }
2144  };
2145 
2146  template<class A>
2147  class const_pow_expression : public shell_expression<const_pow_expression<A> >
2148  {
2149  const A & right;
2150  INMOST_DATA_REAL_TYPE value, rdmult;
2151  public:
2152  const_pow_expression(INMOST_DATA_REAL_TYPE pleft, const shell_expression<A> & pright) : right(pright)
2153  {
2154  INMOST_DATA_REAL_TYPE rval = right.GetValue();
2155  value = ::pow(pleft,rval);
2156  if( pleft != 0 )
2157  rdmult = value * ::log(pleft);
2158  else
2159  rdmult = 0;
2160  }
2161  const_pow_expression(const const_pow_expression & other)
2162  :right(other.right), value(other.value), rdmult(other.rdmult) {}
2163  const_pow_expression(const const_pow_expression & other, const A & pright)
2164  :right(pright), value(other.value), rdmult(other.rdmult) {}
2165  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2166  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2167  {
2168  right.GetJacobian(mult*rdmult,r);
2169  }
2170  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2171  {
2172  right.GetJacobian(mult*rdmult,r);
2173  }
2174  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2175  {
2176  throw NotImplemented;
2177  }
2178  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return right.GetCount(); }
2179  };
2180 
2181  template<class A, class B, class C>
2182  class condition_expression : public shell_expression<condition_expression<A,B,C> >
2183  {
2184  const A & cond;
2185  const B & left;
2186  const C & right;
2187  INMOST_DATA_REAL_TYPE value, cond_value;
2188  public:
2189  condition_expression(const shell_expression<A> & pcond, const shell_expression<B> & pleft, const shell_expression<C> & pright) : cond(pcond), left(pleft), right(pright)
2190  {
2191  cond_value = cond.GetValue();
2192  value = cond_value >= 0.0 ? left.GetValue() : right.GetValue();
2193  }
2194  condition_expression(const condition_expression & other)
2195  :cond(other.cond), left(other.left), right(other.right),
2196  value(other.value), cond_value(other.cond_value) {}
2197  condition_expression(const condition_expression & other, const A & pcond, const B & pleft, const C & pright)
2198  :cond(pcond), left(pleft), right(pright),
2199  value(other.value), cond_value(other.cond_value) {}
2200  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2201  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2202  {
2203  if( cond_value >= 0.0 )
2204  left.GetJacobian(mult,r);
2205  else
2206  right.GetJacobian(mult,r);
2207  }
2208  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2209  {
2210  if( cond_value >= 0.0 )
2211  left.GetJacobian(mult,r);
2212  else
2213  right.GetJacobian(mult,r);
2214  }
2215  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2216  {
2217  if( cond_value >= 0.0 )
2218  left.GetHessian(multJ,J,multH,H);
2219  else
2220  right.GetHessian(multJ,J,multH,H);
2221  }
2222  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return cond_value >= 0.0 ? left.GetCount() : right.GetCount(); }
2223  };
2224 
2225  template<class A, class B>
2226  class branch_expression : public shell_expression<branch_expression<A,B> >
2227  {
2228  bool cond;
2229  const A & left;
2230  const B & right;
2231  INMOST_DATA_REAL_TYPE value;
2232  public:
2233  branch_expression(const shell_expression<A> & pleft, const shell_expression<B> & pright) : left(pleft), right(pright)
2234  {
2235  value = 0;
2236  }
2237  branch_expression(bool pcond, const shell_expression<A> & pleft, const shell_expression<B> & pright) : cond(pcond), left(pleft), right(pright)
2238  {
2239  value = cond ? left.GetValue() : right.GetValue();
2240  }
2241  branch_expression(const branch_expression & other)
2242  :cond(other.cond), left(other.left), right(other.right),
2243  value(other.value) {}
2244  branch_expression(const branch_expression & other, bool pcond, const A & pleft, const B & pright)
2245  :cond(pcond), left(pleft), right(pright),
2246  value(other.value) {}
2247  __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; }
2248  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2249  {
2250  if( cond )
2251  left.GetJacobian(mult,r);
2252  else
2253  right.GetJacobian(mult,r);
2254  }
2255  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2256  {
2257  if( cond )
2258  left.GetJacobian(mult,r);
2259  else
2260  right.GetJacobian(mult,r);
2261  }
2262  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2263  {
2264  if( cond )
2265  left.GetHessian(multJ,J,multH,H);
2266  else
2267  right.GetHessian(multJ,J,multH,H);
2268  }
2269  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return cond ? left.GetCount() : right.GetCount(); }
2270  void SetCondition(bool _cond)
2271  {
2272  cond = _cond;
2273  value = cond ? left.GetValue() : right.GetValue();
2274  }
2275  };
2276 
2277 
2278  template<class A>
2279  class function_expression : public shell_expression< function_expression<A> >
2280  {
2281  const A &arg;
2282  INMOST_DATA_REAL_TYPE value, dmult, ddmult;
2283  public:
2284  function_expression(const shell_expression<A> &_arg)
2285  : arg(_arg), value(1), dmult(0) {}
2286 
2287  function_expression(const shell_expression<A> &_arg, INMOST_DATA_REAL_TYPE pvalue, INMOST_DATA_REAL_TYPE pdmult,
2288  INMOST_DATA_REAL_TYPE pddmult = 0)
2289  : arg(_arg), value(pvalue), dmult(pdmult), ddmult(pddmult) {}
2290 
2291  function_expression(const function_expression &other)
2292  : arg(other.arg), value(other.value), dmult(other.dmult), ddmult(other.ddmult) {}
2293  function_expression(const function_expression &other, const A & parg)
2294  : arg(parg), value(other.value), dmult(other.dmult), ddmult(other.ddmult) {}
2295 
2296  function_expression &operator=(function_expression const &b)
2297  {
2298  arg = b.arg;
2299  value = b.value;
2300  dmult = b.dmult;
2301  ddmult = b.ddmult;
2302  return *this;
2303  }
2304  __INLINE INMOST_DATA_REAL_TYPE GetValue() const {return value;}
2305  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const
2306  {
2307  arg.GetJacobian(mult*dmult,r);
2308  }
2309  __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const
2310  {
2311  arg.GetJacobian(mult*dmult,r);
2312  }
2313  __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const
2314  {
2315  arg.GetHessian(multJ*dmult,J,multH*ddmult,H);
2316  }
2317  __INLINE INMOST_DATA_ENUM_TYPE GetCount() const { return arg.GetCount(); }
2318  void SetFunctionValue(INMOST_DATA_REAL_TYPE val) {value = val;}
2319  void SetFunctionDerivative(INMOST_DATA_REAL_TYPE val) {dmult = val;}
2320  };
2321 
2322 
2323  class keyval_table
2324  {
2325  std::string name;
2326  INMOST_DATA_REAL_TYPE * vals;
2327  INMOST_DATA_REAL_TYPE * args;
2328  INMOST_DATA_ENUM_TYPE size;
2329  INMOST_DATA_ENUM_TYPE binary_search(INMOST_DATA_REAL_TYPE arg) const
2330  {
2331  int l = 0, r = static_cast<int>(size)-1, mid = 0;
2332  while (r >= l)
2333  {
2334  mid = (l + r) / 2;
2335  if (args[mid] > arg) r = mid - 1;
2336  else if (args[mid] < arg) l = mid + 1;
2337  else break;
2338  }
2339  mid = (l + r) / 2;
2340  if (mid > static_cast<int>(size - 2)) mid = static_cast<int>(size - 2);
2341  return static_cast<INMOST_DATA_ENUM_TYPE>(mid);
2342  }
2343  public:
2344  INMOST_DATA_REAL_TYPE GetValue(INMOST_DATA_REAL_TYPE arg) const
2345  {
2346  if (arg < args[0]) return vals[0];
2347  INMOST_DATA_ENUM_TYPE i = binary_search(arg);
2348  return vals[i] + (vals[i + 1] - vals[i]) * (arg - args[i]) / (args[i + 1] - args[i]);
2349  }
2350  INMOST_DATA_REAL_TYPE GetDerivative(INMOST_DATA_REAL_TYPE arg) const
2351  {
2352  if (arg < args[0]) return 0.0;
2353  INMOST_DATA_ENUM_TYPE i = binary_search(arg);
2354  return (vals[i + 1] - vals[i]) / (args[i + 1] - args[i]);
2355  }
2356  std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> GetBoth(INMOST_DATA_REAL_TYPE arg) const
2357  {
2358  if (arg < args[0]) return std::make_pair(vals[0], 0.0);
2359  INMOST_DATA_ENUM_TYPE i = binary_search(arg);
2360  INMOST_DATA_REAL_TYPE der = (vals[i + 1] - vals[i]) / (args[i + 1] - args[i]);
2361  return std::make_pair(vals[i] + der * (arg - args[i]), der);
2362  }
2363  keyval_table() :name(""), vals(NULL), args(NULL), size(0) {}
2364  keyval_table(std::string _name, const INMOST_DATA_REAL_TYPE * _args, const INMOST_DATA_REAL_TYPE * _vals, INMOST_DATA_ENUM_TYPE _size)
2365  {
2366  name = _name;
2367  size = _size;
2368  vals = new INMOST_DATA_REAL_TYPE[size];
2369  memcpy(vals,_vals,sizeof(INMOST_DATA_REAL_TYPE)*size);
2370  //std::copy(_vals,_vals+size,vals);
2371  args = new INMOST_DATA_REAL_TYPE[size];
2372  memcpy(args,_args,sizeof(INMOST_DATA_REAL_TYPE)*size);
2373  //std::copy(_args,_args+size,args);
2374  }
2375  keyval_table(const keyval_table & other)
2376  {
2377  name = other.name;
2378  size = other.size;
2379  vals = new INMOST_DATA_REAL_TYPE[size];
2380  memcpy(vals,other.vals,sizeof(INMOST_DATA_REAL_TYPE)*size);
2381  //std::copy(other.vals,other.vals+size,vals);
2382  args = new INMOST_DATA_REAL_TYPE[size];
2383  memcpy(args,other.args,sizeof(INMOST_DATA_REAL_TYPE)*size);
2384  //std::copy(other.args,other.args+size,args);
2385  }
2386  keyval_table & operator = (keyval_table const & other)
2387  {
2388  Clear();
2389  name = other.name;
2390  size = other.size;
2391  vals = new INMOST_DATA_REAL_TYPE[size];
2392  memcpy(vals,other.vals,sizeof(INMOST_DATA_REAL_TYPE)*size);
2393  //std::copy(other.vals,other.vals+size,vals);
2394  args = new INMOST_DATA_REAL_TYPE[size];
2395  memcpy(args,other.args,sizeof(INMOST_DATA_REAL_TYPE)*size);
2396  //std::copy(other.args,other.args+size,args);
2397  return * this;
2398  }
2399  ~keyval_table()
2400  {
2401  Clear();
2402  }
2403  void Clear()
2404  {
2405  name = "";
2406  if( args ) {delete [] args; args = NULL;}
2407  if( vals ) {delete [] vals; vals = NULL;}
2408  size = 0;
2409  }
2410  bool Empty() const
2411  {
2412  return size == 0;
2413  }
2414  const INMOST_DATA_REAL_TYPE * GetTableArguments() const {return args;}
2415  const INMOST_DATA_REAL_TYPE * GetTableValues() const {return vals;}
2416  INMOST_DATA_REAL_TYPE * GetTableArguments() {return args;}
2417  INMOST_DATA_REAL_TYPE * GetTableValues() {return vals;}
2418  INMOST_DATA_ENUM_TYPE GetSize() const {return size;}
2419  };
2420 
2421  typedef multivar_expression variable;
2422  typedef hessian_multivar_expression hessian_variable;
2423  typedef var_expression unknown;
2424 }
2425 
2426 
2427 __INLINE bool check_nans(INMOST_DATA_REAL_TYPE val) {return val != val;}
2428 __INLINE bool check_nans(INMOST::var_expression const & e) {return e.check_nans();}
2429 __INLINE bool check_nans(INMOST::multivar_expression const & e) {return e.check_nans();}
2430 __INLINE bool check_nans(INMOST::multivar_expression_reference const & e) {return e.check_nans();}
2431 __INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return __isinf__(val);}
2432 __INLINE bool check_infs(INMOST::var_expression const & e) {return e.check_infs();}
2433 __INLINE bool check_infs(INMOST::multivar_expression const & e) {return e.check_infs();}
2434 __INLINE bool check_infs(INMOST::multivar_expression_reference const & e) {return e.check_infs();}
2435 __INLINE bool check_nans_infs(INMOST_DATA_REAL_TYPE val) {return check_nans(val) || check_infs(val);}
2436 __INLINE bool check_nans_infs(INMOST::var_expression const & e) {return e.check_nans() || e.check_infs();}
2437 __INLINE bool check_nans_infs(INMOST::multivar_expression const & e) {return e.check_nans() || e.check_infs();}
2438 __INLINE bool check_nans_infs(INMOST::multivar_expression_reference const & e) {return e.check_nans() || e.check_infs();}
2439 __INLINE INMOST_DATA_ENUM_TYPE GetCount(INMOST_DATA_REAL_TYPE val) { (void)val; return 0; }
2440 __INLINE INMOST_DATA_ENUM_TYPE GetCount(INMOST_DATA_CPLX_TYPE val) { (void)val; return 0; }
2441 __INLINE INMOST_DATA_ENUM_TYPE GetCount(const INMOST::basic_expression& expr) { return expr.GetCount(); }
2442 
2443 template<class A, class B, class C> __INLINE INMOST::condition_expression<A,B,C> condition(INMOST::shell_expression<A> const & control, INMOST::shell_expression<B> const & if_ge_zero, INMOST::shell_expression<C> const & if_lt_zero) { return INMOST::condition_expression<A,B,C>(control,if_ge_zero,if_lt_zero); }
2444 __INLINE INMOST_DATA_REAL_TYPE condition(INMOST_DATA_REAL_TYPE control, INMOST_DATA_REAL_TYPE if_ge_zero, INMOST_DATA_REAL_TYPE if_lt_zero) {return control >= 0.0 ? if_ge_zero : if_lt_zero;}
2445 template<class A> __INLINE INMOST::unary_minus_expression<A> operator-(INMOST::shell_expression<A> const & Arg) { return INMOST::unary_minus_expression<A>(Arg); }
2446 template<class A> __INLINE INMOST::unary_plus_expression<A> operator+(INMOST::shell_expression<A> const & Arg) { return INMOST::unary_plus_expression<A>(Arg); }
2447 template<class A> __INLINE INMOST::abs_expression<A> fabs(INMOST::shell_expression<A> const & Arg) { return INMOST::abs_expression<A>(Arg); }
2448 template<class A> __INLINE INMOST::exp_expression<A> exp(INMOST::shell_expression<A> const & Arg) { return INMOST::exp_expression<A> (Arg); }
2449 template<class A> __INLINE INMOST::log_expression<A> log(INMOST::shell_expression<A> const & Arg) { return INMOST::log_expression<A> (Arg); }
2450 template<class A> __INLINE INMOST::sin_expression<A> sin(INMOST::shell_expression<A> const & Arg) { return INMOST::sin_expression<A> (Arg ); }
2451 template<class A> __INLINE INMOST::cos_expression<A> cos(INMOST::shell_expression<A> const & Arg) { return INMOST::cos_expression<A> (Arg); }
2452 template<class A> __INLINE INMOST::sqrt_expression<A> sqrt(INMOST::shell_expression<A> const & Arg) { return INMOST::sqrt_expression<A> (Arg); }
2453 template<class A> __INLINE INMOST::variation_multiplication_expression<A> variation(INMOST::shell_expression<A> const & Arg,INMOST_DATA_REAL_TYPE Mult) {return INMOST::variation_multiplication_expression<A>(Arg,Mult);}
2454 __INLINE INMOST_DATA_REAL_TYPE variation(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE) {return Arg;}
2455 template<class A> __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST::shell_expression<A> const & Arg) { return Arg.GetValue(); }
2456  __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE Arg) {return Arg;}
2457  __INLINE INMOST_DATA_CPLX_TYPE get_value(INMOST_DATA_CPLX_TYPE Arg) { return Arg; }
2458  __INLINE INMOST_DATA_CPLX_TYPE get_value(std::complex<INMOST::var_expression> const& Arg) { return INMOST_DATA_CPLX_TYPE(Arg.real().GetValue(), Arg.imag().GetValue()); }
2459  __INLINE INMOST_DATA_CPLX_TYPE get_value(std::complex<INMOST::multivar_expression> const& Arg) { return INMOST_DATA_CPLX_TYPE(Arg.real().GetValue(), Arg.imag().GetValue()); }
2460  __INLINE INMOST_DATA_CPLX_TYPE get_value(std::complex<INMOST::hessian_multivar_expression> const& Arg) { return INMOST_DATA_CPLX_TYPE(Arg.real().GetValue(), Arg.imag().GetValue()); }
2461  __INLINE const INMOST::var_expression& conj(const INMOST::var_expression& Arg) { return Arg; }
2462  __INLINE const INMOST::multivar_expression& conj(const INMOST::multivar_expression& Arg) { return Arg; }
2463  __INLINE const INMOST::multivar_expression_reference& conj(const INMOST::multivar_expression_reference& Arg) { return Arg; }
2464  __INLINE const INMOST::hessian_multivar_expression& conj(const INMOST::hessian_multivar_expression& Arg) { return Arg; }
2465  __INLINE const INMOST::hessian_multivar_expression_reference& conj(const INMOST::hessian_multivar_expression_reference& Arg) { return Arg; }
2466  __INLINE void set_value(INMOST::var_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg.SetValue(Val); }
2467  __INLINE void set_value(INMOST::multivar_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg.SetValue(Val); }
2468  __INLINE void set_value(INMOST::multivar_expression_reference & Arg, INMOST_DATA_REAL_TYPE Val) {Arg.SetValue(Val); }
2469  __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;}
2470  __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, const INMOST::var_expression & Val) {Arg = Val.GetValue();}
2471  __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression & Val) {Arg = Val.GetValue();}
2472  __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val.GetValue();}
2473  __INLINE void set_value(INMOST::multivar_expression & Arg, const INMOST::var_expression & Val) {Arg.SetValue(Val.GetValue()); }
2474  __INLINE void set_value(INMOST::multivar_expression & Arg, const INMOST::multivar_expression & Val) {Arg.SetValue(Val.GetValue()); }
2475  __INLINE void set_value(INMOST::multivar_expression & Arg, const INMOST::multivar_expression_reference & Val) {Arg.SetValue(Val.GetValue()); }
2476  __INLINE void set_value(INMOST::multivar_expression_reference & Arg, const INMOST::var_expression & Val) {Arg.SetValue(Val.GetValue()); }
2477  __INLINE void set_value(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression & Val) {Arg.SetValue(Val.GetValue()); }
2478  __INLINE void set_value(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression_reference & Val) {Arg.SetValue(Val.GetValue()); }
2479  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val;}
2480  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val;}
2481  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::var_expression & Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val.GetValue();}
2482  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::multivar_expression & Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val.GetValue();}
2483  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::multivar_expression_reference & Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val.GetValue();}
2484  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val.GetValue(); }
2485  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::hessian_multivar_expression_reference & Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val.GetValue(); }
2486  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = (INMOST_DATA_REAL_TYPE)Val;}
2487  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;}
2488  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::var_expression & Val) {Arg = Val.GetValue();}
2489  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression & Val) {Arg = Val.GetValue();}
2490  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val.GetValue();}
2491  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = Val.GetValue(); }
2492  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::hessian_multivar_expression_reference & Val) {Arg = Val.GetValue(); }
2493 // __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::value_reference& Val) { Arg = Val.GetValue(); }
2494  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST_DATA_CPLX_TYPE& Val) { Arg = Val; }
2495  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, INMOST_DATA_INTEGER_TYPE Val) { Arg = (INMOST_DATA_REAL_TYPE)Val; }
2496  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, INMOST_DATA_REAL_TYPE Val) { Arg = Val; }
2497  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::var_expression& Val) { Arg = Val.GetValue(); }
2498  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::multivar_expression& Val) { Arg = Val.GetValue(); }
2499  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::multivar_expression_reference& Val) { Arg = Val.GetValue(); }
2500  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::hessian_multivar_expression& Val) { Arg = Val.GetValue(); }
2501  __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::hessian_multivar_expression_reference& Val) { Arg = Val.GetValue(); }
2502 // __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, const INMOST::value_reference& Val) { Arg = Val.GetValue(); }
2503  __INLINE void assign(INMOST::var_expression & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2504  __INLINE void assign(INMOST::var_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val; }
2505  __INLINE void assign(INMOST::var_expression & Arg, const INMOST::var_expression & Val) {Arg = Val; }
2506  __INLINE void assign(INMOST::multivar_expression & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2507  __INLINE void assign(INMOST::multivar_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val; }
2508 // __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::var_expression & Val) {Arg = Val; }
2509 // __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::multivar_expression & Val) {Arg = Val; }
2510 // __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val; }
2511 // __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = Val; }
2512 // __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::hessian_multivar_expression_reference & Val) {Arg = Val; }
2513  __INLINE void assign(INMOST::multivar_expression_reference & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2514  __INLINE void assign(INMOST::multivar_expression_reference & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val; }
2515 // __INLINE void assign(INMOST::multivar_expression_reference & Arg, const INMOST::var_expression & Val) {Arg = Val; }
2516 // __INLINE void assign(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression & Val) {Arg = Val; }
2517 // __INLINE void assign(INMOST::multivar_expression_reference & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val; }
2518 // __INLINE void assign(INMOST::multivar_expression_reference & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = Val; }
2519  __INLINE void assign(INMOST::hessian_multivar_expression & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2520  __INLINE void assign(INMOST::hessian_multivar_expression & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val; }
2521 // __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::var_expression & Val) {Arg = Val; }
2522 // __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::multivar_expression & Val) {Arg = Val; }
2523 // __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val; }
2524 // __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = Val; }
2525 // __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::hessian_multivar_expression_reference & Val) {Arg = Val; }
2526  __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val; }
2527  __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val; }
2528 // __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::var_expression & Val) {Arg = Val; }
2529 // __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::multivar_expression & Val) {Arg = Val; }
2530 // __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::multivar_expression_reference & Val) {Arg = Val; }
2531 // __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::hessian_multivar_expression & Val) {Arg = Val; }
2532 // __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::hessian_multivar_expression_reference & Val) {Arg = Val; }
2533 template<class A> __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, const INMOST::shell_expression<A> & Val) {Arg = (INMOST_DATA_REAL_TYPE)Val.GetValue();}
2534 template<class A> __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, const INMOST::shell_expression<A> & Val) {Arg = Val.GetValue();}
2535 template<class A> __INLINE void assign(INMOST::multivar_expression & Arg, const INMOST::shell_expression<A> & Val) {Arg = Val;}
2536 //template<class A> __INLINE void assign(INMOST::value_reference& Arg, const INMOST::shell_expression<A>& Val) {Arg.SetValue(Val.GetValue());}
2537 template<class A> __INLINE void assign(INMOST::multivar_expression_reference & Arg, const INMOST::shell_expression<A> & Val) {Arg = Val;}
2538 template<class A> __INLINE void assign(INMOST::hessian_multivar_expression & Arg, const INMOST::shell_expression<A> & Val) {Arg = Val;}
2539 template<class A> __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, const INMOST::shell_expression<A> & Val) {Arg = Val;}
2540 // __INLINE void assign(INMOST::value_reference& Arg, INMOST_DATA_REAL_TYPE Val) { Arg = (INMOST_DATA_REAL_TYPE)Val; }
2541 #if defined(USE_FP64)
2542  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, float Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val; }
2543  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2544  __INLINE void assign(float& Arg, float Val) {Arg = (float)Val; }
2545  __INLINE void assign(INMOST::var_expression & Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2546  __INLINE void assign(INMOST::multivar_expression & Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2547  __INLINE void assign(INMOST::multivar_expression_reference & Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2548 // __INLINE void assign(INMOST::value_reference& Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2549  __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, float Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2550 #else //USE_FP64
2551  __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, double Val) {Arg = (INMOST_DATA_INTEGER_TYPE)Val; }
2552  __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2553  __INLINE void assign(double& Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2554  __INLINE void assign(INMOST::var_expression & Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2555  __INLINE void assign(INMOST::multivar_expression & Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2556  __INLINE void assign(INMOST::multivar_expression_reference & Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2557  __INLINE void assign(INMOST::hessian_multivar_expression_reference & Arg, double Val) {Arg = (INMOST_DATA_REAL_TYPE)Val; }
2558 #endif //USE_FP64
2559 template<class A> __INLINE INMOST::soft_abs_expression<A> soft_fabs(INMOST::shell_expression<A> const & Arg, INMOST_DATA_REAL_TYPE tol = 0) { return INMOST::soft_abs_expression<A>(Arg,tol); }
2560 __INLINE INMOST_DATA_REAL_TYPE soft_fabs(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol = 0) {return ::sqrt(Arg*Arg+tol*tol);}
2561 template<class A> __INLINE INMOST::soft_sign_expression<A> soft_sign(INMOST::shell_expression<A> const & Arg, INMOST_DATA_REAL_TYPE tol = 0) { return INMOST::soft_sign_expression<A>(Arg,tol); }
2562 __INLINE INMOST_DATA_REAL_TYPE soft_sign(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol = 0) {return Arg/::sqrt(Arg*Arg+tol*tol);}
2563 template<class A, class B> __INLINE INMOST::multiplication_expression<A, B> operator*(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::multiplication_expression<A, B> (Left, Right); }
2564 template<class A, class B> __INLINE INMOST::division_expression<A, B> operator/(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::division_expression<A, B> (Left, Right); }
2565 template<class A, class B> __INLINE INMOST::addition_expression<A, B> operator+(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::addition_expression<A, B> (Left, Right); }
2566 template<class A, class B> __INLINE INMOST::subtraction_expression<A, B> operator-(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::subtraction_expression<A, B> (Left, Right); }
2567 template<class A, class B> __INLINE INMOST::pow_expression<A, B> pow(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::pow_expression<A, B> (Left, Right); }
2568 template<class A, class B> __INLINE INMOST::atan2_expression<A, B> atan2(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) { return INMOST::atan2_expression<A, B> (Left, Right); }
2569 template<class A, class B> __INLINE INMOST::soft_max_expression<A, B> soft_max(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right ,INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_max_expression<A, B> (Left, Right,tol); }
2570 template<class A> __INLINE INMOST::soft_max_const_expression<A> soft_max(INMOST::shell_expression<A> const& Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_max_const_expression<A>(Left, Right, tol); }
2571 template<class B> __INLINE INMOST::soft_max_const_expression<B> soft_max(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const& Right, INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_max_const_expression<B>(Right, Left, tol); }
2572  __INLINE INMOST_DATA_REAL_TYPE soft_max(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol = 0.0) {return 0.5*(Left+Right+::sqrt((Left-Right)*(Left-Right)+tol*tol));}
2573 template<class A, class B> __INLINE INMOST::soft_min_expression<A, B> soft_min(INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right ,INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_min_expression<A, B> (Left, Right,tol); }
2574 template<class A> __INLINE INMOST::soft_min_const_expression<A> soft_min(INMOST::shell_expression<A> const& Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_min_const_expression<A>(Left, Right, tol); }
2575 template<class B> __INLINE INMOST::soft_min_const_expression<B> soft_min(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const& Right, INMOST_DATA_REAL_TYPE tol = 0.0) { return INMOST::soft_min_const_expression<B>(Right, Left, tol); }
2576  __INLINE INMOST_DATA_REAL_TYPE soft_min(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol = 0.0) {return 0.5*(Left+Right-::sqrt((Left-Right)*(Left-Right)+tol*tol));}
2577 template<class B> __INLINE INMOST::const_pow_expression<B> pow(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) { return INMOST::const_pow_expression<B> (Left, Right); }
2578 template<class A> __INLINE INMOST::pow_const_expression<A> pow(INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::pow_const_expression<A> (Left, Right); }
2579 template<class B> __INLINE INMOST::const_multiplication_expression<B> operator*(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) { return INMOST::const_multiplication_expression<B>(Right,Left); }
2580 template<class A> __INLINE INMOST::const_multiplication_expression<A> operator*(INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_multiplication_expression<A>(Left,Right); }
2581 template<class B> __INLINE INMOST::reciprocal_expression<B> operator/(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) { return INMOST::reciprocal_expression<B>(Right,Left); }
2582 template<class A> __INLINE INMOST::const_division_expression<A> operator/(INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_division_expression<A>(Left, Right); }
2583 template<class B> __INLINE INMOST::const_addition_expression<B> operator+(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) { return INMOST::const_addition_expression<B>(Right,Left); }
2584 template<class A> __INLINE INMOST::const_addition_expression<A> operator+(INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_addition_expression<A>(Left,Right); }
2585 template<class B> __INLINE INMOST::const_subtraction_expression<B> operator-(INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) { return INMOST::const_subtraction_expression<B>(Right, Left); }
2586 template<class A> __INLINE INMOST::const_addition_expression<A> operator-(INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) { return INMOST::const_addition_expression<A>(Left, -Right); }
2587 template<class A, class B> __INLINE bool operator == (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() == Right.GetValue();}
2588 template<class A, class B> __INLINE bool operator != (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() != Right.GetValue();}
2589 template<class A, class B> __INLINE bool operator < (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() < Right.GetValue();}
2590 template<class A, class B> __INLINE bool operator > (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() > Right.GetValue();}
2591 template<class A, class B> __INLINE bool operator <= (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() <= Right.GetValue();}
2592 template<class A, class B> __INLINE bool operator >= (INMOST::shell_expression<A> const & Left, INMOST::shell_expression<B> const & Right) {return Left.GetValue() >= Right.GetValue();}
2593 template<class A> __INLINE bool operator == (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() == Right;}
2594 template<class A> __INLINE bool operator != (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() != Right;}
2595 template<class A> __INLINE bool operator < (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() < Right;}
2596 template<class A> __INLINE bool operator > (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() > Right;}
2597 template<class A> __INLINE bool operator <= (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() <= Right;}
2598 template<class A> __INLINE bool operator >= (INMOST::shell_expression<A> const & Left, INMOST_DATA_REAL_TYPE Right) {return Left.GetValue() >= Right;}
2599 template<class B> __INLINE bool operator == (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left == Right.GetValue();}
2600 template<class B> __INLINE bool operator != (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left != Right.GetValue();}
2601 template<class B> __INLINE bool operator < (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left < Right.GetValue();}
2602 template<class B> __INLINE bool operator > (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left > Right.GetValue();}
2603 template<class B> __INLINE bool operator <= (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left <= Right.GetValue();}
2604 template<class B> __INLINE bool operator >= (INMOST_DATA_REAL_TYPE Left, INMOST::shell_expression<B> const & Right) {return Left >= Right.GetValue();}
2605 template<class A> __INLINE INMOST::function_expression<A> get_table(INMOST::shell_expression<A> const & Arg, const INMOST::keyval_table & Table)
2606 {
2607  std::pair<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> both = Table.GetBoth(Arg.GetValue());
2608  return INMOST::function_expression<A>(Arg,both.first,both.second);
2609 }
2610 __INLINE INMOST_DATA_REAL_TYPE get_table(INMOST_DATA_REAL_TYPE Arg, const INMOST::keyval_table & Table) {return Table.GetValue(Arg);}
2611 #else //USE_AUTODIFF
2612 __INLINE INMOST_DATA_ENUM_TYPE GetCount(INMOST_DATA_REAL_TYPE val) { (void)val; return 0; }
2613 __INLINE INMOST_DATA_ENUM_TYPE GetCount(INMOST_DATA_CPLX_TYPE val) { (void)val; return 0; }
2614 __INLINE bool check_nans(INMOST_DATA_REAL_TYPE val) {return val != val;}
2615 __INLINE bool check_infs(INMOST_DATA_REAL_TYPE val) {return __isinf__(val);}
2616 __INLINE bool check_nans_infs(INMOST_DATA_REAL_TYPE val) {return check_nans(val) || check_infs(val);}
2617 __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = Val;}
2618 __INLINE void assign(INMOST_DATA_INTEGER_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = static_cast<INMOST_DATA_INTEGER_TYPE>(Val);}
2619 __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_INTEGER_TYPE Val) {Arg = static_cast<INMOST_DATA_REAL_TYPE>(Val);}
2620 __INLINE void assign(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;}
2621 __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, INMOST_DATA_REAL_TYPE Val) { Arg = Val; }
2622 __INLINE void assign(INMOST_DATA_CPLX_TYPE& Arg, INMOST_DATA_CPLX_TYPE Val) { Arg = Val; }
2623 __INLINE void set_value(INMOST_DATA_REAL_TYPE & Arg, INMOST_DATA_REAL_TYPE Val) {Arg = Val;}
2624 __INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE Arg) {return Arg;}
2625 __INLINE INMOST_DATA_CPLX_TYPE get_value(INMOST_DATA_CPLX_TYPE Arg) { return Arg; }
2626 __INLINE INMOST_DATA_REAL_TYPE variation(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE) {return Arg;}
2627 __INLINE INMOST_DATA_REAL_TYPE soft_fabs(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol = 0) {return ::sqrt(Arg*Arg+tol*tol);}
2628 __INLINE INMOST_DATA_REAL_TYPE soft_sign(INMOST_DATA_REAL_TYPE Arg, INMOST_DATA_REAL_TYPE tol = 0) {return Arg/::sqrt(Arg*Arg+tol*tol);}
2629 __INLINE INMOST_DATA_REAL_TYPE soft_max(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right+::sqrt((Left-Right)*(Left-Right)+tol*tol));}
2630 __INLINE INMOST_DATA_REAL_TYPE soft_min(INMOST_DATA_REAL_TYPE Left, INMOST_DATA_REAL_TYPE Right, INMOST_DATA_REAL_TYPE tol) {return 0.5*(Left+Right-::sqrt((Left-Right)*(Left-Right)+tol*tol));}
2631 #endif //USE_AUTODIFF
2632 #endif //INMOST_AUTODIFF_ETEXPR_H_INCLUDED
INMOST::Sparse::HessianRow
Class to store the compressed symmetric matrix of a hessian row.
Definition: inmost_sparse.h:334
INMOST::Sparse::HessianRow::Clear
void Clear()
Clear all data of the current row.
Definition: inmost_sparse.h:402
INMOST::Sparse::HessianRow::MergeSortedRows
static void MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const HessianRow &left, INMOST_DATA_REAL_TYPE beta, const HessianRow &right, HessianRow &output)
output = alpha * left + beta *right
INMOST::Sparse::HessianRow::Size
INMOST_DATA_ENUM_TYPE Size() const
The size of the sparse row, i.e. the total number of nonzero elements.
Definition: inmost_sparse.h:405
INMOST::Sparse::HessianRow::index
struct INMOST::Sparse::HessianRow::hessian_index_s index
Entry of the sparse matrix row.
INMOST::Sparse::RowMerger
This class may be used to sum multiple sparse rows.
Definition: inmost_sparse.h:638
INMOST::Sparse::RowMerger::AddRow
void AddRow(INMOST_DATA_REAL_TYPE coef, const Row &r)
Add a row with a coefficient into non-empty linked list.
INMOST::Sparse::Row
Class to store the sparse matrix row.
Definition: inmost_sparse.h:162
INMOST::Sparse::Row::Swap
void Swap(Row &other)
Exchange all the data with another row.
Definition: inmost_sparse.h:251
INMOST::Sparse::Row::get_safe
INMOST_DATA_REAL_TYPE get_safe(INMOST_DATA_ENUM_TYPE i) const
Finds and returns value with specified index.
Definition: inmost_sparse.h:242
INMOST::Sparse::Row::Clear
void Clear()
Clear all data of the current row.
Definition: inmost_sparse.h:248
INMOST::Sparse::Row::GetValue
INMOST_DATA_REAL_TYPE & GetValue(INMOST_DATA_ENUM_TYPE k)
Retrieve a value corresponding to certain position in the array of pairs of index and value.
Definition: inmost_sparse.h:263
INMOST::Sparse::Row::Print
void Print(double eps=-1, std::ostream &sout=std::cout) const
Output all entries of the row.
Definition: inmost_sparse.h:310
INMOST::Sparse::Row::iterator
Entries::iterator iterator
Iterator over pairs of index and value.
Definition: inmost_sparse.h:190
INMOST::Sparse::Row::isSorted
bool isSorted() const
Check whether the row is sorted.
INMOST::Sparse::Row::Resize
void Resize(INMOST_DATA_ENUM_TYPE size)
Resize row to specified size.
Definition: inmost_sparse.h:308
INMOST::Sparse::Row::MergeSortedRows
static void MergeSortedRows(INMOST_DATA_REAL_TYPE alpha, const Row &left, INMOST_DATA_REAL_TYPE beta, const Row &right, Row &output)
Add up two rows.
INMOST::Sparse::Row::End
iterator End()
An iterator pointing behind the last position in the array of pairs of index and value.
Definition: inmost_sparse.h:277
INMOST::Sparse::Row::Size
INMOST_DATA_ENUM_TYPE Size() const
The size of the sparse row, i.e. the total number of nonzero elements.
Definition: inmost_sparse.h:253
INMOST::Sparse::Row::const_iterator
Entries::const_iterator const_iterator
Iterator over constant pairs of index and value.
Definition: inmost_sparse.h:191
INMOST::Sparse::Row::Begin
iterator Begin()
Retrive interval of nonzeroes.
Definition: inmost_sparse.h:275
INMOST::Sparse::Row::GetIndex
INMOST_DATA_ENUM_TYPE & GetIndex(INMOST_DATA_ENUM_TYPE k)
Retrieve an index corresponding to certain position in the array of pairs of index and value.
Definition: inmost_sparse.h:259
INMOST::Sparse::Row::Push
void Push(INMOST_DATA_ENUM_TYPE ind, INMOST_DATA_REAL_TYPE val)
Push specified element into sparse row.
Definition: inmost_sparse.h:301
INMOST::abs_expression
Definition: inmost_expression.h:1407
INMOST::addition_expression
Definition: inmost_expression.h:1933
INMOST::atan2_expression
Definition: inmost_expression.h:2064
INMOST::basic_expression
Definition: inmost_expression.h:37
INMOST::branch_expression
Definition: inmost_expression.h:2227
INMOST::condition_expression
Definition: inmost_expression.h:2183
INMOST::const_addition_expression
Definition: inmost_expression.h:1263
INMOST::const_division_expression
Definition: inmost_expression.h:1234
INMOST::const_expression
Definition: inmost_expression.h:67
INMOST::const_multiplication_expression
Definition: inmost_expression.h:1177
INMOST::const_pow_expression
Definition: inmost_expression.h:2148
INMOST::const_subtraction_expression
Definition: inmost_expression.h:1291
INMOST::cos_expression
Definition: inmost_expression.h:1546
INMOST::division_expression
Definition: inmost_expression.h:1877
INMOST::exp_expression
Definition: inmost_expression.h:1440
INMOST::function_expression
Definition: inmost_expression.h:2280
INMOST::hessian_multivar_expression_reference
Definition: inmost_expression.h:956
INMOST::hessian_multivar_expression_reference::GetValue
__INLINE INMOST_DATA_REAL_TYPE GetValue() const
Retrieve value.
Definition: inmost_expression.h:971
INMOST::hessian_multivar_expression_reference::hessian_multivar_expression_reference
hessian_multivar_expression_reference(hessian_multivar_expression &other)
Copy constructor from hessian_multivar_expression, sets links to the same reference of value and entr...
Definition: inmost_expression.h:968
INMOST::hessian_multivar_expression_reference::GetJacobian
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger &r) const
Retrieve derivatives with multiplier into Sparse::RowMerger structure.
Definition: inmost_expression.h:975
INMOST::hessian_multivar_expression_reference::hessian_multivar_expression_reference
hessian_multivar_expression_reference(INMOST_DATA_REAL_TYPE &_value, Sparse::Row *_entries, Sparse::HessianRow *_hentries)
Constructor, set links to the provided value and entries.
Definition: inmost_expression.h:962
INMOST::hessian_multivar_expression_reference::GetJacobian
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row &r) const
Retrieve derivatives with multiplier into Sparse::Row structure.
Definition: inmost_expression.h:982
INMOST::hessian_multivar_expression_reference::hessian_multivar_expression_reference
hessian_multivar_expression_reference(const hessian_multivar_expression_reference &other)
Copy constructor, sets links to the same reference of value and entries.
Definition: inmost_expression.h:965
INMOST::hessian_multivar_expression_reference::SetValue
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val)
Set value without changing derivatives.
Definition: inmost_expression.h:973
INMOST::hessian_multivar_expression
A class that represents a variable with multiple first order and second order variations.
Definition: inmost_expression.h:408
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue)
Sets value and no first or second order variations.
Definition: inmost_expression.h:422
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, Sparse::Row &pentries, Sparse::HessianRow &phentries)
Sets value and all first and second order variations.
Definition: inmost_expression.h:431
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, INMOST_DATA_ENUM_TYPE pindex, INMOST_DATA_REAL_TYPE pdmult=1.0)
Sets value and it's variation index.
Definition: inmost_expression.h:434
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(const basic_expression &expr)
Evaluates argument expression to get the value and all variations of variable.
Definition: inmost_expression.h:440
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression()
Sets zero value and no first or second order variations.
Definition: inmost_expression.h:420
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(INMOST_DATA_REAL_TYPE pvalue, Sparse::Row &pentries)
Copy value and all first variations from variable, no second order variations.
Definition: inmost_expression.h:428
INMOST::hessian_multivar_expression::hessian_multivar_expression
hessian_multivar_expression(const hessian_multivar_expression &other)
Copy value and all first and second order variations from hessian_variable.
Definition: inmost_expression.h:424
INMOST::keyval_table
Definition: inmost_expression.h:2324
INMOST::log_expression
Definition: inmost_expression.h:1474
INMOST::multiplication_expression
Definition: inmost_expression.h:1826
INMOST::multivar_expression_reference
Definition: inmost_expression.h:642
INMOST::multivar_expression_reference::GetJacobian
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger &r) const
Retrieve derivatives with multiplier into Sparse::RowMerger structure.
Definition: inmost_expression.h:662
INMOST::multivar_expression_reference::multivar_expression_reference
multivar_expression_reference(INMOST_DATA_REAL_TYPE &_value, Sparse::Row *_entries)
Constructor, set links to the provided value and entries.
Definition: inmost_expression.h:649
INMOST::multivar_expression_reference::multivar_expression_reference
multivar_expression_reference()
Default constructor.
Definition: inmost_expression.h:647
INMOST::multivar_expression_reference::GetJacobian
__INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row &r) const
Retrieve derivatives with multiplier into Sparse::Row structure.
Definition: inmost_expression.h:672
INMOST::multivar_expression_reference::GetValue
__INLINE INMOST_DATA_REAL_TYPE GetValue() const
Retrieve value.
Definition: inmost_expression.h:658
INMOST::multivar_expression_reference::multivar_expression_reference
multivar_expression_reference(multivar_expression &other)
Copy constructor from multivar_expression, sets links to the same reference of value and entries.
Definition: inmost_expression.h:655
INMOST::multivar_expression_reference::SetValue
__INLINE void SetValue(INMOST_DATA_REAL_TYPE val)
Set value without changing derivatives.
Definition: inmost_expression.h:660
INMOST::multivar_expression_reference::multivar_expression_reference
multivar_expression_reference(const multivar_expression_reference &other)
Copy constructor, sets links to the same reference of value and entries.
Definition: inmost_expression.h:652
INMOST::multivar_expression
A class that represents a variable with multiple first order variations.
Definition: inmost_expression.h:131
INMOST::multivar_expression::RetrieveSize
static INMOST_DATA_ENUM_TYPE RetrieveSize(const Sparse::Row::entry *v)
Number of entries used.
Definition: inmost_expression.h:387
INMOST::multivar_expression::Retrieve
INMOST_DATA_ENUM_TYPE Retrieve(const Sparse::Row::entry *v)
Retrieve variable from array of entries.
Definition: inmost_expression.h:372
INMOST::multivar_expression::Record
INMOST_DATA_ENUM_TYPE Record(Sparse::Row::entry *v) const
Write variable into array of entries.
Definition: inmost_expression.h:349
INMOST::multivar_expression::RecordSize
INMOST_DATA_ENUM_TYPE RecordSize() const
Number of entries required to record the variable.
Definition: inmost_expression.h:364
INMOST::pow_const_expression
Definition: inmost_expression.h:2103
INMOST::pow_expression
Definition: inmost_expression.h:2016
INMOST::reciprocal_expression
(c/x)' = -c dx / (x*x) (c/x)'' = 2 c dx dx / (x*x*x)
Definition: inmost_expression.h:1321
INMOST::shell_expression
Definition: inmost_expression.h:52
INMOST::sin_expression
Definition: inmost_expression.h:1510
INMOST::soft_abs_expression
Definition: inmost_expression.h:1618
INMOST::soft_max_const_expression
Definition: inmost_expression.h:1725
INMOST::soft_max_expression
Definition: inmost_expression.h:1682
INMOST::soft_min_const_expression
Definition: inmost_expression.h:1797
INMOST::soft_min_expression
Definition: inmost_expression.h:1754
INMOST::soft_sign_expression
Definition: inmost_expression.h:1649
INMOST::sqrt_expression
Definition: inmost_expression.h:1580
INMOST::subtraction_expression
Definition: inmost_expression.h:1977
INMOST::thread_private
Definition: container.hpp:1569
INMOST::unary_minus_expression
Definition: inmost_expression.h:1357
INMOST::unary_plus_expression
Definition: inmost_expression.h:1382
INMOST::value_reference
Definition: inmost_expression.h:871
INMOST::value_reference::value_reference
value_reference(const value_reference &other)
Copy constructor, sets links to the same reference of value and entries.
Definition: inmost_expression.h:880
INMOST::value_reference::value_reference
value_reference(INMOST_DATA_REAL_TYPE &_value)
Constructor, set links to the provided value and entries.
Definition: inmost_expression.h:877
INMOST::value_reference::value_reference
value_reference()
Default constructor.
Definition: inmost_expression.h:875
INMOST::var_expression
Definition: inmost_expression.h:86
INMOST::variation_multiplication_expression
Definition: inmost_expression.h:1205
INMOST::Sparse::Row::entry_s
Entry of the sparse matrix row.
Definition: inmost_sparse.h:167
INMOST::Sparse::Row::entry_s::second
INMOST_DATA_REAL_TYPE second
the real value of the row element.
Definition: inmost_sparse.h:169
INMOST::Sparse::Row::entry_s::first
INMOST_DATA_ENUM_TYPE first
the column number of the row element.
Definition: inmost_sparse.h:168