func (s *S) TestSetRowColumn(c *check.C) {
	for _, as := range [][][]float64{
		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
		{{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}},
	} {
		for ri, row := range as {
			a := NewDense(flatten(as))
			t := &Dense{}
			t.Clone(a)
			a.SetRow(ri, make([]float64, a.mat.Cols))
			t.Sub(t, a)
			c.Check(t.Norm(0), check.Equals, floats.Norm(row, 2))
		}

		for ci := range as[0] {
			a := NewDense(flatten(as))
			t := &Dense{}
			t.Clone(a)
			a.SetCol(ci, make([]float64, a.mat.Rows))
			col := make([]float64, a.mat.Rows)
			for j := range col {
				col[j] = float64(ci + 1 + j*a.mat.Cols)
			}
			t.Sub(t, a)
			c.Check(t.Norm(0), check.Equals, floats.Norm(col, 2))
		}
	}
}

func (b *BFGS) InitDirection(loc *Location, dir []float64) (stepSize float64) {
	dim := len(loc.X)
	b.dim = dim

	b.x = resize(b.x, dim)
	copy(b.x, loc.X)
	b.grad = resize(b.grad, dim)
	copy(b.grad, loc.Gradient)

	b.y = resize(b.y, dim)
	b.s = resize(b.s, dim)
	b.tmp = resize(b.tmp, dim)
	b.yVec = mat64.NewVector(dim, b.y)
	b.sVec = mat64.NewVector(dim, b.s)
	b.tmpVec = mat64.NewVector(dim, b.tmp)

	if b.invHess == nil || cap(b.invHess.RawSymmetric().Data) < dim*dim {
		b.invHess = mat64.NewSymDense(dim, nil)
	} else {
		b.invHess = mat64.NewSymDense(dim, b.invHess.RawSymmetric().Data[:dim*dim])
	}

	// The values of the hessian are initialized in the first call to NextDirection

	// initial direcion is just negative of gradient because the hessian is 1
	copy(dir, loc.Gradient)
	floats.Scale(-1, dir)

	b.first = true

	return 1 / floats.Norm(dir, 2)
}

// SetCurr sets the current value of the float
// Assumes that the length does not change per iteration.
func (f *Floats) SetCurrent(val []float64) {

	copy(f.previous, f.current)
	copy(f.current, val)
	floats.SubTo(f.diff, f.current, f.previous)
	f.norm = floats.Norm(f.current, 2)
}

func TestMinimalSurface(t *testing.T) {
	for _, size := range [][2]int{
		{20, 30},
		{30, 30},
		{50, 40},
	} {
		f := NewMinimalSurface(size[0], size[1])
		x0 := f.InitX()
		grad := make([]float64, len(x0))
		f.Grad(grad, x0)
		fdGrad := fd.Gradient(nil, f.Func, x0, &fd.Settings{Formula: fd.Central})

		// Test that the numerical and analytical gradients agree.
		dist := floats.Distance(grad, fdGrad, math.Inf(1))
		if dist > 1e-9 {
			t.Errorf("grid %v x %v: numerical and analytical gradient do not match. |fdGrad - grad|_∞ = %v",
				size[0], size[1], dist)
		}

		// Test that the gradient at the minimum is small enough.
		// In some sense this test is not completely correct because ExactX
		// returns the exact solution to the continuous problem projected on the
		// grid, not the exact solution to the discrete problem which we are
		// solving. This is the reason why a relatively loose tolerance 1e-4
		// must be used.
		xSol := f.ExactX()
		f.Grad(grad, xSol)
		norm := floats.Norm(grad, math.Inf(1))
		if norm > 1e-4 {
			t.Errorf("grid %v x %v: gradient at the minimum not small enough. |grad|_∞ = %v",
				size[0], size[1], norm)
		}
	}
}

func Solve(a sparse.Matrix, b, xInit []float64, settings *Settings, method Method) (result Result, err error) {
	stats := Stats{
		StartTime: time.Now(),
	}

	dim := len(xInit)
	if dim == 0 {
		panic("iterative: invalid dimension")
	}

	r, c := a.Dims()
	if r != c {
		panic("iterative: matrix is not square")
	}
	if c != dim {
		panic("iterative: mismatched size of the matrix")
	}
	if len(b) != dim {
		panic("iterative: mismatched size of the right-hand side vector")
	}

	if settings == nil {
		settings = DefaultSettings(dim)
	}

	ctx := Context{
		X:        make([]float64, dim),
		Residual: make([]float64, dim),
	}
	copy(ctx.X, xInit)
	copy(ctx.Residual, b)
	if floats.Norm(ctx.X, math.Inf(1)) > 0 {
		sparse.MulMatVec(-1, false, a, ctx.X, 1, 1, ctx.Residual, 1)
		stats.MatVecMultiplies++
	}

	if floats.Norm(ctx.Residual, 2) >= settings.Tolerance {
		err = iterate(method, a, b, settings, &ctx, &stats)
	}

	result = Result{
		X:       ctx.X,
		Stats:   stats,
		Runtime: time.Since(stats.StartTime),
	}
	return result, err
}

