search

ampl / asl

BSD 3-Clause "New" or "Revised" License
18 stars 14 forks source link

SCIP driver do not set ".sstatus" on solution file #3

Open mingodadampl opened 4 years ago

mingodadampl commented 4 years ago

While testing the snapshot command with scipampl I noticed that no ".sstatus" is returned on the solution file and looking through the source code I found this functions in "scip/src/lpi/lpi.h":

/*
 * LP Basis Methods
 */

/**@name LP Basis Methods */
/**@{ */

/** gets current basis status for columns and rows; arrays must be large enough to store the basis status */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBase(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int*                  cstat,              /**< array to store column basis status, or NULL */
   int*                  rstat               /**< array to store row basis status, or NULL */
   );

/** sets current basis status for columns and rows */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiSetBase(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   const int*            cstat,              /**< array with column basis status */
   const int*            rstat               /**< array with row basis status */
   );

/** returns the indices of the basic columns and rows; basic column n gives value n, basic row m gives value -1-m */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBasisInd(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int*                  bind                /**< pointer to store basis indices ready to keep number of rows entries */
   );

/** get row of inverse basis matrix B^-1
 *
 *  @note The LP interface defines slack variables to have coefficient +1. This means that if, internally, the LP solver
 *        uses a -1 coefficient, then rows associated with slacks variables whose coefficient is -1, should be negated;
 *        see also the explanation in lpi.h.
 */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBInvRow(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int                   r,                  /**< row number */
   SCIP_Real*            coef,               /**< pointer to store the coefficients of the row */
   int*                  inds,               /**< array to store the non-zero indices, or NULL */
   int*                  ninds               /**< pointer to store the number of non-zero indices, or NULL
                                              *   (-1: if we do not store sparsity information) */
   );

/** get column of inverse basis matrix B^-1
 *
 *  @note The LP interface defines slack variables to have coefficient +1. This means that if, internally, the LP solver
 *        uses a -1 coefficient, then rows associated with slacks variables whose coefficient is -1, should be negated;
 *        see also the explanation in lpi.h.
 */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBInvCol(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int                   c,                  /**< column number of B^-1; this is NOT the number of the column in the LP;
                                              *   you have to call SCIPlpiGetBasisInd() to get the array which links the
                                              *   B^-1 column numbers to the row and column numbers of the LP!
                                              *   c must be between 0 and nrows-1, since the basis has the size
                                              *   nrows * nrows */
   SCIP_Real*            coef,               /**< pointer to store the coefficients of the column */
   int*                  inds,               /**< array to store the non-zero indices, or NULL */
   int*                  ninds               /**< pointer to store the number of non-zero indices, or NULL
                                              *   (-1: if we do not store sparsity information) */
   );

/** get row of inverse basis matrix times constraint matrix B^-1 * A
 *
 *  @note The LP interface defines slack variables to have coefficient +1. This means that if, internally, the LP solver
 *        uses a -1 coefficient, then rows associated with slacks variables whose coefficient is -1, should be negated;
 *        see also the explanation in lpi.h.
 */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBInvARow(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int                   r,                  /**< row number */
   const SCIP_Real*      binvrow,            /**< row in (A_B)^-1 from prior call to SCIPlpiGetBInvRow(), or NULL */
   SCIP_Real*            coef,               /**< vector to return coefficients of the row */
   int*                  inds,               /**< array to store the non-zero indices, or NULL */
   int*                  ninds               /**< pointer to store the number of non-zero indices, or NULL
                                              *   (-1: if we do not store sparsity information) */
   );

/** get column of inverse basis matrix times constraint matrix B^-1 * A
 *
 *  @note The LP interface defines slack variables to have coefficient +1. This means that if, internally, the LP solver
 *        uses a -1 coefficient, then rows associated with slacks variables whose coefficient is -1, should be negated;
 *        see also the explanation in lpi.h.
 */
