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-2024 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file exprinterpret_cppad.cpp
26  * @brief methods to interpret (evaluate) an expression "fast" using CppAD
27  * @ingroup DEFPLUGINS_EXPRINT
28  * @author Stefan Vigerske
29  */
30 
31 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
32 
33 #include "scip/exprinterpret.h"
34 #include "scip/def.h"
35 #include "scip/intervalarith.h"
36 #include "scip/pub_expr.h"
37 #include "scip/scip_expr.h"
38 #include "scip/expr_pow.h"
39 #include "scip/expr_exp.h"
40 #include "scip/expr_log.h"
41 #include "scip/expr_varidx.h"
42 
43 #include <cmath>
44 #include <cstring>
45 #include <algorithm>
46 #include <vector>
47 using std::vector;
48 
49 /* Turn off lint warning "747: Significant prototype coercion" and "732: Loss of sign".
50  * The first warning is generated for expressions like t[0], where t is a vector, since 0 is an integer constant, but a
51  * size_t is expected (usually long unsigned). The second is generated for expressions like t[n], where n is an
52  * integer. Both code pieces are likely to be correct. It seems to be impossible to inhibit these messages for
53  * vector<*>::operator[] only. */
54 /*lint --e{747,732}*/
55 
56 /* defining EVAL_USE_EXPRHDLR_ALWAYS uses the exprhdlr methods for the evaluation and differentiation of all expression (except for var and value)
57  * instead of only those that are not recognized in this code
58  * this is incomplete (no hessians) and slow (only forward-mode gradients), but can be useful for testing
59  */
60 // #define EVAL_USE_EXPRHDLR_ALWAYS
61 
62 /* defining NO_CPPAD_USER_ATOMIC disables the use of our own implementation of derivatives of power operators
63  * via CppAD's user-atomic function feature
64  * our customized implementation should give better results (tighter intervals) for the interval data type
65  * but we don't use interval types here anymore and the atomic operator does not look threadsafe (might be better in newer CppAD version)
66  */
67 #define NO_CPPAD_USER_ATOMIC
68 
69 /* fallback to non-thread-safe version if C++ is too old to have std::atomic */
70 #if __cplusplus < 201103L && defined(SCIP_THREADSAFE)
71 #undef SCIP_THREADSAFE
72 #endif
73 
74 /* CppAD needs to know a fixed upper bound on the number of threads at compile time.
75  * 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.
76  */
77 #ifndef CPPAD_MAX_NUM_THREADS
78 #ifdef SCIP_THREADSAFE
79 #define CPPAD_MAX_NUM_THREADS 64
80 #else
81 #define CPPAD_MAX_NUM_THREADS 1
82 #endif
83 #endif
84 
85 /* disable -Wshadow warnings for upcoming includes of CppAD if using some old GCC
86  * -Wshadow was too strict with some versions of GCC 4 (https://stackoverflow.com/questions/2958457/gcc-wshadow-is-too-strict)
87  * disable -Wimplicit-fallthrough as I don't want to maintain extra comments in CppAD code to suppress these
88  */
89 #ifdef __GNUC__
90 #if __GNUC__ == 4
91 #pragma GCC diagnostic ignored "-Wshadow"
92 #endif
93 #if __GNUC__ >= 7
94 #pragma GCC diagnostic ignored "-Wimplicit-fallthrough"
95 #endif
96 #endif
97 
98 // MS compiler doesn't have log2() - define an expensive workaround
99 #ifdef _MSC_VER
100 #define log2(x) (log((double)x) / log(2.0))
101 #endif
102 
103 #include <cppad/cppad.hpp>
104 #include <cppad/utility/error_handler.hpp>
105 
106 /* CppAD is not thread-safe by itself, but uses some static datastructures
107  * To run it in a multithreading environment, a special CppAD memory allocator that is aware of the multiple threads has to be used.
108  * This allocator requires to know the number of threads and a thread number for each thread.
109  * To implement this, we follow the team_pthread example of CppAD, which uses pthread's thread-specific data management.
110  */
111 #ifdef SCIP_THREADSAFE
112 
113 #include <atomic>
114 
115 /** currently registered number of threads */
116 static std::atomic_size_t ncurthreads{0};
117 static thread_local int thread_number{-1};
118 
119 /** CppAD callback function that indicates whether we are running in parallel mode */
120 static
121 bool in_parallel(void)
122 {
123  return ncurthreads > 1;
124 }
125 
126 /** CppAD callback function that returns the number of the current thread
127  *
128  * assigns a new number to the thread if new
129  */
130 static
131 size_t thread_num(void)
132 {
133  size_t threadnum;
134 
135  /* if no thread_number for this thread yet, then assign a new thread number to the current thread
136  */
137  if( thread_number == -1 )
138  {
139  thread_number = static_cast<int>(ncurthreads.fetch_add(1, std::memory_order_relaxed));
140  }
141 
142  threadnum = static_cast<size_t>(thread_number);
143 
144  return threadnum;
145 }
146 
147 /** sets up CppAD's datastructures for running in multithreading mode
148  *
149  * It must be called once before multithreading is started.
150  * For GCC-compatible compilers, this will happen automatically.
151  */
152 extern "C" SCIP_EXPORT char SCIPexprintCppADInitParallel(void);
153 #ifdef __GNUC__
154 __attribute__((constructor))
155 #endif
156 char SCIPexprintCppADInitParallel(void)
157 {
158  CppAD::thread_alloc::parallel_setup(CPPAD_MAX_NUM_THREADS, in_parallel, thread_num);
159  CppAD::parallel_ad<double>();
160 
161  return 0;
162 }
163 
164 #if !defined(__GNUC__)
165 /** a dummy variable that is initialized to the result of init_parallel
166  *
167  * The purpose is to make sure that init_parallel() is called before any multithreading is started.
168  */
169 static char init_parallel_return = SCIPexprintCppADInitParallel();
170 #endif
171 
172 #endif // SCIP_THREADSAFE
173 
174 using CppAD::AD;
175 
176 class atomic_userexpr;
177 
178 /** expression specific interpreter data */
179 struct SCIP_ExprIntData
180 {
181 public:
182  /** constructor */
183  SCIP_ExprIntData()
184  : val(0.0),
185  need_retape(true),
186  need_retape_always(false),
187  userevalcapability(SCIP_EXPRINTCAPABILITY_ALL),
188  hesrowidxs(NULL),
189  hescolidxs(NULL),
190  hesnnz(0),
191  hesconstant(false)
192  { }
193 
194  /** destructor */
195  ~SCIP_ExprIntData()
196  { }/*lint --e{1540}*/
197 
198  /** gives position of index in varidxs vector */
199  int getVarPos(
200  int varidx /**< variable index to look for */
201  ) const
202  {
203  // varidxs is sorted, so can use binary search functions
204  assert(std::binary_search(varidxs.begin(), varidxs.end(), varidx));
205  return std::lower_bound(varidxs.begin(), varidxs.end(), varidx) - varidxs.begin();
206  }
207 
208  vector< int > varidxs; /**< variable indices used in expression (unique and sorted) */
209  vector< AD<double> > X; /**< vector of dependent variables (same size as varidxs) */
210  vector< AD<double> > Y; /**< result vector (size 1) */
211  CppAD::ADFun<double> f; /**< the function to evaluate as CppAD object */
212 
213  vector<double> x; /**< current values of dependent variables (same size as varidxs) */
214  double val; /**< current function value */
215  bool need_retape; /**< will retaping be required for the next point evaluation? */
216  bool need_retape_always; /**< will retaping be always required? */
217 
218  vector<atomic_userexpr*> userexprs; /**< vector of atomic_userexpr that are created during eval() and need to be kept as long as the tape exists */
219  SCIP_EXPRINTCAPABILITY userevalcapability;/**< (intersection of) capabilities of evaluation routines of user expressions */
220 
221  int* hesrowidxs; /**< row indices of Hessian sparsity: indices are the variables used in expression */
222  int* hescolidxs; /**< column indices of Hessian sparsity: indices are the variables used in expression */
223  vector<SCIP_Real> hesvalues; /**< values of Hessian (a vector<> so we can pass it to CppAD) */
224  int hesnnz; /**< number of nonzeros in Hessian */
225  SCIP_Bool hesconstant; /**< whether Hessian is constant (because expr is at most quadratic) */
226 
227  // Hessian data in CppAD style: indices are 0..n-1 and elements on both lower and upper diagonal of matrix are considered
228  CppAD::local::internal_sparsity<bool>::pattern_type hessparsity_pattern; /**< packed sparsity pattern of Hessian in CppAD-internal form */
229  CppAD::vector<size_t> hessparsity_row; /**< row indices of sparsity pattern of Hessian in CppAD-internal form */
230  CppAD::vector<size_t> hessparsity_col; /**< column indices of sparsity pattern of Hessian in CppAD-internal form */
231  CppAD::sparse_hessian_work heswork; /**< work memory of CppAD for sparse Hessians */
232 };
233 
234 #ifndef NO_CPPAD_USER_ATOMIC
235 
236 /** computes sparsity of jacobian for a univariate function during a forward sweep
237  *
238  * 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.
239  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
240  */
241 static
242 bool univariate_for_sparse_jac(
243  size_t q, /**< number of columns in R */
244  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
245  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
246  )
247 {
248  assert(r.size() == q);
249  assert(s.size() == q);
250 
251  s = r;
252 
253  return true;
254 }
255 
256 /** Computes sparsity of jacobian during a reverse sweep
257  *
258  * 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).
259  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
260  */
261 static
262 bool univariate_rev_sparse_jac(
263  size_t q, /**< number of rows in R */
264  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
265  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
266  )
267 {
268  assert(r.size() == q);
269  assert(s.size() == q);
270 
271  s = r;
272 
273  return true;
274 }
275 
276 /** computes sparsity of hessian during a reverse sweep
277  *
278  * Assume V(x) = (g(f(x)))'' R with f(x) = x^p for a function g:R->R and a matrix R.
279  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
280  */ /*lint -e715*/
281 static
282 bool univariate_rev_sparse_hes(
283  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
284  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
285  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
286  size_t q, /**< number of columns in R, U, and V */
287  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
288  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
289  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
290  )
291 { /*lint --e{439,715}*/ /* @todo take vx into account */
292  assert(r.size() == q);
293  assert(s.size() == 1);
294  assert(t.size() == 1);
295  assert(u.size() == q);
296  assert(v.size() == q);
297 
298  // T(x) = g'(f(x)) * f'(x) = S * f'(x), and f' is not identically 0
299  t[0] = s[0];
300 
301  // V(x) = g''(f(x)) f'(x) f'(x) R + g'(f(x)) f''(x) R466
302  // = f'(x) U + S f''(x) R, with f'(x) and f''(x) not identically 0
303  v = u;
304  if( s[0] )
305  for( size_t j = 0; j < q; ++j )
306  if( r[j] )
307  v[j] = true;
308 
309  return true;
310 }
311 
312 
313 /** Automatic differentiation of x -> x^p, p>=2 integer, as CppAD user-atomic function.
314  *
315  * This class implements forward and reverse operations for the function x -> x^p for use within CppAD.
316  * 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.
317  *
318  * @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)
319  */
320 template<class Type>
321 class atomic_posintpower : public CppAD::atomic_base<Type>
322 {
323 public:
324  atomic_posintpower()
325  : CppAD::atomic_base<Type>("posintpower"),
326  exponent(0)
327  {
328  /* indicate that we want to use bool-based sparsity pattern */
329  this->option(CppAD::atomic_base<Type>::bool_sparsity_enum);
330  }
331 
332 private:
333  /** exponent value for next call to forward or reverse */
334  int exponent;
335 
336  /** stores exponent value corresponding to next call to forward or reverse
337  *
338  * how is this supposed to be threadsafe? (we use only one global instantiation of this class)
339  * TODO using this function is deprecated, do as with atomic_userexpr instead
340  */
341  virtual void set_old(size_t id)
342  {
343  exponent = (int) id;
344  }
345 
346  /** forward sweep of positive integer power
347  *
348  * Given the taylor coefficients for x, we have to compute the taylor coefficients for f(x),
349  * that is, given tx = (x, x', x'', ...), we compute the coefficients ty = (y, y', y'', ...)
350  * in the taylor expansion of f(x) = x^p.
351  * Thus, y = x^p
352  * = tx[0]^p,
353  * y' = p * x^(p-1) * x'
354  * = p * tx[0]^(p-1) * tx[1],
355  * y'' = 1/2 * p * (p-1) * x^(p-2) * x'^2 + p * x^(p-1) * x''
356  * = 1/2 * p * (p-1) * tx[0]^(p-2) * tx[1]^2 + p * tx[0]^(p-1) * tx[2]
357  */
358  bool forward(
359  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
360  size_t p, /**< highest order Taylor coefficient that we are evaluating */
361  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
362  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
363  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
364  CppAD::vector<Type>& ty /**< vector to store taylor coefficients of y */
365  )
366  {
367  assert(exponent > 1);
368  assert(tx.size() >= p+1);
369  assert(ty.size() >= p+1);
370  assert(q <= p);
371 
372  if( vx.size() > 0 )
373  {
374  assert(vx.size() == 1);
375  assert(vy.size() == 1);
376  assert(p == 0);
377 
378  vy[0] = vx[0];
379  }
380 
381  if( q == 0 /* q <= 0 && 0 <= p */ )
382  {
383  ty[0] = CppAD::pow(tx[0], exponent);
384  }
385 
386  if( q <= 1 && 1 <= p )
387  {
388  ty[1] = CppAD::pow(tx[0], exponent-1) * tx[1];
389  ty[1] *= double(exponent);
390  }
391 
392  if( q <= 2 && 2 <= p )
393  {
394  if( exponent > 2 )
395  {
396  // ty[2] = 1/2 * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1] * tx[1] + exponent * pow(tx[0], exponent-1) * tx[2];
397  ty[2] = CppAD::pow(tx[0], exponent-2) * tx[1] * tx[1];
398  ty[2] *= (exponent-1) / 2.0;
399  ty[2] += CppAD::pow(tx[0], exponent-1) * tx[2];
400  ty[2] *= exponent;
401  }
402  else
403  {
404  assert(exponent == 2);
405  // ty[2] = 1/2 * exponent * tx[1] * tx[1] + exponent * tx[0] * tx[2];
406  ty[2] = tx[1] * tx[1] + 2.0 * tx[0] * tx[2];
407  }
408  }
409 
410  /* higher order derivatives not implemented */
411  if( p > 2 )
412  return false;
413 
414  return true;
415  }
416 
417  /** reverse sweep of positive integer power
418  *
419  * Assume y(x) is a function of the taylor coefficients of f(x) = x^p for x, i.e.,
420  * y(x) = [ x^p, p * x^(p-1) * x', p * (p-1) * x^(p-2) * x'^2 + p * x^(p-1) * x'', ... ].
421  * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
422  * 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.
423  * That is, we have to compute
424  *\f$
425  * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
426  * = \sum_{i=0}^k (\partial g / \partial y_i) * (\partial y_i / \partial x^[l])
427  * = \sum_{i=0}^k py[i] * (\partial y_i / \partial x^[l])
428  * \f$
429  *
430  * For k = 0, this means
431  *\f$
432  * px[0] = py[0] * (\partial y_0 / \partial x^[0])
433  * = py[0] * (\partial x^p / \partial x)
434  * = py[0] * p * tx[0]^(p-1)
435  *\f$
436  *
437  * For k = 1, this means
438  * \f$
439  * px[0] = py[0] * (\partial y_0 / \partial x^[0]) + py[1] * (\partial y_1 / \partial x^[0])
440  * = py[0] * (\partial x^p / \partial x) + py[1] * (\partial (p * x^(p-1) * x') / \partial x)
441  * = py[0] * p * tx[0]^(p-1) + py[1] * p * (p-1) * tx[0]^(p-2) * tx[1]
442  * px[1] = py[0] * (\partial y_0 / \partial x^[1]) + py[1] * (\partial y_1 / \partial x^[1])
443  * = py[0] * (\partial x^p / \partial x') + py[1] * (\partial (p * x^(p-1) x') / \partial x')
444  * = py[0] * 0 + py[1] * p * tx[0]^(p-1)
445  * \f$
446  */ /*lint -e715*/
447  bool reverse(
448  size_t p, /**< highest order Taylor coefficient that we are evaluating */
449  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
450  const CppAD::vector<Type>& ty, /**< values for taylor coefficients of y */
451  CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
452  const CppAD::vector<Type>& py /**< values for partial derivatives of g(x) w.r.t. y */
453  )
454  { /*lint --e{715}*/
455  assert(exponent > 1);
456  assert(px.size() >= p+1);
457  assert(py.size() >= p+1);
458  assert(tx.size() >= p+1);
459 
460  switch( p )
461  {
462  case 0:
463  // px[0] = py[0] * exponent * pow(tx[0], exponent-1);
464  px[0] = py[0] * CppAD::pow(tx[0], exponent-1);
465  px[0] *= exponent;
466  break;
467 
468  case 1:
469  // px[0] = py[0] * exponent * pow(tx[0], exponent-1) + py[1] * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1];
470  px[0] = py[1] * tx[1] * CppAD::pow(tx[0], exponent-2);
471  px[0] *= exponent-1;
472  px[0] += py[0] * CppAD::pow(tx[0], exponent-1);
473  px[0] *= exponent;
474  // px[1] = py[1] * exponent * pow(tx[0], exponent-1);
475  px[1] = py[1] * CppAD::pow(tx[0], exponent-1);
476  px[1] *= exponent;
477  break;
478 
479  default:
480  return false;
481  }
482 
483  return true;
484  }
485 
486  using CppAD::atomic_base<Type>::for_sparse_jac;
487 
488  /** computes sparsity of jacobian during a forward sweep
489  *
490  * 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.
491  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
492  */
493  bool for_sparse_jac(
494  size_t q, /**< number of columns in R */
495  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
496  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
497  )
498  {
499  return univariate_for_sparse_jac(q, r, s);
500  }
501 
502  using CppAD::atomic_base<Type>::rev_sparse_jac;
503 
504  /** computes sparsity of jacobian during a reverse sweep
505  *
506  * 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).
507  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
508  */
509  bool rev_sparse_jac(
510  size_t q, /**< number of rows in R */
511  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
512  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
513  )
514  {
515  return univariate_rev_sparse_jac(q, r, s);
516  }
517 
518  using CppAD::atomic_base<Type>::rev_sparse_hes;
519 
520  /** computes sparsity of hessian during a reverse sweep
521  *
522  * Assume V(x) = (g(f(x)))'' R with f(x) = x^p for a function g:R->R and a matrix R.
523  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
524  */
525  bool rev_sparse_hes(
526  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
527  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
528  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
529  size_t q, /**< number of columns in R, U, and V */
530  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
531  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
532  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
533  )
534  {
535  return univariate_rev_sparse_hes(vx, s, t, q, r, u, v);
536  }
537 };
538 
539 /** power function with natural exponents */
540 template<class Type>
541 static
542 void posintpower(
543  const vector<Type>& in, /**< vector which first argument is base */
544  vector<Type>& out, /**< vector where to store result in first argument */
545  size_t exponent /**< exponent */
546  )
547 {
548  static atomic_posintpower<typename Type::value_type> pip;
549  pip(in, out, exponent);
550 }
551 
552 #else
553 
554 /** power function with natural exponents */
555 template<class Type>
557  const vector<Type>& in, /**< vector which first argument is base */
558  vector<Type>& out, /**< vector where to store result in first argument */
559  size_t exponent /**< exponent */
560  )
561 {
562  out[0] = pow(in[0], (int)exponent);
563 }
564 
565 #endif
566 
567 
568 #ifndef NO_CPPAD_USER_ATOMIC
569 
570 /** sign of a value (-1 or +1)
571  *
572  * 0.0 has sign +1
573  */
574 #define SIGN(x) ((x) >= 0.0 ? 1.0 : -1.0)
575 
576 /** Automatic differentiation of x -> sign(x)abs(x)^p, p>=1, as CppAD user-atomic function.
577  *
578  * This class implements forward and reverse operations for the function x -> sign(x)abs(x)^p for use within CppAD.
579  * While we otherwise would have to use discontinuous sign and abs functions, our own implementation allows to provide
580  * a continuously differentiable function.
581  *
582  * @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)
583  */
584 template<class Type>
585 class atomic_signpower : public CppAD::atomic_base<Type>
586 {
587 public:
588  atomic_signpower()
589  : CppAD::atomic_base<Type>("signpower"),
590  exponent(0.0)
591  {
592  /* indicate that we want to use bool-based sparsity pattern */
593  this->option(CppAD::atomic_base<Type>::bool_sparsity_enum);
594  }
595 
596 private:
597  /** exponent for use in next call to forward or reverse */
598  SCIP_Real exponent;
599 
600  /** stores exponent corresponding to next call to forward or reverse
601  *
602  * How is this supposed to be threadsafe? (we use only one global instantiation of this class)
603  * TODO using this function is deprecated, do as with atomic_userexpr instead
604  */
605  virtual void set_old(size_t id)
606  {
607  exponent = SCIPgetExponentExprPow((SCIP_EXPR*)(void*)id);
608  }
609 
610  /** forward sweep of signpower
611  *
612  * Given the taylor coefficients for x, we have to compute the taylor coefficients for f(x),
613  * that is, given tx = (x, x', x'', ...), we compute the coefficients ty = (y, y', y'', ...)
614  * in the taylor expansion of f(x) = sign(x)abs(x)^p.
615  * Thus, y = sign(x)abs(x)^p
616  * = sign(tx[0])abs(tx[0])^p,
617  * y' = p * abs(x)^(p-1) * x'
618  * = p * abs(tx[0])^(p-1) * tx[1],
619  * y'' = 1/2 * p * (p-1) * sign(x) * abs(x)^(p-2) * x'^2 + p * abs(x)^(p-1) * x''
620  * = 1/2 * p * (p-1) * sign(tx[0]) * abs(tx[0])^(p-2) * tx[1]^2 + p * abs(tx[0])^(p-1) * tx[2]
621  */
622  bool forward(
623  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
624  size_t p, /**< highest order Taylor coefficient that we are evaluating */
625  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
626  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
627  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
628  CppAD::vector<Type>& ty /**< vector to store taylor coefficients of y */
629  )
630  {
631  assert(exponent > 0.0);
632  assert(tx.size() >= p+1);
633  assert(ty.size() >= p+1);
634  assert(q <= p);
635 
636  if( vx.size() > 0 )
637  {
638  assert(vx.size() == 1);
639  assert(vy.size() == 1);
640  assert(p == 0);
641 
642  vy[0] = vx[0];
643  }
644 
645  if( q == 0 /* q <= 0 && 0 <= p */ )
646  {
647  ty[0] = SIGN(tx[0]) * pow(REALABS(tx[0]), exponent);
648  }
649 
650  if( q <= 1 && 1 <= p )
651  {
652  ty[1] = pow(REALABS(tx[0]), exponent - 1.0) * tx[1];
653  ty[1] *= exponent;
654  }
655 
656  if( q <= 2 && 2 <= p )
657  {
658  if( exponent != 2.0 )
659  {
660  ty[2] = SIGN(tx[0]) * pow(REALABS(tx[0]), exponent - 2.0) * tx[1] * tx[1];
661  ty[2] *= (exponent - 1.0) / 2.0;
662  ty[2] += pow(REALABS(tx[0]), exponent - 1.0) * tx[2];
663  ty[2] *= exponent;
664  }
665  else
666  {
667  // y'' = 2 (1/2 * sign(x) * x'^2 + |x|*x'') = sign(tx[0]) * tx[1]^2 + 2 * abs(tx[0]) * tx[2]
668  ty[2] = SIGN(tx[0]) * tx[1] * tx[1];
669  ty[2] += 2.0 * REALABS(tx[0]) * tx[2];
670  }
671  }
672 
673  /* higher order derivatives not implemented */
674  if( p > 2 )
675  return false;
676 
677  return true;
678  }
679 
680  /** reverse sweep of signpower
681  *
682  * Assume y(x) is a function of the taylor coefficients of f(x) = sign(x)|x|^p for x, i.e.,
683  * y(x) = [ f(x), f'(x), f''(x), ... ].
684  * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
685  * 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.
686  * That is, we have to compute
687  *\f$
688  * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
689  * = \sum_{i=0}^k (\partial g / \partial y_i) * (\partial y_i / \partial x^[l])
690  * = \sum_{i=0}^k py[i] * (\partial y_i / \partial x^[l])
691  *\f$
692  *
693  * For k = 0, this means
694  *\f$
695  * px[0] = py[0] * (\partial y_0 / \partial x^[0])
696  * = py[0] * (\partial f(x) / \partial x)
697  * = py[0] * p * abs(tx[0])^(p-1)
698  * \f$
699  *
700  * For k = 1, this means
701  *\f$
702  * px[0] = py[0] * (\partial y_0 / \partial x^[0]) + py[1] * (\partial y_1 / \partial x^[0])
703  * = py[0] * (\partial f(x) / \partial x) + py[1] * (\partial f'(x) / \partial x)
704  * = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1]
705  * px[1] = py[0] * (\partial y_0 / \partial x^[1]) + py[1] * (\partial y_1 / \partial x^[1])
706  * = py[0] * (\partial f(x) / \partial x') + py[1] * (\partial f'(x) / \partial x')
707  * = py[0] * 0 + py[1] * p * abs(tx[0])^(p-1)
708  * \f$
709  */
710  bool reverse(
711  size_t p, /**< highest order Taylor coefficient that we are evaluating */
712  const CppAD::vector<Type>& tx, /**< values for taylor coefficients of x */
713  const CppAD::vector<Type>& ty, /**< values for taylor coefficients of y */
714  CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
715  const CppAD::vector<Type>& py /**< values for partial derivatives of g(x) w.r.t. y */
716  )
717  { /*lint --e{715}*/
718  assert(exponent > 1);
719  assert(px.size() >= p+1);
720  assert(py.size() >= p+1);
721  assert(tx.size() >= p+1);
722 
723  switch( p )
724  {
725  case 0:
726  // px[0] = py[0] * p * pow(abs(tx[0]), p-1);
727  px[0] = py[0] * pow(REALABS(tx[0]), exponent - 1.0);
728  px[0] *= p;
729  break;
730 
731  case 1:
732  if( exponent != 2.0 )
733  {
734  // 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]
735  px[0] = py[1] * tx[1] * pow(REALABS(tx[0]), exponent - 2.0) * SIGN(tx[0]);
736  px[0] *= exponent - 1.0;
737  px[0] += py[0] * pow(REALABS(tx[0]), exponent - 1.0);
738  px[0] *= exponent;
739  // px[1] = py[1] * p * abs(tx[0])^(p-1)
740  px[1] = py[1] * pow(REALABS(tx[0]), exponent - 1.0);
741  px[1] *= exponent;
742  }
743  else
744  {
745  // px[0] = py[0] * 2.0 * abs(tx[0]) + py[1] * 2.0 * sign(tx[0]) * tx[1]
746  px[0] = py[1] * tx[1] * SIGN(tx[0]);
747  px[0] += py[0] * REALABS(tx[0]);
748  px[0] *= 2.0;
749  // px[1] = py[1] * 2.0 * abs(tx[0])
750  px[1] = py[1] * REALABS(tx[0]);
751  px[1] *= 2.0;
752  }
753  break;
754 
755  default:
756  return false;
757  }
758 
759  return true;
760  }
761 
762  using CppAD::atomic_base<Type>::for_sparse_jac;
763 
764  /** computes sparsity of jacobian during a forward sweep
765  *
766  * 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.
767  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
768  */
769  bool for_sparse_jac(
770  size_t q, /**< number of columns in R */
771  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
772  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
773  )
774  {
775  return univariate_for_sparse_jac(q, r, s);
776  }
777 
778  using CppAD::atomic_base<Type>::rev_sparse_jac;
779 
780  /** computes sparsity of jacobian during a reverse sweep
781  *
782  * 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).
783  * Since f'(x) is dense, the sparsity of S will be the sparsity of R.
784  */
785  bool rev_sparse_jac(
786  size_t q, /**< number of rows in R */
787  const CppAD::vector<bool>& r, /**< sparsity of R, rowwise */
788  CppAD::vector<bool>& s /**< vector to store sparsity of S, rowwise */
789  )
790  {
791  return univariate_rev_sparse_jac(q, r, s);
792  }
793 
794  using CppAD::atomic_base<Type>::rev_sparse_hes;
795 
796  /** computes sparsity of hessian during a reverse sweep
797  *
798  * 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.
799  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
800  */
801  bool rev_sparse_hes(
802  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
803  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
804  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
805  size_t q, /**< number of columns in S and R */
806  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
807  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
808  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
809  )
810  {
811  return univariate_rev_sparse_hes(vx, s, t, q, r, u, v);
812  }
813 
814 };
815 
816 /** template for evaluation for signpower operator */
817 template<class Type>
818 static
819 void evalSignPower(
820  Type& resultant, /**< resultant */
821  const Type& arg, /**< operand */
822  SCIP_EXPR* expr /**< expression that holds the exponent */
823  )
824 {
825  vector<Type> in(1, arg);
826  vector<Type> out(1);
827 
828  static atomic_signpower<typename Type::value_type> sp;
829  sp(in, out, (size_t)(void*)expr);
830 
831  resultant = out[0];
832  return;
833 }
834 
835 #else
836 
837 /** specialization of signpower evaluation for real numbers */
838 template<class Type>
839 static
841  CppAD::AD<Type>& resultant, /**< resultant */
842  const CppAD::AD<Type>& arg, /**< operand */
843  SCIP_EXPR* expr /**< expression that holds the exponent */
844  )
845 {
846  AD<Type> adzero(0.);
847  SCIP_Real exponent;
848 
849  exponent = SCIPgetExponentExprPow(expr);
850  assert(exponent >= 1.0);
851 
852  if( EPSISINT(exponent, 0.0) )
853  {
854  resultant = CppAD::CondExpGe(arg, adzero, pow(arg, (int)exponent), -pow(-arg, (int)exponent));
855  }
856  else
857  {
858  /* pow(0,fractional>1) is not differential in the CppAD world (https://github.com/coin-or/CppAD/discussions/93)
859  * this works around this by evaluating pow(eps,exponent) in this case
860  */
861  resultant = CppAD::CondExpEq(arg, adzero, pow(arg+std::numeric_limits<SCIP_Real>::epsilon(), exponent)-pow(std::numeric_limits<SCIP_Real>::epsilon(), exponent),
862  CppAD::CondExpGe(arg, adzero, pow(arg, exponent), -pow(-arg, exponent)));
863  }
864 }
865 #endif
866 
867 /** Automatic differentiation of user expression as CppAD user-atomic function.
868  *
869  * This class implements forward and reverse operations for a function given by a user expression for use within CppAD.
870  */
871 class atomic_userexpr : public CppAD::atomic_base<SCIP_Real>
872 {
873 public:
875  SCIP* scip_, /**< SCIP data structure */
876  SCIP_EXPR* expr_ /**< expression to use */
877  )
878  : CppAD::atomic_base<SCIP_Real>(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr_)), CppAD::atomic_base<SCIP_Real>::bool_sparsity_enum),
879  scip(scip_),
880  expr(expr_)
881  { }
882 
883 private:
884  /** SCIP data structure */
885  SCIP* scip;
886 
887  /** expression */
888  SCIP_EXPR* expr;
889 
890  /** forward sweep of userexpr
891  *
892  * We follow http://www.coin-or.org/CppAD/Doc/atomic_forward.xml
893  * Note, that p and q are interchanged!
894  *
895  * For a scalar variable t, let
896  * Y(t) = f(X(t))
897  * X(t) = x^0 + x^1 t^1 + ... + x^p t^p
898  * where for x^i the i an index, while for t^i the i is an exponent.
899  * Thus, x^k = 1/k! X^(k) (0), where X^(k)(.) denotes the k-th derivative.
900  *
901  * Next, let y^k = 1/k! Y^(k)(0) be the k'th taylor coefficient of Y. Thus,
902  * y^0 = Y^(0)(0) = Y(0) = f(X(0)) = f(x^0)
903  * y^1 = Y^(1)(0) = Y'(0) = f'(X(0)) * X'(0) = f'(x^0) * x^1
904  * y^2 = 1/2 Y^(2)(0) = 1/2 Y''(0) = 1/2 X'(0) * f''(X(0)) X'(0) + 1/2 * f'(X(0)) * X''(0) = 1/2 x^1 * f''(x^0) * x^1 + f'(x^0) * x^2
905  *
906  * As x^k = (tx[k], tx[(p+1)+k], tx[2*(p+1)+k], ..., tx[n*(p+1)+k], we get
907  * ty[0] = y^0 = f(x^0) = f(tx[{1..n}*(p+1)])
908  * ty[1] = y^1 = f'(x^0) * tx[{1..n}*(p+1)+1] = sum(i=1..n, grad[i] * tx[i*(p+1)+1]), where grad = f'(x^0)
909  * ty[2] = 1/2 sum(i,j=1..n, x[i*(p+1)+1] * x[j*(p+1)+q] * hessian[i,j]) + sum(i=1..n, grad[i] * x[i*(p+1)+2])
910  */
911  /* cppcheck-suppress unusedPrivateFunction */
912  bool forward(
913  size_t q, /**< lowest order Taylor coefficient that we are evaluating */
914  size_t p, /**< highest order Taylor coefficient that we are evaluating */
915  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
916  CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
917  const CppAD::vector<SCIP_Real>& tx, /**< values for taylor coefficients of x */
918  CppAD::vector<SCIP_Real>& ty /**< vector to store taylor coefficients of y */
919  )
920  {
921  assert(scip != NULL);
922  assert(expr != NULL);
923  assert(ty.size() == p+1);
924  assert(q <= p);
925 
926  size_t n = tx.size() / (p+1);
927  assert(n == (size_t)SCIPexprGetNChildren(expr)); /*lint !e571*/
928  assert(n >= 1);
929 
930  SCIPdebugMsg(scip, "expr_%s:forward, q=%zd, p=%zd\n", SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), q, p);
931 
932  if( vx.size() > 0 )
933  {
934  assert(vx.size() == n);
935  assert(vy.size() == 1);
936  assert(p == 0);
937 
938  /* y_0 is a variable if at least one of the x_i is a variable */
939  vy[0] = false;
940  for( size_t i = 0; i < n; ++i )
941  if( vx[i] )
942  {
943  vy[0] = true;
944  break;
945  }
946  }
947 
948  if( p == 0 )
949  {
950  /* only eval requested
951  * arguments are in tx[i * (p+1) + 0], i == 0..n-1, that is tx[0]..tx[n-1]
952  */
953  if( SCIPcallExprEval(scip, expr, const_cast<double*>(tx.data()), &ty[0]) != SCIP_OKAY )
954  return false;
955 
956  if( ty[0] == SCIP_INVALID )
958 
959  return true;
960  }
961 
962  if( p == 1 )
963  {
964  /* eval and forward-diff requested
965  *
966  * point to eval is in tx[i * (p+1) + 0], i == 0..n-1
967  * direction is in tx[i * (p+1) + 1], i == 0..n-1
968  */
969  SCIP_Real* x = new SCIP_Real[n];
970  SCIP_Real* dir = new SCIP_Real[n];
971  for( size_t i = 0; i < n; ++i )
972  {
973  x[i] = tx[i * (p+1) + 0]; /*lint !e835*/
974  dir[i] = tx[i * (p+1) + 1]; /*lint !e835*/
975  }
976 
977  SCIP_RETCODE rc = SCIPcallExprEvalFwdiff(scip, expr, x, dir, &ty[0], &ty[1]);
978 
979  if( ty[0] == SCIP_INVALID )
981  if( ty[1] == SCIP_INVALID )
983 
984  delete[] dir;
985  delete[] x;
986 
987  return rc == SCIP_OKAY;
988  }
989 
990  /* higher order derivatives not implemented yet (TODO) */
991  SCIPABORT();
992  return false;
993 
994 #if SCIP_DISABLED_CODE
995  /* cppcheck-suppress unreachableCode */
996  SCIP_Real* gradient = NULL;
997  SCIP_Real* hessian = NULL;
998 
999  if( q <= 2 && 1 <= p )
1000  gradient = new SCIP_Real[n];
1001  if( q <= 2 && 2 <= p )
1002  hessian = new SCIP_Real[n*n];
1003 
1004  if( exprEvalUser(expr, x, ty[0], gradient, hessian) != SCIP_OKAY )
1005  {
1006  delete[] x;
1007  delete[] gradient;
1008  delete[] hessian;
1009  return false;
1010  }
1011 
1012  if( gradient != NULL )
1013  {
1014  ty[1] = 0.0;
1015  for( size_t i = 0; i < n; ++i )
1016  ty[1] += gradient[i] * tx[i * (p+1) + 1];
1017  }
1018 
1019  if( hessian != NULL )
1020  {
1021  assert(gradient != NULL);
1022 
1023  ty[2] = 0.0;
1024  for( size_t i = 0; i < n; ++i )
1025  {
1026  for( size_t j = 0; j < n; ++j )
1027  ty[2] += 0.5 * hessian[i*n+j] * tx[i * (p+1) + 1] * tx[j * (p+1) + 1];
1028 
1029  ty[2] += gradient[i] * tx[i * (p+1) + 2];
1030  }
1031  }
1032 
1033  delete[] x;
1034  delete[] gradient;
1035  delete[] hessian;
1036 
1037  /* higher order derivatives not implemented */
1038  if( p > 2 )
1039  return false;
1040 
1041  return true;
1042 #endif
1043  }
1044 
1045  /** reverse sweep of userexpr
1046  *
1047  * We follow http://www.coin-or.org/CppAD/Doc/atomic_reverse.xml
1048  * Note, that there q is our p.
1049  *
1050  * For a scalar variable t, let
1051  * Y(t) = f(X(t))
1052  * X(t) = x^0 + x^1 t^1 + ... + x^p t^p
1053  * where for x^i the i an index, while for t^i the i is an exponent.
1054  * Thus, x^k = 1/k! X^(k) (0), where X^(k)(.) denotes the k-th derivative.
1055  *
1056  * Next, let y^k = 1/k! Y^(k)(0) be the k'th taylor coefficient of Y. Thus,
1057  * Y(t) = y^0 + y^1 t^1 + y^2 t^2 + ...
1058  * y^0, y^1, ... are the taylor coefficients of f(x).
1059  *
1060  * Further, let F(x^0,..,x^p) by given as F^k(x) = y^k. Thus,
1061  * F^0(x) = y^0 = Y^(0)(0) = f(x^0)
1062  * F^1(x) = y^1 = Y^(1)(0) = f'(x^0) * x^1
1063  * F^2(x) = y^2 = 1/2 Y''(0) = 1/2 x^1 f''(x^0) x^1 + f'(x^0) x^2
1064  *
1065  * Given functions G: R^(p+1) -> R and H: R^(n*(p+1)) -> R, where H(x^0, x^1, .., x^p) = G(F(x^0,..,x^p)),
1066  * we have to return the value of \f$\partial H / \partial x^l, l = 0..p,\f$ in px. Therefor,
1067  * \f$
1068  * px^l = \partial H / \partial x^l
1069  * = sum(k=0..p, (\partial G / \partial y^k) * (\partial y^k / \partial x^l)
1070  * = sum(k=0..p, py[k] * (\partial F^k / \partial x^l)
1071  * \f$
1072  *
1073  * For p = 0, this means
1074  * \f$
1075  * px^0 = py[0] * \partial F^0 / \partial x^0
1076  * = py[0] * \partial f(x^0) / \partial x^0
1077  * = py[0] * f'(x^0)
1078  * \f$
1079  *
1080  * For p = 1, this means (for l = 0):
1081  * \f[
1082  * px^0 = py[0] * \partial F^0 / \partial x^0 + py[1] * \partial F^1 / \partial x^0
1083  * = py[0] * \partial f(x^0) / \partial x^0 + py[1] * \partial (f'(x^0) * x^1) / \partial x^0
1084  * = py[0] * f'(x^0) + py[1] * f''(x^0) * x^1
1085  * \f]
1086  * and (for l=1):
1087  * \[
1088  * px^1 = py[0] * \partial F^0 / \partial x^1 + py[1] * \partial F^1 / \partial x^1
1089  * = py[0] * \partial f(x^0) / \partial x^1 + py[1] * \partial (f'(x^0) * x^1) / \partial x^0
1090  * = py[0] * 0 + py[1] * f'(x^0)
1091  * \f]
1092  *
1093  * As x^k = (tx[k], tx[(p+1)+k], tx[2*(p+1)+k], ..., tx[n*(p+1)+k] and
1094  * px^k = (px[k], px[(p+1)+k], px[2*(p+1)+k], ..., px[n*(p+1)+k], we get
1095  * for p = 0:
1096  * px[i] = (px^0)_i = py[0] * grad[i]
1097  * for p = 1:
1098  * px[i*2+0] = (px^0)_i = py[0] * grad[i] + py[1] * sum(j, hessian[j,i] * tx[j*2+1])
1099  * px[i*2+1] = (px^1)_i = py[1] * grad[i]
1100  */
1101  /* cppcheck-suppress unusedPrivateFunction */
1102  bool reverse(
1103  size_t p, /**< highest order Taylor coefficient that we are evaluating */
1104  const CppAD::vector<SCIP_Real>& tx, /**< values for taylor coefficients of x */
1105  const CppAD::vector<SCIP_Real>& ty, /**< values for taylor coefficients of y */
1106  CppAD::vector<SCIP_Real>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
1107  const CppAD::vector<SCIP_Real>& py /**< values for partial derivatives of g(x) w.r.t. y */
1108  )
1109  {
1110  assert(expr != NULL);
1111  assert(px.size() == tx.size());
1112  assert(py.size() == p+1);
1113 
1114  // TODO implement this again and then remove check for userexprs in SCIPexprintGrad
1115  SCIPABORT();
1116  return false;
1117 
1118 #ifdef SCIP_DISABLED_CODE
1119  /* cppcheck-suppress unreachableCode */
1120  size_t n = tx.size() / (p+1);
1121  assert(n == (size_t)SCIPexprGetNChildren(expr)); /*lint !e571*/
1122  assert(n >= 1);
1123 
1124  SCIP_Real* x = new SCIP_Real[n];
1125  SCIP_Real funcval;
1126  SCIP_Real* gradient = new SCIP_Real[n];
1127  SCIP_Real* hessian = NULL;
1128 
1129  if( p == 1 )
1130  hessian = new SCIP_Real[n*n];
1131 
1132  for( size_t i = 0; i < n; ++i )
1133  x[i] = tx[i * (p+1) + 0]; /*lint !e835*/
1134 
1135  if( exprEvalUser(expr, x, funcval, gradient, hessian) != SCIP_OKAY )
1136  {
1137  delete[] x;
1138  delete[] gradient;
1139  delete[] hessian;
1140  return false;
1141  }
1142 
1143  switch( p )
1144  {
1145  case 0:
1146  // px[j] = (px^0)_j = py[0] * grad[j]
1147  for( size_t i = 0; i < n; ++i )
1148  px[i] = py[0] * gradient[i];
1149  break;
1150 
1151  case 1:
1152  // px[i*2+0] = (px^0)_i = py[0] * grad[i] + py[1] * sum(j, hessian[j,i] * tx[j*2+1])
1153  // px[i*2+1] = (px^1)_i = py[1] * grad[i]
1154  assert(hessian != NULL);
1155  for( size_t i = 0; i < n; ++i )
1156  {
1157  px[i*2+0] = py[0] * gradient[i]; /*lint !e835*/
1158  for( size_t j = 0; j < n; ++j )
1159  px[i*2+0] += py[1] * hessian[i+n*j] * tx[j*2+1]; /*lint !e835*/
1160 
1161  px[i*2+1] = py[1] * gradient[i];
1162  }
1163  break;
1164 
1165  default:
1166  return false;
1167  }
1168 
1169  return true;
1170 #endif
1171  } /*lint !e715*/
1172 
1173  using CppAD::atomic_base<SCIP_Real>::for_sparse_jac;
1174 
1175  /** computes sparsity of jacobian during a forward sweep
1176  * 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.
1177  * Since we assume f'(x) to be dense, the sparsity of S will be the sparsity of R.
1178  */
1179  /* cppcheck-suppress unusedPrivateFunction */
1180  bool for_sparse_jac(
1181  size_t q, /**< number of columns in R */
1182  const CppAD::vector<bool>& r, /**< sparsity of R, columnwise */
1183  CppAD::vector<bool>& s /**< vector to store sparsity of S, columnwise */
1184  )
1185  {
1186  assert(expr != NULL);
1187  assert(s.size() == q);
1188 
1189  size_t n = r.size() / q;
1190  assert(n == (size_t)SCIPexprGetNChildren(expr)); /*lint !e571*/
1191 
1192  // sparsity for S(x) = f'(x) * R
1193  for( size_t j = 0; j < q; j++ )
1194  {
1195  s[j] = false;
1196  for( size_t i = 0; i < n; i++ )
1197  s[j] |= (bool)r[i * q + j]; /*lint !e1786*/
1198  }
1199 
1200  return true;
1201  }
1202 
1203  using CppAD::atomic_base<SCIP_Real>::rev_sparse_jac;
1204 
1205  /** computes sparsity of jacobian during a reverse sweep
1206  * For a q x 1 matrix S, we have to return the sparsity pattern of the q x 1 matrix R(x) = S * f'(x).
1207  * Since we assume f'(x) to be dense, the sparsity of R will be the sparsity of S.
1208  */
1209  /* cppcheck-suppress unusedPrivateFunction */
1210  bool rev_sparse_jac(
1211  size_t q, /**< number of rows in R */
1212  const CppAD::vector<bool>& rt, /**< sparsity of R, rowwise */
1213  CppAD::vector<bool>& st /**< vector to store sparsity of S, rowwise */
1214  )
1215  {
1216  assert(expr != NULL);
1217  assert(rt.size() == q);
1218 
1219  size_t n = st.size() / q;
1220  assert(n == (size_t)SCIPexprGetNChildren(expr));
1221 
1222  // sparsity for S(x)^T = f'(x)^T * R^T
1223  for( size_t j = 0; j < q; j++ )
1224  for( size_t i = 0; i < n; i++ )
1225  st[i * q + j] = rt[j];
1226 
1227  return true;
1228  }
1229 
1230  using CppAD::atomic_base<SCIP_Real>::rev_sparse_hes;
1231 
1232  /** computes sparsity of hessian during a reverse sweep
1233  * Assume V(x) = (g(f(x)))'' R for a function g:R->R and a matrix R.
1234  * we have to specify the sparsity pattern of V(x) and T(x) = (g(f(x)))'.
1235  */
1236  /* cppcheck-suppress unusedPrivateFunction */
1237  bool rev_sparse_hes(
1238  const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
1239  const CppAD::vector<bool>& s, /**< sparsity pattern of S = g'(y) */
1240  CppAD::vector<bool>& t, /**< vector to store sparsity pattern of T(x) = (g(f(x)))' */
1241  size_t q, /**< number of columns in S and R */
1242  const CppAD::vector<bool>& r, /**< sparsity pattern of R */
1243  const CppAD::vector<bool>& u, /**< sparsity pattern of U(x) = g''(f(x)) f'(x) R */
1244  CppAD::vector<bool>& v /**< vector to store sparsity pattern of V(x) = (g(f(x)))'' R */
1245  )
1246  {
1247  assert(expr != NULL);
1248  size_t n = vx.size();
1249  assert((size_t)SCIPexprGetNChildren(expr) == n);
1250  assert(s.size() == 1);
1251  assert(t.size() == n);
1252  assert(r.size() == n * q);
1253  assert(u.size() == q);
1254  assert(v.size() == n * q);
1255 
1256  size_t i, j, k;
1257 
1258  // sparsity for T(x) = S(x) * f'(x)
1259  for( i = 0; i < n; ++i )
1260  t[i] = s[0];
1261 
1262  // V(x) = f'(x)^T * g''(y) * f'(x) * R + g'(y) * f''(x) * R
1263  // U(x) = g''(y) * f'(x) * R
1264  // S(x) = g'(y)
1265 
1266  // back propagate the sparsity for U
1267  for( j = 0; j < q; j++ )
1268  for( i = 0; i < n; i++ )
1269  v[ i * q + j] = u[j];
1270 
1271  // include forward Jacobian sparsity in Hessian sparsity
1272  // sparsity for g'(y) * f''(x) * R (Note f''(x) is assumed to be dense)
1273  if( s[0] )
1274  for( j = 0; j < q; j++ )
1275  for( i = 0; i < n; i++ )
1276  for( k = 0; k < n; ++k )
1277  v[ i * q + j] |= (bool) r[ k * q + j];
1278 
1279  return true;
1280  }
1281 };
1282 
1283 
1284 /** integer power operation for arbitrary integer exponents */
1285 template<class Type>
1286 static
1288  Type& resultant, /**< resultant */
1289  const Type& arg, /**< operand */
1290  const int exponent /**< exponent */
1291  )
1292 {
1293  if( exponent > 1 )
1294  {
1295  vector<Type> in(1, arg);
1296  vector<Type> out(1);
1297 
1298  posintpower(in, out, exponent);
1299 
1300  resultant = out[0];
1301  return;
1302  }
1303 
1304  if( exponent < -1 )
1305  {
1306  vector<Type> in(1, arg);
1307  vector<Type> out(1);
1308 
1309  posintpower(in, out, -exponent);
1310 
1311  resultant = Type(1.0)/out[0];
1312  return;
1313  }
1314 
1315  if( exponent == 1 )
1316  {
1317  resultant = arg;
1318  return;
1319  }
1320 
1321  if( exponent == 0 )
1322  {
1323  resultant = Type(1.0);
1324  return;
1325  }
1326 
1327  assert(exponent == -1);
1328  resultant = Type(1.0)/arg;
1329 }
1330 
1331 /** CppAD compatible evaluation of an expression for given arguments */
1332 template<class Type>
1333 static
1335  SCIP* scip, /**< SCIP data structure */
1336  SCIP_EXPR* expr, /**< expression */
1337  SCIP_EXPRINTDATA* exprintdata, /**< interpreter data for root expression */
1338  const vector<Type>& x, /**< values of variables */
1339  Type& val /**< buffer to store expression value */
1340  )
1341 {
1342  Type* buf = NULL;
1343 
1344  assert(expr != NULL);
1345 
1346  // TODO this should iterate instead of using recursion
1347  // but the iterdata wouldn't work to hold Type at the moment
1348  // they could hold Type*, but then we need to alloc small portions all the time
1349  // or we have a big Type-array outside and point to it in iterdata
1350 
1351  if( SCIPisExprVaridx(scip, expr) )
1352  {
1353  val = x[exprintdata->getVarPos(SCIPgetIndexExprVaridx(expr))];
1354  return SCIP_OKAY;
1355  }
1356  if( SCIPisExprValue(scip, expr) )
1357  {
1358  val = SCIPgetValueExprValue(expr);
1359  return SCIP_OKAY;
1360  }
1361 
1362  SCIP_CALL( SCIPallocBufferArray(scip, &buf, SCIPexprGetNChildren(expr)) );
1363  for( int i = 0; i < SCIPexprGetNChildren(expr); ++i )
1364  {
1365  SCIP_CALL( eval(scip, SCIPexprGetChildren(expr)[i], exprintdata, x, buf[i]) );
1366  }
1367 
1368 #ifndef EVAL_USE_EXPRHDLR_ALWAYS
1369  if( SCIPisExprSum(scip, expr) )
1370  {
1371  val = SCIPgetConstantExprSum(expr);
1372  for( int i = 0; i < SCIPexprGetNChildren(expr); ++i )
1373  val += SCIPgetCoefsExprSum(expr)[i] * buf[i];
1374  }
1375  else if( SCIPisExprProduct(scip, expr) )
1376  {
1377  val = SCIPgetCoefExprProduct(expr);
1378  for( int i = 0; i < SCIPexprGetNChildren(expr); ++i )
1379  val *= buf[i];
1380  }
1381  else if( SCIPisExprPower(scip, expr) )
1382  {
1383  SCIP_Real exponent = SCIPgetExponentExprPow(expr);
1384  if( EPSISINT(exponent, 0.0) )
1385  evalIntPower(val, buf[0], (int)SCIPgetExponentExprPow(expr));
1386  else if( exponent == 0.5 )
1387  val = sqrt(buf[0]);
1388  else if( exponent < 1.0 )
1389  val = CppAD::pow(buf[0], SCIPgetExponentExprPow(expr));
1390  else
1391  {
1392  // workaround bug where CppAD claims pow(x,fractional>0) is nondiff at x=0
1393  // https://github.com/coin-or/CppAD/discussions/93#discussioncomment-327876
1394  AD<double> adzero(0.);
1395  val = CppAD::CondExpEq(buf[0], adzero, pow(buf[0]+std::numeric_limits<SCIP_Real>::epsilon(), exponent)-pow(std::numeric_limits<SCIP_Real>::epsilon(), exponent),
1396  pow(buf[0], exponent));
1397  }
1398  }
1399  else if( SCIPisExprSignpower(scip, expr) )
1400  {
1401  evalSignPower(val, buf[0], expr);
1402  }
1403  else if( SCIPisExprExp(scip, expr) )
1404  {
1405  val = exp(buf[0]);
1406  }
1407  else if( SCIPisExprLog(scip, expr) )
1408  {
1409  val = log(buf[0]);
1410  }
1411  else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "sin") == 0 )
1412  {
1413  val = sin(buf[0]);
1414  }
1415  else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") == 0 )
1416  {
1417  val = cos(buf[0]);
1418  }
1419  else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "erf") == 0 )
1420  {
1421  val = erf(buf[0]);
1422  }
1423  else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "abs") == 0 )
1424  {
1425  val = abs(buf[0]);
1426  }
1427  else if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "entropy") == 0 )
1428  {
1429  /* -x*log(x) if x > 0, else 0
1430  * https://coin-or.github.io/CppAD/doc/cond_exp.cpp.htm suggest to use 0 for the x=0 case
1431  * but then derivatives at 0 are 0, while they are actually infinite (see expr_entropy.c:bwdiff)
1432  * so we use -sqrt(x) for the x=0 case, as this also has value 0 and first derivative -inf at 0
1433  */
1434  val = CppAD::CondExpGt(buf[0], AD<double>(0.), -buf[0] * log(buf[0]), -sqrt(buf[0]));
1435  }
1436  else
1437 #endif
1438  {
1439  //SCIPerrorMessage("using derivative methods of exprhdlr %s in CppAD is not implemented yet", SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)));
1440  //CppAD::ErrorHandler::Call(true, __LINE__, __FILE__, "evalUser()", "Error: user expressions in CppAD not possible without CppAD user atomic facility");
1441  vector<Type> in(buf, buf + SCIPexprGetNChildren(expr));
1442  vector<Type> out(1);
1443 
1444  exprintdata->userexprs.push_back(new atomic_userexpr(scip, expr));
1445  (*exprintdata->userexprs.back())(in, out);
1446 
1447  val = out[0];
1448  }
1449 
1450  SCIPfreeBufferArray(scip, &buf);
1451 
1452  return SCIP_OKAY;
1453 }
1454 
1455 /** replacement for CppAD's default error handler
1456  *
1457  * In debug mode, CppAD gives an error when an evaluation contains a nan.
1458  * 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.
1459  * Since we cannot ignore this particular error, we ignore all.
1460  * @todo find a way to check whether the error corresponds to a nan and communicate this back
1461  */
1462 static
1464  bool known, /**< is the error from a known source? */
1465  int line, /**< line where error occured */
1466  const char* file, /**< file where error occured */
1467  const char* cond, /**< error condition */
1468  const char* msg /**< error message */
1469  )
1470 {
1471  SCIPdebugMessage("ignore CppAD error from %sknown source %s:%d: msg: %s exp: %s\n", known ? "" : "un", file, line, msg, cond);
1472 }
1473 
1474 /* install our error handler */
1475 static CppAD::ErrorHandler errorhandler(cppaderrorcallback);
1476 
1477 /** gets name and version of expression interpreter */
1478 const char* SCIPexprintGetName(void)
1479 {
1480  return CPPAD_PACKAGE_STRING;
1481 }
1482 
1483 /** gets descriptive text of expression interpreter */
1484 const char* SCIPexprintGetDesc(void)
1485 {
1486  return "Algorithmic Differentiation of C++ algorithms developed by B. Bell (github.com/coin-or/CppAD)";
1487 }
1488 
1489 /** gets capabilities of expression interpreter (using bitflags) */
1491  void
1492  )
1493 {
1495 }
1496 
1497 /** creates an expression interpreter object */
1499  SCIP* scip, /**< SCIP data structure */
1500  SCIP_EXPRINT** exprint /**< buffer to store pointer to expression interpreter */
1501  )
1502 {
1503  assert(exprint != NULL);
1504 
1505  *exprint = (SCIP_EXPRINT*)1u; /* some code checks that a non-NULL pointer is returned here, even though it may not point anywhere */
1506 
1507  return SCIP_OKAY;
1508 }
1509 
1510 /** frees an expression interpreter object */
1512  SCIP* scip, /**< SCIP data structure */
1513  SCIP_EXPRINT** exprint /**< expression interpreter that should be freed */
1514  )
1515 {
1516  assert(exprint != NULL);
1517 
1518  *exprint = NULL;
1519 
1520  return SCIP_OKAY;
1521 }
1522 
1523 /** compiles an expression and returns interpreter-specific data for expression
1524  *
1525  * can be called again with existing exprintdata if expression has been changed
1526  *
1527  * @attention *exprintdata needs to be initialized to NULL at first call
1528  * @attention the expression is assumed to use varidx expressions instead of var expressions
1529  */
1531  SCIP* scip, /**< SCIP data structure */
1532  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1533  SCIP_EXPR* rootexpr, /**< expression */
1534  SCIP_EXPRINTDATA** exprintdata /**< buffer to store pointer to compiled data */
1535  )
1536 {
1537  SCIP_EXPRITER* it;
1538  SCIP_EXPR* expr;
1539 
1540  assert(rootexpr != NULL);
1541  assert(exprintdata != NULL);
1542 
1543  if( *exprintdata == NULL )
1544  {
1545  *exprintdata = new SCIP_EXPRINTDATA();
1546  assert(*exprintdata != NULL);
1547  }
1548  else
1549  {
1550  (*exprintdata)->need_retape = true;
1551  (*exprintdata)->need_retape_always = false;
1552  (*exprintdata)->hesconstant = false;
1553  }
1554 
1555  SCIP_CALL( SCIPcreateExpriter(scip, &it) );
1557 
1558  std::set<int> varidxs;
1559  for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
1560  {
1561  /* cannot handle var-expressions in exprint so far, should be varidx expressions */
1562  assert(!SCIPisExprVar(scip, expr));
1563 
1564  if( SCIPisExprVaridx(scip, expr) )
1565  {
1566  varidxs.insert(SCIPgetIndexExprVaridx(expr));
1567  continue;
1568  }
1569 
1570  /* check whether expr is of a type we don't recognize */
1571 #ifndef EVAL_USE_EXPRHDLR_ALWAYS
1572  if( !SCIPisExprSum(scip, expr) &&
1573  !SCIPisExprProduct(scip, expr) &&
1574  !SCIPisExprPower(scip, expr) &&
1575  !SCIPisExprSignpower(scip, expr) &&
1576  !SCIPisExprExp(scip, expr) &&
1577  !SCIPisExprLog(scip, expr) &&
1578  strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "abs") != 0 &&
1579  strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "sin") != 0 &&
1580  strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "cos") != 0 &&
1581  strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "entropy") != 0 &&
1582  strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "erf") != 0 &&
1583  !SCIPisExprValue(scip, expr) )
1584 #endif
1585  {
1586  /* an expression for which we have no taping implemented in eval()
1587  * so we will try to use CppAD's atomic operand and the expression handler methods
1588  * however, only atomic_userexpr:forward() is implemented for p=0,1 at the moment, so we cannot do Hessians
1589  * also the exprhdlr needs to implement the fwdiff callback for derivatives
1590  */
1592  (*exprintdata)->userevalcapability &= SCIP_EXPRINTCAPABILITY_FUNCVALUE | SCIP_EXPRINTCAPABILITY_GRADIENT;
1593  else
1594  (*exprintdata)->userevalcapability &= SCIP_EXPRINTCAPABILITY_FUNCVALUE;
1595  }
1596 
1597  //if( strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), "xyz") == 0 )
1598  // (*exprintdata)->need_retape_always = true;
1599  }
1600 
1601  SCIPfreeExpriter(&it);
1602 
1603  (*exprintdata)->varidxs.reserve(varidxs.size());
1604  (*exprintdata)->varidxs.insert((*exprintdata)->varidxs.begin(), varidxs.begin(), varidxs.end());
1605 
1606  size_t n = (*exprintdata)->varidxs.size();
1607  (*exprintdata)->X.resize(n);
1608  (*exprintdata)->x.resize(n);
1609  (*exprintdata)->Y.resize(1);
1610 
1611  // check whether we are quadratic (or linear), so we can save on Hessian time (so
1612  // assumes simplified and skips over x^2 and x*y cases
1613  // not using SCIPcheckExprQuadratic(), because we don't need the quadratic form
1614  if( SCIPisExprSum(scip, rootexpr) )
1615  {
1616  (*exprintdata)->hesconstant = true;
1617  for( int i = 0; i < SCIPexprGetNChildren(rootexpr); ++i )
1618  {
1619  SCIP_EXPR* child;
1620  child = SCIPexprGetChildren(rootexpr)[i];
1621  /* linear term is ok */
1622  if( SCIPisExprVaridx(scip, child) )
1623  continue;
1624  /* square term is ok */
1625  if( SCIPisExprPower(scip, child) && SCIPgetExponentExprPow(child) == 2.0 && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[0]) )
1626  continue;
1627  /* bilinear term is ok */
1628  if( SCIPisExprProduct(scip, child) && SCIPexprGetNChildren(child) == 2 && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[0]) && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[1]) )
1629  continue;
1630  /* everything else means not quadratic (or not simplified) */
1631  (*exprintdata)->hesconstant = false;
1632  break;
1633  }
1634  SCIPdebugMsg(scip, "Hessian found %sconstant\n", (*exprintdata)->hesconstant ? "" : "not ");
1635  }
1636 
1637  return SCIP_OKAY;
1638 }
1639 
1640 /** frees interpreter data for expression */
1642  SCIP* scip, /**< SCIP data structure */
1643  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1644  SCIP_EXPR* expr, /**< expression */
1645  SCIP_EXPRINTDATA** exprintdata /**< pointer to pointer to compiled data to be freed */
1646  )
1647 {
1648  assert( exprintdata != NULL);
1649  assert(*exprintdata != NULL);
1650 
1651  SCIPfreeBlockMemoryArrayNull(scip, &(*exprintdata)->hesrowidxs, (*exprintdata)->hesnnz);
1652  SCIPfreeBlockMemoryArrayNull(scip, &(*exprintdata)->hescolidxs, (*exprintdata)->hesnnz);
1653 
1654  delete *exprintdata;
1655  *exprintdata = NULL;
1656 
1657  return SCIP_OKAY;
1658 }
1659 
1660 /** gives the capability to evaluate an expression by the expression interpreter
1661  *
1662  * In cases of user-given expressions, higher order derivatives may not be available for the user-expression,
1663  * even if the expression interpreter could handle these. This method allows to recognize that, e.g., the
1664  * Hessian for an expression is not available because it contains a user expression that does not provide
1665  * Hessians.
1666  */
1668  SCIP* scip, /**< SCIP data structure */
1669  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1670  SCIP_EXPR* expr, /**< expression */
1671  SCIP_EXPRINTDATA* exprintdata /**< interpreter-specific data for expression */
1672  )
1673 {
1674  assert(exprintdata != NULL);
1675 
1676  return exprintdata->userevalcapability;
1677 }/*lint !e715*/
1678 
1679 /** evaluates an expression */
1681  SCIP* scip, /**< SCIP data structure */
1682  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1683  SCIP_EXPR* expr, /**< expression */
1684  SCIP_EXPRINTDATA* exprintdata, /**< interpreter-specific data for expression */
1685  SCIP_Real* varvals, /**< values of variables */
1686  SCIP_Real* val /**< buffer to store value of expression */
1687  )
1688 {
1689  assert(expr != NULL);
1690  assert(exprintdata != NULL);
1691  assert(varvals != NULL);
1692  assert(val != NULL);
1693 
1694  size_t n = exprintdata->varidxs.size();
1695 
1696  if( n == 0 )
1697  {
1698  SCIP_CALL( SCIPevalExpr(scip, expr, NULL, 0L) );
1699  exprintdata->val = *val = SCIPexprGetEvalValue(expr);
1700  return SCIP_OKAY;
1701  }
1702 
1703  if( exprintdata->need_retape_always || exprintdata->need_retape )
1704  {
1705  /* free old Hessian data, if any */
1706  SCIPfreeBlockMemoryArrayNull(scip, &exprintdata->hesrowidxs, exprintdata->hesnnz);
1707  SCIPfreeBlockMemoryArrayNull(scip, &exprintdata->hescolidxs, exprintdata->hesnnz);
1708  exprintdata->hesvalues.clear();
1709  exprintdata->hesnnz = 0;
1710 
1711  for( size_t i = 0; i < n; ++i )
1712  {
1713  int idx = exprintdata->varidxs[i];
1714  exprintdata->X[i] = varvals[idx];
1715  exprintdata->x[i] = varvals[idx]; /* need this for a following grad or hessian eval with new_x = false */
1716  }
1717 
1718  // delete old atomic_userexprs before we start collecting new ones
1719  for( vector<atomic_userexpr*>::iterator it(exprintdata->userexprs.begin()); it != exprintdata->userexprs.end(); ++it )
1720  delete *it;
1721  exprintdata->userexprs.clear();
1722 
1723  CppAD::Independent(exprintdata->X);
1724 
1725  SCIP_CALL( eval(scip, expr, exprintdata, exprintdata->X, exprintdata->Y[0]) );
1726 
1727  exprintdata->f.Dependent(exprintdata->X, exprintdata->Y);
1728 
1729  exprintdata->val = Value(exprintdata->Y[0]);
1730  SCIPdebugMessage("Eval retaped and computed value %g\n", exprintdata->val);
1731 
1732  // the following is required if the gradient shall be computed by a reverse sweep later
1733  // exprintdata->val = exprintdata->f.Forward(0, exprintdata->x)[0];
1734 
1735  // https://coin-or.github.io/CppAD/doc/optimize.htm
1736  exprintdata->f.optimize();
1737 
1738  exprintdata->need_retape = false;
1739  }
1740  else
1741  {
1742  assert(exprintdata->x.size() >= n);
1743  for( size_t i = 0; i < n; ++i )
1744  exprintdata->x[i] = varvals[exprintdata->varidxs[i]];
1745 
1746  exprintdata->val = exprintdata->f.Forward(0, exprintdata->x)[0];
1747  SCIPdebugMessage("Eval used forward sweep to compute value %g\n", exprintdata->val);
1748  }
1749 
1750  *val = exprintdata->val;
1751 
1752  return SCIP_OKAY;
1753 }
1754 
1755 /** computes value and gradient of an expression */
1757  SCIP* scip, /**< SCIP data structure */
1758  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1759  SCIP_EXPR* expr, /**< expression */
1760  SCIP_EXPRINTDATA* exprintdata, /**< interpreter-specific data for expression */
1761  SCIP_Real* varvals, /**< values of variables, can be NULL if new_varvals is FALSE */
1762  SCIP_Bool new_varvals, /**< have variable values changed since last call to a point evaluation routine? */
1763  SCIP_Real* val, /**< buffer to store expression value */
1764  SCIP_Real* gradient /**< buffer to store expression gradient */
1765  )
1766 {
1767  assert(expr != NULL);
1768  assert(exprintdata != NULL);
1769  assert(varvals != NULL || new_varvals == FALSE);
1770  assert(val != NULL);
1771  assert(gradient != NULL);
1772 
1773  if( new_varvals )
1774  {
1775  SCIP_CALL( SCIPexprintEval(scip, exprint, expr, exprintdata, varvals, val) );
1776  }
1777  else
1778  *val = exprintdata->val;
1779 
1780  size_t n = exprintdata->varidxs.size();
1781 
1782  if( n == 0 )
1783  return SCIP_OKAY;
1784 
1785  vector<double> jac;
1786  if( exprintdata->userexprs.empty() )
1787  jac = exprintdata->f.Jacobian(exprintdata->x);
1788  else
1789  {
1790  // atomic_userexpr:reverse not implemented yet (even for p=0)
1791  // so force use of forward-eval for gradient (though that seems pretty inefficient)
1792  exprintdata->f.Forward(0, exprintdata->x);
1793  jac.resize(n);
1794  CppAD::JacobianFor(exprintdata->f, exprintdata->x, jac);
1795  }
1796 
1797  for( size_t i = 0; i < n; ++i )
1798  gradient[exprintdata->varidxs[i]] = jac[i];
1799 
1800 #ifdef SCIP_DEBUG
1801  SCIPinfoMessage(scip, NULL, "Grad for ");
1802  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
1803  SCIPinfoMessage(scip, NULL, "\nx = ");
1804  for( size_t i = 0; i < n; ++i )
1805  {
1806  SCIPinfoMessage(scip, NULL, "\t %g", exprintdata->x[i]);
1807  }
1808  SCIPinfoMessage(scip, NULL, "\ngrad = ");
1809  for( size_t i = 0; i < n; ++i )
1810  {
1811  SCIPinfoMessage(scip, NULL, "\t %g", jac[i]);
1812  }
1813  SCIPinfoMessage(scip, NULL, "\n");
1814 #endif
1815 
1816  return SCIP_OKAY;
1817 }
1818 
1819 /** gives sparsity pattern of lower-triangular part of Hessian
1820  *
1821  * Since the AD code might need to do a forward sweep, variable values need to be passed in here.
1822  *
1823  * Result will have `(*colidxs)[i] <= (*rowidixs)[i]` for `i=0..*nnz`.
1824  */
1826  SCIP* scip, /**< SCIP data structure */
1827  SCIP_EXPRINT* exprint, /**< interpreter data structure */
1828  SCIP_EXPR* expr, /**< expression */
1829  SCIP_EXPRINTDATA* exprintdata, /**< interpreter-specific data for expression */
1830  SCIP_Real* varvals, /**< values of variables */
1831  int** rowidxs, /**< buffer to return array with row indices of Hessian elements */
1832  int** colidxs, /**< buffer to return array with column indices of Hessian elements */
1833  int* nnz /**< buffer to return length of arrays */
1834  )
1835 {
1836  assert(expr != NULL);
1837  assert(exprintdata != NULL);
1838  assert(varvals != NULL);
1839  assert(rowidxs != NULL);
1840  assert(colidxs != NULL);
1841  assert(nnz != NULL);
1842 
1843  if( exprintdata->hesrowidxs == NULL )
1844  {
1845  assert(exprintdata->hescolidxs == NULL);
1846  assert(exprintdata->hesvalues.size() == 0);
1847  assert(exprintdata->hesnnz == 0);
1848 
1849  size_t n = exprintdata->varidxs.size();
1850  if( n == 0 )
1851  {
1852  *nnz = 0;
1853  return SCIP_OKAY;
1854  }
1855 
1856  size_t nn = n*n;
1857 
1858  if( exprintdata->need_retape_always )
1859  {
1860  // pretend dense
1861  exprintdata->hesnnz = (n * (n+1))/2;
1862  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &exprintdata->hesrowidxs, exprintdata->hesnnz) );
1863  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &exprintdata->hescolidxs, exprintdata->hesnnz) );
1864 
1865  int k = 0;
1866  for( size_t i = 0; i < n; ++i )
1867  for( size_t j = 0; j <= i; ++j )
1868  {
1869  exprintdata->hesrowidxs[k] = exprintdata->varidxs[i];
1870  exprintdata->hescolidxs[k] = exprintdata->varidxs[j];
1871  k++;
1872  }
1873 
1874 #ifdef SCIP_DEBUG
1875  SCIPinfoMessage(scip, NULL, "HessianSparsity for ");
1876  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
1877  SCIPinfoMessage(scip, NULL, ": all elements, due to discontinuities\n");
1878 #endif
1879 
1880  *rowidxs = exprintdata->hesrowidxs;
1881  *colidxs = exprintdata->hescolidxs;
1882  *nnz = exprintdata->hesnnz;
1883 
1884  return SCIP_OKAY;
1885  }
1886 
1887  if( exprintdata->need_retape )
1888  {
1889  SCIP_Real val;
1890  SCIP_CALL( SCIPexprintEval(scip, exprint, expr, exprintdata, varvals, &val) );
1891  } /*lint !e438*/
1892 
1893  // following https://github.com/coin-or/CppAD/blob/20180000.0/cppad_ipopt/src/vec_fun_pattern.cpp
1894 
1895  SCIPdebugMessage("calling ForSparseJac\n");
1896 
1897  vector<bool> r(nn, false);
1898  for( size_t i = 0; i < n; ++i )
1899  r[i*n+i] = true;
1900  (void) exprintdata->f.ForSparseJac(n, r); // need to compute sparsity for Jacobian first
1901 
1902  SCIPdebugMessage("calling RevSparseHes\n");
1903 
1904  // this was originally
1905  // vector<bool> s(1, true);
1906  // hessparsity = exprintdata->f.RevSparseHes(n, s);
1907  // RevSparseHes is just calling RevSparseHesCase
1908  // to avoid copying hessparsity, call RevSparseHesCase directly
1909  vector<bool> hessparsity;
1910  exprintdata->f.RevSparseHesCase(true, false, n, vector<bool>(1, true), hessparsity);
1911 
1912  // count number of hessian elements and setup hessparsity_pattern
1913  exprintdata->hessparsity_pattern.resize(0, 0); // clear old data, if any
1914  exprintdata->hessparsity_pattern.resize(n, n);
1915  size_t hesnnz_full = 0; // number of nonzeros in full matrix, that is, not only lower-diagonal
1916  for( size_t i = 0; i < nn; ++i )
1917  if( hessparsity[i] )
1918  {
1919  size_t row = i / n;
1920  size_t col = i % n;
1921 
1922  ++hesnnz_full;
1923 
1924  if( col > row )
1925  continue;
1926  ++exprintdata->hesnnz;
1927  }
1928 
1929  // hessian sparsity in sparse form
1930  // - hesrowidxs,hescolidxs are nonzero entries in the lower-diagonal of the Hessian and are returned to the caller; indices are variable indices as in expression
1931  // - hessparsity_row,hessparsity_col are nonzero entries in the full Hessian and are used in SCIPexprintHessian(); indices are 0..n-1 (n=dimension of f)
1932  // - hessparsity_pattern is the same information as hessparsity_row,hessparsity_col, but stored in different form
1933  // originally we were calling CppAD::local::sparsity_user2internal(exprintdata->hessparsity_pattern, exprintdata->hessparsity, n, n, false, ""),
1934  // which was again looping over hessparsity, so this is included into one loop here
1935  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &exprintdata->hesrowidxs, exprintdata->hesnnz) );
1936  SCIP_CALL( SCIPallocBlockMemoryArray(scip, &exprintdata->hescolidxs, exprintdata->hesnnz) );
1937 
1938  exprintdata->hessparsity_row.resize(hesnnz_full);
1939  exprintdata->hessparsity_col.resize(hesnnz_full);
1940 
1941  for( size_t i = 0, j = 0, k = 0; i < nn; ++i )
1942  if( hessparsity[i] )
1943  {
1944  size_t row = i / n;
1945  size_t col = i % n;
1946 
1947  if( (size_t)exprintdata->hesnnz <= nn/4 )
1948  {
1949  // we need hessparsity_pattern only if SCIPexprintHessian is doing a sparse hessian eval; nn/4 is the threshold for that
1950  exprintdata->hessparsity_pattern.add_element(row, col);
1951  }
1952 
1953  if( col > row )
1954  {
1955  // add above-diagonal elements into end part of hessparsity_row/col, so we have a 1:1 correspondence
1956  // between hessparsity_row/col and hesrowidxs/hescolidxs
1957  assert(exprintdata->hesnnz + k < hesnnz_full);
1958  exprintdata->hessparsity_row[exprintdata->hesnnz + k] = row;
1959  exprintdata->hessparsity_col[exprintdata->hesnnz + k] = col;
1960  ++k;
1961  continue;
1962  }
1963 
1964  exprintdata->hessparsity_row[j] = row;
1965  exprintdata->hessparsity_col[j] = col;
1966 
1967  assert(j < (size_t)exprintdata->hesnnz);
1968  exprintdata->hesrowidxs[j] = exprintdata->varidxs[row];
1969  exprintdata->hescolidxs[j] = exprintdata->varidxs[col];
1970  ++j;
1971  }
1972 
1973 #ifdef SCIP_DEBUG
1974  SCIPinfoMessage(scip, NULL, "HessianSparsity for ");
1975  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
1976  SCIPinfoMessage(scip, NULL, ":");
1977  for( int i = 0; i < exprintdata->hesnnz; ++i )
1978  {
1979  SCIPinfoMessage(scip, NULL, " (%d,%d)", exprintdata->hesrowidxs[i], exprintdata->hescolidxs[i]);
1980  }
1981  SCIPinfoMessage(scip, NULL, "\n");
1982 #endif
1983  }
1984 
1985  *rowidxs = exprintdata->hesrowidxs;
1986  *colidxs = exprintdata->hescolidxs;
1987  *nnz = exprintdata->hesnnz;
1988 
1989  return SCIP_OKAY;
1990 }
1991 
1992 /** computes value and Hessian of an expression
1993  *
1994  * Returned arrays `rowidxs` and `colidxs` and number of elements `nnz` are the same as given by SCIPexprintHessianSparsity().
1995  * Returned array `hessianvals` will contain the corresponding Hessian elements.
1996  */
1998  SCIP* scip, /**< SCIP data structure */
1999  SCIP_EXPRINT* exprint, /**< interpreter data structure */
2000  SCIP_EXPR* expr, /**< expression */
2001  SCIP_EXPRINTDATA* exprintdata, /**< interpreter-specific data for expression */
2002  SCIP_Real* varvals, /**< values of variables, can be NULL if new_varvals is FALSE */
2003  SCIP_Bool new_varvals, /**< have variable values changed since last call to an evaluation routine? */
2004  SCIP_Real* val, /**< buffer to store function value */
2005  int** rowidxs, /**< buffer to return array with row indices of Hessian elements */
2006  int** colidxs, /**< buffer to return array with column indices of Hessian elements */
2007  SCIP_Real** hessianvals, /**< buffer to return array with Hessian elements */
2008  int* nnz /**< buffer to return length of arrays */
2009  )
2010 {
2011  assert(expr != NULL);
2012  assert(exprintdata != NULL);
2013 
2014  if( exprintdata->hesrowidxs == NULL )
2015  {
2016  /* setup sparsity if not done yet */
2017  int dummy1;
2018  int* dummy2;
2019 
2020  assert(exprintdata->hescolidxs == NULL);
2021  assert(exprintdata->hesvalues.size() == 0);
2022  assert(exprintdata->hesnnz == 0);
2023 
2024  SCIP_CALL( SCIPexprintHessianSparsity(scip, exprint, expr, exprintdata, varvals, &dummy2, &dummy2, &dummy1) );
2025 
2026  new_varvals = FALSE;
2027  }
2028 
2029  if( new_varvals )
2030  {
2031  SCIP_CALL( SCIPexprintEval(scip, exprint, expr, exprintdata, varvals, val) );
2032  }
2033  else
2034  *val = exprintdata->val;
2035 
2036  size_t n = exprintdata->varidxs.size();
2037 
2038  // eval hessian; if constant, then only if not evaluated yet (hesvalues has size 0, but hesnnz > 0)
2039  if( n > 0 && (!exprintdata->hesconstant || exprintdata->hesvalues.size() < (size_t)exprintdata->hesnnz) )
2040  {
2041  size_t nn = n*n;
2042 
2043  if( exprintdata->hesvalues.size() == 0 )
2044  {
2045  // using here hessparsity_row.size() instead of hesnnz as we want to pass hesvalues to CppAD and it prefers to calculate a full Hessian
2046  exprintdata->hesvalues.resize(exprintdata->hessparsity_row.size());
2047  }
2048 
2049  // these use reverse mode
2050  if( (size_t)exprintdata->hesnnz > nn/4 )
2051  {
2052  // use dense Hessian if there are many nonzero elements (approx more than half (recall only lower (or upper)-triangular is considered))
2053  // because it seems faster than the sparse Hessian if there isn't actually much sparsity
2054  vector<double> hess = exprintdata->f.Hessian(exprintdata->x, 0);
2055  for( int i = 0; i < exprintdata->hesnnz; ++i )
2056  exprintdata->hesvalues[i] = hess[exprintdata->hessparsity_row[i] * n + exprintdata->hessparsity_col[i]];
2057  }
2058  else
2059  {
2060  // originally, this was hess = exprintdata->f.SparseHessian(exprintdata->x, vector<double>(1, 1.0), exprintdata->hessparsity),
2061  // where hess was a dense nxn matrix as in the case above
2062  // to reuse the coloring of the sparsity pattern and use also a sparse matrix for the Hessian values, we now call SparseHessianCompute directly
2063  exprintdata->f.SparseHessianCompute(exprintdata->x, vector<double>(1, 1.0), exprintdata->hessparsity_pattern, exprintdata->hessparsity_row, exprintdata->hessparsity_col, exprintdata->hesvalues, exprintdata->heswork);
2064 
2065 #ifndef NDEBUG
2066  // hessparsity_row, hessparsity_col, hessparsity_pattern are no longer used by SparseHessianCompute after coloring has been computed in first call, except for some asserts, so we free some mem
2067  exprintdata->hessparsity_row.clear();
2068  exprintdata->hessparsity_col.clear();
2069  exprintdata->hessparsity_pattern.resize(0,0);
2070 #endif
2071  }
2072  }
2073 
2074 #ifdef SCIP_DEBUG
2075  SCIPinfoMessage(scip, NULL, "Hessian for ");
2076  SCIP_CALL( SCIPprintExpr(scip, expr, NULL) );
2077  SCIPinfoMessage(scip, NULL, "\nat x = ");
2078  for( size_t i = 0; i < n; ++i )
2079  {
2080  SCIPinfoMessage(scip, NULL, "\t %g", exprintdata->x[i]);
2081  }
2082  SCIPinfoMessage(scip, NULL, "\nis ");
2083  for( int i = 0; i < exprintdata->hesnnz; ++i )
2084  {
2085  SCIPinfoMessage(scip, NULL, " (%d,%d)=%g", exprintdata->hesrowidxs[i], exprintdata->hescolidxs[i], exprintdata->hesvalues[i]);
2086  }
2087  SCIPinfoMessage(scip, NULL, "\n");
2088 #endif
2089 
2090  *rowidxs = exprintdata->hesrowidxs;
2091  *colidxs = exprintdata->hescolidxs;
2092  *hessianvals = exprintdata->hesvalues.data();
2093  *nnz = exprintdata->hesnnz;
2094 
2095  return SCIP_OKAY;
2096 }
static SCIP_RETCODE eval(SCIP *scip, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, const vector< Type > &x, Type &val)
#define NULL
Definition: def.h:267
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:501
const char * SCIPexprintGetName(void)
#define SCIPallocBlockMemoryArray(scip, ptr, num)
Definition: scip_mem.h:93
methods to interpret (evaluate) an expression "fast"
SCIP_Bool SCIPexprhdlrHasFwdiff(SCIP_EXPRHDLR *exprhdlr)
Definition: expr.c:605
SCIP_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition: scip_expr.c:1486
SCIP_RETCODE SCIPexprintCompile(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *rootexpr, SCIP_EXPRINTDATA **exprintdata)
int SCIPexprGetNChildren(SCIP_EXPR *expr)
Definition: expr.c:3854
#define infinity
Definition: gastrans.c:80
const char * SCIPexprhdlrGetName(SCIP_EXPRHDLR *exprhdlr)
Definition: expr.c:545
void posintpower(const vector< Type > &in, vector< Type > &out, size_t exponent)
SCIP_RETCODE SCIPexprintFreeData(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA **exprintdata)
static void evalSignPower(CppAD::AD< Type > &resultant, const CppAD::AD< Type > &arg, SCIP_EXPR *expr)
SCIP_Real SCIPgetConstantExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1567
SCIP_RETCODE SCIPexprintHessianSparsity(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, int **rowidxs, int **colidxs, int *nnz)
#define FALSE
Definition: def.h:94
#define EPSISINT(x, eps)
Definition: def.h:210
SCIP_Real SCIPgetExponentExprPow(SCIP_EXPR *expr)
Definition: expr_pow.c:3457
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
const char * SCIPexprintGetDesc(void)
SCIP_Bool SCIPisExprLog(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_log.c:648
#define CPPAD_MAX_NUM_THREADS
#define SCIPdebugMessage
Definition: pub_message.h:96
SCIP_RETCODE SCIPexprintGrad(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, SCIP_Real *gradient)
#define SCIPfreeBufferArray(scip, ptr)
Definition: scip_mem.h:136
unsigned int SCIP_EXPRINTCAPABILITY
#define SCIPdebugMsg
Definition: scip_message.h:78
SCIP_EXPR * SCIPexpriterGetCurrent(SCIP_EXPRITER *iterator)
Definition: expriter.c:683
SCIP_VAR ** x
Definition: circlepacking.c:63
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
#define SIGN(x)
Definition: expr_pow.c:72
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1464
SCIP_RETCODE SCIPevalExpr(SCIP *scip, SCIP_EXPR *expr, SCIP_SOL *sol, SCIP_Longint soltag)
Definition: scip_expr.c:1635
public functions to work with algebraic expressions
SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr)
Definition: expr.c:3864
SCIP_Real * SCIPgetCoefsExprSum(SCIP_EXPR *expr)
Definition: expr_sum.c:1552
SCIP_RETCODE SCIPcallExprEvalFwdiff(SCIP *scip, SCIP_EXPR *expr, SCIP_Real *childrenvalues, SCIP_Real *direction, SCIP_Real *val, SCIP_Real *dot)
Definition: scip_expr.c:2212
SCIP_Real SCIPexprGetEvalValue(SCIP_EXPR *expr)
Definition: expr.c:3928
handler for variable index expressions
interval arithmetics for provable bounds
static void cppaderrorcallback(bool known, int line, const char *file, const char *cond, const char *msg)
public functions to work with algebraic expressions
SCIP_Bool SCIPisExprValue(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1442
power and signed power expression handlers
#define REALABS(x)
Definition: def.h:197
SCIP_RETCODE SCIPexprintCreate(SCIP *scip, SCIP_EXPRINT **exprint)
SCIP_Bool SCIPisExprSum(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1453
#define SCIP_CALL(x)
Definition: def.h:380
SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void)
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2337
SCIP_Real SCIPgetCoefExprProduct(SCIP_EXPR *expr)
SCIP_RETCODE SCIPexprintEval(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, SCIP_Real *val)
SCIP_RETCODE SCIPcallExprEval(SCIP *scip, SCIP_EXPR *expr, SCIP_Real *childrenvalues, SCIP_Real *val)
Definition: scip_expr.c:2185
#define SCIPallocBufferArray(scip, ptr, num)
Definition: scip_mem.h:124
#define SCIP_EXPRINTCAPABILITY_ALL
#define SCIP_EXPRINTCAPABILITY_HESSIAN
SCIP_EXPRINTCAPABILITY SCIPexprintGetExprCapability(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata)
logarithm expression handler
#define SCIP_Bool
Definition: def.h:91
unsigned short Type
Definition: cons_xor.c:132
SCIP_RETCODE SCIPexprintFree(SCIP *scip, SCIP_EXPRINT **exprint)
struct SCIP_ExprInt SCIP_EXPRINT
SCIP_EXPR * SCIPexpriterGetNext(SCIP_EXPRITER *iterator)
Definition: expriter.c:858
SCIP_EXPRHDLR * SCIPexprGetHdlr(SCIP_EXPR *expr)
Definition: expr.c:3877
SCIP_Bool SCIPisExprPower(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1475
SCIP_RETCODE SCIPexprintHessian(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, int **rowidxs, int **colidxs, SCIP_Real **hessianvals, int *nnz)
static CppAD::ErrorHandler errorhandler(cppaderrorcallback)
static void evalIntPower(Type &resultant, const Type &arg, const int exponent)
void SCIPfreeExpriter(SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2351
SCIP_Real * r
Definition: circlepacking.c:59
#define SCIP_EXPRINTCAPABILITY_FUNCVALUE
SCIP_Bool SCIPisExprSignpower(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_pow.c:3243
SCIP_Bool SCIPisExprVaridx(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_varidx.c:252
SCIP_Bool SCIPisExprVar(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1431
SCIP_Bool SCIPisExprExp(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_exp.c:528
#define SCIP_Real
Definition: def.h:173
#define SCIP_INVALID
Definition: def.h:193
SCIP_Real SCIPgetValueExprValue(SCIP_EXPR *expr)
Definition: expr_value.c:294
atomic_userexpr(SCIP *scip_, SCIP_EXPR *expr_)
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:111
common defines and data types used in all packages of SCIP
SCIP_Bool SCIPexpriterIsEnd(SCIP_EXPRITER *iterator)
Definition: expriter.c:969
#define SCIPABORT()
Definition: def.h:352
int SCIPgetIndexExprVaridx(SCIP_EXPR *expr)
Definition: expr_varidx.c:267
struct SCIP_ExprIntData SCIP_EXPRINTDATA
exponential expression handler
#define SCIP_EXPRINTCAPABILITY_GRADIENT