func checkConvergence(loc *Location, iterType IterationType, stats *Stats, settings *Settings) Status {
	if iterType == MajorIteration || iterType == InitIteration {
		if loc.Gradient != nil {
			norm := floats.Norm(loc.Gradient, math.Inf(1))
			if norm < settings.GradientThreshold {
				return GradientThreshold
			}
		}
		if loc.F < settings.FunctionThreshold {
			return FunctionThreshold
		}
	}

	if iterType == MajorIteration && settings.FunctionConverge != nil {
		status := settings.FunctionConverge.FunctionConverged(loc.F)
		if status != NotTerminated {
			return status
		}
	}

	// Check every step for negative infinity because it could break the
	// linesearches and -inf is the best you can do anyway.
	if math.IsInf(loc.F, -1) {
		return FunctionNegativeInfinity
	}

	if settings.FuncEvaluations > 0 {
		if stats.FuncEvaluations >= settings.FuncEvaluations {
			return FunctionEvaluationLimit
		}
	}

	if settings.GradEvaluations > 0 {
		if stats.GradEvaluations >= settings.GradEvaluations {
			return GradientEvaluationLimit
		}
	}

	if settings.HessEvaluations > 0 {
		if stats.HessEvaluations >= settings.HessEvaluations {
			return HessianEvaluationLimit
		}
	}

	if settings.Runtime > 0 {
		// TODO(vladimir-ch): It would be nice to update Runtime here.
		if stats.Runtime >= settings.Runtime {
			return RuntimeLimit
		}
	}

	if iterType == MajorIteration && settings.MajorIterations > 0 {
		if stats.MajorIterations >= settings.MajorIterations {
			return IterationLimit
		}
	}
	return NotTerminated
}

// Init sets the initial value of the variable to be
// used by the optimizer. (for example, initial location, initial
// function value, etc.)
func (f *Floats) SetInit(val []float64) {
	if len(f.init) > len(val) {
		f.init = f.init[:len(val)]
	} else {
		f.init = make([]float64, len(val))
	}
	copy(f.init, val)
	f.norm = floats.Norm(val, 2)
}

// SetCurr sets a new value  for the current location and updates the
// delta from the last value
func (o *Objective) SetCurr(f []float64) {
	// Find the current delta in the values
	o.delta = math.Abs(o.Floats.norm - floats.Norm(f, 2))
	// Set the initial delta if necessary
	if math.IsNaN(o.initDelta) {
		o.initDelta = o.delta
	}
	// Set the new current value
	o.Floats.SetCurr(f)
}

func (o OneNorm) LossDeriv(parameters, derivative []float64) float64 {
	loss := o.Gamma * floats.Norm(parameters, 2)
	for i, p := range parameters {
		derivative[i] = o.Gamma
		if p < 0 {
			derivative[i] = -o.Gamma
		}
	}
	return loss
}

func (o OneNorm) LossAddDeriv(parameters, derivative []float64) float64 {
	loss := o.Gamma * floats.Norm(parameters, 2)
	for i, p := range parameters {
		add := o.Gamma
		if p < 0 {
			add = -o.Gamma
		}
		derivative[i] += add
	}
	return loss
}

// Status checks if either of the tolerances have converged
func (f *Floats) Status() common.Status {
	s := f.AbsTol.Status(f.norm)
	if s != common.Continue {
		return s
	}
	s = f.ChangeTol.Status(floats.Norm(f.diff, 2))
	if s != common.Continue {
		return s
	}
	return s
}

// Initialize initializes the Float to be ready to optimize by
// setting the history slice to have length zero, and setting
// the current value equal to the initial value
// This should be called by the optimizer at the beginning of
// the optimization
func (f *Floats) Initialize() error {
	f.hist = f.hist[:0]

	if f.init == nil {
		return errors.New("multivariate: inial slice is nil")
	}

	f.curr = make([]float64, len(f.init))
	copy(f.curr, f.init)
	f.opt = nil
	f.normInit = floats.Norm(f.curr, 2)
	return nil
}

