func RegularizerTest(t *testing.T, r Regularizer, name string, parameters []float64, trueLoss float64, trueDeriv []float64) {
	// Test that Loss works
	loss := r.Loss(parameters)
	if math.Abs(loss-trueLoss) > 1e-14 {
		t.Errorf("Loss doesn't match for case %v. Expected: %v, Found: %v", name, trueLoss, loss)
	}
	// Test that LossDeriv works
	derivative := make([]float64, len(trueDeriv))
	lossDeriv := r.LossDeriv(parameters, derivative)

	if math.Abs(lossDeriv-trueLoss) > 1e-14 {
		t.Errorf("Loss doesn't match from LossDeriv for case %v. Expected: %v, Found: %v", name, trueLoss, lossDeriv)
	}
	if !floats.EqualApprox(trueDeriv, derivative, 1e-14) {
		t.Errorf("Derivative doesn't match from LossDeriv for case %v", name)
	}

	for i := range derivative {
		derivative[i] = float64(i)
	}

	lossAddDeriv := r.LossAddDeriv(parameters, derivative)
	if math.Abs(lossAddDeriv-trueLoss) > 1e-14 {
		t.Errorf("Loss doesn't match from LossAddDeriv for case %v. Expected: %v, Found: %v", name, trueLoss, lossAddDeriv)
	}
	for i := range derivative {
		derivative[i] -= float64(i)
	}
	if !floats.EqualApprox(trueDeriv, derivative, 1e-14) {
		t.Errorf("Derivative doesn't match from LossAddDeriv for case %v", name)
	}
}

func (s *S) TestVectorMul(c *check.C) {

	for i, test := range []struct {
		m int
		n int
	}{
		{
			m: 10,
			n: 5,
		},
		{
			m: 5,
			n: 5,
		},
		{
			m: 5,
			n: 10,
		},
	} {
		vData := make([]float64, test.n)
		for i := range vData {
			vData[i] = rand.Float64()
		}
		vDataCopy := make([]float64, test.n)
		copy(vDataCopy, vData)
		v := NewVector(test.n, vData)
		aData := make([]float64, test.n*test.m)
		for i := range aData {
			aData[i] = rand.Float64()
		}
		a := NewDense(test.m, test.n, aData)
		var v2 Vector
		v2.MulVec(a, false, v)
		var v2M Dense
		v2M.Mul(a, v)
		same := floats.EqualApprox(v2.mat.Data, v2M.mat.Data, 1e-14)
		c.Check(same, check.Equals, true, check.Commentf("Test %d", i))

		var aT Dense
		aT.TCopy(a)
		v2.MulVec(&aT, true, v)
		same = floats.EqualApprox(v2.mat.Data, v2M.mat.Data, 1e-14)
		c.Check(same, check.Equals, true, check.Commentf("Test %d", i))

		/*
			v.MulVec(&aT, true, v)
			same = floats.EqualApprox(v.mat.Data, v2M.mat.Data, 1e-14)
			c.Check(same, check.Equals, true, check.Commentf("Test %d", i))
		*/
	}
}

func TestLogSquared(t *testing.T) {
	prediction := []float64{1, -2, 3}
	truth := []float64{1.1, -2.2, 2.7}
	trueloss := (math.Log(.1*.1+1) + math.Log(.2*.2+1) + math.Log(.3*.3+1)) / 3
	derivative := []float64{0, 0, 0}

	sq := LogSquared{}
	loss := sq.Loss(prediction, truth)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from Loss(). Expected %v, Found: %v", trueloss, loss)
	}

	loss = sq.LossDeriv(prediction, truth, derivative)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from LossDeriv()")
	}
	derivative, fdDerivative := finiteDifferenceLosser(sq, prediction, truth)
	if !floats.EqualApprox(derivative, fdDerivative, FDTol) {
		t.Errorf("Derivative doesn't match. \n deriv: %v \n fdDeriv: %v ", derivative, fdDerivative)
	}
	err := common.InterfaceTestMarshalAndUnmarshal(sq)
	if err != nil {
		t.Errorf("Error marshaling and unmarshaling")
	}
}

