Scippy

SCIP

Solving Constraint Integer Programs

exprinterpret_cppad.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2014 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not email to scip@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file exprinterpret_cppad.cpp
17  * @brief methods to interpret (evaluate) an expression tree "fast" using CppAD
18  * @ingroup EXPRINTS
19  * @author Stefan Vigerske
20  */
21 
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
23 
24 #include "scip/def.h"
25 #include "blockmemshell/memory.h"
26 #include "nlpi/pub_expr.h"
27 #include "nlpi/exprinterpret.h"
28 
29 #include <cmath>
30 #include <vector>
31 using std::vector;
32 
33 /* Turn off lint warning "747: Significant prototype coercion" and "732: Loss of sign".
34  * The first warning is generated for expressions like t[0], where t is a vector, since 0 is an integer constant, but a
35  * size_t is expected (usually long unsigned). The second is generated for expressions like t[n], where n is an
36  * integer. Both code pieces are likely to be correct. It seems to be impossible to inhibit these messages for
37  * vector<*>::operator[] only. */
38 /*lint --e{747,732}*/
39 
40 /* Turn off lint info "1702 operator '...' is both an ordinary function 'CppAD::operator...' and a member function 'CppAD::SCIPInterval::operator...'.
41  * However, the functions have different signatures (the CppAD working on double, the SCIPInterval member
42  * function working on SCIPInterval's.
43  */
44 /*lint --e{1702}*/
45 
46 /* defining NO_CPPAD_USER_ATOMIC disables the use of our own implementation of derivaties of power operators
47  * via CppAD's user-atomic function feature
48  * our customized implementation should give better results (tighter intervals) for the interval data type
49  */
50 /* #define NO_CPPAD_USER_ATOMIC */
51 
52 /** sign of a value (-1 or +1)
53  *
54  * 0.0 has sign +1
55  */
56 #define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
57 
58 /* in order to use intervals as operands in CppAD,
59  * we need to include the intervalarithext.h very early and require the interval operations to be in the CppAD namespace */
60 #define SCIPInterval_NAMESPACE CppAD
61 #include "nlpi/intervalarithext.h"
62 
63 SCIP_Real CppAD::SCIPInterval::infinity = SCIP_DEFAULT_INFINITY;
64 using CppAD::SCIPInterval;
65 
66 /* CppAD needs to know a fixed upper bound on the number of threads at compile time.
67  * It is wise to set it to a power of 2, so that if the tape id overflows, it is likely to start at 0 again, which avoids difficult to debug errors.
68  */
69 #ifndef CPPAD_MAX_NUM_THREADS
70 #ifndef NPARASCIP
71 #define CPPAD_MAX_NUM_THREADS 64
72 #else
73 #define CPPAD_MAX_NUM_THREADS 1
74 #endif
75 #endif
76 
77 #include <cppad/cppad.hpp>
78 #include <cppad/error_handler.hpp>
79 
80 /* CppAD is not thread-safe by itself, but uses some static datastructures
81  * To run it in a multithreading environment, a special CppAD memory allocator that is aware of the multiple threads has to be used.
82  * This allocator requires to know the number of threads and a thread number for each thread.
83  * To implement this, we follow the team_pthread example of CppAD, which uses pthread's thread-specific data management.
84  */
85 #ifndef NPARASCIP
86 #include <pthread.h>
87 
88 /* workaround error message regarding missing implementation of tanh during initialization of static variables (see cppad/local/erf.hpp) */
89 namespace CppAD
90 {
91 template <> SCIPInterval erf_template(
92  const SCIPInterval &x
93  )
94 {
95  CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
96  return SCIPInterval();
97 }
98 template <> AD<SCIPInterval> erf_template(
99  const AD<SCIPInterval> &x
100  )
101 {
102  CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
103  return AD<SCIPInterval>();
104 }
105 }
106 
107 /** mutex for locking in pthread case */
108 static pthread_mutex_t cppadmutex = PTHREAD_MUTEX_INITIALIZER;
109 
110 /** key for accessing thread specific information */
111 static pthread_key_t thread_specific_key;
112 
113 /** currently registered number of threads */
114 static size_t ncurthreads = 0;
115 
116 /** CppAD callback function that indicates whether we are running in parallel mode */
117 static
118 bool in_parallel(void)
119 {
120  return ncurthreads > 1;
121 }
122 
123 /** CppAD callback function that returns the number of the current thread
124  *
125  * assigns a new number to the thread if new
126  */
127 static
128 size_t thread_num(void)
129 {
130  size_t threadnum;
131  void* specific;
132 
133  specific = pthread_getspecific(thread_specific_key);
134 
135  /* if no data for this thread yet, then assign a new thread number to the current thread
136  * we store the thread number incremented by one, to distinguish the absence of data (=0) from existing data
137  */
138  if( specific == NULL )
139  {
140  pthread_mutex_lock(&cppadmutex);
141 
142  SCIPdebugMessage("Assigning thread number %lu to thread %p.\n", (long unsigned int)ncurthreads, (void*)pthread_self());
143 
144  pthread_setspecific(thread_specific_key, (void*)(ncurthreads + 1));
145 
146  threadnum = ncurthreads;
147 
148  ++ncurthreads;
149 
150  pthread_mutex_unlock(&cppadmutex);
151 
152  assert(pthread_getspecific(thread_specific_key) != NULL);
153  assert((size_t)pthread_getspecific(thread_specific_key) == threadnum + 1);
154  }
155  else
156  {
157  threadnum = (size_t)(specific) - 1;
158  }
159 
160  assert(threadnum < ncurthreads);
161 
162  return threadnum;
163 }
164 
165 /** sets up CppAD's datastructures for running in multithreading mode
166  *
167  * It must be called once before multithreading is started.
168  */
169 static
170 char init_parallel(void)
171 {
172  pthread_key_create(&thread_specific_key, NULL);
173 
174  CppAD::thread_alloc::parallel_setup(CPPAD_MAX_NUM_THREADS, in_parallel, thread_num);
175  CppAD::parallel_ad<double>();
176  CppAD::parallel_ad<SCIPInterval>();
177 
178  return 0;
179 }
180 
181 /** a dummy variable that can is initialized to the result of init_parallel
182  *
183  * The purpose is to make sure that init_parallel() is called before any multithreading is started.
184  */
185 static char init_parallel_return = init_parallel();
186 
187 #endif // NPARASCIP
188 
189 /** definition of CondExpOp for SCIPInterval (required by CppAD) */
190 inline
192  enum CppAD::CompareOp cop,
193  const SCIPInterval& left,
194  const SCIPInterval& right,
195  const SCIPInterval& trueCase,
196  const SCIPInterval& falseCase)
197 { /*lint --e{715}*/
198  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
199  "SCIPInterval CondExpOp(...)",
200  "Error: cannot use CondExp with an interval type"
201  );
202 
203  return SCIPInterval();
204 }
205 
206 /** another function that returns whether two intervals are the same (required by CppAD) */
207 inline
209  const SCIPInterval& x, /**< first operand */
210  const SCIPInterval& y /**< second operand */
211  )
212 {
213  return x == y;
214 }
215 
216 /** another function required by CppAD */
217 inline
219  const SCIPInterval& x /**< operand */
220  )
221 { /*lint --e{715}*/
222  return true;
223 }
224 
225 /** returns whether the interval equals [0,0] */
226 inline
228  const SCIPInterval& x /**< operand */
229  )
230 {
231  return (x == 0.0);
232 }
233 
234 /** returns whether the interval equals [1,1] */
235 inline
237  const SCIPInterval& x /**< operand */
238  )
239 {
240  return (x == 1.0);
241 }
242 
243 /** yet another function that checks whether two intervals are equal */
244 inline
246  const SCIPInterval& x, /**< first operand */
247  const SCIPInterval& y /**< second operand */
248  )
249 {
250  return (x == y);
251 }
252 
253 /** greater than zero not defined for intervals */
254 inline
256  const SCIPInterval& x /**< operand */
257  )
258 { /*lint --e{715}*/
259  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
260  "GreaterThanZero(x)",
261  "Error: cannot use GreaterThanZero with interval"
262  );
263 
264  return false;
265 }
266 
267 /** greater than or equal zero not defined for intervals */
268 inline
270  const SCIPInterval& x /**< operand */
271  )
272 { /*lint --e{715}*/
273  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__ ,
274  "GreaterThanOrZero(x)",
275  "Error: cannot use GreaterThanOrZero with interval"
276  );
277 
278  return false;
279 }
280 
281 /** less than not defined for intervals */
282 inline
284  const SCIPInterval& x /**< operand */
285  )
286 { /*lint --e{715}*/
287  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
288  "LessThanZero(x)",
289  "Error: cannot use LessThanZero with interval"
290  );
291 
292  return false;
293 }
294 
295 /** less than or equal not defined for intervals */
296 inline
298  const SCIPInterval& x /**< operand */
299  )
300 { /*lint --e{715}*/
301  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
302  "LessThanOrZero(x)",
303  "Error: cannot use LessThanOrZero with interval"
304  );
305 
306  return false;
307 }
308 
309 /** conversion to integers not defined for intervals */
310 inline
312  const SCIPInterval& x /**< operand */
313  )
314 { /*lint --e{715}*/
315  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
316  "Integer(x)",
317  "Error: cannot use Integer with interval"
318  );
319 
320  return 0;
321 }
322 
323 /** printing of an interval (required by CppAD) */
324 inline
325 std::ostream& operator<<(std::ostream& out, const SCIP_INTERVAL& x)
326 {
327  out << '[' << x.inf << ',' << x.sup << ']';
328  return out;
329 }
330 
331 using CppAD::AD;
332 
333 /** expression interpreter */
334 struct SCIP_ExprInt
335 {
336  BMS_BLKMEM* blkmem; /**< block memory data structure */
337 };
338 
339 /** expression specific interpreter data */
340 struct SCIP_ExprIntData
341 {
342 public:
343  /** constructor */
344  SCIP_ExprIntData()
345  : val(0.0), need_retape(true), int_need_retape(true), need_retape_always(false), blkmem(NULL), root(NULL)
346  { }
347 
348  /** destructor */
349  ~SCIP_ExprIntData()
350  { }/*lint --e{1540}*/
351 
352  vector< AD<double> > X; /**< vector of dependent variables */
353  vector< AD<double> > Y; /**< result vector */
354  CppAD::ADFun<double> f; /**< the function to evaluate as CppAD object */
355 
356  vector<double> x; /**< current values of dependent variables */
357  double val; /**< current function value */
358  bool need_retape; /**< will retaping be required for the next point evaluation? */
359 
360  vector< AD<SCIPInterval> > int_X; /**< interval vector of dependent variables */
361  vector< AD<SCIPInterval> > int_Y; /**< interval result vector */
362  CppAD::ADFun<SCIPInterval> int_f; /**< the function to evaluate on intervals as CppAD object */
363 
364  vector<SCIPInterval> int_x; /**< current interval values of dependent variables */
365  SCIPInterval int_val; /**< current interval function value */
366  bool int_need_retape; /**< will retaping be required for the next interval evaluation? */
367 
368  bool need_retape_always; /**< will retaping be always required? */
369 
370  BMS_BLKMEM* blkmem; /**< block memory used to allocate expresstion tree */
371  SCIP_EXPR* root; /**< copy of expression tree; @todo we should not need to make a copy */
372 };
373 
374 
375 #ifndef NO_CPPAD_USER_ATOMIC
376 
377 /** computes sparsity of jacobian for a univariate function during a forward sweep
378  *
379  * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R.
380  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
381  */
382 static
383 bool univariate_for_sparse_jac(
384  size_t q, /**< number of columns in R */
385  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
386  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
387  )
388 {
389  assert(r.size() == q);
390  assert(s.size() == q);
391 
392  s = r;
393 
394  return true;
395 }
396 
397 /** Computes sparsity of jacobian during a reverse sweep
398  *
399  * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x).
400  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
401  */
402 static
403 bool univariate_rev_sparse_jac(
404  size_t q, /**< number of rows in R */
405  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
406  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
407  )
408 {
409  assert(r.size() == q);
410  assert(s.size() == q);
411 
412  s = r;
413 
414  return true;
415 }
416 
417 /** computes sparsity of hessian during a reverse sweep
418  *
419  * Assume V(x) = (g(f(x)))'' R with f(x) = x^p for a function g:R->R and a matrix R.
420  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
421  */
422 static
423 bool univariate_rev_sparse_hes(
424  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
425  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
426  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
427  size_t q, /**< number of columns in R, U, and V */
428  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
429  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
430  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
431  )
432 { /*lint --e{439,715}*/ /* @todo take vx into account */
433  assert(r.size() == q);
434  assert(s.size() == 1);
435  assert(t.size() == 1);
436  assert(u.size() == q);
437  assert(v.size() == q);
438 
439  // T(x) = g'(f(x)) * f'(x) = S * f'(x), and f' is not identically 0
440  t[0] = s[0];
441 
442  // V(x) = g''(f(x)) f'(x) f'(x) R + g'(f(x)) f''(x) R466
443  // = f'(x) U + S f''(x) R, with f'(x) and f''(x) not identically 0
444  v = u;
445  if( s[0] )
446  for( size_t j = 0; j < q; ++j )
447  if( r[j] )
448  v[j] = true;
449 
450  return true;
451 }
452 
453 
454 /** Automatic differentiation of x -> x^p, p>=2 integer, as CppAD user-atomic function.
455  *
456  * This class implements forward and reverse operations for the function x -> x^p for use within CppAD.
457  * While CppAD would implement integer powers as a recursion of multiplications, we still use pow functions as they allow us to avoid overestimation in interval arithmetics.
458  *
459  * @todo treat the exponent as a (variable) argument to the function, with the assumption that we never differentiate w.r.t. it (this should make the approach threadsafe again)
460  */
461 template<class Type>
462 class atomic_posintpower : public CppAD::atomic_base<Type>
463 {
464 public:
466  : CppAD::atomic_base<Type>("posintpower"),
467  exponent(0)
468  {
469  /* indicate that we want to use bool-based sparsity pattern */
470  this->option(CppAD::atomic_base<Type>::bool_sparsity_enum);
471  }
472 
473 private:
474  /** exponent value for next call to forward or reverse */
475  int exponent;
476 
477  /** stores exponent value corresponding to next call to forward or reverse
478  *
479  * how is this supposed to be threadsafe? (we use only one global instantiation of this class)
480  */
481  virtual void set_id(size_t id)
482  {
483  exponent = (int) id;
484  }
485 
486  /** forward sweep of positive integer power
487  *
488  * Given the taylor coefficients for x, we have to compute the taylor coefficients for f(x),
489  * that is, given tx = (x, x', x'', ...), we compute the coefficients ty = (y, y', y'', ...)
490  * in the taylor expansion of f(x) = x^p.
491  * Thus, y = x^p
492  * = tx[0]^p,
493  * y' = p * x^(p-1) * x'
494  * = p * tx[0]^(p-1) * tx[1],
495  * y'' = 1/2 * p * (p-1) * x^(p-2) * x'^2 + p * x^(p-1) * x''
496  * = 1/2 * p * (p-1) * tx[0]^(p-2) * tx[1]^2 + p * tx[0]^(p-1) * tx[2]
497  */
498  bool forward(
499  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
500  size_t p, /**< highest order Taylor coefficient that we are evaluating */
501  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
502  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
503  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
504  CppAD::vector<Type>& ty /**< vector to store taylor coefficients of y */
505  )
506  {
507  assert(exponent > 1);
508  assert(tx.size() >= p+1);
509  assert(ty.size() >= p+1);
510  assert(q <= p);
511 
512  if( vx.size() > 0 )
513  {
514  assert(vx.size() == 1);
515  assert(vy.size() == 1);
516  assert(p == 0);
517 
518  vy[0] = vx[0];
519  }
520 
521  if( q == 0 /* q <= 0 && 0 <= p */ )
522  {
523  ty[0] = CppAD::pow(tx[0], exponent);
524  }
525 
526  if( q <= 1 && 1 <= p )
527  {
528  ty[1] = CppAD::pow(tx[0], exponent-1) * tx[1];
529  ty[1] *= double(exponent);
530  }
531 
532  if( q <= 2 && 2 <= p )
533  {
534  if( exponent > 2 )
535  {
536  // ty[2] = 1/2 * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1] * tx[1] + exponent * pow(tx[0], exponent-1) * tx[2];
537  ty[2] = CppAD::pow(tx[0], exponent-2) * tx[1] * tx[1];
538  ty[2] *= (exponent-1) / 2.0;
539  ty[2] += CppAD::pow(tx[0], exponent-1) * tx[2];
540  ty[2] *= exponent;
541  }
542  else
543  {
544  assert(exponent == 2);
545  // ty[2] = 1/2 * exponent * tx[1] * tx[1] + exponent * tx[0] * tx[2];
546  ty[2] = tx[1] * tx[1] + 2.0 * tx[0] * tx[2];
547  }
548  }
549 
550  /* higher order derivatives not implemented */
551  if( p > 2 )
552  return false;
553 
554  return true;
555  }
556 
557  /** reverse sweep of positive integer power
558  *
559  * Assume y(x) is a function of the taylor coefficients of f(x) = x^p for x, i.e.,
560  * y(x) = [ x^p, p * x^(p-1) * x', p * (p-1) * x^(p-2) * x'^2 + p * x^(p-1) * x'', ... ].
561  * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
562  * where x^[l] is the l'th taylor coefficient (x, x', x'', ...) and h(x) = g(y(x)) for some function g:R^k -> R.
563  * That is, we have to compute
564  *\f$
565  * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
566  * = \sum_{i=0}^k (\partial g / \partial y_i) * (\partial y_i / \partial x^[l])
567  * = \sum_{i=0}^k py[i] * (\partial y_i / \partial x^[l])
568  * \f$
569  *
570  * For k = 0, this means
571  *\f$
572  * px[0] = py[0] * (\partial y_0 / \partial x^[0])
573  * = py[0] * (\partial x^p / \partial x)
574  * = py[0] * p * tx[0]^(p-1)
575  *\f$
576  *
577  * For k = 1, this means
578  * \f$
579  * px[0] = py[0] * (\partial y_0 / \partial x^[0]) + py[1] * (\partial y_1 / \partial x^[0])
580  * = py[0] * (\partial x^p / \partial x) + py[1] * (\partial (p * x^(p-1) * x') / \partial x)
581  * = py[0] * p * tx[0]^(p-1) + py[1] * p * (p-1) * tx[0]^(p-2) * tx[1]
582  * px[1] = py[0] * (\partial y_0 / \partial x^[1]) + py[1] * (\partial y_1 / \partial x^[1])
583  * = py[0] * (\partial x^p / \partial x') + py[1] * (\partial (p * x^(p-1) x') / \partial x')
584  * = py[0] * 0 + py[1] * p * tx[0]^(p-1)
585  * \f$
586  */
587  bool reverse(
588  size_t p, /**< highest order Taylor coefficient that we are evaluating */
589  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
590  const CppAD::vector<Type>& ty, /**< values for taylor coefficients of y */
591  CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
592  const CppAD::vector<Type>& py /**< values for partial derivatives of g(x) w.r.t. y */
593  )
594  { /*lint --e{715}*/
595  assert(exponent > 1);
596  assert(px.size() >= p+1);
597  assert(py.size() >= p+1);
598  assert(tx.size() >= p+1);
599 
600  switch( p )
601  {
602  case 0:
603  // px[0] = py[0] * exponent * pow(tx[0], exponent-1);
604  px[0] = py[0] * CppAD::pow(tx[0], exponent-1);
605  px[0] *= exponent;
606  break;
607 
608  case 1:
609  // px[0] = py[0] * exponent * pow(tx[0], exponent-1) + py[1] * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1];
610  px[0] = py[1] * tx[1] * CppAD::pow(tx[0], exponent-2);
611  px[0] *= exponent-1;
612  px[0] += py[0] * CppAD::pow(tx[0], exponent-1);
613  px[0] *= exponent;
614  // px[1] = py[1] * exponent * pow(tx[0], exponent-1);
615  px[1] = py[1] * CppAD::pow(tx[0], exponent-1);
616  px[1] *= exponent;
617  break;
618 
619  default:
620  return false;
621  }
622 
623  return true;
624  }
625 
626  using CppAD::atomic_base<Type>::for_sparse_jac;
627 
628  /** computes sparsity of jacobian during a forward sweep
629  *
630  * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R.
631  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
632  */
633  bool for_sparse_jac(
634  size_t q, /**< number of columns in R */
635  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
636  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
637  )
638  {
639  return univariate_for_sparse_jac(q, r, s);
640  }
641 
642  using CppAD::atomic_base<Type>::rev_sparse_jac;
643 
644  /** computes sparsity of jacobian during a reverse sweep
645  *
646  * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x).
647  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
648  */
649  bool rev_sparse_jac(
650  size_t q, /**< number of rows in R */
651  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
652  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
653  )
654  {
655  return univariate_rev_sparse_jac(q, r, s);
656  }
657 
658  using CppAD::atomic_base<Type>::rev_sparse_hes;
659 
660  /** computes sparsity of hessian during a reverse sweep
661  *
662  * Assume V(x) = (g(f(x)))'' R with f(x) = x^p for a function g:R->R and a matrix R.
663  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
664  */
665  bool rev_sparse_hes(
666  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
667  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
668  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
669  size_t q, /**< number of columns in R, U, and V */
670  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
671  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
672  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
673  )
674  {
675  return univariate_rev_sparse_hes(vx, s, t, q, r, u, v);
676  }
677 };
678 
679 /** power function with natural exponents */
680 template<class Type>
681 static
682 void posintpower(
683  const vector<Type>& in, /**< vector which first argument is base */
684  vector<Type>& out, /**< vector where to store result in first argument */
685  size_t exponent /**< exponent */
686  )
687 {
689  pip(in, out, exponent);
690 }
691 
692 #else
693 
694 /** power function with natural exponents */
695 template<class Type>
696 void posintpower(
697  const vector<Type>& in, /**< vector which first argument is base */
698  vector<Type>& out, /**< vector where to store result in first argument */
699  size_t exponent /**< exponent */
700  )
701 {
702  out[0] = pow(in[0], (int)exponent);
703 }
704 
705 #endif
706 
707 
708 #ifndef NO_CPPAD_USER_ATOMIC
709 
710 /** Automatic differentiation of x -> sign(x)abs(x)^p, p>=1, as CppAD user-atomic function.
711  *
712  * This class implements forward and reverse operations for the function x -> sign(x)abs(x)^p for use within CppAD.
713  * While we otherwise would have to use discontinuous sign and abs functions, our own implementation allows to provide
714  * a continuously differentiable function.
715  *
716  * @todo treat the exponent as a (variable) argument to the function, with the assumption that we never differentiate w.r.t. it (this should make the approach threadsafe again)
717  */
718 template<class Type>
719 class atomic_signpower : public CppAD::atomic_base<Type>
720 {
721 public:
723  : CppAD::atomic_base<Type>("signpower"),
724  exponent(0.0)
725  {
726  /* indicate that we want to use bool-based sparsity pattern */
727  this->option(CppAD::atomic_base<Type>::bool_sparsity_enum);
728  }
729 
730 private:
731  /** exponent for use in next call to forward or reverse */
732  SCIP_Real exponent;
733 
734  /** stores exponent corresponding to next call to forward or reverse
735  *
736  * How is this supposed to be threadsafe? (we use only one global instantiation of this class)
737  */
738  virtual void set_id(size_t id)
739  {
740  exponent = SCIPexprGetSignPowerExponent((SCIP_EXPR*)(void*)id);
741  }
742 
743  /** forward sweep of signpower
744  *
745  * Given the taylor coefficients for x, we have to compute the taylor coefficients for f(x),
746  * that is, given tx = (x, x', x'', ...), we compute the coefficients ty = (y, y', y'', ...)
747  * in the taylor expansion of f(x) = sign(x)abs(x)^p.
748  * Thus, y = sign(x)abs(x)^p
749  * = sign(tx[0])abs(tx[0])^p,
750  * y' = p * abs(x)^(p-1) * x'
751  * = p * abs(tx[0])^(p-1) * tx[1],
752  * y'' = 1/2 * p * (p-1) * sign(x) * abs(x)^(p-2) * x'^2 + p * abs(x)^(p-1) * x''
753  * = 1/2 * p * (p-1) * sign(tx[0]) * abs(tx[0])^(p-2) * tx[1]^2 + p * abs(tx[0])^(p-1) * tx[2]
754  */
755  bool forward(
756  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
757  size_t p, /**< highest order Taylor coefficient that we are evaluating */
758  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
759  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
760  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
761  CppAD::vector<Type>& ty /**< vector to store taylor coefficients of y */
762  )
763  {
764  assert(exponent > 0.0);
765  assert(tx.size() >= p+1);
766  assert(ty.size() >= p+1);
767  assert(q <= p);
768 
769  if( vx.size() > 0 )
770  {
771  assert(vx.size() == 1);
772  assert(vy.size() == 1);
773  assert(p == 0);
774 
775  vy[0] = vx[0];
776  }
777 
778  if( q == 0 /* q <= 0 && 0 <= p */ )
779  {
780  ty[0] = SIGN(tx[0]) * pow(REALABS(tx[0]), exponent);
781  }
782 
783  if( q <= 1 && 1 <= p )
784  {
785  ty[1] = pow(REALABS(tx[0]), exponent - 1.0) * tx[1];
786  ty[1] *= exponent;
787  }
788 
789  if( q <= 2 && 2 <= p )
790  {
791  if( exponent != 2.0 )
792  {
793  ty[2] = SIGN(tx[0]) * pow(REALABS(tx[0]), exponent - 2.0) * tx[1] * tx[1];
794  ty[2] *= (exponent - 1.0) / 2.0;
795  ty[2] += pow(REALABS(tx[0]), exponent - 1.0) * tx[2];
796  ty[2] *= exponent;
797  }
798  else
799  {
800  // y'' = 2 (1/2 * sign(x) * x'^2 + |x|*x'') = sign(tx[0]) * tx[1]^2 + 2 * abs(tx[0]) * tx[2]
801  ty[2] = SIGN(tx[0]) * tx[1] * tx[1];
802  ty[2] += 2.0 * REALABS(tx[0]) * tx[2];
803  }
804  }
805 
806  /* higher order derivatives not implemented */
807  if( p > 2 )
808  return false;
809 
810  return true;
811  }
812 
813  /** reverse sweep of signpower
814  *
815  * Assume y(x) is a function of the taylor coefficients of f(x) = sign(x)|x|^p for x, i.e.,
816  * y(x) = [ f(x), f'(x), f''(x), ... ].
817  * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
818  * where x^[l] is the l'th taylor coefficient (x, x', x'', ...) and h(x) = g(y(x)) for some function g:R^k -> R.
819  * That is, we have to compute
820  *\f$
821  * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
822  * = \sum_{i=0}^k (\partial g / \partial y_i) * (\partial y_i / \partial x^[l])
823  * = \sum_{i=0}^k py[i] * (\partial y_i / \partial x^[l])
824  *\f$
825  *
826  * For k = 0, this means
827  *\f$
828  * px[0] = py[0] * (\partial y_0 / \partial x^[0])
829  * = py[0] * (\partial f(x) / \partial x)
830  * = py[0] * p * abs(tx[0])^(p-1)
831  * \f$
832  *
833  * For k = 1, this means
834  *\f$
835  * px[0] = py[0] * (\partial y_0 / \partial x^[0]) + py[1] * (\partial y_1 / \partial x^[0])
836  * = py[0] * (\partial f(x) / \partial x) + py[1] * (\partial f'(x) / \partial x)
837  * = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1]
838  * px[1] = py[0] * (\partial y_0 / \partial x^[1]) + py[1] * (\partial y_1 / \partial x^[1])
839  * = py[0] * (\partial f(x) / \partial x') + py[1] * (\partial f'(x) / \partial x')
840  * = py[0] * 0 + py[1] * p * abs(tx[0])^(p-1)
841  * \f$
842  */
843  bool reverse(
844  size_t p, /**< highest order Taylor coefficient that we are evaluating */
845  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
846  const CppAD::vector<Type>& ty, /**< values for taylor coefficients of y */
847  CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
848  const CppAD::vector<Type>& py /**< values for partial derivatives of g(x) w.r.t. y */
849  )
850  { /*lint --e{715}*/
851  assert(exponent > 1);
852  assert(px.size() >= p+1);
853  assert(py.size() >= p+1);
854  assert(tx.size() >= p+1);
855 
856  switch( p )
857  {
858  case 0:
859  // px[0] = py[0] * p * pow(abs(tx[0]), p-1);
860  px[0] = py[0] * pow(REALABS(tx[0]), exponent - 1.0);
861  px[0] *= p;
862  break;
863 
864  case 1:
865  if( exponent != 2.0 )
866  {
867  // px[0] = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1]
868  px[0] = py[1] * tx[1] * pow(REALABS(tx[0]), exponent - 2.0) * SIGN(tx[0]);
869  px[0] *= exponent - 1.0;
870  px[0] += py[0] * pow(REALABS(tx[0]), exponent - 1.0);
871  px[0] *= exponent;
872  // px[1] = py[1] * p * abs(tx[0])^(p-1)
873  px[1] = py[1] * pow(REALABS(tx[0]), exponent - 1.0);
874  px[1] *= exponent;
875  }
876  else
877  {
878  // px[0] = py[0] * 2.0 * abs(tx[0]) + py[1] * 2.0 * sign(tx[0]) * tx[1]
879  px[0] = py[1] * tx[1] * SIGN(tx[0]);
880  px[0] += py[0] * REALABS(tx[0]);
881  px[0] *= 2.0;
882  // px[1] = py[1] * 2.0 * abs(tx[0])
883  px[1] = py[1] * REALABS(tx[0]);
884  px[1] *= 2.0;
885  }
886  break;
887 
888  default:
889  return false;
890  }
891 
892  return true;
893  }
894 
895  using CppAD::atomic_base<Type>::for_sparse_jac;
896 
897  /** computes sparsity of jacobian during a forward sweep
898  *
899  * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R.
900  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
901  */
902  bool for_sparse_jac(
903  size_t q, /**< number of columns in R */
904  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
905  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
906  )
907  {
908  return univariate_for_sparse_jac(q, r, s);
909  }
910 
911  using CppAD::atomic_base<Type>::rev_sparse_jac;
912 
913  /** computes sparsity of jacobian during a reverse sweep
914  *
915  * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x).
916  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
917  */
918  bool rev_sparse_jac(
919  size_t q, /**< number of rows in R */
920  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
921  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
922  )
923  {
924  return univariate_rev_sparse_jac(q, r, s);
925  }
926 
927  using CppAD::atomic_base<Type>::rev_sparse_hes;
928 
929  /** computes sparsity of hessian during a reverse sweep
930  *
931  * Assume V(x) = (g(f(x)))'' R with f(x) = sign(x)abs(x)^p for a function g:R->R and a matrix R.
932  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
933  */
934  bool rev_sparse_hes(
935  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
936  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
937  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
938  size_t q, /**< number of columns in S and R */
939  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
940  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
941  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
942  )
943  {
944  return univariate_rev_sparse_hes(vx, s, t, q, r, u, v);
945  }
946 
947 };
948 
949 /** Specialization of atomic_signpower template for intervals */
950 template<>
951 class atomic_signpower<SCIPInterval> : public CppAD::atomic_base<SCIPInterval>
952 {
953 public:
955  : CppAD::atomic_base<SCIPInterval>("signpowerint"),
956  exponent(0.0)
957  {
958  /* indicate that we want to use bool-based sparsity pattern */
959  this->option(CppAD::atomic_base<SCIPInterval>::bool_sparsity_enum);
960  }
961 
962 private:
963  /** exponent for use in next call to forward or reverse */
964  SCIP_Real exponent;
965 
966  /** stores exponent corresponding to next call to forward or reverse
967  *
968  * How is this supposed to be threadsafe? (we use only one global instantiation of this class)
969  */
970  virtual void set_id(size_t id)
971  {
972  exponent = SCIPexprGetSignPowerExponent((SCIP_EXPR*)(void*)id);
973  }
974 
975  /** specialization of atomic_signpower::forward template for SCIPinterval
976  *
977  * @todo try to compute tighter resultants
978  */
979  bool forward(
980  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
981  size_t p, /**< highest order Taylor coefficient that we are evaluating */
982  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
983  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
984  const CppAD::vector<SCIPInterval>& tx, /**< values for taylor coefficients of x */
985  CppAD::vector<SCIPInterval>& ty /**< vector to store taylor coefficients of y */
986  )
987  {
988  assert(exponent > 0.0);
989  assert(tx.size() >= p+1);
990  assert(ty.size() >= p+1);
991  assert(q <= p);
992 
993  if( vx.size() > 0 )
994  {
995  assert(vx.size() == 1);
996  assert(vy.size() == 1);
997  assert(p == 0);
998 
999  vy[0] = vx[0];
1000  }
1001 
1002  if( q == 0 /* q <= 0 && 0 <= p */ )
1003  {
1004  ty[0] = CppAD::signpow(tx[0], exponent);
1005  }
1006 
1007  if( q <= 1 && 1 <= p )
1008  {
1009  ty[1] = CppAD::pow(CppAD::abs(tx[0]), exponent - 1.0) * tx[1];
1010  ty[1] *= p;
1011  }
1012 
1013  if( q <= 2 && 2 <= p )
1014  {
1015  if( exponent != 2.0 )
1016  {
1017  ty[2] = CppAD::signpow(tx[0], exponent - 2.0) * CppAD::square(tx[1]);
1018  ty[2] *= (exponent - 1.0) / 2.0;
1019  ty[2] += CppAD::pow(CppAD::abs(tx[0]), exponent - 1.0) * tx[2];
1020  ty[2] *= exponent;
1021  }
1022  else
1023  {
1024  // y'' = 2 (1/2 * sign(x) * x'^2 + |x|*x'') = sign(tx[0]) * tx[1]^2 + 2 * abs(tx[0]) * tx[2]
1025  ty[2] = CppAD::sign(tx[0]) * CppAD::square(tx[1]);
1026  ty[2] += 2.0 * CppAD::abs(tx[0]) * tx[2];
1027  }
1028  }
1029 
1030  /* higher order derivatives not implemented */
1031  if( p > 2 )
1032  return false;
1033 
1034  return true;
1035  }
1036 
1037  /** specialization of atomic_signpower::reverse template for SCIPinterval
1038  *
1039  * @todo try to compute tighter resultants
1040  */
1041  bool reverse(
1042  size_t p, /**< highest order Taylor coefficient that we are evaluating */
1043  const CppAD::vector<SCIPInterval>& tx, /**< values for taylor coefficients of x */
1044  const CppAD::vector<SCIPInterval>& ty, /**< values for taylor coefficients of y */
1045  CppAD::vector<SCIPInterval>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
1046  const CppAD::vector<SCIPInterval>& py /**< values for partial derivatives of g(x) w.r.t. y */
1047  )
1048  { /*lint --e{715} */
1049  assert(exponent > 1);
1050  assert(px.size() >= p+1);
1051  assert(py.size() >= p+1);
1052  assert(tx.size() >= p+1);
1053 
1054  switch( p )
1055  {
1056  case 0:
1057  // px[0] = py[0] * p * pow(abs(tx[0]), p-1);
1058  px[0] = py[0] * CppAD::pow(CppAD::abs(tx[0]), exponent - 1.0);
1059  px[0] *= exponent;
1060  break;
1061 
1062  case 1:
1063  if( exponent != 2.0 )
1064  {
1065  // px[0] = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1]
1066  px[0] = py[1] * tx[1] * CppAD::signpow(tx[0], exponent - 2.0);
1067  px[0] *= exponent - 1.0;
1068  px[0] += py[0] * CppAD::pow(CppAD::abs(tx[0]), exponent - 1.0);
1069  px[0] *= exponent;
1070  // px[1] = py[1] * p * abs(tx[0])^(p-1)
1071  px[1] = py[1] * CppAD::pow(CppAD::abs(tx[0]), exponent - 1.0);
1072  px[1] *= exponent;
1073  }
1074  else
1075  {
1076  // px[0] = py[0] * 2.0 * abs(tx[0]) + py[1] * 2.0 * sign(tx[0]) * tx[1]
1077  px[0] = py[1] * tx[1] * CppAD::sign(tx[0]);
1078  px[0] += py[0] * CppAD::abs(tx[0]);
1079  px[0] *= 2.0;
1080  // px[1] = py[1] * 2.0 * abs(tx[0])
1081  px[1] = py[1] * CppAD::abs(tx[0]);
1082  px[1] *= 2.0;
1083  }
1084  break;
1085 
1086  default:
1087  return false;
1088  }
1089 
1090  return true;
1091  }
1092 
1093  using CppAD::atomic_base<SCIPInterval>::for_sparse_jac;
1094 
1095  /** computes sparsity of jacobian during a forward sweep
1096  *
1097  * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R.
1098  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
1099  */
1100  bool for_sparse_jac(
1101  size_t q, /**< number of columns in R */
1102  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
1103  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
1104  )
1105  {
1106  return univariate_for_sparse_jac(q, r, s);
1107  }
1108 
1109  using CppAD::atomic_base<SCIPInterval>::rev_sparse_jac;
1110 
1111  /** computes sparsity of jacobian during a reverse sweep
1112  *
1113  * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x).
1114  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
1115  */
1116  bool rev_sparse_jac(
1117  size_t q, /**< number of rows in R */
1118  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
1119  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
1120  )
1121  {
1122  return univariate_rev_sparse_jac(q, r, s);
1123  }
1124 
1125  using CppAD::atomic_base<SCIPInterval>::rev_sparse_hes;
1126 
1127  /** computes sparsity of hessian during a reverse sweep
1128  *
1129  * Assume V(x) = (g(f(x)))'' R with f(x) = sign(x)abs(x)^p for a function g:R->R and a matrix R.
1130  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
1131  */
1132  bool rev_sparse_hes(
1133  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
1134  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
1135  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
1136  size_t q, /**< number of columns in S and R */
1137  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
1138  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
1139  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
1140  )
1141  {
1142  return univariate_rev_sparse_hes(vx, s, t, q, r, u, v);
1143  }
1144 };
1145 
1146 /** template for evaluation for signpower operator */
1147 template<class Type>
1148 static
1149 void evalSignPower(
1150  Type& resultant, /**< resultant */
1151  const Type& arg, /**< operand */
1152  SCIP_EXPR* expr /**< expression that holds the exponent */
1153  )
1154 {
1155  vector<Type> in(1, arg);
1156  vector<Type> out(1);
1157 
1159  sp(in, out, (size_t)(void*)expr);
1160 
1161  resultant = out[0];
1162  return;
1163 }
1164 
1165 #else
1166 
1167 /** template for evaluation for signpower operator
1168  *
1169  * Only implemented for real numbers, thus gives error by default.
1170  */
1171 template<class Type>
1172 static
1173 void evalSignPower(
1174  Type& resultant, /**< resultant */
1175  const Type& arg, /**< operand */
1176  SCIP_EXPR* expr /**< expression that holds the exponent */
1177  )
1178 { /*lint --e{715}*/
1179  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
1180  "evalSignPower()",
1181  "Error: SignPower not implemented for this value type"
1182  );
1183 }
1184 
1185 /** specialization of signpower evaluation for real numbers */
1186 template<>
1187 void evalSignPower(
1188  CppAD::AD<double>& resultant, /**< resultant */
1189  const CppAD::AD<double>& arg, /**< operand */
1190  SCIP_EXPR* expr /**< expression that holds the exponent */
1191  )
1192 {
1193  SCIP_Real exponent;
1194 
1195  exponent = SCIPexprGetSignPowerExponent(expr);
1196 
1197  if( arg == 0.0 )
1198  resultant = 0.0;
1199  else if( arg > 0.0 )
1200  resultant = pow( arg, exponent);
1201  else
1202  resultant = -pow(-arg, exponent);
1203 }
1204 
1205 #endif
1206 
1207 /** template for evaluation for minimum operator
1208  *
1209  * Only implemented for real numbers, thus gives error by default.
1210  * @todo implement own userad function
1211  */
1212 template<class Type>
1213 static
1214 void evalMin(
1215  Type& resultant, /**< resultant */
1216  const Type& arg1, /**< first operand */
1217  const Type& arg2 /**< second operand */
1218  )
1219 { /*lint --e{715,1764}*/
1220  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
1221  "evalMin()",
1222  "Error: Min not implemented for this value type"
1223  );
1224 }
1225 
1226 /** specialization of minimum evaluation for real numbers */
1227 template<>
1228 void evalMin(
1229  CppAD::AD<double>& resultant, /**< resultant */
1230  const CppAD::AD<double>& arg1, /**< first operand */
1231  const CppAD::AD<double>& arg2 /**< second operand */
1232  )
1233 {
1234  resultant = MIN(arg1, arg2);
1235 }
1236 
1237 /** template for evaluation for maximum operator
1238  *
1239  * Only implemented for real numbers, thus gives error by default.
1240  * @todo implement own userad function
1241  */
1242 template<class Type>
1243 static
1244 void evalMax(
1245  Type& resultant, /**< resultant */
1246  const Type& arg1, /**< first operand */
1247  const Type& arg2 /**< second operand */
1248  )
1249 { /*lint --e{715,1764}*/
1250  CppAD::ErrorHandler::Call(true, __LINE__, __FILE__,
1251  "evalMax()",
1252  "Error: Max not implemented for this value type"
1253  );
1254 }
1255 
1256 /** specialization of maximum evaluation for real numbers */
1257 template<>
1258 void evalMax(
1259  CppAD::AD<double>& resultant, /**< resultant */
1260  const CppAD::AD<double>& arg1, /**< first operand */
1261  const CppAD::AD<double>& arg2 /**< second operand */
1262  )
1263 {
1264  resultant = MAX(arg1, arg2);
1265 }
1266 
1267 /** template for evaluation for square-root operator
1268  *
1269  * Default is to use the standard sqrt-function.
1270  */
1271 template<class Type>
1272 static
1273 void evalSqrt(
1274  Type& resultant, /**< resultant */
1275  const Type& arg /**< operand */
1276  )
1277 {
1278  resultant = sqrt(arg);
1279 }
1280 
1281 /** specialization of square-root operator for numbers
1282  *
1283  * We perturb the function a little bit so that it's derivatives are defined in 0.0.
1284  */
1285 template<>
1287  CppAD::AD<double>& resultant, /**< resultant */
1288  const CppAD::AD<double>& arg /**< operand */
1289  )
1290 {
1291  resultant = sqrt(arg + 1e-20) - 1e-10;
1292 }
1293 
1294 /** template for evaluation for absolute value operator */
1295 template<class Type>
1296 static
1297 void evalAbs(
1298  Type& resultant, /**< resultant */
1299  const Type& arg /**< operand */
1300  )
1301 {
1302  resultant = abs(arg);
1303 }
1304 
1305 /** specialization of absolute value evaluation for intervals
1306  *
1307  * Use sqrt(x^2) for now @todo implement own userad function.
1308  */
1309 template<>
1310 void evalAbs(
1311  CppAD::AD<SCIPInterval>& resultant, /**< resultant */
1312  const CppAD::AD<SCIPInterval>& arg /**< operand */
1313  )
1314 {
1315  vector<CppAD::AD<SCIPInterval> > in(1, arg);
1316  vector<CppAD::AD<SCIPInterval> > out(1);
1317 
1318  posintpower(in, out, 2);
1319 
1320  resultant = sqrt(out[0]);
1321 }
1322 
1323 /** integer power operation for arbitrary integer exponents */
1324 template<class Type>
1325 static
1326 void evalIntPower(
1327  Type& resultant, /**< resultant */
1328  const Type& arg, /**< operand */
1329  const int exponent /**< exponent */
1330  )
1331 {
1332  if( exponent > 1 )
1333  {
1334  vector<Type> in(1, arg);
1335  vector<Type> out(1);
1336 
1337  posintpower(in, out, exponent);
1338 
1339  resultant = out[0];
1340  return;
1341  }
1342 
1343  if( exponent < -1 )
1344  {
1345  vector<Type> in(1, arg);
1346  vector<Type> out(1);
1347 
1348  posintpower(in, out, -exponent);
1349 
1350  resultant = Type(1.0)/out[0];
1351  return;
1352  }
1353 
1354  if( exponent == 1 )
1355  {
1356  resultant = arg;
1357  return;
1358  }
1359 
1360  if( exponent == 0 )
1361  {
1362  resultant = Type(1.0);
1363  return;
1364  }
1365 
1366  assert(exponent == -1);
1367  resultant = Type(1.0)/arg;
1368 }
1369 
1370 /** CppAD compatible evaluation of an expression for given arguments and parameters */
1371 template<class Type>
1372 static
1373 SCIP_RETCODE eval(
1374  SCIP_EXPR* expr, /**< expression */
1375  const vector<Type>& x, /**< values of variables */
1376  SCIP_Real* param, /**< values of parameters */
1377  Type& val /**< buffer to store expression value */
1378  )
1379 {
1380  Type* buf = 0;
1381 
1382  assert(expr != NULL);
1383 
1384  /* todo use SCIP_MAXCHILD_ESTIMATE as in expression.c */
1385 
1386  if( SCIPexprGetNChildren(expr) )
1387  {
1388  if( BMSallocMemoryArray(&buf, SCIPexprGetNChildren(expr)) == NULL ) /*lint !e666*/
1389  return SCIP_NOMEMORY;
1390 
1391  for( int i = 0; i < SCIPexprGetNChildren(expr); ++i )
1392  {
1393  SCIP_CALL( eval(SCIPexprGetChildren(expr)[i], x, param, buf[i]) );
1394  }
1395  }
1396 
1397  switch(SCIPexprGetOperator(expr))
1398  {
1399  case SCIP_EXPR_VARIDX:
1400  assert(SCIPexprGetOpIndex(expr) < (int)x.size());
1401  val = x[SCIPexprGetOpIndex(expr)];
1402  break;
1403 
1404  case SCIP_EXPR_CONST:
1405  val = SCIPexprGetOpReal(expr);
1406  break;
1407 
1408  case SCIP_EXPR_PARAM:
1409  assert(param != NULL);
1410  val = param[SCIPexprGetOpIndex(expr)];
1411  break;
1412 
1413  case SCIP_EXPR_PLUS:
1414  assert( buf != 0 );
1415  val = buf[0] + buf[1];
1416  break;
1417 
1418  case SCIP_EXPR_MINUS:
1419  assert( buf != 0 );
1420  val = buf[0] - buf[1];
1421  break;
1422 
1423  case SCIP_EXPR_MUL:
1424  assert( buf != 0 );
1425  val = buf[0] * buf[1];
1426  break;
1427 
1428  case SCIP_EXPR_DIV:
1429  assert( buf != 0 );
1430  val = buf[0] / buf[1];
1431  break;
1432 
1433  case SCIP_EXPR_SQUARE:
1434  assert( buf != 0 );
1435  evalIntPower(val, buf[0], 2);
1436  break;
1437 
1438  case SCIP_EXPR_SQRT:
1439  assert( buf != 0 );
1440  evalSqrt(val, buf[0]);
1441  break;
1442 
1443  case SCIP_EXPR_REALPOWER:
1444  assert( buf != 0 );
1445  val = CppAD::pow(buf[0], SCIPexprGetRealPowerExponent(expr));
1446  break;
1447 
1448  case SCIP_EXPR_INTPOWER:
1449  assert( buf != 0 );
1450  evalIntPower(val, buf[0], SCIPexprGetIntPowerExponent(expr));
1451  break;
1452 
1453  case SCIP_EXPR_SIGNPOWER:
1454  assert( buf != 0 );
1455  evalSignPower(val, buf[0], expr);
1456  break;
1457 
1458  case SCIP_EXPR_EXP:
1459  assert( buf != 0 );
1460  val = exp(buf[0]);
1461  break;
1462 
1463  case SCIP_EXPR_LOG:
1464  assert( buf != 0 );
1465  val = log(buf[0]);
1466  break;
1467 
1468  case SCIP_EXPR_SIN:
1469  assert( buf != 0 );
1470  val = sin(buf[0]);
1471  break;
1472 
1473  case SCIP_EXPR_COS:
1474  assert( buf != 0 );
1475  val = cos(buf[0]);
1476  break;
1477 
1478  case SCIP_EXPR_TAN:
1479  assert( buf != 0 );
1480  val = tan(buf[0]);
1481  break;
1482 #ifdef SCIP_DISABLED_CODE /* these operators are currently disabled */
1483  case SCIP_EXPR_ERF:
1484  assert( buf != 0 );
1485  val = erf(buf[0]);
1486  break;
1487 
1488  case SCIP_EXPR_ERFI:
1489  return SCIP_ERROR;
1490 #endif
1491  case SCIP_EXPR_MIN:
1492  assert( buf != 0 );
1493  evalMin(val, buf[0], buf[1]);
1494  break;
1495 
1496  case SCIP_EXPR_MAX:
1497  assert( buf != 0 );
1498  evalMax(val, buf[0], buf[1]);
1499  break;
1500 
1501  case SCIP_EXPR_ABS:
1502  assert( buf != 0 );
1503  evalAbs(val, buf[0]);
1504  break;
1505 
1506  case SCIP_EXPR_SIGN:
1507  assert( buf != 0 );
1508  val = sign(buf[0]);
1509  break;
1510 
1511  case SCIP_EXPR_SUM:
1512  assert( buf != 0 );
1513  val = 0.0;
1514  for (int i = 0; i < SCIPexprGetNChildren(expr); ++i)
1515  val += buf[i];
1516  break;
1517 
1518  case SCIP_EXPR_PRODUCT:
1519  assert( buf != 0 );
1520  val = 1.0;
1521  for (int i = 0; i < SCIPexprGetNChildren(expr); ++i)
1522  val *= buf[i];
1523  break;
1524 
1525  case SCIP_EXPR_LINEAR:
1526  {
1527  SCIP_Real* coefs;
1528 
1529  coefs = SCIPexprGetLinearCoefs(expr);
1530  assert(coefs != NULL || SCIPexprGetNChildren(expr) == 0);
1531 
1532  assert( buf != 0 );
1533  val = SCIPexprGetLinearConstant(expr);
1534  for (int i = 0; i < SCIPexprGetNChildren(expr); ++i)
1535  val += coefs[i] * buf[i]; /*lint !e613*/
1536  break;
1537  }
1538 
1539  case SCIP_EXPR_QUADRATIC:
1540  {
1541  SCIP_Real* lincoefs;
1542  SCIP_QUADELEM* quadelems;
1543  int nquadelems;
1544  SCIP_Real sqrcoef;
1545  Type lincoef;
1546  vector<Type> in(1);
1547  vector<Type> out(1);
1548 
1549  assert( buf != 0 );
1550 
1551  lincoefs = SCIPexprGetQuadLinearCoefs(expr);
1552  nquadelems = SCIPexprGetNQuadElements(expr);
1553  quadelems = SCIPexprGetQuadElements(expr);
1554  assert(quadelems != NULL || nquadelems == 0);
1555 
1556  SCIPexprSortQuadElems(expr);
1557 
1558  val = SCIPexprGetQuadConstant(expr);
1559 
1560  /* for each argument, we collect it's linear index from lincoefs, it's square coefficients and all factors from bilinear terms
1561  * then we compute the interval sqrcoef*x^2 + lincoef*x and add it to result */
1562  int i = 0;
1563  for( int argidx = 0; argidx < SCIPexprGetNChildren(expr); ++argidx )
1564  {
1565  if( i == nquadelems || quadelems[i].idx1 > argidx ) /*lint !e613*/
1566  {
1567  /* there are no quadratic terms with argidx in its first argument, that should be easy to handle */
1568  if( lincoefs != NULL )
1569  val += lincoefs[argidx] * buf[argidx];
1570  continue;
1571  }
1572 
1573  sqrcoef = 0.0;
1574  lincoef = lincoefs != NULL ? lincoefs[argidx] : 0.0;
1575 
1576  assert(i < nquadelems && quadelems[i].idx1 == argidx); /*lint !e613*/
1577  do
1578  {
1579  if( quadelems[i].idx2 == argidx ) /*lint !e613*/
1580  sqrcoef += quadelems[i].coef; /*lint !e613*/
1581  else
1582  lincoef += quadelems[i].coef * buf[quadelems[i].idx2]; /*lint !e613*/
1583  ++i;
1584  } while( i < nquadelems && quadelems[i].idx1 == argidx ); /*lint !e613*/
1585  assert(i == nquadelems || quadelems[i].idx1 > argidx); /*lint !e613*/
1586 
1587  /* this is not as good as what we can get from SCIPintervalQuad, but easy to implement */
1588  if( sqrcoef != 0.0 )
1589  {
1590  in[0] = buf[argidx];
1591  posintpower(in, out, 2);
1592  val += sqrcoef * out[0];
1593  }
1594 
1595  val += lincoef * buf[argidx];
1596  }
1597  assert(i == nquadelems);
1598 
1599  break;
1600  }
1601 
1602  case SCIP_EXPR_POLYNOMIAL:
1603  {
1604  SCIP_EXPRDATA_MONOMIAL** monomials;
1605  Type childval;
1606  Type monomialval;
1607  SCIP_Real exponent;
1608  int nmonomials;
1609  int nfactors;
1610  int* childidxs;
1611  SCIP_Real* exponents;
1612  int i;
1613  int j;
1614 
1615  assert( buf != 0 );
1616 
1617  val = SCIPexprGetPolynomialConstant(expr);
1618 
1619  nmonomials = SCIPexprGetNMonomials(expr);
1620  monomials = SCIPexprGetMonomials(expr);
1621 
1622  for( i = 0; i < nmonomials; ++i )
1623  {
1624  nfactors = SCIPexprGetMonomialNFactors(monomials[i]);
1625  childidxs = SCIPexprGetMonomialChildIndices(monomials[i]);
1626  exponents = SCIPexprGetMonomialExponents(monomials[i]);
1627  monomialval = SCIPexprGetMonomialCoef(monomials[i]);
1628 
1629  for( j = 0; j < nfactors; ++j )
1630  {
1631  assert(childidxs[j] >= 0);
1632  assert(childidxs[j] < SCIPexprGetNChildren(expr));
1633 
1634  childval = buf[childidxs[j]];
1635  exponent = exponents[j];
1636 
1637  /* cover some special exponents separately to avoid calling expensive pow function */
1638  if( exponent == 0.0 )
1639  continue;
1640  if( exponent == 1.0 )
1641  {
1642  monomialval *= childval;
1643  continue;
1644  }
1645  if( (int)exponent == exponent )
1646  {
1647  Type tmp;
1648  evalIntPower(tmp, childval, (int)exponent);
1649  monomialval *= tmp;
1650  continue;
1651  }
1652  if( exponent == 0.5 )
1653  {
1654  Type tmp;
1655  evalSqrt(tmp, childval);
1656  monomialval *= tmp;
1657  continue;
1658  }
1659  monomialval *= pow(childval, exponent);
1660  }
1661 
1662  val += monomialval;
1663  }
1664 
1665  break;
1666  }
1667 
1668  default:
1669  return SCIP_ERROR;
1670  }
1671 
1672  BMSfreeMemoryArrayNull(&buf);
1673 
1674  return SCIP_OKAY;
1675 }
1676 
1677 /** analysis an expression tree whether it requires retaping on every evaluation
1678  *
1679  * This may be the case if the evaluation sequence depends on values of operands (e.g., in case of abs, sign, signpower, ...).
1680  */
1681 static
1682 bool needAlwaysRetape(SCIP_EXPR* expr)
1683 {
1684  assert(expr != NULL);
1685  assert(SCIPexprGetChildren(expr) != NULL || SCIPexprGetNChildren(expr) == 0);
1686 
1687  for( int i = 0; i < SCIPexprGetNChildren(expr); ++i )
1688  {
1689  if( needAlwaysRetape(SCIPexprGetChildren(expr)[i]) )
1690  return true;
1691  }
1692 
1693  switch( SCIPexprGetOperator(expr) )
1694  {
1695  case SCIP_EXPR_MIN:
1696  case SCIP_EXPR_MAX:
1697  case SCIP_EXPR_ABS:
1698 #ifdef NO_CPPAD_USER_ATOMIC
1699  case SCIP_EXPR_SIGNPOWER:
1700 #endif
1701  return true;
1702 
1703  default: ;
1704  } /*lint !e788*/
1705 
1706  return false;
1707 }
1708 
1709 /** replacement for CppAD's default error handler
1710  *
1711  * In debug mode, CppAD gives an error when an evaluation contains a nan.
1712  * We do not want to stop execution in such a case, since the calling routine should check for nan's and decide what to do.
1713  * Since we cannot ignore this particular error, we ignore all.
1714  * @todo find a way to check whether the error corresponds to a nan and communicate this back
1715  */
1716 static
1717 void cppaderrorcallback(
1718  bool known, /**< is the error from a known source? */
1719  int line, /**< line where error occured */
1720  const char* file, /**< file where error occured */
1721  const char* cond, /**< error condition */
1722  const char* msg /**< error message */
1723  )
1724 {
1725  SCIPdebugMessage("ignore CppAD error from %sknown source %s:%d: msg: %s exp: %s\n", known ? "" : "un", file, line, msg, cond);
1726 }
1727 
1728 /* install our error handler */
1729 static CppAD::ErrorHandler errorhandler(cppaderrorcallback);
1730 
1731 /** gets name and version of expression interpreter */
1732 const char* SCIPexprintGetName(void)
1733 {
1734  return CPPAD_PACKAGE_STRING;
1735 }
1736 
1737 /** gets descriptive text of expression interpreter */
1738 const char* SCIPexprintGetDesc(void)
1739 {
1740  return "Algorithmic Differentiation of C++ algorithms developed by B. Bell (www.coin-or.org/CppAD)";
1741 }
1742 
1743 /** gets capabilities of expression interpreter (using bitflags) */
1745  void
1746  )
1747 {
1751 }
1752 
1753 /** creates an expression interpreter object */
1755  BMS_BLKMEM* blkmem, /**< block memory data structure */
1756  SCIP_EXPRINT** exprint /**< buffer to store pointer to expression interpreter */
1757  )
1758 {
1759  assert(blkmem != NULL);
1760  assert(exprint != NULL);
1761 
1762  if( BMSallocMemory(exprint) == NULL )
1763  return SCIP_NOMEMORY;
1764 
1765  (*exprint)->blkmem = blkmem;
1766 
1767  return SCIP_OKAY;
1768 }
1769 
1770 /** frees an expression interpreter object */
1772  SCIP_EXPRINT** exprint /**< expression interpreter that should be freed */
1773  )
1774 {
1775  assert( exprint != NULL);
1776  assert(*exprint != NULL);
1777 
1778  BMSfreeMemory(exprint);
1779 
1780  return SCIP_OKAY;
1781 }
1782 
1783 /** compiles an expression tree and stores compiled data in expression tree */
1785  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1786  SCIP_EXPRTREE* tree /**< expression tree */
1787  )
1788 { /*lint --e{429} */
1789  assert(tree != NULL);
1790 
1792  if (!data)
1793  {
1794  data = new SCIP_EXPRINTDATA();
1795  assert( data != NULL );
1796  SCIPexprtreeSetInterpreterData(tree, data);
1797  SCIPdebugMessage("set interpreter data in tree %p to %p\n", (void*)tree, (void*)data);
1798  }
1799  else
1800  {
1801  data->need_retape = true;
1802  data->int_need_retape = true;
1803  }
1804 
1805  int n = SCIPexprtreeGetNVars(tree);
1806 
1807  data->X.resize(n);
1808  data->x.resize(n);
1809  data->Y.resize(1);
1810 
1811  data->int_X.resize(n);
1812  data->int_x.resize(n);
1813  data->int_Y.resize(1);
1814 
1815  if( data->root != NULL )
1816  {
1817  SCIPexprFreeDeep(exprint->blkmem, &data->root);
1818  }
1819 
1820  SCIP_EXPR* root = SCIPexprtreeGetRoot(tree);
1821 
1822  SCIP_CALL( SCIPexprCopyDeep(exprint->blkmem, &data->root, root) );
1823 
1824  data->need_retape_always = needAlwaysRetape(SCIPexprtreeGetRoot(tree));
1825 
1826  data->blkmem = exprint->blkmem;
1827 
1828  return SCIP_OKAY;
1829 }
1830 
1831 /** frees interpreter data */
1833  SCIP_EXPRINTDATA** interpreterdata /**< interpreter data that should freed */
1834  )
1835 {
1836  assert( interpreterdata != NULL);
1837  assert(*interpreterdata != NULL);
1838 
1839  if( (*interpreterdata)->root != NULL )
1840  SCIPexprFreeDeep((*interpreterdata)->blkmem, &(*interpreterdata)->root);
1841 
1842  delete *interpreterdata;
1843  *interpreterdata = NULL;
1844 
1845  return SCIP_OKAY;
1846 }
1847 
1848 /** notify expression interpreter that a new parameterization is used
1849  *
1850  * This probably causes retaping by AD algorithms.
1851  */
1853  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1854  SCIP_EXPRTREE* tree /**< expression tree */
1855  )
1856 {
1857  assert(exprint != NULL);
1858  assert(tree != NULL);
1859 
1861  if( data != NULL )
1862  {
1863  data->need_retape = true;
1864  data->int_need_retape = true;
1865  }
1866 
1867  return SCIP_OKAY;
1868 }
1869 
1870 /** evaluates an expression tree */
1872  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1873  SCIP_EXPRTREE* tree, /**< expression tree */
1874  SCIP_Real* varvals, /**< values of variables */
1875  SCIP_Real* val /**< buffer to store value */
1876  )
1877 {
1878  SCIP_EXPRINTDATA* data;
1879 
1880  assert(exprint != NULL);
1881  assert(tree != NULL);
1882  assert(varvals != NULL);
1883  assert(val != NULL);
1884 
1885  data = SCIPexprtreeGetInterpreterData(tree);
1886  assert(data != NULL);
1887  assert(SCIPexprtreeGetNVars(tree) == (int)data->X.size());
1888  assert(SCIPexprtreeGetRoot(tree) != NULL);
1889 
1890  int n = SCIPexprtreeGetNVars(tree);
1891 
1892  if( data->need_retape_always || data->need_retape )
1893  {
1894  for( int i = 0; i < n; ++i )
1895  {
1896  data->X[i] = varvals[i];
1897  data->x[i] = varvals[i];
1898  }
1899 
1900  CppAD::Independent(data->X);
1901 
1902  if( data->root != NULL )
1903  SCIP_CALL( eval(data->root, data->X, SCIPexprtreeGetParamVals(tree), data->Y[0]) );
1904  else
1905  data->Y[0] = 0.0;
1906 
1907  data->f.Dependent(data->X, data->Y);
1908 
1909  data->val = Value(data->Y[0]);
1910  SCIPdebugMessage("Eval retaped and computed value %g\n", data->val);
1911 
1912  // the following is required if the gradient shall be computed by a reverse sweep later
1913  // data->val = data->f.Forward(0, data->x)[0];
1914 
1915  data->need_retape = false;
1916  }
1917  else
1918  {
1919  assert((int)data->x.size() >= n);
1920  for( int i = 0; i < n; ++i )
1921  data->x[i] = varvals[i];
1922 
1923  data->val = data->f.Forward(0, data->x)[0]; /*lint !e1793*/
1924  SCIPdebugMessage("Eval used forward sweep to compute value %g\n", data->val);
1925  }
1926 
1927  *val = data->val;
1928 
1929  return SCIP_OKAY;
1930 }
1931 
1932 /** evaluates an expression tree on intervals */
1934  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1935  SCIP_EXPRTREE* tree, /**< expression tree */
1936  SCIP_Real infinity, /**< value for infinity */
1937  SCIP_INTERVAL* varvals, /**< interval values of variables */
1938  SCIP_INTERVAL* val /**< buffer to store interval value of expression */
1939  )
1940 {
1941  SCIP_EXPRINTDATA* data;
1942 
1943  assert(exprint != NULL);
1944  assert(tree != NULL);
1945  assert(varvals != NULL);
1946  assert(val != NULL);
1947 
1948  data = SCIPexprtreeGetInterpreterData(tree);
1949  assert(data != NULL);
1950  assert(SCIPexprtreeGetNVars(tree) == (int)data->int_X.size());
1951  assert(SCIPexprtreeGetRoot(tree) != NULL);
1952 
1953  int n = SCIPexprtreeGetNVars(tree);
1954 
1955  SCIPInterval::infinity = infinity;
1956 
1957  if( data->int_need_retape || data->need_retape_always )
1958  {
1959  for( int i = 0; i < n; ++i )
1960  {
1961  data->int_X[i] = varvals[i];
1962  data->int_x[i] = varvals[i];
1963  }
1964 
1965  CppAD::Independent(data->int_X);
1966 
1967  if( data->root != NULL )
1968  SCIP_CALL( eval(data->root, data->int_X, SCIPexprtreeGetParamVals(tree), data->int_Y[0]) );
1969  else
1970  data->int_Y[0] = 0.0;
1971 
1972  data->int_f.Dependent(data->int_X, data->int_Y);
1973 
1974  data->int_val = Value(data->int_Y[0]);
1975 
1976  data->int_need_retape = false;
1977  }
1978  else
1979  {
1980  assert((int)data->int_x.size() >= n);
1981  for( int i = 0; i < n; ++i )
1982  data->int_x[i] = varvals[i];
1983 
1984  data->int_val = data->int_f.Forward(0, data->int_x)[0]; /*lint !e1793*/
1985  }
1986 
1987  *val = data->int_val;
1988 
1989  return SCIP_OKAY;
1990 }
1991 
1992 /** computes value and gradient of an expression tree */
1994  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1995  SCIP_EXPRTREE* tree, /**< expression tree */
1996  SCIP_Real* varvals, /**< values of variables, can be NULL if new_varvals is FALSE */
1997  SCIP_Bool new_varvals, /**< have variable values changed since last call to a point evaluation routine? */
1998  SCIP_Real* val, /**< buffer to store expression value */
1999  SCIP_Real* gradient /**< buffer to store expression gradient, need to have length at least SCIPexprtreeGetNVars(tree) */
2000  )
2001 {
2002  assert(exprint != NULL);
2003  assert(tree != NULL);
2004  assert(varvals != NULL || new_varvals == FALSE);
2005  assert(val != NULL);
2006  assert(gradient != NULL);
2007 
2009  assert(data != NULL);
2010 
2011  if( new_varvals )
2012  {
2013  SCIP_CALL( SCIPexprintEval(exprint, tree, varvals, val) );
2014  }
2015  else
2016  *val = data->val;
2017 
2018  int n = SCIPexprtreeGetNVars(tree);
2019 
2020  vector<double> jac(data->f.Jacobian(data->x));
2021 
2022  for( int i = 0; i < n; ++i )
2023  gradient[i] = jac[i];
2024 
2025 /* disable debug output since we have no message handler here
2026 #ifdef SCIP_DEBUG
2027  SCIPdebugMessage("Grad for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2028  SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n");
2029  SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t %g", gradient[i]); printf("\n");
2030 #endif
2031 */
2032 
2033  return SCIP_OKAY;
2034 }
2035 
2036 /** computes interval value and interval gradient of an expression tree */
2038  SCIP_EXPRINT* exprint, /**< interpreter data structure */
2039  SCIP_EXPRTREE* tree, /**< expression tree */
2040  SCIP_Real infinity, /**< value for infinity */
2041  SCIP_INTERVAL* varvals, /**< interval values of variables, can be NULL if new_varvals is FALSE */
2042  SCIP_Bool new_varvals, /**< have variable interval values changed since last call to an interval evaluation routine? */
2043  SCIP_INTERVAL* val, /**< buffer to store expression interval value */
2044  SCIP_INTERVAL* gradient /**< buffer to store expression interval gradient, need to have length at least SCIPexprtreeGetNVars(tree) */
2045  )
2046 {
2047  assert(exprint != NULL);
2048  assert(tree != NULL);
2049  assert(varvals != NULL || new_varvals == false);
2050  assert(val != NULL);
2051  assert(gradient != NULL);
2052 
2054  assert(data != NULL);
2055 
2056  if (new_varvals)
2057  SCIP_CALL( SCIPexprintEvalInt(exprint, tree, infinity, varvals, val) );
2058  else
2059  *val = data->int_val;
2060 
2061  int n = SCIPexprtreeGetNVars(tree);
2062 
2063  vector<SCIPInterval> jac(data->int_f.Jacobian(data->int_x));
2064 
2065  for (int i = 0; i < n; ++i)
2066  gradient[i] = jac[i];
2067 
2068 /* disable debug output since we have no message handler here
2069 #ifdef SCIP_DEBUG
2070  SCIPdebugMessage("GradInt for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2071  SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t [%g,%g]", SCIPintervalGetInf(data->int_x[i]), SCIPintervalGetSup(data->int_x[i])); printf("\n");
2072  SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t [%g,%g]", SCIPintervalGetInf(gradient[i]), SCIPintervalGetSup(gradient[i])); printf("\n");
2073 #endif
2074 */
2075 
2076  return SCIP_OKAY;
2077 }
2078 
2079 /** gives sparsity pattern of hessian
2080  *
2081  * NOTE: this function might be replaced later by something nicer.
2082  * Since the AD code might need to do a forward sweep, you should pass variable values in here.
2083  */
2085  SCIP_EXPRINT* exprint, /**< interpreter data structure */
2086  SCIP_EXPRTREE* tree, /**< expression tree */
2087  SCIP_Real* varvals, /**< values of variables */
2088  SCIP_Bool* sparsity /**< buffer to store sparsity pattern of Hessian, sparsity[i+n*j] indicates whether entry (i,j) is nonzero in the hessian */
2089  )
2090 {
2091  assert(exprint != NULL);
2092  assert(tree != NULL);
2093  assert(varvals != NULL);
2094  assert(sparsity != NULL);
2095 
2097  assert(data != NULL);
2098 
2099  int n = SCIPexprtreeGetNVars(tree);
2100  int nn = n*n;
2101 
2102  if( data->need_retape_always )
2103  {
2104  // @todo can we do something better here, e.g., by looking at the expression tree by ourself?
2105 
2106  for( int i = 0; i < nn; ++i )
2107  sparsity[i] = TRUE;
2108 
2109 /* disable debug output since we have no message handler here
2110 #ifdef SCIP_DEBUG
2111  SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2112  SCIPdebugMessage("sparsity = all elements, due to discontinuouities\n");
2113 #endif
2114 */
2115 
2116  return SCIP_OKAY;
2117  }
2118 
2119  if( data->need_retape )
2120  {
2121  SCIP_Real val;
2122  SCIP_CALL( SCIPexprintEval(exprint, tree, varvals, &val) );
2123  }
2124 
2125  SCIPdebugMessage("calling ForSparseJac\n");
2126 
2127  vector<bool> r(nn, false);
2128  for (int i = 0; i < n; ++i)
2129  r[i*n+i] = true; /*lint !e647 !e1793*/
2130  (void) data->f.ForSparseJac(n, r); // need to compute sparsity for Jacobian first
2131 
2132  SCIPdebugMessage("calling RevSparseHes\n");
2133 
2134  vector<bool> s(1, true);
2135  vector<bool> sparsehes(data->f.RevSparseHes(n, s));
2136 
2137  for( int i = 0; i < nn; ++i )
2138  sparsity[i] = sparsehes[i];
2139 
2140 /* disable debug output since we have no message handler here
2141 #ifdef SCIP_DEBUG
2142  SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2143  SCIPdebugMessage("sparsity ="); for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) if (sparsity[i*n+j]) printf(" (%d,%d)", i, j); printf("\n");
2144 #endif
2145 */
2146 
2147  return SCIP_OKAY;
2148 }
2149 
2150 /** computes value and dense hessian of an expression tree
2151  *
2152  * The full hessian is computed (lower left and upper right triangle).
2153  */
2155  SCIP_EXPRINT* exprint, /**< interpreter data structure */
2156  SCIP_EXPRTREE* tree, /**< expression tree */
2157  SCIP_Real* varvals, /**< values of variables, can be NULL if new_varvals is FALSE */
2158  SCIP_Bool new_varvals, /**< have variable values changed since last call to an evaluation routine? */
2159  SCIP_Real* val, /**< buffer to store function value */
2160  SCIP_Real* hessian /**< buffer to store hessian values, need to have size at least n*n */
2161  )
2162 {
2163  assert(exprint != NULL);
2164  assert(tree != NULL);
2165  assert(varvals != NULL || new_varvals == FALSE);
2166  assert(val != NULL);
2167  assert(hessian != NULL);
2168 
2170  assert(data != NULL);
2171 
2172  if( new_varvals )
2173  {
2174  SCIP_CALL( SCIPexprintEval(exprint, tree, varvals, val) );
2175  }
2176  else
2177  *val = data->val;
2178 
2179  int n = SCIPexprtreeGetNVars(tree);
2180 
2181 #if 1
2182  /* this one uses reverse mode */
2183  vector<double> hess(data->f.Hessian(data->x, 0));
2184 
2185  int nn = n*n;
2186  for (int i = 0; i < nn; ++i)
2187  hessian[i] = hess[i];
2188 
2189 #else
2190  /* this one uses forward mode */
2191  for( int i = 0; i < n; ++i )
2192  for( int j = 0; j < n; ++j )
2193  {
2194  vector<int> ii(1,i);
2195  vector<int> jj(1,j);
2196  hessian[i*n+j] = data->f.ForTwo(data->x, ii, jj)[0];
2197  }
2198 #endif
2199 
2200 /* disable debug output since we have no message handler here
2201 #ifdef SCIP_DEBUG
2202  SCIPdebugMessage("HessianDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2203  SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n");
2204  SCIPdebugMessage("hess ="); for (int i = 0; i < n*n; ++i) printf("\t %g", hessian[i]); printf("\n");
2205 #endif
2206 */
2207  return SCIP_OKAY;
2208 }
SCIP_RETCODE SCIPexprintEval(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Real *val)
SCIPInterval square(const SCIPInterval &x)
SCIP_RETCODE SCIPexprintHessianDense(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, SCIP_Real *hessian)
bool LessThanZero(const SCIPInterval &x)
const char * SCIPexprintGetName(void)
bool IdenticalEqualPar(const SCIPInterval &x, const SCIPInterval &y)
methods to interpret (evaluate) an expression tree "fast"
int SCIPexprGetNChildren(SCIP_EXPR *expr)
void SCIPexprtreeSetInterpreterData(SCIP_EXPRTREE *tree, SCIP_EXPRINTDATA *interpreterdata)
struct SCIP_QuadElement SCIP_QUADELEM
Definition: type_expr.h:105
SCIP_QUADELEM * SCIPexprGetQuadElements(SCIP_EXPR *expr)
const char * SCIPexprintGetDesc(void)
SCIPInterval pow(const SCIPInterval &x, const SCIPInterval &y)
SCIP_Real SCIPexprGetMonomialCoef(SCIP_EXPRDATA_MONOMIAL *monomial)
#define NULL
Definition: lpi_spx.cpp:129
struct SCIP_ExprData_Monomial SCIP_EXPRDATA_MONOMIAL
Definition: type_expr.h:108
#define SCIP_EXPRINTCAPABILITY_INTGRADIENT
SCIPInterval cos(const SCIPInterval &x)
SCIP_RETCODE SCIPexprintFree(SCIP_EXPRINT **exprint)
SCIP_Real * SCIPexprGetMonomialExponents(SCIP_EXPRDATA_MONOMIAL *monomial)
int SCIPexprtreeGetNVars(SCIP_EXPRTREE *tree)
#define FALSE
Definition: def.h:52
SCIP_RETCODE SCIPexprintGrad(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, SCIP_Real *gradient)
#define TRUE
Definition: def.h:51
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
SCIP_Real SCIPexprGetSignPowerExponent(SCIP_EXPR *expr)
static SCIP_Real infinity
SCIPInterval exp(const SCIPInterval &x)
SCIP_RETCODE SCIPexprCopyDeep(BMS_BLKMEM *blkmem, SCIP_EXPR **targetexpr, SCIP_EXPR *sourceexpr)
#define CPPAD_MAX_NUM_THREADS
#define SCIPdebugMessage
Definition: pub_message.h:77
SCIP_EXPRINTDATA * SCIPexprtreeGetInterpreterData(SCIP_EXPRTREE *tree)
int SCIPexprGetNMonomials(SCIP_EXPR *expr)
SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void)
unsigned int SCIP_EXPRINTCAPABILITY
int SCIPexprGetNQuadElements(SCIP_EXPR *expr)
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
SCIPInterval abs(const SCIPInterval &x)
public methods for expressions, expression trees, expression graphs, and related stuff ...
#define SCIP_EXPRINTCAPABILITY_INTFUNCVALUE
void SCIPexprSortQuadElems(SCIP_EXPR *expr)
SCIPInterval CondExpOp(enum CppAD::CompareOp cop, const SCIPInterval &left, const SCIPInterval &right, const SCIPInterval &trueCase, const SCIPInterval &falseCase)
SCIPInterval signpow(const SCIPInterval &x, const SCIP_Real p)
bool IdenticalOne(const SCIPInterval &x)
int * SCIPexprGetMonomialChildIndices(SCIP_EXPRDATA_MONOMIAL *monomial)
SCIP_Real SCIPexprGetQuadConstant(SCIP_EXPR *expr)
SCIP_RETCODE SCIPexprintFreeData(SCIP_EXPRINTDATA **interpreterdata)
void evalAbs(CppAD::AD< SCIPInterval > &resultant, const CppAD::AD< SCIPInterval > &arg)
SCIP_EXPR * SCIPexprtreeGetRoot(SCIP_EXPRTREE *tree)
SCIPInterval sqrt(const SCIPInterval &x)
void evalMin(CppAD::AD< double > &resultant, const CppAD::AD< double > &arg1, const CppAD::AD< double > &arg2)
SCIPInterval sign(const SCIPInterval &x)
#define SIGN(x)
#define REALABS(x)
Definition: def.h:146
SCIP_RETCODE SCIPexprintEvalInt(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_INTERVAL *val)
#define SCIP_CALL(x)
Definition: def.h:258
SCIP_EXPRDATA_MONOMIAL ** SCIPexprGetMonomials(SCIP_EXPR *expr)
SCIP_RETCODE SCIPexprintHessianSparsityDense(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool *sparsity)
SCIP_Real SCIPexprGetOpReal(SCIP_EXPR *expr)
struct SCIP_ExprTree SCIP_EXPRTREE
Definition: type_expr.h:91
SCIP_Real SCIPexprGetRealPowerExponent(SCIP_EXPR *expr)
int SCIPexprGetMonomialNFactors(SCIP_EXPRDATA_MONOMIAL *monomial)
#define SCIP_EXPRINTCAPABILITY_HESSIAN
bool IdenticalZero(const SCIPInterval &x)
#define SCIP_Bool
Definition: def.h:49
void SCIPexprFreeDeep(BMS_BLKMEM *blkmem, SCIP_EXPR **expr)
SCIPInterval sin(const SCIPInterval &x)
int SCIPexprGetOpIndex(SCIP_EXPR *expr)
#define MAX(x, y)
Definition: tclique_def.h:75
SCIP_Real SCIPexprGetPolynomialConstant(SCIP_EXPR *expr)
struct SCIP_ExprInt SCIP_EXPRINT
SCIP_Real * SCIPexprtreeGetParamVals(SCIP_EXPRTREE *tree)
SCIP_RETCODE SCIPexprintNewParametrization(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree)
SCIPInterval log(const SCIPInterval &x)
SCIP_RETCODE SCIPexprintGradInt(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_Bool new_varvals, SCIP_INTERVAL *val, SCIP_INTERVAL *gradient)
bool GreaterThanZero(const SCIPInterval &x)
bool EqualOpSeq(const SCIPInterval &x, const SCIPInterval &y)
void evalMax(CppAD::AD< double > &resultant, const CppAD::AD< double > &arg1, const CppAD::AD< double > &arg2)
bool IdenticalPar(const SCIPInterval &x)
#define SCIP_DEFAULT_INFINITY
Definition: def.h:128
#define SCIP_EXPRINTCAPABILITY_FUNCVALUE
SCIP_RETCODE SCIPexprintCreate(BMS_BLKMEM *blkmem, SCIP_EXPRINT **exprint)
std::ostream & operator<<(std::ostream &out, const SCIP_INTERVAL &x)
int Integer(const SCIPInterval &x)
SCIP_EXPROP SCIPexprGetOperator(SCIP_EXPR *expr)
#define SCIP_Real
Definition: def.h:123
int SCIPexprGetIntPowerExponent(SCIP_EXPR *expr)
SCIP_Real SCIPexprGetLinearConstant(SCIP_EXPR *expr)
bool LessThanOrZero(const SCIPInterval &x)
SCIP_Real * SCIPexprGetLinearCoefs(SCIP_EXPR *expr)
common defines and data types used in all packages of SCIP
void evalSqrt(CppAD::AD< double > &resultant, const CppAD::AD< double > &arg)
struct SCIP_ExprIntData SCIP_EXPRINTDATA
bool GreaterThanOrZero(const SCIPInterval &x)
#define SCIP_EXPRINTCAPABILITY_GRADIENT
SCIP_Real * SCIPexprGetQuadLinearCoefs(SCIP_EXPR *expr)
C++ extensions to interval arithmetics for provable bounds.
SCIP_RETCODE SCIPexprintCompile(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree)
struct SCIP_Expr SCIP_EXPR
Definition: type_expr.h:90