func iterate(method Method, a sparse.Matrix, b []float64, settings *Settings, ctx *Context, stats *Stats) error {
	dim := len(ctx.X)
	bNorm := floats.Norm(b, 2)
	if bNorm == 0 {
		bNorm = 1
	}

	request := method.Init(ctx)
	for {
		switch request {
		case NoOperation:

		case ComputeAp:
			ctx.Ap = resize(ctx.Ap, dim)
			sparse.MulMatVec(1, false, a, ctx.P, 1, 0, ctx.Ap, 1)
			stats.MatVecMultiplies++

		case SolvePreconditioner:
			ctx.Z = resize(ctx.Z, dim)
			copy(ctx.Z, ctx.Residual)
			stats.PrecondionerSolves++

		case CheckConvergence:
			stats.Iterations++
			stats.Residual = floats.Norm(ctx.Residual, 2) / bNorm
			fmt.Println(stats.Residual)
			if stats.Residual < settings.Tolerance {
				return nil
			}
			if stats.Iterations == settings.Iterations {
				return errors.New("iterative: reached iteration limit")
			}
		}

		request = method.Iterate(ctx)
	}
}

func TestSetRowColumn(t *testing.T) {
	for _, as := range [][][]float64{
		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}},
		{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
		{{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}},
	} {
		for ri, row := range as {
			a := NewDense(flatten(as))
			m := &Dense{}
			m.Clone(a)
			a.SetRow(ri, make([]float64, a.mat.Cols))
			m.Sub(m, a)
			nt := Norm(m, 2)
			nr := floats.Norm(row, 2)
			if math.Abs(nt-nr) > 1e-14 {
				t.Errorf("Row %d norm mismatch, want: %g, got: %g", ri, nr, nt)
			}
		}

		for ci := range as[0] {
			a := NewDense(flatten(as))
			m := &Dense{}
			m.Clone(a)
			a.SetCol(ci, make([]float64, a.mat.Rows))
			col := make([]float64, a.mat.Rows)
			for j := range col {
				col[j] = float64(ci + 1 + j*a.mat.Cols)
			}
			m.Sub(m, a)
			nt := Norm(m, 2)
			nc := floats.Norm(col, 2)
			if math.Abs(nt-nc) > 1e-14 {
				t.Errorf("Column %d norm mismatch, want: %g, got: %g", ci, nc, nt)
			}
		}
	}
}

func (p *Printer) Record(loc *Location, _ EvaluationType, iter IterationType, stats *Stats) error {
	if iter != MajorIteration && iter != InitIteration && iter != PostIteration {
		return nil
	}

	// Print values always on PostIteration or when ValueInterval has elapsed.
	printValues := time.Since(p.lastValue) > p.ValueInterval || iter == PostIteration
	if !printValues {
		// Return early if not printing anything.
		return nil
	}

	// Print heading when HeadingInterval lines have been printed, but never on PostIteration.
	printHeading := p.lastHeading >= p.HeadingInterval && iter != PostIteration
	if printHeading {
		p.lastHeading = 1
	} else {
		p.lastHeading++
	}

	if printHeading {
		headings := "\n" + fmt.Sprintf(printerBaseTmpl, printerHeadings[0], printerHeadings[1], printerHeadings[2], printerHeadings[3])
		if loc.Gradient != nil {
			headings += fmt.Sprintf(printerGradTmpl, printerHeadings[4], printerHeadings[5])
		}
		if loc.Hessian != nil {
			headings += fmt.Sprintf(printerHessTmpl, printerHeadings[6])
		}
		_, err := fmt.Fprintln(p.Writer, headings)
		if err != nil {
			return err
		}
	}

	values := fmt.Sprintf(printerBaseTmpl, stats.MajorIterations, stats.Runtime, stats.FuncEvaluations, loc.F)
	if loc.Gradient != nil {
		values += fmt.Sprintf(printerGradTmpl, stats.GradEvaluations, floats.Norm(loc.Gradient, math.Inf(1)))
	}
	if loc.Hessian != nil {
		values += fmt.Sprintf(printerHessTmpl, stats.HessEvaluations)
	}
	_, err := fmt.Fprintln(p.Writer, values)
	if err != nil {
		return err
	}

	p.lastValue = time.Now()
	return nil
}

