docs: def adaptive(u_kn, N_k, f_k, tol = 1.0e-12, options = None):
  
      """
      Determine dimensionless free energies by a combination of Newton-Raphson iteration and self-consistent iteration.
      Picks whichever method gives the lowest gradient.
      Is slower than NR since it calculates the log norms twice each iteration.
  
      OPTIONAL ARGUMENTS
      tol (float between 0 and 1) - relative tolerance for convergence (default 1.0e-12)
  
      options: dictionary of options
          gamma (float between 0 and 1) - incrementor for NR iterations (default 1.0).  Usually not changed now, since adaptively switch.
          maximum_iterations (int) - maximum number of Newton-Raphson iterations (default 250: either NR converges or doesn't, pretty quickly)
          verbose (boolean) - verbosity level for debug output
  
      NOTES
  
  
      This method determines the dimensionless free energies by
      minimizing a convex function whose solution is the desired
      estimator.  The original idea came from the construction of a
      likelihood function that independently reproduced the work of
      Geyer (see [1] and Section 6 of [2]).  This can alternatively be
      formulated as a root-finding algorithm for the Z-estimator.  More
      details of this procedure will follow in a subsequent paper.  Only
      those states with nonzero counts are include in the estimation
      procedure.
  
      REFERENCES
      See Appendix C.2 of [1].
  
      """
      # put the defaults here in case we get passed an 'options' dictionary that is only partial
      options.setdefault('verbose',False)
      options.setdefault('maximum_iterations',250)
      options.setdefault('print_warning',False)
      options.setdefault('gamma',1.0)
  
      gamma = options['gamma']
      doneIterating = False
      if options['verbose'] == True:
          print("Determining dimensionless free energies by Newton-Raphson / self-consistent iteration.")
  
      if tol < 1.5e-15:
          print("Tolerance may be too close to machine precision to converge.")
      # keep track of Newton-Raphson and self-consistent iterations
      nr_iter = 0
      sci_iter = 0
  
      f_sci = np.zeros(len(f_k), dtype=np.float64)
      f_nr = np.zeros(len(f_k), dtype=np.float64)
  
      # Perform Newton-Raphson iterations (with sci computed on the way)
      for iteration in range(0, options['maximum_iterations']):
          g = mbar_gradient(u_kn, N_k, f_k)  # Objective function gradient
          H = mbar_hessian(u_kn, N_k, f_k)  # Objective function hessian
          Hinvg = np.linalg.lstsq(H, g, rcond=-1)[0]
          Hinvg -= Hinvg[0]
          f_nr = f_k - gamma * Hinvg
  
          # self-consistent iteration gradient norm and saved log sums.
          f_sci = self_consistent_update(u_kn, N_k, f_k)
          f_sci = f_sci -  f_sci[0]   # zero out the minimum
          g_sci = mbar_gradient(u_kn, N_k, f_sci)
          gnorm_sci = np.dot(g_sci, g_sci)
  
          # newton raphson gradient norm and saved log sums.
          g_nr = mbar_gradient(u_kn, N_k, f_nr)
          gnorm_nr = np.dot(g_nr, g_nr)
  
          # we could save the gradient, for the next round, but it's not too expensive to
          # compute since we are doing the Hessian anyway.
  
          if options['verbose']:
              print("self consistent iteration gradient norm is %10.5g, Newton-Raphson gradient norm is %10.5g" % (gnorm_sci, gnorm_nr))
          # decide which directon to go depending on size of gradient norm
          f_old = f_k
          if (gnorm_sci < gnorm_nr or sci_iter < 2):
              f_k = f_sci
              sci_iter += 1
              if options['verbose']:
                  if sci_iter < 2:
                      print("Choosing self-consistent iteration on iteration %d" % iteration)
                  else:
                      print("Choosing self-consistent iteration for lower gradient on iteration %d" % iteration)
          else:
              f_k = f_nr
              nr_iter += 1
              if options['verbose']:
                  print("Newton-Raphson used on iteration %d" % iteration)
  
          div = np.abs(f_k[1:]) # what we will divide by to get relative difference
          zeroed = np.abs(f_k[1:])< np.min([10**-8,tol]) # check which values are near enough to zero, hard coded max for now.
          div[zeroed] = 1.0  # for these values, use absolute values.
          max_delta = np.max(np.abs(f_k[1:]-f_old[1:])/div)
          if np.isnan(max_delta) or (max_delta < tol):
              doneIterating = True
              break
  
      if doneIterating:
          if options['verbose']:
              print('Converged to tolerance of {:e} in {:d} iterations.'.format(max_delta, iteration + 1))
              print('Of {:d} iterations, {:d} were Newton-Raphson iterations and {:d} were self-consistent iterations'.format(iteration + 1, nr_iter, sci_iter))
              if np.all(f_k == 0.0):
                  # all f_k appear to be zero
                  print('WARNING: All f_k appear to be zero.')
      else:
          print('WARNING: Did not converge to within specified tolerance.')
          if options['maximum_iterations'] <= 0:
              print("No iterations ran be cause maximum_iterations was <= 0 ({})!".format(options['maximum_iterations']))
          else:
              print('max_delta = {:e}, tol = {:e}, maximum_iterations = {:d}, iterations completed = {:d}'.format(max_delta,tol, options['maximum_iterations'], iteration))
      return f_k