docs: def estimate_skeleton(indep_test_func, data_matrix, alpha, **kwargs):
      """Estimate a skeleton graph from the statistis information.
  
      Args:
          indep_test_func: the function name for a conditional
              independency test.
          data_matrix: data (as a numpy array).
          alpha: the significance level.
          kwargs:
              'max_reach': maximum value of l (see the code).  The
                  value depends on the underlying distribution.
              'method': if 'stable' given, use stable-PC algorithm
                  (see [Colombo2014]).
              'init_graph': initial structure of skeleton graph
                  (as a networkx.Graph). If not specified,
                  a complete graph is used.
              other parameters may be passed depending on the
                  indep_test_func()s.
      Returns:
          g: a skeleton graph (as a networkx.Graph).
          sep_set: a separation set (as an 2D-array of set()).
  
      [Colombo2014] Diego Colombo and Marloes H Maathuis. Order-independent
      constraint-based causal structure learning. In The Journal of Machine
      Learning Research, Vol. 15, pp. 3741-3782, 2014.
      """
  
      def method_stable(kwargs):
          return ('method' in kwargs) and kwargs['method'] == "stable"
  
      node_ids = range(data_matrix.shape[1])
      node_size = data_matrix.shape[1]
      sep_set = [[set() for i in range(node_size)] for j in range(node_size)]
      if 'init_graph' in kwargs:
          g = kwargs['init_graph']
          if not isinstance(g, nx.Graph):
              raise ValueError
          elif not g.number_of_nodes() == len(node_ids):
              raise ValueError('init_graph not matching data_matrix shape')
          for (i, j) in combinations(node_ids, 2):
              if not g.has_edge(i, j):
                  sep_set[i][j] = None
                  sep_set[j][i] = None
      else:
          g = _create_complete_graph(node_ids)
  
      l = 0
      while True:
          cont = False
          remove_edges = []
          for (i, j) in permutations(node_ids, 2):
              adj_i = list(g.neighbors(i))
              if j not in adj_i:
                  continue
              else:
                  adj_i.remove(j)
              if len(adj_i) >= l:
                  _logger.debug('testing %s and %s' % (i,j))
                  _logger.debug('neighbors of %s are %s' % (i, str(adj_i)))
                  if len(adj_i) < l:
                      continue
                  for k in combinations(adj_i, l):
                      _logger.debug('indep prob of %s and %s with subset %s'
                                    % (i, j, str(k)))
                      p_val = indep_test_func(data_matrix, i, j, set(k),
                                              **kwargs)
                      _logger.debug('p_val is %s' % str(p_val))
                      if p_val > alpha:
                          if g.has_edge(i, j):
                              _logger.debug('p: remove edge (%s, %s)' % (i, j))
                              if method_stable(kwargs):
                                  remove_edges.append((i, j))
                              else:
                                  g.remove_edge(i, j)
                          sep_set[i][j] |= set(k)
                          sep_set[j][i] |= set(k)
                          break
                  cont = True
          l += 1
          if method_stable(kwargs):
              g.remove_edges_from(remove_edges)
          if cont is False:
              break
          if ('max_reach' in kwargs) and (l > kwargs['max_reach']):
              break
  
      return (g, sep_set)