Problem:
 Score : 700 points    Problem Statement We have a graph with N vertices and M edges, and there are two people on the graph: Takahashi and Aoki.  The i -th edge connects Vertex U_i and Vertex V_i .
The time it takes to traverse this edge is D_i minutes, regardless of direction and who traverses the edge (Takahashi or Aoki).  Takahashi departs Vertex S and Aoki departs Vertex T at the same time. Takahashi travels to Vertex T and Aoki travels to Vertex S , both in the shortest time possible.
Find the number of the pairs of ways for Takahashi and Aoki to choose their shortest paths such that they never meet (at a vertex or on an edge) during the travel, modulo 10^9 + 7 .      Constraints  1 \leq N \leq 100  000  1 \leq M \leq 200  000  1 \leq S, T \leq N  S \neq T  1 \leq U_i, V_i \leq N ( 1 \leq i \leq M )  1 \leq D_i \leq 10^9 ( 1 \leq i \leq M )  If i \neq j , then (U_i, V_i) \neq (U_j, V_j) and (U_i, V_i) \neq (V_j, U_j) .  U_i \neq V_i ( 1 \leq i \leq M )  D_i are integers.  The given graph is connected.         Input Input is given from Standard Input in the following format:  N  M  S  T  U_1  V_1  D_1  U_2  V_2  D_2  :  U_M  V_M  D_M       Output Print the answer.        Sample Input 1 4 4
1 3
1 2 1
2 3 1
3 4 1
4 1 1      Sample Output 1 2  There are two ways to choose shortest paths that satisfies the condition:   Takahashi chooses the path 1 \rightarrow 2 \rightarrow 3 , and Aoki chooses the path 3 \rightarrow 4 \rightarrow 1 .  Takahashi chooses the path 1 \rightarrow 4 \rightarrow 3 , and Aoki chooses the path 3 \rightarrow 2 \rightarrow 1 .        Sample Input 2 3 3
1 3
1 2 1
2 3 1
3 1 2      Sample Output 2 2       Sample Input 3 3 3
1 3
1 2 1
2 3 1
3 1 2      Sample Output 3 2       Sample Input 4 8 13
4 2
7 3 9
6 2 3
1 6 4
7 6 9
3 8 9
1 2 2
2 8 12
8 6 9
2 5 5
4 2 18
5 3 7
5 1 515371567
4 8 6      Sample Output 4 6

Solution (Python):
import sys
readline = sys.stdin.readline

from heapq import heappop as hpp, heappush as hp
def dijkstra(N, s, Edge):
    inf = geta
    dist = [inf] * N
    dist[s] = 0
    Q = [(0, s)]
    dp = [0]*N
    dp[s] = 1
    while Q:
        dn, vn = hpp(Q)
        if dn > dist[vn]:
            continue
        for df, vf in Edge[vn]:
            if dist[vn] + df < dist[vf]:
                dist[vf] = dist[vn] + df
                dp[vf] = dp[vn]
                hp(Q, (dn + df,vf))
            elif dist[vn] + df == dist[vf]:
                dp[vf] = (dp[vf] + dp[vn]) % MOD
    return dist, dp

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

geta = 10**15
Edge = [[] for _ in range(N)]
for _ in range(M):
    u, v, d = map(int, readline().split())
    u -= 1
    v -= 1
    Edge[u].append((d, v))
    Edge[v].append((d, u))
    

dists, dps = dijkstra(N, S, Edge)
distt, dpt = dijkstra(N, T, Edge)
st = dists[T]

onpath = [i for i in range(N) if dists[i] + distt[i] == st]

ans = dps[T]**2%MOD
for i in onpath:
    if 2*dists[i] == 2*distt[i] == st:
        ans = (ans - pow(dps[i]*dpt[i], 2, MOD))%MOD
    
    for cost, vf in Edge[i]:
        if dists[i] + cost + distt[vf] == st:
            if 2*dists[i] < st < 2*dists[vf]:
                ans = (ans - pow(dps[i]*dpt[vf], 2, MOD))%MOD
print(ans)