// Explicitly forms vectors and computes normalized dot product.
func cosCorrMultiNaive(f, g *rimg64.Multi) *rimg64.Image {
	h := rimg64.New(f.Width-g.Width+1, f.Height-g.Height+1)
	n := g.Width * g.Height * g.Channels
	a := make([]float64, n)
	b := make([]float64, n)
	for i := 0; i < h.Width; i++ {
		for j := 0; j < h.Height; j++ {
			a = a[:0]
			b = b[:0]
			for u := 0; u < g.Width; u++ {
				for v := 0; v < g.Height; v++ {
					for p := 0; p < g.Channels; p++ {
						a = append(a, f.At(i+u, j+v, p))
						b = append(b, g.At(u, v, p))
					}
				}
			}
			floats.Scale(1/floats.Norm(a, 2), a)
			floats.Scale(1/floats.Norm(b, 2), b)
			h.Set(i, j, floats.Dot(a, b))
		}
	}
	return h
}

func (l *LBFGS) InitDirection(loc *Location, dir []float64) (stepSize float64) {
	dim := len(loc.X)
	l.dim = dim

	if l.Store == 0 {
		l.Store = 15
	}

	l.oldest = l.Store - 1 // the first vector will be put in at 0

	l.x = resize(l.x, dim)
	l.grad = resize(l.grad, dim)
	copy(l.x, loc.X)
	copy(l.grad, loc.Gradient)

	l.y = resize(l.y, dim)
	l.s = resize(l.s, dim)
	l.a = resize(l.a, l.Store)
	l.rhoHist = resize(l.rhoHist, l.Store)

	if cap(l.yHist) < l.Store {
		n := make([][]float64, l.Store-cap(l.yHist))
		l.yHist = append(l.yHist, n...)
	}
	if cap(l.sHist) < l.Store {
		n := make([][]float64, l.Store-cap(l.sHist))
		l.sHist = append(l.sHist, n...)
	}
	l.yHist = l.yHist[:l.Store]
	l.sHist = l.sHist[:l.Store]
	for i := range l.sHist {
		l.sHist[i] = resize(l.sHist[i], dim)
		for j := range l.sHist[i] {
			l.sHist[i][j] = 0
		}
	}
	for i := range l.yHist {
		l.yHist[i] = resize(l.yHist[i], dim)
		for j := range l.yHist[i] {
			l.yHist[i][j] = 0
		}
	}

	copy(dir, loc.Gradient)
	floats.Scale(-1, dir)

	return 1 / floats.Norm(dir, 2)
}

// checkConvergence returns NotTerminated if the Location does not satisfy the
// convergence criteria given by settings. Otherwise a corresponding status is
// returned.
// Unlike checkLimits, checkConvergence is called by Local only at MajorIterations.
func checkConvergence(loc *Location, settings *Settings) Status {
	if loc.Gradient != nil {
		norm := floats.Norm(loc.Gradient, math.Inf(1))
		if norm < settings.GradientThreshold {
			return GradientThreshold
		}
	}

	if loc.F < settings.FunctionThreshold {
		return FunctionThreshold
	}

	if settings.FunctionConverge != nil {
		return settings.FunctionConverge.FunctionConverged(loc.F)
	}

	return NotTerminated
}

func (l *LBFGS) InitDirection(loc *Location, dir []float64) (stepSize float64) {
	dim := len(loc.X)
	l.dim = dim
	l.oldest = 0

	l.a = resize(l.a, l.Store)
	l.rho = resize(l.rho, l.Store)
	l.y = l.initHistory(l.y)
	l.s = l.initHistory(l.s)

	l.x = resize(l.x, dim)
	copy(l.x, loc.X)

	l.grad = resize(l.grad, dim)
	copy(l.grad, loc.Gradient)

	copy(dir, loc.Gradient)
	floats.Scale(-1, dir)
	return 1 / floats.Norm(dir, 2)
}

// Initialize initializes the Float to be ready to optimize by
// setting the history slice to have length zero, and setting
// the current value equal to the initial value
// This should be called by the optimizer at the beginning of
// the optimization
func (f *Floats) Init() error {
	f.Hist = f.Hist[:0]

	if f.Initial == nil {
		return errors.New("multivariate: initial slice is nil")
	}

	f.length = len(f.Initial)

	f.diff = make([]float64, len(f.Initial))
	f.current = make([]float64, len(f.Initial))
	f.previous = make([]float64, len(f.Initial))
	for i := range f.previous {
		f.previous[i] = math.Inf(1)
	}
	copy(f.current, f.Initial)
	f.norm = floats.Norm(f.current, 2)
	floats.SubTo(f.diff, f.current, f.previous)

	f.AddToHist(f.Initial)
	//f.prevNorm = math.Inf(1)
	return nil
}