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
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].