Python DSP

[MiftariB]

Hi,

I am looking to pass a structured MILP problem to DSP via python ctypes. However, I don't fully understand the workflow that needs to be implemented. I tried to add a single block to the DspApiEnv but I get a SEGFAULT each time.

My python code :

from ctypes import *
import numpy as np

libDSP = cdll.LoadLibrary("/DSP/build/src/libDSP.so")

pointer_to_model = c_void_p(libDSP.createEnv())
print(type(pointer_to_model))

row = [0, 1, 2, 3]
row = (c_double * len(row))(*row)

nrow = 4

col = [0, 1, 2, 3]
col = (c_int * len(col))(*col)

ncol = 4

val = [1, 2, 3, 4]
val = (c_double * len(val))(*val)

nval = 4

col_lb = [0, 0, 0, 0]
col_lb = (c_double* len(col_lb))(*col_lb)

col_up = [10, 10, 10, 10]
col_up = (c_double* len(col_up))(*col_up)

number = c_int(0)

coltypes = [73,73,73,73] # int corresponding to "I"
coltypes = (c_char* len(coltypes))(*coltypes)

libDSP.setNumberOfScenarios(pointer_to_model,1)

libDSP.setDimensions(pointer_to_model, 4, 4, 4, 4)

libDSP.loadBlockProblem(pointer_to_model, 0, ncol, nrow, nval, byref(number), 
						col, val, col_lb, col_up, coltypes, row, row, row)

Could you help me understand the different functions I need to call to pass a block structured MILP to DSP ?

Thanks in advance for your help.

Best Regards,

Miftari B

[kibaekkim]

In case you don't know yet, we have a modeling interface in Julia: https://github.com/kibaekkim/DSPopt.jl With the package, you can simply model problems in algebraic form.

---

If you want/need to use Python, you have to use the C-library functions available in DSP, as you are already trying.

I think the problem data is not correctly defined in your example. The problem data should be compatible to the function arguments defined here:

DSP/src/DspCInterface.cpp