func TestRelativeLog(t *testing.T) {
	tol := 1e-2
	prediction := []float64{1, -2, 3}
	truth := []float64{1.1, -2.2, 2.7}
	trueloss := ((.1/(1.1+tol))*(.1/(1.1+tol)) + (.2/(2.2+tol))*(.2/(2.2+tol)) + (.3/(2.7+tol))*(.3/(2.7+tol))) / 3
	trueloss = math.Log(trueloss + 1)
	derivative := []float64{0, 0, 0}

	sq := RelativeLog(tol)
	loss := sq.Loss(prediction, truth)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from Loss(). Expected %v, Found: %v", trueloss, loss)
	}

	loss = sq.LossDeriv(prediction, truth, derivative)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from LossDeriv()")
	}
	derivative, fdDerivative := finiteDifferenceLosser(sq, prediction, truth)
	if !floats.EqualApprox(derivative, fdDerivative, FDTol) {
		t.Errorf("Derivative doesn't match. \n deriv: %v \n fdDeriv: %v ", derivative, fdDerivative)
	}

	err := common.InterfaceTestMarshalAndUnmarshal(sq)
	if err != nil {
		t.Errorf("Error marshaling and unmarshaling: " + err.Error())
	}
}

func DrsclTest(t *testing.T, impl Drscler) {
	for _, test := range []struct {
		x []float64
		a float64
	}{
		{
			x: []float64{1, 2, 3, 4, 5},
			a: 4,
		},
		{
			x: []float64{1, 2, 3, 4, 5},
			a: math.MaxFloat64,
		},
		{
			x: []float64{1, 2, 3, 4, 5},
			a: 1e-307,
		},
	} {
		xcopy := make([]float64, len(test.x))
		copy(xcopy, test.x)

		// Cannot test the scaling directly because of floating point scaling issues
		// (the purpose of Drscl). Instead, check that scaling and scaling back
		// yeilds approximately x. If overflow or underflow occurs then the scaling
		// won't match.
		impl.Drscl(len(test.x), test.a, xcopy, 1)
		if floats.Equal(xcopy, test.x) {
			t.Errorf("x unchanged during call to drscl. a = %v, x = %v.", test.a, test.x)
		}
		impl.Drscl(len(test.x), 1/test.a, xcopy, 1)
		if !floats.EqualApprox(xcopy, test.x, 1e-14) {
			t.Errorf("x not equal after scaling and unscaling. a = %v, x = %v.", test.a, test.x)
		}
	}
}

func testDerivParam(t *testing.T, d derivParamTester) {
	// Tests that the derivative matches for a number of different quantiles
	// along the distribution.
	nTest := 10
	quantiles := make([]float64, nTest)
	floats.Span(quantiles, 0.1, 0.9)

	deriv := make([]float64, d.NumParameters())
	fdDeriv := make([]float64, d.NumParameters())

	initParams := d.parameters(nil)
	init := make([]float64, d.NumParameters())
	for i, v := range initParams {
		init[i] = v.Value
	}
	for _, v := range quantiles {
		d.setParameters(initParams)
		x := d.Quantile(v)
		d.DLogProbDParam(x, deriv)
		f := func(p []float64) float64 {
			params := d.parameters(nil)
			for i, v := range p {
				params[i].Value = v
			}
			d.setParameters(params)
			return d.LogProb(x)
		}
		fd.Gradient(fdDeriv, f, init, nil)
		if !floats.EqualApprox(deriv, fdDeriv, 1e-6) {
			t.Fatal("Derivative mismatch. Want", fdDeriv, ", got", deriv, ".")
		}
	}
}