Solution (Python):
def main():
    import sys
    sys.setrecursionlimit(100000)
    input = sys.stdin.readline
    from heapq import heappop, heappush
    mod = 10**9+7
    INF = 1<<60
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    E = [[] for _ in range(N+1)]
    edges = []
    
    
    
    
    
    for u, v, d in zip(*[iter(map(int, sys.stdin.read().split()))]*3):
        E[u].append((v, d))
        E[v].append((u, d))
        edges.append((u, v, d))

    def dijksrtra(start):
        Dist = [INF] * (N+1)
        Dist[start] = 0
        q = [0<<17 | start]
        mask = (1<<17)-1
        while q:
            dist_v = heappop(q)
            v = dist_v & mask
            dist = dist_v >> 17
            if Dist[v] != dist:
                continue
            for u, d in E[v]:
                new_dist = dist + d
                if Dist[u] > new_dist:
                    Dist[u] = new_dist
                    heappush(q, new_dist<<17 | u)
        return Dist

    Dist_S = dijksrtra(S)
    Dist_T = dijksrtra(T)
    dist_st = Dist_S[T]

    DAG_edges = []  
    DAG = [[] for _ in range(N+1)]
    DAG_rev = [[] for _ in range(N+1)]
    for u, v, d in edges:
        if Dist_S[u] + Dist_T[v] + d == dist_st:
            DAG_edges.append((u, v))
            DAG[u].append(v)
            DAG_rev[v].append(u)
        elif Dist_T[u] + Dist_S[v] + d == dist_st:
            DAG_edges.append((v, u))
            DAG[v].append(u)
            DAG_rev[u].append(v)


    
    V = []
    rem = [0] * (N+1)
    for _, v in DAG_edges:
        rem[v] += 1
    q = [S]
    while q:
        v = q.pop()
        V.append(v)
        for u in DAG[v]:
            rem[u] -= 1
            if rem[u] == 0:
                q.append(u)


    
    
    n_paths_S = [0] * (N+1)
    n_paths_T = [0] * (N+1)
    n_paths_S[S] = 1
    n_paths_T[T] = 1
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    

    
    

    for v in V:
        n = n_paths_S[v]
        for u in DAG[v]:
            n_paths_S[u] += n
    for v in V[::-1]:
        n = n_paths_T[v]
        for u in DAG_rev[v]:
            n_paths_T[u] += n
    ans = n_paths_S[T]
    ans = ans * ans % mod

    for v, u in DAG_edges:  
        if Dist_S[v]*2 < dist_st and dist_st < Dist_S[u]*2:
            ans = (ans - (n_paths_S[v]*n_paths_T[u])**2) % mod

    
    for v, (dist, ns, nt) in enumerate(zip(Dist_S, n_paths_S, n_paths_T)):
        if dist*2 == dist_st:
            ans = (ans - (ns*nt)**2) % mod

    print(ans)

main()

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from functools import reduce


sys.setrecursionlimit(2147483647)
INF = 10 ** 15
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007


n, m = LI()
s, t = LI()
G = [[] for _ in range(n)]
for i in range(m):
    a, b, c = LI()
    G[a - 1] += [(b - 1, c)]
    G[b - 1] += [(a - 1, c)]


def dijkstra(graph, start=0):
    route_cnt = [start] * n
    route_cnt[start] = 1
    que = [(0, start)]
    dist = [INF] * n
    dist[start] = 0
    while que:
        min_dist, u = heappop(que)
        if min_dist > dist[u]:
            continue
        for v, c in graph[u]:
            if dist[u] + c < dist[v]:
                dist[v] = dist[u] + c
                route_cnt[v] = route_cnt[u]
                heappush(que, (dist[u] + c  , v))
            elif dist[u] + c == dist[v]:
                route_cnt[v] = (route_cnt[v] + route_cnt[u]) % mod
    return route_cnt, dist


route_cnt, dist = dijkstra(G, s - 1)
r_route_cnt, r_dist = dijkstra(G, t - 1)
total_dist = dist[t - 1]
ret = 0
for i in range(n):
    if dist[i] == r_dist[i]:
        ret = ret + (route_cnt[i] * r_route_cnt[i] % mod) ** 2 % mod


