# Source code for pyproximal.projection.Intersection

import numpy as np
from pyproximal.projection import HyperPlaneBoxProj

[docs]class IntersectionProj():
r"""Intersection of multiple convex sets

Parameters
----------
k : :obj:int
Size of vector to be projected
n : :obj:int
Number of vectors to be projected simultaneously
sigma : :obj:np.ndarray
Matrix of distances of size :math:k \times k
k : :obj:int, optional
Number of iterations
tol : :obj:float, optional
Tolerance of update

Notes
-----
Given an Intersection of simple sets defined as:

.. math::

K = \bigcap_{1 \leq i_1 < i_2 \leq k} K_{i_1,i_2}, \quad
K_{i_1,i_2}= \{ \mathbf{x}: |x_{i_2} - x_{i_1}| \leq \sigma_{i1, i2} \}

its orthogonal projection can be obtained using the Dykstra's
algorithm [1]_.

.. [1] A., Chambolle, D., Cremers, and T., Pock, "A Convex Approach to
Minimal Partitions", Journal of Mathematical, 2011.

"""
def __init__(self, k, n, sigma, niter=100, tol=1e-5):
self.k, self.n = k, n
if isinstance(sigma, np.ndarray):
self.sigma = sigma
else:
self.sigma = sigma * np.ones((k, k))
self.niter = niter
self.tol = tol

def __call__(self, x):
x = x.reshape(self.k, self.n)
x12 = np.zeros((self.k, self.k, self.n))
for iiter in range(self.niter):
xold = x.copy()
for i1 in range(self.k - 1):
for i2 in range(i1 + 1, self.k):
xtilde = x[i2] - x[i1] + x12[i1, i2]
xtildeabs = np.abs(xtilde)
xdtilde = \
np.maximum(0, xtildeabs - self.sigma[i1, i2]) * \
xtilde / (xtildeabs + 1e-10)
x[i1] = x[i1] + 0.5 * (xdtilde - x12[i1, i2])
x[i2] = x[i2] - 0.5 * (xdtilde - x12[i1, i2])
x12[i1, i2] = xdtilde
if max(np.sum(np.abs(x - xold), axis=0)) < self.tol:
break
return x.ravel()