Buggy python submission (verdict: Time Limit Exceeded):
import scipy.sparse as sps
import scipy.sparse.csgraph as spsc
import numpy as np
from sys import stdin
 
N,M = map(int,stdin.readline().split())
 
S,T = map(int,stdin.readline().split())
S,T = S-1,T-1
 
E = sps.dok_matrix((N,N),dtype=int)
 
def it():
  lines = stdin.readlines()
  for line in lines:
    u,v,d = map(int,line.split())
    u,v = u-1,v-1
    yield ((u,v),d)
E.update(it())
 
 
dS,dT = spsc.dijkstra(E,directed=False, indices=[S,T]).astype(int)
D = dS[T]
D2 = D//2
 
E = tuple((u,v) for d,u,v in ((d,u,v) if dS[u] <= dS[v] else (d,v,u) for (u,v),d in E.items()) if dS[u]+d+dT[v] == D)
 
MOD = 10**9+7
 
 
def helper(start, edges):
  cnts = [0]*N
  cnts[start] = 1
  for d,u,v in edges:
    cnts[v] = (cnts[v] + cnts[u]) % MOD
  return cnts
 
 
Sn = helper(S,sorted((dS[u],u,v) for u,v in E))
Tn = helper(T,sorted((dT[v],v,u) for u,v in E))
 
from functools import reduce
 
def it():
  for v in range(N):
    if dS[v] == D2:
      yield (Sn[v]*Tn[v])**2
 
c0 = reduce(lambda x,y: (x+y)%MOD, it(), 0) if D%2==0 else 0
 
diff = dS-dT
 
def it():
  for u,v in E:
    if diff[u] < 0 < diff[v]:
      yield (Sn[u]*Tn[v])**2
 
c1 = reduce(lambda x,y: (x+y)%MOD, it(), 0)
 
total = Sn[T]*Tn[S]
cnt = c0+c1
 
print((total-cnt)%MOD)

Fixed submission:
from collections import defaultdict
from functools import reduce
import heapq as hq
from sys import stdin
 
MOD = 10**9+7
N,M = map(int,stdin.readline().split())
S,T = map(int,stdin.readline().split())
S,T = S-1,T-1
 
E = defaultdict(list)
lines = stdin.readlines()
for line in lines:
  u,v,d = map(int,line.split())
  E[u-1].append((v-1,d))
  E[v-1].append((u-1,d))
 
dS = [None]*N
dT = [None]*N
 
Q = [(0,S,False),(0,T,True)]
 
D = float('inf')
 
 
while Q:
  d,v,f = hq.heappop(Q)
  if d*2 > D:
    break

  dM = dT if f else dS
  dN = dS if f else dT
  if dM[v] is None:
    dM[v] = d
    for u,dd in E[v]:
      if dM[u] is None and (d+dd)*2 <= D:
        hq.heappush(Q,(d+dd,u,f))
 
  if dN[v] is not None:
    D = min(dN[v]+d,D)
 
del Q

def helper(start,dist):
  cnts = [0]*N
  cnts[start] = 1

  for d,v in sorted((d,v) for v,d in enumerate(dist) if d is not None):
    c = cnts[v]
    for u,dd in E[v]:
      if dd+d == dist[u]:
        cnts[u] = (cnts[u]+c) % MOD

  return cnts

Sn = helper(S,dS)
Tn = helper(T,dT)


def it():
  for v in range(N):
    if dS[v] is not None:
      ds = dS[v]
      if dT[v] == ds and ds*2 == D:
        yield Sn[v]*Tn[v] % MOD
      else:
        for u,d in E[v]:
          if dT[u] is not None and abs(dT[u]-ds) < d and ds+d+dT[u] == D:
            yield Sn[v]*Tn[u] % MOD


X = list(it())

 
total = reduce(lambda a,b:(a+b)%MOD, X,0)
total = total*total % MOD
 
