[Feature request] WebAssembly target

MuTsunTsai commented 4 months ago

I want to incorporate a MINLP solver for my web app, and it seems that SCIP is a good choice, but the architecture is a bit too complicated for me to figure out how to build it into WebAssembly with Emscripten. It will be wonderful if the SCIP team can add this target to the release!

MuTsunTsai commented 4 months ago

OK, I've got some initial success in compiling SCIP + Soplex into a static library for WASM, but I still need some help. I've written a basic example of my use case:

#include "scip/scip.h"
#include "scip/scipdefplugins.h"
#include "soplex/exceptions.h"

#include <iostream>

using namespace std;

struct PT {
    double x;
    double y;
};

struct POINT {
    SCIP_VAR *x;
    SCIP_VAR *y;
};

struct FLAP {
    POINT pt;
    double width;
    double height;
};

SCIP_RETCODE setLessThan(SCIP *scip, SCIP_VAR *left, SCIP_VAR *right, double by, const char *name) {
    SCIP_CONS *cons;
    SCIP_CALL(SCIPcreateConsBasicLinear(scip, &cons, name, 0, NULL, NULL, -SCIPinfinity(scip), -by));
    SCIP_CALL(SCIPaddCoefLinear(scip, cons, left, 1));
    SCIP_CALL(SCIPaddCoefLinear(scip, cons, right, -1));
    SCIP_CALL(SCIPaddCons(scip, cons));
    SCIP_CALL(SCIPreleaseCons(scip, &cons));
    return SCIP_OKAY;
}

SCIP_RETCODE setPointMinDist(SCIP *scip, POINT p1, POINT p2, double dist, const char *name) {
    SCIP_Real qcf[] = {1, -2, 1, 1, -2, 1};
    SCIP_VAR *qv1[] = {p1.x, p1.x, p2.x, p1.y, p1.y, p2.y};
    SCIP_VAR *qv2[] = {p1.x, p2.x, p2.x, p1.y, p2.y, p2.y};
    SCIP_CONS *cons;
    SCIP_CALL(SCIPcreateConsBasicQuadraticNonlinear(
        scip, &cons, name,
        0, NULL, NULL, 6, qv1, qv2, qcf,
        dist, SCIPinfinity(scip)));
    SCIP_CALL(SCIPaddCons(scip, cons));
    SCIP_CALL(SCIPreleaseCons(scip, &cons));
    return SCIP_OKAY;
}

/**
 * Create an expression expr = max(child, 0) by using (x + abs(x)) / 2.
 */
SCIP_RETCODE SCIPcreateExprMaxWithZero(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child) {
    SCIP_EXPR *abs;
    SCIP_Real cf[] = {0.5, 0.5};
    SCIP_CALL(SCIPcreateExprAbs(scip, &abs, child, NULL, NULL));
    SCIP_EXPR *exprs[] = {child, abs};
    SCIP_CALL(SCIPcreateExprSum(scip, expr, 2, exprs, cf, 0, NULL, NULL));
    SCIP_CALL(SCIPreleaseExpr(scip, &abs));
    return SCIP_OKAY;
}

SCIP_RETCODE SCIPcreateExprOrdinateDiff(SCIP *scip, SCIP_EXPR **expr, SCIP_VAR *v1, double u1, SCIP_VAR *v2, double u2) {
    SCIP_EXPR *e1, *e2, *abs, *sum1, *sum2;
    SCIP_CALL(SCIPcreateExprVar(scip, &e1, v1, NULL, NULL));
    SCIP_CALL(SCIPcreateExprVar(scip, &e2, v2, NULL, NULL));
    SCIP_EXPR *exprs[] = {e1, e2};
    SCIP_Real cf[] = {1, -1};
    SCIP_CALL(SCIPcreateExprSum(scip, &sum1, 2, exprs, cf, (u1 - u2) / 2, NULL, NULL));
    SCIP_CALL(SCIPcreateExprAbs(scip, &abs, sum1, NULL, NULL));
    SCIP_CALL(SCIPcreateExprSum(scip, &sum2, 1, &abs, cf, (-u1 - u2) / 2, NULL, NULL));
    SCIP_CALL(SCIPcreateExprMaxWithZero(scip, expr, sum2));
    SCIP_CALL(SCIPreleaseExpr(scip, &e1));
    SCIP_CALL(SCIPreleaseExpr(scip, &e2));
    SCIP_CALL(SCIPreleaseExpr(scip, &abs));
    SCIP_CALL(SCIPreleaseExpr(scip, &sum1));
    SCIP_CALL(SCIPreleaseExpr(scip, &sum2));
    return SCIP_OKAY;
}