for u in range(n):
    for v, c in G[u]:
        if dist[u] < total_dist / 2 < dist[v] and dist[u] + c + r_dist[v] == total_dist:
            ret = (ret + route_cnt[u] ** 2 % mod * r_route_cnt[v] ** 2 % mod) % mod


print((route_cnt[t - 1] ** 2 - ret) % mod)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**6)
from heapq import *

"""
とりあえず何通りあるかを数える：各頂点を通る時点での場合の数
最後に衝突する方法を除外する：頂点での衝突、辺での衝突
辺での衝突：距離を見れば分かる
"""

MOD = 10 ** 9 + 7

N,M = map(int,input().split())
S,T = map(int,input().split())
graph = [[] for _ in range(N+1)]
for _ in range(M):
    u,v,d = map(int,input().split())
    graph[u].append((v,d))
    graph[v].append((u,d))

def dijkstra(S):
    
    INF = 10 ** 15
    dist = [INF] * (N+1)
    cnt_path = [0] * (N+1)
    cnt_path[S] = 1
    q = [(0,S)] 
    while q:
        d,x = heappop(q)
        if dist[x] < d:
            continue
        dist[x] = d
        for y,dy in graph[x]:
            if d + dy > dist[y]:
                continue
            if d + dy == dist[y]:
                cnt_path[y] += cnt_path[x]
                continue
            cnt_path[y] = cnt_path[x]
            dist[y] = d + dy
            heappush(q,(d+dy,y))
    return dist, cnt_path

dist_S, cnt_path_S = dijkstra(S)
dist_T, cnt_path_T = dijkstra(T)

D = dist_S[T]

answer = cnt_path_S[T] * cnt_path_T[S]

for x in range(1,N+1):
    if dist_S[x] == dist_T[x]:
        t = cnt_path_S[x] * cnt_path_T[x] % MOD
        answer -= t * t % MOD

for x in range(1,N+1):
    for y,d in graph[x]:
        dx = dist_S[x]
        dy = dist_T[y]
        if dx + dy + d == D and dx + d > dy and dy + d > dx:
            t = cnt_path_S[x] * cnt_path_T[y] % MOD
            answer -= t * t % MOD

answer %= MOD
print(answer)

Solution (Python):
from collections import defaultdict, deque
from heapq import heappop, heappush
def inpl(): return list(map(int, input().split()))

N, M = inpl()
S, T = inpl()
G = [[] for _ in range(N+1)]
MOD = 10**9 + 7
for _ in range(M):
    u, v, d = inpl()
    G[u].append((v, d))
    G[v].append((u, d))

cost = [10**18]*(N+1)
appear = [0]*(N+1)

Q = [(0, S)]
cost[S] = 0
appear[S] = 1

while Q:
    c, p = heappop(Q)
    if cost[p] < c:
        continue
    for q, d in G[p]:
        if cost[p] + d == cost[q]:
            appear[q] = (appear[q] + appear[p])%MOD
        elif cost[p] + d < cost[q]:
            cost[q] = cost[p] + d
            appear[q] = appear[p] % MOD
            heappush(Q, (cost[q], q))

end = cost[T]*1
ans = pow(appear[T], 2, MOD)

Q = [(0, T)]
cost2 = [10**18]*(N+1)
cost2[T] = 0
appear2 = [0]*(N+1)
appear2[T] = 1
while Q:
    c, p = heappop(Q)
    if cost2[p] < c:
        continue
    for q, d in G[p]:
        if cost[p] - d != cost[q]:
            continue
        if cost2[p] + d == cost2[q]:
            appear2[q] = (appear2[q] + appear2[p])%MOD
        elif cost2[p] + d < cost2[q]:
            cost2[q] = cost2[p] + d
            appear2[q] = appear2[p] % MOD
            heappush(Q, (cost2[q], q))

Q = deque([S])
searched = [False]*(N+1)
searched[S] = True