sub = reduce(lambda a,b:(a+b)%MOD, map(lambda x: x*x % MOD, X),0)
print((total-sub)%MOD)

Buggy python submission (verdict: Time Limit Exceeded):
from collections import defaultdict
from functools import reduce
import heapq as hq

MOD = 10**9+7
N,M = map(int,input().split())
S,T = map(int,input().split())
S,T = S-1,T-1

E = defaultdict(list)

for _ in range(M):
  u,v,d = map(int,input().split())
  E[u-1].append((v-1,d))
  E[v-1].append((u-1,d))

dist = [None]*N
flag = [None]*N

Q = [(0,S,False),(0,T,True)]

D = float('inf')

while Q:
  d,v,f = hq.heappop(Q)
  if d*2 > D:
    break
  if flag[v] is None:
    flag[v] = f
    dist[v] = d
    for u,dd in E[v]:
      if flag[u] != f:
        hq.heappush(Q,(d+dd,u,f))

  elif flag[v] != f:
    D = min(dist[v]+d,D)

del Q


cnts = [0]*N
cnts[S] = 1
cnts[T] = 1

for d,v in sorted((d,v) for v,d in enumerate(dist) if d is not None):
  f = flag[v]
  c = cnts[v]
  for u,dd in E[v]:
    if f == flag[u] and d+dd == dist[u]:
      cnts[u] = (cnts[u]+c)%MOD

def it():
  for v in range(N):
    for u,d in E[v]:
      if flag[v] == True and flag[u] == False and dist[v]+d+dist[u] == D:
        yield cnts[u]*cnts[v] % MOD
X = list(it())

total = reduce(lambda a,b:(a+b)%MOD, X)
total = total*total % MOD

sub = reduce(lambda a,b:(a+b)%MOD, map(lambda x: x*x % MOD, X))

print((total-sub)%MOD)

Fixed submission:
from collections import defaultdict
from functools import reduce
import heapq as hq
from sys import stdin
 
MOD = 10**9+7
N,M = map(int,stdin.readline().split())
S,T = map(int,stdin.readline().split())
S,T = S-1,T-1
 
E = defaultdict(list)
lines = stdin.readlines()
for line in lines:
  u,v,d = map(int,line.split())
  E[u-1].append((v-1,d))
  E[v-1].append((u-1,d))
 
dS = [None]*N
dT = [None]*N
 
Q = [(0,S,False),(0,T,True)]
 
D = float('inf')
 
 
while Q:
  d,v,f = hq.heappop(Q)
  if d*2 > D:
    break

  dM = dT if f else dS
  dN = dS if f else dT
  if dM[v] is None:
    dM[v] = d
    for u,dd in E[v]:
      if dM[u] is None and (d+dd)*2 <= D:
        hq.heappush(Q,(d+dd,u,f))
 
  if dN[v] is not None:
    D = min(dN[v]+d,D)
 
del Q

def helper(start,dist):
  cnts = [0]*N
  cnts[start] = 1

  for d,v in sorted((d,v) for v,d in enumerate(dist) if d is not None):
    c = cnts[v]
    for u,dd in E[v]:
      if dd+d == dist[u]:
        cnts[u] = (cnts[u]+c) % MOD

  return cnts

Sn = helper(S,dS)
Tn = helper(T,dT)


def it():
  for v in range(N):
    if dS[v] is not None:
      ds = dS[v]
      if dT[v] == ds and ds*2 == D:
        yield Sn[v]*Tn[v] % MOD
      else:
        for u,d in E[v]:
          if dT[u] is not None and abs(dT[u]-ds) < d and ds+d+dT[u] == D:
            yield Sn[v]*Tn[u] % MOD


X = list(it())

 
total = reduce(lambda a,b:(a+b)%MOD, X,0)
total = total*total % MOD
 
sub = reduce(lambda a,b:(a+b)%MOD, map(lambda x: x*x % MOD, X),0)
print((total-sub)%MOD)