SCIP_RETCODE setFlapMinDist(SCIP *scip, FLAP f1, FLAP f2, double dist, const char *name) {
    SCIP_CONS *cons;
    SCIP_EXPR *diff_x, *diff_y, *sqx, *sqy, *sum;
    SCIP_CALL(SCIPcreateExprOrdinateDiff(scip, &diff_x, f1.pt.x, f1.width, f2.pt.x, f2.width));
    SCIP_CALL(SCIPcreateExprPow(scip, &sqx, diff_x, 2, NULL, NULL));
    SCIP_CALL(SCIPcreateExprOrdinateDiff(scip, &diff_y, f1.pt.y, f1.height, f2.pt.y, f2.height));
    SCIP_CALL(SCIPcreateExprPow(scip, &sqy, diff_y, 2, NULL, NULL));
    SCIP_EXPR *exprs[] = {sqx, sqy};
    SCIP_Real cf[] = {1, 1};
    SCIP_CALL(SCIPcreateExprSum(scip, &sum, 2, exprs, cf, 0, NULL, NULL));
    SCIP_CALL(SCIPcreateConsBasicNonlinear(scip, &cons, name, sum, dist * dist, SCIPinfinity(scip)));
    SCIP_CALL(SCIPaddCons(scip, cons));
    SCIP_CALL(SCIPreleaseCons(scip, &cons));
    SCIP_CALL(SCIPreleaseExpr(scip, &diff_x));
    SCIP_CALL(SCIPreleaseExpr(scip, &diff_y));
    SCIP_CALL(SCIPreleaseExpr(scip, &sqx));
    SCIP_CALL(SCIPreleaseExpr(scip, &sqy));
    SCIP_CALL(SCIPreleaseExpr(scip, &sum));
    return SCIP_OKAY;
}

SCIP_RETCODE setupOffset(SCIP *scip, int n, PT *r, SCIP_VAR *s, POINT *p, SCIP_VAR *o) {
    SCIP_CONS *cons;
    SCIP_EXPR *quad, *prod, *es, *eo, *sum;
    size_t size = n * 6;

    SCIP_Real *qcf;
    SCIP_VAR **qv1;
    SCIP_VAR **qv2;
    SCIP_CALL(SCIPallocBufferArray(scip, &qcf, size));
    SCIP_CALL(SCIPallocBufferArray(scip, &qv1, size));
    SCIP_CALL(SCIPallocBufferArray(scip, &qv2, size));
    for (int i = 0; i < n; i++) {
        int base = i * 6;
        double values_qcf[] = {1, -2 * r[i].x, r[i].x * r[i].x, 1, -2 * r[i].y, r[i].y * r[i].y};
        SCIP_VAR *values_qv1[] = {p[i].x, p[i].x, s, p[i].y, p[i].y, s};
        SCIP_VAR *values_qv2[] = {p[i].x, s, s, p[i].y, s, s};
        for (int j = 0; j < 6; j++) {
            qcf[base + j] = values_qcf[j];
            qv1[base + j] = values_qv1[j];
            qv2[base + j] = values_qv2[j];
        }
    }
    SCIP_CALL(SCIPcreateExprQuadratic(scip, &quad, 0, NULL, NULL, size, qv1, qv2, qcf, NULL, NULL));

    SCIP_CALL(SCIPcreateExprVar(scip, &es, s, NULL, NULL));
    SCIP_CALL(SCIPcreateExprVar(scip, &eo, o, NULL, NULL));
    SCIP_EXPR *prod_exprs[] = {eo, es, es};
    SCIP_CALL(SCIPcreateExprProduct(scip, &prod, 3, prod_exprs, -n, NULL, NULL));

    SCIP_EXPR *sum_exprs[] = {quad, prod};
    SCIP_Real scf[] = {1, 1};
    SCIP_CALL(SCIPcreateExprSum(scip, &sum, 2, sum_exprs, scf, 0, NULL, NULL));

    SCIP_CALL(SCIPcreateConsBasicNonlinear(scip, &cons, "offset", sum, 0, 0));
    SCIP_CALL(SCIPaddCons(scip, cons));
    SCIP_CALL(SCIPreleaseCons(scip, &cons));

    SCIPfreeBufferArray(scip, &qcf);
    SCIPfreeBufferArray(scip, &qv1);
    SCIPfreeBufferArray(scip, &qv2);
    SCIP_CALL(SCIPreleaseExpr(scip, &quad));
    SCIP_CALL(SCIPreleaseExpr(scip, &prod));
    SCIP_CALL(SCIPreleaseExpr(scip, &es));
    SCIP_CALL(SCIPreleaseExpr(scip, &eo));
    SCIP_CALL(SCIPreleaseExpr(scip, &sum));
    return SCIP_OKAY;
}

