diff --git a/AGENTS.md b/AGENTS.md index 99ff52fc..994c5871 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -44,7 +44,7 @@ when ready. | `spr.cpp` | 7 KB | SPR distance approximation | | `spr_lookup.cpp` | — | SPR lookup-table implementation | | `nni_distance.cpp` | 16 KB | NNI distance approximations; HybridBuffer allocation | -| `li_diameters.h` | 30 KB | Precomputed NNI diameter lookup tables | +| `li_diameters.h` | <1 KB | Exact NNI diameters for trees with ≤ 12 leaves | | `information.h` | 6 KB | log₂ / factorial lookup tables (max 8192); cached at startup | | `binary_entropy_counts.cpp` | — | Binary entropy calculations | | `day_1985.cpp` | 10 KB | Consensus tree information | diff --git a/NEWS.md b/NEWS.md index 4db05c58..ae53d755 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,17 @@ ## Bug fixes +- `NNIDist()`: `best_lower` was not a lower bound. It added the diameter of + each region of unmatched edges, which bounds the distance from above; it + could exceed the true distance (#220). `lower` added a lower bound on the + diameter, which exceeded the number of unmatched edges in regions of 37 or + more edges. Both now count unmatched edges, the only known lower bound, and + are equal. + +- `NNIDist()` and `NNIDiameter()` use theorem 3 of Fack _et al._ (2002) for + `fack_upper` and `fackMax`, tightening `loose_upper` and `best_upper`. + All seven fields of `NNIDist()` are now documented. + - Fix `MASTSize()` under-counting on some trees: the root of each tree was never compared with the subtrees of the other, so a maximum agreement subtree confined to one side of a root could be missed. Affects both diff --git a/R/tree_distance_nni.R b/R/tree_distance_nni.R index d47d4ac6..432021df 100644 --- a/R/tree_distance_nni.R +++ b/R/tree_distance_nni.R @@ -31,21 +31,46 @@ #' \tab 15 \tab 18 \tab 21 \tab ? \cr #' } #' +#' Upper bounds treat each region separately: transforming each region in turn +#' is a valid path between the trees, so the sum of the diameters of the regions +#' bounds the distance from above. +#' A diameter is the distance between the _most_ different pair of trees with +#' that many leaves, so it says nothing about how close a particular pair might +#' be: the only known lower bound is the number of unmatched edges +#' \insertCite{Li1996}{TreeDist}. +#' An optimal path need not leave matched edges untouched +#' \insertCite{Li1996}{TreeDist}, so even the exact distance between the trees +#' that each region induces is not a lower bound. +#' #' @template tree12Params #' #' @return `NNIDist()` returns, for each pair of trees, a named vector -#' containing three integers: -#' -#' - `lower` is a lower bound on the NNI distance, and corresponds -#' to the RF distance between the trees. -#' +#' containing seven integers: +#' +#' - `lower` is a lower bound on the NNI distance: the number of edges in +#' one tree that do not occur in the other, i.e. half the Robinson-Foulds +#' distance between the trees. +#' +#' - `best_lower` is the tightest available lower bound. +#' No bound tighter than `lower` is known, so the two values are equal. +#' #' - `tight_upper` is an upper bound on the distance, based on calculated #' maximum diameters for trees with < 13 leaves. _NA_ is returned if trees are #' too different to employ this approach. -#' +#' +#' - `best_upper` is the tightest available upper bound: it uses the exact +#' diameter for regions with < 13 leaves, and `loose_upper` for larger +#' regions. +#' #' - `loose_upper` is a looser upper bound on the distance, using _n_ log _n_ + -#' O(_n_). -#' +#' O(_n_): the smaller of `fack_upper` and `li_upper`. +#' +#' - `fack_upper` is an upper bound based on the maximum diameter given by +#' theorem 3 of \insertCite{Fack2002;textual}{TreeDist}. +#' +#' - `li_upper` is an upper bound based on the sorting approach of +#' \insertCite{Li1996;textual}{TreeDist}. +#' #' #' @references #' \insertAllCited{} @@ -108,9 +133,9 @@ NNIDist <- function(tree1, tree2 = tree1) { #' - `liMax`: Upper bound on diameter following Li _et al._ (1996), i.e. #' \deqn{n \log2{n} + \textrm{O}(n)}{n log₂ _n_ + O(_n_)}; #' -#' - `fackMax`: Upper bound on diameter following Fack _et al_. (2002), i.e. -#' (\deqn{N - 2}{_N_ − 2}) ceiling(\deqn{\log2{n}}{log₂ _N_}) -#' + _N_; +#' - `fackMax`: Upper bound on diameter following theorem 3 of +#' Fack _et al_. (2002), i.e. _N_ _c_ − 2^_c_ + 1 + 2(_N_ − +#' ceiling(log₂ (_N_ + 1))), where _c_ = ceiling(log₂ _N_); #' #' - `max`: The smaller of `liMax` and `fackMax`; #' @@ -145,7 +170,8 @@ NNIDiameter.numeric <- function(tree) { fackMin <- ceiling(lfactorial(n) / 4L / log(2)) liMax <- .SortingNumber(tree) + (2 * .DegenerateDistance(tree)) - fackMax <- ((n - 2L) * ceiling(log2(n))) + n + lgN <- ceiling(log2(n)) + fackMax <- (n * lgN) - (2 ^ lgN) + 1 + (2 * (n - ceiling(log2(n + 1)))) cbind(liMin = liMin, fackMin = fackMin, min = pmax(liMin, fackMin), diff --git a/man/NNIDist.Rd b/man/NNIDist.Rd index 75fadd1d..d782b683 100644 --- a/man/NNIDist.Rd +++ b/man/NNIDist.Rd @@ -18,15 +18,25 @@ or an integer specifying the number of leaves in a tree/trees.} } \value{ \code{NNIDist()} returns, for each pair of trees, a named vector -containing three integers: +containing seven integers: \itemize{ -\item \code{lower} is a lower bound on the NNI distance, and corresponds -to the RF distance between the trees. +\item \code{lower} is a lower bound on the NNI distance: the number of edges in +one tree that do not occur in the other, i.e. half the Robinson-Foulds +distance between the trees. +\item \code{best_lower} is the tightest available lower bound. +No bound tighter than \code{lower} is known, so the two values are equal. \item \code{tight_upper} is an upper bound on the distance, based on calculated maximum diameters for trees with < 13 leaves. \emph{NA} is returned if trees are too different to employ this approach. +\item \code{best_upper} is the tightest available upper bound: it uses the exact +diameter for regions with < 13 leaves, and \code{loose_upper} for larger +regions. \item \code{loose_upper} is a looser upper bound on the distance, using \emph{n} log \emph{n} + -O(\emph{n}). +O(\emph{n}): the smaller of \code{fack_upper} and \code{li_upper}. +\item \code{fack_upper} is an upper bound based on the maximum diameter given by +theorem 3 of \insertCite{Fack2002;textual}{TreeDist}. +\item \code{li_upper} is an upper bound based on the sorting approach of +\insertCite{Li1996;textual}{TreeDist}. } \code{NNIDiameter()} returns a matrix specifying (bounds on) the diameter @@ -41,11 +51,9 @@ Columns correspond to: \item \code{exact}: The exact value of the diameter, where \emph{n} < 13; \item \code{liMax}: Upper bound on diameter following Li \emph{et al.} (1996), i.e. \deqn{n \log2{n} + \textrm{O}(n)}{n log₂ _n_ + O(_n_)}; -\item \code{fackMax}: Upper bound on diameter following Fack \emph{et al}. (2002), i.e. -(\deqn{N - 2}{_N_ − 2}) ceiling(\deqn{\log2{n}}{log₂ _N_}) -\itemize{ -\item \emph{N}; -} +\item \code{fackMax}: Upper bound on diameter following theorem 3 of +Fack \emph{et al}. (2002), i.e. \emph{N} \emph{c} − 2^\emph{c} + 1 + 2(\emph{N} − +ceiling(log₂ (\emph{N} + 1))), where \emph{c} = ceiling(log₂ \emph{N}); \item \code{max}: The smaller of \code{liMax} and \code{fackMax}; } @@ -83,6 +91,17 @@ Leaves: \tab 1 \tab 2 \tab 3 \tab 4 \tab 5 \tab 6 \tab 7 \tab 8 \tab 9 Diameter: \tab 0 \tab 0 \tab 0 \tab 1 \tab 3 \tab 5 \tab 7 \tab 10 \tab 12 \tab 15 \tab 18 \tab 21 \tab ? \cr } + +Upper bounds treat each region separately: transforming each region in turn +is a valid path between the trees, so the sum of the diameters of the regions +bounds the distance from above. +A diameter is the distance between the \emph{most} different pair of trees with +that many leaves, so it says nothing about how close a particular pair might +be: the only known lower bound is the number of unmatched edges +\insertCite{Li1996}{TreeDist}. +An optimal path need not leave matched edges untouched +\insertCite{Li1996}{TreeDist}, so even the exact distance between the trees +that each region induces is not a lower bound. } \examples{ library("TreeTools", quietly = TRUE) diff --git a/src/li_diameters.h b/src/li_diameters.h index a4c45853..b71bfda6 100644 --- a/src/li_diameters.h +++ b/src/li_diameters.h @@ -1,10 +1,12 @@ #ifndef _TREEDIST_LI_DIAMETERS_H #define _TREEDIST_LI_DIAMETERS_H -// Exact value of diameter for trees with 0..N_EXACT edges, -// calculated by Li et al. 1996. +// Exact diameter of NNI tree space on 0..N_EXACT leaves: the largest NNI +// distance between any two trees with that many leaves. +// Values for <= 11 leaves from Li et al. (1996, table 1); 12 leaves from +// Fack et al. (2002, table 2, n = 10 internal nodes). +// A diameter bounds the distance between a given pair from above, not below. constexpr int32_t N_EXACT = 12; constexpr int32_t exact_diameter[] = {0, 0, 0, 0, 1, 3, 5, 7, 10, 12, 15, 18, 21}; -constexpr int32_t min_diameter[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 38, 39, 40, 42, 43, 44, 46, 47, 48, 50, 51, 53, 54, 55, 57, 58, 60, 61, 63, 64, 66, 67, 69, 70, 72, 73, 74, 76, 78, 79, 81, 82, 84, 85, 87, 88, 90, 91, 93, 94, 96, 98, 99, 101, 102, 104, 106, 107, 109, 110, 112, 114, 115, 117, 119, 120, 122, 123, 125, 127, 128, 130, 132, 133, 135, 137, 138, 140, 142, 143, 145, 147, 148, 150, 152, 154, 155, 157, 159, 160, 162, 164, 166, 167, 169, 171, 173, 174, 176, 178, 180, 181, 183, 185, 187, 188, 190, 192, 194, 195, 197, 199, 201, 203, 204, 206, 208, 210, 212, 213, 215, 217, 219, 221, 222, 224, 226, 228, 230, 231, 233, 235, 237, 239, 241, 242, 244, 246, 248, 250, 252, 253, 255, 257, 259, 261, 263, 265, 267, 268, 270, 272, 274, 276, 278, 280, 281, 283, 285, 287, 289, 291, 293, 295, 297, 299, 300, 302, 304, 306, 308, 310, 312, 314, 316, 318, 320, 321, 323, 325, 327, 329, 331, 333, 335, 337, 339, 341, 343, 345, 347, 348, 350, 352, 354, 356, 358, 360, 362, 364, 366, 368, 370, 372, 374, 376, 378, 380, 382, 384, 386, 388, 390, 392, 394, 396, 398, 400, 402, 404, 406, 408, 410, 412, 414, 416, 418, 419, 421, 424, 426, 428, 430, 432, 434, 436, 438, 440, 442, 444, 446, 448, 450, 452, 454, 456, 458, 460, 462, 464, 466, 468, 470, 472, 474, 476, 478, 480, 482, 484, 486, 488, 490, 492, 494, 496, 498, 501, 503, 505, 507, 509, 511, 513, 515, 517, 519, 521, 523, 525, 527, 529, 531, 534, 536, 538, 540, 542, 544, 546, 548, 550, 552, 554, 556, 558, 561, 563, 565, 567, 569, 571, 573, 575, 577, 579, 581, 584, 586, 588, 590, 592, 594, 596, 598, 600, 602, 605, 607, 609, 611, 613, 615, 617, 619, 621, 624, 626, 628, 630, 632, 634, 636, 638, 641, 643, 645, 647, 649, 651, 653, 655, 658, 660, 662, 664, 666, 668, 670, 673, 675, 677, 679, 681, 683, 685, 688, 690, 692, 694, 696, 698, 700, 703, 705, 707, 709, 711, 713, 716, 718, 720, 722, 724, 726, 729, 731, 733, 735, 737, 739, 742, 744, 746, 748, 750, 752, 755, 757, 759, 761, 763, 765, 768, 770, 772, 774, 776, 779, 781, 783, 785, 787, 789, 792, 794, 796, 798, 800, 803, 805, 807, 809, 811, 814, 816, 818, 820, 822, 825, 827, 829, 831, 833, 836, 838, 840, 842, 844, 847, 849, 851, 853, 855, 858, 860, 862, 864, 867, 869, 871, 873, 875, 878, 880, 882, 884, 887, 889, 891, 893, 895, 898, 900, 902, 904, 907, 909, 911, 913, 915, 918, 920, 922, 924, 927, 929, 931, 933, 936, 938, 940, 942, 945, 947, 949, 951, 954, 956, 958, 960, 963, 965, 967, 969, 972, 974, 976, 978, 981, 983, 985, 987, 990, 992, 994, 996, 999, 1001, 1003, 1005, 1008, 1010, 1012, 1014, 1017, 1019, 1021, 1023, 1026, 1028, 1030, 1033, 1035, 1037, 1039, 1042, 1044, 1046, 1048, 1051, 1053, 1055, 1058, 1060, 1062, 1064, 1067, 1069, 1071, 1074, 1076, 1078, 1080, 1083, 1085, 1087, 1090, 1092, 1094, 1096, 1099, 1101, 1103, 1106, 1108, 1110, 1112, 1115, 1117, 1119, 1122, 1124, 1126, 1128, 1131, 1133, 1135, 1138, 1140, 1142, 1145, 1147, 1149, 1151, 1154, 1156, 1158, 1161, 1163, 1165, 1168, 1170, 1172, 1175, 1177, 1179, 1181, 1184, 1186, 1188, 1191, 1193, 1195, 1198, 1200, 1202, 1205, 1207, 1209, 1212, 1214, 1216, 1218, 1221, 1223, 1225, 1228, 1230, 1232, 1235, 1237, 1239, 1242, 1244, 1246, 1249, 1251, 1253, 1256, 1258, 1260, 1263, 1265, 1267, 1270, 1272, 1274, 1277, 1279, 1281, 1284, 1286, 1288, 1291, 1293, 1295, 1298, 1300, 1302, 1305, 1307, 1309, 1312, 1314, 1316, 1319, 1321, 1323, 1326, 1328, 1330, 1333, 1335, 1338, 1340, 1342, 1345, 1347, 1349, 1352, 1354, 1356, 1359, 1361, 1363, 1366, 1368, 1370, 1373, 1375, 1378, 1380, 1382, 1385, 1387, 1389, 1392, 1394, 1396, 1399, 1401, 1404, 1406, 1408, 1411, 1413, 1415, 1418, 1420, 1422, 1425, 1427, 1430, 1432, 1434, 1437, 1439, 1441, 1444, 1446, 1448, 1451, 1453, 1456, 1458, 1460, 1463, 1465, 1467, 1470, 1472, 1475, 1477, 1479, 1482, 1484, 1487, 1489, 1491, 1494, 1496, 1498, 1501, 1503, 1506, 1508, 1510, 1513, 1515, 1518, 1520, 1522, 1525, 1527, 1529, 1532, 1534, 1537, 1539, 1541, 1544, 1546, 1549, 1551, 1553, 1556, 1558, 1561, 1563, 1565, 1568, 1570, 1573, 1575, 1577, 1580, 1582, 1585, 1587, 1589, 1592, 1594, 1597, 1599, 1601, 1604, 1606, 1609, 1611, 1613, 1616, 1618, 1621, 1623, 1625, 1628, 1630, 1633, 1635, 1637, 1640, 1642, 1645, 1647, 1650, 1652, 1654, 1657, 1659, 1662, 1664, 1666, 1669, 1671, 1674, 1676, 1678, 1681, 1683, 1686, 1688, 1691, 1693, 1695, 1698, 1700, 1703, 1705, 1708, 1710, 1712, 1715, 1717, 1720, 1722, 1725, 1727, 1729, 1732, 1734, 1737, 1739, 1742, 1744, 1746, 1749, 1751, 1754, 1756, 1759, 1761, 1763, 1766, 1768, 1771, 1773, 1776, 1778, 1780, 1783, 1785, 1788, 1790, 1793, 1795, 1797, 1800, 1802, 1805, 1807, 1810, 1812, 1815, 1817, 1819, 1822, 1824, 1827, 1829, 1832, 1834, 1837, 1839, 1841, 1844, 1846, 1849, 1851, 1854, 1856, 1859, 1861, 1863, 1866, 1868, 1871, 1873, 1876, 1878, 1881, 1883, 1886, 1888, 1890, 1893, 1895, 1898, 1900, 1903, 1905, 1908, 1910, 1913, 1915, 1917, 1920, 1922, 1925, 1927, 1930, 1932, 1935, 1937, 1940, 1942, 1945, 1947, 1949, 1952, 1954, 1957, 1959, 1962, 1964, 1967, 1969, 1972, 1974, 1977, 1979, 1982, 1984, 1986, 1989, 1991, 1994, 1996, 1999, 2001, 2004, 2006, 2009, 2011, 2014, 2016, 2019, 2021, 2024, 2026, 2029, 2031, 2033, 2036, 2038, 2041, 2043, 2046, 2048, 2051, 2053, 2056, 2058, 2061, 2063, 2066, 2068, 2071, 2073, 2076, 2078, 2081, 2083, 2086, 2088, 2091, 2093, 2096, 2098, 2100, 2103, 2105, 2108, 2110, 2113, 2115, 2118, 2120, 2123, 2125, 2128, 2130, 2133, 2135, 2138, 2140, 2143, 2145, 2148, 2150, 2153, 2155, 2158, 2160, 2163, 2165, 2168, 2170, 2173, 2175, 2178, 2180, 2183, 2185, 2188, 2190, 2193, 2195, 2198, 2200, 2203, 2205, 2208, 2210, 2213, 2215, 2218, 2220, 2223, 2225, 2228, 2230, 2233, 2235, 2238, 2240, 2243, 2245, 2248, 2250, 2253, 2255, 2258, 2260, 2263, 2265, 2268, 2270, 2273, 2275, 2278, 2280, 2283, 2285, 2288, 2291, 2293, 2296, 2298, 2301, 2303, 2306, 2308, 2311, 2313, 2316, 2318, 2321, 2323, 2326, 2328, 2331, 2333, 2336, 2338, 2341, 2343, 2346, 2348, 2351, 2353, 2356, 2359, 2361, 2364, 2366, 2369, 2371, 2374, 2376, 2379, 2381, 2384, 2386, 2389, 2391, 2394, 2396, 2399, 2401, 2404, 2407, 2409, 2412, 2414, 2417, 2419, 2422, 2424, 2427, 2429, 2432, 2434, 2437, 2439, 2442, 2444, 2447, 2450, 2452, 2455, 2457, 2460, 2462, 2465, 2467, 2470, 2472, 2475, 2477, 2480, 2483, 2485, 2488, 2490, 2493, 2495, 2498, 2500, 2503, 2505, 2508, 2510, 2513, 2516, 2518, 2521, 2523, 2526, 2528, 2531, 2533, 2536, 2538, 2541, 2544, 2546, 2549, 2551, 2554, 2556, 2559, 2561, 2564, 2566, 2569, 2572, 2574, 2577, 2579, 2582, 2584, 2587, 2589, 2592, 2595, 2597, 2600, 2602, 2605, 2607, 2610, 2612, 2615, 2617, 2620, 2623, 2625, 2628, 2630, 2633, 2635, 2638, 2641, 2643, 2646, 2648, 2651, 2653, 2656, 2658, 2661, 2664, 2666, 2669, 2671, 2674, 2676, 2679, 2681, 2684, 2687, 2689, 2692, 2694, 2697, 2699, 2702, 2705, 2707, 2710, 2712, 2715, 2717, 2720, 2723, 2725, 2728, 2730, 2733, 2735, 2738, 2740, 2743, 2746, 2748, 2751, 2753, 2756, 2758, 2761, 2764, 2766, 2769, 2771, 2774, 2776, 2779, 2782, 2784, 2787, 2789, 2792, 2794, 2797, 2800, 2802, 2805, 2807, 2810, 2813, 2815, 2818, 2820, 2823, 2825, 2828, 2831, 2833, 2836, 2838, 2841, 2843, 2846, 2849, 2851, 2854, 2856, 2859, 2862, 2864, 2867, 2869, 2872, 2874, 2877, 2880, 2882, 2885, 2887, 2890, 2893, 2895, 2898, 2900, 2903, 2905, 2908, 2911, 2913, 2916, 2918, 2921, 2924, 2926, 2929, 2931, 2934, 2937, 2939, 2942, 2944, 2947, 2950, 2952, 2955, 2957, 2960, 2962, 2965, 2968, 2970, 2973, 2975, 2978, 2981, 2983, 2986, 2988, 2991, 2994, 2996, 2999, 3001, 3004, 3007, 3009, 3012, 3014, 3017, 3020, 3022, 3025, 3027, 3030, 3033, 3035, 3038, 3040, 3043, 3046, 3048, 3051, 3053, 3056, 3059, 3061, 3064, 3066, 3069, 3072, 3074, 3077, 3079, 3082, 3085, 3087, 3090, 3092, 3095, 3098, 3100, 3103, 3106, 3108, 3111, 3113, 3116, 3119, 3121, 3124, 3126, 3129, 3132, 3134, 3137, 3139, 3142, 3145, 3147, 3150, 3153, 3155, 3158, 3160, 3163, 3166, 3168, 3171, 3173, 3176, 3179, 3181, 3184, 3186, 3189, 3192, 3194, 3197, 3200, 3202, 3205, 3207, 3210, 3213, 3215, 3218, 3221, 3223, 3226, 3228, 3231, 3234, 3236, 3239, 3241, 3244, 3247, 3249, 3252, 3255, 3257, 3260, 3262, 3265, 3268, 3270, 3273, 3276, 3278, 3281, 3283, 3286, 3289, 3291, 3294, 3297, 3299, 3302, 3304, 3307, 3310, 3312, 3315, 3318, 3320, 3323, 3325, 3328, 3331, 3333, 3336, 3339, 3341, 3344, 3347, 3349, 3352, 3354, 3357, 3360, 3362, 3365, 3368, 3370, 3373, 3375, 3378, 3381, 3383, 3386, 3389, 3391, 3394, 3397, 3399, 3402, 3404, 3407, 3410, 3412, 3415, 3418, 3420, 3423, 3426, 3428, 3431, 3434, 3436, 3439, 3441, 3444, 3447, 3449, 3452, 3455, 3457, 3460, 3463, 3465, 3468, 3470, 3473, 3476, 3478, 3481, 3484, 3486, 3489, 3492, 3494, 3497, 3500, 3502, 3505, 3507, 3510, 3513, 3515, 3518, 3521, 3523, 3526, 3529, 3531, 3534, 3537, 3539, 3542, 3545, 3547, 3550, 3553, 3555, 3558, 3560, 3563, 3566, 3568, 3571, 3574, 3576, 3579, 3582, 3584, 3587, 3590, 3592, 3595, 3598, 3600, 3603, 3606, 3608, 3611, 3614, 3616, 3619, 3621, 3624, 3627, 3629, 3632, 3635, 3637, 3640, 3643, 3645, 3648, 3651, 3653, 3656, 3659, 3661, 3664, 3667, 3669, 3672, 3675, 3677, 3680, 3683, 3685, 3688, 3691, 3693, 3696, 3699, 3701, 3704, 3707, 3709, 3712, 3715, 3717, 3720, 3723, 3725, 3728, 3731, 3733, 3736, 3739, 3741, 3744, 3747, 3749, 3752, 3755, 3757, 3760, 3763, 3765, 3768, 3771, 3773, 3776, 3779, 3781, 3784, 3787, 3789, 3792, 3795, 3797, 3800, 3803, 3805, 3808, 3811, 3813, 3816, 3819, 3821, 3824, 3827, 3829, 3832, 3835, 3837, 3840, 3843, 3845, 3848, 3851, 3853, 3856, 3859, 3861, 3864, 3867, 3869, 3872, 3875, 3877, 3880, 3883, 3886, 3888, 3891, 3894, 3896, 3899, 3902, 3904, 3907, 3910, 3912, 3915, 3918, 3920, 3923, 3926, 3928, 3931, 3934, 3936, 3939, 3942, 3944, 3947, 3950, 3953, 3955, 3958, 3961, 3963, 3966, 3969, 3971, 3974, 3977, 3979, 3982, 3985, 3987, 3990, 3993, 3995, 3998, 4001, 4004, 4006, 4009, 4012, 4014, 4017, 4020, 4022, 4025, 4028, 4030, 4033, 4036, 4038, 4041, 4044, 4047, 4049, 4052, 4055, 4057, 4060, 4063, 4065, 4068, 4071, 4073, 4076, 4079, 4082, 4084, 4087, 4090, 4092, 4095, 4098, 4100, 4103, 4106, 4108, 4111, 4114, 4117, 4119, 4122, 4125, 4127, 4130, 4133, 4135, 4138, 4141, 4144, 4146, 4149, 4152, 4154, 4157, 4160, 4162, 4165, 4168, 4171, 4173, 4176, 4179, 4181, 4184, 4187, 4189, 4192, 4195, 4198, 4200, 4203, 4206, 4208, 4211, 4214, 4216, 4219, 4222, 4225, 4227, 4230, 4233, 4235, 4238, 4241, 4244, 4246, 4249, 4252, 4254, 4257, 4260, 4262, 4265, 4268, 4271, 4273, 4276, 4279, 4281, 4284, 4287, 4290, 4292, 4295, 4298, 4300, 4303, 4306, 4308, 4311, 4314, 4317, 4319, 4322, 4325, 4327, 4330, 4333, 4336, 4338, 4341, 4344, 4346, 4349, 4352, 4355, 4357, 4360, 4363, 4365, 4368, 4371, 4374, 4376, 4379, 4382, 4384, 4387, 4390, 4393, 4395, 4398, 4401, 4403, 4406, 4409, 4412, 4414, 4417, 4420, 4422, 4425, 4428, 4431, 4433, 4436, 4439, 4442, 4444, 4447, 4450, 4452, 4455, 4458, 4461, 4463, 4466, 4469, 4471, 4474, 4477, 4480, 4482, 4485, 4488, 4491, 4493, 4496, 4499, 4501, 4504, 4507, 4510, 4512, 4515, 4518, 4520, 4523, 4526, 4529, 4531, 4534, 4537, 4540, 4542, 4545, 4548, 4550, 4553, 4556, 4559, 4561, 4564, 4567, 4570, 4572, 4575, 4578, 4580, 4583, 4586, 4589, 4591, 4594, 4597, 4600, 4602, 4605, 4608, 4611, 4613, 4616, 4619, 4621, 4624, 4627, 4630, 4632, 4635, 4638, 4641, 4643, 4646, 4649, 4651, 4654, 4657, 4660, 4662, 4665, 4668, 4671, 4673, 4676, 4679, 4682, 4684, 4687, 4690, 4693, 4695, 4698, 4701, 4703, 4706, 4709, 4712, 4714, 4717, 4720, 4723, 4725, 4728, 4731, 4734, 4736, 4739, 4742, 4745, 4747, 4750, 4753, 4756, 4758, 4761, 4764, 4766, 4769, 4772, 4775, 4777, 4780, 4783, 4786, 4788, 4791, 4794, 4797, 4799, 4802, 4805, 4808, 4810, 4813, 4816, 4819, 4821, 4824, 4827, 4830, 4832, 4835, 4838, 4841, 4843, 4846, 4849, 4852, 4854, 4857, 4860, 4863, 4865, 4868, 4871, 4874, 4876, 4879, 4882, 4885, 4887, 4890, 4893, 4896, 4898, 4901, 4904, 4907, 4909, 4912, 4915, 4918, 4920, 4923, 4926, 4929, 4931, 4934, 4937, 4940, 4942, 4945, 4948, 4951, 4953, 4956, 4959, 4962, 4964, 4967, 4970, 4973, 4975, 4978, 4981, 4984, 4986, 4989, 4992, 4995, 4997, 5000, 5003, 5006, 5008, 5011, 5014, 5017, 5019, 5022, 5025, 5028, 5031, 5033, 5036, 5039, 5042, 5044, 5047, 5050, 5053, 5055, 5058, 5061, 5064, 5066, 5069, 5072, 5075, 5077, 5080, 5083, 5086, 5088, 5091, 5094, 5097, 5100, 5102, 5105, 5108, 5111, 5113, 5116, 5119, 5122, 5124, 5127, 5130, 5133, 5135, 5138, 5141, 5144, 5147, 5149, 5152, 5155, 5158, 5160, 5163, 5166, 5169, 5171, 5174, 5177, 5180, 5182, 5185, 5188, 5191, 5194, 5196, 5199, 5202, 5205, 5207, 5210, 5213, 5216, 5218, 5221, 5224, 5227, 5230, 5232, 5235, 5238, 5241, 5243, 5246, 5249, 5252, 5255, 5257, 5260, 5263, 5266, 5268, 5271, 5274, 5277, 5279, 5282, 5285, 5288, 5291, 5293, 5296, 5299, 5302, 5304, 5307, 5310, 5313, 5316, 5318, 5321, 5324, 5327, 5329, 5332, 5335, 5338, 5341, 5343, 5346, 5349, 5352, 5354, 5357, 5360, 5363, 5366, 5368, 5371, 5374, 5377, 5379, 5382, 5385, 5388, 5391, 5393, 5396, 5399, 5402, 5404, 5407, 5410, 5413, 5416, 5418, 5421, 5424, 5427, 5429, 5432, 5435, 5438, 5441, 5443, 5446, 5449, 5452, 5455, 5457, 5460, 5463, 5466, 5468, 5471, 5474, 5477, 5480, 5482, 5485, 5488, 5491, 5494, 5496, 5499, 5502, 5505, 5507, 5510, 5513, 5516, 5519, 5521, 5524, 5527, 5530, 5533, 5535, 5538, 5541, 5544, 5547, 5549, 5552, 5555, 5558, 5560, 5563, 5566, 5569, 5572, 5574, 5577, 5580, 5583, 5586, 5588, 5591, 5594, 5597, 5600, 5602, 5605, 5608, 5611, 5613, 5616, 5619, 5622, 5625, 5627, 5630, 5633, 5636, 5639, 5641, 5644, 5647, 5650, 5653, 5655, 5658, 5661, 5664, 5667, 5669, 5672, 5675, 5678, 5681, 5683, 5686, 5689, 5692, 5695, 5697, 5700, 5703, 5706, 5709, 5711, 5714, 5717, 5720, 5723, 5725, 5728, 5731, 5734, 5737, 5739, 5742, 5745, 5748, 5751, 5753, 5756, 5759, 5762, 5765, 5767, 5770, 5773, 5776, 5779, 5781, 5784, 5787, 5790, 5793, 5795, 5798, 5801, 5804, 5807, 5809, 5812, 5815, 5818, 5821, 5823, 5826, 5829, 5832, 5835, 5837, 5840, 5843, 5846, 5849, 5851, 5854, 5857, 5860, 5863, 5865, 5868, 5871, 5874, 5877, 5880, 5882, 5885, 5888, 5891, 5894, 5896, 5899, 5902, 5905, 5908, 5910, 5913, 5916, 5919, 5922, 5924, 5927, 5930, 5933, 5936, 5939, 5941, 5944, 5947, 5950, 5953, 5955, 5958, 5961, 5964, 5967, 5969, 5972, 5975, 5978, 5981, 5984, 5986, 5989, 5992, 5995, 5998, 6000, 6003, 6006, 6009, 6012, 6014, 6017, 6020, 6023, 6026, 6029, 6031, 6034, 6037, 6040, 6043, 6045, 6048, 6051, 6054, 6057, 6060, 6062, 6065, 6068, 6071, 6074, 6076, 6079, 6082, 6085, 6088, 6091, 6093, 6096, 6099, 6102, 6105, 6107, 6110, 6113, 6116, 6119, 6122, 6124, 6127, 6130, 6133, 6136, 6138, 6141, 6144, 6147, 6150, 6153, 6155, 6158, 6161, 6164, 6167, 6169, 6172, 6175, 6178, 6181, 6184, 6186, 6189, 6192, 6195, 6198, 6201, 6203, 6206, 6209, 6212, 6215, 6217, 6220, 6223, 6226, 6229, 6232, 6234, 6237, 6240, 6243, 6246, 6249, 6251, 6254, 6257, 6260, 6263, 6266, 6268, 6271, 6274, 6277, 6280, 6283, 6285, 6288, 6291, 6294, 6297, 6299, 6302, 6305, 6308, 6311, 6314, 6316, 6319, 6322, 6325, 6328, 6331, 6333, 6336, 6339, 6342, 6345, 6348, 6350, 6353, 6356, 6359, 6362, 6365, 6367, 6370, 6373, 6376, 6379, 6382, 6384, 6387, 6390, 6393, 6396, 6399, 6401, 6404, 6407, 6410, 6413, 6416, 6418, 6421, 6424, 6427, 6430, 6433, 6435, 6438, 6441, 6444, 6447, 6450, 6452, 6455, 6458, 6461, 6464, 6467, 6469, 6472, 6475, 6478, 6481, 6484, 6487, 6489, 6492, 6495, 6498, 6501, 6504, 6506, 6509, 6512, 6515, 6518, 6521, 6523, 6526, 6529, 6532, 6535, 6538, 6540, 6543, 6546, 6549, 6552, 6555, 6558, 6560, 6563, 6566, 6569, 6572, 6575, 6577, 6580, 6583, 6586, 6589, 6592, 6594, 6597, 6600, 6603, 6606, 6609, 6612, 6614, 6617, 6620, 6623, 6626, 6629, 6631, 6634, 6637, 6640, 6643, 6646, 6649, 6651, 6654, 6657, 6660, 6663, 6666, 6668, 6671, 6674, 6677, 6680, 6683, 6686, 6688, 6691, 6694, 6697, 6700, 6703, 6705, 6708, 6711, 6714, 6717, 6720, 6723, 6725, 6728, 6731, 6734, 6737, 6740, 6743, 6745, 6748, 6751, 6754, 6757, 6760, 6762, 6765, 6768, 6771, 6774, 6777, 6780, 6782, 6785, 6788, 6791, 6794, 6797, 6800, 6802, 6805, 6808, 6811, 6814, 6817, 6820, 6822, 6825, 6828, 6831, 6834, 6837, 6840, 6842, 6845, 6848, 6851, 6854, 6857, 6860, 6862, 6865, 6868, 6871, 6874, 6877, 6880, 6882, 6885, 6888, 6891, 6894, 6897, 6900, 6902, 6905, 6908, 6911, 6914, 6917, 6920, 6922, 6925, 6928, 6931, 6934, 6937, 6940, 6942, 6945, 6948, 6951, 6954, 6957, 6960, 6962, 6965, 6968, 6971, 6974, 6977, 6980, 6982, 6985, 6988, 6991, 6994, 6997, 7000, 7002, 7005, 7008, 7011, 7014, 7017, 7020, 7023, 7025, 7028, 7031, 7034, 7037, 7040, 7043, 7045, 7048, 7051, 7054, 7057, 7060, 7063, 7066, 7068, 7071, 7074, 7077, 7080, 7083, 7086, 7088, 7091, 7094, 7097, 7100, 7103, 7106, 7109, 7111, 7114, 7117, 7120, 7123, 7126, 7129, 7131, 7134, 7137, 7140, 7143, 7146, 7149, 7152, 7154, 7157, 7160, 7163, 7166, 7169, 7172, 7174, 7177, 7180, 7183, 7186, 7189, 7192, 7195, 7197, 7200, 7203, 7206, 7209, 7212, 7215, 7218, 7220, 7223, 7226, 7229, 7232, 7235, 7238, 7241, 7243, 7246, 7249, 7252, 7255, 7258, 7261, 7264, 7266, 7269, 7272, 7275, 7278, 7281, 7284, 7287, 7289, 7292, 7295, 7298, 7301, 7304, 7307, 7310, 7312, 7315, 7318, 7321, 7324, 7327, 7330, 7333, 7335, 7338, 7341, 7344, 7347, 7350, 7353, 7356, 7358, 7361, 7364, 7367, 7370, 7373, 7376, 7379, 7381, 7384, 7387, 7390, 7393, 7396, 7399, 7402, 7405, 7407, 7410, 7413, 7416, 7419, 7422, 7425, 7428, 7430, 7433, 7436, 7439, 7442, 7445, 7448, 7451, 7453, 7456, 7459, 7462, 7465, 7468, 7471, 7474, 7477, 7479, 7482, 7485, 7488, 7491, 7494, 7497, 7500, 7503, 7505, 7508, 7511, 7514, 7517, 7520, 7523, 7526, 7528, 7531, 7534, 7537, 7540, 7543, 7546, 7549, 7552, 7554, 7557, 7560, 7563, 7566, 7569, 7572, 7575, 7578, 7580, 7583, 7586, 7589, 7592, 7595, 7598, 7601, 7604, 7606, 7609, 7612, 7615, 7618, 7621, 7624, 7627, 7630, 7632, 7635, 7638, 7641, 7644, 7647, 7650, 7653, 7656, 7658, 7661, 7664, 7667, 7670, 7673, 7676, 7679, 7682, 7684, 7687, 7690, 7693, 7696, 7699, 7702, 7705, 7708, 7710, 7713, 7716, 7719, 7722, 7725, 7728, 7731, 7734, 7737, 7739, 7742, 7745, 7748, 7751, 7754, 7757, 7760, 7763, 7765, 7768, 7771, 7774, 7777, 7780, 7783, 7786, 7789, 7792, 7794, 7797, 7800, 7803, 7806, 7809, 7812, 7815, 7818, 7821, 7823, 7826, 7829, 7832, 7835, 7838, 7841, 7844, 7847, 7849, 7852, 7855, 7858, 7861, 7864, 7867, 7870, 7873, 7876, 7878, 7881, 7884, 7887, 7890, 7893, 7896, 7899, 7902, 7905, 7907, 7910, 7913, 7916, 7919, 7922, 7925, 7928, 7931, 7934, 7937, 7939, 7942, 7945, 7948, 7951, 7954, 7957, 7960, 7963, 7966, 7968, 7971, 7974, 7977, 7980, 7983, 7986, 7989, 7992, 7995, 7997, 8000, 8003, 8006, 8009, 8012, 8015, 8018, 8021, 8024, 8027, 8029, 8032, 8035, 8038, 8041, 8044, 8047, 8050, 8053, 8056, 8058, 8061, 8064, 8067, 8070, 8073, 8076, 8079, 8082, 8085, 8088, 8090, 8093, 8096, 8099, 8102, 8105, 8108, 8111, 8114, 8117, 8120, 8122, 8125, 8128, 8131, 8134, 8137, 8140, 8143, 8146, 8149, 8152, 8154, 8157, 8160, 8163, 8166, 8169, 8172, 8175, 8178, 8181, 8184, 8187, 8189, 8192, 8195, 8198, 8201, 8204, 8207, 8210, 8213, 8216, 8219, 8221, 8224, 8227, 8230, 8233, 8236, 8239, 8242, 8245, 8248, 8251, 8254, 8256, 8259, 8262, 8265, 8268, 8271, 8274, 8277, 8280, 8283, 8286, 8288, 8291, 8294, 8297, 8300, 8303, 8306, 8309, 8312, 8315, 8318, 8321, 8323, 8326, 8329, 8332, 8335, 8338, 8341, 8344, 8347, 8350, 8353, 8356, 8359, 8361, 8364, 8367, 8370, 8373, 8376, 8379, 8382, 8385, 8388, 8391, 8394, 8396, 8399, 8402, 8405, 8408, 8411, 8414, 8417, 8420, 8423, 8426, 8429, 8432, 8434, 8437, 8440, 8443, 8446, 8449, 8452, 8455, 8458, 8461, 8464, 8467, 8469, 8472, 8475, 8478, 8481, 8484, 8487, 8490, 8493, 8496, 8499, 8502, 8505, 8507, 8510, 8513, 8516, 8519, 8522, 8525, 8528, 8531, 8534, 8537, 8540, 8543, 8546, 8548, 8551, 8554, 8557, 8560, 8563, 8566, 8569, 8572, 8575, 8578, 8581, 8584, 8586, 8589, 8592, 8595, 8598, 8601, 8604, 8607, 8610, 8613, 8616, 8619, 8622, 8625, 8627, 8630, 8633, 8636, 8639, 8642, 8645, 8648, 8651, 8654, 8657, 8660, 8663, 8666, 8668, 8671, 8674, 8677, 8680, 8683, 8686, 8689, 8692, 8695, 8698, 8701, 8704, 8707, 8710, 8712, 8715, 8718, 8721, 8724, 8727, 8730, 8733, 8736, 8739, 8742, 8745, 8748, 8751, 8753, 8756, 8759, 8762, 8765, 8768, 8771, 8774, 8777, 8780, 8783, 8786, 8789, 8792, 8795, 8798, 8800, 8803, 8806, 8809, 8812, 8815, 8818, 8821, 8824, 8827, 8830, 8833, 8836, 8839, 8842, 8844, 8847, 8850, 8853, 8856, 8859, 8862, 8865, 8868, 8871, 8874, 8877, 8880, 8883, 8886, 8889, 8891, 8894, 8897, 8900, 8903, 8906, 8909, 8912, 8915, 8918, 8921, 8924, 8927, 8930, 8933, 8936, 8938, 8941, 8944, 8947, 8950, 8953, 8956, 8959, 8962, 8965, 8968, 8971, 8974, 8977, 8980, 8983, 8986, 8988, 8991, 8994, 8997, 9000, 9003, 9006, 9009, 9012, 9015, 9018, 9021, 9024, 9027, 9030, 9033, 9036, 9038, 9041, 9044, 9047, 9050, 9053, 9056, 9059, 9062, 9065, 9068, 9071, 9074, 9077, 9080, 9083, 9086, 9089, 9091, 9094, 9097, 9100, 9103, 9106, 9109, 9112, 9115, 9118, 9121, 9124, 9127, 9130, 9133, 9136, 9139, 9142, 9145, 9147, 9150, 9153, 9156, 9159, 9162, 9165, 9168, 9171, 9174, 9177, 9180, 9183, 9186, 9189, 9192, 9195, 9198, 9201, 9203, 9206, 9209, 9212, 9215, 9218, 9221, 9224, 9227, 9230, 9233, 9236, 9239, 9242, 9245, 9248, 9251, 9254, 9257, 9260, 9262, 9265, 9268, 9271, 9274, 9277, 9280, 9283, 9286, 9289, 9292, 9295, 9298, 9301, 9304, 9307, 9310, 9313, 9316, 9319, 9322, 9324, 9327, 9330, 9333, 9336, 9339, 9342, 9345, 9348, 9351, 9354, 9357, 9360, 9363, 9366, 9369, 9372, 9375, 9378, 9381, 9384, 9387, 9389, 9392, 9395, 9398, 9401, 9404, 9407, 9410, 9413, 9416, 9419, 9422, 9425, 9428, 9431, 9434, 9437, 9440, 9443, 9446, 9449, 9452, 9455, 9457, 9460, 9463, 9466, 9469, 9472, 9475, 9478, 9481, 9484, 9487, 9490, 9493, 9496, 9499, 9502, 9505, 9508, 9511, 9514, 9517, 9520, 9523, 9526, 9528, 9531, 9534, 9537, 9540, 9543, 9546, 9549, 9552, 9555, 9558, 9561, 9564, 9567, 9570, 9573, 9576, 9579, 9582, 9585, 9588, 9591, 9594, 9597, 9600, 9603, 9605, 9608, 9611, 9614, 9617, 9620, 9623, 9626, 9629, 9632, 9635, 9638, 9641, 9644, 9647, 9650, 9653, 9656, 9659, 9662, 9665, 9668, 9671, 9674, 9677, 9680, 9683, 9685, 9688, 9691, 9694, 9697, 9700, 9703, 9706, 9709, 9712, 9715, 9718, 9721, 9724, 9727, 9730, 9733, 9736, 9739, 9742, 9745, 9748, 9751, 9754, 9757, 9760, 9763, 9766, 9769, 9772, 9774, 9777, 9780, 9783, 9786, 9789, 9792, 9795, 9798, 9801, 9804, 9807, 9810, 9813, 9816, 9819, 9822, 9825, 9828, 9831, 9834, 9837, 9840, 9843, 9846, 9849, 9852, 9855, 9858, 9861, 9864, 9867, 9870, 9872, 9875, 9878, 9881, 9884, 9887, 9890, 9893, 9896, 9899, 9902, 9905, 9908, 9911, 9914, 9917, 9920, 9923, 9926, 9929, 9932, 9935, 9938, 9941, 9944, 9947, 9950, 9953, 9956, 9959, 9962, 9965, 9968, 9971, 9974, 9977, 9980, 9982, 9985, 9988, 9991, 9994, 9997, 10000, 10003, 10006, 10009, 10012, 10015, 10018, 10021, 10024, 10027, 10030, 10033, 10036, 10039, 10042, 10045, 10048, 10051, 10054, 10057, 10060, 10063, 10066, 10069, 10072, 10075, 10078, 10081, 10084, 10087, 10090, 10093, 10096, 10099, 10102, 10105, 10108, 10110, 10113, 10116, 10119, 10122, 10125, 10128, 10131, 10134, 10137, 10140, 10143, 10146, 10149, 10152, 10155, 10158, 10161, 10164, 10167, 10170, 10173, 10176, 10179, 10182, 10185, 10188, 10191, 10194, 10197, 10200, 10203, 10206, 10209, 10212, 10215, 10218, 10221, 10224, 10227, 10230, 10233, 10236, 10239, 10242, 10245, 10248, 10251, 10254, 10257, 10260, 10263, 10266, 10268, 10271, 10274, 10277, 10280, 10283, 10286, 10289, 10292, 10295, 10298, 10301, 10304, 10307, 10310, 10313, 10316, 10319, 10322, 10325, 10328, 10331, 10334, 10337, 10340, 10343, 10346, 10349, 10352, 10355, 10358, 10361, 10364, 10367, 10370, 10373, 10376, 10379, 10382, 10385, 10388, 10391, 10394, 10397, 10400, 10403, 10406, 10409, 10412, 10415, 10418, 10421, 10424, 10427, 10430, 10433, 10436, 10439, 10442, 10445, 10448, 10451, 10454, 10457, 10460, 10463, 10466, 10469, 10472, 10475, 10478, 10481, 10484, 10487, 10490, 10493, 10496, 10498, 10501, 10504, 10507, 10510, 10513, 10516, 10519, 10522, 10525, 10528, 10531, 10534, 10537, 10540, 10543, 10546, 10549, 10552, 10555, 10558, 10561, 10564, 10567, 10570, 10573, 10576, 10579, 10582, 10585, 10588, 10591, 10594, 10597, 10600, 10603, 10606, 10609, 10612, 10615, 10618, 10621, 10624, 10627, 10630, 10633, 10636, 10639, 10642, 10645, 10648, 10651, 10654, 10657, 10660, 10663, 10666, 10669, 10672, 10675, 10678, 10681, 10684, 10687, 10690, 10693, 10696, 10699, 10702, 10705, 10708, 10711, 10714, 10717, 10720, 10723, 10726, 10729, 10732, 10735, 10738, 10741, 10744, 10747, 10750, 10753, 10756, 10759, 10762, 10765, 10768, 10771, 10774, 10777, 10780, 10783, 10786, 10789, 10792, 10795, 10798, 10801, 10804, 10807, 10810, 10813, 10816, 10819, 10822, 10825, 10828, 10831, 10834, 10837, 10840, 10843, 10846, 10849, 10852, 10855, 10858, 10861, 10864, 10867, 10870, 10873, 10876, 10879, 10882, 10885, 10888, 10891, 10894, 10897, 10900, 10903, 10906, 10909, 10912, 10915, 10918, 10921, 10924, 10927, 10930, 10933, 10936, 10939, 10942, 10945, 10948, 10951, 10954, 10957, 10960, 10963, 10966, 10969, 10972, 10975, 10978, 10981, 10984, 10987, 10990, 10993, 10996, 10999, 11002, 11005, 11008, 11011, 11014, 11017, 11020, 11023, 11026, 11029, 11032, 11035, 11038, 11041, 11044, 11047, 11050, 11053, 11056, 11059, 11062, 11065, 11068, 11071, 11074, 11077, 11080, 11083, 11086, 11089, 11092, 11095, 11098, 11101, 11104, 11107, 11110, 11113, 11116, 11119, 11122, 11125, 11128, 11132, 11135, 11138, 11141, 11144, 11147, 11150, 11153, 11156, 11159, 11162, 11165, 11168, 11171, 11174, 11177, 11180, 11183, 11186, 11189, 11192, 11195, 11198, 11201, 11204, 11207, 11210, 11213, 11216, 11219, 11222, 11225, 11228, 11231, 11234, 11237, 11240, 11243, 11246, 11249, 11252, 11255, 11258, 11261, 11264, 11267, 11270, 11273, 11276, 11279, 11282, 11285, 11288, 11291, 11294, 11297, 11300, 11303, 11306, 11309, 11312, 11315, 11318, 11321, 11324, 11327, 11330, 11333, 11336, 11339, 11342, 11345, 11348, 11351, 11354, 11357, 11360, 11363, 11366, 11370, 11373, 11376, 11379, 11382, 11385, 11388, 11391, 11394, 11397, 11400, 11403, 11406, 11409, 11412, 11415, 11418, 11421, 11424, 11427, 11430, 11433, 11436, 11439, 11442, 11445, 11448, 11451, 11454, 11457, 11460, 11463, 11466, 11469, 11472, 11475, 11478, 11481, 11484, 11487, 11490, 11493, 11496, 11499, 11502, 11505, 11508, 11511, 11514, 11517, 11520, 11523, 11526, 11529, 11532, 11536, 11539, 11542, 11545, 11548, 11551, 11554, 11557, 11560, 11563, 11566, 11569, 11572, 11575, 11578, 11581, 11584, 11587, 11590, 11593, 11596, 11599, 11602, 11605, 11608, 11611, 11614, 11617, 11620, 11623, 11626, 11629, 11632, 11635, 11638, 11641, 11644, 11647, 11650, 11653, 11656, 11659, 11662, 11665, 11668, 11672, 11675, 11678, 11681, 11684, 11687, 11690, 11693, 11696, 11699, 11702, 11705, 11708, 11711, 11714, 11717, 11720, 11723, 11726, 11729, 11732, 11735, 11738, 11741, 11744, 11747, 11750, 11753, 11756, 11759, 11762, 11765, 11768, 11771, 11774, 11777, 11780, 11783, 11787, 11790, 11793, 11796, 11799, 11802, 11805, 11808, 11811, 11814, 11817, 11820, 11823, 11826, 11829, 11832, 11835, 11838, 11841, 11844, 11847, 11850, 11853, 11856, 11859, 11862, 11865, 11868, 11871, 11874, 11877, 11880, 11883, 11886, 11889, 11893, 11896, 11899, 11902, 11905, 11908, 11911, 11914, 11917, 11920, 11923, 11926, 11929, 11932, 11935, 11938, 11941, 11944, 11947, 11950, 11953, 11956, 11959, 11962, 11965, 11968, 11971, 11974, 11977, 11980, 11983, 11986, 11990, 11993, 11996, 11999, 12002, 12005, 12008, 12011, 12014, 12017, 12020, 12023, 12026, 12029, 12032, 12035, 12038, 12041, 12044, 12047, 12050, 12053, 12056, 12059, 12062, 12065, 12068, 12071, 12074, 12078, 12081, 12084, 12087, 12090, 12093, 12096, 12099, 12102, 12105, 12108, 12111, 12114, 12117, 12120, 12123, 12126, 12129, 12132, 12135, 12138, 12141, 12144, 12147, 12150, 12153, 12156, 12159, 12163, 12166, 12169, 12172, 12175, 12178, 12181, 12184, 12187, 12190, 12193, 12196, 12199, 12202, 12205, 12208, 12211, 12214, 12217, 12220, 12223, 12226, 12229, 12232, 12235, 12238, 12242, 12245, 12248, 12251, 12254, 12257, 12260, 12263, 12266, 12269, 12272, 12275, 12278, 12281, 12284, 12287, 12290, 12293, 12296, 12299, 12302, 12305, 12308, 12311, 12314, 12318, 12321, 12324, 12327, 12330, 12333, 12336, 12339, 12342, 12345, 12348, 12351, 12354, 12357, 12360, 12363, 12366, 12369, 12372, 12375, 12378, 12381, 12384, 12388, 12391, 12394, 12397, 12400, 12403, 12406, 12409, 12412, 12415, 12418, 12421, 12424, 12427, 12430, 12433, 12436, 12439, 12442, 12445, 12448, 12451, 12454, 12458, 12461, 12464, 12467, 12470, 12473, 12476, 12479, 12482, 12485, 12488, 12491, 12494, 12497, 12500, 12503, 12506, 12509, 12512, 12515, 12518, 12521, 12525, 12528, 12531, 12534, 12537, 12540, 12543, 12546, 12549, 12552, 12555, 12558, 12561, 12564, 12567, 12570, 12573, 12576, 12579, 12582, 12585, 12589, 12592, 12595, 12598, 12601, 12604, 12607, 12610, 12613, 12616, 12619, 12622, 12625, 12628, 12631, 12634, 12637, 12640, 12643, 12646, 12650, 12653, 12656, 12659, 12662, 12665, 12668, 12671, 12674, 12677, 12680, 12683, 12686, 12689, 12692, 12695, 12698, 12701, 12704, 12707, 12711, 12714, 12717, 12720, 12723, 12726, 12729, 12732, 12735, 12738, 12741, 12744, 12747, 12750, 12753, 12756, 12759, 12762, 12765, 12769, 12772, 12775, 12778, 12781, 12784, 12787, 12790, 12793, 12796, 12799, 12802, 12805, 12808, 12811, 12814, 12817, 12820, 12823, 12827, 12830, 12833, 12836, 12839, 12842, 12845, 12848, 12851, 12854, 12857, 12860, 12863, 12866, 12869, 12872, 12875, 12878, 12882, 12885, 12888, 12891, 12894, 12897, 12900, 12903, 12906, 12909, 12912, 12915, 12918, 12921, 12924, 12927, 12930, 12933, 12937, 12940, 12943, 12946, 12949, 12952, 12955, 12958, 12961, 12964, 12967, 12970, 12973, 12976, 12979, 12982, 12985, 12989, 12992, 12995, 12998, 13001, 13004, 13007, 13010, 13013, 13016, 13019, 13022, 13025, 13028, 13031, 13034, 13037, 13041, 13044, 13047, 13050, 13053, 13056, 13059, 13062, 13065, 13068, 13071, 13074, 13077, 13080, 13083, 13086, 13089, 13093, 13096, 13099, 13102, 13105, 13108, 13111, 13114, 13117, 13120, 13123, 13126, 13129, 13132, 13135, 13138, 13142, 13145, 13148, 13151, 13154, 13157, 13160, 13163, 13166, 13169, 13172, 13175, 13178, 13181, 13184, 13187, 13191, 13194, 13197, 13200, 13203, 13206, 13209, 13212, 13215, 13218, 13221, 13224, 13227, 13230, 13233, 13236, 13240, 13243, 13246, 13249, 13252, 13255, 13258, 13261, 13264, 13267, 13270, 13273, 13276, 13279, 13282, 13286, 13289, 13292, 13295, 13298, 13301, 13304, 13307, 13310, 13313, 13316, 13319, 13322, 13325, 13328, 13332, 13335, 13338, 13341, 13344, 13347, 13350, 13353, 13356, 13359, 13362, 13365, 13368, 13371, 13374, 13378, 13381, 13384, 13387, 13390, 13393, 13396, 13399, 13402, 13405, 13408, 13411, 13414, 13417, 13420, 13424, 13427, 13430, 13433, 13436, 13439, 13442, 13445, 13448, 13451, 13454, 13457, 13460, 13463, 13467, 13470, 13473, 13476, 13479, 13482, 13485, 13488, 13491, 13494, 13497, 13500, 13503, 13506, 13509, 13513, 13516, 13519, 13522, 13525, 13528, 13531, 13534, 13537, 13540, 13543, 13546, 13549, 13552}; - + #endif diff --git a/src/nni_distance.cpp b/src/nni_distance.cpp index db3a5b3c..7b53ee84 100644 --- a/src/nni_distance.cpp +++ b/src/nni_distance.cpp @@ -99,8 +99,14 @@ static void initialize_cache() { for (int32_t i = 8192 + 1; i != 16384 + 1; i++) lg2_ceiling[i] = 14; for (int32_t i = 16384 + 1; i != 32768 + 1; i++) lg2_ceiling[i] = 15; + /* Fack et al. (2002), theorem 3: with N = n - 2 internal nodes, + * diam <= N ceil(lg N) - 2^ceil(lg N) + 1 + 2 (N - ceil(lg(N + 1))). + * Never looser than the bound in their theorem 1. */ for (int32_t i = 4; i != NNI_MAX_TIPS + 1; i++) { - fack_lookup[i] = int32_t(((i - 2 - 2) * lg2_ceiling[i - 2]) + i - 2); + const int32_t n_node = i - 2; + const int32_t lg_n = lg2_ceiling[n_node]; + fack_lookup[i] = n_node * lg_n - (int32_t(1) << lg_n) + 1 + + 2 * (n_node - lg2_ceiling[n_node + 1]); } for (int32_t i = 4; i != NNI_MAX_TIPS + 1; i++) { const int32_t log_ceiling = lg2_ceiling[i]; @@ -138,14 +144,16 @@ inline void ensure_cache() { // Score subtree, add to score, and reset subtree size void update_score(const int32_t subtree_edges, int32_t *lower_bound, - int32_t *closest_lower_bound, int32_t *tight_bound, int32_t *closest_upper_bound, int32_t *loose_bound, int32_t *li_bound, int32_t *fack_bound) { if (subtree_edges) { const int32_t subtree_tips = subtree_edges + 3; - const int32_t lower = min_diameter[subtree_edges]; + // Each unmatched edge must undergo at least one NNI (Li et al. 1996). + // Diameters bound a single pair's distance from above only, so they + // cannot tighten this lower bound. + const int32_t lower = subtree_edges; const int32_t this_li = li[subtree_tips], fack = fack_lookup[subtree_tips], upper = this_li < fack ? this_li : fack; @@ -157,13 +165,11 @@ void update_score(const int32_t subtree_edges, int32_t *lower_bound, if (subtree_tips <= N_EXACT) { const int32_t exact = exact_diameter[subtree_tips]; - *closest_lower_bound += exact; if (*tight_bound != NA_INT32) { *tight_bound += exact; } *closest_upper_bound += exact; } else { - *closest_lower_bound += lower; *tight_bound = NA_INT32; *closest_upper_bound += upper; } @@ -296,7 +302,6 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, const int32_t n_edge = int32_t(edge1.nrow()); int32_t lower_bound = 0; - int32_t best_lower_bound = 0; int32_t tight_score_bound = 0; int32_t loose_score_bound = 0; @@ -379,8 +384,7 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, unmatched_below[parent_i - node_0] += unmatched_below[child_i - node_0]; } else { update_score(unmatched_below[child_i - node_0], &lower_bound, - &best_lower_bound, &tight_score_bound, - &best_upper_bound, &loose_score_bound, + &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); } } @@ -400,31 +404,31 @@ IntegerVector cpp_nni_distance(const IntegerMatrix& edge1, if (!matched_1[root_child_2 - node_0]) { update_score(int32_t(unmatched_below[root_node] + unmatched_1 - + unmatched_2), &lower_bound, &best_lower_bound, + + unmatched_2), &lower_bound, &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); } else { - update_score(unmatched_1, &lower_bound, &best_lower_bound, + update_score(unmatched_1, &lower_bound, &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); - update_score(unmatched_2, &lower_bound, &best_lower_bound, + update_score(unmatched_2, &lower_bound, &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); } } else { - update_score(unmatched_1, &lower_bound, &best_lower_bound, - &tight_score_bound, &best_upper_bound, &loose_score_bound, + update_score(unmatched_1, &lower_bound, + &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); } } else { - update_score(unmatched_below[root_node], &lower_bound, &best_lower_bound, + update_score(unmatched_below[root_node], &lower_bound, &tight_score_bound, &best_upper_bound, &loose_score_bound, &li_score_bound, &fack_score_bound); } return IntegerVector::create( Named("lower") = lower_bound, - _["best_lower"] = best_lower_bound, + _["best_lower"] = lower_bound, // No tighter lower bound is known _["tight_upper"] = tight_score_bound == NA_INT32 ? NA_INTEGER : tight_score_bound, _["best_upper"] = best_upper_bound, diff --git a/tests/testthat/test-tree_distance_nni.R b/tests/testthat/test-tree_distance_nni.R index 2259f4e8..80b62afa 100644 --- a/tests/testthat/test-tree_distance_nni.R +++ b/tests/testthat/test-tree_distance_nni.R @@ -57,7 +57,12 @@ test_that("Simple NNI approximations", { edge1 <- Postorder(tree1$edge) edge2 <- Postorder(tree2$edge) - Fack <- function(n) ((n - 2) * ceiling(log2(n))) + n + # Fack et al. 2002, theorem 3; N = internal nodes = n_edges + 1 + Fack <- function(n_edges) { + N <- n_edges + 1 + lgN <- ceiling(log2(N)) + N * lgN - 2 ^ lgN + 1 + 2 * (N - ceiling(log2(N + 1))) + } Sorting <- function(n) { lc <- ceiling(log2(n)) @@ -76,10 +81,10 @@ test_that("Simple NNI approximations", { best_upper = 0L, loose_upper = 0L, fack_upper = 0L, li_upper = 0L) oneUnmatched <- c(lower = 1L, best_lower = 1L, tight_upper = 1L, - best_upper = 1L, loose_upper = 2L, fack_upper = 2L, + best_upper = 1L, loose_upper = 1L, fack_upper = Fack(1), li_upper = Li(1)) - fiveUnmatched <- c(lower = 5L, best_lower = 10L, tight_upper = 10L, - best_upper = 10L, loose_upper = 18L, fack_upper = 18L, + fiveUnmatched <- c(lower = 5L, best_lower = 5L, tight_upper = 10L, + best_upper = 10L, loose_upper = 17L, fack_upper = Fack(5), li_upper = Li(5)) Test <- function(tree, expect) { @@ -87,8 +92,11 @@ test_that("Simple NNI approximations", { names(expectation) <- c("lower", "best_lower", "tight_upper", "best_upper", "loose_upper", "fack_upper", "li_upper") expectation[names(expect)] <- expect - if (is.na(expectation["best_lower"]) && !is.na(expect["tight_upper"])) { - expectation[c("best_lower", "best_upper")] <- expect["tight_upper"] + if (is.na(expectation["best_lower"])) { + expectation["best_lower"] <- expect["lower"] + } + if (is.na(expectation["best_upper"]) && !is.na(expect["tight_upper"])) { + expectation["best_upper"] <- expect["tight_upper"] } if (is.na(expect["loose_upper"])) { expectation["loose_upper"] <- min(expect[c("fack_upper", "li_upper")]) @@ -130,11 +138,11 @@ test_that("Simple NNI approximations", { # One region of three unmatched edges Test(read.tree(text="(((a, e), (c, d)), ((b, f), (g, h)));"), - c(lower = 3L, tight_upper = 5L, fack_upper = 8L, li_upper = Li(3))) + c(lower = 3L, tight_upper = 5L, fack_upper = Fack(3), li_upper = Li(3))) # One region of four unmatched edges Test(ape::read.tree(text="(((a, e), (f, d)), ((b, c), (g, h)));"), - c(lower = 4L, tight_upper = 7L, fack_upper = 14L, li_upper = Li(4))) + c(lower = 4L, tight_upper = 7L, fack_upper = Fack(4), li_upper = Li(4))) # One region of five unmatched edges Test(ape::read.tree(text="(((a, e), (f, d)), ((b, g), (c, h)));"), @@ -154,6 +162,57 @@ test_that("Simple NNI approximations", { }) +test_that("NNIDist() lower bounds are valid", { + # Exact NNI distance is 2 (ms609/TreeDist issue 220) + a <- read.tree(text = "(t1,(t2,(t4,t5)),(t3,t6));") + b <- read.tree(text = "(t6,(t1,(t4,t5)),(t3,t2));") + expect_equal(NNIDist(a, b)[c("lower", "best_lower")], + c(lower = 2L, best_lower = 2L)) +}) + +.ExpectBoundsBracketExact <- function(nTip) { + # Exact distances by breadth-first search over all trees, which are an NNI + # apart iff they differ in a single split + trees <- as.phylo(seq_len(NUnrooted(nTip)) - 1, nTip) + adjacent <- as.matrix(RobinsonFoulds(trees)) == 2 + exact <- matrix(NA_integer_, length(trees), length(trees)) + for (origin in seq_along(trees)) { + dist <- rep(NA_integer_, length(trees)) + dist[origin] <- 0L + frontier <- origin + step <- 0L + while (length(frontier)) { + step <- step + 1L + reached <- which(colSums(adjacent[frontier, , drop = FALSE]) > 0 & + is.na(dist)) + dist[reached] <- step + frontier <- reached + } + exact[origin, ] <- dist + } + testthat::expect_equal(max(exact), NNIDiameter(nTip)[, "exact"], + ignore_attr = TRUE) + + bounds <- lapply(NNIDist(trees), as.matrix) + exact <- exact[lower.tri(exact)] + .Bound <- function(field) bounds[[field]][lower.tri(bounds[[field]])] + testthat::expect_true(all(.Bound("lower") <= exact)) + testthat::expect_equal(.Bound("best_lower"), .Bound("lower")) + for (upper in c("tight_upper", "best_upper", "loose_upper", "fack_upper", + "li_upper")) { + testthat::expect_true(all(.Bound(upper) >= exact)) + } +} + +test_that("NNIDist() bounds bracket exact distances: 6 leaves", { + .ExpectBoundsBracketExact(6) +}) + +test_that("NNIDist() bounds bracket exact distances: 7 leaves", { + skip_on_cran() # c. 70 s + .ExpectBoundsBracketExact(7) +}) + test_that("NNI with lists of trees", { tree1 <- BalancedTree(1:8) list1 <- list(tree1, PectinateTree(as.character(1:8)), @@ -192,7 +251,10 @@ test_that("NNIDiameter() is sane", { FackMin <- function(n) ceiling(0.25 * lfactorial(n) / log(2)) exacts <- c(0, 0, 0, 1, 3, 5, 7, 10, 12, 15, 18, 21) liMaxes <- c(0, 1, 3, 5, 8, 13, 16, 21, 25, 31, 37, 43, 47, 53, 59, 65) - FackMax <- function(n) n * ceiling(log2(n)) + n - (2 * ceiling(log2(n))) + FackMax <- function(n) { + lgN <- ceiling(log2(n)) + n * lgN - 2 ^ lgN + 1 + 2 * (n - ceiling(log2(n + 1))) + } n <- 4:8 expect_equal(NNIDiameter(n), cbind( liMin = n - 3L,