some comments in readme

README.rst

@@ -0,0 +1,96 @@
+Project EuroMir
+===============
+
+Project codename `EuroMir <https://rcdb.com/972.htm>`_
+`(also, song) <https://open.spotify.com/track/3ffkbz5OvPjXjOsYTsEjKu>`_
+is a prototype of a conic programming solver that builds on the work done for
+the older `conic program refinement
+<https://github.com/cvxgrp/cone_prog_refine>`_ research project, of which it
+may inherit the name, once it's completed.
+
+Compared to that 2018 Stanford research project, it has a completely new
+codebase (written from scratch) and it removes various unnecessary
+dependencies. It also has a modified algorithm, which is guaranteed to preserve
+strong convexity (unlike many similar attempts). It uses a simplified version
+of the 2018 algorithm only for the final polishing.
+
+Algorithm (draft)
+=================
+
+The algorithm is under development. This is the current model, see the
+`scs <https://web.stanford.edu/~boyd/papers/pdf/scs.pdf>`_ and
+`conic refinement
+<https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf>`_ papers to
+understand the notation:
+
+.. math::
+
+    \begin{array}{ll}
+
+        \text{minimize} & \|Q u - v \|_2^2 + \| u - \Pi u \|_2^2  + \| v - \Pi^\star v \|_2^2
+
+    \end{array}
+
+The system matrix :math:`Q` from the homogeneous self-dual embedding is skew
+symmetric, so at convergence it is guaranteed that :math:`u` and :math:`v` are
+orthogonal. The objective function is clearly convex and has continuous
+derivative. If we drop the zero cone variables (as we do), the projections
+onto :math:`\mathbf{R}`, and any all-zero rows and columns of :math:`Q`, it is
+also strongly convex. The conditioning and Lipschitz constant are all dependent
+on the conditioning of :math:`Q`, we apply by default standard `Ruiz diagonal
+pre-conditioning <https://web.stanford.edu/~takapoui/preconditioning.pdf>`_.
+
+We use `limited-memory BFGS
+<https://doi.org/10.1090/S0025-5718-1980-0572855-7>`_, as it's implemented in
+the `variable-metric limited memory module of MINPACK-2
+<https://ftp.mcs.anl.gov/pub/MINPACK-2/>`_, for the minimization routine.
+
+We then use the 2018 conic refinement algorithm (simplified, without the
+normalization step), for the final refinement of the solution obtained from the
+BFGS loop.
+
+The approach is clearly globally convergent, and it can probably be showed to
+converge super-linearly. However, theoretical bounds on convergence have very
+limited applicability in real mathematical engineering; smart implementation
+choices and good code quality are by far more important. Only experimentation
+can tell how this approach compares to mature alternative ones, such as
+the various flavors of interior point and operator splitting algorithms.
+
+
+Development
+===========
+
+We mirror development in Python and C (and, in the future, C with OpenMP). We
+set up the build with CMake, patching ``setuptools`` to use it. All testing and
+profiling is done in Python. Binding of C to Python is done using ``ctypes``:
+we link at runtime. This guarantees that pre-built wheels are agnostic to the
+Python version, the Numpy ABI, ....
+Memory communication between Python (mostly Numpy arrays) and C is still
+zero-copy and the GIL is released by ``ctypes``; there are no appreciable
+overheads compared to building with ``Python.h``, since the user only calls a
+``ctypes``-linked function once when solving a conic program.
+
+Memory
+======
+
+All memory is pre-allocated by the caller, which in typical usage is the CVXPY
+interface. There are no ``malloc`` in our C code. We require to allocate space
+for a copy of the problem data (for rescaling, that may be prevented if we
+change the internal logic, adding some computational overhead), and the size of
+the primal and dual variables times the L-BFGS memory, which is a small number
+like 5 or 10.
+
+In layman terms, this is the least memory any conic solver needs, and
+dramatically less than interior-point solvers.
+
+Usage
+=====
+
+We will provide a CVXPY and raw Python interface as part of our packages. The
+single C function the user interacts with will be also documented, for usage
+from other runtime environments.
+
+Installation
+============
+
+Pre-built wheels will be available on PyPi soon.

project_euromir/__init__.py

@@ -141,3 +141,12 @@ def fun(**kwargs):
         ('mult', 'double'),
     )
 )
+
+###
+# CVXPY interface
+###
+
+try:
+    from .cvxpy_solver import Solver
+except ImportError:
+    pass

project_euromir/cvxpy_solver.py

@@ -0,0 +1,65 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""CVXPY solver interface.
+
+Documentation:
+    https://www.cvxpy.org/tutorial/solvers/index.html#custom-solvers
+
+Relevant base class:
+    https://github.com/cvxpy/cvxpy/blob/master/cvxpy/reductions/solvers/conic_solvers/conic_solver.py
+
+Model class:
+    https://github.com/cvxpy/cvxpy/blob/master/cvxpy/reductions/solvers/conic_solvers/scs_conif.py
+"""
+
+try:
+    from cvxpy.reductions.solvers.conic_solvers.conic_solver import (
+        ConicSolver, NonNeg, Zero)
+except ImportError as exc: # pragma: no cover
+    raise ImportError(
+        "Can't use CVXPY interface if CVXPY is not installed!") from exc
+
+
+class Solver(ConicSolver):
+    """CVXPY solver interface.
+
+    We follow SCS conventions for simplicity: cones ordering, in the future
+    exp and sdp cone definition, names in the dims dict (SCS 3), ...
+    """
+
+    MIP_CAPABLE = False
+    SUPPORTED_CONSTRAINTS = [Zero, NonNeg]
+    REQUIRES_CONSTR = False
+
+    def name(self):
+        return "PROJECT_EUROMIR"
+
+    def solve_via_data(self, *args, **kwargs):
+        print("Solving with a custom QP solver!")
+        super().solve_via_data(*args, **kwargs)

project_euromir/linear_algebra.py

@@ -0,0 +1,57 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""Basic linear algebra operations."""
+
+from __future__ import annotations
+
+import numpy as np
+import scipy as sp
+
+
+def Q_matvec(
+    matrix: sp.csc_matrix,
+    b: np.array,
+    c: np.array,
+    input: np.array) -> np.array:
+    """Python implementation of HSDE's Q matvec."""
+
+    m = len(b)
+    n = len(c)
+    output = np.zeros(m + n + 1)
+
+    output[:n] += matrix.T @ input[n:n+m]
+    output[:n] += c * input[-1]
+
+    output[n:n+m] -= matrix @ input[:n]
+    output[n:n+m] += b * input[-1]
+
+    output[-1] -= np.dot(c, input[:n])
+    output[-1] -= np.dot(b, input[n:n+m])
+
+    return output

project_euromir/loss.py

@@ -0,0 +1,29 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""Define loss function and gradient."""

project_euromir/tests/test_linear_algebra.py

@@ -34,6 +34,9 @@
 import scipy.sparse as sp
 
 import project_euromir as lib
+from project_euromir.linear_algebra import Q_matvec
+
+from .test_equilibrate import _make_Q
 
 
 class TestLinearAlgebra(TestCase):