func TestPrivatePredictsMatch(t *testing.T) {
	for i, test := range netIniters {
		for j := 0; j < nRandSamp; j++ {
			n := testNets[i]
			input := make([]float64, test.inputDim)
			floats.Fill(rand.NormFloat64, input)
			outputSimple := make([]float64, test.outputDim)
			floats.Fill(rand.NormFloat64, outputSimple)
			outputCache := make([]float64, test.outputDim)
			floats.Fill(rand.NormFloat64, outputCache)

			// predict using uncached method
			tmp1, tmp2 := newPredictMemory(n.neurons)
			predict(input, n.neurons, n.parameters, tmp1, tmp2, outputSimple)

			// predict using cached method
			combinations := newPerNeuronMemory(n.neurons)
			outputs := newPerNeuronMemory(n.neurons)
			cachePredict(input, n.neurons, n.parameters, combinations, outputs, outputCache)

			if !floats.EqualApprox(outputSimple, outputCache, 1e-14) {
				t.Errorf("test %v: output mismatch between simple and cached predict. Simple: %v, Cached: %v", test.name, outputSimple, outputCache)
				break
			}
		}
	}
}

func DspmvTest(t *testing.T, blasser Dspmver) {
	for i, test := range []struct {
		ul    blas.Uplo
		n     int
		a     [][]float64
		x     []float64
		y     []float64
		alpha float64
		beta  float64
		ans   []float64
	}{
		{
			ul: blas.Upper,
			n:  3,
			a: [][]float64{
				{5, 6, 7},
				{0, 8, 10},
				{0, 0, 13},
			},
			x:     []float64{3, 4, 5},
			y:     []float64{6, 7, 8},
			alpha: 2.1,
			beta:  -3,
			ans:   []float64{137.4, 189, 240.6},
		},
		{
			ul: blas.Lower,
			n:  3,
			a: [][]float64{
				{5, 0, 0},
				{6, 8, 0},
				{7, 10, 13},
			},
			x:     []float64{3, 4, 5},
			y:     []float64{6, 7, 8},
			alpha: 2.1,
			beta:  -3,
			ans:   []float64{137.4, 189, 240.6},
		},
	} {
		incTest := func(incX, incY, extra int) {
			x := makeIncremented(test.x, incX, extra)
			y := makeIncremented(test.y, incY, extra)
			aFlat := flattenTriangular(test.a, test.ul)
			ans := makeIncremented(test.ans, incY, extra)

			blasser.Dspmv(test.ul, test.n, test.alpha, aFlat, x, incX, test.beta, y, incY)
			if !floats.EqualApprox(ans, y, 1e-14) {
				t.Errorf("Case %v, incX=%v, incY=%v: Want %v, got %v.", i, incX, incY, ans, y)
			}
		}
		incTest(1, 1, 0)
		incTest(2, 3, 0)
		incTest(3, 2, 0)
		incTest(-3, 2, 0)
		incTest(-2, 4, 0)
		incTest(2, -1, 0)
		incTest(-3, -4, 3)
	}
}

func Dlasv2Test(t *testing.T, impl Dlasv2er) {
	rnd := rand.New(rand.NewSource(1))
	for i := 0; i < 100; i++ {
		f := rnd.NormFloat64()
		g := rnd.NormFloat64()
		h := rnd.NormFloat64()

		ssmin, ssmax, snr, csr, snl, csl := impl.Dlasv2(f, g, h)

		// tmp =
		// [ csl snl] [f g]
		// [-snl csl] [0 h]
		tmp11 := csl * f
		tmp12 := csl*g + snl*h
		tmp21 := -snl * f
		tmp22 := -snl*g + csl*h
		// lhs =
		// [tmp11 tmp12] [csr -snr]
		// [tmp21 tmp22] [snr  csr]
		ans11 := tmp11*csr + tmp12*snr
		ans12 := tmp11*-snr + tmp12*csr
		ans21 := tmp21*csr + tmp22*snr
		ans22 := tmp21*-snr + tmp22*csr

		lhs := []float64{ans11, ans12, ans21, ans22}
		rhs := []float64{ssmax, 0, 0, ssmin}
		if !floats.EqualApprox(rhs, lhs, 1e-12) {
			t.Errorf("SVD mismatch. f = %v, g = %v, h = %v.\nLHS: %v\nRHS: %v", f, g, h, lhs, rhs)
		}
	}
}

