Add windows build on CI (#105)

* Add beginning of nnqp implementation.

* Add current work

* Add current work refactoring nnqp and nnls.

* Add constraint integration with working tests for unconstrained case.

* Add test with constraints

* Change constraint API and refactor test with pytest.mark.parametrize

* Fix interaction test

* Add current work on lower bound constraints

* Add separate build for windows

* Change versioning to be localized to __init__.py with partial import.

* Add starting code for qcp sood

* Add glinternet docs

* Add windows into build workflow

* Add a few more flags for windows

* Add changes to pass clang on stringent build

* Add MAMBA_ROOT_PREFIX to GITHUB_ENV

* Add separate conda handling for windows

* Specify bash for windows

* Add debug on finding conda

* Add find on $CONDA

* Add debug on MAMBA_ROOT_PREFIX within CIBW

* Check MAMBA_ROOT_PREFIX from non-micromamba shell

* Update cibuildwheel

* Add environment variable to CIBW_ENVIRONMENT

* Help Linux build with MAMBA_ROOT_PREFIX

* Add debug on include for eigen3

* Find eigen3 in prefix

* Add conda create in linux

* Add to conda to path

* Add fresh install of conda

* Install wget omg

* Change to apt-get

* Access host file system instead

* Try adding as environment to CIBW

* Ignore warnings for msvc around all pybind11 related things

* Make paths general and warnings compiler specific

* Add debug on prefix for windows

* Find eigen3

* Add new setup to adjust for windows mamba pathing

* Remove stringent build on windows

* Change all openmp loops to be signed and u_char -> uint8_t

* Change all u_char

* Change 0/0 to inf to bypass msvc error

* Set warning to 0 not 1 and set MACOSX_DEPLOYMENT_TARGET to 10.9

* Remove cibw_before_build and runtime_library_dirs

* Add no deps and add back runtime dirs

* Add relaxed msvc warnings

* Ignore C4566

* Add correct syntax

* Ignore conversion or truncation

* Ignore more conversions and helps with left-shift warning

* Ignore openmp clauses being ignored

* Add final build fixes and update version

* Clean-up from review 1

* Add more warnings about constraints

* Add a cleaner hack for reading __version__

* Fix build error on unused params

* Fix build errors 2

.github/workflows/wheels.yml

@@ -14,6 +14,7 @@ jobs:
           - [ubuntu-latest, manylinux_x86_64]
           - [macos-14, macosx_arm64]
           - [macos-13, macosx_x86_64]
+          - [windows-latest, win_amd64]
         python: 
           - ["cp39", '3.9']
           - ["cp310", '3.10']
@@ -24,20 +25,27 @@ jobs:
     steps:
     - uses: actions/checkout@v4
     - uses: actions/setup-python@v3
-    - name: Install micromamba
+    - name: Setup micromamba
       uses: mamba-org/setup-micromamba@v1
       with:
         environment-name: adelie
         init-shell: bash
         generate-run-shell: true
+        # Install with no dependencies to avoid installing libcxx.
+        # libcxx seems to mess with delocator.
         create-args: >-
+          --no-deps
           llvm-openmp==11.1.0
           eigen==3.4.0
     - name: Build wheels
-      uses: pypa/cibuildwheel@v2.18.1
+      uses: pypa/cibuildwheel@v2.19.0
       env:
-        CIBW_BEFORE_BUILD: bash tools/wheels/cibw_before_build.sh
+        # Linux build is fully self-contained, so extra care is needed.
+        # The host file system is mounted as /host (https://cibuildwheel.pypa.io/en/stable/faq/).
+        # We can access MAMBA_ROOT_PREFIX of host for the conda environment.
+        CIBW_ENVIRONMENT_LINUX: MAMBA_ROOT_PREFIX="/host/home/runner/micromamba"
         CIBW_BUILD: ${{ matrix.python[0] }}-${{ matrix.buildplat[1] }}
+        CIBW_REPAIR_WHEEL_COMMAND_WINDOWS: bash tools/github/repair_windows_wheels.sh {wheel} {dest_dir}
     - name: Store the distribution packages
       uses: actions/upload-artifact@v4
       with:

MANIFEST.in

@@ -1,4 +1,3 @@
-include VERSION
 graft adelie
 
 global-exclude *.so

adelie/__init__.py

@@ -4,6 +4,8 @@
 #   - Leda Guin
 #   - Ginnie Guin
 
+__version__ = "1.1.35"
+
 # Set environment flags before loading adelie_core
 import os
 
@@ -25,6 +27,7 @@
 from . import adelie_core
 from . import bcd
 from . import configs
+from . import constraint
 from . import cv
 from . import data
 from . import diagnostic
@@ -41,4 +44,4 @@
 from .solver import (
     gaussian_cov, 
     grpnet,
-)
+)