int main() {
    SCIP *scip;
    SCIP_CALL(SCIPcreate(&scip));
    SCIPprintVersion(scip, NULL);

    SCIP_CALL(SCIPsetIntParam(scip, "display/freq", 1));

    SCIP_CALL(SCIPincludeDefaultPlugins(scip));
    SCIP_CALL(SCIPcreateProbBasic(scip, "Test"));

    const int n = 2;

    // Create variables
    SCIP_VAR *s;
    SCIP_VAR *o;
    POINT *p;
    SCIP_CALL(SCIPallocMemoryArray(scip, &p, (size_t)n));
    char name[SCIP_MAXSTRLEN];

    double infinity = SCIPinfinity(scip);
    SCIP_CALL(SCIPcreateVarBasic(scip, &s, "s", 1, infinity, 1, SCIP_VARTYPE_INTEGER));
    SCIP_CALL(SCIPaddVar(scip, s));
    SCIP_CALL(SCIPcreateVarBasic(scip, &o, "o", 0, 1, 0.5, SCIP_VARTYPE_CONTINUOUS));
    SCIP_CALL(SCIPaddVar(scip, o));
    for (int i = 0; i < n; i++) {
        SCIPsnprintf(name, SCIP_MAXSTRLEN, "x%d", i);
        SCIP_CALL(SCIPcreateVarBasic(scip, &p[i].x, name, 0, infinity, 0, SCIP_VARTYPE_INTEGER));
        SCIP_CALL(SCIPaddVar(scip, p[i].x));
        SCIPsnprintf(name, SCIP_MAXSTRLEN, "y%d", i);
        SCIP_CALL(SCIPcreateVarBasic(scip, &p[i].y, name, 0, infinity, 0, SCIP_VARTYPE_INTEGER));
        SCIP_CALL(SCIPaddVar(scip, p[i].y));
    }

    // Setup constraints
    FLAP f[] = {
        {p[0], 1, 0},
        {p[1], 1, 0}};
    for (int i = 0; i < n; i++) {
        SCIPsnprintf(name, SCIP_MAXSTRLEN, "bound_x%d", i);
        SCIP_CALL(setLessThan(scip, f[i].pt.x, s, f[i].width, name));
        SCIPsnprintf(name, SCIP_MAXSTRLEN, "bound_y%d", i);
        SCIP_CALL(setLessThan(scip, f[i].pt.y, s, f[i].height, name));
    }
    PT r[2] = {{1 / 3.0, 1 / 3.0}, {2 / 3.0, 2 / 3.0}};
    SCIP_CALL(setupOffset(scip, n, r, s, p, o));

    SCIP_CALL(setFlapMinDist(scip, f[0], f[1], 7, "dist_1_2"));

    // Print problem
    SCIPinfoMessage(scip, NULL, "\nOriginal problem:\n");
    SCIP_CALL(SCIPprintOrigProblem(scip, NULL, "cip", FALSE));
    SCIPinfoMessage(scip, NULL, "\n\n");

    // Solve and print solution
    try {
        SCIP_CALL(SCIPsolve(scip));
    } catch (soplex::SPxException ex) {
        std::cout << ex.what() << '\n';
        return SCIP_ERROR;
    }
    SCIP_CALL(SCIPprintBestSol(scip, NULL, TRUE));

    // Release memory
    SCIP_CALL(SCIPreleaseVar(scip, &s));
    SCIP_CALL(SCIPreleaseVar(scip, &o));
    for (int i = 0; i < n; i++) {
        SCIP_CALL(SCIPreleaseVar(scip, &p[i].x));
        SCIP_CALL(SCIPreleaseVar(scip, &p[i].y));
    }
    SCIPfreeMemoryArray(scip, &p);
    SCIP_CALL(SCIPfree(&scip));

    return SCIP_OKAY;
}