while Q:
    p = Q.pop()
    if 2*cost2[p] == end:
        ans = (ans - (appear[p]*appear2[p])**2)%MOD

    for q, d in G[p]:
        if cost2[p] - d != cost2[q]:
            continue
        if 2*cost2[q] < end < 2*cost2[p]:
            ans = (ans - (appear[p]*appear2[q])**2)%MOD
        if not searched[q]:
            Q.append(q)
            searched[q] = True
print(ans)

Solution (Python):
from heapq import heappush, heappop
import sys
input = sys.stdin.readline

INF = float('inf')
MOD = 10**9 + 7

def Dijkstra(adjList, vSt):
    numV = len(adjList)
    numUsed = 0
    costs = [INF] * numV
    costs[vSt] = 0
    nums = [0] * numV
    nums[vSt] = 1
    PQ = []
    heappush(PQ, (0, vSt))
    while PQ:
        cNow, vNow = heappop(PQ)
        if cNow > costs[vNow]: continue
        numUsed += 1
        if numUsed == numV: break
        for v2, wt in adjList[vNow]:
            c2 = cNow + wt
            if c2 < costs[v2]:
                costs[v2] = c2
                nums[v2] = nums[vNow]
                heappush(PQ, (c2, v2))
            elif c2 == costs[v2]:
                nums[v2] = (nums[v2]+nums[vNow]) % MOD
    return (costs, nums)

def solve():
    N, M = map(int, input().split())
    S, T = map(int, input().split())
    S, T = S-1, T-1
    edges = [tuple(map(int, input().split())) for _ in range(M)]
    adjL = [[] for _ in range(N)]
    for u, v, d in edges:
        adjL[u-1].append((v-1, d))
        adjL[v-1].append((u-1, d))

    costSs, numSs = Dijkstra(adjL, S)
    costTs, numTs = Dijkstra(adjL, T)
    minCost = costSs[T]

    ans = (numSs[T]**2) % MOD
    for v in range(N):
        if costSs[v] == costTs[v] == minCost/2:
            k = (numSs[v]*numTs[v]) % MOD
            ans -= (k**2) % MOD
            ans %= MOD
    for u, v, d in edges:
        u, v = u-1, v-1
        if costSs[u] > costSs[v]:
            u, v = v, u
        if costSs[u] < minCost/2 and costTs[v] < minCost/2 and costSs[u]+d+costTs[v] == minCost:
            k = (numSs[u]*numTs[v]) % MOD
            ans -= (k**2) % MOD
            ans %= MOD

    print(ans)

solve()

Solution (Python):
import heapq
import sys
def input():
	return sys.stdin.readline()[:-1]

MOD = 10**9+7
class DijkstraList():
	
	
	def __init__(self, adj, start):
		self.list = adj
		self.start = start
		self.size = len(adj)

	def solve(self):
		self.dist = [float("inf") for _ in range(self.size)]
		self.dist[self.start] = 0
		self.dp = [0 for _ in range(self.size)]
		self.dp[self.start] = 1
		
		self.q = []
		self.cnt = 0

		heapq.heappush(self.q, (0, self.start))

		while self.q and self.cnt < self.size:
			u_dist, u = heapq.heappop(self.q)
			if self.dist[u] < u_dist:
				continue
			for v, w in self.list[u]:
				if self.dist[v] > u_dist + w:
					self.dist[v] = u_dist + w
					self.dp[v] = self.dp[u]
					
					heapq.heappush(self.q, (self.dist[v], v))
				elif self.dist[v] == u_dist + w:
					self.dp[v] += self.dp[u]
					self.dp[v] %= MOD

			self.cnt += 1
		
		return

	def distance(self, goal):
		return self.dist[goal]

	def get_dp(self, x):
		return self.dp[x]

n, m = map(int, input().split())
s, t = map(int, input().split())
s, t = s-1, t-1
adj = [[] for _ in range(n)]
for _ in range(m):
	u, v, d = map(int, input().split())
	adj[u-1].append([v-1, d])
	adj[v-1].append([u-1, d])
S, T = DijkstraList(adj, s), DijkstraList(adj, t)
S.solve()
T.solve()
ans = pow(S.get_dp(t), 2, MOD)
goal = S.distance(t)


