exprinterpret_cppad.cpp
Go to the documentation of this file.
31/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
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
56/* defining EVAL_USE_EXPRHDLR_ALWAYS uses the exprhdlr methods for the evaluation and differentiation of all expression (except for var and value)
58 * this is incomplete (no hessians) and slow (only forward-mode gradients), but can be useful for testing
62/* defining NO_CPPAD_USER_ATOMIC disables the use of our own implementation of derivatives of power operators
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)
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.
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
107 * To run it in a multithreading environment, a special CppAD memory allocator that is aware of the multiple threads has to be used.
109 * To implement this, we follow the team_pthread example of CppAD, which uses pthread's thread-specific data management.
135 /* if no thread_number for this thread yet, then assign a new thread number to the current thread
167 * The purpose is to make sure that init_parallel() is called before any multithreading is started.
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 */
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 */
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 */
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.
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).
283 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
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.
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)
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 */
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];
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.
425 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
469 // px[0] = py[0] * exponent * pow(tx[0], exponent-1) + py[1] * exponent * (exponent-1) * pow(tx[0], exponent-2) * tx[1];
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.
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).
526 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
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
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)
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 */
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.
688 * px[l] = \partial h / \partial x^[l] = (\partial g / \partial y) * (\partial y / \partial x^[l])
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]
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.
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).
802 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
858 /* pow(0,fractional>1) is not differential in the CppAD world (https://github.com/coin-or/CppAD/discussions/93)
861 resultant = CppAD::CondExpEq(arg, adzero, pow(arg+std::numeric_limits<SCIP_Real>::epsilon(), exponent)-pow(std::numeric_limits<SCIP_Real>::epsilon(), exponent),
869 * This class implements forward and reverse operations for a function given by a user expression for use within CppAD.
878 : CppAD::atomic_base<SCIP_Real>(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr_)), CppAD::atomic_base<SCIP_Real>::bool_sparsity_enum),
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
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])
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 */
930 SCIPdebugMsg(scip, "expr_%s:forward, q=%zd, p=%zd\n", SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), q, p);
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)),
1106 CppAD::vector<SCIP_Real>& px, /**< vector to store partial derivatives of h(x) = g(y(x)) w.r.t. x */
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.
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).
1238 const CppAD::vector<bool>& vx, /**< indicates whether argument is a variable, or empty vector */
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),
1431 * but then derivatives at 0 are 0, while they are actually infinite (see expr_entropy.c:bwdiff)
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");
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.
1471 SCIPdebugMessage("ignore CppAD error from %sknown source %s:%d: msg: %s exp: %s\n", known ? "" : "un", file, line, msg, cond);
1486 return "Algorithmic Differentiation of C++ algorithms developed by B. Bell (github.com/coin-or/CppAD)";
1494 return SCIP_EXPRINTCAPABILITY_FUNCVALUE | SCIP_EXPRINTCAPABILITY_GRADIENT | SCIP_EXPRINTCAPABILITY_HESSIAN;
1505 *exprint = (SCIP_EXPRINT*)1u; /* some code checks that a non-NULL pointer is returned here, even though it may not point anywhere */
1559 for( expr = SCIPexpriterGetCurrent(it); !SCIPexpriterIsEnd(it); expr = SCIPexpriterGetNext(it) )
1588 * however, only atomic_userexpr:forward() is implemented for p=0,1 at the moment, so we cannot do Hessians
1592 (*exprintdata)->userevalcapability &= SCIP_EXPRINTCAPABILITY_FUNCVALUE | SCIP_EXPRINTCAPABILITY_GRADIENT;
1604 (*exprintdata)->varidxs.insert((*exprintdata)->varidxs.begin(), varidxs.begin(), varidxs.end());
1625 if( SCIPisExprPower(scip, child) && SCIPgetExponentExprPow(child) == 2.0 && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[0]) )
1628 if( SCIPisExprProduct(scip, child) && SCIPexprGetNChildren(child) == 2 && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[0]) && SCIPisExprVaridx(scip, SCIPexprGetChildren(child)[1]) )
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
1715 exprintdata->x[i] = varvals[idx]; /* need this for a following grad or hessian eval with new_x = false */
1719 for( vector<atomic_userexpr*>::iterator it(exprintdata->userexprs.begin()); it != exprintdata->userexprs.end(); ++it )
1762 SCIP_Bool new_varvals, /**< have variable values changed since last call to a point evaluation routine? */
1821 * Since the AD code might need to do a forward sweep, variable values need to be passed in here.
1893 // following https://github.com/coin-or/CppAD/blob/20180000.0/cppad_ipopt/src/vec_fun_pattern.cpp
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, ""),
1949 // we need hessparsity_pattern only if SCIPexprintHessian is doing a sparse hessian eval; nn/4 is the threshold for that
1955 // add above-diagonal elements into end part of hessparsity_row/col, so we have a 1:1 correspondence
1979 SCIPinfoMessage(scip, NULL, " (%d,%d)", exprintdata->hesrowidxs[i], exprintdata->hescolidxs[i]);
1994 * Returned arrays `rowidxs` and `colidxs` and number of elements `nnz` are the same as given by SCIPexprintHessianSparsity().
2003 SCIP_Bool new_varvals, /**< have variable values changed since last call to an evaluation routine? */
2024 SCIP_CALL( SCIPexprintHessianSparsity(scip, exprint, expr, exprintdata, varvals, &dummy2, &dummy2, &dummy1) );
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) )
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
2052 // use dense Hessian if there are many nonzero elements (approx more than half (recall only lower (or upper)-triangular is considered))
2056 exprintdata->hesvalues[i] = hess[exprintdata->hessparsity_row[i] * n + exprintdata->hessparsity_col[i]];
2060 // originally, this was hess = exprintdata->f.SparseHessian(exprintdata->x, vector<double>(1, 1.0), exprintdata->hessparsity),
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);
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
2085 SCIPinfoMessage(scip, NULL, " (%d,%d)=%g", exprintdata->hesrowidxs[i], exprintdata->hescolidxs[i], exprintdata->hesvalues[i]);
Definition: exprinterpret_cppad.cpp:872
atomic_userexpr(SCIP *scip_, SCIP_EXPR *expr_)
Definition: exprinterpret_cppad.cpp:874
common defines and data types used in all packages of SCIP
exponential expression handler
logarithm expression handler
power and signed power expression handlers
handler for variable index expressions
methods to interpret (evaluate) an expression "fast"
void posintpower(const vector< Type > &in, vector< Type > &out, size_t exponent)
Definition: exprinterpret_cppad.cpp:556
static CppAD::ErrorHandler errorhandler(cppaderrorcallback)
static void cppaderrorcallback(bool known, int line, const char *file, const char *cond, const char *msg)
Definition: exprinterpret_cppad.cpp:1463
static void evalIntPower(Type &resultant, const Type &arg, const int exponent)
Definition: exprinterpret_cppad.cpp:1287
static void evalSignPower(CppAD::AD< Type > &resultant, const CppAD::AD< Type > &arg, SCIP_EXPR *expr)
Definition: exprinterpret_cppad.cpp:840
static SCIP_RETCODE eval(SCIP *scip, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, const vector< Type > &x, Type &val)
Definition: exprinterpret_cppad.cpp:1334
SCIP_Bool SCIPisExprVaridx(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_varidx.c:252
SCIP_Bool SCIPisExprSignpower(SCIP *scip, SCIP_EXPR *expr)
Definition: expr_pow.c:3242
SCIP_RETCODE SCIPexprintCompile(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *rootexpr, SCIP_EXPRINTDATA **exprintdata)
Definition: exprinterpret_cppad.cpp:1530
SCIP_RETCODE SCIPexprintFreeData(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA **exprintdata)
Definition: exprinterpret_cppad.cpp:1641
SCIP_EXPRINTCAPABILITY SCIPexprintGetCapability(void)
Definition: exprinterpret_cppad.cpp:1490
SCIP_RETCODE SCIPexprintHessianSparsity(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, int **rowidxs, int **colidxs, int *nnz)
Definition: exprinterpret_cppad.cpp:1825
SCIP_RETCODE SCIPexprintFree(SCIP *scip, SCIP_EXPRINT **exprint)
Definition: exprinterpret_cppad.cpp:1511
SCIP_RETCODE SCIPexprintEval(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata, SCIP_Real *varvals, SCIP_Real *val)
Definition: exprinterpret_cppad.cpp:1680
SCIP_EXPRINTCAPABILITY SCIPexprintGetExprCapability(SCIP *scip, SCIP_EXPRINT *exprint, SCIP_EXPR *expr, SCIP_EXPRINTDATA *exprintdata)
Definition: exprinterpret_cppad.cpp:1667
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)
Definition: exprinterpret_cppad.cpp:1997
SCIP_RETCODE SCIPexprintCreate(SCIP *scip, SCIP_EXPRINT **exprint)
Definition: exprinterpret_cppad.cpp:1498
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)
Definition: exprinterpret_cppad.cpp:1756
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
SCIP_RETCODE SCIPevalExpr(SCIP *scip, SCIP_EXPR *expr, SCIP_SOL *sol, SCIP_Longint soltag)
Definition: scip_expr.c:1635
SCIP_Bool SCIPisExprProduct(SCIP *scip, SCIP_EXPR *expr)
Definition: scip_expr.c:1464
SCIP_Real SCIPgetCoefExprProduct(SCIP_EXPR *expr)
Definition: expr_product.c:2300
SCIP_EXPR * SCIPexpriterGetCurrent(SCIP_EXPRITER *iterator)
Definition: expriter.c:683
SCIP_RETCODE SCIPcallExprEval(SCIP *scip, SCIP_EXPR *expr, SCIP_Real *childrenvalues, SCIP_Real *val)
Definition: scip_expr.c:2185
SCIP_RETCODE SCIPcreateExpriter(SCIP *scip, SCIP_EXPRITER **iterator)
Definition: scip_expr.c:2337
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_RETCODE SCIPprintExpr(SCIP *scip, SCIP_EXPR *expr, FILE *file)
Definition: scip_expr.c:1486
SCIP_EXPR * SCIPexpriterGetNext(SCIP_EXPRITER *iterator)
Definition: expriter.c:858
SCIP_RETCODE SCIPexpriterInit(SCIP_EXPRITER *iterator, SCIP_EXPR *expr, SCIP_EXPRITER_TYPE type, SCIP_Bool allowrevisit)
Definition: expriter.c:501
#define SCIPfreeBlockMemoryArrayNull(scip, ptr, num)
Definition: scip_mem.h:111
interval arithmetics for provable bounds
Definition: objbenders.h:44
public functions to work with algebraic expressions
public functions to work with algebraic expressions
Definition: struct_expr.h:204
Definition: struct_expr.h:106
Definition: struct_scip.h:70
#define SCIP_EXPRINTCAPABILITY_GRADIENT
Definition: type_exprinterpret.h:51
#define SCIP_EXPRINTCAPABILITY_HESSIAN
Definition: type_exprinterpret.h:52
#define SCIP_EXPRINTCAPABILITY_FUNCVALUE
Definition: type_exprinterpret.h:50