func enetCoordinateDescent(w *mat.VecDense, alpha, beta float64, X *mat.Dense, Y *mat.VecDense, maxIter int, tol float64, rng base.Intner, random, positive bool) *CDResult {
	gap := tol + 1.
	dwtol := tol
	_, NFeatures := X.Dims()
	R := &mat.VecDense{}
	rmX := X.RawMatrix()
	normColsX := make([]float64, NFeatures)
	for r, rp := 0, 0; r < rmX.Rows; r, rp = r+1, rp+rmX.Stride {
		for c := 0; c < rmX.Cols; c++ {
			e := rmX.Data[rp+c]
			normColsX[c] += e * e
		}
	}
	R.MulVec(X, w)
	R.SubVec(Y, R)
	{
		Y2 := &mat.VecDense{}
		Y2.MulElemVec(Y, Y)
		tol *= mat.Sum(Y2)
	}
	var nIter int
	var wii, dwii, wmax, dwmax, dualNormXtA, cons, l1norm float64
	tmp := &mat.VecDense{}
	XtA := &mat.VecDense{}
	XtArv := XtA.RawVector()
	RNorm2, wNorm2 := &mat.VecDense{}, &mat.VecDense{}
	RNorm2rv, wNorm2rv := RNorm2.RawVector(), wNorm2.RawVector()
	RY := &mat.VecDense{}
	fsign := func(x float64) float64 {
		if x > 0 {
			return 1.
		} else if x == 0. {
			return 0.
		}
		return -1
	}
	for nIter = 0; nIter < maxIter; nIter++ {
		wmax, dwmax = 0., 0.
		var ii int
		for fIter := 0; fIter < NFeatures; fIter++ {
			if random {
				if rng != nil {
					ii = rng.Intn(NFeatures)
				} else {
					ii = rand.Intn(NFeatures)
				}
			} else {
				ii = fIter
			}
			if normColsX[ii] == 0. {
				continue
			}
			wii = w.AtVec(ii)  
			if wii != 0. {
				R.AddScaledVec(R, wii, X.ColView(ii))
			}
			tmp.MulVec(X.ColView(ii).T(), R)
			if positive && tmp.AtVec(0) < 0. {
				w.SetVec(ii, 0)
			} else {
				w.SetVec(ii, (fsign(tmp.AtVec(0))*math.Max(math.Abs(tmp.AtVec(0))-alpha, 0))/(normColsX[ii]+beta))
			}
			if w.AtVec(ii) != 0. {
				R.AddScaledVec(R, -w.AtVec(ii), X.ColView(ii))
			}
			dwii = math.Abs(w.AtVec(ii) - wii)
			if dwii > dwmax {
				dwmax = dwii
			}
			if math.Abs(w.AtVec(ii)) > wmax {
				wmax = math.Abs(w.AtVec(ii))
			}
		}
		if wmax == 0. || dwmax/wmax < dwtol || nIter == maxIter-1 {
			XtA.MulVec(X.T(), R)
			XtArv = XtA.RawVector()
			XtA.AddScaledVec(XtA, -beta, w)
			if positive {
				dualNormXtA = mat.Max(XtA)
			} else {
				dualNormXtA = math.Abs(XtA.AtVec(blas64.Iamax(XtArv)))
			}
			RNorm2.MulVec(R.T(), R)
			RNorm2rv = RNorm2.RawVector()
			wNorm2.MulVec(w.T(), w)
			wNorm2rv = wNorm2.RawVector()
			if dualNormXtA > alpha {
				cons = alpha / dualNormXtA
				ANorm2 := RNorm2rv.Data[0] * cons * cons
				gap = .5 * (RNorm2rv.Data[0] + ANorm2)
			} else {
				cons = 1.
				gap = RNorm2rv.Data[0]
			}
			l1norm = blas64.Asum(w.RawVector())
			RY.MulElemVec(R, Y)
			gap += (alpha*l1norm - cons*mat.Sum(RY)) + .5*beta*(1.+cons*cons)*wNorm2rv.Data[0]
			if gap < tol {
				break
			}
		}
	}
	return &CDResult{Gap: gap, Eps: tol, NIter: nIter + 1}
}
coordinate descent algorithm for Elastic-Net
v https://github.com/scikit-learn/scikit-learn/blob/a24c8b464d094d2c468a16ea9f8bf8d42d949f84/sklearn/linear_model/cd_fast.pyx
package linearmodel
import (
	"math"
	"github.com/pa-m/sklearn/base"
	"golang.org/x/exp/rand"
	"gonum.org/v1/gonum/blas/blas64"
	"gonum.org/v1/gonum/mat"
)