func testDpotf2(t *testing.T, impl Dpotf2er, testPos bool, a, ans [][]float64, stride int, ul blas.Uplo) {
	aFlat := flattenTri(a, stride, ul)
	ansFlat := flattenTri(ans, stride, ul)
	pos := impl.Dpotf2(ul, len(a[0]), aFlat, stride)
	if pos != testPos {
		t.Errorf("Positive definite mismatch: Want %v, Got %v", testPos, pos)
		return
	}
	if testPos && !floats.EqualApprox(ansFlat, aFlat, 1e-14) {
		t.Errorf("Result mismatch: Want %v, Got  %v", ansFlat, aFlat)
	}
}

func TestNormRand(t *testing.T) {
	for _, test := range []struct {
		mean []float64
		cov  []float64
	}{
		{
			mean: []float64{0, 0},
			cov: []float64{
				1, 0,
				0, 1,
			},
		},
		{
			mean: []float64{0, 0},
			cov: []float64{
				1, 0.9,
				0.9, 1,
			},
		},
		{
			mean: []float64{6, 7},
			cov: []float64{
				5, 0.9,
				0.9, 2,
			},
		},
	} {
		dim := len(test.mean)
		cov := mat64.NewSymDense(dim, test.cov)
		n, ok := NewNormal(test.mean, cov, nil)
		if !ok {
			t.Errorf("bad covariance matrix")
		}

		nSamples := 1000000
		samps := mat64.NewDense(nSamples, dim, nil)
		for i := 0; i < nSamples; i++ {
			n.Rand(samps.RawRowView(i))
		}
		estMean := make([]float64, dim)
		for i := range estMean {
			estMean[i] = stat.Mean(mat64.Col(nil, i, samps), nil)
		}
		if !floats.EqualApprox(estMean, test.mean, 1e-2) {
			t.Errorf("Mean mismatch: want: %v, got %v", test.mean, estMean)
		}
		estCov := stat.CovarianceMatrix(nil, samps, nil)
		if !mat64.EqualApprox(estCov, cov, 1e-2) {
			t.Errorf("Cov mismatch: want: %v, got %v", cov, estCov)
		}
	}
}

func testLinear(t *testing.T, kind linearTest) {
	u := &Linear{}

	data := flatten(kind.data)
	err := u.SetScale(data)

	if err != nil {
		if kind.eqDim != true {
			t.Errorf("Error where there shouldn't be for case " + kind.name + ": " + err.Error())
		}
	}
	if !floats.EqualApprox(u.Min, kind.min, 1e-14) {
		t.Errorf("Min doesn't match for case " + kind.name)
	}
	if !floats.EqualApprox(u.Max, kind.max, 1e-14) {
		t.Errorf("Max doesn't match for case " + kind.name)
	}
	scaledData := flatten(kind.scaledData)
	testScaling(t, u, data, scaledData, kind.name)
	u2 := &Linear{}
	testGob(u, u2, t)
}

func denseEqualApprox(a *Dense, acomp matComp, tol float64) bool {
	ar2, ac2 := a.Dims()
	if ar2 != acomp.r {
		return false
	}
	if ac2 != acomp.c {
		return false
	}
	if !floats.EqualApprox(a.mat.Data, acomp.data, tol) {
		return false
	}
	return true
}

