| 78 | |
| 79 | @pytest.mark.slow |
| 80 | def test_diophantine_fuzz(): |
| 81 | # Fuzz test the diophantine solver |
| 82 | rng = np.random.RandomState(1234) |
| 83 | |
| 84 | max_int = np.iinfo(np.intp).max |
| 85 | |
| 86 | for ndim in range(10): |
| 87 | feasible_count = 0 |
| 88 | infeasible_count = 0 |
| 89 | |
| 90 | min_count = 500 // (ndim + 1) |
| 91 | |
| 92 | while min(feasible_count, infeasible_count) < min_count: |
| 93 | # Ensure big and small integer problems |
| 94 | A_max = 1 + rng.randint(0, 11, dtype=np.intp)**6 |
| 95 | U_max = rng.randint(0, 11, dtype=np.intp)**6 |
| 96 | |
| 97 | A_max = min(max_int, A_max) |
| 98 | U_max = min(max_int - 1, U_max) |
| 99 | |
| 100 | A = tuple(int(rng.randint(1, A_max + 1, dtype=np.intp)) |
| 101 | for j in range(ndim)) |
| 102 | U = tuple(int(rng.randint(0, U_max + 2, dtype=np.intp)) |
| 103 | for j in range(ndim)) |
| 104 | |
| 105 | b_ub = min(max_int - 2, sum(a * ub for a, ub in zip(A, U))) |
| 106 | b = int(rng.randint(-1, b_ub + 2, dtype=np.intp)) |
| 107 | |
| 108 | if ndim == 0 and feasible_count < min_count: |
| 109 | b = 0 |
| 110 | |
| 111 | X = solve_diophantine(A, U, b) |
| 112 | |
| 113 | if X is None: |
| 114 | # Check the simplified decision problem agrees |
| 115 | X_simplified = solve_diophantine(A, U, b, simplify=1) |
| 116 | assert_(X_simplified is None, (A, U, b, X_simplified)) |
| 117 | |
| 118 | # Check no solution exists (provided the problem is |
| 119 | # small enough so that brute force checking doesn't |
| 120 | # take too long) |
| 121 | ranges = tuple(range(0, a * ub + 1, a) for a, ub in zip(A, U)) |
| 122 | |
| 123 | size = 1 |
| 124 | for r in ranges: |
| 125 | size *= len(r) |
| 126 | if size < 100000: |
| 127 | assert_(not any(sum(w) == b for w in itertools.product(*ranges))) |
| 128 | infeasible_count += 1 |
| 129 | else: |
| 130 | # Check the simplified decision problem agrees |
| 131 | X_simplified = solve_diophantine(A, U, b, simplify=1) |
| 132 | assert_(X_simplified is not None, (A, U, b, X_simplified)) |
| 133 | |
| 134 | # Check validity |
| 135 | assert_(sum(a * x for a, x in zip(A, X)) == b) |
| 136 | assert_(all(0 <= x <= ub for x, ub in zip(X, U))) |
| 137 | feasible_count += 1 |