adelie/constraint.py

@@ -0,0 +1,68 @@
+from .adelie_core.constraint import (
+    ConstraintBase32,
+    ConstraintBase64,
+)
+from .glm import _coerce_dtype
+from . import adelie_core as core
+from typing import Union
+import numpy as np
+
+
+def lower(
+    b: np.ndarray,
+    *,
+    max_iters: int =100,
+    tol: float =1e-7,
+    nnls_max_iters: int =10000,
+    nnls_tol: float =1e-7,
+    dtype: Union[np.float32, np.float64] =None,
+):
+    """Creates a lower bound constraint.
+
+    The lower bound constraint is given by :math:`x \\geq -b` where :math:`b \\geq 0`.
+
+    Parameters
+    ----------
+    b : (d,) np.ndarray
+        Lower bound.
+    max_iters: int, optional
+        Maximum number of proximal Newton iterations.
+        Default is ``100``.
+    tol : float, optional
+        Convergence tolerance for proximal Newton.
+        Default is ``1e-7``.
+    nnls_max_iters: int, optional
+        Maximum number of non-negative least squares iterations.
+        Default is ``10000``.
+    nnls_tol: float, optional
+        Maximum number of non-negative least squares iterations.
+        Default is ``1e-7``.
+    dtype : Union[np.float32, np.float64], optional
+        The underlying data type.
+        If ``None``, it is inferred from ``b``,
+        in which case ``b`` must have an underlying data type of
+        ``np.float32`` or ``np.float64``.
+        Default is ``None``.
+
+    Returns
+    -------
+    wrap
+        Wrapper constraint object.
+
+    See Also
+    --------
+    adelie.adelie_core.constraint.ConstraintLowerUpper64
+    """
+    b, dtype = _coerce_dtype(b, dtype)
+
+    core_base = {
+        np.float32: core.constraint.ConstraintLowerUpper32,
+        np.float64: core.constraint.ConstraintLowerUpper64,
+    }[dtype]
+
+    class _lower(core_base):
+        def __init__(self):
+            self._b = np.array(b, copy=True, dtype=dtype)
+            core_base.__init__(self, -1, self._b, max_iters, tol, nnls_max_iters, nnls_tol)
+
+    return _lower()

adelie/diagnostic.py