Lines 292 to 306 in 6675d30

	void loadBlockProblem(
	DspApiEnv * env, /**< pointer to API object */
	int id, /**< block ID */
	int ncols, /**< number of columns */
	int nrows, /**< number of rows */
	int numels, /**< number of elements in the matrix */
	const CoinBigIndex * start, /**< start index for each row */
	const int * index, /**< column indices */
	const double * value, /**< constraint elements */
	const double * clbd, /**< column lower bounds */
	const double * cubd, /**< column upper bounds */
	const char * ctype, /**< column types */
	const double * obj, /**< objective coefficients */
	const double * rlbd, /**< row lower bounds */
	const double * rubd /**< row upper bounds */) {

In particular, these three arguments should represent the constraint matrix in row-wise sparse format:

DSP/src/DspCInterface.cpp

Lines 298 to 300 in 6675d30

	const CoinBigIndex * start, /**< start index for each row */
	const int * index, /**< column indices */
	const double * value, /**< constraint elements */

What you are trying to do has been coded in Julia like this:
https://github.com/kibaekkim/DSPopt.jl/blob/0812121ea8031c3f67a8cf93a8f59ecdfe0ba897/src/Core.jl#L138-L183
with the row-wise sparse matrix from
https://github.com/kibaekkim/DSPopt.jl/blob/0812121ea8031c3f67a8cf93a8f59ecdfe0ba897/src/Model.jl#L172-L226

[MiftariB]

First, Thanks for you time and your quick response. All my models are written in Python so I would rather implement my own interface with DSP directly rather than going through Julia.

I see my mistake in the argument CoinBigIndex * start, I should have passed an int type row there. However, I still get a segfault when I do so. If possible, could you explain with a toy c++ model how you would instanciate a block MILP directly from handwritten vectors and solve it ?

Thanks for your help again,

MIFTARI B

[kibaekkim]

@MiftariB I am working on a c++ example, but this will take some time.

[MiftariB]

Dear @kibaekkim,

Thanks for the answer and sorry for my previous enquiry, it came too early in my using of DSP. However, I made some progress in the python front. This is the current state:

from ctypes import *
import numpy as np
import pickle 
from scipy.sparse import coo_matrix, csr_matrix


libDSP = cdll.LoadLibrary("/DSP/build/src/libDSP.so")

#
# LOAD THE 5 PICKLE FILES
#
row, col, data = pickle.load( open( "constraint_matrix.pkl", "rb" ) ) # full constraint matrix
index_constraints = pickle.load( open( "index_constraints.pkl", "rb" ) ) # tuples  [start, end] of block's constraint index
index_variables = pickle.load( open( "index_var.pkl", "rb" ) ) # tuples  [start, end] of block's variables
b = pickle.load(open("b.pkl", "rb")) # right-hand side of constraints
C = pickle.load(open("C.pkl", "rb")) # objective vector


# 
# CONVERT THEIR TYPE AND INDEX FROM 0
#
row = np.array(row, dtype = int)-1
col = np.array(col, dtype = int)-1
data = np.array(data)
b = np.array(b)
C = np.array(C)
nb_lines = max(row)+1
nb_col = max(col)+1
index_constraints = np.array(index_constraints, dtype=int)-1
index_var = np.array(index_variables, dtype=int)-1


#
# CREATE THE FULL CONSTRAINT MATRIX 
#
A = coo_matrix((data, (row, col)), shape=(nb_lines, nb_col)).todense()


#
# CREATE THE BLOCK 0 -> MASTER BLOCK
#

# first its constraint matrix
start_master_constraint, end_master_constraint = index_constraints[-1]
master_row = []
master_col = []
master_data = []
for i in range(start_master_constraint, end_master_constraint+1):
	for j in range(nb_col):
		if A[i,j] != 0:
			master_row.append(i-start_master_constraint)
			master_col.append(j)
			master_data.append(A[i,j])

master_nb_lines = end_master_constraint-start_master_constraint+1
master_row = np.array(master_row, dtype=int)
master_col = np.array(master_col, dtype=int)
master_data = np.array(master_data)
master_csr = csr_matrix(coo_matrix((master_data,(master_row, master_col)), 
	shape=(master_nb_lines, nb_col)))

master_val, master_row, master_col = master_csr.data, master_csr.indptr, master_csr.indices

# arguments for block adding
master_val = (c_double*len(master_val))(*master_val)
master_row = (c_int*len(master_row))(*master_row)
master_col = (c_int*len(master_col))(*master_col)
master_C = (c_double*len(C))(*C)

minus_infinity_col = [-float("inf")]*nb_col
minus_infinity_col = (c_double*len(minus_infinity_col))(*minus_infinity_col)

plus_infinity_col =  [float("inf")]*nb_col
plus_infinity_col = (c_double*len(plus_infinity_col))(*plus_infinity_col)

master_coltypes = [67]*nb_col
master_coltypes = (c_char* len(master_coltypes))(*master_coltypes)

row_lb = [-float("inf")]*master_nb_lines
row_lb = (c_double*len(row_lb))(*row_lb)

row_up = b[start_master_constraint:end_master_constraint+1]
row_up = (c_double*len(row_up))(*row_up)

pointer_to_model = c_void_p(libDSP.createEnv())

libDSP.loadBlockProblem(pointer_to_model, 0, c_int(nb_col), c_int(master_nb_lines), len(master_val), 
	master_row, master_col, master_val, minus_infinity_col, plus_infinity_col, master_coltypes, master_C, row_lb, row_up)

print("Master loaded ")

#
# ADDING THE BLOCKS
#
objective_nul = [0]*len(C)
objective_nul = (c_double*len(objective_nul))(*objective_nul)
nb_blocks, _ = index_var.shape

for b_number in range(nb_blocks):
	block_row = []
	block_col = []
	block_data = []

	block_start_constr, block_end_constr = index_constraints[b_number]
	block_start_var, block_end_var = index_var[b_number]

	for i in range(block_start_constr, block_end_constr+1):
		for j in range(block_start_var, block_end_var+1):
			if A[i,j] != 0:
				block_row.append(i-block_start_constr)
				block_col.append(j)
				block_data.append(A[i,j])

	block_nb_col = block_end_var-block_start_var+1
	block_nb_lines = block_end_constr-block_start_constr+1

	block_row = np.array(block_row, dtype=int)
	block_col = np.array(block_col, dtype=int)
	block_data = np.array(block_data)

	block_csr = csr_matrix(coo_matrix((block_data,(block_row, block_col)), 
		shape=(block_nb_lines, nb_col)))

	block_val, block_row, block_col = block_csr.data, block_csr.indptr, block_csr.indices

	block_val = (c_double*len(block_val))(*block_val)
	block_row = (c_int*len(block_row))(*block_row)
	block_col = (c_int*len(block_col))(*block_col)

	block_coltypes = [67]*block_nb_col
	block_coltypes = (c_char* len(block_coltypes))(*block_coltypes)

	row_up = b[block_start_constr:block_end_constr+1]
	row_up = (c_double*len(row_up))(*row_up)

	libDSP.loadBlockProblem(pointer_to_model, b_number+1, c_int(block_nb_col), c_int(block_nb_lines), len(block_val), 
		block_row, block_col, block_val, minus_infinity_col, plus_infinity_col, block_coltypes, objective_nul, row_lb, row_up)

	print("Added block "+str(b_number+1))

print("Loaded ")
libDSP.printModel(pointer_to_model)
#libDSP.writeMps(pointer_to_model, "salut.mps")

libDSP.solveDe(pointer_to_model)

print(libDSP.getStatus(pointer_to_model))

#libDSP.solveDe(pointer_to_model)
libDSP.freeEnv(byref(pointer_to_model))
print("finished")

Adding the different blocks seems to work. However, the model printing yields a segfault after printing the "parent/master". I added the corresponding log file python_print.log. I tried not printing and going straight to the solving but the extended form also segfaults (python_de.log), dantzig-wolf yields a more complicated python error that i am yet to figure completely out (python_dw.log). I am not sure to understand where these segfaults could come from. Is there something else that i am using wrongly ? Is the work-flow more or less the correct one now ?
I also join the pickle files pickle.zip.

Best regards,

MIFTARI B

[kibaekkim]

By quickly looking at your example, I see the example is using coo_matrix and as a dense. DSP takes the matrix in CSR format, so I believe you should use csr_matrix (not as a dense). Please take a look at an work-in-progress example: https://github.com/Argonne-National-Laboratory/DSP/blob/12cdd1f6cf7c0f7e6fbc8dda3be5705fbe3dd7e3/examples/c/farmer.c

I also found one potential issue causing segfault from the print function. I will fix that shortly.

[MiftariB]

Hi,

I only use the COO form to match the pickle files. CSR matrices are sent to DSP corresponding to the lines :

master_csr = csr_matrix(coo_matrix((master_data,(master_row, master_col)), 
	shape=(master_nb_lines, nb_col)))

master_val, master_row, master_col = master_csr.data, master_csr.indptr, master_csr.indices

and the lines

block_csr = csr_matrix(coo_matrix((block_data,(block_row, block_col)), 
		shape=(block_nb_lines, nb_col)))

block_val, block_row, block_col = block_csr.data, block_csr.indptr, block_csr.indices

I will definitively have a look at your C++ example. Thanks.

Regards,

MIFTARI B

[kibaekkim]

Oh, I see. Thanks for the clarification. I must have overlooked that. I will try to run your example on my side.

[kibaekkim]

I am not familiar with python. But, this simple code is not running on my mac or linux, returning segfault.

from ctypes import *

libDSP = cdll.LoadLibrary("libDsp.dylib")
pointer_to_model = c_void_p(libDSP.createEnv())
# libDSP.freeEnv(byref(pointer_to_model)) # This does not run either.
libDSP.freeEnv(pointer_to_model)

But, this simple c++ code is running fine:

#include <stdlib.h>
#include <stdio.h>
#include <dlfcn.h>

int main(int argc, char **argv)
{
    void *handle;
    char *error;

    handle = dlopen ("libDsp.dylib", RTLD_LAZY);
    if (!handle) {
        fputs (dlerror(), stderr);
        exit(1);
    }

    void* env = NULL;
    env = ((void *(*)())dlsym(handle, "createEnv"))();
    if (env == NULL) {
        printf("Failed to create DSP environment.\n");
    } else {
        printf("Successfully created DSP environment.\n");
    }
   
    if (env != NULL) {
        ((void (*)(void*&))dlsym(handle, "freeEnv"))(env);
        env = NULL;
        printf("Successfully freed DSP environment.\n");
    }
    dlclose(handle);

    return 0;
}

@MiftariB Any idea from this?

[MiftariB]

The following code works perfectly on my end

from ctypes import *

libDSP = cdll.LoadLibrary("DSP/build/src/libDSP.so")
pointer_to_model = c_void_p(libDSP.createEnv())
print("hi")
libDSP.freeEnv(byref(pointer_to_model))
print("finished")

producing the following output test.log. I ran this code on Windows subsystem Linux Ubuntu 18.04 with python 3.7.5 and DSP installed with Cplex and Gurobi.

Can you maybe try the following :

from ctypes import *

libDSP = cdll.LoadLibrary("/DSP/build/src/libDSP.so")
 
libDSP.createEnv.restype = c_void_p
pointer_to_model = c_void_p(libDSP.createEnv())
# libDSP.freeEnv(byref(pointer_to_model)) # This does not run either.
print("hi")

libDSP.freeEnv.argtypes = [c_void_p]
libDSP.freeEnv(byref(pointer_to_model))
print("finished")

The problem might come from not knowing the return and argument types respectively. Can you tell me whether it's freeing or allocating the model that creates the segfault ? (I guess the freeing does.)

[kibaekkim]

Yes, giving explicit types fixed the python issue. I have also fixed some bugs in the DSP code. Can you try again with branch issue/205?

[MiftariB]

Thanks. I changed slightly the file by adding the arguments and return types :

#
# SETTING THE ARGSTYPE AND RETURN TYPES
#

libDSP.createEnv.restype = c_void_p

libDSP.loadBlockProblem.argtypes = [c_void_p, c_int, c_int, c_int, c_int, 
	POINTER(c_int), POINTER(c_int), POINTER(c_double), POINTER(c_double), 
	POINTER(c_double), POINTER(c_char), POINTER(c_double), POINTER(c_double), POINTER(c_double)]

libDSP.solveDe.argtypes = [c_void_p]

libDSP.solveDw.argtypes = [c_void_p]

libDSP.getStatus.argtypes = [c_void_p]

libDSP.getPrimalSolution.argtypes = [c_void_p, c_int, POINTER(c_double)]

libDSP.freeEnv.argtypes = [c_void_p]

libDSP.printModel.argtypes = [c_void_p]

libDSP.getPrimalBound.argtypes = [c_void_p]

libDSP.getPrimalBound.restype = c_double

I also pulled the branch (and rebuild DSP). The model printing and the extensive method do not segfault anymore. The result from DE is slightly different from the one i get with CPLEX (DE : 154.4 - CPLEX : 146.6000000003857). I still have to figure out why. Might come from the encoding. DW still produces the error :

python3: malloc.c:2401: sysmalloc: Assertion `(old_top == initial_top (av) && old_size == 0) || ((unsigned long) (old_size) >= MINSIZE && prev_inuse (old_top) && ((unsigned long) old_end & (pagesize - 1)) == 0)' failed.
Aborted

Here is also the full file, I used linking_test.zip.

[kibaekkim]

linking_test.py needs the modification:

# libDSP.loadBlockProblem(pointer_to_model, b_number+1, c_int(block_nb_col), c_int(block_nb_lines), len(block_val), 
# 	block_row, block_col, block_val, minus_infinity_col, plus_infinity_col, block_coltypes, objective_nul, row_lb, row_up)
libDSP.loadBlockProblem(pointer_to_model, b_number+1, c_int(nb_col), c_int(block_nb_lines), len(block_val), 
	block_row, block_col, block_val, minus_infinity_col, plus_infinity_col, master_coltypes, master_C, row_lb, row_up)

But, it does not seem to give the correct solution. I will look into that more.

[MiftariB]

Hi @kibaekkim,

I made a few changes. My lower bound in blocks was incorrect (not enough values if master block had less constraints than another block). One of the pickle files was incorrect too. I corrected both issues and now I have the correct objective : 146.6. Thanks for the help ! I still need to figure out what is going on with DW.
Here are all the files you will need to run the test : pythonDSP_test.zip

However, I don't fully understand why do we need the changes referenced in the previous comment : #205 (comment).
I will try it and let you know if I get the same results.

Regards,

MIFTARI B

[MiftariB]

Sorry for the double comment. I tried your modifications. I obtain the exact same output for the extensive form python_de.log. However, when trying with DW, I get the following output python_dw.log which ends with a segfault when trying to retrieve the solution (I guess) as it classifies the problem as unfeasible.

[kibaekkim]

However, I don't fully understand why do we need the changes referenced in the previous comment : #205 (comment).
I will try it and let you know if I get the same results.

The reason is based on how it should be used, for example:
https://github.com/kibaekkim/DSPopt.jl/blob/0812121ea8031c3f67a8cf93a8f59ecdfe0ba897/src/Core.jl#L168-L169

However, when trying with DW, I get the following output python_dw.log which ends with a segfault when trying to retrieve the solution (I guess) as it classifies the problem as unfeasible.

Yes, the subproblem is infeasible. See sub1.mps.zip.

[MiftariB]

There is something I am not getting. In the problem I am trying to solve, I put every variable as free. However in your mps file, they are fixed to 0. And of course, the constraint x49 <= -6.9 becomes unfeasible in that case.
Can you help me understand what am I missing here ?

[kibaekkim]

I found that these need to be added after adding blocks.

libDSP.updateBlocks.argtypes = [c_void_p]
libDSP.updateBlocks(pointer_to_model)

In future, I will see if I can prevent solveXX without calling this function.

Known bug

After adding the lines above, I ran again and found that one of your subproblems becomes unbounded and terminates the algorithm with error. This is a known bug I just found today (#211).

Workaround

For now you can avoid this by adding some bounds on your variables. For example, I was able to find the objective value of 154.4 with the variable bounds [-20, 20].

[MiftariB]

Thanks for the update. I will keep up with the issue 211 then. Thanks for your help ! I will be implementing a clean-up ( and more general) version of the python linking with DSP.

Thanks for your help,

Best regards,

MIFTARI B

[kibaekkim]

@MiftariB Awesome! I would really appreciate if you can share your python example with us in this repo.

[MiftariB]

Sure ! We will keep in touch ! Thanks again !

[MiftariB]

Dear @kibaekkim,

I am writing to follow up on our earlier exchange about a Python interface for DSP. To give some context, @MathiasPBerger and I are researchers at the University of Liège (Belgium) and we have been developing a new optimization modelling language for structured mixed-integer linear programming problems called GBOML. In short, GBOML makes it possible to exploit the underlying problem structure in order to simplify model encoding, enable model re-use, speed up model generation (e.g., by parallelizing matrix generation) and facilitate the application of structure-exploiting algorithms. The GBOML parser is written in Python and we decided to develop an interface to DSP to promote the use of structure-exploiting algorithms. DSP's Extensive Form and Dantzig-Wolfe algorithms have been tested via the Python interface so far, and we plan on testing the dual decomposition algorithm shortly. To learn more about GBOML, here is the documentation.

Our understanding is that you were interested in having a look at the Python interface, and we wanted to let you know that we would be happy to share the code with you. We plan on maintaining an interface to DSP as part of GBOML but the code could also be made stand-alone for other people to use. We would love to hear your thoughts on this.

Best regards,

Miftari B

[kibaekkim]

@MiftariB Thanks for the update. The package looks very interesting. I am interested in trying your package for some of our examples. Please feel free to let us know if there is anything we can do on the DSP side. I would also like to draw your attention on a Julia package that may be relevant to this direction of your research: https://github.com/zavalab/Plasmo.jl

Does your team also have a plan for modeling beyond linear? For example, quadratic, second-order cone, etc..? Given that one of your examples considers power system, it would be great to have some nonlinear supports.

[MiftariB]

Dear @kibaekkim,

Sorry with the late reply, we were a little overwhelmed by work lately. First, thanks for trying our package, any feedback we can get is very valuable to us. Do not hesitate to ask questions on our gitlab. Thanks for the link, we already knew about Plasmo and were asked about it in a very interesting discussion with one of our reviewers for our JOSS submission.

Thank you for offering your help, we will recontact you in a near future to talk about our experimenting of DSP. Concerning modeling beyond linear, it is not in our mid-term goals to do so. We would like to focus our attention to graph partitioning and structure exploiting methods in MILPs. If you have any interesting ideas or references on these subjects, we would love to hear them.

Best Regards,

Miftari B