func testNormal(t *testing.T, kind normalTest) {
	u := &Normal{}
	data := flatten(kind.data)
	err := u.SetScale(data)

	if err != nil {
		if kind.eqDim != true {
			t.Errorf("Error where there shouldn't be for case " + kind.name + ": " + err.Error())
		}
	}
	if !floats.EqualApprox(u.Mu, kind.mu, 1e-14) {
		t.Errorf("Mu doesn't match for case "+kind.name+". Expected: %v, Found: %v", kind.mu, u.Mu)
	}
	if !floats.EqualApprox(u.Sigma, kind.sigma, 1e-14) {
		t.Errorf("Sigma doesn't match for case "+kind.name+". Expected: %v, Found: %v", kind.sigma, u.Sigma)
	}
	scaledData := flatten(kind.scaledData)
	testScaling(t, u, data, scaledData, kind.name)

	u2 := &Normal{}
	testGob(u, u2, t)
}

func DorgqlTest(t *testing.T, impl Dorgqler) {
	rnd := rand.New(rand.NewSource(1))
	for _, test := range []struct {
		m, n, k, lda int
	}{
		{5, 4, 3, 0},
		{100, 100, 100, 0},
		{200, 100, 50, 0},
		{200, 200, 50, 0},
	} {
		m := test.m
		n := test.n
		k := test.k
		lda := test.lda
		if lda == 0 {
			lda = n
		}
		a := make([]float64, m*lda)
		for i := range a {
			a[i] = rnd.NormFloat64()
		}
		tau := nanSlice(min(m, n))
		work := nanSlice(max(m, n))

		impl.Dgeql2(m, n, a, lda, tau, work)

		aCopy := make([]float64, len(a))
		copy(aCopy, a)

		impl.Dorg2l(m, n, k, a, lda, tau, work)
		ans := make([]float64, len(a))
		copy(ans, a)

		impl.Dorgql(m, n, k, a, lda, tau, work, -1)
		work = make([]float64, int(work[0]))
		copy(a, aCopy)
		impl.Dorgql(m, n, k, a, lda, tau, work, len(work))

		if !floats.EqualApprox(a, ans, 1e-8) {
			t.Errorf("Answer mismatch. m = %v, n = %v, k = %v", m, n, k)
		}
	}
}

// TestDeriv uses finite difference to test that the prediction from Deriv
// is correct, and tests that computing the loss in parallel works properly
// Only does finite difference for the first nTest to save time
func TestDeriv(t *testing.T, trainable DerivTester, inputs, trueOutputs common.RowMatrix, name string) {

	// Set the parameters to something random
	trainable.RandomizeParameters()

	// Compute the loss and derivative
	losser := loss.SquaredDistance{}
	regularizer := regularize.TwoNorm{}

	batchGrad := train.NewBatchGradBased(trainable, true, inputs, trueOutputs, losser, regularizer)

	derivative := make([]float64, trainable.NumParameters())
	parameters := trainable.Parameters(nil)
	// Don't need to check loss, because if predict is right and losser is right then loss must be correct
	_ = batchGrad.ObjGrad(parameters, derivative)

	fdDerivative := make([]float64, trainable.NumParameters())

	wg := &sync.WaitGroup{}
	wg.Add(trainable.NumParameters())
	for i := 0; i < trainable.NumParameters(); i++ {
		go func(i int) {
			newParameters := make([]float64, trainable.NumParameters())
			tmpDerivative := make([]float64, trainable.NumParameters())
			copy(newParameters, parameters)
			newParameters[i] += fdStep
			loss1 := batchGrad.ObjGrad(newParameters, tmpDerivative)
			newParameters[i] -= 2 * fdStep
			loss2 := batchGrad.ObjGrad(newParameters, tmpDerivative)
			newParameters[i] += fdStep
			fdDerivative[i] = (loss1 - loss2) / (2 * fdStep)
			wg.Done()
		}(i)
	}
	wg.Wait()
	if !floats.EqualApprox(derivative, fdDerivative, 1e-6) {
		t.Errorf("%v: deriv doesn't match: Finite Difference: %v, Analytic: %v", name, fdDerivative, derivative)
	}

}

