From 27542a6d591d891378bcb184c7da6e6019163d0a Mon Sep 17 00:00:00 2001 From: Brendan O'Donoghue Date: Tue, 4 Apr 2023 19:37:49 +0100 Subject: [PATCH] New exp cone projection (#247) * New exponential cone projection, about 50x faster than the original one. * Based on "Projection onto the exponential cone: a univariate root-finding problem", Friberg, 2021. * Add test covering exp cone projection. --- CMakeLists.txt | 1 + Makefile | 3 +- src/cones.c | 172 +--------- src/exp_cone.c | 399 +++++++++++++++++++++++ test/problems/max_ent.h | 2 +- test/problems/random_prob.h | 2 +- test/problems/test_exp_cone.h | 84 +++++ test/problems/test_prob_from_data_file.h | 4 +- test/run_tests.c | 2 + 9 files changed, 506 insertions(+), 163 deletions(-) create mode 100644 src/exp_cone.c create mode 100644 test/problems/test_exp_cone.h diff --git a/CMakeLists.txt b/CMakeLists.txt index e40a8e58..26d6be17 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -161,6 +161,7 @@ set(${PROJECT_NAME}_SRC src/aa.c src/cones.c src/ctrlc.c + src/exp_cone.c src/linalg.c src/normalize.c src/rw.c diff --git a/Makefile b/Makefile index 7576b07a..06f49efb 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ # MAKEFILE for scs include scs.mk -SCS_OBJECTS = src/util.o src/cones.o src/aa.o src/rw.o src/linalg.o src/ctrlc.o src/scs_version.o src/normalize.o +SCS_OBJECTS = src/util.o src/cones.o src/exp_cone.o src/aa.o src/rw.o src/linalg.o src/ctrlc.o src/scs_version.o src/normalize.o SCS_O = src/scs.o SCS_INDIR_O = src/scs_indir.o @@ -41,6 +41,7 @@ $(SCS_INDIR_O): src/scs.c $(INC_FILES) src/util.o : src/util.c $(INC_FILES) src/cones.o : src/cones.c $(INC_FILES) +src/exp_cone.o : src/exp_cone.c $(INC_FILES) src/aa.o : src/aa.c $(INC_FILES) src/rw.o : src/rw.c $(INC_FILES) src/linalg.o: src/linalg.c $(INC_FILES) diff --git a/src/cones.c b/src/cones.c index cf09d0ab..42eb03fc 100644 --- a/src/cones.c +++ b/src/cones.c @@ -4,10 +4,8 @@ #include "scs_blas.h" /* contains BLAS(X) macros and type info */ #include "util.h" -#define CONE_TOL (1e-9) -#define CONE_THRESH (1e-8) -#define EXP_CONE_MAX_ITERS (100) #define BOX_CONE_MAX_ITERS (25) +#define POW_CONE_TOL (1e-9) #define POW_CONE_MAX_ITERS (20) /* Box cone limits (+ or -) taken to be INF */ @@ -35,6 +33,11 @@ void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx, #endif +/* Forward declare exponential cone projection routine. + * Implemented in exp_cone.c. + */ +scs_float SCS(proj_pd_exp_cone)(scs_float *v0, scs_int primal); + void SCS(free_cone)(ScsCone *k) { if (k) { if (k->bu) @@ -325,118 +328,6 @@ char *SCS(get_cone_header)(const ScsCone *k) { return tmp; } -static scs_float exp_newton_one_d(scs_float rho, scs_float y_hat, - scs_float z_hat, scs_float w) { - scs_float t_prev, t = MAX(w - z_hat, MAX(-z_hat, 1e-9)); - scs_float f = 1., fp = 1.; - scs_int i; - for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) { - t_prev = t; - f = t * (t + z_hat) / rho / rho - y_hat / rho + log(t / rho) + 1; - fp = (2 * t + z_hat) / rho / rho + 1 / t; - - t = t - f / fp; - - if (t <= -z_hat) { - t = -z_hat; - break; - } else if (t <= 0) { - t = 0; - break; - } else if (ABS(t - t_prev) < CONE_TOL) { - break; - } else if (SQRTF(f * f / fp) < CONE_TOL) { - break; - } - } - if (i == EXP_CONE_MAX_ITERS) { - scs_printf("warning: exp cone newton step hit maximum %i iters\n", (int)i); - scs_printf("rho=%1.5e; y_hat=%1.5e; z_hat=%1.5e; w=%1.5e; f=%1.5e, " - "fp=%1.5e, t=%1.5e, t_prev= %1.5e\n", - rho, y_hat, z_hat, w, f, fp, t, t_prev); - } - return t + z_hat; -} - -static void exp_solve_for_x_with_rho(const scs_float *v, scs_float *x, - scs_float rho, scs_float w) { - x[2] = exp_newton_one_d(rho, v[1], v[2], w); - x[1] = (x[2] - v[2]) * x[2] / rho; - x[0] = v[0] - rho; -} - -static scs_float exp_calc_grad(const scs_float *v, scs_float *x, scs_float rho, - scs_float w) { - exp_solve_for_x_with_rho(v, x, rho, w); - if (x[1] <= 1e-12) { - return x[0]; - } - return x[0] + x[1] * log(x[1] / x[2]); -} - -static void exp_get_rho_ub(const scs_float *v, scs_float *x, scs_float *ub, - scs_float *lb) { - *lb = 0; - *ub = 0.125; - while (exp_calc_grad(v, x, *ub, v[1]) > 0) { - *lb = *ub; - (*ub) *= 2; - } -} - -/* project onto the exponential cone, v has dimension *exactly* 3 */ -static scs_int proj_exp_cone(scs_float *v) { - scs_int i; - scs_float ub, lb, rho, g, x[3]; - scs_float r = v[0], s = v[1], t = v[2]; - - /* v in cl(Kexp) */ - if ((s * exp(r / s) - t <= CONE_THRESH && s > 0) || - (r <= 0 && s == 0 && t >= 0)) { - return 0; - } - - /* -v in Kexp^* */ - if ((r > 0 && r * exp(s / r) + exp(1) * t <= CONE_THRESH) || - (r == 0 && s <= 0 && t <= 0)) { - memset(v, 0, 3 * sizeof(scs_float)); - return 0; - } - - /* special case with analytical solution */ - if (r < 0 && s < 0) { - v[1] = 0.0; - v[2] = MAX(v[2], 0); - return 0; - } - - /* iterative procedure to find projection, bisects on dual variable: */ - exp_get_rho_ub(v, x, &ub, &lb); /* get starting upper and lower bounds */ - for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) { - rho = (ub + lb) / 2; /* halfway between upper and lower bounds */ - g = exp_calc_grad(v, x, rho, x[1]); /* calculates gradient wrt dual var */ - if (g > 0) { - lb = rho; - } else { - ub = rho; - } - if (ub - lb < CONE_TOL) { - break; - } - } -#if VERBOSITY > 10 - scs_printf("exponential cone proj iters %i\n", (int)i); -#endif - if (i == EXP_CONE_MAX_ITERS) { - scs_printf("warning: exp cone outer step hit maximum %i iters\n", (int)i); - scs_printf("r=%1.5e; s=%1.5e; t=%1.5e\n", r, s, t); - } - v[0] = x[0]; - v[1] = x[1]; - v[2] = x[2]; - return 0; -} - static scs_int set_up_sd_cone_work_space(ScsConeWork *c, const ScsCone *k) { scs_int i; #ifdef USE_LAPACK @@ -741,13 +632,13 @@ static void proj_power_cone(scs_float *v, scs_float a) { scs_int i; /* v in K_a */ if (xh >= 0 && yh >= 0 && - CONE_THRESH + POWF(xh, a) * POWF(yh, (1 - a)) >= rh) { + POW_CONE_TOL + POWF(xh, a) * POWF(yh, (1 - a)) >= rh) { return; } /* -v in K_a^* */ if (xh <= 0 && yh <= 0 && - CONE_THRESH + POWF(-xh, a) * POWF(-yh, 1 - a) >= + POW_CONE_TOL + POWF(-xh, a) * POWF(-yh, 1 - a) >= rh * POWF(a, a) * POWF(1 - a, 1 - a)) { v[0] = v[1] = v[2] = 0; return; @@ -760,7 +651,7 @@ static void proj_power_cone(scs_float *v, scs_float a) { y = pow_calc_x(r, yh, rh, 1 - a); f = pow_calc_f(x, y, r, a); - if (ABS(f) < CONE_TOL) { + if (ABS(f) < POW_CONE_TOL) { break; } @@ -829,51 +720,16 @@ static scs_int proj_cone(scs_float *x, const ScsCone *k, ScsConeWork *c, } } - if (k->ep) { /* doesn't use r_y */ - /* - * exponential cone is not self dual, if s \in K - * then y \in K^* and so if K is the primal cone - * here we project onto K^*, via Moreau - * \Pi_C^*(y) = y + \Pi_C(-y) - */ + if (k->ep || k->ed) { /* doesn't use r_y */ #ifdef _OPENMP #pragma omp parallel for #endif - for (i = 0; i < k->ep; ++i) { - proj_exp_cone(&(x[count + 3 * i])); + for (i = 0; i < k->ep + k->ed; ++i) { + /* provided in exp_cone.c */ + SCS(proj_pd_exp_cone)(&(x[count + 3 * i]), i < k->ep); } - count += 3 * k->ep; + count += 3 * (k->ep + k->ed); } - - /* dual exponential cone */ - if (k->ed) { /* doesn't use r_y */ - /* - * exponential cone is not self dual, if s \in K - * then y \in K^* and so if K is the primal cone - * here we project onto K^*, via Moreau - * \Pi_C^*(y) = y + \Pi_C(-y) - */ - scs_int idx; - scs_float r, s, t; - SCS(scale_array)(&(x[count]), -1, 3 * k->ed); /* x = -x; */ -#ifdef _OPENMP -#pragma omp parallel for private(r, s, t, idx) -#endif - for (i = 0; i < k->ed; ++i) { - idx = count + 3 * i; - r = x[idx]; - s = x[idx + 1]; - t = x[idx + 2]; - - proj_exp_cone(&(x[idx])); - - x[idx] -= r; - x[idx + 1] -= s; - x[idx + 2] -= t; - } - count += 3 * k->ed; - } - if (k->psize && k->p) { /* doesn't use r_y */ scs_float v[3]; scs_int idx; diff --git a/src/exp_cone.c b/src/exp_cone.c new file mode 100644 index 00000000..96c5f081 --- /dev/null +++ b/src/exp_cone.c @@ -0,0 +1,399 @@ +#include "cones.h" +#include "glbopts.h" +#include "linalg.h" +#include "scs.h" + +#define EXP_CONE_INFINITY_VALUE (1E15) + +/* + * Exponential cone projection routines, from: + * + * Projection onto the exponential cone: a univariate root-finding problem, + * by Henrik A. Friberg, 2021. + * + */ + +static inline scs_int _isfinite(scs_float x) { + return ABS(x) < EXP_CONE_INFINITY_VALUE; +} + +static inline scs_float _clip(scs_float x, scs_float l, scs_float u) { + return MAX(l, MIN(u, x)); +} + +/* As defined in Friberg, 2021 (multiplied by positive polynomial) */ +static void hfun(const scs_float *v0, scs_float rho, scs_float *f, + scs_float *df) { + scs_float t0 = v0[2], s0 = v0[1], r0 = v0[0]; + scs_float exprho = exp(rho); + scs_float expnegrho = exp(-rho); + /* function value at v0 */ + *f = ((rho - 1) * r0 + s0) * exprho - (r0 - rho * s0) * expnegrho - + (rho * (rho - 1) + 1) * t0; + /* gradient of function at v0 */ + *df = (rho * r0 + s0) * exprho + (r0 - (rho - 1) * s0) * expnegrho - + (2 * rho - 1) * t0; +} + +/* Binary search for the root of the hfun function */ +static scs_float root_search_binary(const scs_float *v0, scs_float xl, + scs_float xu, scs_float x) { +#if VERBOSITY > 0 + scs_printf("Exp cone: Newton method failed, resorting to binary search.\n"); +#endif + const scs_float EPS = 1e-12; /* expensive so loosen tol */ + const scs_int MAXITER = 40; + scs_int i; + scs_float x_plus = x, f, df; + for (i = 0; i < MAXITER; i++) { + hfun(v0, x, &f, &df); + if (f < 0.0) { + xl = x; + } else { + xu = x; + } + /* binary search step */ + x_plus = 0.5 * (xl + xu); + if (ABS(x_plus - x) <= EPS * MAX(1., ABS(x_plus)) || (x_plus == xl) || + (x_plus == xu)) { + break; + } + x = x_plus; + } + return x_plus; +} + +/* Use damped Newton's to find the root of the hfun function */ +static scs_float root_search_newton(const scs_float *v0, scs_float xl, + scs_float xu, scs_float x) { + /* params taken from Friberg code */ + const scs_float EPS = 1e-15; + const scs_float DFTOL = 1e-13; /* pow(EPS, 6.0 / 7.0) */ + const scs_int MAXITER = 20; + const scs_float LODAMP = 0.05; + const scs_float HIDAMP = 0.95; + + scs_float x_plus, f, df; + scs_int i; + + for (i = 0; i < MAXITER; i++) { + hfun(v0, x, &f, &df); + + if (ABS(f) <= EPS) { /* root found */ + break; + } + + if (f < 0.0) { + xl = x; + } else { + xu = x; + } + + if (xu <= xl) { + xu = 0.5 * (xu + xl); + xl = xu; + break; + } + + if (!_isfinite(f) || df < DFTOL) { + break; + } + + /* Newton step */ + x_plus = x - f / df; + + if (ABS(x_plus - x) <= EPS * MAX(1., ABS(x_plus))) { + break; + } + + if (x_plus >= xu) { + x = MIN(LODAMP * x + HIDAMP * xu, xu); + } else if (x_plus <= xl) { + x = MAX(LODAMP * x + HIDAMP * xl, xl); + } else { + x = x_plus; + } + } + if (i < MAXITER) { /* Newton's method converged */ +#if VERBOSITY > 0 + scs_printf("Exp cone: Newton iters:%i, f:%.4e, df:%.4e\n", (int)i, f, df); +#endif + return _clip(x, xl, xu); + } + /* Fall back to binary search if Newton failed */ + return root_search_binary(v0, xl, xu, x); +} + +/* try heuristic (cheap) projection */ +static scs_float proj_primal_exp_cone_heuristic(const scs_float *v0, + scs_float *vp) { + scs_float t0 = v0[2], s0 = v0[1], r0 = v0[0]; + scs_float dist, tp, newdist; + /* perspective boundary */ + vp[2] = MAX(t0, 0); + vp[1] = 0.0; + vp[0] = MIN(r0, 0); + dist = SCS(norm_diff)(v0, vp, 3); + + /* perspective interior */ + if (s0 > 0.0) { + tp = MAX(t0, s0 * exp(r0 / s0)); + newdist = tp - t0; + if (newdist < dist) { + vp[2] = tp; + vp[1] = s0; + vp[0] = r0; + dist = newdist; + } + } + return dist; +} + +/* try heuristic (cheap) projection */ +static scs_float proj_polar_exp_cone_heuristic(const scs_float *v0, + scs_float *vd) { + scs_float t0 = v0[2], s0 = v0[1], r0 = v0[0]; + scs_float dist, td, newdist; + /* perspective boundary */ + vd[2] = MIN(t0, 0); + vd[1] = MIN(s0, 0); + vd[0] = 0.0; + dist = SCS(norm_diff)(v0, vd, 3); + + /* perspective interior */ + if (r0 > 0.0) { + td = MIN(t0, -r0 * exp(s0 / r0 - 1)); + newdist = t0 - td; + if (newdist < dist) { + vd[2] = td; + vd[1] = s0; + vd[0] = r0; + dist = newdist; + } + } + return dist; +} + +static scs_float ppsi(const scs_float *v0) { + scs_float s0 = v0[1], r0 = v0[0]; + scs_float psi; + + if (r0 > s0) { + psi = (r0 - s0 + sqrt(r0 * r0 + s0 * s0 - r0 * s0)) / r0; + } else { + psi = -s0 / (r0 - s0 - sqrt(r0 * r0 + s0 * s0 - r0 * s0)); + } + + return ((psi - 1) * r0 + s0) / (psi * (psi - 1) + 1); +} + +static scs_float pomega(scs_float rho) { + scs_float val = exp(rho) / (rho * (rho - 1) + 1); + + if (rho < 2.0) { + val = MIN(val, exp(2.0) / 3); + } + + return val; +} + +static scs_float dpsi(const scs_float *v0) { + scs_float s0 = v0[1], r0 = v0[0]; + scs_float psi; + + if (s0 > r0) { + psi = (r0 - sqrt(r0 * r0 + s0 * s0 - r0 * s0)) / s0; + } else { + psi = (r0 - s0) / (r0 + sqrt(r0 * r0 + s0 * s0 - r0 * s0)); + } + + return (r0 - psi * s0) / (psi * (psi - 1) + 1); +} + +static scs_float domega(scs_float rho) { + scs_float val = -exp(-rho) / (rho * (rho - 1) + 1); + + if (rho > -1.0) { + val = MAX(val, -exp(1.0) / 3); + } + + return val; +} + +/* Generate upper and lower search bounds for root of hfun */ +static void exp_search_bracket(const scs_float *v0, scs_float pdist, + scs_float ddist, scs_float *low_out, + scs_float *upr_out) { + scs_float t0 = v0[2], s0 = v0[1], r0 = v0[0]; + scs_float baselow = -EXP_CONE_INFINITY_VALUE, + baseupr = EXP_CONE_INFINITY_VALUE; + scs_float low = -EXP_CONE_INFINITY_VALUE, upr = EXP_CONE_INFINITY_VALUE; + + scs_float Dp = SQRTF(pdist * pdist - MIN(s0, 0) * MIN(s0, 0)); + scs_float Dd = SQRTF(ddist * ddist - MIN(r0, 0) * MIN(r0, 0)); + + scs_float curbnd, fl, fu, df, tpu, tdl; + + if (t0 > 0) { + curbnd = log(t0 / ppsi(v0)); + low = MAX(low, curbnd); + } else if (t0 < 0) { + curbnd = -log(-t0 / dpsi(v0)); + upr = MIN(upr, curbnd); + } + + if (r0 > 0) { + baselow = 1 - s0 / r0; + low = MAX(low, baselow); + tpu = MAX(1e-12, MIN(Dd, Dp + t0)); + curbnd = MAX(low, baselow + tpu / r0 / pomega(low)); + upr = MIN(upr, curbnd); + } + + if (s0 > 0) { + baseupr = r0 / s0; + upr = MIN(upr, baseupr); + tdl = -MAX(1e-12, MIN(Dp, Dd - t0)); + curbnd = MIN(upr, baseupr - tdl / s0 / domega(upr)); + low = MAX(low, curbnd); + } + + /* Guarantee valid bracket */ + /* TODO do we need these 2 lines? */ + low = _clip(MIN(low, upr), baselow, baseupr); + upr = _clip(MAX(low, upr), baselow, baseupr); + + if (low != upr) { + hfun(v0, low, &fl, &df); + hfun(v0, upr, &fu, &df); + + if (fl * fu > 0) { + if (ABS(fl) < ABS(fu)) { + upr = low; + } else { + low = upr; + } + } + } + + *low_out = low; + *upr_out = upr; +} + +/* convert from rho to primal projection */ +static scs_float proj_sol_primal_exp_cone(const scs_float *v0, scs_float rho, + scs_float *vp) { + scs_float linrho = (rho - 1) * v0[0] + v0[1]; + scs_float exprho = exp(rho); + scs_float quadrho, dist; + if (linrho > 0 && _isfinite(exprho)) { + quadrho = rho * (rho - 1) + 1; + vp[2] = exprho * linrho / quadrho; + vp[1] = linrho / quadrho; + vp[0] = rho * linrho / quadrho; + dist = SCS(norm_diff)(vp, v0, 3); + } else { + vp[2] = EXP_CONE_INFINITY_VALUE; + vp[1] = 0.0; + vp[0] = 0.0; + dist = EXP_CONE_INFINITY_VALUE; + } + return dist; +} + +/* convert from rho to polar projection */ +static scs_float proj_sol_polar_exp_cone(const scs_float *v0, scs_float rho, + scs_float *vd) { + scs_float linrho = v0[0] - rho * v0[1]; + scs_float exprho = exp(-rho); + scs_float quadrho, dist; + if (linrho > 0 && _isfinite(exprho)) { + quadrho = rho * (rho - 1) + 1; + vd[2] = -exprho * linrho / quadrho; + vd[1] = (1 - rho) * linrho / quadrho; + vd[0] = linrho / quadrho; + dist = SCS(norm_diff)(v0, vd, 3); + } else { + vd[2] = -EXP_CONE_INFINITY_VALUE; + vd[1] = 0.0; + vd[0] = 0.0; + dist = EXP_CONE_INFINITY_VALUE; + } + return dist; +} + +static inline void _copy(scs_float *dest, scs_float *src) { + dest[0] = src[0]; + dest[1] = src[1]; + dest[2] = src[2]; +} + +/* Project onto primal or dual exponential cone, performed in-place. + * If `primal=0` then project on the dual cone, otherwise project + * onto primal. Taken from algorithm in Friberg, 2021. + */ +scs_float SCS(proj_pd_exp_cone)(scs_float *v0, scs_int primal) { + scs_float TOL = 1e-8; /* pow(1e-10, 2.0 / 3.0); */ + scs_float xl, xh, pdist, ddist, err, rho, dist_hat; + scs_float vp[3], vd[3], v_hat[3]; + scs_int opt; + if (!primal) { + /* This routine actually projects onto primal and polar cones + * simultaneously. So to make it project onto dual, use this: + * Pi_{C^*}(v0) = -Pi_{C^polar}(-v0) + */ + v0[0] *= -1.; + v0[1] *= -1.; + v0[2] *= -1.; + } + + pdist = proj_primal_exp_cone_heuristic(v0, vp); + ddist = proj_polar_exp_cone_heuristic(v0, vd); + + err = ABS(vp[0] + vd[0] - v0[0]); + err = MAX(err, ABS(vp[1] + vd[1] - v0[1])); + err = MAX(err, ABS(vp[2] + vd[2] - v0[2])); + + /* Skip root search if presolve rules apply + * or optimality conditions are satisfied + */ + opt = (v0[1] <= 0 && v0[0] <= 0); + opt |= (MIN(pdist, ddist) <= TOL); + opt |= (err <= TOL && SCS(dot)(vp, vd, 3) <= TOL); + if (opt) { + if (primal) { + _copy(v0, vp); + return pdist; + } + /* polar cone -> dual cone */ + v0[0] = -vd[0]; + v0[1] = -vd[1]; + v0[2] = -vd[2]; + return ddist; + } + + exp_search_bracket(v0, pdist, ddist, &xl, &xh); + rho = root_search_newton(v0, xl, xh, 0.5 * (xl + xh)); + + if (primal) { + /* primal cone projection */ + dist_hat = proj_sol_primal_exp_cone(v0, rho, v_hat); + if (dist_hat <= pdist) { + _copy(vp, v_hat); + pdist = dist_hat; + } + _copy(v0, vp); + return pdist; + } + /* polar cone projection */ + dist_hat = proj_sol_polar_exp_cone(v0, rho, v_hat); + if (dist_hat <= ddist) { + _copy(vd, v_hat); + ddist = dist_hat; + } + /* polar cone -> dual cone */ + v0[0] = -vd[0]; + v0[1] = -vd[1]; + v0[2] = -vd[2]; + return ddist; +} diff --git a/test/problems/max_ent.h b/test/problems/max_ent.h index 7a9b4b69..408cffea 100644 --- a/test/problems/max_ent.h +++ b/test/problems/max_ent.h @@ -1,6 +1,6 @@ #include "glbopts.h" -#include "scs.h" #include "problems/test_prob_from_data_file.h" +#include "scs.h" static const char *max_ent(void) { scs_float OPT = -6.067087663361563; /* from ecos */ diff --git a/test/problems/random_prob.h b/test/problems/random_prob.h index 91f67747..6d5544ba 100644 --- a/test/problems/random_prob.h +++ b/test/problems/random_prob.h @@ -1,6 +1,6 @@ #include "glbopts.h" -#include "scs.h" #include "problems/test_prob_from_data_file.h" +#include "scs.h" static const char *random_prob(void) { scs_float OPT = 5.751458006385587; diff --git a/test/problems/test_exp_cone.h b/test/problems/test_exp_cone.h new file mode 100644 index 00000000..c86dda47 --- /dev/null +++ b/test/problems/test_exp_cone.h @@ -0,0 +1,84 @@ +#include "glbopts.h" +#include "linalg.h" +#include "minunit.h" +#include "scs.h" + +/* Forward declare exponential cone projection routine. + * Implemented in exp_cone.c. + */ +scs_float SCS(proj_pd_exp_cone)(scs_float *v0, scs_int primal); + +static scs_int _run_exp_cone_test(scs_float *v0, scs_float *vp_true, + scs_float *vd_true) { + scs_int success = 1; + scs_float vp[3], vd[3]; + const scs_float TOL = 1e-6; + + memcpy(vp, v0, 3 * sizeof(scs_float)); + memcpy(vd, v0, 3 * sizeof(scs_float)); + + /* inefficient, but just for testing */ + SCS(proj_pd_exp_cone)(vp, 1); + SCS(proj_pd_exp_cone)(vd, 0); + + scs_printf("*******************************\n"); + scs_printf("v0: (%f, %f, %f)\n", v0[0], v0[1], v0[2]); + scs_printf("vp: (%f, %f, %f)\n", vp[0], vp[1], vp[2]); + scs_printf("vp_true: (%f, %f, %f)\n", vp_true[0], vp_true[1], vp_true[2]); + scs_printf("vd: (%f, %f, %f)\n", vd[0], vd[1], vd[2]); + scs_printf("vd_true: (%f, %f, %f)\n", vd_true[0], vd_true[1], vd_true[2]); + + success &= (SCS(norm_diff)(vp, vp_true, 3) <= TOL); + success &= (SCS(norm_diff)(vd, vd_true, 3) <= TOL); + /* Moreau decomposition holds only for polar */ + /* + success &= (SCS(dot)(vp, vd, 3) <= TOL); + success &= (ABS(v0[0] - vp[0] + vd[0]) <= TOL); + success &= (ABS(v0[1] - vp[1] + vd[1]) <= TOL); + success &= (ABS(v0[2] - vp[2] + vd[2]) <= TOL); + */ + + if (!success) { + scs_printf("Failed.\n"); + } + + return success; +} + +static const char *test_exp_cone(void) { + scs_int success = 1; + scs_int i; + /* test points */ + scs_float v0[6][3] = { + {1, 2, 3}, + {0.14814832, 1.04294573, 0.67905585}, + {-0.78301134, 1.82790084, -1.05417044}, + {1.3282585, -0.43277314, 1.7468072}, + {0.67905585, 0.14814832, 1.04294573}, + {0.50210027, 0.12314491, -1.77568921}, + }; + /* primal projections */ + scs_float vp[6][3] = { + {0.8899428, 1.94041881, 3.06957226}, + {-0.02001571, 0.8709169, 0.85112944}, + {-1.17415616, 0.9567094, 0.280399}, + {0.53160512, 0.2804836, 1.86652094}, + {0.38322814, 0.27086569, 1.11482228}, + {0., 0., 0.}, + }; + /* dual projections */ + scs_float vd[6][3] = { + {-0., 2., 3.}, + {-0., 1.04294573, 0.67905585}, + {-0.68541419, 1.85424082, 0.01685653}, + {-0.02277033, -0.12164823, 1.75085347}, + {-0., 0.14814832, 1.04294573}, + {-0., 0.12314491, -0.}, + }; + + for (i = 0; i < 6; ++i) { + success &= _run_exp_cone_test(v0[i], vp[i], vd[i]); + } + mu_assert("test_exp_cone: Failure", success); + return 0; +} diff --git a/test/problems/test_prob_from_data_file.h b/test/problems/test_prob_from_data_file.h index 655d9275..5e972bdb 100644 --- a/test/problems/test_prob_from_data_file.h +++ b/test/problems/test_prob_from_data_file.h @@ -1,4 +1,4 @@ -#ifndef _SCS_FILE_TEST_CHASSIS +#ifndef _SCS_FILE_TEST_CHASSIS #define _SCS_FILE_TEST_CHASSIS #include "glbopts.h" @@ -8,7 +8,7 @@ #include "scs.h" #include "util.h" -static const char *_test_prob_from_data(const char * file, scs_float OPT) { +static const char *_test_prob_from_data(const char *file, scs_float OPT) { scs_int read_status; ScsData *d; ScsCone *k; diff --git a/test/run_tests.c b/test/run_tests.c index b61c28bd..337b5d7a 100644 --- a/test/run_tests.c +++ b/test/run_tests.c @@ -12,6 +12,7 @@ #include "problems/qafiro_tiny_qp.h" #include "problems/small_lp.h" #include "problems/small_qp.h" +#include "problems/test_exp_cone.h" #include "problems/unbounded_tiny_qp.h" int tests_run = 0; @@ -60,6 +61,7 @@ static const char *all_tests(void) { mu_run_test(unbounded_tiny_qp); mu_run_test(random_prob); mu_run_test(max_ent); + mu_run_test(test_exp_cone); return 0; } int main(void) {