for i in range(n):
	SI, TI = S.distance(i), T.distance(i)
	if SI + TI != goal:
		continue
	if SI*2 == goal:
		ans -= S.get_dp(i) * T.get_dp(i) * S.get_dp(i) * T.get_dp(i)
		ans %= MOD

	else:
		for j, d in adj[i]:
			if i>j:
				continue
			SJ, TJ = S.distance(j), T.distance(j)
			if SJ + TJ != goal or abs(SI - SJ) != d:
				continue
			if SI*2 < goal < SJ*2:
				
				ans -= S.get_dp(i) * T.get_dp(j) * S.get_dp(i) * T.get_dp(j)
				ans %= MOD
			elif SJ*2 < goal < SI*2:
				
				ans -= S.get_dp(j) * T.get_dp(i) * S.get_dp(j) * T.get_dp(i)
				ans %= MOD

print(ans)

Solution (Python):
from heapq import heapify,heappop,heappush
import sys
input=sys.stdin.readline
N,M=map(int,input().split())
S,T=map(int,input().split())
Graph=[set() for i in range(N)]
Edge=[]
mod=10**9+7
inf=float("inf")
for i in range(M):
  a,b,c=map(int,input().split())
  Graph[a-1].add((c,b-1))
  Graph[b-1].add((c,a-1))
  Edge.append((a-1,b-1,c))
  Edge.append((b-1,a-1,c))
def dijkstra(s,n,links):
  Visited=[False]*n
  cost=[inf]*n
  P=[0]*n
  P[s]=1
  cost[s]=0
  heap=[(0,s)]
  heapify(heap)
  while heap:
    hc,hp=heappop(heap)
    if Visited[hp]:
      continue
    Visited[hp]=True
    for c,p in links[hp]:
      if Visited[p]:
        continue
      if c+hc<cost[p]:
        cost[p]=c+hc
        P[p]=P[hp]
        heappush(heap,(cost[p],p))
      elif c+hc==cost[p]:
        P[p]=(P[p]+P[hp])%mod
  return cost,P

def dp(s,n,links,d):
  Visited=[False]*n
  cost=[inf]*n
  P=[0]*n
  P[s]=1
  cost[s]=0
  heap=[(0,s)]
  heapify(heap)
  while heap:
    hc,hp=heappop(heap)
    if 2*hc>d:
      break
    if Visited[hp]:
      continue
    Visited[hp]=True
    for c,p in links[hp]:
      if Visited[p]:
        continue
      if c+hc<cost[p]:
        cost[p]=c+hc
        P[p]=P[hp]
        heappush(heap,(cost[p],p))
      elif c+hc==cost[p]:
        P[p]=(P[p]+P[hp])%mod
  return cost,P

cost_S,P_S=dijkstra(S-1,N,Graph)
D=cost_S[T-1]
cost_T,P_T=dp(T-1,N,Graph,D)
Pat=P_S[T-1]
ans=(Pat**2)%mod
for a,b,c in Edge:
  if cost_S[a]+cost_T[b]+c==D:
    if 2*cost_S[a]<D and 2*cost_T[b]<D:
      ans-=((P_S[a]*P_T[b])**2)%mod
for i in range(N):
  if 2*cost_S[i]==2*cost_T[i]==D:
    ans-=((P_S[i]*P_T[i])**2)%mod
print(ans%mod)

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n,m = LI()
    s,t = LI()
    e = collections.defaultdict(list)
    for _ in range(m):
        u,v,d = LI()
        e[u].append((v,d))
        e[v].append((u,d))

    def search(s,t):
        d = collections.defaultdict(lambda: inf)
        dc = collections.defaultdict(int)
        d[s] = 0
        dc[s] = 1
        q = []
        heapq.heappush(q, (0, s))
        v = collections.defaultdict(bool)
        while len(q):
            k, u = heapq.heappop(q)
            if v[u]:
                continue
            v[u] = True

            dc[u] %= mod
            if u == t:
                return (d,dc)

            uc = dc[u]

            for uv, ud in e[u]:
                if v[uv]:
                    continue
                vd = k + ud
                if d[uv] < vd:
                    continue
                if d[uv] > vd:
                    d[uv] = vd
                    dc[uv] = uc
                    heapq.heappush(q, (vd, uv))
                elif d[uv] == vd:
                    dc[uv] += uc

        return (d,dc)

    d1,dc1 = search(s,t)
    d2,dc2 = search(t,s)
    rd = d1[t]
    kk = rd / 2.0

    r = dc1[t] ** 2 % mod
    for k in list(d1.keys()):
        t = d1[k]
        c = dc1[k]
        if t > kk:
            continue
        if t == kk:
            if d2[k] == t:
                r -= pow(c,2,mod) * pow(dc2[k],2,mod) % mod
            continue

        for uv, ud in e[k]:
            if d2[uv] >= kk or t + ud + d2[uv] != rd:
                continue
            r -= pow(c,2,mod) * pow(dc2[uv],2,mod) % mod

    return r % mod