func TestPredictFeaturized(t *testing.T) {
	for _, test := range []struct {
		z              []float64
		featureWeights [][]float64
		output         []float64
		Name           string
	}{
		{
			Name: "General",
			z:    []float64{1, 2, 3},
			featureWeights: [][]float64{
				{3, 4},
				{1, 2},
				{0.5, 0.4},
			},
			output: []float64{6.5, 9.2},
		},
	} {
		zCopy := make([]float64, len(test.z))
		copy(zCopy, test.z)
		fwMat := flatten(test.featureWeights)
		fwMatCopy := &mat64.Dense{}
		fwMatCopy.Clone(fwMat)

		output := make([]float64, len(test.output))

		predictFeaturized(zCopy, fwMat, output)

		// Test that z wasn't changed
		if !floats.Equal(test.z, zCopy) {
			t.Errorf("z changed during call")
		}

		if !floats.EqualApprox(output, test.output, 1e-14) {
			t.Errorf("output doesn't match for test %v. Expected %v, found %v", test.Name, test.output, output)
		}
	}
}

func TestManhattanDistance(t *testing.T) {
	prediction := []float64{1, 2, 3}
	truth := []float64{1.1, 2.2, 2.7}
	trueloss := (.1 + .2 + .3) / 3
	derivative := []float64{0, 0, 0}

	sq := ManhattanDistance{}
	loss := sq.Loss(prediction, truth)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from Loss()")
	}

	loss = sq.LossDeriv(prediction, truth, derivative)
	if math.Abs(loss-trueloss) > TOL {
		t.Errorf("loss doesn't match from LossDeriv()")
	}
	derivative, fdDerivative := finiteDifferenceLosser(sq, prediction, truth)
	if !floats.EqualApprox(derivative, fdDerivative, FDTol) {
		t.Errorf("Derivative doesn't match. \n deriv: %v \n fdDeriv: %v ", derivative, fdDerivative)
	}

	err := common.InterfaceTestMarshalAndUnmarshal(sq)
	if err != nil {
		t.Errorf("Error marshaling and unmarshaling")
	}

	truth = []float64{1, 2, 3}
	loss = sq.LossDeriv(prediction, truth, derivative)
	if loss != 0 {
		t.Errorf("Non-zero loss for equal pred and truth")
	}
	for _, val := range derivative {
		if val != 0 {
			t.Errorf("Non-zero derivative for equal pred and truth")
		}
	}
}

