SCIP Doxygen Documentation
 
Loading...
Searching...
No Matches
exprcurv.c
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-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 exprcurv.c
26 * @ingroup OTHER_CFILES
27 * @brief functions to work with curvature (convex, concave, etc)
28 * @author Stefan Vigerske
29 *
30 * Declarations are in pub_expr.h
31 */
32
33#include "scip/pub_expr.h"
34
35/** curvature names as strings */
36static
37const char* curvnames[4] =
38 {
39 "unknown",
40 "convex",
41 "concave",
42 "linear"
43 };
44
45#ifdef NDEBUG
46#undef SCIPexprcurvAdd
47#undef SCIPexprcurvNegate
48#undef SCIPexprcurvMultiply
49#endif
50
51/** gives curvature for a sum of two functions with given curvature */
53 SCIP_EXPRCURV curv1, /**< curvature of first summand */
54 SCIP_EXPRCURV curv2 /**< curvature of second summand */
55 )
56{
57 return (SCIP_EXPRCURV) (curv1 & curv2);
58}
59
60/** gives the curvature for the negation of a function with given curvature */
62 SCIP_EXPRCURV curvature /**< curvature of function */
63 )
64{
65 switch( curvature )
66 {
69
72
75 /* can return curvature, do this below */
76 break;
77
78 default:
79 SCIPerrorMessage("unknown curvature status.\n");
80 SCIPABORT();
81 }
82
83 return curvature;
84}
85
86/** gives curvature for a functions with given curvature multiplied by a constant factor */
88 SCIP_Real factor, /**< constant factor */
89 SCIP_EXPRCURV curvature /**< curvature of other factor */
90 )
91{
92 if( factor == 0.0 )
94 if( factor > 0.0 )
95 return curvature;
96 return SCIPexprcurvNegate(curvature);
97}
98
99/** gives curvature for base^exponent for given bounds and curvature of base-function and constant exponent */
101 SCIP_INTERVAL basebounds, /**< bounds on base function */
102 SCIP_EXPRCURV basecurv, /**< curvature of base function */
103 SCIP_Real exponent /**< exponent */
104 )
105{
106 SCIP_Bool expisint;
107
108 assert(basebounds.inf <= basebounds.sup);
109
110 if( exponent == 0.0 )
112
113 if( exponent == 1.0 )
114 return basecurv;
115
116 expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
117
118 /* if exponent is fractional, then power is not defined for a negative base
119 * thus, consider only positive part of basebounds
120 */
121 if( !expisint && basebounds.inf < 0.0 )
122 {
123 basebounds.inf = 0.0;
124 if( basebounds.sup < 0.0 )
126 }
127
128 /* if basebounds contains 0.0, consider negative and positive interval separately, if possible */
129 if( basebounds.inf < 0.0 && basebounds.sup > 0.0 )
130 {
133
134 /* something like x^(-2) may look convex on each side of zero, but is not convex on the whole interval
135 * due to the singularity at 0.0 */
136 if( exponent < 0.0 )
138
141
143 }
144 assert(basebounds.inf >= 0.0 || basebounds.sup <= 0.0);
145
146 /* (base^exponent)'' = exponent * ( (exponent-1) base^(exponent-2) (base')^2 + base^(exponent-1) base'' )
147 *
148 * if base'' is positive, i.e., base is convex, then
149 * - for base > 0.0 and exponent > 1.0, the second deriv. is positive -> convex
150 * - for base < 0.0 and exponent > 1.0, we can't say (first and second summand opposite signs)
151 * - for base > 0.0 and 0.0 < exponent < 1.0, we can't say (first sommand negative, second summand positive)
152 * - for base > 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
153 * - for base < 0.0 and exponent < 0.0 and even, the second deriv. is positive -> convex
154 * - for base < 0.0 and exponent < 0.0 and odd, the second deriv. is negative -> concave
155 *
156 * if base'' is negative, i.e., base is concave, then
157 * - for base > 0.0 and exponent > 1.0, we can't say (first summand positive, second summand negative)
158 * - for base < 0.0 and exponent > 1.0 and even, the second deriv. is positive -> convex
159 * - for base < 0.0 and exponent > 1.0 and odd, the second deriv. is negative -> concave
160 * - for base > 0.0 and 0.0 < exponent < 1.0, the second deriv. is negative -> concave
161 * - for base > 0.0 and exponent < 0.0, the second deriv. is positive -> convex
162 * - for base < 0.0 and exponent < 0.0, we can't say (first and second summand opposite signs)
163 *
164 * if base'' is zero, i.e., base is linear, then
165 * (base^exponent)'' = exponent * (exponent-1) base^(exponent-2) (base')^2
166 * - just multiply signs
167 */
168
170 {
171 SCIP_Real sign;
172
173 /* base^(exponent-2) is negative, if base < 0.0 and exponent is odd */
174 sign = exponent * (exponent - 1.0);
175 assert(basebounds.inf >= 0.0 || expisint);
176 if( basebounds.inf < 0.0 && ((int)exponent)%2 != 0 )
177 sign *= -1.0;
178 assert(sign != 0.0);
179
181 }
182
184 {
185 if( basebounds.sup <= 0.0 && exponent < 0.0 && expisint )
186 return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
187 if( basebounds.inf >= 0.0 && exponent > 1.0 )
188 return SCIP_EXPRCURV_CONVEX ;
190 }
191
193 {
194 if( basebounds.sup <= 0.0 && exponent > 1.0 && expisint )
195 return ((int)exponent)%2 == 0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
196 if( basebounds.inf >= 0.0 && exponent < 1.0 )
197 return exponent < 0.0 ? SCIP_EXPRCURV_CONVEX : SCIP_EXPRCURV_CONCAVE;
199 }
200
202}
203
204/** gives required curvature for base so that base^exponent has given curvature under given bounds on base and constant exponent
205 *
206 * returns curvature unknown if expected curvature cannot be obtained
207 */
209 SCIP_INTERVAL basebounds, /**< bounds on base function */
210 SCIP_Real exponent, /**< exponent, must not be 0 */
211 SCIP_EXPRCURV powercurv /**< expected curvature for power */
212 )
213{
214 SCIP_Bool expisint;
215
216 assert(basebounds.inf <= basebounds.sup);
217 assert(exponent != 0.0);
219
220 if( exponent == 1.0 )
221 return powercurv;
222
223 /* power is usually never linear, now that exponent != 1 */
226
227 expisint = EPSISINT(exponent, 0.0); /*lint !e835*/
228
229 /* if exponent is fractional, then power is only defined for a non-negative base
230 * boundtightening should have ensured this before calling this function,
231 * but sometimes this does not work and so we correct this here for us
232 */
233 if( !expisint && basebounds.inf < 0.0 )
234 {
235 basebounds.inf = 0.0;
236 if( basebounds.sup < 0.0 )
238 }
239
240 /* if basebounds contains 0.0, consider negative and positive interval separately, if possible */
241 if( basebounds.inf < 0.0 && basebounds.sup > 0.0 )
242 {
247
248 /* something like x^(-2) may look convex on each side of zero, but is not convex on the whole
249 * interval due to the singularity at 0.0 */
250 if( exponent < 0.0 )
252
255
258
259 /* now need to intersect */
261 return rightcurv;
263 return leftcurv;
269 }
270 assert(basebounds.inf >= 0.0 || basebounds.sup <= 0.0);
271
272 /* inverting the logic from SCIPexprcurvPower here */
274 {
275 SCIP_Real sign;
276
277 if( basebounds.sup <= 0.0 && exponent < 0.0 && expisint && ((int)exponent)%2 == 0 )
279 if( basebounds.inf >= 0.0 && exponent > 1.0 )
281 if( basebounds.sup <= 0.0 && exponent > 1.0 && expisint && ((int)exponent)%2 == 0 )
283 if( basebounds.inf >= 0.0 && exponent < 0.0 )
285
286 /* base^(exponent-2) is negative, if base < 0.0 and exponent is odd */
287 sign = exponent * (exponent - 1.0);
288 assert(basebounds.inf >= 0.0 || expisint);
289 if( basebounds.inf < 0.0 && ((int)exponent)%2 != 0 )
290 sign *= -1.0;
291 assert(sign != 0.0);
292
293 if( sign > 0.0 )
295 }
296 else
297 {
298 SCIP_Real sign;
299
300 assert(powercurv == SCIP_EXPRCURV_CONCAVE); /* linear handled at top, unknown should not be the case */
301
302 if( basebounds.sup <= 0.0 && exponent < 0.0 && expisint && ((int)exponent)%2 != 0 )
304 if( basebounds.sup <= 0.0 && exponent > 1.0 && expisint && ((int)exponent)%2 != 0 )
306 if( basebounds.inf >= 0.0 && exponent < 1.0 && exponent >= 0.0 )
308
309 /* base^(exponent-2) is negative, if base < 0.0 and exponent is odd */
310 sign = exponent * (exponent - 1.0);
311 assert(basebounds.inf >= 0.0 || expisint);
312 if( basebounds.inf < 0.0 && ((int)exponent)%2 != 0 )
313 sign *= -1.0;
314 assert(sign != 0.0);
315
316 if( sign < 0.0 )
318 }
319
321}
322
323/** gives curvature for a monomial with given curvatures and bounds for each factor
324 *
325 * See Maranas and Floudas, Finding All Solutions of Nonlinearly Constrained Systems of Equations, JOGO 7, 1995
326 * for the categorization in the case that all factors are linear.
327 *
328 * Exponents can also be negative or rational.
329 */
331 int nfactors, /**< number of factors in monomial */
332 SCIP_Real* exponents, /**< exponents in monomial, or NULL if all 1.0 */
333 int* factoridxs, /**< indices of factors (but not exponents), or NULL if identity mapping */
334 SCIP_EXPRCURV* factorcurv, /**< curvature of each factor */
335 SCIP_INTERVAL* factorbounds /**< bounds of each factor */
336 )
337{
338 SCIP_Real mult;
339 SCIP_Real e;
340 SCIP_INTERVAL bounds;
341 SCIP_EXPRCURV curv;
343 int nnegative;
344 int npositive;
345 SCIP_Real sum;
346 SCIP_Bool expcurvpos;
347 SCIP_Bool expcurvneg;
348 int j;
349 int f;
350
351 assert(nfactors >= 0);
352 assert(factorcurv != NULL || nfactors == 0);
353 assert(factorbounds != NULL || nfactors == 0);
354
355 if( nfactors == 0 )
357
358 if( nfactors == 1 )
359 {
360 f = factoridxs != NULL ? factoridxs[0] : 0;
361 e = exponents != NULL ? exponents[0] : 1.0;
362 /* SCIPdebugMessage("monomial [%g,%g]^%g is %s\n",
363 factorbounds[f].inf, factorbounds[f].sup, e,
364 SCIPexprcurvGetName(SCIPexprcurvPower(factorbounds[f], factorcurv[f], e))); */
365 return SCIPexprcurvPower(factorbounds[f], factorcurv[f], e); /*lint !e613*/
366 }
367
368 mult = 1.0;
369
370 nnegative = 0; /* number of negative exponents */
371 npositive = 0; /* number of positive exponents */
372 sum = 0.0; /* sum of exponents */
373 expcurvpos = TRUE; /* whether exp_j * f_j''(x) >= 0 for all factors (assuming f_j >= 0) */
374 expcurvneg = TRUE; /* whether exp_j * f_j''(x) <= 0 for all factors (assuming f_j >= 0) */
375
376 for( j = 0; j < nfactors; ++j )
377 {
378 f = factoridxs != NULL ? factoridxs[j] : j;
379 if( factorcurv[f] == SCIP_EXPRCURV_UNKNOWN ) /*lint !e613*/
381
382 e = exponents != NULL ? exponents[j] : 1.0;
383 bounds = factorbounds[f]; /*lint !e613*/
384
385 /* if argument is negative, then exponent should be integer; correct bounds if that doesn't hold */
386 if( !EPSISINT(e, 0.0) && bounds.inf < 0.0 ) /*lint !e835*/
387 {
388 bounds.inf = 0.0;
389 if( bounds.sup < 0.0 )
391 }
392
393 if( bounds.inf < 0.0 && bounds.sup > 0.0 )
395
396 if( e < 0.0 )
397 ++nnegative;
398 else
399 ++npositive;
400 sum += e;
401
402 if( bounds.inf < 0.0 )
403 {
404 /* flip j'th argument: (f_j)^(exp_j) = (-1)^(exp_j) (-f_j)^(exp_j) */
405
406 /* -f_j has negated curvature of f_j */
407 fcurv = SCIPexprcurvNegate(factorcurv[f]); /*lint !e613*/
408
409 /* negate monomial, if exponent is odd, i.e., (-1)^(exp_j) = -1 */
410 if( (int)e % 2 != 0 )
411 mult *= -1.0;
412 }
413 else
414 {
415 fcurv = factorcurv[f]; /*lint !e613*/
416 }
417
418 /* check if exp_j * fcurv is convex (>= 0) and/or concave */
420 if( !(fcurv & SCIP_EXPRCURV_CONVEX) )
424 }
425
426 /* if all factors are linear, then a product f_j^exp_j with f_j >= 0 is convex if
427 * - all exponents are negative, or
428 * - all except one exponent j* are negative and exp_j* >= 1 - sum_{j!=j*}exp_j, but the latter is equivalent to sum_j exp_j >= 1
429 * further, the product is concave if
430 * - all exponents are positive and the sum of exponents is <= 1.0
431 *
432 * if factors are nonlinear, then we require additionally, that for convexity
433 * - each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0
434 * and for concavity, we require that
435 * - all factors are concave, i.e., exp_j*f_j'' <= 0
436 */
437
438 if( nnegative == nfactors && expcurvpos )
440 else if( nnegative == nfactors-1 && EPSGE(sum, 1.0, 1e-9) && expcurvpos )
442 else if( npositive == nfactors && EPSLE(sum, 1.0, 1e-9) && expcurvneg )
444 else
446 curv = SCIPexprcurvMultiply(mult, curv);
447
448 return curv;
449}
450
451/** for a monomial with given bounds for each factor, gives condition on the curvature of each factor,
452 * so that monomial has a requested curvature, if possible
453 *
454 * @return whether `monomialcurv` can be achieved
455 */
457 SCIP_EXPRCURV monomialcurv, /**< desired curvature */
458 int nfactors, /**< number of factors in monomial */
459 SCIP_Real* exponents, /**< exponents in monomial, or NULL if all 1.0 */
460 SCIP_INTERVAL* factorbounds, /**< bounds of each factor */
461 SCIP_EXPRCURV* factorcurv /**< buffer to store required curvature of each factor */
462 )
463{
464 int nnegative;
465 int npositive;
466 SCIP_INTERVAL bounds;
467 SCIP_Real e;
468 SCIP_Real sum;
469 int j;
470
472 assert(nfactors >= 1);
475
476 if( nfactors == 1 )
477 {
480 }
481
482 /* any decent monomial with at least 2 factors is not linear */
484 return FALSE;
485
486 /* count positive and negative exponents, sum of exponents; flip negative factors */
487 nnegative = 0; /* number of negative exponents */
488 npositive = 0; /* number of positive exponents */
489 sum = 0.0; /* sum of exponents */
490 for( j = 0; j < nfactors; ++j )
491 {
492 e = exponents != NULL ? exponents[j] : 1.0;
493 assert(e != 0.0); /* should have been simplified away */
494
495 bounds = factorbounds[j];
496
497 /* if argument is negative, then exponent should be integer
498 * if that didn't happen, consider argument as if non-negative
499 */
500 if( !EPSISINT(e, 0.0) && bounds.inf < 0.0 ) /*lint !e835*/
501 {
502 bounds.inf = 0.0;
503 if( bounds.sup < 0.0 )
504 return FALSE;
505 }
506
507 /* mixed signs are bad */
508 if( bounds.inf < 0.0 && bounds.sup > 0.0 )
509 return FALSE;
510
511 if( e < 0.0 )
512 ++nnegative;
513 else
514 ++npositive;
515 sum += e;
516
517 if( bounds.inf < 0.0 )
518 {
519 /* flip j'th argument: (f_j)^(exp_j) = (-1)^(exp_j) (-f_j)^(exp_j)
520 * thus, negate monomial, if exponent is odd, i.e., (-1)^(exp_j) = -1
521 */
522 if( (int)e % 2 != 0 )
524 }
525 }
526
527 /* if all factors are linear, then a product f_j^exp_j with f_j >= 0 is convex if
528 * - all exponents are negative, or
529 * - all except one exponent j* are negative and exp_j* >= 1 - sum_{j!=j*}exp_j, but the latter is equivalent to sum_j exp_j >= 1
530 * further, the product is concave if
531 * - all exponents are positive and the sum of exponents is <= 1.0
532 *
533 * if factors are nonlinear, then we require additionally, that for convexity
534 * - each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0
535 * and for concavity, we require that
536 * - all factors are concave, i.e., exp_j*f_j'' <= 0
537 */
538
540 {
541 if( nnegative < nfactors-1 ) /* at least two positive exponents */
542 return FALSE;
543 if( nnegative < nfactors && !EPSGE(sum, 1.0, 1e-9) ) /* one negative exponent, but sum is not >= 1 */
544 return FALSE;
545
546 /* monomial will be convex, if each factor is convex if exp_j >= 0, or concave if exp_j <= 0, i.e., exp_j*f_j'' >= 0 */
547 for( j = 0; j < nfactors; ++j )
548 {
549 e = exponents != NULL ? exponents[j] : 1.0;
550
551 /* if factor is negative, then factorcurv[j] need to be flipped, which we can also get by flipping e */
552 if( factorbounds[j].inf < 0.0 && EPSISINT(e, 0.0) ) /*lint !e835*/
553 e = -e;
554 if( e >= 0.0 )
556 else
558 }
559 }
560 else
561 {
563 if( npositive < nfactors ) /* at least one negative exponent */
564 return FALSE;
565 if( !EPSLE(sum, 1.0, 1e-9) ) /* sum is not <= 1 */
566 return FALSE;
567
568 /* monomial will be concave, if each factor is concave */
569 for( j = 0; j < nfactors; ++j )
570 {
571 e = exponents != NULL ? exponents[j] : 1.0;
572
573 /* if factor is negative, then factorcurv[j] need to be flipped, i.e. convex */
574 if( factorbounds[j].inf < 0.0 && EPSISINT(e, 0.0) ) /*lint !e835*/
576 else
578 }
579 }
580
581 return TRUE;
582}
583
584/** gives name as string for a curvature */
586 SCIP_EXPRCURV curv /**< curvature */
587 )
588{
589 assert(0 <= curv && curv <= SCIP_EXPRCURV_LINEAR); /*lint !e685 !e2650 !e587 !e831 !e641 !e568*/
590
591 return curvnames[curv];
592}
#define EPSGE(x, y, eps)
Definition def.h:215
#define EPSISINT(x, eps)
Definition def.h:223
#define EPSLE(x, y, eps)
Definition def.h:213
#define TRUE
Definition def.h:95
#define FALSE
Definition def.h:96
#define SCIPABORT()
Definition def.h:360
static const char * curvnames[4]
Definition exprcurv.c:37
const char * SCIPexprcurvGetName(SCIP_EXPRCURV curv)
Definition exprcurv.c:585
SCIP_EXPRCURV SCIPexprcurvMonomial(int nfactors, SCIP_Real *exponents, int *factoridxs, SCIP_EXPRCURV *factorcurv, SCIP_INTERVAL *factorbounds)
Definition exprcurv.c:330
SCIP_EXPRCURV SCIPexprcurvPower(SCIP_INTERVAL basebounds, SCIP_EXPRCURV basecurv, SCIP_Real exponent)
Definition exprcurv.c:100
SCIP_EXPRCURV SCIPexprcurvPowerInv(SCIP_INTERVAL basebounds, SCIP_Real exponent, SCIP_EXPRCURV powercurv)
Definition exprcurv.c:208
SCIP_Bool SCIPexprcurvMonomialInv(SCIP_EXPRCURV monomialcurv, int nfactors, SCIP_Real *exponents, SCIP_INTERVAL *factorbounds, SCIP_EXPRCURV *factorcurv)
Definition exprcurv.c:456
SCIP_EXPRCURV SCIPexprcurvMultiply(SCIP_Real factor, SCIP_EXPRCURV curvature)
Definition exprcurv.c:87
SCIP_EXPRCURV SCIPexprcurvAdd(SCIP_EXPRCURV curv1, SCIP_EXPRCURV curv2)
Definition exprcurv.c:52
SCIP_EXPRCURV SCIPexprcurvNegate(SCIP_EXPRCURV curvature)
Definition exprcurv.c:61
void SCIPintervalSetBounds(SCIP_INTERVAL *resultant, SCIP_Real inf, SCIP_Real sup)
assert(minobj< SCIPgetCutoffbound(scip))
#define NULL
Definition lpi_spx1.cpp:161
public functions to work with algebraic expressions
#define SCIPerrorMessage
Definition pub_message.h:64
SCIP_Real sup
SCIP_Real inf
SCIP_EXPRCURV
Definition type_expr.h:58
@ SCIP_EXPRCURV_CONVEX
Definition type_expr.h:60
@ SCIP_EXPRCURV_LINEAR
Definition type_expr.h:62
@ SCIP_EXPRCURV_UNKNOWN
Definition type_expr.h:59
@ SCIP_EXPRCURV_CONCAVE
Definition type_expr.h:61