Buggy python submission (verdict: Time Limit Exceeded):
import scipy.sparse as sps
import scipy.sparse.csgraph as spsc
import numpy as np
from time import process_time as timer

N,M = map(int,input().split())

S,T = map(int,input().split())
S,T = S-1,T-1

E = sps.dok_matrix((N,N),dtype=int)

def it():
  for _ in range(M):
    u,v,d = map(int,input().split())
    u,v = u-1,v-1
    yield ((u,v),d)
E.update(it())


dS,dT = spsc.dijkstra(E,directed=False, indices=[S,T]).astype(int)
D = dS[T]
D2 = D//2

E = tuple((u,v) for d,u,v in ((d,u,v) if dS[u] <= dS[v] else (d,v,u) for (u,v),d in E.items()) if dS[u]+d+dT[v] == D)

MOD = 10**9+7


def helper(start, edges):
  cnts = [0]*N
  cnts[start] = 1
  for d,u,v in edges:
    cnts[v] = (cnts[v] + cnts[u]) % MOD
  return cnts


Sn = helper(S,sorted((dS[u],u,v) for u,v in E))
Tn = helper(T,sorted((dT[v],v,u) for u,v in E))

from functools import reduce

def it():
  for v in range(N):
    if dS[v] == D2:
      yield (Sn[v]*Tn[v])**2

c0 = reduce(lambda x,y: (x+y)%MOD, it(), 0) if D%2==0 else 0

diff = dS-dT

def it():
  for u,v in E:
    if diff[u] < 0 < diff[v]:
      yield (Sn[u]*Tn[v])**2

c1 = reduce(lambda x,y: (x+y)%MOD, it(), 0)

total = Sn[T]*Tn[S]
cnt = c0+c1

print((total-cnt)%MOD)


Fixed submission:
from collections import defaultdict
from functools import reduce
import heapq as hq
from sys import stdin
 
MOD = 10**9+7
N,M = map(int,stdin.readline().split())
S,T = map(int,stdin.readline().split())
S,T = S-1,T-1
 
E = defaultdict(list)
lines = stdin.readlines()
for line in lines:
  u,v,d = map(int,line.split())
  E[u-1].append((v-1,d))
  E[v-1].append((u-1,d))
 
dS = [None]*N
dT = [None]*N
 
Q = [(0,S,False),(0,T,True)]
 
D = float('inf')
 
 
while Q:
  d,v,f = hq.heappop(Q)
  if d*2 > D:
    break

  dM = dT if f else dS
  dN = dS if f else dT
  if dM[v] is None:
    dM[v] = d
    for u,dd in E[v]:
      if dM[u] is None and (d+dd)*2 <= D:
        hq.heappush(Q,(d+dd,u,f))
 
  if dN[v] is not None:
    D = min(dN[v]+d,D)
 
del Q

def helper(start,dist):
  cnts = [0]*N
  cnts[start] = 1

  for d,v in sorted((d,v) for v,d in enumerate(dist) if d is not None):
    c = cnts[v]
    for u,dd in E[v]:
      if dd+d == dist[u]:
        cnts[u] = (cnts[u]+c) % MOD

  return cnts

Sn = helper(S,dS)
Tn = helper(T,dT)


def it():
  for v in range(N):
    if dS[v] is not None:
      ds = dS[v]
      if dT[v] == ds and ds*2 == D:
        yield Sn[v]*Tn[v] % MOD
      else:
        for u,d in E[v]:
          if dT[u] is not None and abs(dT[u]-ds) < d and ds+d+dT[u] == D:
            yield Sn[v]*Tn[u] % MOD


X = list(it())

 
total = reduce(lambda a,b:(a+b)%MOD, X,0)
total = total*total % MOD
 
sub = reduce(lambda a,b:(a+b)%MOD, map(lambda x: x*x % MOD, X),0)
print((total-sub)%MOD)

