Commit 450725b9 authored by Elvira Clement's avatar Elvira Clement
Browse files

Importing code

parent 3532e829
......@@ -23,9 +23,8 @@ setup(
packages=find_packages(exclude=('tests', 'docs')),
ext_modules=cythonize([
"src/screening/utils.pyx",
"src/screening/gap_ptest.pyx",
"src/screening/gap_rqtest.pyx",
"src/screening/kappa_test.pyx",
"src/screening/gap_test_all.pyx",
"src/screening/gap_test_p_1.pyx",
]),
)
# -*- coding: utf-8 -*-
__version__ = "0.9.0"
\ No newline at end of file
__version__ = "2.0.0"
\ No newline at end of file
......@@ -12,18 +12,16 @@ class AbstractGapScreening(ABC):
self.name = "NONAME"
def get_name(self):
return self.name
@abstractmethod
def apply_test(self, Atu_abs, gap, lbd, vec_gammas, offset_radius=0, index=None) -> np.ndarray:
def apply_test(self, Atu_abs, gap, lbd, vec_gammas, coeff_dual_scaling=1., offset_radius=0, index=None) -> np.ndarray:
""" GAP Sphere screening test detailed in Cédric's node
Parameters
----------
Atu_abs : np.ndarray
vector |matA.T @ vecu| where vecu is dual admissible
vector |matA.T @ vecu|
vecu is not necessarily dual admissible
In that case, the dual scaling factor must be provided
size [n]
gap : positive float
duality gap
......@@ -32,6 +30,11 @@ class AbstractGapScreening(ABC):
vec_gammas : np.ndarray
slope parameters
size [n,]
coeff_dual_scaling : positif float
If coeff_dual_scaling is not feasible, dual scaling factor
such taht vecu / coeff_dual_scaling os dual feasible
Here for code optimization purposes
Default value is 1. (vecu is feasible)
offset_radius : float
additive term added to the redius
index : np.ndarray
......
This diff is collapsed.
# -*- coding: utf-8 -*-
# Based on https://cython.readthedocs.io/en/latest/src/tutorial/numpy.html
# https://cython.readthedocs.io/en/latest/src/tutorial/memory_allocation.html
import numpy as np
from src.screening.ascreening import AbstractGapScreening
# "cimport" is used to import special compile-time information
# about the numpy module (this is stored in a file numpy.pxd which is
# currently part of the Cython distribution).
cimport numpy as np
# from libcpp cimport bool as bool_t
cimport cython
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
ctypedef np.npy_bool bool_t
#(AbstractGapScreening)
cdef class GAP_RQtest:
""" Generalized test
"""
cdef double *vec_cumsum_gammas
cdef long *vec_rq
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
def __init__(
self,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
np.ndarray[long, ndim=1] vec_rq,
):
cdef int n = vec_gammas.shape[0]
cdef int i
self.vec_cumsum_gammas = <double*> PyMem_Malloc((n+1) * sizeof(double))
self.vec_cumsum_gammas[0] = 0
self.vec_rq = <long*> PyMem_Malloc(n * sizeof(long))
for i in range(n):
self.vec_cumsum_gammas[i+1] = self.vec_cumsum_gammas[i] + vec_gammas[i]
self.vec_rq[i] = vec_rq[i]
def __dealloc__(self):
PyMem_Free(self.vec_cumsum_gammas) # no-op if self.data is NULL
PyMem_Free(self.vec_rq)
def get_name(self):
return "p-test cython"
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
def apply_test(self,
np.ndarray[np.npy_double, ndim=1] Atcabs,
double gap,
double lbd,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
double offset_radius=0.,
np.ndarray[long, ndim=1] index=None
):
""" Apply the Gap safe sphere screening test
Implementation of the safe screening test
Parameters
----------
Atcabs : np.ndarray
vector matA.T @ vecc where vec is dual admissible in absolute value
size [n]
gap : positive float
duality gap
lbd : positive float
regularization parameters
vec_gammas : np.ndarray
slope parameters
size [n,]
offset_radius : float
additive term added to the redius
default is 0
index : np.ndarray
Array of indices that sort Atu in absolute value
default is None
Returns
-------
calI_screen : np.ndarray
vector of boolean
True if screening test passes and False otherwise
size [n,]
"""
cdef double radius = np.sqrt(2 * gap) + offset_radius
cdef int n = Atcabs.shape[0]
cdef int k, l, q, r_q
cdef double lhs
# right-hand side
cdef np.ndarray[np.npy_double, ndim=1] vec_kappaqr
vec_kappaqr = np.zeros(n+1, dtype=np.double)
for q in range(1, n+1):
r_q = self.vec_rq[q-1]
vec_kappaqr[q] += lbd * (self.vec_cumsum_gammas[q] - self.vec_cumsum_gammas[r_q-1])
vec_kappaqr[q] -= (q - r_q + 1) * radius
# output
cdef np.ndarray[np.npy_bool, ndim=1] calI_screen = np.ones(n, dtype=bool)
# 1. Sort in descending order
if index is None:
index = np.argsort(Atcabs)[::-1]
cdef np.ndarray[np.npy_double, ndim=1] vec_cumsumAtcabs = np.zeros(n+1, dtype=np.double)
vec_cumsumAtcabs[1:] = np.cumsum(Atcabs[index])
for l in range(n):
index_l = index[l]
for q in range(1, n+1): # Index are shifted by 1
r_q = self.vec_rq[q-1]
# lhs = 0
if r_q >= q:
lhs = Atcabs[index_l]
# elif index_l < index[r_q-1] or index_l >= index[q-1]:
elif l < r_q or l > q-1:
#In this calse |A_{\l}c|_[k] =|Atc|_[k] for the range of k
lhs = Atcabs[index_l] + vec_cumsumAtcabs[q-1] - vec_cumsumAtcabs[r_q-2]
# for k in range(r_q-1, q-1):
# lhs += Atcabs[index[k]]
else:
# print("ici else")
lhs = vec_cumsumAtcabs[q] - vec_cumsumAtcabs[r_q-2]
# for k in range(r_q-1, q):
# lhs += Atcabs[index[k]]
if lhs >= vec_kappaqr[q]:
calI_screen[index_l] = False
break
return calI_screen
\ No newline at end of file
......@@ -18,9 +18,7 @@ from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
ctypedef np.npy_bool bool_t
#(AbstractGapScreening)
cdef class GAP_Ptest:
cdef class GapTestAll:
""" Generalized test
"""
......@@ -47,22 +45,21 @@ cdef class GAP_Ptest:
# self.vec_cumsum_gammas = &buf[0]
self.name = "GAP test-all (cython)"
def __dealloc__(self):
PyMem_Free(self.vec_cumsum_gammas) # no-op if self.data is NULL
def get_name(self):
return "p-test cython"
@cython.boundscheck(False) # turn off bounds-checking for entire function
# @cython.wraparound(False) # turn off negative index wrapping for entire function
def apply_test(self,
np.ndarray[np.npy_double, ndim=1] Atcabs,
double gap,
double lbd,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
double coeff_dual_scaling =1.,
double offset_radius=0.,
np.ndarray[long, ndim=1] index=None
):
......@@ -82,6 +79,11 @@ cdef class GAP_Ptest:
vec_gammas : np.ndarray
slope parameters
size [n,]
coeff_dual_scaling : positif float
If coeff_dual_scaling is not feasible, dual scaling factor
such taht vecu / coeff_dual_scaling os dual feasible
Here for code optimization purposes
Default value is 1. (vecu is feasible)
offset_radius : float
additive term added to the redius
default is 0
......@@ -101,7 +103,9 @@ cdef class GAP_Ptest:
# all indexes below range from 1 to n (instead of 0 to n-1)
# in order to match the paper indexation rules
cdef double radius = np.sqrt(2 * gap) + offset_radius
cdef double radius = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius
# cdef double lbd_aug = coeff_dual_scaling * lbd
cdef np.ndarray[np.npy_double, ndim=1] coeff_lbd_gamma = coeff_dual_scaling * lbd * vec_gammas
cdef int n = Atcabs.shape[0]
cdef int k, q, r
......@@ -121,10 +125,10 @@ cdef class GAP_Ptest:
# 2. Precomputing quantities
cdef np.ndarray[np.npy_double, ndim=1] vec_f = np.cumsum(
lbd * vec_gammas[::-1] - Atcabs[index[::-1]] - radius
coeff_lbd_gamma[::-1] - Atcabs[index[::-1]] - radius
)[::-1]
cdef double best_bound_q = vec_f[0] - lbd * vec_gammas[0]
cdef double best_bound_q = vec_f[0] - coeff_lbd_gamma[0]
cdef double curr_bound_q = 0.
cdef double best_bound_p = vec_f[0]
cdef double curr_bound_p = 0.
......@@ -140,7 +144,7 @@ cdef class GAP_Ptest:
vec_p_star[k] = vec_p_star[k-1]
# 1. Evaluate q*
curr_bound_q = vec_f[k-1] - lbd * vec_gammas[k-1]
curr_bound_q = vec_f[k-1] - coeff_lbd_gamma[k-1]
if curr_bound_q > best_bound_q:
best_bound_q = curr_bound_q
vec_q_star[k] = k
......@@ -155,7 +159,7 @@ cdef class GAP_Ptest:
p = vec_p_star[q]
# Evaluation of the threshold
tau = vec_f[p-1] - vec_f[q-1] + (lbd * vec_gammas[q-1] - radius)
tau = vec_f[p-1] - vec_f[q-1] + (coeff_lbd_gamma[q-1] - radius)
# Test
if Atcabs[index[l-1]] >= tau:
......
......@@ -20,7 +20,7 @@ from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
#(AbstractGapScreening)
cdef class Kappa_test:
cdef class GapTestPequalOne:
""" Generalized test
"""
......@@ -72,14 +72,13 @@ cdef class Kappa_test:
# (1. + vec_coherence_function) * np.arange(1, n+1, dtype=np.double)
# ).copy()
self.name = "Gap test p=1 (cython)"
def __dealloc__(self):
PyMem_Free(self.vec_cumsum_gammas) # no-op if self.data is NULL
PyMem_Free(self.vec_kappa_q) # no-op if self.data is NULL
def get_name(self):
return "kappa-test cython"
@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False) # turn off negative index wrapping for entire function
......@@ -87,7 +86,8 @@ cdef class Kappa_test:
np.ndarray[np.npy_double, ndim=1] Atcabs,
double gap,
double lbd,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
np.ndarray[np.npy_double, ndim=1] vec_gammas,
double coeff_dual_scaling = 1.,
double offset_radius=0.,
np.ndarray[long, ndim=1] index=None
):
......@@ -107,6 +107,14 @@ cdef class Kappa_test:
vec_gammas : np.ndarray
slope parameters
size [n,]
vec_gammas : np.ndarray
slope parameters
size [n,]
coeff_dual_scaling : positif float
If coeff_dual_scaling is not feasible, dual scaling factor
such taht vecu / coeff_dual_scaling os dual feasible
Here for code optimization purposes
Default value is 1. (vecu is feasible)
offset_radius : float
additive term added to the redius
default is 0
......@@ -127,7 +135,8 @@ cdef class Kappa_test:
cdef int l_min = 0
cdef double bound
cdef double radius = np.sqrt(2 * gap) + offset_radius
cdef double radius = coeff_dual_scaling * np.sqrt(2 * gap) + offset_radius
cdef double coeff_lbd = coeff_dual_scaling * lbd
cdef np.ndarray[np.npy_double, ndim=1] sigmas = np.zeros(n+1)
......@@ -142,7 +151,7 @@ cdef class Kappa_test:
sigmas[1:] = np.cumsum(Atcabs[index])
for q in range(n-1, -1, -1):
bound = lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
bound = coeff_lbd * self.vec_cumsum_gammas[q] - self.vec_kappa_q[q] * radius - sigmas[q]
if Atcabs[index[q]] < bound:
l_q = 0
......
......@@ -4,15 +4,15 @@ import numpy as np
from src.screening.ascreening import AbstractGapScreening
class GapSphereSingleTest(AbstractGapScreening):
class GapTestPequalQ(AbstractGapScreening):
""" Single test safe sphere screening test for Slope
"""
def __init__(self):
super(GapSphereSingleTest, self).__init__()
super(GapTestPequalQ, self).__init__()
self.name = "single test"
self.name = "Gap test p=q"
def apply_test(self, Atu_abs, gap, lbd, vec_gammas, offset_radius=0, index=None) -> np.ndarray:
def apply_test(self, Atu_abs, gap, lbd, vec_gammas, coeff_dual_scaling=1., offset_radius=0, index=None) -> np.ndarray:
""" GAP Sphere screening test detailed in Cédric's node
Parameters
......@@ -27,6 +27,11 @@ class GapSphereSingleTest(AbstractGapScreening):
vec_gammas : np.ndarray
slope parameters
size [n,]
coeff_dual_scaling : positif float
If coeff_dual_scaling is not feasible, dual scaling factor
such taht vecu / coeff_dual_scaling os dual feasible
Here for code optimization purposes
Default value is 1. (vecu is feasible)
offset_radius : float
additive term added to the redius
index : np.ndarray
......@@ -43,4 +48,4 @@ class GapSphereSingleTest(AbstractGapScreening):
"""
radius = np.sqrt(2 * gap)
return Atu_abs < lbd * vec_gammas[-1] - radius - offset_radius
return Atu_abs < coeff_dual_scaling * (lbd * vec_gammas[-1] - radius) - offset_radius
......@@ -8,31 +8,53 @@ from src.screening.ascreening import AbstractGapScreening
class EnumLipchitzOptions(enum.Enum):
EXACT = enum.auto()
GERSHGORIN = enum.auto()
EXACT = enum.auto()
GERSHGORIN = enum.auto()
BACKTRACKING = enum.auto()
class DualScalingOptions(enum.Enum):
EXACT = enum.auto()
BAO_ET_AL = enum.auto()
class MyList:
def __getitem__(self, key):
return float(key)
def __getitem__(self, key):
return float(key)
@dataclass
class SlopeParameters(object):
max_it: int = 1e3
gap_stopping: float = 1e-9
time_stopping: float = 1e2
accelerated: bool = False
vecx_init: np.ndarray = None
verbose: bool = True
max_it: int = 1e3
gap_stopping: float = 1e-9
time_stopping: float = 1e2
accelerated: bool = False
vecx_init: np.ndarray = None
verbose: bool = True
dual_scaling: DualScalingOptions = DualScalingOptions.EXACT
eval_gap: bool = True
lipchitz_constant: float = None
lipchitz_update: EnumLipchitzOptions = EnumLipchitzOptions.EXACT
coherence: float = 1.
coherence_function: np.ndarray = MyList()
screening1: AbstractGapScreening = None
screening2: AbstractGapScreening = None
screening_it_div: float = 2.
save_nactive: bool = False
save_pfunc: bool = False
save_dfunc: bool = False
# @dataclass
# class SlopeParameters(AlgParameters):
lipchitz_constant: float = None
lipchitz_update: EnumLipchitzOptions = EnumLipchitzOptions.EXACT
coherence: float = 1.
coherence_function: np.ndarray = MyList()
screening1: AbstractGapScreening = None
screening2: AbstractGapScreening = None
screening_it_div: float = 2.
save_nactive = False
......@@ -3,8 +3,8 @@ import time
import numpy as np
from src.solver.parameters import SlopeParameters, EnumLipchitzOptions
from src.solver.prox import prox_owl
from src.parameters import SlopeParameters, EnumLipchitzOptions, DualScalingOptions
from src.prox import prox_owl
def primal_func(vecy, Ax, x, lbd, vec_gammas) -> float:
......@@ -85,6 +85,7 @@ def slope_gp(vecy, matA, lbd, vecgamma, algParameters=SlopeParameters()):
# Precomputed quantities
Aty = matA.T @ vecy
normy2 = np.linalg.norm(vecy)**2
gamma_returned = vecgamma[::-1]
# Coherence
mu = algParameters.coherence
......@@ -124,6 +125,8 @@ def slope_gp(vecy, matA, lbd, vecgamma, algParameters=SlopeParameters()):
screening_test_1 = algParameters.screening1
screening_test_2 = algParameters.screening2
vec_cumsum_gammas = np.cumsum(vecgamma)
is_test_1 = False if screening_test_1 is None else True
is_test_2 = False if screening_test_2 is None else True
......@@ -145,10 +148,13 @@ def slope_gp(vecy, matA, lbd, vecgamma, algParameters=SlopeParameters()):
best_costfunc = .5 * normy2
best_vecx = np.copy(vecx_hat)
best_vecu = np.zeros(m)
best_dualfunc = 0
it = 0
# Security to prevent error with negative gap
gap_stop = max(algParameters.gap_stopping, 5e-16)
time_starting = time.time()
ellapsed_time = lambda: time.time() - time_starting
while True:
......@@ -160,34 +166,55 @@ def slope_gp(vecy, matA, lbd, vecgamma, algParameters=SlopeParameters()):
vec_res = vecy - Ax
neg_grad = matA[:, ind_active].T @ vec_res
# Dual scaling
index_sort_neg_grad = np.argsort(np.abs(neg_grad))[::-1]
beta_dual = np.abs(neg_grad[index_sort_neg_grad])
beta_dual = np.cumsum(beta_dual) / np.cumsum(vecgamma[:n_active])
max_beta_dual = max(np.max(beta_dual) / lbd, 1.)
vecu_hat = vec_res / max_beta_dual
# -- Evaluating primal function
# I drop the idea below. it is the "best strategy" only when evaluating the dual function
# ytAx = vecy.dot(Ax)
# normAx2 = np.linalg.norm(Ax, 2)**2
half_norm_res_2 = .5 * np.linalg.norm(vec_res, 2)**2
lbd_slope_norm = lbd * gamma_returned.dot(np.sort(np.abs(vecx_hat[ind_active])))
vec_cost_func.append(primal_func(vecy, Ax, vecx_hat[ind_active], lbd, vecgamma[:n_active]))
vec_dual_func.append(dual_func(vecy, normy2, vecu_hat))
vec_cost_func.append(half_norm_res_2 + lbd_slope_norm)
if vec_cost_func[-1] < best_costfunc:
best_costfunc = vec_cost_func[-1]
best_vecx = np.copy(vecx_hat)
if vec_dual_func[-1] > best_dualfunc:
best_dualfunc = vec_dual_func[-1]
best_vecu = np.copy(vecu_hat)
gap = best_costfunc
if (it % 20 == 0) or algParameters.eval_gap:
# -- Evaluating dual function (if needed)
# if algParameters.eval_gap:
index_sort_neg_grad = np.argsort(np.abs(neg_grad))[::-1]
# print(index_sort_neg_grad)
# Dual scaling
if algParameters.dual_scaling == DualScalingOptions.EXACT:
beta_dual = np.cumsum(np.abs(neg_grad[index_sort_neg_grad])) \
/ vec_cumsum_gammas[:n_active]
coeff_dual_scaling = max(np.max(beta_dual) / lbd, 1.)
elif algParameters.dual_scaling == DualScalingOptions.BAO_ET_AL:
beta_dual = np.abs(neg_grad[index_sort_neg_grad]) / vecgamma[:n_active]
coeff_dual_scaling = max(np.max(beta_dual) / lbd, 1.)
else:
raise NotImplemented
gap = np.abs(best_costfunc - best_dualfunc)
# print(f"it {it} gap {gap}")
yT_res = vecy.dot(vec_res)
d_func_val = (yT_res - half_norm_res_2 / coeff_dual_scaling) / coeff_dual_scaling
if algParameters.save_dfunc:
vec_dual_func.append(d_func_val)
gap = np.abs(best_costfunc - d_func_val)
# ------------------
# Stopping rules
# ------------------
if (it >= algParameters.max_it) or (gap <= algParameters.gap_stopping) or (ellapsed_time() > algParameters.time_stopping):
if (it >= algParameters.max_it) or (gap <= gap_stop) or (ellapsed_time() > algParameters.time_stopping):
break
......@@ -195,16 +222,45 @@ def slope_gp(vecy, matA, lbd, vecgamma, algParameters=SlopeParameters()):
# Screening
# ------------------
if is_test_1 or is_test_2:
if (it % 20 == 0) and (is_test_1 or is_test_2):
# -- Evaluating dual function (if needed)
# if algParameters.eval_gap:
index_sort_neg_grad = np.argsort(np.abs(neg_grad))[::-1]
# print(index_sort_neg_grad)
# Dual scaling
if algParameters.dual_scaling == DualScalingOptions.EXACT:
beta_dual = np.cumsum(np.abs(neg_grad[index_sort_neg_grad])) \
/ vec_cumsum_gammas[:n_active]
coeff_dual_scaling = max(np.max(beta_dual) / lbd, 1.)