SCIP_EXPORT
SCIP_RETCODE SCIPlpiGetBInvACol(
   SCIP_LPI*             lpi,                /**< LP interface structure */
   int                   c,                  /**< column number */
   SCIP_Real*            coef,               /**< vector to return coefficients of the column */
   int*                  inds,               /**< array to store the non-zero indices, or NULL */
   int*                  ninds               /**< pointer to store the number of non-zero indices, or NULL
                                              *   (-1: if we do not store sparsity information) */
   );

/**@} */

And looking at other drivers it's clear that who wrote the actual scipampl interface didn't bothered to implement it.

mingodadampl commented 4 years ago

After add some code to scipampl driver:

--- /home/mingo/dev/lp/scipoptsuite-6.0.1/scip/interfaces/ampl/src/reader_nl.c
+++ /home/mingo/dev/lp/scipoptsuite-7.0.0/scip/interfaces/ampl/src/reader_nl.c
@@ -1496,6 +1496,8 @@
 }
 #endif

+static SufDecl suftable[17];
+
 /** problem reading method of reader */
 static
 SCIP_DECL_READERREAD(readerReadNl)
@@ -1505,7 +1507,6 @@
    const char* filebasename;
    FILE* nl;
    ASL* asl;
-   SufDecl suftable[15];

    assert(scip != NULL);
    assert(reader != NULL);
@@ -1581,7 +1582,14 @@
    suftable[14].kind = ASL_Sufkind_var;
    suftable[14].nextra = 0;

-   suf_declare(suftable, 15);
+   suftable[15].name = (char*)"sstatus";
+   suftable[15].kind = ASL_Sufkind_var;
+   suftable[15].nextra = 0;
+
+   suftable[16].name = (char*)"sstatus";
+   suftable[16].kind = ASL_Sufkind_con;
+   suftable[16].nextra = 0;
+

    /* let ASL read .nl file
     * jac0dim will do exit(1) if the file is not found
@@ -1593,6 +1601,8 @@
       SCIPerrorMessage("error processing <%s> file\n", filename);
       return SCIP_READERROR;
    }
+
+   suf_declare(suftable, 17);

    want_derivs = 0;
    want_xpi0 = 1;
@@ -1896,6 +1906,22 @@
          assert(y[c] != SCIP_INVALID); /*lint !e777*/
       }
    }
+   
+   if (solve_result_num < 500 && scip->lp && scip->lp->lpi) {
+        int *consstat = (int*)M1alloc((n_con + n_var)*sizeof(int));
+        int *varstat = consstat + n_con;
+        SufDesc *csd = suf_iput("sstatus", ASL_Sufkind_con, consstat);
+        SufDesc *vsd = suf_iput("sstatus", ASL_Sufkind_var, varstat );
+        SCIP_CALL( SCIPlpiGetBase(scip->lp->lpi, varstat, consstat) );
+        for(int i=0; i < n_con; ++i) {
+            if(!consstat[i]) consstat[i] = 3;
+            //printf("consstat[%d] = %d %s\n", i, consstat[i], con_name(i));
+        }
+        for(int i=0; i < n_var; ++i) {
+            if(!varstat[i]) varstat[i] = 3;
+            //printf("varstat[%d] = %d %s\n", i, varstat[i], var_name(i));
+        }
+   }

    write_sol((char*)msg, x, y, NULL);

And test with this ampl script:

option precision 0;
model diet.mod;
data diet.dat;
problem DietProb: Buy, Total_Cost, Diet;
problem DietProb;

option auxfiles rc;
write gmain;

option solver scipampl;
#option solver cplex;
#option solver gurobi;
#option solver xpress;
option solver cbc;

solve;
#display Diet;
#snapshot > snapshot-main0.ampl;

display Buy.astatus, Buy.defeqn, Buy.dual, Buy.init, Buy.init0, Buy.lb, Buy.ub, Buy.lb0, Buy.ub0, Buy.lb1, Buy.ub1, Buy.lb2, Buy.ub2, Buy.lrc, Buy.urc, Buy.lslack, Buy.uslack, Buy.rc, Buy.slack, Buy.sstatus, Buy.status, Buy.val;