//.........这里部分代码省略.........
					k := test.k
					ldc := test.ldc
					if ldc == 0 {
						ldc = n
					}
					nq := n
					if side == blas.Left {
						nq = m
					}

					// Compute a decomposition.
					var ma, na int
					var a []float64
					if vect == lapack.ApplyQ {
						ma = nq
						na = k
					} else {
						ma = k
						na = nq
					}
					lda := test.lda
					if lda == 0 {
						lda = na
					}
					a = make([]float64, ma*lda)
					for i := range a {
						a[i] = rnd.NormFloat64()
					}
					nTau := min(nq, k)
					tauP := make([]float64, nTau)
					tauQ := make([]float64, nTau)
					d := make([]float64, nTau)
					e := make([]float64, nTau)
					lwork := -1
					work := make([]float64, 1)
					impl.Dgebrd(ma, na, a, lda, d, e, tauQ, tauP, work, lwork)
					work = make([]float64, int(work[0]))
					lwork = len(work)

					impl.Dgebrd(ma, na, a, lda, d, e, tauQ, tauP, work, lwork)

					// Apply and compare update.
					c := make([]float64, m*ldc)
					for i := range c {
						c[i] = rnd.NormFloat64()
					}

					cCopy := make([]float64, len(c))
					copy(cCopy, c)

					if vect == lapack.ApplyQ {
						impl.Dormbr(vect, side, trans, m, n, k, a, lda, tauQ, c, ldc, work, lwork)
					} else {
						impl.Dormbr(vect, side, trans, m, n, k, a, lda, tauP, c, ldc, work, lwork)
					}

					// Check that the multiplication was correct.
					cOrig := blas64.General{
						Rows:   m,
						Cols:   n,
						Stride: ldc,
						Data:   make([]float64, len(cCopy)),
					}
					copy(cOrig.Data, cCopy)
					cAns := blas64.General{
						Rows:   m,
						Cols:   n,
						Stride: ldc,
						Data:   make([]float64, len(cCopy)),
					}
					copy(cAns.Data, cCopy)
					nb := min(ma, na)
					var mulMat blas64.General
					if vect == lapack.ApplyQ {
						mulMat = constructQPBidiagonal(lapack.ApplyQ, ma, na, nb, a, lda, tauQ)
					} else {
						mulMat = constructQPBidiagonal(lapack.ApplyP, ma, na, nb, a, lda, tauP)
					}

					mulTrans := trans

					if side == blas.Left {
						bi.Dgemm(mulTrans, blas.NoTrans, m, n, m, 1, mulMat.Data, mulMat.Stride, cOrig.Data, cOrig.Stride, 0, cAns.Data, cAns.Stride)
					} else {
						bi.Dgemm(blas.NoTrans, mulTrans, m, n, n, 1, cOrig.Data, cOrig.Stride, mulMat.Data, mulMat.Stride, 0, cAns.Data, cAns.Stride)
					}

					if !floats.EqualApprox(cAns.Data, c, 1e-8) {
						isApplyQ := vect == lapack.ApplyQ
						isLeft := side == blas.Left
						isTrans := trans == blas.Trans

						t.Errorf("C mismatch. isApplyQ: %v, isLeft: %v, isTrans: %v, m = %v, n = %v, k = %v, lda = %v, ldc = %v",
							isApplyQ, isLeft, isTrans, m, n, k, lda, ldc)
					}
				}
			}
		}
	}
}

// testFunction checks that the function can evaluate itself (and its gradient)
// correctly.
func testFunction(f function, ftests []funcTest, t *testing.T) {
	// Make a copy of tests because we may append to the slice.
	tests := make([]funcTest, len(ftests))
	copy(tests, ftests)

	// Get information about the function.
	fMinima, isMinimumer := f.(minimumer)
	fGradient, isGradient := f.(gradient)

	// If the function is a Minimumer, append its minima to the tests.
	if isMinimumer {
		for _, minimum := range fMinima.Minima() {
			// Allocate gradient only if the function can evaluate it.
			var grad []float64
			if isGradient {
				grad = make([]float64, len(minimum.X))
			}
			tests = append(tests, funcTest{
				X:        minimum.X,
				F:        minimum.F,
				Gradient: grad,
			})
		}
	}

	for i, test := range tests {
		F := f.Func(test.X)

		// Check that the function value is as expected.
		if math.Abs(F-test.F) > defaultTol {
			t.Errorf("Test #%d: function value given by Func is incorrect. Want: %v, Got: %v",
				i, test.F, F)
		}

		if test.Gradient == nil {
			continue
		}

		// Evaluate the finite difference gradient.
		fdGrad := fd.Gradient(nil, f.Func, test.X, nil)

		// Check that the finite difference and expected gradients match.
		if !floats.EqualApprox(fdGrad, test.Gradient, defaultFDGradTol) {
			dist := floats.Distance(fdGrad, test.Gradient, math.Inf(1))
			t.Errorf("Test #%d: numerical and expected gradients do not match. |fdGrad - WantGrad|_∞ = %v",
				i, dist)
		}

		// If the function is a Gradient, check that it computes the gradient correctly.
		if isGradient {
			grad := make([]float64, len(test.Gradient))
			fGradient.Grad(grad, test.X)

			if !floats.EqualApprox(grad, test.Gradient, defaultGradTol) {
				dist := floats.Distance(grad, test.Gradient, math.Inf(1))
				t.Errorf("Test #%d: gradient given by Grad is incorrect. |grad - WantGrad|_∞ = %v",
					i, dist)
			}
		}
	}
}