City Construction Solution times out for large input Problem


  The country of Hackerland has n cities connected by m
  uni-directional roads. The cities are numbered from 1 to n.
  Recently the government decided to build new cities in Hackerland.
  
  Your task is to simulate q events. An event can have two types as
  described below:
  
  
  1 x d : A new city n + 1 is constructed in the Hackerland and it is connected to city x. If d = 0 the direction of the new road
  is from x to n + 1. If d = 1 the direction of the new road is
  from n + 1 to x.
  2 x y : Print Yes if it's possible to move from city x to city y, print No otherwise.
  


Constraints


  0 &lt; n, m &lt;= 50000
  
  0 &lt; q &lt;= 10000
  
  x, y always correspond to an existing city in the Hackerland.
  
  Total number of cities in the Hackerland wont exceed 50000


You can find a full description of the problem here

My Effort

I tried to use a BFS here but all the test cases but the default are failing due to a timeout. For those test cases n &gt;= 4000. m &gt;= 20000.

Code

from collections import defaultdict, deque

n, m = map(int, input().split())

adjacentNodes = defaultdict(set)
for _ in range(m):
    u, v = map(int, input().split())
    adjacentNodes[u].add(v)

q = int(input())
for _ in range(q):
    c, x, d = map(int, input().split())

    if c == 1:
        n += 1
        y = n

        if d == 0:
            adjacentNodes[x].add(y)
        else:
            adjacentNodes[y].add(x)
    else:
        y = d
        visited = set()
        current = deque()
        current.append(x)
        answer = "No"
        while current:
            parent = current.popleft()
            if parent == y:
                answer = "Yes"
                break
            if parent in adjacentNodes:    
                for i in adjacentNodes[parent]:
                    if i not in visited:
                        visited.add(i)
                        current.append(i)
        print (answer) 


Question

Does anyone have a suggestion how to optimise the solution?
It was a tough one but I managed to solve it and pass all the test cases. Below are main points.

Build a full directed graph first

Save the queries because there can be up to 100,000 of them and it's more efficient to build a full graph first and then run all the queries on it.

Find all strongly connected components

Find all strongly connected components in the graph and then build a reachability map based on them. It has a lower space complexity since each strong component can be represented by a single node and so we don't need to have a reachability set for each node but only for its representative.

Use DFS for finding reachability sets

Modify a DFS so that after running it each connected component contains representatives (so-called low-link nodes) from other connected components that are reachable. But then, if a low link node is reachable so is any node from its connected component.

Use memory efficiently for storing strong components

Create an array of (n + 1) integer values, where n - the number of nodes in the graph. Index identifies a node in the graph and a value - all connected components (low link node) reachable from it. Value can have up to n + 1 binary digits (if all nodes are reachable from the current one).

Increase recursion limit 

As the number of nodes can reach 50000 so the recursion call stack must be set to this value.

Solution Code

from collections import defaultdict, deque
import sys

sys.setrecursionlimit(50001)

n, m = map(int, input().split())
adjacentNodes = defaultdict(set)
for _ in range(m):
    u, v = map(int, input().split())
    adjacentNodes[u].add(v)

q = int(input())
# read all queries first
queries = []
for i in range(q):
    c, x, d = map(int, input().split())
    if c == 1:
        n += 1
        y = n
        if d == 0:
            adjacentNodes[x].add(y)
        else:
            adjacentNodes[y].add(x)
    else:
        y = d
        queries.append((x, y))

# each node with number i corresponds to 2^i
nodeDpPosition = {i: (1 &lt;&lt; i) for i in range(n + 1)}
# Implementation using the Tarjan algorithm
# https://en.wikipedia.org/wiki/Tarjan%27s_strongly_connected_components_algorithm
def findAllConnectedComponents(edges):
    connectedComponents = [0 for _ in range(n + 1)]
    index = {'value' : 0}
    # map node to its connected component low link node
    connectedComponentLink = {}

    nodeStack = deque()
    indices = [-1] * (n + 1)
    lowLinks = [0] * (n + 1)
    onStack = [False] * (n + 1)

    def connectedComponent(node):
        indices[node] = index['value']
        lowLinks[node] = index['value']
        index['value'] += 1
        nodeStack.append(node)
        onStack[node] = True

        for adjacentNode in edges[node]:
            if indices[adjacentNode] &lt; 0:
                connectedComponent(adjacentNode)
                lowLinks[node] = min(lowLinks[node], lowLinks[adjacentNode])
            elif onStack[adjacentNode]:
                # we have reached the previously visited node or come back in the call stack
                lowLinks[node] = min(lowLinks[node], indices[adjacentNode])

        if lowLinks[node] == indices[node]:
            while nodeStack:
                current = nodeStack.pop()
                connectedComponentLink[current] = node
                onStack[current] = False
                connectedComponents[node] |= nodeDpPosition[current]

                if current == node:
                    break

    for i in range(1, n + 1):
        if indices[i] &lt; 0:
            connectedComponent(i)

    return (connectedComponents, connectedComponentLink)

connectedComponents, connectedComponentLink = findAllConnectedComponents(adjacentNodes)

def advancedDfs(node):
    lowLinkNode = connectedComponentLink[node]

    neighbours = adjacentNodes.get(node)
    if neighbours is not None:
        for neighbour in neighbours:
            lowLinkNeighbourNode = connectedComponentLink[neighbour]
            if (connectedComponents[lowLinkNode] &amp; nodeDpPosition[lowLinkNeighbourNode]) == 0:
                connectedComponents[lowLinkNode] |= nodeDpPosition[lowLinkNeighbourNode]
                advancedDfs(lowLinkNeighbourNode)

            connectedComponents[lowLinkNode] |= connectedComponents[lowLinkNeighbourNode]

for i in range(1, n + 1):
    advancedDfs(i)

for query in queries:
     x, y = query
     lowLinkXNode = connectedComponentLink[x]
     lowLinkYNode = connectedComponentLink[y]
     if connectedComponents[lowLinkXNode] &amp; nodeDpPosition[lowLinkYNode]:
         print ("Yes")
     else:
         print ("No")