print(main())

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n,m = LI()
    s,t = LI()
    e = collections.defaultdict(list)
    for _ in range(m):
        u,v,d = LI()
        e[u].append((v,d))
        e[v].append((u,d))

    def search(s,t):
        d = collections.defaultdict(lambda: inf)
        dc = collections.defaultdict(int)
        d[s] = 0
        dc[s] = 1
        q = []
        heapq.heappush(q, (0, s))
        v = collections.defaultdict(bool)
        while len(q):
            k, u = heapq.heappop(q)
            if v[u]:
                continue
            v[u] = True

            dc[u] %= mod
            if u == t:
                return (d,dc)

            uc = dc[u]

            for uv, ud in e[u]:
                if v[uv]:
                    continue
                vd = k + ud
                if d[uv] < vd:
                    continue
                if d[uv] > vd:
                    d[uv] = vd
                    dc[uv] = uc
                    heapq.heappush(q, (vd, uv))
                elif d[uv] == vd:
                    dc[uv] += uc

        return (d,dc)

    d1,dc1 = search(s,t)
    d2,dc2 = search(t,s)
    rd = d1[t]
    kk = rd / 2.0

    r = dc1[t] ** 2 % mod
    for k in list(d1.keys()):
        t = d1[k]
        c = dc1[k]
        if t > kk:
            continue
        if t == kk:
            if d2[k] == t:
                r -= pow(c,2,mod) * pow(dc2[k],2,mod) % mod
            continue

        for uv, ud in e[k]:
            if d2[uv] >= kk or t + ud + d2[uv] != rd:
                continue
            r -= pow(c,2,mod) * pow(dc2[uv],2,mod) % mod

    return r % mod




print(main())

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n,m = LI()
    s,t = LI()
    e = collections.defaultdict(list)
    for _ in range(m):
        u,v,d = LI()
        e[u].append((v,d))
        e[v].append((u,d))

    def search(s,t):
        d = collections.defaultdict(lambda: inf)
        dc = collections.defaultdict(int)
        d[s] = 0
        dc[s] = 1
        q = []
        heapq.heappush(q, (0, s))
        v = collections.defaultdict(bool)
        while len(q):
            k, u = heapq.heappop(q)
            if v[u]:
                continue
            v[u] = True

            if u == t:
                return (d,dc)

            uc = dc[u]

            for uv, ud in e[u]:
                if v[uv]:
                    continue
                vd = k + ud
                if d[uv] < vd:
                    continue
                if d[uv] > vd:
                    d[uv] = vd
                    dc[uv] = uc
                    heapq.heappush(q, (vd, uv))
                elif d[uv] == vd:
                    dc[uv] += uc

        return (d,dc)

    d1,dc1 = search(s,t)
    d2,dc2 = search(t,s)
    rd = d1[t]
    kk = rd / 2.0

    r = dc1[t] ** 2 % mod
    for k in list(d1.keys()):
        t = d1[k]
        c = dc1[k]
        if t > kk:
            continue
        if t == kk:
            if d2[k] == t:
                r -= (c**2 * dc2[k]**2) % mod
            continue

        for uv, ud in e[k]:
            if d2[uv] >= kk or t + ud + d2[uv] != rd:
                continue
            r -= (c**2 * dc2[uv]**2) % mod

    return r % mod




