Formula for dual gap of elastic net in coordinate descent solver

### Describe the bug
The computation of the dual gap for the elastic net in the coordinate descent solver (`enet_coordinate_descent`) might be wrong.

The elastic net minimizes
```
Primal(w) = (1/2) * ||y - X w||_2^2 + alpha * ||w||_1 + beta/2 * ||w||_2^2
```

### Lasso
For the pure Lasso, i.e. `beta=0`, the dual becomes [1]
```
Dual(nu) = -1/2 * ||nu||_2^2 - nu'y   if ||X'nu||_∞ = max(abs(X'nu)) <= alpha   else -∞
```
with `nu = Xw - y` (=minus the residual `R`) as possible dual point. This yields the Lasso dual gap
```
Primal(w) - Dual(nu)
```
In case of `||X'nu||_∞ > alpha`, one uses a down-scaled variable `nu` in the dual.

### Elastic Net
For the elastic net, the dual, see section 5.2.3 in [2], becomes
```
Dual(nu) = -1/2 * ||nu||_2^2 - nu'y - 1/(2*beta) * sum_j (|X_j'nu| - alpha)_+^2
```
with `X_j` the j-th column of `X` and `x_+ = max(0, x)` the positive part.

### References
[1] Kim, Seung-Jean et al. “An Interior-Point Method for Large-Scale $\ell_1$-Regularized Least Squares.” IEEE Journal of Selected Topics in Signal Processing 1 (2007): 606-617. [pdf link](http://stanford.edu/~boyd/papers/pdf/l1_ls.pdf)
[2] Mendler-Dünner, Celestine et al. “Primal-Dual Rates and Certificates.” ICML (2016). [arxiv link](https://arxiv.org/pdf/1602.05205.pdf)

### Expected Results

```python
if beta > =0:
    R = y - X @ w  # note: R = -nu
    R_norm2 = R @ R
    gap = R_norm2 - R @ y + alpha * l1_norm + beta/2 * l2_norm
    gap += 1/(2*beta) * np.sum(np.max(0, np.abs(X.T @ R) - alpha)**2)
```
At least a test for it, similar to the pure Lasso case in https://github.com/scikit-learn/scikit-learn/blob/2c2f31d68c21b7647557b2f776e86c05954d80bf/sklearn/linear_model/tests/test_coordinate_descent.py#L292

### Actual Results

I start to be rattled by the `-beta * w` term in line 205.
https://github.com/scikit-learn/scikit-learn/blob/2c2f31d68c21b7647557b2f776e86c05954d80bf/sklearn/linear_model/_cd_fast.pyx#L205-L236

COMMENT:
Here's an implementation if you want to play with it: 
```python
import numpy as np
from numba import njit
from numpy.linalg import norm
import matplotlib.pyplot as plt

from benchopt.datasets import make_correlated_data


@njit
def st(x, mu):
    if x > mu:
        return x - mu
    if x < - mu:
        return x + mu
    return 0.


def st_vec(x, mu):
    return np.sign(x) * np.maximum(0, np.abs(x) - mu)


@njit
def cd_epoch(w, X, R, lc, alpha, beta):
    n_features = X.shape[1]
    for j in range(n_features):
        old = w[j]
        w[j] = st(w[j] + X[:, j] @ R / lc[j],
                  alpha / lc[j]) / (1 + beta / lc[j])
        diff = old - w[j]
        if diff != 0:
            R += diff * X[:, j]


def cd(X, y, alpha, beta, n_iter):
    lc = norm(X, axis=0) ** 2
    n_samples, n_features = X.shape
    R = np.copy(y)
    w = np.zeros(n_features)
    ny2 = norm(y) ** 2
    primals = np.zeros(n_iter)
    duals1 = np.zeros(n_iter)
    duals2 = np.zeros(n_iter)

    for it in range(n_iter):
        cd_epoch(w, X, R, lc, alpha, beta)
        primals[it] = (R ** 2).sum() / 2 + alpha * \
            norm(w, ord=1) + beta / 2 * norm(w) ** 2

        # unconstrained approach
        theta1 = - R
        duals1[it] = (ny2 / 2 - norm(y + theta1) ** 2 / 2
                      - 1 / (2 * beta) * norm(st_vec(X.T @ theta1, alpha)) ** 2)

        # constrained approach:
        scaling = max(1, norm(X.T @ (- R) + beta * w, ord=np.inf) / alpha)
        duals2[it] = (ny2 / 2 - norm(y - R / scaling) ** 2 / 2
                      - beta * norm(w / scaling) ** 2 / 2)

    return w, primals, duals1, duals2


X, y, _ = make_correlated_data(n_samples=500, n_features=2000, random_state=0)

# normalize y to get roughly 1e-14 as machine prec in gap:
y /= norm(y)
alpha_max = norm(X.T @ y, ord=np.inf)
alpha = alpha_max / 30
beta = 0.1 * alpha

w, primals, duals1, duals2 = cd(X, y, alpha, beta, 1000)

plt.semilogy(primals - np.min(primals), label='primal subopt')
plt.semilogy(primals - duals1, label='gap with unconstrained formulation')
plt.semilogy(primals - duals2, label='gap with Lasso formulation')
plt.legend()
plt.show(block=False)
```

Labels: Bug, module:linear_model