Skip to content

Commit 3c596ea

Browse files
lkirklkirk
committed
revert changes to r2_ij, perhaps revisit these at another point
1 parent c5bdbc7 commit 3c596ea

File tree

3 files changed

+54
-56
lines changed

3 files changed

+54
-56
lines changed

‎c/tests/test_stats.c‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2637,11 +2637,10 @@ test_paper_ex_two_site(void)
26372637

26382638
tsk_treeseq_from_text(&ts, 10, paper_ex_nodes, paper_ex_edges, NULL, paper_ex_sites,
26392639
paper_ex_mutations, paper_ex_individuals, NULL, 0);
2640-
double truth_three_index_tuples[27] = { 1, 1, 0, 0.1111111111111111,
2641-
0.1111111111111111, 0, 0.1111111111111111, 0.1111111111111111, 0,
2642-
0.1111111111111111, 0.1111111111111111, 0, 1, 1, 0.94921874999999978, 1, 1,
2643-
0.94921874999999978, 0.1111111111111111, 0.1111111111111111, 0, 1, 1,
2644-
0.94921874999999978, 1, 1, 0.94921874999999978 };
2640+
double truth_three_index_tuples[27] = { 1, 1, NAN, 0.1111111111111111,
2641+
0.1111111111111111, NAN, 0.1111111111111111, 0.1111111111111111, NAN,
2642+
0.1111111111111111, 0.1111111111111111, NAN, 1, 1, 1, 1, 1, 1,
2643+
0.1111111111111111, 0.1111111111111111, NAN, 1, 1, 1, 1, 1, 1 };
26452644

26462645
tsk_size_t sample_set_sizes[3], num_index_tuples;
26472646
tsk_id_t sample_sets[ts.num_samples * 3], index_tuples[2 * 3] = { 0, 1, 0, 0, 0, 2 };

‎c/tskit/trees.c‎

Lines changed: 19 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -4944,35 +4944,33 @@ r2_ij_summary_func(tsk_size_t TSK_UNUSED(state_dim), const double *state,
49444944
const double *state_row;
49454945
tsk_size_t k;
49464946
tsk_id_t i, j;
4947-
double ni, w_AB_i, w_Ab_i, w_aB_i, w_A_i, w_B_i, D_i;
4948-
double nj, w_AB_j, w_Ab_j, w_aB_j, w_A_j, w_B_j, D_j;
4949-
double p_A, p_B;
4947+
double n, pAB, pAb, paB, pA, pB, D_i, D_j, denom_i, denom_j;
49504948

49514949
for (k = 0; k < result_dim; k++) {
49524950
i = args.set_indexes[2 * k];
49534951
j = args.set_indexes[2 * k + 1];
49544952

4955-
ni = (double) args.sample_set_sizes[i];
4953+
n = (double) args.sample_set_sizes[i];
49564954
state_row = GET_2D_ROW(state, 3, i);
4957-
w_AB_i = state_row[0];
4958-
w_Ab_i = state_row[1];
4959-
w_aB_i = state_row[2];
4960-
w_A_i = w_AB_i + w_Ab_i;
4961-
w_B_i = w_AB_i + w_aB_i;
4962-
D_i = (ni * w_AB_i - (w_A_i * w_B_i)) / (ni * ni);
4955+
pAB = state_row[0] / n;
4956+
pAb = state_row[1] / n;
4957+
paB = state_row[2] / n;
4958+
pA = pAB + pAb;
4959+
pB = pAB + paB;
4960+
D_i = pAB - (pA * pB);
4961+
denom_i = sqrt(pA * (1 - pA) * pB * (1 - pB));
49634962

4964-
nj = (double) args.sample_set_sizes[j];
4963+
n = (double) args.sample_set_sizes[j];
49654964
state_row = GET_2D_ROW(state, 3, j);
4966-
w_AB_j = state_row[0];
4967-
w_Ab_j = state_row[1];
4968-
w_aB_j = state_row[2];
4969-
w_A_j = w_AB_j + w_Ab_j;
4970-
w_B_j = w_AB_j + w_aB_j;
4971-
D_j = (nj * w_AB_j - (w_A_j * w_B_j)) / (nj * nj);
4972-
4973-
p_A = (w_A_i + w_A_j) / (ni + nj);
4974-
p_B = (w_B_i + w_B_j) / (ni + nj);
4975-
result[k] = (D_i * D_j) / (p_A * (1 - p_A) * p_B * (1 - p_B));
4965+
pAB = state_row[0] / n;
4966+
pAb = state_row[1] / n;
4967+
paB = state_row[2] / n;
4968+
pA = pAB + pAb;
4969+
pB = pAB + paB;
4970+
D_j = pAB - (pA * pB);
4971+
denom_j = sqrt(pA * (1 - pA) * pB * (1 - pB));
4972+
4973+
result[k] = (D_i * D_j) / (denom_i * denom_j);
49764974
}
49774975
return 0;
49784976
}

‎python/tests/test_ld_matrix.py‎

