func (mlp *BaseMultilayerPerceptron64) Predict(X mat.Matrix, Y Mutable) {
	var xb, yb General64
	if xg, ok := X.(RawMatrixer64); ok {
		if yg, ok := Y.(RawMatrixer64); ok {
			xb, yb = General64(xg.RawMatrix()), General64(yg.RawMatrix())
		}
	} else {
		xb.Copy(X)
		yb.Copy(Y)
	}
	mlp.predict(xb.RawMatrix(), yb.RawMatrix())
	FromDense64(Y, yb)
}
Predict do forward pass and fills Y (Y must be Mutable)
package neuralnetwork
import (
	"encoding/json"
	"fmt"
	"log"
	"reflect"
	"runtime"
	"sort"
	"strings"
	"sync"
	"time"
	"gonum.org/v1/gonum/blas/blas64"
	"github.com/pa-m/sklearn/base"
	"golang.org/x/exp/rand"
	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/optimize"
)
type BaseMultilayerPerceptron64 struct {
	Activation         string  `json:"activation"`
	Solver             string  `json:"solver"`
	Alpha              float64 `json:"alpha"`
	WeightDecay        float64 `json:"weight_decay"`
	BatchSize          int     `json:"batch_size"`
	BatchNormalize     bool
	LearningRate       string           `json:"learning_rate"`
	LearningRateInit   float64          `json:"learning_rate_init"`
	PowerT             float64          `json:"power_t"`
	MaxIter            int              `json:"max_iter"`
	LossFuncName       string           `json:"loss_func_name"`
	HiddenLayerSizes   []int            `json:"hidden_layer_sizes"`
	Shuffle            bool             `json:"shuffle"`
	RandomState        base.RandomState `json:"random_state"`
	Tol                float64          `json:"tol"`
	Verbose            bool             `json:"verbose"`
	WarmStart          bool             `json:"warm_start"`
	Momentum           float64          `json:"momentum"`
	NesterovsMomentum  bool             `json:"nesterovs_momentum"`
	EarlyStopping      bool             `json:"early_stopping"`
	ValidationFraction float64          `json:"validation_fraction"`
	Beta1              float64          `json:"beta_1"`
	Beta2              float64          `json:"beta_2"`
	Epsilon            float64          `json:"epsilon"`
	NIterNoChange      int              `json:"n_iter_no_change"`
	NLayers       int
	NIter         int
	NOutputs      int
	Intercepts    [][]float64     `json:"intercepts_"`
	Coefs         []blas64General `json:"coefs_"`
	OutActivation string          `json:"out_activation_"`
	Loss          float64
	t                   int
	LossCurve           []float64
	ValidationScores    []float64
	BestValidationScore float64
	BestLoss            float64
	NoImprovementCount  int
	optimizer           Optimizer64
	packedParameters    []float64
	packedGrads         []float64  
	bestParameters      []float64
	batchNorm           [][]float64
	lb                  *LabelBinarizer64
	beforeMinimize func(optimize.Problem, []float64)
}
var Activations64 = map[string]func(z blas64General){
	"identity": func(z blas64General) {},
	"logistic": func(z blas64General) {
		for row, zpos := 0, 0; row < z.Rows; row, zpos = row+1, zpos+z.Stride {
			for col := 0; col < z.Cols; col++ {
				z.Data[zpos+col] = 1 / (1 + M64.Exp(-z.Data[zpos+col]))
			}
		}
	},
	"tanh": func(z blas64General) {
		for row, zpos := 0, 0; row < z.Rows; row, zpos = row+1, zpos+z.Stride {
			for col := 0; col < z.Cols; col++ {
				z.Data[zpos+col] = M64.Tanh(-z.Data[zpos+col])
			}
		}
	},
	"relu": func(z blas64General) {
		for row, zpos := 0, 0; row < z.Rows; row, zpos = row+1, zpos+z.Stride {
			for col := 0; col < z.Cols; col++ {
				if z.Data[zpos+col] < 0 {
					z.Data[zpos+col] = 0
				}
			}
		}
	},
	"softmax": func(z blas64General) {
		for row, zpos := 0, 0; row < z.Rows; row, zpos = row+1, zpos+z.Stride {
			sum := float64(0)
			for col := 0; col < z.Cols; col++ {
				z.Data[zpos+col] = M64.Exp(z.Data[zpos+col])
				sum += z.Data[zpos+col]
			}
			for col := 0; col < z.Cols; col++ {
				z.Data[zpos+col] /= sum
			}
		}
	},
}
var Derivatives64 = map[string]func(Z, deltas blas64General){
	"identity": func(Z, deltas blas64General) {
	},
	"logistic": func(Z, deltas blas64General) {
		for row, zpos, dpos := 0, 0, 0; row < Z.Rows; row, zpos, dpos = row+1, zpos+Z.Stride, dpos+deltas.Stride {
			for col := 0; col < Z.Cols; col++ {
				z := Z.Data[zpos+col]
				deltas.Data[dpos+col] *= z * (1 - z)
			}
		}
	},
	"tanh": func(Z, deltas blas64General) {
		for row, zpos, dpos := 0, 0, 0; row < Z.Rows; row, zpos, dpos = row+1, zpos+Z.Stride, dpos+deltas.Stride {
			for col := 0; col < Z.Cols; col++ {
				z := Z.Data[zpos+col]
				deltas.Data[dpos+col] *= 1 - z*z
			}
		}
	},
	"relu": func(Z, deltas blas64General) {
		for row, zpos, dpos := 0, 0, 0; row < Z.Rows; row, zpos, dpos = row+1, zpos+Z.Stride, dpos+deltas.Stride {
			for col := 0; col < Z.Cols; col++ {
				if Z.Data[zpos+col] == 0 {
					deltas.Data[dpos+col] = 0
				}
			}
		}
	},
}
var LossFunctions64 = map[string]func(y, h blas64General) float64{
	"square_loss": func(y, h blas64General) float64 {
		sum := float64(0)
		for row, hpos, ypos := 0, 0, 0; row < y.Rows; row, hpos, ypos = row+1, hpos+h.Stride, ypos+y.Stride {
			for col := 0; col < y.Cols; col++ {
				e := h.Data[hpos+col] - y.Data[ypos+col]
				sum += e * e
			}
		}
		return sum / 2 / float64(h.Rows)
	},
	"log_loss": func(y, h blas64General) float64 {
		sum := float64(0)
		hmin, hmax := M64.Nextafter(0, 1), M64.Nextafter(1, 0)
		for row, hpos, ypos := 0, 0, 0; row < y.Rows; row, hpos, ypos = row+1, hpos+h.Stride, ypos+y.Stride {
			for col := 0; col < y.Cols; col++ {
				hval := h.Data[hpos+col]
				if hval < hmin {
					hval = hmin
				} else if hval > hmax {
					hval = hmax
				}
				if y.Data[ypos+col] != 0 {
					sum += -y.Data[ypos+col] * M64.Log(hval)
				}
			}
		}
		return sum / float64(h.Rows)
	},
	"binary_log_loss": func(y, h blas64General) float64 {
		sum := float64(0)
		hmin, hmax := M64.Nextafter(0, 1), M64.Nextafter(1, 0)
		for row, hpos, ypos := 0, 0, 0; row < y.Rows; row, hpos, ypos = row+1, hpos+h.Stride, ypos+y.Stride {
			for col := 0; col < y.Cols; col++ {
				hval := h.Data[hpos+col]
				if hval < hmin {
					hval = hmin
				} else if hval > hmax {
					hval = hmax
				}
				sum += -y.Data[ypos+col]*M64.Log(hval) - (1-y.Data[ypos+col])*M64.Log1p(-hval)
			}
		}
		return sum / float64(h.Rows)
	},
}
type Optimizer64 interface {
	iterationEnds(timeStep float64)
	triggerStopping(msg string, verbose bool) bool
	updateParams(grads []float64)
}
func addIntercepts64(a blas64General, b []float64) {
	for arow, apos := 0, 0; arow < a.Rows; arow, apos = arow+1, apos+a.Stride {
		for c := 0; c < a.Cols; c++ {
			a.Data[apos+c] += b[c]
		}
	}
}
func matRowMean64(a blas64General, b []float64) {
	for c := 0; c < a.Cols; c++ {
		b[c] = 0
	}
	for arow, apos := 0, 0; arow < a.Rows; arow, apos = arow+1, apos+a.Stride {
		for c := 0; c < a.Cols; c++ {
			b[c] += a.Data[apos+c]
		}
	}
	for c := 0; c < a.Cols; c++ {
		b[c] /= float64(a.Rows)
	}
}
func NewBaseMultilayerPerceptron64() *BaseMultilayerPerceptron64 {
	return &BaseMultilayerPerceptron64{
		Activation:       "relu",
		Solver:           "adam",
		Alpha:            0.0001,
		BatchSize:        200,
		LearningRate:     "constant",
		LearningRateInit: 0.001,
		PowerT:           .5,
		MaxIter:          200,
		HiddenLayerSizes: []int{100},
		Shuffle:          true,
		Tol:                1e-4,
		Verbose:            false,
		WarmStart:          false,
		Momentum:           .9,
		NesterovsMomentum:  true,
		EarlyStopping:      false,
		ValidationFraction: .1,
		Beta1:              .9,
		Beta2:              .999,
		Epsilon:            1e-8,
		NIterNoChange:      10,
	}
}
func (mlp *BaseMultilayerPerceptron64) forwardPass(activations []blas64General) {
	hiddenActivation := Activations64[mlp.Activation]
	var i int
	for i = 0; i < mlp.NLayers-1; i++ {
		gemm64(blas.NoTrans, blas.NoTrans, 1, activations[i], mlp.Coefs[i], 0, activations[i+1])
		addIntercepts64(activations[i+1], mlp.Intercepts[i])
		if (i + 1) != (mlp.NLayers - 1) {
			hiddenActivation(activations[i+1])
		}
	}
	i = mlp.NLayers - 2
	outputActivation := Activations64[mlp.OutActivation]
	outputActivation(activations[i+1])
}
func (mlp *BaseMultilayerPerceptron64) batchNormalize(activations []blas64General) {
	for i := 0; i < mlp.NLayers-2; i++ {
		activation := activations[i+1]
		batchNorm := mlp.batchNorm[i]
		for o := 0; o < activation.Cols; o++ {
			M := float64(0)
			for r, rpos := 0, 0; r < activation.Rows; r, rpos = r+1, rpos+activation.Stride {
				a := M64.Abs(activation.Data[rpos+o])
				if M < a {
					M = a
				}
			}
			if M > 0 {
				for r, rpos := 0, 0; r < activation.Rows; r, rpos = r+1, rpos+activation.Stride {
					activation.Data[rpos+o] /= M
				}
			}
			batchNorm[o] = M
		}
	}
}
func (mlp *BaseMultilayerPerceptron64) batchNormalizeDeltas(deltas blas64General, batchNorm []float64) {
	for r, rpos := 0, 0; r < deltas.Rows; r, rpos = r+1, rpos+deltas.Stride {
		for o := 0; o < deltas.Cols; o++ {
			deltas.Data[rpos+o] /= batchNorm[o]
		}
	}
}
func (mlp *BaseMultilayerPerceptron64) sumCoefSquares() float64 {
	s := float64(0)
	for _, c := range mlp.Coefs {
		for _, co := range c.Data {
			s += co * co
		}
	}
	return s
}
func (mlp *BaseMultilayerPerceptron64) computeLossGrad(layer, NSamples int, activations []blas64General, deltas []blas64General, coefGrads []blas64General, interceptGrads [][]float64) {
	gemm64(blas.Trans, blas.NoTrans, 1/float64(NSamples), activations[layer], deltas[layer], 0, coefGrads[layer])
	axpy64(len(coefGrads[layer].Data), mlp.Alpha/float64(NSamples), mlp.Coefs[layer].Data, coefGrads[layer].Data)
	matRowMean64(deltas[layer], interceptGrads[layer])
}
func (mlp *BaseMultilayerPerceptron64) backprop(X, y blas64General, activations, deltas, coefGrads []blas64General, interceptGrads [][]float64) float64 {
	nSamples := X.Rows
	if mlp.WeightDecay > 0 {
		for iw := range mlp.packedParameters {
			mlp.packedParameters[iw] *= (1 - mlp.WeightDecay)
		}
	}
	mlp.forwardPass(activations)
	if mlp.BatchNormalize {
		mlp.batchNormalize(activations)
	}
	lossFuncName := mlp.LossFuncName
	if strings.EqualFold(lossFuncName, "log_loss") && strings.EqualFold(mlp.OutActivation, "logistic") {
		lossFuncName = "binary_log_loss"
	}
	loss := LossFunctions64[lossFuncName](y, activations[len(activations)-1])
	loss += (0.5 * mlp.Alpha) * mlp.sumCoefSquares() / float64(nSamples)
	last := mlp.NLayers - 2
	{
		H := activations[len(activations)-1]
		D := deltas[last]
		for r, pos := 0, 0; r < y.Rows; r, pos = r+1, pos+y.Stride {
			for o, posc := 0, pos; o < y.Cols; o, posc = o+1, posc+1 {
				D.Data[posc] = H.Data[posc] - y.Data[posc]
			}
		}
	}
	mlp.computeLossGrad(
		last, nSamples, activations, deltas, coefGrads, interceptGrads)
	for i := mlp.NLayers - 2; i >= 1; i-- {
		gemm64(blas.NoTrans, blas.Trans, 1, deltas[i], mlp.Coefs[i], 0, deltas[i-1])
		inplaceDerivative := Derivatives64[mlp.Activation]
		inplaceDerivative(activations[i], deltas[i-1])
		if mlp.BatchNormalize {
			mlp.batchNormalizeDeltas(deltas[i-1], mlp.batchNorm[i-1])
		}
		mlp.computeLossGrad(
			i-1, nSamples, activations, deltas, coefGrads,
			interceptGrads)
	}
	return loss
}
func (mlp *BaseMultilayerPerceptron64) initialize(yCols int, layerUnits []int, isClassifier, isMultiClass bool) {
	mlp.NIter = 0
	mlp.t = 0
	mlp.NOutputs = yCols
	mlp.NLayers = len(layerUnits)
	if !isClassifier {
		mlp.OutActivation = "identity"
		mlp.LossFuncName = "square_loss"
	} else if isMultiClass {
		mlp.OutActivation = "softmax"
		mlp.LossFuncName = "log_loss"
	} else {
		mlp.OutActivation = "logistic"
		mlp.LossFuncName = "binary_log_loss"
	}
	mlp.Coefs = make([]blas64General, mlp.NLayers-1)
	mlp.Intercepts = make([][]float64, mlp.NLayers-1)
	off := 0
	for i := 0; i < mlp.NLayers-1; i++ {
		off += (1 + layerUnits[i]) * layerUnits[i+1]
	}
	mem := make([]float64, off)
	mlp.packedParameters = mem[0:off]
	if mlp.BatchNormalize {
		mlp.batchNorm = make([][]float64, mlp.NLayers-2)
	}
	off = 0
	if mlp.RandomState == (base.RandomState)(nil) {
		mlp.RandomState = base.NewLockedSource(uint64(time.Now().UnixNano()))
	}
	type Float64er interface {
		Float64() float64
	}
	var rndFloat64 func() float64
	if float64er, ok := mlp.RandomState.(Float64er); ok {
		rndFloat64 = float64er.Float64
	} else {
		rndFloat64 = rand.New(mlp.RandomState).Float64
	}
	for i := 0; i < mlp.NLayers-1; i++ {
		prevOff := off
		mlp.Intercepts[i] = mem[off : off+layerUnits[i+1]]
		off += layerUnits[i+1]
		mlp.Coefs[i] = blas64General{Rows: layerUnits[i], Cols: layerUnits[i+1], Stride: layerUnits[i+1], Data: mem[off : off+layerUnits[i]*layerUnits[i+1]]}
		off += layerUnits[i] * layerUnits[i+1]
		factor := float64(6.)
		fanIn, fanOut := layerUnits[i], layerUnits[i+1]
		if strings.EqualFold(mlp.Activation, "logistic") {
			factor = 2.
		}
		initBound := M64.Sqrt(factor / float64(fanIn+fanOut))
		for pos := prevOff; pos < off; pos++ {
			mem[pos] = rndFloat64() * initBound
		}
		if mlp.BatchNormalize && i < mlp.NLayers-2 {
			mlp.batchNorm[i] = make([]float64, layerUnits[i+1])
		}
	}
	mlp.BestLoss = M64.Inf(1)
}
func (mlp *BaseMultilayerPerceptron64) fit(X, y blas64General, incremental bool) {
	mlp.validateHyperparameters()
	for _, s := range mlp.HiddenLayerSizes {
		if s < 0 {
			log.Panicf("hiddenLayerSizes must be > 0, got %v.", mlp.HiddenLayerSizes)
		}
	}
	X, y = mlp.validateInput(X, y, incremental)
	nSamples, nFeatures := X.Rows, X.Cols
	mlp.NOutputs = y.Cols
	layerUnits := append([]int{nFeatures}, mlp.HiddenLayerSizes...)
	layerUnits = append(layerUnits, mlp.NOutputs)
	if mlp.RandomState == nil {
		mlp.RandomState = rand.New(base.NewLockedSource(uint64(time.Now().UnixNano())))
	}
	if !mlp.WarmStart && !incremental {
		var isClassifier, isMulticlass = true, y.Cols > 1
		for _, yval := range y.Data {
			if yval != 0 && yval != 1 {
				isClassifier = false
				break
			}
		}
		mlp.initialize(y.Cols, layerUnits, isClassifier, isMulticlass)
	}
	if strings.EqualFold(mlp.Solver, "lbfgs") {
		mlp.BatchSize = nSamples
	} else if mlp.BatchSize <= 0 {
		mlp.BatchSize = nSamples
		if mlp.BatchSize > 200 {
			mlp.BatchSize = 200
		}
	} else {
		if mlp.BatchSize > nSamples {
			log.Printf("Got batchsize larger than sample size.  It is going to be clipped.\n")
			mlp.BatchSize = nSamples
		}
	}
	activations := make([]blas64.General, 1, len(layerUnits))
	activations[0] = X
	deltas := make([]blas64.General, 0, len(layerUnits)-1)
	off := 0
	for _, nFanOut := range layerUnits[1:] {
		size := mlp.BatchSize * nFanOut
		off += size + size
	}
	mem := make([]float64, off)
	off = 0
	for _, nFanOut := range layerUnits[1:] {
		size := mlp.BatchSize * nFanOut
		activations = append(activations, blas64General{Rows: mlp.BatchSize, Cols: nFanOut, Stride: nFanOut, Data: mem[off : off+size]})
		off += size
		deltas = append(deltas, blas64General{Rows: mlp.BatchSize, Cols: nFanOut, Stride: nFanOut, Data: mem[off : off+size]})
		off += size
	}
	off = len(mlp.packedParameters)
	packedGrads := make([]float64, off)
	CoefsGrads := make([]blas64General, mlp.NLayers-1)
	InterceptsGrads := make([][]float64, mlp.NLayers-1)
	off = 0
	for i := 0; i < mlp.NLayers-1; i++ {
		InterceptsGrads[i] = packedGrads[off : off+layerUnits[i+1]]
		off += layerUnits[i+1]
		CoefsGrads[i] = blas64General{Rows: layerUnits[i], Cols: layerUnits[i+1], Stride: layerUnits[i+1], Data: packedGrads[off : off+layerUnits[i]*layerUnits[i+1]]}
		off += layerUnits[i] * layerUnits[i+1]
	}
	if strings.EqualFold(mlp.Solver, "lbfgs") {
		mlp.fitLbfgs(X, y, activations, deltas, CoefsGrads,
			InterceptsGrads, packedGrads, layerUnits)
	} else {
		mlp.fitStochastic(X, y, activations, deltas, CoefsGrads,
			InterceptsGrads, packedGrads, layerUnits, incremental)
	}
	mlp.packedGrads = packedGrads
}
func (mlp *BaseMultilayerPerceptron64) IsClassifier() bool {
	return mlp.LossFuncName != "square_loss"
}
func (mlp *BaseMultilayerPerceptron64) Fit(X, Y Matrix) {
	var xb, yb blas64.General
	if xg, ok := X.(RawMatrixer64); ok && !mlp.Shuffle {
		if yg, ok := Y.(RawMatrixer64); ok {
			xb, yb = xg.RawMatrix(), yg.RawMatrix()
		}
	} else {
		var tmp General64
		tmp = General64(xb)
		tmp.Copy(X)
		xb = tmp.RawMatrix()
		tmp = General64(yb)
		tmp.Copy(Y)
		yb = tmp.RawMatrix()
	}
	if mlp.IsClassifier() && !isBinarized64(yb) {
		mlp.lb = NewLabelBinarizer64(0, 1)
		xbin, ybin := mlp.lb.FitTransform(General64(xb), General64(yb))
		xb, yb = blas64.General(xbin), blas64.General(ybin)
	}
	mlp.fit(xb, yb, false)
}
func (mlp *BaseMultilayerPerceptron64) GetNOutputs() int {
	if mlp.lb != nil {
		return len(mlp.lb.Classes)
	}
	return mlp.NOutputs
}