exprinterpret_cppad.cpp
Go to the documentation of this file.
22 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
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
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
46 /* defining NO_CPPAD_USER_ATOMIC disables the use of our own implementation of derivaties of power operators
48 * our customized implementation should give better results (tighter intervals) for the interval data type
59 * we need to include the intervalarithext.h very early and require the interval operations to be in the CppAD namespace */
70 * 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.
81 * -Wshadow was too strict with some versions of GCC 4 (https://stackoverflow.com/questions/2958457/gcc-wshadow-is-too-strict)
93 * To run it in a multithreading environment, a special CppAD memory allocator that is aware of the multiple threads has to be used.
95 * To implement this, we follow the team_pthread example of CppAD, which uses pthread's thread-specific data management.
121 /* if no thread_number for this thread yet, then assign a new thread number to the current thread
154 * The purpose is to make sure that init_parallel() is called before any multithreading is started.
287 * @return [0,0] if first argument is [0,0] independent of whether the second argument is an empty interval or not
322 : val(0.0), need_retape(true), int_need_retape(true), need_retape_always(false), userevalcapability(SCIP_EXPRINTCAPABILITY_ALL), blkmem(NULL), root(NULL)
346 SCIP_EXPRINTCAPABILITY userevalcapability; /**< (intersection of) capabilities of evaluation rountines of user expressions */
356 * 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.
376 * 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).
401 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
433 * This class implements forward and reverse operations for the function x -> x^p for use within CppAD.
434 * 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.
436 * @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)
457 * TODO according to the CppAD 2018 docu, using this function is deprecated; what is the modern way to do this?
479 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
480 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
514 // ty[2] = 1/2 * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1] * tx[1] + exponent * pow(tx[0], exponent-1) * tx[2];
539 * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
540 * 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.
543 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
587 // px[0] = py[0] * exponent * pow(tx[0], exponent-1) + py[1] * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1];
608 * 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.
624 * 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).
644 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
690 * This class implements forward and reverse operations for the function x -> sign(x)abs(x)^p for use within CppAD.
691 * While we otherwise would have to use discontinuous sign and abs functions, our own implementation allows to provide
694 * @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)
715 * TODO according to the CppAD 2018 docu, using this function is deprecated; what is the modern way to do this?
737 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
738 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
796 * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$
797 * 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.
800 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
846 // 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]
878 * 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.
894 * 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).
914 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
948 * TODO according to the CppAD 2018 docu, using this function is deprecated; what is the modern way to do this?
962 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
963 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
1025 CppAD::vector<SCIPInterval>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
1045 // 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]
1077 * 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.
1093 * 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).
1109 * 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.
1113 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
1216 * This class implements forward and reverse operations for a function given by a user expression for use within CppAD.
1237 * TODO according to the CppAD 2018 docu, using this function is deprecated; what is the modern way to do this?
1259 * 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
1263 * 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)
1264 * 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])
1269 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
1270 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */
1371 * 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)),
1411 CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
1474 * 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.
1503 * 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).
1533 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
1956 /* for each argument, we collect it's linear index from lincoefs, it's square coefficients and all factors from bilinear terms
1963 /* there are no quadratic terms with argidx in its first argument, that should be easy to handle */
2081 * This may be the case if the evaluation sequence depends on values of operands (e.g., in case of abs, sign, signpower, ...).
2118 * 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.
2131 SCIPdebugMessage("ignore CppAD error from %sknown source %s:%d: msg: %s exp: %s\n", known ? "" : "un", file, line, msg, cond);
2146 return "Algorithmic Differentiation of C++ algorithms developed by B. Bell (www.coin-or.org/CppAD)";
2240 * In cases of user-given expressions, higher order derivatives may not be available for the user-expression,
2241 * even if the expression interpreter could handle these. This method allows to recognize that, e.g., the
2242 * Hessian for an expression is not available because it contains a user expression that does not provide
2436 SCIP_Bool new_varvals, /**< have variable values changed since last call to a point evaluation routine? */
2438 SCIP_Real* gradient /**< buffer to store expression gradient, need to have length at least SCIPexprtreeGetNVars(tree) */
2470 SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n");
2471 SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t %g", gradient[i]); printf("\n");
2483 SCIP_INTERVAL* varvals, /**< interval values of variables, can be NULL if new_varvals is FALSE */
2484 SCIP_Bool new_varvals, /**< have variable interval values changed since last call to an interval evaluation routine? */
2486 SCIP_INTERVAL* gradient /**< buffer to store expression interval gradient, need to have length at least SCIPexprtreeGetNVars(tree) */
2516 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");
2517 SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t [%g,%g]", SCIPintervalGetInf(gradient[i]), SCIPintervalGetSup(gradient[i])); printf("\n");
2533 SCIP_Bool* sparsity /**< buffer to store sparsity pattern of Hessian, sparsity[i+n*j] indicates whether entry (i,j) is nonzero in the hessian */
2559 SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2590 SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2591 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");
2606 SCIP_Bool new_varvals, /**< have variable values changed since last call to an evaluation routine? */
2653 SCIPdebugMessage("HessianDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n");
2654 SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n");
2655 SCIPdebugMessage("hess ="); for (int i = 0; i < n*n; ++i) printf("\t %g", hessian[i]); printf("\n");
Definition: exprinterpret_cppad.cpp:439
static void evalMax(Type &resultant, const Type &arg1, const Type &arg2)
Definition: exprinterpret_cppad.cpp:1653
Definition: type_expr.h:57
bool IdenticalEqualPar(const SCIPInterval &x, const SCIPInterval &y)
Definition: exprinterpret_cppad.cpp:207
SCIP_RETCODE SCIPexprintNewParametrization(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree)
Definition: exprinterpret_cppad.cpp:2279
static SCIP_RETCODE eval(SCIP_EXPR *expr, const vector< Type > &x, SCIP_Real *param, Type &val)
Definition: exprinterpret_cppad.cpp:1769
SCIP_RETCODE SCIPexprtreeEvalInt(SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_INTERVAL *val)
Definition: expr.c:8741
Definition: intervalarith.h:37
Definition: type_expr.h:59
methods to interpret (evaluate) an expression tree "fast"
int * SCIPexprGetMonomialChildIndices(SCIP_EXPRDATA_MONOMIAL *monomial)
Definition: expr.c:5924
Definition: type_expr.h:72
static void evalUser(Type &resultant, const Type *args, SCIP_EXPR *expr)
Definition: exprinterpret_cppad.cpp:1581
Definition: type_expr.h:65
static void evalSignPower(Type &resultant, const Type &arg, SCIP_EXPR *expr)
Definition: exprinterpret_cppad.cpp:1129
Definition: type_expr.h:74
SCIP_Real * SCIPexprtreeGetParamVals(SCIP_EXPRTREE *tree)
Definition: expr.c:8634
#define SCIP_EXPRINTCAPABILITY_INTGRADIENT
Definition: type_exprinterpret.h:39
SCIP_EXPRINTCAPABILITY SCIPexprintGetExprtreeCapability(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree)
Definition: exprinterpret_cppad.cpp:2245
SCIP_Real SCIPexprGetRealPowerExponent(SCIP_EXPR *expr)
Definition: expr.c:5760
SCIP_Real SCIPexprGetPolynomialConstant(SCIP_EXPR *expr)
Definition: expr.c:5892
Definition: type_expr.h:62
Definition: type_expr.h:73
SCIP_RETCODE SCIPexprintGrad(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, SCIP_Real *gradient)
Definition: exprinterpret_cppad.cpp:2432
static void evalAbs(Type &resultant, const Type &arg)
Definition: exprinterpret_cppad.cpp:1693
SCIP_RETCODE SCIPexprintEvalInt(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_INTERVAL *val)
Definition: exprinterpret_cppad.cpp:2366
Definition: type_expr.h:100
Definition: type_expr.h:40
static void analyzeTree(SCIP_EXPRINTDATA *data, SCIP_EXPR *expr)
Definition: exprinterpret_cppad.cpp:2084
Definition: exprinterpret_cppad.cpp:63
Definition: type_expr.h:55
static void evalMin(Type &resultant, const Type &arg1, const Type &arg2)
Definition: exprinterpret_cppad.cpp:1623
Definition: type_expr.h:46
SCIP_EXPRDATA_MONOMIAL ** SCIPexprGetMonomials(SCIP_EXPR *expr)
Definition: expr.c:5868
public methods for expressions, expression trees, expression graphs, and related stuff ...
int SCIPexprGetMonomialNFactors(SCIP_EXPRDATA_MONOMIAL *monomial)
Definition: expr.c:5914
#define SCIP_EXPRINTCAPABILITY_INTFUNCVALUE
Definition: type_exprinterpret.h:37
SCIPInterval CondExpOp(enum CppAD::CompareOp cop, const SCIPInterval &left, const SCIPInterval &right, const SCIPInterval &trueCase, const SCIPInterval &falseCase)
Definition: exprinterpret_cppad.cpp:163
Definition: type_expr.h:47
SCIPInterval signpow(const SCIPInterval &x, const SCIP_Real p)
Definition: intervalarithext.h:319
Definition: type_expr.h:49
Definition: type_expr.h:76
Definition: type_expr.h:53
static bool univariate_for_sparse_jac(size_t q, const CppAD::vector< bool > &r, CppAD::vector< bool > &s)
Definition: exprinterpret_cppad.cpp:360
Definition: type_expr.h:51
SCIP_Real * SCIPexprGetQuadLinearCoefs(SCIP_EXPR *expr)
Definition: expr.c:5844
static void cppaderrorcallback(bool known, int line, const char *file, const char *cond, const char *msg)
Definition: exprinterpret_cppad.cpp:2123
SCIP_RETCODE SCIPexprCopyDeep(BMS_BLKMEM *blkmem, SCIP_EXPR **targetexpr, SCIP_EXPR *sourceexpr)
Definition: expr.c:6145
Definition: type_expr.h:52
Definition: exprinterpret_cppad.cpp:1219
SCIP_RETCODE SCIPexprintFree(SCIP_EXPRINT **exprint)
Definition: exprinterpret_cppad.cpp:2177
Definition: type_retcode.h:33
SCIP_RETCODE SCIPexprEvalUser(SCIP_EXPR *expr, SCIP_Real *argvals, SCIP_Real *val, SCIP_Real *gradient, SCIP_Real *hessian)
Definition: expr.c:7973
SCIP_RETCODE exprEvalUser(SCIP_EXPR *expr, Type *x, Type &funcval, Type *gradient, Type *hessian)
Definition: exprinterpret_cppad.cpp:1191
SCIP_RETCODE SCIPexprintCreate(BMS_BLKMEM *blkmem, SCIP_EXPRINT **exprint)
Definition: exprinterpret_cppad.cpp:2160
Definition: type_retcode.h:34
static void evalSqrt(Type &resultant, const Type &arg)
Definition: exprinterpret_cppad.cpp:1682
Definition: struct_expr.h:46
SCIP_EXPRINTCAPABILITY SCIPexprGetUserEvalCapability(SCIP_EXPR *expr)
Definition: expr.c:5966
#define SCIP_EXPRINTCAPABILITY_HESSIAN
Definition: type_exprinterpret.h:40
SCIP_Real SCIPexprGetSignPowerExponent(SCIP_EXPR *expr)
Definition: expr.c:5782
SCIP_RETCODE SCIPexprintEval(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Real *val)
Definition: exprinterpret_cppad.cpp:2298
SCIPInterval azmul(const SCIPInterval &x, const SCIPInterval &y)
Definition: exprinterpret_cppad.cpp:290
Definition: type_expr.h:39
SCIP_RETCODE SCIPexprintFreeData(SCIP_EXPRINTDATA **interpreterdata)
Definition: exprinterpret_cppad.cpp:2259
SCIP_EXPRINTDATA * SCIPexprtreeGetInterpreterData(SCIP_EXPRTREE *tree)
Definition: expr.c:8659
SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void)
Definition: exprinterpret_cppad.cpp:2150
static bool univariate_rev_sparse_jac(size_t q, const CppAD::vector< bool > &r, CppAD::vector< bool > &s)
Definition: exprinterpret_cppad.cpp:380
Definition: type_expr.h:58
Definition: type_expr.h:64
Definition: struct_expr.h:89
SCIP_RETCODE SCIPexprintHessianDense(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool new_varvals, SCIP_Real *val, SCIP_Real *hessian)
Definition: exprinterpret_cppad.cpp:2602
SCIP_Real SCIPexprGetMonomialCoef(SCIP_EXPRDATA_MONOMIAL *monomial)
Definition: expr.c:5904
static CppAD::ErrorHandler errorhandler(cppaderrorcallback)
static void evalIntPower(Type &resultant, const Type &arg, const int exponent)
Definition: exprinterpret_cppad.cpp:1722
Definition: type_expr.h:63
Definition: exprinterpret_cppad.cpp:697
SCIP_RETCODE SCIPexprEvalIntUser(SCIP_EXPR *expr, SCIP_Real infinity, SCIP_INTERVAL *argvals, SCIP_INTERVAL *val, SCIP_INTERVAL *gradient, SCIP_INTERVAL *hessian)
Definition: expr.c:7996
Definition: type_expr.h:79
#define SCIP_EXPRINTCAPABILITY_FUNCVALUE
Definition: type_exprinterpret.h:36
SCIP_QUADELEM * SCIPexprGetQuadElements(SCIP_EXPR *expr)
Definition: expr.c:5819
std::ostream & operator<<(std::ostream &out, const SCIP_INTERVAL &x)
Definition: exprinterpret_cppad.cpp:302
Definition: type_expr.h:71
Definition: type_expr.h:54
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)
Definition: exprinterpret_cppad.cpp:2479
SCIP_RETCODE SCIPexprintHessianSparsityDense(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool *sparsity)
Definition: exprinterpret_cppad.cpp:2529
SCIP_EXPORT char SCIPexprintCppADInitParallel(void)
Definition: exprinterpret_cppad.cpp:142
SCIP_Real * SCIPexprGetMonomialExponents(SCIP_EXPRDATA_MONOMIAL *monomial)
Definition: expr.c:5934
Definition: type_expr.h:75
Definition: type_expr.h:56
static bool univariate_rev_sparse_hes(const CppAD::vector< bool > &vx, const CppAD::vector< bool > &s, CppAD::vector< bool > &t, size_t q, const CppAD::vector< bool > &r, const CppAD::vector< bool > &u, CppAD::vector< bool > &v)
Definition: exprinterpret_cppad.cpp:400
Definition: struct_expr.h:55
Definition: type_retcode.h:35
common defines and data types used in all packages of SCIP
Definition: type_expr.h:50
SCIP_Real SCIPexprGetLinearConstant(SCIP_EXPR *expr)
Definition: expr.c:5806
void SCIPexprtreeSetInterpreterData(SCIP_EXPRTREE *tree, SCIP_EXPRINTDATA *interpreterdata)
Definition: expr.c:8669
bool GreaterThanOrZero(const SCIPInterval &x)
Definition: exprinterpret_cppad.cpp:231
Definition: exprinterpret_cppad.cpp:311
#define SCIP_EXPRINTCAPABILITY_GRADIENT
Definition: type_exprinterpret.h:38
SCIP_RETCODE SCIPexprintCompile(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree)
Definition: exprinterpret_cppad.cpp:2190
SCIP_RETCODE SCIPexprtreeEval(SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Real *val)
Definition: expr.c:8725
static void posintpower(const vector< Type > &in, vector< Type > &out, size_t exponent)
Definition: exprinterpret_cppad.cpp:660
C++ extensions to interval arithmetics for provable bounds.
Definition: type_expr.h:48
Definition: type_expr.h:41
memory allocation routines