|
| 1 | +#!/usr/bin/env python3 |
| 2 | +""" |
| 3 | +Test for rounding issues when values come from different sources. |
| 4 | +
|
| 5 | +The hypothesis: When Box A's max and Box B's min are computed or stored |
| 6 | +independently, float32 rounding can create gaps that don't exist in float64. |
| 7 | +""" |
| 8 | +import numpy as np |
| 9 | +from python_prtree import PRTree2D |
| 10 | + |
| 11 | + |
| 12 | +def test_computed_vs_literal(): |
| 13 | + """ |
| 14 | + Test when coordinates come from computations vs literals. |
| 15 | +
|
| 16 | + A computed value might round differently than a literal value |
| 17 | + due to intermediate precision. |
| 18 | + """ |
| 19 | + print("\n=== Test: Computed vs Literal Values ===\n") |
| 20 | + |
| 21 | + # Computed value (with intermediate float64 precision) |
| 22 | + computed_f64 = np.float64(1.0) / np.float64(3.0) * np.float64(300.0) # = 100.0 |
| 23 | + |
| 24 | + # Literal value |
| 25 | + literal_f64 = np.float64(100.0) |
| 26 | + |
| 27 | + print(f"Computed (f64): {computed_f64:.20f}") |
| 28 | + print(f"Literal (f64): {literal_f64:.20f}") |
| 29 | + print(f"Equal in f64: {computed_f64 == literal_f64}") |
| 30 | + |
| 31 | + computed_f32 = np.float32(computed_f64) |
| 32 | + literal_f32 = np.float32(literal_f64) |
| 33 | + |
| 34 | + print(f"\nComputed (f32): {computed_f32:.20f}") |
| 35 | + print(f"Literal (f32): {literal_f32:.20f}") |
| 36 | + print(f"Equal in f32: {computed_f32 == literal_f32}") |
| 37 | + |
| 38 | + # Create boxes |
| 39 | + boxes = np.array([ |
| 40 | + [0.0, 0.0, computed_f64, 100.0], # Ends at computed value |
| 41 | + [literal_f64, 0.0, 200.0, 100.0], # Starts at literal value |
| 42 | + ], dtype=np.float64) |
| 43 | + |
| 44 | + boxes_f32 = boxes.astype(np.float32) |
| 45 | + |
| 46 | + print(f"\nOverlap (f64): {boxes[0,2] >= boxes[1,0]}") |
| 47 | + print(f"Overlap (f32): {boxes_f32[0,2] >= boxes_f32[1,0]}") |
| 48 | + |
| 49 | + idx = np.array([0, 1], dtype=np.int64) |
| 50 | + |
| 51 | + tree_f64 = PRTree2D(idx, boxes) |
| 52 | + pairs_f64 = tree_f64.query_intersections() |
| 53 | + |
| 54 | + tree_f32 = PRTree2D(idx, boxes_f32) |
| 55 | + pairs_f32 = tree_f32.query_intersections() |
| 56 | + |
| 57 | + print(f"\nFloat64 tree: {len(pairs_f64)} pairs") |
| 58 | + print(f"Float32 tree: {len(pairs_f32)} pairs") |
| 59 | + |
| 60 | + if len(pairs_f64) != len(pairs_f32): |
| 61 | + print("\n❌ FALSE NEGATIVE!") |
| 62 | + return True |
| 63 | + return False |
| 64 | + |
| 65 | + |
| 66 | +def test_accumulated_computation(): |
| 67 | + """ |
| 68 | + Test when values are accumulated through multiple operations. |
| 69 | + """ |
| 70 | + print("\n=== Test: Accumulated Computation ===\n") |
| 71 | + |
| 72 | + # Create a value through accumulation |
| 73 | + accumulated_f64 = np.float64(0.0) |
| 74 | + step = np.float64(0.1) |
| 75 | + for i in range(1000): |
| 76 | + accumulated_f64 += step |
| 77 | + |
| 78 | + # Direct value |
| 79 | + direct_f64 = np.float64(100.0) |
| 80 | + |
| 81 | + print(f"Accumulated (f64): {accumulated_f64:.20f}") |
| 82 | + print(f"Direct (f64): {direct_f64:.20f}") |
| 83 | + print(f"Difference: {abs(accumulated_f64 - direct_f64):.20e}") |
| 84 | + |
| 85 | + accumulated_f32 = np.float32(accumulated_f64) |
| 86 | + direct_f32 = np.float32(direct_f64) |
| 87 | + |
| 88 | + print(f"\nAccumulated (f32): {accumulated_f32:.20f}") |
| 89 | + print(f"Direct (f32): {direct_f32:.20f}") |
| 90 | + print(f"Equal in f32: {accumulated_f32 == direct_f32}") |
| 91 | + |
| 92 | + # Create boxes with these values |
| 93 | + boxes_f64 = np.array([ |
| 94 | + [0.0, 0.0, accumulated_f64, 100.0], |
| 95 | + [direct_f64, 0.0, 200.0, 100.0], |
| 96 | + ], dtype=np.float64) |
| 97 | + |
| 98 | + boxes_f32 = boxes_f64.astype(np.float32) |
| 99 | + |
| 100 | + print(f"\nOverlap (f64): {boxes_f64[0,2] >= boxes_f64[1,0]}") |
| 101 | + print(f"Overlap (f32): {boxes_f32[0,2] >= boxes_f32[1,0]}") |
| 102 | + |
| 103 | + idx = np.array([0, 1], dtype=np.int64) |
| 104 | + |
| 105 | + tree_f64 = PRTree2D(idx, boxes_f64) |
| 106 | + pairs_f64 = tree_f64.query_intersections() |
| 107 | + |
| 108 | + tree_f32 = PRTree2D(idx, boxes_f32) |
| 109 | + pairs_f32 = tree_f32.query_intersections() |
| 110 | + |
| 111 | + print(f"\nFloat64 tree: {len(pairs_f64)} pairs") |
| 112 | + print(f"Float32 tree: {len(pairs_f32)} pairs") |
| 113 | + |
| 114 | + if len(pairs_f64) != len(pairs_f32): |
| 115 | + print("\n❌ FALSE NEGATIVE!") |
| 116 | + return True |
| 117 | + return False |
| 118 | + |
| 119 | + |
| 120 | +def test_separate_float32_arrays(): |
| 121 | + """ |
| 122 | + Test when float32 values are created in separate arrays. |
| 123 | +
|
| 124 | + Key insight: If two float32 values are created independently, |
| 125 | + they might have different representations even if they should be equal. |
| 126 | + """ |
| 127 | + print("\n=== Test: Separate Float32 Arrays ===\n") |
| 128 | + |
| 129 | + # Create a problematic float64 value |
| 130 | + problematic = np.float64(100.0) + np.float64(1e-7) |
| 131 | + |
| 132 | + print(f"Problematic value (f64): {problematic:.20f}") |
| 133 | + |
| 134 | + # Create first array with this value as max |
| 135 | + array1_f32 = np.array([0.0, 0.0, problematic, 100.0], dtype=np.float32) |
| 136 | + |
| 137 | + # Create second array with this value as min |
| 138 | + array2_f32 = np.array([problematic, 0.0, 200.0, 100.0], dtype=np.float32) |
| 139 | + |
| 140 | + print(f"\nArray1[2] (max): {array1_f32[2]:.20f}") |
| 141 | + print(f"Array2[0] (min): {array2_f32[0]:.20f}") |
| 142 | + print(f"Equal: {array1_f32[2] == array2_f32[0]}") |
| 143 | + print(f"Overlap: {array1_f32[2] >= array2_f32[0]}") |
| 144 | + |
| 145 | + # Combine into boxes |
| 146 | + boxes_f32 = np.vstack([array1_f32.reshape(1, -1), array2_f32.reshape(1, -1)]) |
| 147 | + |
| 148 | + print(f"\nCombined boxes (f32):\n{boxes_f32}") |
| 149 | + |
| 150 | + idx = np.array([0, 1], dtype=np.int64) |
| 151 | + |
| 152 | + tree = PRTree2D(idx, boxes_f32) |
| 153 | + pairs = tree.query_intersections() |
| 154 | + |
| 155 | + print(f"\nIntersections found: {len(pairs)}") |
| 156 | + print(f"Pairs: {pairs}") |
| 157 | + |
| 158 | + # Expected: should find intersection since they touch |
| 159 | + if len(pairs) == 0: |
| 160 | + print("\n❌ FALSE NEGATIVE: Touching boxes not detected!") |
| 161 | + return True |
| 162 | + |
| 163 | + return False |
| 164 | + |
| 165 | + |
| 166 | +def test_binary_representation(): |
| 167 | + """ |
| 168 | + Test values that have identical decimal representation but different binary. |
| 169 | + """ |
| 170 | + print("\n=== Test: Binary Representation ===\n") |
| 171 | + |
| 172 | + # Create a value that cannot be exactly represented in binary |
| 173 | + decimal_val = 0.1 |
| 174 | + |
| 175 | + # In float64 |
| 176 | + val_f64 = np.float64(decimal_val) |
| 177 | + print(f"0.1 in float64: {val_f64:.60f}") |
| 178 | + print(f"Hex: {val_f64.hex()}") |
| 179 | + |
| 180 | + # In float32 |
| 181 | + val_f32 = np.float32(decimal_val) |
| 182 | + print(f"\n0.1 in float32: {val_f32:.60f}") |
| 183 | + print(f"Hex: {val_f32.hex()}") |
| 184 | + |
| 185 | + # Scale up |
| 186 | + scale = 1000 |
| 187 | + scaled_f64 = val_f64 * scale |
| 188 | + scaled_f32 = val_f32 * scale |
| 189 | + |
| 190 | + print(f"\nScaled (f64): {scaled_f64:.60f}") |
| 191 | + print(f"Scaled (f32): {scaled_f32:.60f}") |
| 192 | + |
| 193 | + # Now use these as coordinates |
| 194 | + boxes_f64 = np.array([ |
| 195 | + [0.0, 0.0, scaled_f64, 100.0], |
| 196 | + [scaled_f64, 0.0, 200.0, 100.0], |
| 197 | + ], dtype=np.float64) |
| 198 | + |
| 199 | + # Create float32 version two ways: |
| 200 | + # 1. Convert from float64 |
| 201 | + boxes_f32_converted = boxes_f64.astype(np.float32) |
| 202 | + |
| 203 | + # 2. Create directly with float32 |
| 204 | + boxes_f32_direct = np.array([ |
| 205 | + [0.0, 0.0, val_f32 * scale, 100.0], |
| 206 | + [val_f32 * scale, 0.0, 200.0, 100.0], |
| 207 | + ], dtype=np.float32) |
| 208 | + |
| 209 | + print(f"\nBoxes f32 (converted):\n{boxes_f32_converted}") |
| 210 | + print(f"\nBoxes f32 (direct):\n{boxes_f32_direct}") |
| 211 | + print(f"\nAre they equal? {np.array_equal(boxes_f32_converted, boxes_f32_direct)}") |
| 212 | + |
| 213 | + idx = np.array([0, 1], dtype=np.int64) |
| 214 | + |
| 215 | + tree_f64 = PRTree2D(idx, boxes_f64) |
| 216 | + pairs_f64 = tree_f64.query_intersections() |
| 217 | + |
| 218 | + tree_f32_conv = PRTree2D(idx, boxes_f32_converted) |
| 219 | + pairs_f32_conv = tree_f32_conv.query_intersections() |
| 220 | + |
| 221 | + tree_f32_dir = PRTree2D(idx, boxes_f32_direct) |
| 222 | + pairs_f32_dir = tree_f32_dir.query_intersections() |
| 223 | + |
| 224 | + print(f"\nFloat64 tree: {len(pairs_f64)} pairs") |
| 225 | + print(f"Float32 (converted): {len(pairs_f32_conv)} pairs") |
| 226 | + print(f"Float32 (direct): {len(pairs_f32_dir)} pairs") |
| 227 | + |
| 228 | + if len(pairs_f64) != len(pairs_f32_conv) or len(pairs_f64) != len(pairs_f32_dir): |
| 229 | + print("\n❌ FALSE NEGATIVE!") |
| 230 | + return True |
| 231 | + |
| 232 | + return False |
| 233 | + |
| 234 | + |
| 235 | +if __name__ == "__main__": |
| 236 | + print("=" * 70) |
| 237 | + print("Testing Different Sources of Rounding") |
| 238 | + print("=" * 70) |
| 239 | + |
| 240 | + issue_found = False |
| 241 | + |
| 242 | + issue_found |= test_computed_vs_literal() |
| 243 | + issue_found |= test_accumulated_computation() |
| 244 | + issue_found |= test_separate_float32_arrays() |
| 245 | + issue_found |= test_binary_representation() |
| 246 | + |
| 247 | + print("\n" + "=" * 70) |
| 248 | + if issue_found: |
| 249 | + print("❌ FALSE NEGATIVE CONFIRMED!") |
| 250 | + else: |
| 251 | + print("⚠️ No false negatives in these tests") |
| 252 | + print("=" * 70) |
0 commit comments