display Diet.astatus, Diet.body, Diet.defvar, Diet.dinit, Diet.dinit0, Diet.dual, Diet.lb, Diet.ub, Diet.lbs, Diet.ubs, Diet.ldual, Diet.udual, Diet.lslack, Diet.uslack, Diet.slack, Diet.sstatus, Diet.status;

display Total_Cost.astatus, Total_Cost.exitcode, Total_Cost.message, Total_Cost.result, Total_Cost.val;

display DietProb.astatus, DietProb.exitcode, DietProb.message, DietProb.result;

And testing with different solvers I noticed that doesn't seem to have any harmonization of the returning sstatus codes, and the gurobi solver/driver do not send any sstatus for constraints:

MINOS 5.51: optimal solution found.
6 iterations, objective 88.2

Options
3
0
1
0
4
4
8
8
0.0018181818181818108
0.008181818181818278
0.11599999999999991
-8.169327588411304e-18
0
0
0
0
46.66666666666666
1.5761812194954111e-15
8.429823983987501e-15
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 3
1 3
2 3
3 3
4 1
5 1
6 1
7 3
suffix 1 4 8 0 0
sstatus
0 3
1 3
2 3
3 1
CPLEX 12.10.0.0: optimal solution; objective 88.2
1 dual simplex iterations (0 in phase I)

Options
3
0
1
0
4
4
8
8
0
0
0.126
0
0
0
0
0
46.666666666666664
0
0
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 3
1 3
2 3
3 3
4 1
5 3
6 3
7 3
suffix 1 4 8 0 0
sstatus
0 1
1 1
2 3
3 1
Gurobi 9.0.2: optimal solution; objective 88.2
1 simplex iterations

Options
3
0
1
0
4
4
8
8
0
0
0.126
0
0
0
0
0
46.666666666666664
0
0
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 3
1 3
2 3
3 3
4 1
5 3
6 3
7 3
suffix 1 0 8 0 0
sstatus
XPRESS 8.8.0(35.01.01): Optimal solution found
Objective 88.2
1 simplex iteration

Options
3
0
1
0
4
4
8
8
0
0
0.126
0
0
0
0
0
46.666666666666664
0
0
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 3
1 3
2 3
3 3
4 1
5 3
6 3
7 3
suffix 1 4 8 0 0
sstatus
0 1
1 1
2 4
3 1
CBC 2.10.5 optimal, objective 88.2
1 iterations

Options
3
0
1
0
4
4
8
8
0
0
0.126
0
0
0
0
0
46.66666666666667
0
0
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 1
1 1
2 1
3 1
4 3
5 1
6 1
7 1
suffix 1 4 8 0 0
sstatus
0 3
1 3
2 4
3 3

scipampl:

optimal solution found

Options
3
0
1
0
4
0
8
8
0
0
0
0
46.666666666666664
0
0
0
objno 0 0
suffix 0 8 8 0 0
sstatus
0 3
1 3
2 3
3 3
4 1
5 3
6 3
7 3
suffix 1 4 8 0 0
sstatus
0 1
1 1
2 3
3 1
mingodadampl commented 4 years ago

Also note that scipampl do not send the value 0.126 and several calculated fields like reduced costs do not seem to be correct (showing Buy.sstatus). scipampl:

BEEF        3.19  low
CHK     2.59  low
FISH        2.29  low
HAM     2.80  low
MCH     1.80  bas
MTL     1.99  low
SPG     1.99  low
TUR     2.49  low

minos:

BEEF        1.25909         low
CHK     0.0918182       low
FISH        0.992727        low
HAM     1.37091         low
MCH     2.85926e-16     bas
MTL     6.89169e-16     bas
SPG     -4.79712e-16    bas
TUR     1.09818         low

cplex, gurobi, xpress, cbc:

BEEF        1.3
CHK     0.07  
FISH        1.03  
HAM     1.63 
MCH     -2.22045e-16
MTL     0.1  
SPG     0.1  
TUR     1.23

cplex, gurobi and xpress show Buy.sstatus as shown in scipampl, but cbc just inverts then:

BEEF   bas
CHK    bas
FISH   bas 
HAM    bas
MCH    low
MTL    bas
SPG    bas 
TUR    bas