func LogLoss(Ytrue, X mat.Matrix, Theta, Ypred, Ydiff, grad *mat.Dense, Alpha, L1Ratio float64, nSamples int, activation Activation, disableRegularizationOfFirstFeature bool) (J float64) {
	Ypred.Mul(X, Theta)
	Ypred.Apply(func(i, o int, xtheta float64) float64 { return activation.F(xtheta) }, Ypred)
	Ydiff.Sub(Ypred, Ytrue)
	J = 0.
	Ypred.Apply(func(i int, o int, hpred float64) float64 {
		eps := 1e-10
		y := Ytrue.At(i, o)
		h := hpred
		if hpred == 0. {
			h += eps
		} else if h == 1. {
			h -= eps
		}
		J += -y * math.Log(h)
		if math.IsNaN(J) {
			panic(fmt.Errorf("J Nan after -y*math.Log(h) . y=%g h=%g", y, h))
		}
		if math.IsInf(J, 1) {
			panic(fmt.Errorf("LogLoss J Inf h=%g y=%g", h, y))
		}
		return hpred
	}, Ypred)
	if grad != nil {
		grad.Apply(func(j, o int, theta float64) float64 {
			g := 0.
			for i := 0; i < nSamples; i++ {
				h := Ypred.At(i, o)
				g += -Ytrue.At(i, o) * activation.Fprime(h) / h
			}
			return g
		}, Theta)
	}
	if Alpha > 0. {
		J += regularization(Theta, grad, Alpha, L1Ratio, disableRegularizationOfFirstFeature)
	}
	J /= float64(nSamples)
	if grad != nil {
		grad.Scale(1./float64(nSamples), grad)
	}
	if math.IsNaN(J) {
		panic("J Nan")
	}
	return
}
LogLoss for one versus rest classifiers
package linearmodel
import (
	"fmt"
	"math"
	"github.com/pa-m/sklearn/base"
	"gonum.org/v1/gonum/mat"
)
type Loss func(Ytrue, X mat.Matrix, Theta, Ypred, Ydiff, grad *mat.Dense, Alpha, L1Ratio float64, nSamples int, activation Activation, disableRegularizationOfFirstFeature bool) (J float64)
var LossFunctions = map[string]Loss{"square": SquareLoss, "log": LogLoss, "cross-entropy": CrossEntropyLoss}
func SquareLoss(Ytrue, X mat.Matrix, Theta, Ypred, Ydiff, grad *mat.Dense, Alpha, L1Ratio float64, nSamples int, activation Activation, disableRegularizationOfFirstFeature bool) (J float64) {
	Ypred.Mul(X, Theta)
	Ypred.Apply(func(i, o int, xtheta float64) float64 { return activation.F(xtheta) }, Ypred)
	Ydiff.Sub(Ypred, Ytrue)
	J = 0.
	Ydiff.Apply(func(_ int, _ int, ydiff float64) float64 {
		J += ydiff * ydiff
		return ydiff
	}, Ydiff)
	if grad != nil {
		if _, ok := activation.(base.Identity); ok {
			grad.Mul(X.T(), Ydiff)  
		} else {
			grad.Apply(func(j, o int, theta float64) float64 {
				g := 0.
				for i := 0; i < nSamples; i++ {
					h := Ypred.At(i, o)
					g += Ydiff.At(i, o) * X.At(i, j) * activation.Fprime(h)
				}
				return g
			}, Theta)
		}
	}
	if Alpha > 0. {
		J += regularization(Theta, grad, Alpha, L1Ratio, disableRegularizationOfFirstFeature)
	}
	J /= 2. * float64(nSamples)
	if grad != nil {
		grad.Scale(1./float64(nSamples), grad)
	}
	return
}