docs: def _adaptive(self, gamma=1.0, relative_tolerance=1.0e-8, maximum_iterations=1000, verbose=True, print_warning=True):
          """
          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 (approximated, not calculated) since it calculates the log norms twice each iteration.
  
          OPTIONAL ARGUMENTS
            gamma (float between 0 and 1) - incrementor for NR iterations.
            relative_tolerance (float between 0 and 1) - relative tolerance for convergence (default 1.0e-6)
            maximum_iterations (int) - maximum number of Newton-Raphson iterations (default 1000)
            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].
  
          """
  
          if verbose:
              print("Determining dimensionless free energies by Newton-Raphson iteration.")
  
          # keep track of Newton-Raphson and self-consistent iterations
          nr_iter = 0
          sci_iter = 0
  
          N_k = self.N_k[self.states_with_samples]
          K = len(N_k)
  
          f_k_sci = np.zeros([K], dtype=np.float64)
          f_k_new = np.zeros([K], dtype=np.float64)
  
          # Perform Newton-Raphson iterations (with sci computed on the way)
          for iteration in range(0, maximum_iterations):
  
              # Store for new estimate of dimensionless relative free energies.
              f_k = self.f_k[self.states_with_samples].copy()
  
              # compute weights for gradients: the denominators and free energies are from the previous
              # iteration in most cases.
              (W_nk, f_k_sci) = self._computeWeights(
                  recalc_denom=(iteration == 0), return_f_k = True)
  
              # Compute gradient and Hessian of last (K-1) states.
              #
              # gradient (defined by Eq. C6 of [1])
              # g_i(theta) = N_i - \sum_n N_i W_ni
              #
              # Hessian (defined by Eq. C9 of [1])
              # H_ii(theta) = - \sum_n N_i W_ni (1 - N_i W_ni)
              # H_ij(theta) = \sum_n N_i W_ni N_j W_nj
              #
  
              """
        g = np.matrix(np.zeros([K-1,1], dtype=np.float64)) # gradient
        H = np.matrix(np.zeros([K-1,K-1], dtype=np.float64)) # Hessian
        for i in range(1,K):
          g[i-1] = N_k[i] - N_k[i] * W_nk[:,i].sum()
          H[i-1,i-1] = - (N_k[i] * W_nk[:,i] * (1.0 - N_k[i] * W_nk[:,i])).sum() 
          for j in range(1,i):
            H[i-1,j-1] = (N_k[i] * W_nk[:,i] * N_k[j] * W_nk[:,j]).sum()
            H[j-1,i-1] = H[i-1,j-1]
  
        # Update the free energy estimate (Eq. C11 of [1]).
        Hinvg = linalg.lstsq(H,g)[0]      #
        # Hinvg = linalg.solve(H,g)       # This might be faster if we can guarantee full rank.
        for k in range(0,K-1):
          f_k_new[k+1] = f_k[k+1] - gamma*Hinvg[k]
  
        """
              g = N_k - N_k * W_nk.sum(axis=0)
              NW = N_k * W_nk
              H = np.dot(NW.T, NW)
              H += (g.T - N_k) * np.eye(K)
              # Update the free energy estimate (Eq. C11 of [1]).
              # will always have lower rank the way it is set up
              Hinvg = linalg.lstsq(H, g)[0]
              Hinvg -= Hinvg[0]
              f_k_new = f_k - gamma * Hinvg
  
              # self-consistent iteration gradient norm and saved log sums.
              g_sci = self._gradientF(f_k_sci)
              gnorm_sci = np.dot(g_sci, g_sci)
              # save this so we can switch it back in if g_sci is lower.
              log_weight_denom = self.log_weight_denom.copy()
  
              # newton raphson gradient norm and saved log sums.
              g_nr = self._gradientF(f_k_new)
              gnorm_nr = np.dot(g_nr, g_nr)
  
              # we could save the gradient, too, but it's not too expensive to
              # compute since we are doing the Hessian anyway.
  
              if 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
              if (gnorm_sci < gnorm_nr or sci_iter < 2):
                  sci_iter += 1
                  self.log_weight_denom = log_weight_denom.copy()
                  if 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)
  
                  f_k_new = f_k_sci.copy()
              else:
                  nr_iter += 1
                  if verbose:
                      print("Newton-Raphson used on iteration %d" % iteration)
  
              # get rid of big matrices that are not used.
              del(log_weight_denom, NW, W_nk)
  
              # have to set the free energies back in self, since the gradient
              # routine changes them.
              self.f_k[self.states_with_samples] = f_k
              if (self._amIdoneIterating(f_k_new, relative_tolerance, iteration, maximum_iterations, print_warning, verbose)):
                  if verbose:
                      print('Of %d iterations, %d were Newton-Raphson iterations and %d were self-consistent iterations' % (iteration + 1, nr_iter, sci_iter))
                  break
  
          return