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