print(main())

Solution (Python):


import sys,collections,itertools,re,math,fractions,decimal,random,array,bisect,heapq


sys.setrecursionlimit(100000)
MOD = 10**9 + 7

def dijkstra(g, v_size, st):
    """Dijkstra method"""

    d = [float('inf')] * v_size
    d[st] = 0

    que = [(0, st)]     
    heapq.heapify(que)

    while len(que) > 0:
        dist, v0 = heapq.heappop(que)
        if d[v0] < dist: continue
        for v1, dist in g[v0]:
            if d[v1] > d[v0] + dist:
                d[v1] = d[v0] + dist
                heapq.heappush(que, (d[v1], v1))

    return d

def solve(f):
    n, m = f.read_int_list()
    s, t = f.read_int_list()
    g = collections.defaultdict(list)
    for m in xrange(m):
        u, v, d = f.read_int_list()
        g[u-1].append((v-1, d))
        g[v-1].append((u-1, d))

    d = dijkstra(g, n, s-1)

    cand = []
    for i in xrange(n):
        if d[t-1] == d[i]*2: cand.append((i,i))

    memo = [0]*n
    memo[t-1] = 1
    ary = [(d[i], i) for i in xrange(n)]
    ary.sort(reverse=True)

    for _, v in ary:
        for vi, di in g[v]:
            if d[v]-di == d[vi]:
                memo[vi] = (memo[vi]+memo[v])%MOD
                if d[vi]*2 < d[t-1] < (d[vi] + di)*2:
                    cand.append((vi, v))

    memo2 = [0]*n
    memo2[s-1] = 1

    ary.sort()
    for _, v in ary:
        for vi, di in g[v]:
            if d[v]+di == d[vi]:
                memo2[vi] = (memo2[vi]+memo2[v])%MOD

    cnt = []
    for v1, v2 in cand:
        cnt.append((memo2[v1]*memo[v2])%MOD)

    s = sum(cnt)
    ans = 0
    for ci in cnt:
        ans = (ans+ci*(s-ci))%MOD

    return ans


class Reader(object):
    def __init__(self, filename=None):
        self.file = open(filename) if filename is not None else None
        self.case = 1

    def __readline(self):
        return self.file.next().strip() if self.file else raw_input()

    def next_case(self):
        self.file.next()
        self.case += 1

    def read_int(self): return int(self.__readline())
    def read_float(self): return float(self.__readline())
    def read_long(self): return long(self.__readline())
    def read_decimal(self): return decimal.Decimal(self.__readline())
    def read_str(self): return self.__readline()

    def read_int_list(self): return map(int, self.__readline().split())
    def read_float_list(self): return map(float, self.__readline().split())
    def read_long_list(self): return map(long, self.__readline().split())
    def read_decimal_list(self): return map(decimal.Decimal, self.__readline().split())
    def read_str_list(self): return self.__readline().split()

if __name__ == '__main__':
    filename = sys.argv[1] if len(sys.argv) > 1 else None
    f = Reader(filename)
    if f.file:
        while True:
            print "Case 
            try:
                f.next_case()
            except StopIteration:
                break
    else:
        print solve(f)

Solution (Python):
from heapq import heappush, heappop
from collections import deque
N, M = map(int, input().split())
S, T = map(int, input().split())
S -= 1; T -= 1
G = [[] for i in range(N)]
for i in range(M):
    u, v, d = map(int, input().split())
    G[u-1].append((v-1, d))
    G[v-1].append((u-1, d))

INF = 10**18
MOD = 10**9 + 7

dist = [INF]*N
dist[S] = 0
que = [(0, S)]
while que:
    cost, v = heappop(que)
    if dist[v] < cost:
        continue
    for w, c in G[v]:
        if cost + c < dist[w]:
            dist[w] = cost + c
            heappush(que, (cost + c, w))

que = deque()
mi = [0]*N
que.append(T)
mi[T] = 1
H = [[] for i in range(N)]
RH = [[] for i in range(N)]

D = dist[T]