@@ -376,15 +376,24 @@ def gradient_norms(
 
     .. math::
         \\begin{align*}
-            \\hat{h}_g = \\|\\hat{\\gamma}_g - \\lambda (1-\\alpha) \\omega_g \\beta_g\\|_2  \\quad g=1,\\ldots, G
+            \\hat{h}_g = \\|
+                \\hat{\\gamma}_g - 
+                \\lambda (1-\\alpha) \\omega_g \\beta_g
+            \\|_2  \\quad g=1,\\ldots, G
         \\end{align*}
 
     where
     :math:`\\hat{\\gamma}_g` is the gradient as in ``adelie.diagnostic.gradients``,
     :math:`\\omega_g` is the penalty factor,
     :math:`\\lambda` is the regularization,
+    :math:`\\alpha` is the elastic net proportion,
     and :math:`\\beta_g` is the coefficient block for group :math:`g`.
 
+    .. warning::
+        The group-wise gradient norm is primarily used to check the KKT conditions.
+        We *do not* correct for the case of non-trivial constraints 
+        since it may be too time consuming.
+
     Parameters
     ----------
     grads : (L, p) or (L, p, K) np.ndarray

adelie/optimization.py

@@ -1,6 +1,6 @@
 from .adelie_core.optimization import (
     search_pivot,
     symmetric_penalty,
-    nnls_cov_full,
-    nnls_naive,
+    StateNNLS,
+    StateNNQPFull,
 )

adelie/solver.py

@@ -94,6 +94,7 @@ def gaussian_cov(
     A: np.ndarray,
     v: np.ndarray,
     *,
+    constraints: list =None,
     groups: np.ndarray =None,
     alpha: float =1,
     penalty: np.ndarray =None,
@@ -148,6 +149,12 @@ def gaussian_cov(
         It is typically one of the matrices defined in ``adelie.matrix`` submodule.
     v : (p,) np.ndarray
         Linear term.
+    constraints : (G,) list, optional
+        List of constraints for each group.
+        ``constraints[i]`` is the constraint object corresponding to group ``i``.
+        If ``constraints[i]`` is ``None``, then the ``i`` th group is unconstrained.
+        If ``None``, every group is unconstrained.
+        Default is ``None``.
     groups : (G,) np.ndarray, optional
         List of starting indices to each group where `G` is the number of groups.
         ``groups[i]`` is the starting index of the ``i`` th group. 
@@ -319,6 +326,13 @@ def gaussian_cov(
         screen_set = np.arange(G)[(penalty <= 0) | (alpha <= 0)]
         screen_beta = np.zeros(np.sum(group_sizes[screen_set]), dtype=dtype)
         screen_is_active = np.ones(screen_set.shape[0], dtype=bool)
+        screen_dual_size = 0
+        if not (constraints is None):
+            screen_dual_size = np.sum([
+                0 if constraints[k] is None else constraints[k].dual_size
+                for k in screen_set
+            ], dtype=int)
+        screen_dual = np.zeros(screen_dual_size, dtype=dtype)
         active_set_size = screen_set.shape[0]
         active_set = np.empty(G, dtype=int)
         active_set[:active_set_size] = np.arange(active_set_size)
@@ -342,6 +356,7 @@ def gaussian_cov(
         screen_set = warm_start.screen_set
         screen_beta = warm_start.screen_beta
         screen_is_active = warm_start.screen_is_active
+        screen_dual = warm_start.screen_dual
         active_set_size = warm_start.active_set_size
         active_set = warm_start.active_set
         rsq = warm_start.rsq
@@ -350,13 +365,15 @@ def gaussian_cov(
     state = state_gaussian_cov(
         A=A,
         v=v,
+        constraints=constraints,
         groups=groups,
         group_sizes=group_sizes,
         alpha=alpha,
         penalty=penalty,
         screen_set=screen_set,
         screen_beta=screen_beta,
         screen_is_active=screen_is_active,
+        screen_dual=screen_dual,
         active_set_size=active_set_size,
         active_set=active_set,
         rsq=rsq,
@@ -395,6 +412,7 @@ def grpnet(
     X: np.ndarray,
     glm: Union[glm.GlmBase32, glm.GlmBase64],
     *,
+    constraints: list =None,
     groups: np.ndarray =None,
     alpha: float =1,
     penalty: np.ndarray =None,
@@ -484,6 +502,12 @@ def grpnet(
     glm : Union[adelie.glm.GlmBase32, adelie.glm.GlmBase64, adelie.glm.GlmMultiBase32, adelie.glm.GlmMultiBase64]
         GLM object.
         It is typically one of the GLM classes defined in ``adelie.glm`` submodule.
+    constraints : (G,) list, optional
+        List of constraints for each group.
+        ``constraints[i]`` is the constraint object corresponding to group ``i``.
+        If ``constraints[i]`` is ``None``, then the ``i`` th group is unconstrained.
+        If ``None``, every group is unconstrained.
+        Default is ``None``.
     groups : (G,) np.ndarray, optional
         List of starting indices to each group where `G` is the number of groups.
         ``groups[i]`` is the starting index of the ``i`` th group. 
@@ -682,6 +706,7 @@ def grpnet(
 
     solver_args = {
         "X": X_raw,
+        "constraints": constraints,
         "alpha": alpha,
         "offsets": offsets,
         "lmda_path": lmda_path,
@@ -752,6 +777,13 @@ def grpnet(
             screen_set = np.arange(groups.shape[0])[(penalty <= 0) | (alpha <= 0)]
             screen_beta = np.zeros(np.sum(group_sizes[screen_set]), dtype=dtype)
             screen_is_active = np.ones(screen_set.shape[0], dtype=bool)
+            screen_dual_size = 0
+            if not (constraints is None):
+                screen_dual_size = np.sum([
+                    0 if constraints[k] is None else constraints[k].dual_size
+                    for k in screen_set
+                ], dtype=int)
+            screen_dual = np.zeros(screen_dual_size, dtype=dtype)
             active_set_size = screen_set.shape[0]
             active_set = np.empty(groups.shape[0], dtype=int)
             active_set[:active_set_size] = np.arange(active_set_size)
@@ -761,6 +793,7 @@ def grpnet(
             screen_set = warm_start.screen_set
             screen_beta = warm_start.screen_beta
             screen_is_active = warm_start.screen_is_active
+            screen_dual = warm_start.screen_dual
             active_set_size = warm_start.active_set_size
             active_set = warm_start.active_set
 
@@ -772,6 +805,7 @@ def grpnet(
         solver_args["screen_set"] = screen_set
         solver_args["screen_beta"] = screen_beta
         solver_args["screen_is_active"] = screen_is_active
+        solver_args["screen_dual"] = screen_dual
         solver_args["active_set_size"] = active_set_size
         solver_args["active_set"] = active_set
 
@@ -890,6 +924,13 @@ def grpnet(
             screen_set = np.arange(G)[(penalty <= 0) | (alpha <= 0)]
             screen_beta = np.zeros(np.sum(group_sizes[screen_set]), dtype=dtype)
             screen_is_active = np.ones(screen_set.shape[0], dtype=bool)
+            screen_dual_size = 0
+            if not (constraints is None):
+                screen_dual_size = np.sum([
+                    0 if constraints[k] is None else constraints[k].dual_size
+                    for k in screen_set
+                ], dtype=int)
+            screen_dual = np.zeros(screen_dual_size, dtype=dtype)
             active_set_size = screen_set.shape[0]
             active_set = np.empty(groups.shape[0], dtype=int)
             active_set[:active_set_size] = np.arange(active_set_size)
@@ -900,6 +941,7 @@ def grpnet(
             screen_set = warm_start.screen_set
             screen_beta = warm_start.screen_beta
             screen_is_active = warm_start.screen_is_active
+            screen_dual = warm_start.screen_dual
             active_set_size = warm_start.active_set_size
             active_set = warm_start.active_set
 
@@ -911,6 +953,7 @@ def grpnet(
         solver_args["screen_set"] = screen_set
         solver_args["screen_beta"] = screen_beta
         solver_args["screen_is_active"] = screen_is_active
+        solver_args["screen_dual"] = screen_dual
         solver_args["active_set_size"] = active_set_size
         solver_args["active_set"] = active_set
 

adelie/src/adelie_core.cpp

@@ -14,6 +14,9 @@ PYBIND11_MODULE(adelie_core, m) {
     auto m_configs = m.def_submodule("configs", "Configurations submodule.");
     register_configs(m_configs);
 
+    auto m_constraint = m.def_submodule("constraint", "Constraint submodule.");
+    register_constraint(m_constraint);
+
     auto m_glm = m.def_submodule("glm", "GLM submodule.");
     register_glm(m_glm);
 

adelie/src/bcd.cpp

@@ -244,16 +244,16 @@ py::dict constrained_admm_solver(
     double rho,
     size_t max_iters,
     double tol_abs,
-    double tol_rel
+    double tol_rel,
+    ad::util::rowvec_type<double> x,
+    ad::util::rowvec_type<double> z,
+    ad::util::rowvec_type<double> u
 )
 {
     using sw_t = ad::util::Stopwatch;
 
     const auto m = A.rows();
     const auto d = A.cols();
-    ad::util::rowvec_type<double> x(d); x.setZero();
-    ad::util::rowvec_type<double> z(m); z.setZero();
-    ad::util::rowvec_type<double> u(m); u.setZero();
     ad::util::rowvec_type<double> buff(3*m + 4*d);
     size_t iters;
     sw_t sw;
@@ -307,7 +307,7 @@ py::dict constrained_coordinate_descent_solver(
     return dct;
 }
 
-py::dict constrained_proximal_newton_solver(
+py::dict constrained_proximal_newton_general_solver(
     const Eigen::Ref<const ad::util::rowvec_type<double>>& mu0,
     const Eigen::Ref<const ad::util::rowvec_type<double>>& quad,
     const Eigen::Ref<const ad::util::rowvec_type<double>>& linear,
@@ -334,17 +334,16 @@ py::dict constrained_proximal_newton_solver(
     ad::util::rowvec_type<double> x(d);
     ad::util::rowvec_type<double> mu = mu0;
     ad::util::rowvec_type<double> mu_resid = (
-        linear.matrix() - x.matrix() * A
+        linear.matrix() - mu.matrix() * A
     );
     ad::util::rowvec_type<double> AT_vars = (
         A.array().square().rowwise().sum()
     );
-    ad::util::rowmat_type<double> AAT = A * A.transpose();
 
     sw_t sw;
     sw.start();
-    ad::bcd::constrained::proximal_newton_solver(
-        quad, linear, l1, l2, A, b, AT_vars, AAT,
+    ad::bcd::constrained::proximal_newton_general_solver(
+        quad, linear, l1, l2, A, b, AT_vars,
         max_iters, tol, newton_max_iters, newton_tol, nnls_max_iters, nnls_tol, nnls_dtol,
         iters, x, mu, mu_resid, buff
     );
@@ -365,7 +364,7 @@ void register_bcd(py::module_& m)
     /* constrained */
     m.def("constrained_admm_solver", &constrained_admm_solver);
     m.def("constrained_coordinate_descent_solver", &constrained_coordinate_descent_solver);
-    m.def("constrained_proximal_newton_solver", &constrained_proximal_newton_solver);
+    m.def("constrained_proximal_newton_general_solver", &constrained_proximal_newton_general_solver);
 
     /* unconstrained */
     m.def("unconstrained_brent_solver", &unconstrained_brent_solver);