Buggy python submission (verdict: Time Limit Exceeded):
import scipy.sparse as sps
import scipy.sparse.csgraph as spsc
import numpy as np

N,M = map(int,input().split())

S,T = map(int,input().split())
S,T = S-1,T-1

E = sps.dok_matrix((N,N),dtype=int)

adj = list(list() for _ in range(N))

def it():
  for _ in range(M):
    u,v,d = map(int,input().split())
    u,v = u-1,v-1
    adj[u].append((v,d))
    adj[v].append((u,d))
    yield ((u,v),d)
    yield ((v,u),d)
E.update(it())

dS,dT = spsc.dijkstra(E,directed=False, indices=[S,T])
D = dS[T]

V = dS+dT == D
MOD = 10**9+7

def helper(start, dist):
  import heapq as hq
  cnts = [0]*N
  cnts[start] = 1

  edges = []
  for u in range(N):
    if V[u]:
      for v,d in adj[u]:
        if V[v] and dist[v]-dist[u] == d:
          edges.append((dist[u],u,v))

  edges.sort()

  for _,u,v in edges:
    cnts[v] += cnts[u]
    cnts[v] %= MOD

  return cnts


Sn = helper(S,dS)
Tn = helper(T,dT)


c0 = 0
c1 = 0

for v in range(N):
  if V[v]:
    dv = dS[v]

    if dv*2 == D:
      c0 += Sn[v]*Tn[v]
      c0 %= MOD

    elif dv*2 < D:
      for u,d in adj[v]:
        du = dS[u]
        if V[u] and du*2 > D and du-dv == d:
          c1 += Sn[v]*Tn[u]
          c1 %= MOD

total = Sn[T]*Tn[S]
cnt = c0+c1

print((total-cnt)%MOD)


Fixed submission:
from collections import defaultdict
from functools import reduce
import heapq as hq
from sys import stdin
 
MOD = 10**9+7
N,M = map(int,stdin.readline().split())
S,T = map(int,stdin.readline().split())
S,T = S-1,T-1
 
E = defaultdict(list)
lines = stdin.readlines()
for line in lines:
  u,v,d = map(int,line.split())
  E[u-1].append((v-1,d))
  E[v-1].append((u-1,d))
 
dS = [None]*N
dT = [None]*N
 
Q = [(0,S,False),(0,T,True)]
 
D = float('inf')
 
 
while Q:
  d,v,f = hq.heappop(Q)
  if d*2 > D:
    break

  dM = dT if f else dS
  dN = dS if f else dT
  if dM[v] is None:
    dM[v] = d
    for u,dd in E[v]:
      if dM[u] is None and (d+dd)*2 <= D:
        hq.heappush(Q,(d+dd,u,f))
 
  if dN[v] is not None:
    D = min(dN[v]+d,D)
 
del Q

def helper(start,dist):
  cnts = [0]*N
  cnts[start] = 1

  for d,v in sorted((d,v) for v,d in enumerate(dist) if d is not None):
    c = cnts[v]
    for u,dd in E[v]:
      if dd+d == dist[u]:
        cnts[u] = (cnts[u]+c) % MOD

  return cnts

Sn = helper(S,dS)
Tn = helper(T,dT)


def it():
  for v in range(N):
    if dS[v] is not None:
      ds = dS[v]
      if dT[v] == ds and ds*2 == D:
        yield Sn[v]*Tn[v] % MOD
      else:
        for u,d in E[v]:
          if dT[u] is not None and abs(dT[u]-ds) < d and ds+d+dT[u] == D:
            yield Sn[v]*Tn[u] % MOD


X = list(it())

 
total = reduce(lambda a,b:(a+b)%MOD, X,0)
total = total*total % MOD
 
sub = reduce(lambda a,b:(a+b)%MOD, map(lambda x: x*x % MOD, X),0)
print((total-sub)%MOD)