V = []
E = []
while que:
    v = que.popleft()
    d = dist[v]
    if d*2 == D:
        V.append(v)
    for w, c in G[v]:
        if d == dist[w] + c:
            if dist[w]*2 < D < d*2:
                E.append((w, v))
            H[v].append(w)
            RH[w].append(v)
            if mi[w] == 0:
                mi[w] = 1
                que.append(w)

def bfs(G, RG, s):
    *deg, = map(len, RG)
    que = deque([s])
    res = [0]*N
    res[s] = 1
    while que:
        v = que.popleft()
        for w in G[v]:
            deg[w] -= 1
            res[w] += res[v]
            if deg[w] == 0:
                que.append(w)
                res[w] %= MOD
    return res
CT = bfs(H, RH, T)
CS = bfs(RH, H, S)

ans = CS[T] * CT[S] % MOD
for v in V:
    ans -= (CS[v] * CT[v]) ** 2
    ans %= MOD
for v, w in E:
    ans -= (CS[v] * CT[w]) ** 2
    ans %= MOD
print(ans)

Solution (Python):
from heapq import heappop, heappush
import sys

MOD, INF = 1000000007, float('inf')


def solve(s, t, links):
    q = [(0, 0, s, -1, 0), (0, 0, t, -1, 1)]
    visited_fwd, visited_bwd = [INF] * n, [INF] * n
    patterns_fwd, patterns_bwd = [0] * (n + 1), [0] * (n + 1)
    patterns_fwd[-1] = patterns_bwd[-1] = 1
    collision_nodes, collision_links = set(), set()
    limit = 0
    while q:
        cost, cost_a, v, a, is_bwd = heappop(q)
        if is_bwd:
            visited_self = visited_bwd
            visited_opp = visited_fwd
            patterns_self = patterns_bwd
        else:
            visited_self = visited_fwd
            visited_opp = visited_bwd
            patterns_self = patterns_fwd

        relax_flag = False
        cost_preceding = visited_self[v]
        if cost_preceding == INF:
            visited_self[v] = cost
            relax_flag = True
        elif cost > cost_preceding:
            continue

        patterns_self[v] += patterns_self[a]

        cost_opp = visited_opp[v]
        if cost_opp != INF:
            limit = cost + cost_opp
            if cost == cost_opp:
                collision_nodes.add(v)
            else:
                collision_links.add((v, a) if is_bwd else (a, v))
            break

        if relax_flag:
            for u, du in links[v].items():
                nc = cost + du
                if visited_self[u] < nc:
                    continue
                heappush(q, (nc, cost, u, v, is_bwd))

    collision_time = limit / 2
    while q:
        cost, cost_a, v, a, is_bwd = heappop(q)
        if cost > limit:
            break
        visited_self = visited_bwd if is_bwd else visited_fwd
        if visited_self[v] == INF:
            visited_self[v] = cost
        if is_bwd:
            if cost == collision_time:
                patterns_bwd[v] += patterns_bwd[a]
            continue
        if cost_a == collision_time:
            collision_nodes.add(a)
        elif cost == collision_time:
            collision_nodes.add(v)
            patterns_fwd[v] += patterns_fwd[a]
        else:
            collision_links.add((a, v))

    shortest_count = 0
    collision_count = 0

    for v in collision_nodes:
        if visited_fwd[v] == visited_bwd[v]:
            r = patterns_fwd[v] * patterns_bwd[v]
            shortest_count += r
            shortest_count %= MOD
            collision_count += r * r
            collision_count %= MOD

    for u, v in collision_links:
        if visited_fwd[u] + visited_bwd[v] + links[u][v] == limit:
            r = patterns_fwd[u] * patterns_bwd[v]
            shortest_count += r
            shortest_count %= MOD
            collision_count += r * r
            collision_count %= MOD

    return (shortest_count ** 2 - collision_count) % MOD


n, m = map(int, input().split())
s, t = map(int, input().split())
s -= 1
t -= 1
links = [{} for _ in range(n)]
for uvd in sys.stdin.readlines():
    u, v, d = map(int, uvd.split())
    
    
    u -= 1
    v -= 1
    links[u][v] = d
    links[v][u] = d
print(solve(s, t, links))