diff --git a/.gitignore b/.gitignore index b20ea41d..f536e2d9 100644 --- a/.gitignore +++ b/.gitignore @@ -121,3 +121,5 @@ ENV/ .gdb_history .cquery_cached_index/ +pip-wheel-metadata/ +scratch/ diff --git a/.travis.yml b/.travis.yml index be040dbd..7c062f65 100644 --- a/.travis.yml +++ b/.travis.yml @@ -21,9 +21,6 @@ env: matrix: include: - - os: linux - env: - - PYTHON_VERSION=2.7 - os: linux env: - PYTHON_VERSION=3.5 @@ -33,9 +30,6 @@ matrix: - os: linux env: - PYTHON_VERSION=3.7 - - os: osx - env: - - PYTHON_VERSION=2.7 - os: osx env: - PYTHON_VERSION=3.5 @@ -45,14 +39,6 @@ matrix: - os: osx env: - PYTHON_VERSION=3.7 - allow_failures: - - os: linux - env: - - PYTHON_VERSION=2.7 - - os: osx - env: - - PYTHON_VERSION=2.7 - install: - if [[ $TRAVIS_OS_NAME == 'osx' ]]; then @@ -63,6 +49,8 @@ install: - bash miniconda.sh -b -p $HOME/miniconda - source $HOME/miniconda/bin/activate - conda config --set always_yes yes --set changeps1 no --set auto_update_conda no + - conda config --add channels conda-forge + - conda config --set channel_priority strict - conda create -n testenv python=$PYTHON_VERSION pytest numpy - source activate testenv - conda install conda conda-verify conda-build anaconda-client twine diff --git a/appveyor.yml b/appveyor.yml index a84c6a12..6dfe580a 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,6 +1,5 @@ platform: - x64 - - x86 image: Visual Studio 2015 @@ -10,7 +9,6 @@ environment: PYPI_PASSWORD: secure: 3gQtEWf4jXJLJKP8oM22oCShTJ3VI+BtRYczfzO3RSQ= matrix: - - PYTHON_VERSION: 2.7 - PYTHON_VERSION: 3.5 - PYTHON_VERSION: 3.6 - PYTHON_VERSION: 3.7 @@ -36,13 +34,11 @@ install: - if "%PLATFORM%"=="x64" set MINICONDA=C:\Miniconda3-x64 # Load the conda root environment, configure and install some packages - '"%MINICONDA%\Scripts\activate.bat"' - - conda config --append channels conda-forge + - conda config --add channels conda-forge + - conda config --set channel_priority strict - conda config --set always_yes yes --set changeps1 no --set auto_update_conda no - conda create -n testenv python=%PYTHON_VERSION% pytest numpy cmake conda conda-verify conda-build anaconda-client twine - activate testenv - - conda install vs2008_express_vc_python_patch - # Fix for 64-bit Python 2.7 builds, courtesy vs2008_express_vc_python_patch - - call setup_x64 - conda info -a test_script: diff --git a/conda-recipe/conda_build_config.yaml b/conda-recipe/conda_build_config.yaml index b43b49a6..70916fa3 100644 --- a/conda-recipe/conda_build_config.yaml +++ b/conda-recipe/conda_build_config.yaml @@ -1,5 +1,5 @@ -numpy: - - 1.11 +# numpy: +# - 1.11 CONDA_BUILD_SYSROOT: - /var/tmp/MacOSX-SDKs/MacOSX10.9.sdk # [osx] diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index dfb1be8d..5ece40c0 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -23,6 +23,7 @@ requirements: - python {{ python }} - numpy {{ numpy }} - pip + - qdldl-python build: - {{ compiler('c') }} - cmake @@ -30,7 +31,7 @@ requirements: - python - {{ pin_compatible('numpy') }} - scipy >=0.13.2 - - future + - qdldl-python test: requires: diff --git a/extension/include/osqpmodulemethods.h b/extension/include/osqpmodulemethods.h index 219fb43b..8c4ed150 100644 --- a/extension/include/osqpmodulemethods.h +++ b/extension/include/osqpmodulemethods.h @@ -5,297 +5,6 @@ * OSQP methods independently from any object * ***********************************************************************/ -static PyObject * OSQP_module_solve(OSQP *self, PyObject *args, PyObject *kwargs) { - c_int n, m; // Problem dimensions - c_int exitflag_setup, exitflag_solve; - - // Variables for setup - PyOSQPData *pydata; - OSQPData * data; - OSQPSettings * settings; - OSQPWorkspace * workspace; // Pointer to C workspace structure - PyArrayObject *Px, *Pi, *Pp, *q, *Ax, *Ai, *Ap, *l, *u; - - - // Variables for solve - // Create status object - PyObject * status; - // Create obj_val object - PyObject * obj_val; - // Create solution objects - PyObject * x, *y, *prim_inf_cert, *dual_inf_cert; - // Define info related variables - static char *argparse_string_info; - PyObject *info_list; - PyObject *info; - - // Results - PyObject *results_list; - PyObject *results; - - npy_intp nd[1]; - npy_intp md[1]; - - static char *kwlist[] = {"dims", // nvars and ncons - "Px", "Pi", "Pp", "q", // Cost function - "Ax", "Ai", "Ap", "l", "u", // Constraints - "scaling", - "adaptive_rho", "adaptive_rho_interval", - "adaptive_rho_tolerance", "adaptive_rho_fraction", - "rho", "sigma", "max_iter", "eps_abs", "eps_rel", - "eps_prim_inf", "eps_dual_inf", "alpha", "delta", - "linsys_solver", "polish", - "polish_refine_iter", "verbose", - "scaled_termination", - "check_termination", "warm_start", - "time_limit", NULL}; // Settings - -#ifdef DLONG - - // NB: linsys_solver is enum type which is stored as int (regardless on how c_int is defined). -#ifdef DFLOAT - static char * argparse_string_setup = "(LL)O!O!O!O!O!O!O!O!O!|LLLffffLffffffiLLLLLLf"; -#else - static char * argparse_string_setup = "(LL)O!O!O!O!O!O!O!O!O!|LLLddddLddddddiLLLLLLd"; -#endif - -#else - -#ifdef DFLOAT - static char * argparse_string_setup = "(ii)O!O!O!O!O!O!O!O!O!|iiiffffiffffffiiiiiiif"; -#else - static char * argparse_string_setup = "(ii)O!O!O!O!O!O!O!O!O!|iiiddddiddddddiiiiiiid"; -#endif - -#endif - - // Initialize settings - settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings)); - osqp_set_default_settings(settings); - - if( !PyArg_ParseTupleAndKeywords(args, kwargs, argparse_string_setup, kwlist, - &n, &m, - &PyArray_Type, &Px, - &PyArray_Type, &Pi, - &PyArray_Type, &Pp, - &PyArray_Type, &q, - &PyArray_Type, &Ax, - &PyArray_Type, &Ai, - &PyArray_Type, &Ap, - &PyArray_Type, &l, - &PyArray_Type, &u, - &settings->scaling, - &settings->adaptive_rho, - &settings->adaptive_rho_interval, - &settings->adaptive_rho_tolerance, - &settings->adaptive_rho_fraction, - &settings->rho, - &settings->sigma, - &settings->max_iter, - &settings->eps_abs, - &settings->eps_rel, - &settings->eps_prim_inf, - &settings->eps_dual_inf, - &settings->alpha, - &settings->delta, - &settings->linsys_solver, - &settings->polish, - &settings->polish_refine_iter, - &settings->verbose, - &settings->scaled_termination, - &settings->check_termination, - &settings->warm_start, - &settings->time_limit)) { - return (PyObject *) NULL; - } - - // Create Data from parsed vectors - pydata = create_pydata(n, m, Px, Pi, Pp, q, Ax, Ai, Ap, l, u); - data = create_data(pydata); - - // Perform setup and solve - // release the GIL - Py_BEGIN_ALLOW_THREADS; - // Create Workspace object - exitflag_setup = osqp_setup(&(workspace), data, settings); - exitflag_solve = osqp_solve(workspace); - // reacquire the GIL - Py_END_ALLOW_THREADS; - - // Cleanup data and settings - free_data(data, pydata); - c_free(settings); - - // Check successful setup and solve - if (exitflag_setup){ // Workspace allocation error - PyErr_SetString(PyExc_ValueError, "Workspace allocation error!"); - return (PyObject *) NULL; - } - - if(exitflag_solve){ - PyErr_SetString(PyExc_ValueError, "OSQP solve error!"); - return (PyObject *) NULL; - } - - // Temporary solution - nd[0] = (npy_intp)workspace->data->n; // Dimensions in R^n - md[0] = (npy_intp)workspace->data->m; // Dimensions in R^m - - // If problem is not primal or dual infeasible store it - if ((workspace->info->status_val != OSQP_PRIMAL_INFEASIBLE) && - (workspace->info->status_val != OSQP_PRIMAL_INFEASIBLE_INACCURATE) && - (workspace->info->status_val != OSQP_DUAL_INFEASIBLE) && - (workspace->info->status_val != OSQP_DUAL_INFEASIBLE_INACCURATE)){ - - // Primal and dual solutions - x = (PyObject *)PyArrayFromCArray(workspace->solution->x, nd); - y = (PyObject *)PyArrayFromCArray(workspace->solution->y, md); - - // Infeasibility certificates -> None values - prim_inf_cert = PyArray_EMPTY(1, nd, NPY_OBJECT, 0); - dual_inf_cert = PyArray_EMPTY(1, md, NPY_OBJECT, 0); - - } else if (workspace->info->status_val == OSQP_PRIMAL_INFEASIBLE || - workspace->info->status_val == OSQP_PRIMAL_INFEASIBLE_INACCURATE) { - // primal infeasible - - // Primal and dual solution arrays -> None values - x = PyArray_EMPTY(1, nd, NPY_OBJECT, 0); - y = PyArray_EMPTY(1, md, NPY_OBJECT, 0); - - // Primal infeasibility certificate - prim_inf_cert = (PyObject *)PyArrayFromCArray(workspace->delta_y, md); - - // Dual infeasibility certificate -> None values - dual_inf_cert = PyArray_EMPTY(1, nd, NPY_OBJECT, 0); - - // Set objective value to infinity - workspace->info->obj_val = NPY_INFINITY; - - } else { - // dual infeasible - - // Primal and dual solution arrays -> None values - x = PyArray_EMPTY(1, nd, NPY_OBJECT, 0); - y = PyArray_EMPTY(1, md, NPY_OBJECT, 0); - - // Primal infeasibility certificate -> None values - prim_inf_cert = PyArray_EMPTY(1, md, NPY_OBJECT, 0); - - // Dual infeasibility certificate - dual_inf_cert = (PyObject *)PyArrayFromCArray(workspace->delta_x, nd); - - // Set objective value to -infinity - workspace->info->obj_val = -NPY_INFINITY; - } - - /* CREATE INFO OBJECT */ - // Store status string - status = PyUnicode_FromString(workspace->info->status); - - // Store obj_val - if (workspace->info->status_val == OSQP_NON_CVX) { // non convex - obj_val = PyFloat_FromDouble(Py_NAN); - } else { - obj_val = PyFloat_FromDouble(workspace->info->obj_val); - } - -#ifdef PROFILING - -#ifdef DLONG - -#ifdef DFLOAT - argparse_string_info = "LOLLOfffffffLf"; -#else - argparse_string_info = "LOLLOdddddddLd"; -#endif - -#else - -#ifdef DFLOAT - argparse_string_info = "iOiiOfffffffif"; -#else - argparse_string_info = "iOiiOdddddddid"; -#endif - -#endif /* DLONG */ - - info_list = Py_BuildValue(argparse_string_info, - workspace->info->iter, - status, - workspace->info->status_val, - workspace->info->status_polish, - obj_val, - workspace->info->pri_res, - workspace->info->dua_res, - workspace->info->setup_time, - workspace->info->solve_time, - workspace->info->update_time, - workspace->info->polish_time, - workspace->info->run_time, - workspace->info->rho_updates, - workspace->info->rho_estimate - ); -#else /* PROFILING */ - -#ifdef DLONG - -#ifdef DFLOAT - argparse_string = "LOLLOffLf"; -#else - argparse_string = "LOLLOddLd"; -#endif - -#else - -#ifdef DFLOAT - argparse_string = "iOiiOffif"; -#else - argparse_string = "iOiiOddid"; -#endif - -#endif /* DLONG */ - - info_list = Py_BuildValue(argparse_string_info, - workspace->info->iter, - status, - workspace->info->status_val, - workspace->info->status_polish, - obj_val, - workspace->info->pri_res, - workspace->info->dua_res, - workspace->info->rho_updates, - workspace->info->rho_estimate, - ); -#endif /* PROFILING */ - - info = PyObject_CallObject((PyObject *) &OSQP_info_Type, info_list); - - /* Release the info argument list. */ - Py_DECREF(info_list); - - /* CREATE RESULTS OBJECT */ - results_list = Py_BuildValue("OOOOO", x, y, prim_inf_cert, dual_inf_cert, info); - - /* Call the class object. */ - results = PyObject_CallObject((PyObject *) &OSQP_results_Type, results_list); - - // Delete results list - Py_DECREF(results_list); - - // Cleanup workspace - if (osqp_cleanup(workspace)) { - PyErr_SetString(PyExc_ValueError, "Workspace deallocation error!"); - return (PyObject *) NULL; - } - - // Return results - return results; - - -} - - static PyObject *OSQP_constant(OSQP *self, PyObject *args) { char * constant_name; // String less than 32 chars @@ -380,10 +89,9 @@ static PyObject *OSQP_constant(OSQP *self, PyObject *args) { static PyMethodDef OSQP_module_methods[] = { - {"solve", (PyCFunction)OSQP_module_solve,METH_VARARGS|METH_KEYWORDS, PyDoc_STR("Setup solve and cleanup OSQP problem. This function releases GIL.")}, - {"constant", (PyCFunction)OSQP_constant, METH_VARARGS, PyDoc_STR("Return internal OSQP constant")}, - {NULL, NULL} /* sentinel */ + {"constant", (PyCFunction)OSQP_constant, METH_VARARGS, PyDoc_STR("Return internal OSQP constant")}, + {NULL, NULL} /* sentinel */ }; - + #endif diff --git a/extension/include/osqpobjectpy.h b/extension/include/osqpobjectpy.h index 9b01643e..7f80b915 100644 --- a/extension/include/osqpobjectpy.h +++ b/extension/include/osqpobjectpy.h @@ -36,7 +36,7 @@ static c_int OSQP_dealloc(OSQP *self) { // Solve Optimization Problem static PyObject * OSQP_solve(OSQP *self) { c_int exitflag; - + // Create status object PyObject * status; @@ -71,7 +71,12 @@ static PyObject * OSQP_solve(OSQP *self) { /** * Solve QP Problem */ + + // Release the GIL + Py_BEGIN_ALLOW_THREADS; exitflag = osqp_solve(self->workspace); + Py_END_ALLOW_THREADS; + if(exitflag){ PyErr_SetString(PyExc_ValueError, "OSQP solve error!"); return (PyObject *) NULL; @@ -317,7 +322,11 @@ static PyObject * OSQP_setup(OSQP *self, PyObject *args, PyObject *kwargs) { data = create_data(pydata); // Create Workspace object + // Release the GIL + Py_BEGIN_ALLOW_THREADS; exitflag = osqp_setup(&(self->workspace), data, settings); + Py_END_ALLOW_THREADS; + // Cleanup data and settings free_data(data, pydata); diff --git a/extension/src/osqpmodule.c b/extension/src/osqpmodule.c index 24d4da20..b3b36870 100644 --- a/extension/src/osqpmodule.c +++ b/extension/src/osqpmodule.c @@ -9,10 +9,10 @@ /* The PyInt variable is a PyLong in Python3.x. */ -#if PY_MAJOR_VERSION >= 3 -#define PyInt_AsLong PyLong_AsLong -#define PyInt_Check PyLong_Check -#endif +// #if PY_MAJOR_VERSION >= 3 +// #define PyInt_AsLong PyLong_AsLong +// #define PyInt_Check PyLong_Check +// #endif // OSQP Object type @@ -36,8 +36,7 @@ static PyTypeObject OSQP_Type; * Interface Methods * ************************/ - /* Module initialization for Python 3*/ - #if PY_MAJOR_VERSION >= 3 + /* Module initialization*/ static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, "_osqp", /* m_name */ NULL, /* m_doc */ @@ -48,20 +47,14 @@ static PyTypeObject OSQP_Type; NULL, /* m_clear */ NULL, /* m_free */ }; - #endif - static PyObject * moduleinit(void){ PyObject *m; - // Initialize module (no methods. all inside OSQP object) -#if PY_MAJOR_VERSION >= 3 + // Initialize module m = PyModule_Create(&moduledef); -#else - m = Py_InitModule3("_osqp", OSQP_module_methods, NULL); -#endif if (m == NULL) return NULL; @@ -86,20 +79,7 @@ static PyObject * moduleinit(void){ // Init OSQP Internal module -#if PY_MAJOR_VERSION >= 3 PyMODINIT_FUNC PyInit__osqp(void) { -#else -PyMODINIT_FUNC init_osqp(void) { -#endif - import_array(); /* for numpy arrays */ - - // Module initialization is not a global variable in - // Python 3 -#if PY_MAJOR_VERSION >= 3 return moduleinit(); -#else - moduleinit(); -#endif - } diff --git a/module/__init__.py b/module/__init__.py index 497b3ab9..f216e4ce 100644 --- a/module/__init__.py +++ b/module/__init__.py @@ -1,3 +1,2 @@ from osqp.interface import OSQP -from osqp.interface import solve from osqp._osqp import constant diff --git a/module/codegen/files_to_generate/emosqpmodule.c b/module/codegen/files_to_generate/emosqpmodule.c index 645ae077..7859154b 100644 --- a/module/codegen/files_to_generate/emosqpmodule.c +++ b/module/codegen/files_to_generate/emosqpmodule.c @@ -126,10 +126,10 @@ c_float toc(PyTimer* t) { /* The PyInt variable is a PyLong in Python3.x. */ -#if PY_MAJOR_VERSION >= 3 -#define PyInt_AsLong PyLong_AsLong -#define PyInt_Check PyLong_Check -#endif +// #if PY_MAJOR_VERSION >= 3 +// #define PyInt_AsLong PyLong_AsLong +// #define PyInt_Check PyLong_Check +// #endif // Get float type from OSQP setup @@ -664,23 +664,22 @@ static PyObject * OSQP_update_P_A(PyObject *self, PyObject *args) { static PyMethodDef PYTHON_EXT_NAME_methods[] = { - {"solve", (PyCFunction)OSQP_solve, METH_NOARGS, "Solve QP"}, - {"update_lin_cost", (PyCFunction)OSQP_update_lin_cost, METH_VARARGS, "Update linear cost"}, - {"update_lower_bound", (PyCFunction)OSQP_update_lower_bound, METH_VARARGS, "Update lower bound"}, - {"update_upper_bound", (PyCFunction)OSQP_update_upper_bound, METH_VARARGS, "Update upper bound"}, - {"update_bounds", (PyCFunction)OSQP_update_bounds, METH_VARARGS, "Update bounds"}, - #if EMBEDDED != 1 - {"update_P", (PyCFunction)OSQP_update_P, METH_VARARGS, "Update matrix P"}, - {"update_A", (PyCFunction)OSQP_update_A, METH_VARARGS, "Update matrix A"}, - {"update_P_A", (PyCFunction)OSQP_update_P_A, METH_VARARGS, "Update matrices P and A"}, - #endif - {NULL, NULL, 0, NULL} + {"solve", (PyCFunction)OSQP_solve, METH_NOARGS, "Solve QP"}, + {"update_lin_cost", (PyCFunction)OSQP_update_lin_cost, METH_VARARGS, "Update linear cost"}, + {"update_lower_bound", (PyCFunction)OSQP_update_lower_bound, METH_VARARGS, "Update lower bound"}, + {"update_upper_bound", (PyCFunction)OSQP_update_upper_bound, METH_VARARGS, "Update upper bound"}, + {"update_bounds", (PyCFunction)OSQP_update_bounds, METH_VARARGS, "Update bounds"}, +#if EMBEDDED != 1 + {"update_P", (PyCFunction)OSQP_update_P, METH_VARARGS, "Update matrix P"}, + {"update_A", (PyCFunction)OSQP_update_A, METH_VARARGS, "Update matrix A"}, + {"update_P_A", (PyCFunction)OSQP_update_P_A, METH_VARARGS, "Update matrices P and A"}, +#endif + {NULL, NULL, 0, NULL} }; /* Module initialization for Python 3*/ - #if PY_MAJOR_VERSION >= 3 static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, "PYTHON_EXT_NAME", /* m_name */ "Embedded OSQP solver", /* m_doc */ @@ -691,7 +690,6 @@ static struct PyModuleDef moduledef = { NULL, /* m_clear */ NULL, /* m_free */ }; - #endif @@ -699,12 +697,9 @@ static PyObject * moduleinit(void){ PyObject *m; - // Initialize module (no methods. all inside OSQP object) - #if PY_MAJOR_VERSION >= 3 + // Initialize module m = PyModule_Create(&moduledef); - #else - m = Py_InitModule3("PYTHON_EXT_NAME", PYTHON_EXT_NAME_methods, "Embedded OSQP solver"); - #endif + if (m == NULL) return NULL; @@ -715,19 +710,9 @@ static PyObject * moduleinit(void){ // Init Osqp Internal module -#if PY_MAJOR_VERSION >= 3 PyMODINIT_FUNC PyInit_PYTHON_EXT_NAME(void) -#else -PyMODINIT_FUNC initPYTHON_EXT_NAME(void) -#endif { import_array(); /* for numpy arrays */ - // Module initialization is not a global variable in - // Python 3 - #if PY_MAJOR_VERSION >= 3 return moduleinit(); - #else - moduleinit(); - #endif } diff --git a/module/interface.py b/module/interface.py index d551c928..cf0cdde0 100644 --- a/module/interface.py +++ b/module/interface.py @@ -5,10 +5,13 @@ from builtins import object import osqp._osqp as _osqp # Internal low level module import numpy as np +import scipy.sparse as spa +from warnings import warn from platform import system import osqp.codegen as cg import osqp.utils as utils import sys +import qdldl class OSQP(object): @@ -27,6 +30,8 @@ def setup(self, P=None, q=None, A=None, l=None, u=None, **settings): solver settings can be specified as additional keyword arguments """ + # TODO(bart): this will be unnecessary when the derivative will be in C + self._derivative_cache = {'P': P, 'q': q, 'A': A, 'l': l, 'u': u} unpacked_data, settings = utils.prepare_data(P, q, A, l, u, **settings) self._model.setup(*unpacked_data, **settings) @@ -102,6 +107,35 @@ def update(self, q=None, l=None, u=None, if Px is not None and Ax is not None: self._model.update_P_A(Px, Px_idx, len(Px), Ax, Ax_idx, len(Ax)) + + # TODO(bart): this will be unnecessary when the derivative will be in C + # update problem data in self._derivative_cache + if q is not None: + self._derivative_cache["q"] = q + + if l is not None: + self._derivative_cache["l"] = l + + if u is not None: + self._derivative_cache["u"] = u + + if Px is not None: + if Px_idx.size == 0: + self._derivative_cache["P"].data = Px + else: + self._derivative_cache["P"].data[Px_idx] = Px + + if Ax is not None: + if Ax_idx.size == 0: + self._derivative_cache["A"].data = Ax + else: + self._derivative_cache["A"].data[Ax_idx] = Ax + + # delete results from self._derivative_cache to prohibit + # taking the derivative of unsolved problems + if "results" in self._derivative_cache.keys(): + del self._derivative_cache["results"] + def update_settings(self, **kwargs): """ Update OSQP solver settings @@ -198,7 +232,12 @@ def solve(self): Solve QP Problem """ # Solve QP - return self._model.solve() + results = self._model.solve() + + # TODO(bart): this will be unnecessary when the derivative will be in C + self._derivative_cache['results'] = results + + return results def warm_start(self, x=None, y=None): """ @@ -274,18 +313,97 @@ def codegen(self, folder, project_type='', parameters='vectors', cg.codegen(work, folder, python_ext_name, project_type, embedded, force_rewrite, float_flag, long_flag) + def derivative_iterative_refinement(self, rhs, max_iter=20, tol=1e-12): + M = self._derivative_cache['M'] -def solve(P=None, q=None, A=None, l=None, u=None, **settings): - """ - Solve problem of the form + # Prefactor + solver = self._derivative_cache['solver'] - minimize 1/2 x' * P * x + q' * x - subject to l <= A * x <= u + sol = solver.solve(rhs) + for k in range(max_iter): + delta_sol = solver.solve(rhs - M @ sol) + sol = sol + delta_sol - solver settings can be specified as additional keyword arguments. - This function disables the GIL because it internally performs - setup solve and cleanup. + if np.linalg.norm(M @ sol - rhs) < tol: + break + + if k == max_iter - 1: + warn("max_iter iterative refinement reached.") + + return sol + + def adjoint_derivative(self, dx=None, dy_u=None, dy_l=None, + P_idx=None, A_idx=None): + """ + Compute adjoint derivative after solve. """ - unpacked_data, settings = utils.prepare_data(P, q, A, l, u, **settings) - return _osqp.solve(*unpacked_data, **settings) + P, q = self._derivative_cache['P'], self._derivative_cache['q'] + A = self._derivative_cache['A'] + l, u = self._derivative_cache['l'], self._derivative_cache['u'] + + try: + results = self._derivative_cache['results'] + except KeyError: + raise ValueError("Problem has not been solved. " + "You cannot take derivatives. " + "Please call the solve function.") + + if results.info.status != "solved": + raise ValueError("Problem has not been solved to optimality. " + "You cannot take derivatives") + + m, n = A.shape + x = results.x + y = results.y + y_u = np.maximum(y, 0) + y_l = -np.minimum(y, 0) + + if A_idx is None: + A_idx = A.nonzero() + + if P_idx is None: + P_idx = P.nonzero() + + if dy_u is None: + dy_u = np.zeros(m) + if dy_l is None: + dy_l = np.zeros(m) + + # Make sure M matrix exists + if 'M' not in self._derivative_cache: + # Multiply second-third row by diag(y_u)^-1 and diag(y_l)^-1 + # to make the matrix symmetric + inv_dia_y_u = spa.diags(np.reciprocal(y_u + 1e-20)) + inv_dia_y_l = spa.diags(np.reciprocal(y_l + 1e-20)) + M = spa.bmat([ + [P, A.T, -A.T], + [A, spa.diags(A @ x - u) @ inv_dia_y_u, None], + [-A, None, spa.diags(l - A @ x) @ inv_dia_y_l] + ], format='csc') + delta = spa.bmat([[0.001 * spa.eye(n), None], + [None, -0.001 * spa.eye(2 * m)]], format='csc') + self._derivative_cache['M'] = M + self._derivative_cache['solver'] = qdldl.Solver(M + delta) + + rhs = - np.concatenate([dx, dy_u, dy_l]) + + r_sol = self.derivative_iterative_refinement(rhs) + + r_x, r_yu, r_yl = np.split(r_sol, [n, n+m]) + + # Extract derivatives for the constraints + rows, cols = A_idx + dA_vals = (y_u[rows] - y_l[rows]) * r_x[cols] + \ + (r_yu[rows] - r_yl[rows]) * x[cols] + dA = spa.csc_matrix((dA_vals, (rows, cols)), shape=A.shape) + du = - r_yu + dl = r_yl + + # Extract derivatives for the cost (P, q) + rows, cols = P_idx + dP_vals = .5 * (r_x[rows] * x[cols] + r_x[cols] * x[rows]) + dP = spa.csc_matrix((dP_vals, P_idx), shape=P.shape) + dq = r_x + + return (dP, dq, dA, dl, du) diff --git a/module/tests/basic_test.py b/module/tests/basic_test.py index dd526fac..91737992 100644 --- a/module/tests/basic_test.py +++ b/module/tests/basic_test.py @@ -155,20 +155,5 @@ def test_upper_triangular_P(self): # Assert equal nptest.assert_array_almost_equal(res_default.x, res_triu.x) nptest.assert_array_almost_equal(res_default.y, res_triu.y) - nptest.assert_array_almost_equal(res_default.info.obj_val, + nptest.assert_array_almost_equal(res_default.info.obj_val, res_triu.info.obj_val) - - - def test_solve_full_vs_object(self): - - # Solve problem with object - res = self.model.solve() - - # Solve problem with module function - res_module = osqp.solve(P=self.P, q=self.q, A=self.A, l=self.l, u=self.u, - **self.opts) - - # Assert close - nptest.assert_array_almost_equal(res.x, res_module.x) - nptest.assert_array_almost_equal(res.y, res_module.y) - nptest.assert_array_almost_equal(res.info.obj_val, res_module.info.obj_val) diff --git a/module/tests/derivative_test.py b/module/tests/derivative_test.py new file mode 100644 index 00000000..2bff645c --- /dev/null +++ b/module/tests/derivative_test.py @@ -0,0 +1,205 @@ +# Test osqp python module +import osqp +# import osqppurepy as osqp +import numpy as np +import numpy.random as npr +from scipy import sparse +from scipy.optimize import approx_fprime +import numpy.testing as npt + +# Unit Test +import unittest + + +npr.seed(1) + +ATOL = 1e-3 +RTOL = 1e-3 +eps_abs = 1e-10 +eps_rel = 1e-10 +max_iter = 10000 +grad_precision = 1e-08 + + +class derivative_tests(unittest.TestCase): + + def get_prob(self, n=10, m=3, P_scale=1., A_scale=1.): + L = np.random.randn(n, n) + P = sparse.csc_matrix(L.dot(L.T) + 5. * sparse.eye(n)) + x_0 = npr.randn(n) + s_0 = npr.rand(m) + A = sparse.csc_matrix(npr.randn(m, n)) + u = A.dot(x_0) + s_0 + l = -5 - 10 * npr.rand(m) + q = npr.randn(n) + true_x = npr.randn(n) + + return [P, q, A, l, u, true_x] + + def get_grads(self, P, q, A, l, u, true_x): + # Get gradients by solving with osqp + m = osqp.OSQP() + m.setup(P, q, A, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + results = m.solve() + if results.info.status != "solved": + raise ValueError("Problem not solved!") + x = results.x + grads = m.adjoint_derivative(dx=x - true_x) + + return grads + + def test_dl_dq(self, verbose=True): + n, m = 5, 5 + + prob = self.get_prob(n=n, m=m, P_scale=100., A_scale=100.) + P, q, A, l, u, true_x = prob + + def grad(q): + [dP, dq, dA, dl, du] = self.get_grads(P, q, A, l, u, true_x) + return dq + + def f(q): + m = osqp.OSQP() + m.setup(P, q, A, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + res = m.solve() + if res.info.status != "solved": + raise ValueError("Problem not solved!") + x_hat = res.x + + return 0.5 * np.sum(np.square(x_hat - true_x)) + + dq = grad(q) + dq_fd = approx_fprime(q, f, grad_precision) + + if verbose: + print('dq_fd: ', np.round(dq_fd, decimals=4)) + print('dq: ', np.round(dq, decimals=4)) + + npt.assert_allclose(dq_fd, dq, rtol=RTOL, atol=ATOL) + + def test_dl_dP(self, verbose=True): + n, m = 3, 3 + + prob = self.get_prob(n=n, m=m, P_scale=100., A_scale=100.) + P, q, A, l, u, true_x = prob + P_idx = P.nonzero() + + def grad(P_val): + P_qp = sparse.csc_matrix((P_val, P_idx), shape=P.shape) + [dP, dq, dA, dl, du] = self.get_grads(P_qp, q, A, l, u, true_x) + return dP + + def f(P_val): + P_qp = sparse.csc_matrix((P_val, P_idx), shape=P.shape) + m = osqp.OSQP() + m.setup(P_qp, q, A, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + res = m.solve() + if res.info.status != "solved": + raise ValueError("Problem not solved!") + x_hat = res.x + + return 0.5 * np.sum(np.square(x_hat - true_x)) + + dP = grad(P.data) + dP_fd_val = approx_fprime(P.data, f, grad_precision) + dP_fd = sparse.csc_matrix((dP_fd_val, P_idx), shape=P.shape) + dP_fd = (dP_fd + dP_fd.T)/2 + + if verbose: + print('dP_fd: ', np.round(dP_fd.data, decimals=4)) + print('dA: ', np.round(dP.data, decimals=4)) + + npt.assert_allclose(dP.todense(), dP_fd.todense(), rtol=RTOL, atol=ATOL) + + + def test_dl_dA(self, verbose=True): + n, m = 3, 3 + + prob = self.get_prob(n=n, m=m, P_scale=100., A_scale=100.) + P, q, A, l, u, true_x = prob + A_idx = A.nonzero() + + def grad(A_val): + A_qp = sparse.csc_matrix((A_val, A_idx), shape=A.shape) + [dP, dq, dA, dl, du] = self.get_grads(P, q, A_qp, l, u, true_x) + return dA + + def f(A_val): + A_qp = sparse.csc_matrix((A_val, A_idx), shape=A.shape) + m = osqp.OSQP() + m.setup(P, q, A_qp, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + res = m.solve() + if res.info.status != "solved": + raise ValueError("Problem not solved!") + x_hat = res.x + + return 0.5 * np.sum(np.square(x_hat - true_x)) + + dA = grad(A.data) + dA_fd_val = approx_fprime(A.data, f, grad_precision) + dA_fd = sparse.csc_matrix((dA_fd_val, A_idx), shape=A.shape) + + if verbose: + print('dA_fd: ', np.round(dA_fd.data, decimals=4)) + print('dA: ', np.round(dA.data, decimals=4)) + + npt.assert_allclose(dA.todense(), dA_fd.todense(), rtol=RTOL, atol=ATOL) + + def test_dl_dl(self, verbose=True): + n, m = 10, 10 + + prob = self.get_prob(n=n, m=m, P_scale=100., A_scale=100.) + P, q, A, l, u, true_x = prob + + def grad(l): + [dP, dq, dA, dl, du] = self.get_grads(P, q, A, l, u, true_x) + return dl + + def f(l): + m = osqp.OSQP() + m.setup(P, q, A, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + res = m.solve() + if res.info.status != "solved": + raise ValueError("Problem not solved!") + x_hat = res.x + + return 0.5 * np.sum(np.square(x_hat - true_x)) + + dl = grad(l) + dl_fd = approx_fprime(l, f, grad_precision) + + if verbose: + print('dl_fd: ', np.round(dl_fd, decimals=4)) + print('dl: ', np.round(dl, decimals=4)) + + npt.assert_allclose(dl_fd, dl, rtol=RTOL, atol=ATOL) + + def test_dl_du(self, verbose=True): + n, m = 5, 5 + + prob = self.get_prob(n=n, m=m, P_scale=100., A_scale=100.) + P, q, A, l, u, true_x = prob + + def grad(u): + [dP, dq, dA, dl, du] = self.get_grads(P, q, A, l, u, true_x) + return du + + def f(u): + m = osqp.OSQP() + m.setup(P, q, A, l, u, eps_abs=eps_abs, eps_rel=eps_rel, max_iter=max_iter, verbose=False) + res = m.solve() + if res.info.status != "solved": + raise ValueError("Problem not solved!") + x_hat = res.x + + return 0.5 * np.sum(np.square(x_hat - true_x)) + + du = grad(u) + du_fd = approx_fprime(u, f, grad_precision) + + if verbose: + print('du_fd: ', np.round(du_fd, decimals=4)) + print('du: ', np.round(du, decimals=4)) + + npt.assert_allclose(du_fd, du, rtol=RTOL, atol=ATOL) + diff --git a/module/tests/multithread_test.py b/module/tests/multithread_test.py new file mode 100644 index 00000000..b9a3ea2d --- /dev/null +++ b/module/tests/multithread_test.py @@ -0,0 +1,66 @@ +# Test osqp python module +import osqp +# import osqppurepy as osqp +from multiprocessing.pool import ThreadPool +import time +import numpy as np +from scipy import sparse +import scipy as sp + +# Unit Test +import unittest + + +class multithread_tests(unittest.TestCase): + + def test_multithread(self): + data = [] + + n_rep = 20 + + for i in range(n_rep): + m = 1000 + n = 500 + Ad = sparse.random(m, n, density=0.3, format='csc') + b = np.random.randn(m) + + # OSQP data + P = sparse.block_diag([sparse.csc_matrix((n, n)), sparse.eye(m)], format='csc') + q = np.zeros(n+m) + A = sparse.vstack([ + sparse.hstack([Ad, -sparse.eye(m)]), + sparse.hstack([sparse.eye(n), sparse.csc_matrix((n, m))])], format='csc') + l = np.hstack([b, np.zeros(n)]) + u = np.hstack([b, np.ones(n)]) + + data.append((P, q, A, l, u)) + + def f(i): + P, q, A, l, u = data[i] + m = osqp.OSQP() + m.setup(P, q, A, l, u, verbose=False) + m.solve() + + pool = ThreadPool(2) + + tic = time.time() + for i in range(n_rep): + f(i) + t_serial = time.time() - tic + + tic = time.time() + pool.map(f, range(n_rep)) + t_parallel = time.time() - tic + + self.assertLess(t_parallel, t_serial) + + + + + + + + + + + diff --git a/osqp_sources b/osqp_sources index 027d4d69..621d0575 160000 --- a/osqp_sources +++ b/osqp_sources @@ -1 +1 @@ -Subproject commit 027d4d696e0d65ad8bb6e1a41efdd17e391ef134 +Subproject commit 621d0575dcb42168b42a097e70b2628179880654 diff --git a/requirements.txt b/requirements.txt index 777b5a82..a78b5571 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ numpy >= 1.7 scipy >= 0.13.2 -future;python_version<="2.7" +qdldl diff --git a/setup.py b/setup.py index ca33c4d9..e4871dca 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,6 @@ from shutil import copyfile, copy from subprocess import call, check_output -import numpy from setuptools import setup, Extension from setuptools.command.build_ext import build_ext @@ -39,7 +38,7 @@ # necessary to remove OSQP args before passing to setup: if OSQP_ARG_MARK in sys.argv: - sys.argv = sys.argv[0:sys.argv.index(OSQP_ARG_MARK)] + sys.argv = sys.argv[0:sys.argv.index(OSQP_ARG_MARK)] # Add parameters to cmake_args and define_macros cmake_args = ["-DUNITTESTS=OFF"] @@ -49,10 +48,7 @@ # Check if windows linux or mac to pass flag if system() == 'Windows': - if sys.version_info.major == 3: - cmake_args += ['-G', 'Visual Studio 14 2015'] - else: - cmake_args += ['-G', 'Visual Studio 9 2008'] + cmake_args += ['-G', 'Visual Studio 14 2015'] # Differentiate between 32-bit and 64-bit if sys.maxsize // 2 ** 32 > 0: cmake_args[-1] += ' Win64' @@ -75,14 +71,14 @@ "Remove long integers for numpy compatibility. See:\n" + " - https://github.com/numpy/numpy/issues/5906\n" + " - https://github.com/ContinuumIO/anaconda-issues/issues/3823\n" + - "You can reenable long integers by passing: --osqp --long argument.\n") + "You can reenable long integers by passing: " + "--osqp --long argument.\n") cmake_args += ['-DDLONG=OFF'] # Pass python to compiler launched from setup.py define_macros += [('PYTHON', None)] -# Pass python version to cmake -py_version = "%i.%i" % sys.version_info[:2] +# Pass python include dirs to cmake cmake_args += ['-DPYTHON_INCLUDE_DIRS=%s' % sysconfig.get_python_inc()] @@ -92,7 +88,16 @@ osqp_build_dir = os.path.join(osqp_dir, 'build') qdldl_dir = os.path.join(osqp_dir, 'lin_sys', 'direct', 'qdldl') + # Interface files +class get_numpy_include(object): + """Returns Numpy's include path with lazy import. + """ + def __str__(self): + import numpy + return numpy.get_include() + + include_dirs = [ os.path.join(osqp_dir, 'include'), # osqp.h os.path.join(qdldl_dir), # qdldl_interface header to @@ -100,7 +105,7 @@ os.path.join(qdldl_dir, "qdldl_sources", "include"), # qdldl includes for file types os.path.join('extension', 'include'), # auxiliary .h files - numpy.get_include()] # numpy header files + get_numpy_include()] # numpy header files sources_files = glob(os.path.join('extension', 'src', '*.c')) @@ -122,7 +127,7 @@ libraries = [] if system() == 'Linux': libraries += ['rt'] -if system() == 'Windows' and sys.version_info[0] == 3: +if system() == 'Windows': # They moved the stdio library to another place. # We need to include this to fix the dependency libraries += ['legacy_stdio_definitions'] @@ -149,7 +154,8 @@ for f in os.listdir(qdldl_dir) if f.endswith('.c')] cfiles += [os.path.join(qdldl_dir, 'qdldl_sources', 'src', f) - for f in os.listdir(os.path.join(qdldl_dir, 'qdldl_sources', 'src'))] + for f in os.listdir(os.path.join(qdldl_dir, 'qdldl_sources', + 'src'))] osqp_codegen_sources_c_dir = os.path.join(osqp_codegen_sources_dir, 'src') if os.path.exists(osqp_codegen_sources_c_dir): # Create destination directory sh.rmtree(osqp_codegen_sources_c_dir) @@ -160,14 +166,16 @@ # List with OSQP H files hfiles = [os.path.join(osqp_dir, 'include', f) for f in os.listdir(os.path.join(osqp_dir, 'include')) - if f.endswith('.h') and f not in ('qdldl_types.h', 'osqp_configure.h', + if f.endswith('.h') and f not in ('qdldl_types.h', + 'osqp_configure.h', 'cs.h', 'ctrlc.h', 'polish.h', 'lin_sys.h')] hfiles += [os.path.join(qdldl_dir, f) for f in os.listdir(qdldl_dir) if f.endswith('.h')] hfiles += [os.path.join(qdldl_dir, 'qdldl_sources', 'include', f) - for f in os.listdir(os.path.join(qdldl_dir, 'qdldl_sources', 'include')) + for f in os.listdir(os.path.join(qdldl_dir, 'qdldl_sources', + 'include')) if f.endswith('.h')] osqp_codegen_sources_h_dir = os.path.join(osqp_codegen_sources_dir, 'include') if os.path.exists(osqp_codegen_sources_h_dir): # Create destination directory @@ -178,17 +186,21 @@ # List with OSQP configure files configure_files = [os.path.join(osqp_dir, 'configure', 'osqp_configure.h.in'), - os.path.join(qdldl_dir, 'qdldl_sources', 'configure', 'qdldl_types.h.in')] -osqp_codegen_sources_configure_dir = os.path.join(osqp_codegen_sources_dir, 'configure') -if os.path.exists(osqp_codegen_sources_configure_dir): # Create destination directory + os.path.join(qdldl_dir, 'qdldl_sources', 'configure', + 'qdldl_types.h.in')] +osqp_codegen_sources_configure_dir = os.path.join(osqp_codegen_sources_dir, + 'configure') +if os.path.exists(osqp_codegen_sources_configure_dir): sh.rmtree(osqp_codegen_sources_configure_dir) os.makedirs(osqp_codegen_sources_configure_dir) for f in configure_files: # Copy configure files copy(f, osqp_codegen_sources_configure_dir) # Copy cmake files -copy(os.path.join(osqp_dir, 'src', 'CMakeLists.txt'), osqp_codegen_sources_c_dir) -copy(os.path.join(osqp_dir, 'include', 'CMakeLists.txt'), osqp_codegen_sources_h_dir) +copy(os.path.join(osqp_dir, 'src', 'CMakeLists.txt'), + osqp_codegen_sources_c_dir) +copy(os.path.join(osqp_dir, 'include', 'CMakeLists.txt'), + osqp_codegen_sources_h_dir) class build_ext_osqp(build_ext): @@ -243,6 +255,7 @@ def readme(): with open('README.rst') as f: return f.read() + with open('requirements.txt') as f: requirements = f.read().splitlines() @@ -255,7 +268,7 @@ def readme(): package_dir={'osqp': 'module', 'osqppurepy': 'modulepurepy'}, include_package_data=True, # Include package data from MANIFEST.in - setup_requires=["numpy >= 1.7"], + setup_requires=["numpy >= 1.7", "qdldl"], install_requires=requirements, license='Apache 2.0', url="https://osqp.org/",