def bfs(matrix, to_visit, factors, codes, size):
    """Traverse the matrix across connected components.

    Implements breadth first search on an adjacency matrix. In our case, the
    adjacency matrix encodes the statistical relations between factors of
    variation and codes. This is used to traverse the adjacency matrix and
    discover whether a factor is captured in multiple codes and whether there is a
    path in the graph connecting two factors.

    Args:
      matrix: Adjacency matrix for the graph.
      to_visit: Queue with the nodes to visit. We index the factors and codes in
        the adjacency matrix and implement the queue with an array containing the
        nodes that need to  be visited.
      factors: Array of shape (num_factors, ) with flags marking whether factors
        of variation are visited.
      codes: Array of shape (num_codes, ) with flags marking whether codes are
        visited.
      size: Count how many node are in the same connected component.

    Returns:
      factors: Array of shape (num_factors, ) with flags marking whether factors
        of variation are visited.
      codes: Array of shape (num_codes, ) with flags marking whether codes are
        visited.
      size: How many nodes were visited.
    """
    (current_node, flag) = to_visit.pop()
    if flag == 0:
        factors[current_node] = 1
        for i in range(len(matrix[current_node, :])):
            if matrix[current_node, i] != 0:
                if codes[i] == 0:
                    to_visit.append((i, 1))
                    size += 1
                    factors, codes, size = bfs(matrix, to_visit, factors, codes, size)
    else:
        codes[current_node] = 1
        for i in range(len(matrix[:, current_node])):
            if matrix[i, current_node] != 0:
                if factors[i] == 0:
                    to_visit.append((i, 0))
                    size += 1
                    factors, codes, size = bfs(matrix, to_visit, factors, codes, size)
    return factors, codes, size