#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int Maxm = 25;
const ll P = 1e9 + 7;
int mm, OK, must[Maxm], f[Maxm], ok[Maxm], pd[Maxm];
int A, B, c[Maxm], p[Maxm], C[Maxm][Maxm];
vector<int> edge[Maxm];
ll c1[Maxm], c2[Maxm], inv[Maxm];
int _norm(int x)
{
    return x > P ? x - P : x;
}
void adde(int x, int y)
{
    edge[x].push_back(y);
    edge[y].push_back(x);
}
void dfs(int x, int now)
{
    if(!OK) return;
    c[x] = now;
    if(now == 1) A++;
    else B++;
    for(auto &y : edge[x])
    {
        if(must[y] || ok[y]) continue;
        if(c[y])
        {
            if(c[y] != 3 - now)
            {
                OK = 0;
                return;
            }
            continue;
        }
        dfs(y, 3 - now);
        if(!OK) return;
    }
}
ll C1(int x)
{
    if(x < 0) return 0;
    return c1[x];
}
ll C2(int x)
{
    if(x < 0) return 0;
    return c2[x];
}
void solve()
{
    int n1, n2, m, k;
    cin >> n1 >> n2 >> m >> k;
    vector<pair<int, int>> limit;
    int need = 0;
    for(int i = 1; i <= k; i++)
    {
        int x, y;
        cin >> x >> y;
        pd[x] = pd[y] = 1;
        if(x == y)
        {
            if(must[x]) continue;
            must[x] = 1;
            need++;
        }
        else
        {
            limit.push_back({x, y});
        }
    }
    int num = 0;
    mm = 0;
    for(int i = 1; i <= m; i++)
    {
        if(!pd[i])
        {
            num++;
            need++;
            continue;
        }
        if(must[i]) continue;
        f[++mm] = i;
    }
    for(auto &[x, y] : limit)
    {
        adde(x, y);
    }
    int Lim = 1 << mm;
    int ans = 0;
    c1[0] = c2[0] = inv[1] = 1;
    for(int i = 2; i <= m; i++) inv[i] = P - P / i * inv[P % i] % P;
    for(int i = 1; i <= m; i++)
    {
        c1[i] = c1[i - 1] * (n1 + num - i) % P * inv[i] % P;
        c2[i] = c2[i - 1] * (n2 + num - i) % P * inv[i] % P;
    }
    for(int i = 0; i <= m; i++) C[i][0] = 1;
    for(int i = 1; i <= m; i++)
    {
        for(int j = 1; j <= m; j++)
        {
            C[i][j] = _norm(C[i - 1][j] + C[i - 1][j - 1]);
        }
    }
    for(int t = 0; t < Lim; t++)
    {
        int sp = 0;
        OK = 1;
        A = B = 0;
        int now = 0, tmp = need, z = mm;
        for(int i = 1; i <= mm; i++)
        {
            p[i] = c[f[i]] = 0;
            if(t & (1 << (i - 1)))
            {
                ok[f[i]] = 1;
                tmp++;
                z--;
            }
            else ok[f[i]] = 0;
        }
        p[0] = 1;
        for(int i = 1; i <= mm; i++)
        {
            if(!ok[f[i]] && !c[f[i]])
            {
                dfs(f[i], 1);
                if(!OK) break;
                if(A > B) swap(A, B);
                if(A == 0 && B == 1)
                {
                    sp++;
                    continue;
                }
                now = now + B;
                for(int i = now; i >= A; i--)
                {
                    p[i] = p[i - A];
                    if(i >= B) p[i] += p[i - B];
                }
                for(int i = A - 1; i >= 0; i--) p[i] = 0;
                A = B = 0;
            }
        }
        if(!OK) continue;
        now += sp;
        for(int k = now; k >= 0; k--)
        {
            int x = 0;
            for(int i = 0; i <= sp; i++)
            {
                x = _norm(x + p[k - i] * (ll)C[sp][i] % P);
            }
            p[k] = x;
        }
        for(int i = 0; i <= now; i++)
        {
            if(!p[i]) continue;
            ans = _norm(ans + C1(tmp + i - 1) * C2(tmp + (z - i) - 1) % P * p[i] % P);
        }
    }
    cout << ans << "\n";
    for(int i = 1; i <= m; i++)
    {
        must[i] = pd[i] = 0;
        vector<int>().swap(edge[i]);
    }
}
int main()
{
    cin.tie(0)->sync_with_stdio(false);
    int T;
    cin >> T;
    while(T--) solve();
    return 0;
}