Lines changed: 31 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -277,6 +277,7 @@ def norm_hap_weighted_ij(
277277
wAB_i = hap_weights[0, i]
278278
wAB_j = hap_weights[0, j]
279279
result[k] = (wAB_i + wAB_j) / (ni + nj)
280+
# result[k] = (wAB_i / ni / 2) + (wAB_j / nj / 2)
280281

281282

282283
def norm_total_weighted(
@@ -1034,26 +1035,26 @@ def r2_ij_summary_func(
10341035
for k in range(result_dim):
10351036
i = set_indexes[k][0]
10361037
j = set_indexes[k][1]
1037-
ni = sample_set_sizes[i]
1038-
w_AB_i = state[0, i]
1039-
w_Ab_i = state[1, i]
1040-
w_aB_i = state[2, i]
1041-
w_A_i = w_AB_i + w_Ab_i
1042-
w_B_i = w_AB_i + w_aB_i
1043-
D_i = (ni * w_AB_i - (w_A_i * w_B_i)) / (ni * ni)
1038+
n = sample_set_sizes[i]
1039+
pAB = state[0, i] / n
1040+
pAb = state[1, i] / n
1041+
paB = state[2, i] / n
1042+
pA = pAB + pAb
1043+
pB = pAB + paB
1044+
D_i = pAB - pA * pB
1045+
denom_i = np.sqrt(pA * (1 - pA) * pB * (1 - pB))
1046+
1047+
n = sample_set_sizes[j]
1048+
pAB = state[0, j] / n
1049+
pAb = state[1, j] / n
1050+
paB = state[2, j] / n
1051+
pA = pAB + pAb
1052+
pB = pAB + paB
1053+
D_j = pAB - pA * pB
1054+
denom_j = np.sqrt(pA * (1 - pA) * pB * (1 - pB))
10441055

1045-
nj = sample_set_sizes[j]
1046-
w_AB_j = state[0, j]
1047-
w_Ab_j = state[1, j]
1048-
w_aB_j = state[2, j]
1049-
w_A_j = w_AB_j + w_Ab_j
1050-
w_B_j = w_AB_j + w_aB_j
1051-
D_j = (nj * w_AB_j - (w_A_j * w_B_j)) / (nj * nj)
1052-
1053-
p_A = (w_A_i + w_A_j) / (ni + nj)
1054-
p_B = (w_B_i + w_B_j) / (ni + nj)
10551056
with suppress_overflow_div0_warning():
1056-
result[k] = (D_i * D_j) / (p_A * (1 - p_A) * p_B * (1 - p_B))
1057+
result[k] = (D_i * D_j) / (denom_i * denom_j)
10571058

10581059

10591060
def D_summary_func(
@@ -2317,12 +2318,12 @@ def test_two_way_site_ld_matrix(ts, stat):
23172318
(
23182319
correlated,
23192320
(np.array([0, 1, 2, 3, 4, 5]), np.array([6, 7, 8, 9, 10])),
2320-
np.float64(0.9708352229780801),
2321+
np.float64(1.0),
23212322
),
23222323
(
23232324
correlated,
23242325
(np.array([0, 1, 2, 3, 4, 5]), np.array([6, 7, 8, 9])),
2325-
np.float64(0.9526958931720837),
2326+
np.float64(1.0),
23262327
),
23272328
(
23282329
correlated,
@@ -2332,12 +2333,12 @@ def test_two_way_site_ld_matrix(ts, stat):
23322333
(
23332334
correlated,
23342335
(np.array([0, 1, 2, 3, 4, 5]), np.array([6, 7])),
2335-
np.float64(0.7585185185185186),
2336+
np.float64(np.nan),
23362337
),
23372338
(
23382339
correlated,
23392340
(np.array([0, 1, 2, 3, 4, 5]), np.array([6])),
2340-
np.float64(0.0),
2341+
np.float64(np.nan),
23412342
),
23422343
(
23432344
anticorrelated := np.array(
@@ -2355,37 +2356,37 @@ def test_two_way_site_ld_matrix(ts, stat):
23552356
(
23562357
anticorrelated,
23572358
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3, 5, 7, 9, 11, 13])),
2358-
np.float64(0.9798566895766568),
2359+
np.float64(1.0),
23592360
),
23602361
(
23612362
anticorrelated,
23622363
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3, 5, 7, 9, 11])),
2363-
np.float64(0.8574999999999999),
2364+
np.float64(np.nan),
23642365
),
23652366
(
23662367
anticorrelated,
23672368
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3, 5, 7, 9])),
2368-
np.float64(0.8299777777777777),
2369+
np.float64(np.nan),
23692370
),
23702371
(
23712372
anticorrelated,
23722373
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3, 5, 7])),
2373-
np.float64(0.6328124999999999),
2374+
np.float64(np.nan),
23742375
),
23752376
(
23762377
anticorrelated,
23772378
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3, 5])),
2378-
np.float64(0.57179616638322),
2379+
np.float64(np.nan),
23792380
),
23802381
(
23812382
anticorrelated,
23822383
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1, 3])),
2383-
np.float64(0.0),
2384+
np.float64(np.nan),
23842385
),
23852386
(
23862387
anticorrelated,
23872388
(np.array([0, 2, 4, 6, 8, 10, 12, 14]), np.array([1])),
2388-
np.float64(0.0),
2389+
np.float64(np.nan),
23892390
),
23902391
],
23912392
)
@@ -2404,4 +2405,4 @@ def test_multipopulation_r2_varying_unequal_set_sizes(genotypes, sample_sets, ex
24042405
r2_ij_summary_func(state_dim, state, 1, result[i, j], params)
24052406
norm_hap_weighted_ij(1, state, max(a) + 1, max(b) + 1, norm[i, j], params)
24062407

2407-
np.testing.assert_allclose(expected, (result * norm).sum())
2408+
np.testing.assert_allclose((result * norm).sum(), expected)

0 commit comments

Comments
 (0)