in which, long story short, it essentially generate this problem:

STATISTICS
  Problem name     : Test
  Variables        : 6 (0 binary, 5 integer, 0 implicit integer, 1 continuous)
  Constraints      : 0 initial, 6 maximal
OBJECTIVE
  Sense            : minimize
VARIABLES
  [integer] <s>: obj=1, original bounds=[1,+inf]
  [integer] <x0>: obj=0, original bounds=[0,+inf]
  [integer] <y0>: obj=0, original bounds=[0,+inf]
  [integer] <x1>: obj=0, original bounds=[0,+inf]
  [integer] <y1>: obj=0, original bounds=[0,+inf]
  [continuous] <o>: obj=0.5, original bounds=[0,1]
CONSTRAINTS
  [linear] <bound_x0>: <x0>[I] -<s>[I] <= -1;
  [linear] <bound_y0>: <y0>[I] -<s>[I] <= -0;
  [linear] <bound_x1>: <x1>[I] -<s>[I] <= -1;
  [linear] <bound_y1>: <y1>[I] -<s>[I] <= -0;
  [nonlinear] <offset>: ((<x0>)^2-0.666667*<x0>*<s>+0.111111*(<s>)^2+(<y0>)^2-0.666667*<y0>*<s>+0.111111*(<s>)^2+(<x1>)^2-1.33333*<x1>*<s>+0.444444*(<s>)^2+(<y1>)^2-1.33333*<y1>*<s>+0.444444*(<s>)^2)+(-2)*<o>*<s>*<s> == 0;
  [nonlinear] <dist_1_2>: ((0.5*(-1+abs((<x0>-<x1>)))+0.5*abs((-1+abs((<x0>-<x1>))))))^2+((0.5*(abs((<y0>-<y1>)))+0.5*abs((abs((<y0>-<y1>))))))^2 >= 49;
END

Now the weird thing is that running the compiled WASM will result in SPxInternalCodeException("XENTER05 This should never happen.") (found here) as it runs into node 126 during the SCIPsolve(scip) process. However, running the same problem instance with native build SCIP + Soplex will solve the problem correctly.

Also, this issue happens only with some combination of the numbers. For example, if I change {2 / 3.0, 2 / 3.0} in my code to {1.0, 1.0}, then the WASM build will also solve the problem without error.

I can really use some help in figuring this out. I was speculating that the different behavior between WASM build and native build could result from alignment issues etc., but as I add the debugging flags to the WASM building process, it detects no such code.

This is my make.wasm32.web.emcc.opt file:

CC          =   emcc
CXX         =   em++
LINKCC      =   emcc
LINKCXX     =   em++
DCC         =   emcc
DCXX        =   em++
OFLAGS      +=  -O3 -fomit-frame-pointer
LDFLAGS     +=  -lm -m32
FLAGS       +=  -DNDEBUG -DBOOST_NO_CXX98_FUNCTION_BASE -DSCIP_ROUNDING_WASM
CFLAGS      =   -m32 -ffp-contract=off -std=c99 -D_XOPEN_SOURCE=600
CXXFLAGS    =   -m32 -std=c++14 -ffp-contract=off
ARFLAGS     =   crs

