diff --git a/adelie/__init__.py b/adelie/__init__.py index a04be1ff..a6578f86 100644 --- a/adelie/__init__.py +++ b/adelie/__init__.py @@ -2,4 +2,5 @@ from . import bcd from . import matrix from . import state -from . import solver \ No newline at end of file +from . import solver +from .solver import grpnet \ No newline at end of file diff --git a/adelie/solver.py b/adelie/solver.py index b7fca0e9..ca3ec167 100644 --- a/adelie/solver.py +++ b/adelie/solver.py @@ -1,5 +1,6 @@ from . import adelie_core as core from . import logger +from . import matrix import adelie as ad import numpy as np @@ -165,13 +166,138 @@ def solve_basil( return state -def grpnet(): - """ - TODO: - - cap max_strong_size by number of features - - decreasing order of lmdas - - cap max number of lambdas per iter - - alpha <= 0: add all variables to strong_set - - alpha > 0: add all variables with penalty <= 0 to strong_set +def grpnet( + *, + X: np.ndarray | matrix.base, + y: np.ndarray, + groups: np.ndarray, + group_sizes: np.ndarray, + alpha: float, + penalty: np.ndarray, + lmda_path: np.ndarray =None, + max_iters: int =int(1e5), + tol: float =1e-12, + rsq_tol: float =0.9, + rsq_slope_tol: float =1e-3, + rsq_curv_tol: float =1e-3, + newton_tol: float =1e-12, + newton_max_iters: int =1000, + n_threads: int =1, + early_exit: bool =True, + intercept: bool =True, + strong_rule: str ="default", + min_ratio: float =1e-2, + lmda_path_size: int =100, + delta_lmda_path_size: int =5, + delta_strong_size: int =5, + max_strong_size: int =None, + check_state: bool =True, +): + """TODO """ - pass \ No newline at end of file + + if isinstance(X, np.ndarray): + X_raw = X + X = ad.matrix.naive_dense(X_raw, n_threads=n_threads) + _X = X.internal() + else: + assert ( + isinstance(X, matrix.MatrixNaiveBase64) or + isinstance(X, matrix.MatrixNaiveBase32) + ) + _X = X + + dtype = ( + np.float64 + if isinstance(_X, matrix.MatrixNaiveBase64) else + np.float32 + ) + + p = _X.cols() + G = len(groups) + + actual_lmda_path_size = ( + lmda_path_size + if lmda_path is None else + len(lmda_path) + ) + delta_lmda_path_size = np.minimum(delta_lmda_path_size, actual_lmda_path_size) + + max_strong_size = np.minimum(max_strong_size, G) + + X_means = np.empty(p, dtype=dtype) + _X.means(X_means) + + X_group_norms = np.empty(G, dtype=dtype) + _X.group_norms( + groups, + group_sizes, + X_means, + intercept, + X_group_norms, + ) + + y_mean = np.mean(y) + yc = y + if intercept: + yc = yc - y_mean + y_var = np.sum(yc ** 2) + resid = yc + + strong_set = np.arange(G)[(penalty <= 0) | (alpha <= 0)] + strong_beta = np.zeros(np.sum(group_sizes[strong_set]), dtype=dtype) + strong_is_active = np.ones(strong_set.shape[0], dtype=bool) + + rsq = 0 + lmda = np.inf + grad = np.empty(p, dtype=dtype) + _X.bmul(0, p, resid, grad) + + if not (lmda_path is None): + lmda_path = np.flip(np.sort(lmda_path)) + + state = ad.state.basil_naive( + X=X, + X_means=X_means, + X_group_norms=X_group_norms, + y_mean=y_mean, + y_var=y_var, + resid=resid, + groups=groups, + group_sizes=group_sizes, + alpha=alpha, + penalty=penalty, + strong_set=strong_set, + strong_beta=strong_beta, + strong_is_active=strong_is_active, + rsq=rsq, + lmda=lmda, + grad=grad, + lmda_path=lmda_path, + lmda_max=None, + edpp_safe_set=None, + edpp_v1_0=None, + edpp_resid_0=None, + max_iters=max_iters, + tol=tol, + rsq_tol=rsq_tol, + rsq_slope_tol=rsq_slope_tol, + rsq_curv_tol=rsq_curv_tol, + newton_tol=newton_tol, + newton_max_iters=newton_max_iters, + n_threads=n_threads, + early_exit=early_exit, + intercept=intercept, + strong_rule=strong_rule, + min_ratio=min_ratio, + lmda_path_size=lmda_path_size, + delta_lmda_path_size=delta_lmda_path_size, + delta_strong_size=delta_strong_size, + max_strong_size=max_strong_size, + ) + + if check_state: + assert isinstance(X_raw, np.ndarray) + state.check(X_raw, y, method="assert") + + return solve_basil(state) diff --git a/adelie/src/include/adelie_core/matrix/matrix_cov_base.hpp b/adelie/src/include/adelie_core/matrix/matrix_cov_base.hpp index 479d2ffe..501fe78c 100644 --- a/adelie/src/include/adelie_core/matrix/matrix_cov_base.hpp +++ b/adelie/src/include/adelie_core/matrix/matrix_cov_base.hpp @@ -9,8 +9,8 @@ class MatrixCovBase { public: using value_t = ValueType; - using rowvec_t = util::rowvec_type; - using colmat_t = util::colmat_type; + using vec_value_t = util::rowvec_type; + using colmat_value_t = util::colmat_type; virtual ~MatrixCovBase() {} @@ -26,8 +26,8 @@ class MatrixCovBase */ virtual void bmul( int i, int j, int p, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) =0; /** @@ -41,7 +41,7 @@ class MatrixCovBase */ virtual void to_dense( int i, int j, int p, int q, - Eigen::Ref out + Eigen::Ref out ) const =0; /** diff --git a/adelie/src/include/adelie_core/matrix/matrix_cov_dense.hpp b/adelie/src/include/adelie_core/matrix/matrix_cov_dense.hpp index f83e0178..1275a559 100644 --- a/adelie/src/include/adelie_core/matrix/matrix_cov_dense.hpp +++ b/adelie/src/include/adelie_core/matrix/matrix_cov_dense.hpp @@ -12,8 +12,8 @@ class MatrixCovDense: public MatrixCovBase using base_t = MatrixCovBase; using dense_t = DenseType; using typename base_t::value_t; - using typename base_t::rowvec_t; - using typename base_t::colmat_t; + using typename base_t::vec_value_t; + using typename base_t::colmat_value_t; private: const Eigen::Map _mat; // underlying dense matrix @@ -32,8 +32,8 @@ class MatrixCovDense: public MatrixCovBase void bmul( int i, int j, int p, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { auto outm = out.matrix(); @@ -48,7 +48,7 @@ class MatrixCovDense: public MatrixCovBase void to_dense( int i, int j, int p, int q, - Eigen::Ref out + Eigen::Ref out ) const override { out = _mat.block(i, j, p, q); diff --git a/adelie/src/include/adelie_core/matrix/matrix_cov_lazy.hpp b/adelie/src/include/adelie_core/matrix/matrix_cov_lazy.hpp index 61c0b029..34fad3c1 100644 --- a/adelie/src/include/adelie_core/matrix/matrix_cov_lazy.hpp +++ b/adelie/src/include/adelie_core/matrix/matrix_cov_lazy.hpp @@ -14,8 +14,8 @@ class MatrixCovLazy: public: using base_t = MatrixCovBase::Scalar>; using typename base_t::value_t; - using typename base_t::rowvec_t; - using typename base_t::colmat_t; + using typename base_t::vec_value_t; + using typename base_t::colmat_value_t; using dense_t = DenseType; using index_t = Eigen::Index; @@ -72,8 +72,8 @@ class MatrixCovLazy: void bmul( int i, int j, int p, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { if (i < 0 || j < 0 || p <= 0 || q <= 0) { @@ -96,7 +96,7 @@ class MatrixCovLazy: void to_dense( int i, int j, int p, int q, - Eigen::Ref out + Eigen::Ref out ) const override { if (i < 0 || i > static_cast(_index_map.size())) { diff --git a/adelie/src/include/adelie_core/matrix/matrix_naive_base.hpp b/adelie/src/include/adelie_core/matrix/matrix_naive_base.hpp index e13165c8..66f1fe77 100644 --- a/adelie/src/include/adelie_core/matrix/matrix_naive_base.hpp +++ b/adelie/src/include/adelie_core/matrix/matrix_naive_base.hpp @@ -4,13 +4,15 @@ namespace adelie_core { namespace matrix { -template +template class MatrixNaiveBase { public: using value_t = ValueType; - using rowvec_t = util::rowvec_type; - using colmat_t = util::colmat_type; + using index_t = IndexType; + using vec_value_t = util::rowvec_type; + using vec_index_t = util::rowvec_type; + using colmat_value_t = util::colmat_type; virtual ~MatrixNaiveBase() {} @@ -23,7 +25,7 @@ class MatrixNaiveBase */ virtual value_t cmul( int j, - const Eigen::Ref& v + const Eigen::Ref& v ) const =0; /** @@ -36,7 +38,7 @@ class MatrixNaiveBase virtual void ctmul( int j, value_t v, - Eigen::Ref out + Eigen::Ref out ) const =0; /** @@ -49,8 +51,8 @@ class MatrixNaiveBase */ virtual void bmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) =0; /** @@ -63,8 +65,8 @@ class MatrixNaiveBase */ virtual void btmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) =0; /** @@ -76,7 +78,32 @@ class MatrixNaiveBase */ virtual void to_dense( int j, int q, - Eigen::Ref out + Eigen::Ref out + ) const =0; + + /** + * @brief Computes column-wise mean. + * + * @param out resulting column means. + */ + virtual void means( + Eigen::Ref out + ) const =0; + + /** + * @brief Computes group-wise column norms. + * + * @param groups mapping group number to starting position. + * @param group_sizes mapping group number to group size. + * @param center true to compute centered column norms. + * @param out resulting group-wise column norms. + */ + virtual void group_norms( + const Eigen::Ref& groups, + const Eigen::Ref& group_sizes, + const Eigen::Ref& means, + bool center, + Eigen::Ref out ) const =0; /** diff --git a/adelie/src/include/adelie_core/matrix/matrix_naive_dense.hpp b/adelie/src/include/adelie_core/matrix/matrix_naive_dense.hpp index d58c4d9e..986ab7fe 100644 --- a/adelie/src/include/adelie_core/matrix/matrix_naive_dense.hpp +++ b/adelie/src/include/adelie_core/matrix/matrix_naive_dense.hpp @@ -12,8 +12,9 @@ class MatrixNaiveDense: public MatrixNaiveBase using base_t = MatrixNaiveBase; using dense_t = DenseType; using typename base_t::value_t; - using typename base_t::rowvec_t; - using typename base_t::colmat_t; + using typename base_t::vec_value_t; + using typename base_t::vec_index_t; + using typename base_t::colmat_value_t; private: const Eigen::Map _mat; // underlying dense matrix @@ -32,7 +33,7 @@ class MatrixNaiveDense: public MatrixNaiveBase value_t cmul( int j, - const Eigen::Ref& v + const Eigen::Ref& v ) const override { return ddot(_mat.col(j).matrix(), v.matrix(), _n_threads); @@ -41,7 +42,7 @@ class MatrixNaiveDense: public MatrixNaiveBase void ctmul( int j, value_t v, - Eigen::Ref out + Eigen::Ref out ) const override { dax(v, _mat.col(j), _n_threads, out); @@ -49,8 +50,8 @@ class MatrixNaiveDense: public MatrixNaiveBase void bmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { auto outm = out.matrix(); @@ -65,8 +66,8 @@ class MatrixNaiveDense: public MatrixNaiveBase void btmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { auto outm = out.matrix(); @@ -81,7 +82,7 @@ class MatrixNaiveDense: public MatrixNaiveBase void to_dense( int j, int q, - Eigen::Ref out + Eigen::Ref out ) const override { dmmeq( @@ -91,6 +92,46 @@ class MatrixNaiveDense: public MatrixNaiveBase ); } + void means( + Eigen::Ref out + ) const override + { + const size_t p = _mat.cols(); + const int n_blocks = std::min(_n_threads, p); + const int block_size = p / n_blocks; + const int remainder = p % n_blocks; + #pragma omp parallel for schedule(static) num_threads(n_blocks) + for (int t = 0; t < n_blocks; ++t) { + const auto begin = ( + std::min(t, remainder) * (block_size + 1) + + std::max(t-remainder, 0) * block_size + ); + const auto size = block_size + (t < remainder); + out.segment(begin, size) = _mat.middleCols(begin, size).colwise().mean(); + } + } + + void group_norms( + const Eigen::Ref& groups, + const Eigen::Ref& group_sizes, + const Eigen::Ref& means, + bool center, + Eigen::Ref out + ) const override + { + const auto n = _mat.rows(); + const auto G = groups.size(); + const auto n_threads_capped = std::min(_n_threads, G); + #pragma omp parallel for schedule(static) num_threads(n_threads_capped) + for (int i = 0; i < G; ++i) { + const auto g = groups[i]; + const auto gs = group_sizes[i]; + auto Xi_fro = _mat.middleCols(g, gs).squaredNorm(); + if (center) Xi_fro -= n * means.segment(g, gs).matrix().squaredNorm(); + out[i] = std::sqrt(Xi_fro); + } + } + int rows() const override { return _mat.rows(); diff --git a/adelie/src/matrix.cpp b/adelie/src/matrix.cpp index 031eadd9..f0200307 100644 --- a/adelie/src/matrix.cpp +++ b/adelie/src/matrix.cpp @@ -16,13 +16,14 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase /* Inherit the constructors */ using base_t::base_t; using typename base_t::value_t; - using typename base_t::rowvec_t; - using typename base_t::colmat_t; + using typename base_t::vec_value_t; + using typename base_t::vec_index_t; + using typename base_t::colmat_value_t; /* Trampoline (need one for each virtual function) */ value_t cmul( int j, - const Eigen::Ref& v + const Eigen::Ref& v ) const override { PYBIND11_OVERRIDE_PURE( @@ -36,7 +37,7 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase void ctmul( int j, value_t v, - Eigen::Ref out + Eigen::Ref out ) const override { PYBIND11_OVERRIDE_PURE( @@ -49,8 +50,8 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase void bmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { PYBIND11_OVERRIDE_PURE( @@ -63,8 +64,8 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase void btmul( int j, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { PYBIND11_OVERRIDE_PURE( @@ -77,7 +78,7 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase void to_dense( int j, int q, - Eigen::Ref out + Eigen::Ref out ) const override { PYBIND11_OVERRIDE_PURE( @@ -88,6 +89,34 @@ class PyMatrixNaiveBase : public ad::matrix::MatrixNaiveBase ); } + void means( + Eigen::Ref out + ) const override + { + PYBIND11_OVERLOAD_PURE( + void, + base_t, + means, + out + ); + } + + void group_norms( + const Eigen::Ref& groups, + const Eigen::Ref& group_sizes, + const Eigen::Ref& means, + bool center, + Eigen::Ref out + ) const override + { + PYBIND11_OVERLOAD_PURE( + void, + base_t, + group_norms, + groups, group_sizes, means, center, out + ); + } + int rows() const override { PYBIND11_OVERRIDE_PURE( @@ -190,6 +219,37 @@ void matrix_naive_base(py::module_& m, const char* name) out : (n, q) np.ndarray Matrix to store the dense result. )delimiter") + .def("means", &internal_t::means, R"delimiter( + Computes column-wise means. + + Equivalent to ``np.mean(X, axis=0)``. + + Parameters + ---------- + out : (p,) np.ndarray + Vector to store the column-wise means. + )delimiter") + .def("group_norms", &internal_t::group_norms, R"delimiter( + Computes group-wise column norms. + + Equivalent to ``np.linalg.norm(Xc[:, g:g+gs])`` + for every group ``g`` with group size ``gs``. + Note that if ``X`` is to be centered first, + ``Xc`` is the column-wise centered version of ``X``. + + Parameters + ---------- + groups : (G,) np.ndarray + Mapping group number to the starting column index. + group_sizes : (G,) np.ndarray + Mapping group number to the group size. + means : (p,) np.ndarray + Column-wise means. + center : bool + ``True`` if the function should compute centered group-wise column norms. + out : (G,) np.ndarray + Resulting group-wise column norms. + )delimiter") .def("rows", &internal_t::rows, R"delimiter( Number of rows. )delimiter") @@ -207,13 +267,13 @@ class PyMatrixCovBase : public ad::matrix::MatrixCovBase /* Inherit the constructors */ using base_t::base_t; using typename base_t::value_t; - using typename base_t::rowvec_t; - using typename base_t::colmat_t; + using typename base_t::vec_value_t; + using typename base_t::colmat_value_t; void bmul( int i, int j, int p, int q, - const Eigen::Ref& v, - Eigen::Ref out + const Eigen::Ref& v, + Eigen::Ref out ) override { PYBIND11_OVERRIDE_PURE( @@ -226,7 +286,7 @@ class PyMatrixCovBase : public ad::matrix::MatrixCovBase void to_dense( int i, int j, int p, int q, - Eigen::Ref out + Eigen::Ref out ) const override { PYBIND11_OVERRIDE_PURE( diff --git a/adelie/state.py b/adelie/state.py index 1628de22..4bf73aca 100644 --- a/adelie/state.py +++ b/adelie/state.py @@ -1739,7 +1739,7 @@ def basil_naive( lmda_path_size: int =100, delta_lmda_path_size: int =5, delta_strong_size: int =5, - max_strong_size: int =1000, + max_strong_size: int =None, ): """Creates a basil, naive method state object. @@ -1877,13 +1877,44 @@ def basil_naive( Maximum number of strong groups allowed. The function will return a valid state and guaranteed to have strong set size less than or equal to ``max_strong_size``. - Default is ``1000``. + If ``None``, it will be set to the total number of groups. + Default is ``None``. See Also -------- adelie.state.basil_naive_64 adelie.state.basil_naive_32 """ + if max_strong_size is None: + max_strong_size = len(groups) + + if max_iters < 0: + raise ValueError("max_iters must be >= 0.") + if tol <= 0: + raise ValueError("tol must be > 0.") + if rsq_tol < 0 or rsq_tol > 1: + raise ValueError("rsq_tol must be in [0,1].") + if rsq_slope_tol < 0: + raise ValueError("rsq_slope_tol must be >= 0.") + if rsq_curv_tol < 0: + raise ValueError("rsq_curv_tol must be >= 0.") + if newton_tol < 0: + raise ValueError("newton_tol must be >= 0.") + if newton_max_iters < 0: + raise ValueError("newton_max_iters must be >= 0.") + if n_threads < 1: + raise ValueError("n_threads must be >= 1.") + if min_ratio <= 0: + raise ValueError("min_ratio must be > 0.") + if lmda_path_size < 0: + raise ValueError("lmda_path_size must be >= 0.") + if delta_lmda_path_size < 1: + raise ValueError("delta_lmda_path_size must be >= 1.") + if delta_strong_size < 1: + raise ValueError("delta_strong_size must be >= 1.") + if max_strong_size < 0: + raise ValueError("max_strong_size must be >= 0.") + if isinstance(X, matrix.base): X_intr = X.internal() else: diff --git a/research/test.ipynb b/research/test.ipynb index 94bb9700..81ae1d2d 100644 --- a/research/test.ipynb +++ b/research/test.ipynb @@ -12,12 +12,13 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "import adelie as ad\n", - "import numpy as np" + "import numpy as np\n", + "import matplotlib.pyplot as plt" ] }, { @@ -58,46 +59,20 @@ " X /= np.sqrt(n)\n", " y /= np.sqrt(n)\n", "\n", - " X_means = np.mean(X, axis=0)\n", - " X_c = X - intercept * X_means[None]\n", - " X_group_norms = np.array([\n", - " np.linalg.norm(X_c[:, g:g+gs])\n", - " for g, gs in zip(groups, group_sizes)\n", - " ])\n", - " y_mean = np.mean(y)\n", - " y_c = y - y_mean * intercept\n", - " y_var = np.sum(y_c ** 2)\n", - " resid = y_c\n", - " strong_set = np.arange(G)[(penalty <= 0) | (alpha <= 0)]\n", - " strong_beta = np.zeros(np.sum(group_sizes[strong_set]))\n", - " strong_is_active = np.ones(strong_set.shape[0], dtype=bool)\n", - " grad = X_c.T @ y_c\n", - "\n", " return {\n", " \"X\": X, \n", " \"y\": y,\n", - " \"X_means\": X_means,\n", - " \"X_group_norms\": X_group_norms,\n", - " \"y_mean\": y_mean,\n", - " \"y_var\": y_var,\n", - " \"resid\": resid,\n", " \"groups\": groups,\n", " \"group_sizes\": group_sizes,\n", " \"alpha\": alpha,\n", " \"penalty\": penalty,\n", - " \"strong_set\": strong_set,\n", - " \"strong_beta\": strong_beta,\n", - " \"strong_is_active\": strong_is_active,\n", - " \"rsq\": 0,\n", - " \"lmda\": np.inf,\n", - " \"grad\": grad,\n", " \"intercept\": intercept,\n", " }" ] }, { "cell_type": "code", - "execution_count": 70, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -114,78 +89,80 @@ ")\n", "test_data[\"penalty\"] = np.sqrt(test_data[\"group_sizes\"])\n", "X, y = test_data[\"X\"], test_data[\"y\"]\n", - "test_data.pop(\"y\")\n", "Xc = X - np.mean(X, axis=0)[None] * intercept\n", "yc = y - np.mean(y) * intercept" ] }, { "cell_type": "code", - "execution_count": 72, + "execution_count": 32, "metadata": {}, "outputs": [], "source": [ - "n_threads = 8\n", - "test_data[\"X\"] = ad.matrix.naive_dense(X, n_threads=n_threads)\n", - "state = ad.state.basil_naive(\n", + "state = ad.grpnet(\n", " **test_data,\n", " max_strong_size=30000,\n", " strong_rule=\"default\",\n", - " n_threads=n_threads,\n", + " n_threads=8,\n", + " check_state=True,\n", ")" ] }, { "cell_type": "code", - "execution_count": 73, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 2min 44s, sys: 953 ms, total: 2min 45s\n", - "Wall time: 20.7 s\n" - ] - } - ], - "source": [ - "%%time\n", - "next_state = ad.solver.solve_basil(state)" - ] - }, - { - "cell_type": "code", - "execution_count": 74, + "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "6.047451920000001" + "(6.257613330000001, 5.790161783, 0.005346733, 8.978549406)" ] }, - "execution_count": 74, + "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "np.sum(next_state.benchmark_fit)" + "(\n", + " np.sum(state.benchmark_fit),\n", + " np.sum(state.benchmark_screen),\n", + " np.sum(state.benchmark_invariance),\n", + " np.sum(state.benchmark_kkt),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "next_state = state.update_path([state.lmdas[-1] * 0.8])" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "next_state = ad.solver.solve_basil(next_state)" ] }, { "cell_type": "code", - "execution_count": 75, + "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "(6.047451920000001, 5.5086994350000005, 0.005143656, 8.708577768000001)" + "(1.387490109, 0.937606798, 0.000393396, 0.330873943)" ] }, - "execution_count": 75, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } @@ -201,30 +178,30 @@ }, { "cell_type": "code", - "execution_count": 76, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "array([0.32461099, 0.50084198, 0.45982327, 0.46059465, 0.5058221 ,\n", - " 0.440893 , 0.441547 , 0.44808808, 0.47491148, 0.60637055,\n", - " 0.56224193, 0.66115996, 0.70434021, 0.63886951, 0.80449459,\n", - " 0.67396849])" + "array([0.31834238, 0.48176995, 0.44621483, 0.43986034, 0.49144022,\n", + " 0.42791791, 0.43047491, 0.42954165, 0.45668795, 0.57706865,\n", + " 0.54215909, 0.64862796, 0.68794999, 0.61746694, 0.76613428,\n", + " 0.66037768])" ] }, - "execution_count": 76, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "next_state.benchmark_kkt" + "state.benchmark_kkt" ] }, { "cell_type": "code", - "execution_count": 79, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -233,21 +210,21 @@ "(1183, 3003)" ] }, - "execution_count": 79, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(\n", - " len(next_state.strong_set),\n", - " len(next_state.edpp_safe_set),\n", + " len(state.strong_set),\n", + " len(state.edpp_safe_set),\n", ")" ] }, { "cell_type": "code", - "execution_count": 81, + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -257,18 +234,18 @@ "\twith 18460 stored elements in Compressed Sparse Row format>" ] }, - "execution_count": 81, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "next_state.betas" + "state.betas" ] }, { "cell_type": "code", - "execution_count": 82, + "execution_count": 35, "metadata": {}, "outputs": [ { @@ -285,7 +262,27 @@ " 0.89464705, 0.90342317])" ] }, - "execution_count": 82, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "state.rsqs" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.9367745])" + ] + }, + "execution_count": 36, "metadata": {}, "output_type": "execute_result" } diff --git a/tests/test_matrix.py b/tests/test_matrix.py index dc075530..a0d0e1b4 100644 --- a/tests/test_matrix.py +++ b/tests/test_matrix.py @@ -42,6 +42,33 @@ def _test(n, p, dtype, order, seed=0): cX.to_dense(0, p // 2, out) assert np.allclose(X[:, :p//2], out) + # test means + X_means = np.empty(p, dtype=dtype) + cX.means(X_means) + assert np.allclose(np.mean(X, axis=0), X_means) + + # test group_means + groups = np.concatenate([ + [0], + np.random.choice(np.arange(1, p), size=p//2-1, replace=False) + ]) + groups = np.sort(groups).astype(int) + group_sizes = np.concatenate([groups, [p]], dtype=int) + group_sizes = group_sizes[1:] - group_sizes[:-1] + out = np.empty(len(groups), dtype=dtype) + cX.group_norms( + groups, + group_sizes, + np.mean(X, axis=0), + True, + out, + ) + expected = np.array([ + np.linalg.norm(X[:, g:g+gs] - X_means[g:g+gs][None]) + for g, gs in zip(groups, group_sizes) + ]) + assert np.allclose(expected, out) + assert cX.rows() == n assert cX.cols() == p