Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
mingodadampl opened this issue Jun 15, 2020 · 2 comments
Open

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

mingodadampl opened this issue Jun 15, 2020 · 2 comments

Comments

@mingodadampl
Copy link

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
Copy link
Author

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
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant