docs: def nonlinear_optimization_uniform(
      X,
      cells,
      tol,
      max_num_steps,
      verbose=False,
      step_filename_format=None,
      callback=None,
  ):
      """Optimal Delaunay Triangulation smoothing.
  
      This method minimizes the energy
  
          E = int_Omega |u_l(x) - u(x)| rho(x) dx
  
      where u(x) = ||x||^2, u_l is its piecewise linear nodal interpolation and
      rho is the density. Since u(x) is convex, u_l >= u everywhere and
  
          u_l(x) = sum_i phi_i(x) u(x_i)
  
      where phi_i is the hat function at x_i. With rho(x)=1, this gives
  
          E = int_Omega sum_i phi_i(x) u(x_i) - u(x)
            = 1/(d+1) sum_i ||x_i||^2 |omega_i| - int_Omega ||x||^2
  
      where d is the spatial dimension and omega_i is the star of x_i (the set of
      all simplices containing x_i).
      """
      import scipy.optimize
  
      mesh = MeshTri(X, cells)
  
      if step_filename_format:
          mesh.save(
              step_filename_format.format(0),
              show_centroids=False,
              show_coedges=False,
              show_axes=False,
              nondelaunay_edge_color="k",
          )
  
      print("Before:")
      extra_cols = ["energy: {:.5e}".format(energy(mesh))]
      print_stats(mesh, extra_cols=extra_cols)
  
      def f(x):
          mesh.node_coords[mesh.is_interior_node] = x.reshape(-1, X.shape[1])
          mesh.update_values()
          return energy(mesh, uniform_density=True)
  
      # TODO put f and jac together
      def jac(x):
          mesh.node_coords[mesh.is_interior_node] = x.reshape(-1, X.shape[1])
          mesh.update_values()
  
          grad = numpy.zeros(mesh.node_coords.shape)
          cc = mesh.cell_circumcenters
          for mcn in mesh.cells["nodes"].T:
              fastfunc.add.at(
                  grad, mcn, ((mesh.node_coords[mcn] - cc).T * mesh.cell_volumes).T
              )
          gdim = 2
          grad *= 2 / (gdim + 1)
          return grad[mesh.is_interior_node].flatten()
  
      def flip_delaunay(x):
          flip_delaunay.step += 1
          # Flip the edges
          mesh.node_coords[mesh.is_interior_node] = x.reshape(-1, X.shape[1])
          mesh.update_values()
          mesh.flip_until_delaunay()
  
          if step_filename_format:
              mesh.save(
                  step_filename_format.format(flip_delaunay.step),
                  show_centroids=False,
                  show_coedges=False,
                  show_axes=False,
                  nondelaunay_edge_color="k",
              )
          if verbose:
              print("\nStep {}:".format(flip_delaunay.step))
              print_stats(mesh, extra_cols=["energy: {}".format(f(x))])
  
          if callback:
              callback(flip_delaunay.step, mesh)
  
          # mesh.show()
          # exit(1)
          return
  
      flip_delaunay.step = 0
  
      x0 = X[mesh.is_interior_node].flatten()
  
      if callback:
          callback(0, mesh)
  
      out = scipy.optimize.minimize(
          f,
          x0,
          jac=jac,
          # method="Nelder-Mead",
          # method="Powell",
          # method="CG",
          # method="Newton-CG",
          method="BFGS",
          # method="L-BFGS-B",
          # method="TNC",
          # method="COBYLA",
          # method="SLSQP",
          tol=tol,
          callback=flip_delaunay,
          options={"maxiter": max_num_steps},
      )
      # Don't assert out.success; max_num_steps may be reached, that's fine.
  
      # One last edge flip
      mesh.node_coords[mesh.is_interior_node] = out.x.reshape(-1, X.shape[1])
      mesh.update_values()
      mesh.flip_until_delaunay()
  
      print("\nFinal ({} steps):".format(out.nit))
      extra_cols = ["energy: {:.5e}".format(energy(mesh))]
      print_stats(mesh, extra_cols=extra_cols)
      print()
  
      return mesh.node_coords, mesh.cells["nodes"]