1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the program and library */ 4 /* SCIP --- Solving Constraint Integer Programs */ 5 /* */ 6 /* Copyright (c) 2002-2023 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 expr_erf.c 26 * @brief handler for Gaussian error function expressions 27 * @author Benjamin Mueller 28 */ 29 30 /*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/ 31 32 #include "scip/expr_erf.h" 33 #include "scip/expr_value.h" 34 35 /* fundamental expression handler properties */ 36 #define EXPRHDLR_NAME "erf" 37 #define EXPRHDLR_DESC "Gaussian error function" 38 #define EXPRHDLR_PRECEDENCE 79000 39 #define EXPRHDLR_HASHKEY SCIPcalcFibHash(131071.0) 40 41 /* 42 * Data structures 43 */ 44 45 /* 46 * Local methods 47 */ 48 49 /** evaluates the Gaussian error function at a given point */ 50 static 51 SCIP_Real errorf( 52 SCIP_Real x /**< point to evaluate */ 53 ) 54 { 55 SCIP_Real a1 = +0.254829592; 56 SCIP_Real a2 = -0.284496736; 57 SCIP_Real a3 = +1.421413741; 58 SCIP_Real a4 = -1.453152027; 59 SCIP_Real a5 = +1.061405429; 60 SCIP_Real p = +0.3275911; 61 int sign = (x >= 0.0) ? 1 : -1; 62 SCIP_Real t = 1.0 / (1.0 + p * REALABS(x)); 63 SCIP_Real y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * exp(-x*x); 64 65 return sign*y; 66 } 67 68 /* 69 * Callback methods of expression handler 70 */ 71 72 /** expression handler copy callback */ 73 static 74 SCIP_DECL_EXPRCOPYHDLR(copyhdlrErf) 75 { /*lint --e{715}*/ 76 SCIP_CALL( SCIPincludeExprhdlrErf(scip) ); 77 78 return SCIP_OKAY; 79 } 80 81 /** simplifies an erf expression */ 82 static 83 SCIP_DECL_EXPRSIMPLIFY(simplifyErf) 84 { /*lint --e{715}*/ 85 SCIP_EXPR* child; 86 87 assert(scip != NULL); 88 assert(expr != NULL); 89 assert(simplifiedexpr != NULL); 90 assert(SCIPexprGetNChildren(expr) == 1); 91 92 child = SCIPexprGetChildren(expr)[0]; 93 assert(child != NULL); 94 95 /* check for value expression */ 96 if( SCIPisExprValue(scip, child) ) 97 { 98 SCIP_CALL( SCIPcreateExprValue(scip, simplifiedexpr, errorf(SCIPgetValueExprValue(child)), ownercreate, ownercreatedata) ); 99 } 100 else 101 { 102 *simplifiedexpr = expr; 103 104 /* we have to capture it, since it must simulate a "normal" simplified call in which a new expression is created */ 105 SCIPcaptureExpr(*simplifiedexpr); 106 } 107 108 return SCIP_OKAY; 109 } 110 111 /** expression parse callback */ 112 static 113 SCIP_DECL_EXPRPARSE(parseErf) 114 { /*lint --e{715}*/ 115 SCIP_EXPR* childexpr; 116 117 assert(expr != NULL); 118 119 /* parse child expression from remaining string */ 120 SCIP_CALL( SCIPparseExpr(scip, &childexpr, string, endstring, ownercreate, ownercreatedata) ); 121 assert(childexpr != NULL); 122 123 /* create gaussian error function expression */ 124 SCIP_CALL( SCIPcreateExprErf(scip, expr, childexpr, ownercreate, ownercreatedata) ); 125 assert(*expr != NULL); 126 127 /* release child expression since it has been captured by the gaussian error function expression */ 128 SCIP_CALL( SCIPreleaseExpr(scip, &childexpr) ); 129 130 *success = TRUE; 131 132 return SCIP_OKAY; 133 } 134 135 /** expression (point-) evaluation callback */ 136 static 137 SCIP_DECL_EXPREVAL(evalErf) 138 { /*lint --e{715}*/ 139 assert(expr != NULL); 140 assert(SCIPexprGetData(expr) == NULL); 141 assert(SCIPexprGetNChildren(expr) == 1); 142 assert(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0]) != SCIP_INVALID); /*lint !e777*/ 143 144 *val = errorf(SCIPexprGetEvalValue(SCIPexprGetChildren(expr)[0])); 145 146 return SCIP_OKAY; 147 } 148 149 /** expression derivative evaluation callback */ 150 static 151 SCIP_DECL_EXPRBWDIFF(bwdiffErf) 152 { /*lint --e{715}*/ 153 assert(expr != NULL); 154 155 SCIPerrorMessage("method of erf expression handler not implemented yet\n"); 156 SCIPABORT(); /*lint --e{527}*/ 157 158 return SCIP_OKAY; 159 } 160 161 /** expression interval evaluation callback */ 162 static 163 SCIP_DECL_EXPRINTEVAL(intevalErf) 164 { /*lint --e{715}*/ 165 SCIP_INTERVAL childinterval; 166 167 assert(expr != NULL); 168 assert(SCIPexprGetData(expr) == NULL); 169 assert(SCIPexprGetNChildren(expr) == 1); 170 171 childinterval = SCIPexprGetActivity(SCIPexprGetChildren(expr)[0]); 172 173 if( SCIPintervalIsEmpty(SCIP_INTERVAL_INFINITY, childinterval) ) 174 SCIPintervalSetEmpty(interval); 175 else 176 { 177 SCIP_Real childinf = SCIPintervalGetInf(childinterval); 178 SCIP_Real childsup = SCIPintervalGetSup(childinterval); 179 SCIP_Real inf = childinf <= -SCIP_INTERVAL_INFINITY ? -1.0 : errorf(childinf); 180 SCIP_Real sup = childsup >= +SCIP_INTERVAL_INFINITY ? +1.0 : errorf(childsup); 181 assert(inf <= sup); 182 SCIPintervalSetBounds(interval, inf, sup); 183 } 184 185 return SCIP_OKAY; 186 } 187 188 /** erf hash callback */ 189 static 190 SCIP_DECL_EXPRHASH(hashErf) 191 { /*lint --e{715}*/ 192 assert(scip != NULL); 193 assert(expr != NULL); 194 assert(SCIPexprGetNChildren(expr) == 1); 195 assert(hashkey != NULL); 196 assert(childrenhashes != NULL); 197 198 *hashkey = EXPRHDLR_HASHKEY; 199 *hashkey ^= childrenhashes[0]; 200 201 return SCIP_OKAY; 202 } 203 204 /** expression curvature detection callback */ 205 static 206 SCIP_DECL_EXPRCURVATURE(curvatureErf) 207 { /*lint --e{715}*/ 208 assert(scip != NULL); 209 assert(expr != NULL); 210 211 /* expression is 212 * - convex if child is convex and child <= 0 213 * - concave if child is concave and child >= 0 214 */ 215 if( exprcurvature == SCIP_EXPRCURV_CONVEX ) 216 { 217 *success = TRUE; 218 *childcurv = SCIP_EXPRCURV_CONVEX; 219 } 220 else 221 *success = FALSE; 222 223 return SCIP_OKAY; 224 } 225 226 /** expression monotonicity detection callback */ 227 static 228 SCIP_DECL_EXPRMONOTONICITY(monotonicityErf) 229 { /*lint --e{715}*/ 230 assert(scip != NULL); 231 assert(expr != NULL); 232 assert(result != NULL); 233 234 *result = SCIP_MONOTONE_INC; 235 236 return SCIP_OKAY; 237 } 238 239 /** expression integrality detection callback */ 240 static 241 SCIP_DECL_EXPRINTEGRALITY(integralityErf) 242 { /*lint --e{715}*/ 243 assert(scip != NULL); 244 assert(expr != NULL); 245 assert(isintegral != NULL); 246 247 *isintegral = FALSE; 248 249 return SCIP_OKAY; 250 } 251 252 /** creates an erf expression 253 * 254 * @attention The implementation of `erf` expressions is incomplete. 255 * They are not usable for most use cases so far. 256 */ 257 SCIP_RETCODE SCIPcreateExprErf( 258 SCIP* scip, /**< SCIP data structure */ 259 SCIP_EXPR** expr, /**< pointer where to store expression */ 260 SCIP_EXPR* child, /**< child expression */ 261 SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), /**< function to call to create ownerdata */ 262 void* ownercreatedata /**< data to pass to ownercreate */ 263 ) 264 { 265 SCIP_EXPRHDLR* exprhdlr; 266 267 assert(expr != NULL); 268 assert(child != NULL); 269 270 exprhdlr = SCIPfindExprhdlr(scip, EXPRHDLR_NAME); 271 if( exprhdlr == NULL ) 272 { 273 SCIPerrorMessage("could not find %s expression handler -> abort\n", EXPRHDLR_NAME); 274 SCIPABORT(); 275 return SCIP_PLUGINNOTFOUND; 276 } 277 278 /* create expression */ 279 SCIP_CALL( SCIPcreateExpr(scip, expr, exprhdlr, NULL, 1, &child, ownercreate, ownercreatedata) ); 280 281 return SCIP_OKAY; 282 } 283 284 /** indicates whether expression is of erf-type */ /*lint -e{715}*/ 285 SCIP_Bool SCIPisExprErf( 286 SCIP* scip, /**< SCIP data structure */ 287 SCIP_EXPR* expr /**< expression */ 288 ) 289 { /*lint --e{715}*/ 290 assert(expr != NULL); 291 292 return strcmp(SCIPexprhdlrGetName(SCIPexprGetHdlr(expr)), EXPRHDLR_NAME) == 0; 293 } 294 295 /** creates the handler for erf expressions and includes it into SCIP 296 * 297 * @attention The implementation of this expression handler is incomplete. 298 * It is not usable for most use cases so far. 299 */ 300 SCIP_RETCODE SCIPincludeExprhdlrErf( 301 SCIP* scip /**< SCIP data structure */ 302 ) 303 { 304 SCIP_EXPRHDLR* exprhdlr; 305 306 /* include expression handler */ 307 SCIP_CALL( SCIPincludeExprhdlr(scip, &exprhdlr, EXPRHDLR_NAME, EXPRHDLR_DESC, EXPRHDLR_PRECEDENCE, evalErf, NULL) ); 308 assert(exprhdlr != NULL); 309 310 SCIPexprhdlrSetCopyFreeHdlr(exprhdlr, copyhdlrErf, NULL); 311 SCIPexprhdlrSetSimplify(exprhdlr, simplifyErf); 312 SCIPexprhdlrSetParse(exprhdlr, parseErf); 313 SCIPexprhdlrSetIntEval(exprhdlr, intevalErf); 314 SCIPexprhdlrSetHash(exprhdlr, hashErf); 315 SCIPexprhdlrSetDiff(exprhdlr, bwdiffErf, NULL, NULL); 316 SCIPexprhdlrSetCurvature(exprhdlr, curvatureErf); 317 SCIPexprhdlrSetMonotonicity(exprhdlr, monotonicityErf); 318 SCIPexprhdlrSetIntegrality(exprhdlr, integralityErf); 319 320 return SCIP_OKAY; 321 } 322