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 */ 67 * It is wise to set it to a power of 2, so that if the tape id overflows, it is likely to start at 0 again, which avoids difficult to debug errors. 81 * To run it in a multithreading environment, a special CppAD memory allocator that is aware of the multiple threads has to be used. 83 * To implement this, we follow the team_pthread example of CppAD, which uses pthread's thread-specific data management. 88 /* workaround error message regarding missing implementation of tanh during initialization of static variables (see cppad/local/erf.hpp) */ 136 * we store the thread number incremented by one, to distinguish the absence of data (=0) from existing data 142 SCIPdebugMessage("Assigning thread number %lu to thread %p.\n", (long unsigned int)ncurthreads, (void*)pthread_self()); 183 * The purpose is to make sure that init_parallel() is called before any multithreading is started. 345 : val(0.0), need_retape(true), int_need_retape(true), need_retape_always(false), userevalcapability(SCIP_EXPRINTCAPABILITY_ALL), blkmem(NULL), root(NULL) 369 SCIP_EXPRINTCAPABILITY userevalcapability; /**< (intersection of) capabilities of evaluation rountines of user expressions */ 379 * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R. 399 * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x). 424 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 456 * This class implements forward and reverse operations for the function x -> x^p for use within CppAD. 457 * While CppAD would implement integer powers as a recursion of multiplications, we still use pow functions as they allow us to avoid overestimation in interval arithmetics. 459 * @todo treat the exponent as a (variable) argument to the function, with the assumption that we never differentiate w.r.t. it (this should make the approach threadsafe again) 501 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 502 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */ 536 // ty[2] = 1/2 * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1] * tx[1] + exponent * pow(tx[0], exponent-1) * tx[2]; 561 * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$ 562 * where x^[l] is the l'th taylor coefficient (x, x', x'', ...) and h(x) = g(y(x)) for some function g:R^k -> R. 565 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l]) 609 // px[0] = py[0] * exponent * pow(tx[0], exponent-1) + py[1] * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1]; 630 * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R. 646 * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x). 666 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 712 * This class implements forward and reverse operations for the function x -> sign(x)abs(x)^p for use within CppAD. 713 * While we otherwise would have to use discontinuous sign and abs functions, our own implementation allows to provide 716 * @todo treat the exponent as a (variable) argument to the function, with the assumption that we never differentiate w.r.t. it (this should make the approach threadsafe again) 758 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 759 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */ 817 * Then in the reverse sweep we have to compute the elements of \f$\partial h / \partial x^[l], l = 0, ..., k,\f$ 818 * where x^[l] is the l'th taylor coefficient (x, x', x'', ...) and h(x) = g(y(x)) for some function g:R^k -> R. 821 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l]) 867 // px[0] = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1] 899 * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R. 915 * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x). 935 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 982 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 983 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */ 1045 CppAD::vector<SCIPInterval>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */ 1065 // px[0] = py[0] * p * abs(tx[0])^(p-1) + py[1] * p * (p-1) * abs(tx[0])^(p-2) * sign(tx[0]) * tx[1] 1097 * For a 1 x q matrix R, we have to return the sparsity pattern of the 1 x q matrix S(x) = f'(x) * R. 1113 * For a q x 1 matrix R, we have to return the sparsity pattern of the q x 1 matrix S(x) = R * f'(x). 1129 * Assume V(x) = (g(f(x)))'' R with f(x) = sign(x)abs(x)^p for a function g:R->R and a matrix R. 1133 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 1236 * This class implements forward and reverse operations for a function given by a user expression for use within CppAD. 1278 * 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 1282 * 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) 1283 * 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]) 1288 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 1289 CppAD::vector<bool>& vy, /**< vector to store which function values depend on variables, or empty vector */ 1390 * 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)), 1430 CppAD::vector<Type>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */ 1493 * 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. 1522 * 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). 1552 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */ 1975 /* for each argument, we collect it's linear index from lincoefs, it's square coefficients and all factors from bilinear terms 1982 /* there are no quadratic terms with argidx in its first argument, that should be easy to handle */ 2100 * This may be the case if the evaluation sequence depends on values of operands (e.g., in case of abs, sign, signpower, ...). 2137 * 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. 2150 SCIPdebugMessage("ignore CppAD error from %sknown source %s:%d: msg: %s exp: %s\n", known ? "" : "un", file, line, msg, cond); 2165 return "Algorithmic Differentiation of C++ algorithms developed by B. Bell (www.coin-or.org/CppAD)"; 2259 * In cases of user-given expressions, higher order derivatives may not be available for the user-expression, 2260 * even if the expression interpreter could handle these. This method allows to recognize that, e.g., the 2261 * Hessian for an expression is not available because it contains a user expression that does not provide 2455 SCIP_Bool new_varvals, /**< have variable values changed since last call to a point evaluation routine? */ 2457 SCIP_Real* gradient /**< buffer to store expression gradient, need to have length at least SCIPexprtreeGetNVars(tree) */ 2489 SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n"); 2490 SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t %g", gradient[i]); printf("\n"); 2502 SCIP_INTERVAL* varvals, /**< interval values of variables, can be NULL if new_varvals is FALSE */ 2503 SCIP_Bool new_varvals, /**< have variable interval values changed since last call to an interval evaluation routine? */ 2505 SCIP_INTERVAL* gradient /**< buffer to store expression interval gradient, need to have length at least SCIPexprtreeGetNVars(tree) */ 2535 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"); 2536 SCIPdebugMessage("grad ="); for (int i = 0; i < n; ++i) printf("\t [%g,%g]", SCIPintervalGetInf(gradient[i]), SCIPintervalGetSup(gradient[i])); printf("\n"); 2552 SCIP_Bool* sparsity /**< buffer to store sparsity pattern of Hessian, sparsity[i+n*j] indicates whether entry (i,j) is nonzero in the hessian */ 2578 SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n"); 2609 SCIPdebugMessage("HessianSparsityDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n"); 2610 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"); 2625 SCIP_Bool new_varvals, /**< have variable values changed since last call to an evaluation routine? */ 2672 SCIPdebugMessage("HessianDense for "); SCIPexprtreePrint(tree, NULL, NULL, NULL); printf("\n"); 2673 SCIPdebugMessage("x ="); for (int i = 0; i < n; ++i) printf("\t %g", data->x[i]); printf("\n"); 2674 SCIPdebugMessage("hess ="); for (int i = 0; i < n*n; ++i) printf("\t %g", hessian[i]); printf("\n"); SCIP_RETCODE SCIPexprintEval(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Real *val) Definition: exprinterpret_cppad.cpp:2317 Definition: exprinterpret_cppad.cpp:462 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:2621 SCIP_EXPRINTCAPABILITY SCIPexprintGetExprtreeCapability(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree) Definition: exprinterpret_cppad.cpp:2264 Definition: type_expr.h:55 bool IdenticalEqualPar(const SCIPInterval &x, const SCIPInterval &y) Definition: exprinterpret_cppad.cpp:245 Definition: type_expr.h:57 methods to interpret (evaluate) an expression tree "fast" SCIP_RETCODE SCIPexprtreeEval(SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Real *val) Definition: type_expr.h:70 int SCIPexprGetNChildren(SCIP_EXPR *expr) void SCIPexprtreeSetInterpreterData(SCIP_EXPRTREE *tree, SCIP_EXPRINTDATA *interpreterdata) Definition: type_expr.h:63 SCIP_QUADELEM * SCIPexprGetQuadElements(SCIP_EXPR *expr) SCIP_EXPRINTCAPABILITY SCIPexprGetUserEvalCapability(SCIP_EXPR *expr) SCIP_Real SCIPexprGetMonomialCoef(SCIP_EXPRDATA_MONOMIAL *monomial) Definition: type_expr.h:72 struct SCIP_ExprData_Monomial SCIP_EXPRDATA_MONOMIAL Definition: type_expr.h:109 #define SCIP_EXPRINTCAPABILITY_INTGRADIENT Definition: type_exprinterpret.h:39 SCIP_RETCODE SCIPexprintFree(SCIP_EXPRINT **exprint) Definition: exprinterpret_cppad.cpp:2196 SCIP_Real * SCIPexprGetMonomialExponents(SCIP_EXPRDATA_MONOMIAL *monomial) int SCIPexprtreeGetNVars(SCIP_EXPRTREE *tree) 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:2451 Definition: type_expr.h:60 Definition: type_expr.h:71 Definition: intervalarithext.h:42 SCIP_Real SCIPexprGetSignPowerExponent(SCIP_EXPR *expr) SCIP_RETCODE SCIPexprCopyDeep(BMS_BLKMEM *blkmem, SCIP_EXPR **targetexpr, SCIP_EXPR *sourceexpr) SCIP_EXPRINTDATA * SCIPexprtreeGetInterpreterData(SCIP_EXPRTREE *tree) Definition: type_expr.h:38 int SCIPexprGetNMonomials(SCIP_EXPR *expr) SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void) Definition: exprinterpret_cppad.cpp:2169 Definition: exprinterpret_cppad.cpp:89 int SCIPexprGetNQuadElements(SCIP_EXPR *expr) Definition: type_expr.h:53 Definition: type_expr.h:44 SCIP_EXPR ** SCIPexprGetChildren(SCIP_EXPR *expr) public methods for expressions, expression trees, expression graphs, and related stuff ... #define SCIP_EXPRINTCAPABILITY_INTFUNCVALUE Definition: type_exprinterpret.h:37 void SCIPexprSortQuadElems(SCIP_EXPR *expr) SCIPInterval CondExpOp(enum CppAD::CompareOp cop, const SCIPInterval &left, const SCIPInterval &right, const SCIPInterval &trueCase, const SCIPInterval &falseCase) Definition: exprinterpret_cppad.cpp:191 Definition: type_expr.h:45 SCIPInterval signpow(const SCIPInterval &x, const SCIP_Real p) Definition: intervalarithext.h:316 Definition: type_expr.h:47 Definition: type_expr.h:74 Definition: type_expr.h:51 int * SCIPexprGetMonomialChildIndices(SCIP_EXPRDATA_MONOMIAL *monomial) Definition: type_expr.h:49 SCIP_Real SCIPexprGetQuadConstant(SCIP_EXPR *expr) SCIP_RETCODE SCIPexprintFreeData(SCIP_EXPRINTDATA **interpreterdata) Definition: exprinterpret_cppad.cpp:2278 void evalAbs(CppAD::AD< SCIPInterval > &resultant, const CppAD::AD< SCIPInterval > &arg) Definition: exprinterpret_cppad.cpp:1725 SCIP_EXPR * SCIPexprtreeGetRoot(SCIP_EXPRTREE *tree) Definition: type_expr.h:50 void evalMin(CppAD::AD< double > &resultant, const CppAD::AD< double > &arg1, const CppAD::AD< double > &arg2) Definition: exprinterpret_cppad.cpp:1656 SCIP_RETCODE SCIPexprtreeEvalInt(SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_INTERVAL *val) Definition: exprinterpret_cppad.cpp:1239 Definition: type_retcode.h:33 SCIP_RETCODE SCIPexprintEvalInt(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real infinity, SCIP_INTERVAL *varvals, SCIP_INTERVAL *val) Definition: exprinterpret_cppad.cpp:2385 SCIP_RETCODE SCIPexprEvalUser(SCIP_EXPR *expr, SCIP_Real *argvals, SCIP_Real *val, SCIP_Real *gradient, SCIP_Real *hessian) SCIP_EXPRDATA_MONOMIAL ** SCIPexprGetMonomials(SCIP_EXPR *expr) SCIP_RETCODE exprEvalUser(SCIP_EXPR *expr, Type *x, Type &funcval, Type *gradient, Type *hessian) Definition: exprinterpret_cppad.cpp:1211 SCIP_RETCODE SCIPexprintHessianSparsityDense(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree, SCIP_Real *varvals, SCIP_Bool *sparsity) Definition: exprinterpret_cppad.cpp:2548 SCIP_Real SCIPexprGetOpReal(SCIP_EXPR *expr) SCIP_RETCODE SCIPexprEvalIntUser(SCIP_EXPR *expr, SCIP_Real infinity, SCIP_INTERVAL *argvals, SCIP_INTERVAL *val, SCIP_INTERVAL *gradient, SCIP_INTERVAL *hessian) Definition: type_retcode.h:34 SCIP_Real SCIPexprGetRealPowerExponent(SCIP_EXPR *expr) int SCIPexprGetMonomialNFactors(SCIP_EXPRDATA_MONOMIAL *monomial) #define SCIP_EXPRINTCAPABILITY_HESSIAN Definition: type_exprinterpret.h:40 void SCIPexprFreeDeep(BMS_BLKMEM *blkmem, SCIP_EXPR **expr) Definition: type_expr.h:37 int SCIPexprGetOpIndex(SCIP_EXPR *expr) SCIP_Real SCIPexprGetPolynomialConstant(SCIP_EXPR *expr) Definition: type_expr.h:56 SCIP_Real * SCIPexprtreeGetParamVals(SCIP_EXPRTREE *tree) SCIP_RETCODE SCIPexprintNewParametrization(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree) Definition: exprinterpret_cppad.cpp:2298 Definition: type_expr.h:62 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:2498 Definition: type_expr.h:61 Definition: exprinterpret_cppad.cpp:719 bool EqualOpSeq(const SCIPInterval &x, const SCIPInterval &y) Definition: exprinterpret_cppad.cpp:208 void evalMax(CppAD::AD< double > &resultant, const CppAD::AD< double > &arg1, const CppAD::AD< double > &arg2) Definition: exprinterpret_cppad.cpp:1686 Definition: type_expr.h:77 #define SCIP_EXPRINTCAPABILITY_FUNCVALUE Definition: type_exprinterpret.h:36 SCIP_RETCODE SCIPexprintCreate(BMS_BLKMEM *blkmem, SCIP_EXPRINT **exprint) Definition: exprinterpret_cppad.cpp:2179 std::ostream & operator<<(std::ostream &out, const SCIP_INTERVAL &x) Definition: exprinterpret_cppad.cpp:325 SCIP_EXPROP SCIPexprGetOperator(SCIP_EXPR *expr) Definition: type_expr.h:69 int SCIPexprGetIntPowerExponent(SCIP_EXPR *expr) Definition: type_expr.h:52 SCIP_Real SCIPexprGetLinearConstant(SCIP_EXPR *expr) atomic_signpower() Definition: exprinterpret_cppad.cpp:954 Definition: type_expr.h:73 Definition: type_expr.h:54 SCIP_Real * SCIPexprGetLinearCoefs(SCIP_EXPR *expr) Definition: type_retcode.h:35 common defines and data types used in all packages of SCIP Definition: type_expr.h:48 bool GreaterThanOrZero(const SCIPInterval &x) Definition: exprinterpret_cppad.cpp:269 #define SCIP_EXPRINTCAPABILITY_GRADIENT Definition: type_exprinterpret.h:38 SCIP_Real * SCIPexprGetQuadLinearCoefs(SCIP_EXPR *expr) C++ extensions to interval arithmetics for provable bounds. SCIP_RETCODE SCIPexprintCompile(SCIP_EXPRINT *exprint, SCIP_EXPRTREE *tree) Definition: exprinterpret_cppad.cpp:2209 Definition: type_expr.h:46 Definition: type_expr.h:39 |