def compute_tukey_hh_params(l_skewness_and_kurtosis):
  """Computes the H paramesters of a Tukey HH distribution.

  Given the L-skewness and L-kurtosis of a Tukey HH distribution we compute
  the H parameters of the distribution.

  Args:
    l_skewness_and_kurtosis: A np.array with shape (2,) containing L-skewness
    and L-kurtosis.

  Returns:
    An np.array with the same type and shape of the argument containing the
    left and right H parameters of the distribution.
  """

  # Exit conditions for the search loop.
  stop_iter_step = 20  # Max number of iteration for the search loop.
  stop_error_step = 1e-6  # Minimum function variation.
  stop_value_step = 1e-6  # Minimum variable variation.

  dtype = l_skewness_and_kurtosis.dtype

  # Returns zero parameters (i.e. treat as gaussian) if L-kurtosis is smaller
  # than for a gaussian.

  result = np.zeros_like(l_skewness_and_kurtosis)
  if l_skewness_and_kurtosis[1] < 0.1226017:
    return result

  # If L-skewness is negative, swap the parameters.

  swap_params = False
  if l_skewness_and_kurtosis[0] < 0.0:
    l_skewness_and_kurtosis[0] = -l_skewness_and_kurtosis[0]
    swap_params = True

  l_skewness_and_kurtosis[1] = np.minimum(
      l_skewness_and_kurtosis[1], 1.0 - 1.0e-5)

  # If L-skewness is zero, left and right parameters are equal and there is a
  # a closed form to compute them from L-kurtosis. We start from this value
  # and then change them to match simultaneously L-skeweness and L-kurtosis.
  # For that, we parametrize the search space with the array
  # [h_rigth, h_right - h_left], i.e. the value of the right parameter and the
  # difference right minus left paramerters. In the search iteration, we
  # alternate between updates on the first and the second entry of the search
  # parameters.

  initial_h = 3.0 - 1.0 / np.cos(
      np.pi / 15.0 * (l_skewness_and_kurtosis[1] - 6.0))
  search_params = np.array([initial_h, 0.0], dtype=dtype)

  # Current lower and upper bounds for the search parameters.

  min_search_params = np.array([initial_h, 0.0], dtype=dtype)
  max_search_params = np.array([1.0 - 1.0e-7, initial_h], dtype=dtype)

  current_iter = 0
  previous_search_params = np.zeros_like(search_params)
  while current_iter < stop_iter_step:
    # Search for L-skewness at constant h. Increase delta_h.
    error_skewness = lambda x: _params_to_errors(  # pylint: disable=g-long-lambda
        search_params[0], x, l_skewness_and_kurtosis)[0]
    if error_skewness(max_search_params[1]) > 0.0:
      low_delta_h, high_delta_h = _binary_search(
          error_skewness, min_search_params[1], max_search_params[1])
      search_params[1] = high_delta_h
      max_search_params[1] = high_delta_h  # The new delta is an upperbound.
      upperbound_delta_found = True
    else:
      search_params[1] = max_search_params[1]
      min_search_params[1] = max_search_params[1]  # No solution: lowerbound.
      upperbound_delta_found = False

    # Search for L-kurtosis at constant possibly overestimated delta.
    error_kurtosis = lambda x: _params_to_errors(  # pylint: disable=g-long-lambda
        x, search_params[1], l_skewness_and_kurtosis)[1]
    low_h, high_h = _binary_search(
        error_kurtosis, min_search_params[0], max_search_params[0])
    if upperbound_delta_found:
      search_params[0] = high_h
      max_search_params[0] = high_h   # Delta overestimated: upperbound for h.
    else:
      search_params[0] = low_h
      min_search_params[0] = low_h   # Delta underestimated: lowerbound for h.
      max_search_params[1] = low_h  # Delta not found, search on full range.

    if upperbound_delta_found:  # If not found, we repeat the first 2 steps.
      # Otherwise, Search for delta at constant overestimated h.
      error_skewness = lambda x: _params_to_errors(  # pylint: disable=g-long-lambda
          search_params[0], x, l_skewness_and_kurtosis)[0]
      low_delta_h, high_delta_h = _binary_search(
          error_skewness, min_search_params[1], max_search_params[1])
      search_params[1] = low_delta_h
      min_search_params[1] = low_delta_h

      # Search for h at constant delta.
      error_kurtosis = lambda x: _params_to_errors(  # pylint: disable=g-long-lambda
          x, search_params[1], l_skewness_and_kurtosis)[1]
      low_h, high_h = _binary_search(
          error_kurtosis, min_search_params[0], max_search_params[0])
      search_params[0] = low_h
      min_search_params[0] = low_h

    current_error = _params_to_errors(
        search_params[0], search_params[1], l_skewness_and_kurtosis)
    delta_search_params = search_params - previous_search_params
    current_iter += 1
    previous_search_params = search_params.copy()
    if (np.all(np.abs(current_error) < stop_error_step) or
        np.all(np.abs(delta_search_params) < stop_value_step)):
      break

  result[0] = search_params[0] - search_params[1]
  result[1] = search_params[0]
  if swap_params:
    result = result[::-1]
  return result
Computes the H paramesters of a Tukey HH distribution.

  Given the L-skewness and L-kurtosis of a Tukey HH distribution we compute
  the H parameters of the distribution.

  Args:
    l_skewness_and_kurtosis: A np.array with shape (2,) containing L-skewness
    and L-kurtosis.

  Returns:
    An np.array with the same type and shape of the argument containing the
    left and right H parameters of the distribution.
  
tensorflow_transform/gaussianization.py