(Regarding the flag SCIP_ROUNDING_WASM, because WASM doesn't support changing rounding mode, I'll have to add something like this to intervalarith.c to make it work

#ifdef SCIP_ROUNDING_WASM
#define ROUNDING

#define SCIP_ROUND_DOWNWARDS 0
#define SCIP_ROUND_UPWARDS   0
#define SCIP_ROUND_NEAREST   0
#define SCIP_ROUND_ZERO      0

SCIP_Bool SCIPintervalHasRoundingControl() {
   return TRUE;
}

static void intervalSetRoundingMode(SCIP_ROUNDMODE roundmode) {
}

static SCIP_ROUNDMODE intervalGetRoundingMode() {
   return 0;
}

#endif

but I think this is not what causing the issue, since the error is in Soplex and not in SCIP. Besides, I've also tried using this code for the native build, and still the native build runs without error.)

and then I run the compiling with

make scipoptlib OSTYPE=wasm32 ARCH=web COMP=emcc AMPL=false BOOST=false ZLIB=false ZIMPL=false GMP=false READLINE=false EXPRINT=none

Any insights as for why the XENTER05 error might happen will be greatly appreciated!

svigerske commented 4 months ago

Oh dear, maybe should have answered earlier.

The handling of nonlinear constraints in SCIP assumes a properly working rounding mode control, otherwise the domain propagation could wrongly cut off nodes.
It is strongly advised to have a NLP solver available to solve MINLPs, at least when there are continuous nonlinear variables. The three options are closed-source Worhp, difficult-to-obtain FilterSQP (in Fortran), and Ipopt. Ipopt is C++, but requires a linear solver and Blas/Lapack to work, which easily brings in Fortran code. The last I knew, flang wasn't ready yet, so I'm not so sure about the Fortran support of Emscripten.
If you can make a NLP solver available, then also an AD package will be necessary. Currently, the only supported by SCIP is CppAD, so EXPRINT=none should be removed.

Regarding the SoPlex issue, I don't know, but you may want to try solving MIPs first. Missing rounding mode control should play no or only a little role there, so it is more likely that the WebAssembly code runs similar to native code.

MuTsunTsai commented 4 months ago

@svigerske Wow that was really quick and helpful reply! Big thanks. I'm an absolute beginner about this MINLP thing, so please bear with me.

The handling of nonlinear constraints in SCIP assumes a properly working rounding mode control

So I guess that I'm in some sense out of luck here, since WASM simply doesn't have rounding mode control (source, and it probably won't have it in foreseeable future from what I've heard)? Sounds like there's no way I can port SCIP to WASM in general, but perhaps only for my particular use case in which I can use the workaround you suggested?

I'm not so sure about the Fortran support of Emscripten

According to this, people had successfully incorporate Blas/Lapack into their WASM projects, but sounds like a long way to go for me (I mean, I'm so burned out just to get to the current progress!). So if I focus on MIP first, does that mean I can live without an NLP for the moment?

EXPRINT=none should be removed

For now I'm having that just to simplify the matter as much as I can in order to look into the Soplex issue (but maybe I don't have to, hopefully).

Now to give more context to my use case, the only constraint involving <o> (which is the only continuous variable) is <offset>, and its sole purpose is just to select the MIP solution that has the closet relative distribution to an initial guess by the user. If I understand you correctly, you're suggesting that I can just remove <o> and <offset>, find all optimal solutions to the MIP (there's probably not going to be a lot of them in general), and then just pick among them the one with the closet distribution myself? Sounds like a plan to me!

MuTsunTsai commented 4 months ago

@svigerske I think I'm stuck again. I tried to use the approach described here, so I first let it solve essentially this problem:

STATISTICS
  Problem name     : Test
  Variables        : 5 (0 binary, 5 integer, 0 implicit integer, 0 continuous)
  Constraints      : 0 initial, 5 maximal
OBJECTIVE
  Sense            : minimize
VARIABLES
  [integer] <s>: obj=1, original bounds=[1,+inf]
  [integer] <x0>: obj=0, original bounds=[0,+inf]
  [integer] <y0>: obj=0, original bounds=[0,+inf]
  [integer] <x1>: obj=0, original bounds=[0,+inf]
  [integer] <y1>: obj=0, original bounds=[0,+inf]
CONSTRAINTS
  [linear] <bound_x0>: <x0>[I] -<s>[I] <= -1;
  [linear] <bound_y0>: <y0>[I] -<s>[I] <= -0;
  [linear] <bound_x1>: <x1>[I] -<s>[I] <= -1;
  [linear] <bound_y1>: <y1>[I] -<s>[I] <= -0;
  [nonlinear] <dist_1_2>: ((0.5*(-1+abs((<x0>-<x1>)))+0.5*abs((-1+abs((<x0>-<x1>))))))^2+((0.5*(abs((<y0>-<y1>)))+0.5*abs((abs((<y0>-<y1>))))))^2 >= 49;
END

And it correctly find the optimal value for <s> is 6. Next, I changed the bounds of <s> to [6,6] and execute solution counting (I'm expecting four solutions by symmetry), but then it reports that the problem is infeasible with zero solutions. This is the case both for WASM and native builds. What did I do wrong here?

(I got a warning saying a non-trivial solution comes in over <SCIP_DECL_CONSCHECK(consCheckCountsols)>; currently these solutions are ignored.. Maybe it's relevant but I don't know what to do with it.)

scipopt / scip

[Feature request] WebAssembly target #96