Skip to content
Snippets Groups Projects
Commit e2165feb authored by Samuel Grayson's avatar Samuel Grayson
Browse files

Final commit

parent b89bb190
No related branches found
No related tags found
No related merge requests found
......@@ -61,35 +61,58 @@ The quantities in [square brackets] in the last line can be precomputed for all
<p><strong>Algorithm 1</strong> <code>best_order</code> should compute the best order to the tests, using the assumptions above.
<div><pre><code class="language-python">import numpy.typing as npt
<div><pre><code class="language-python">from typing import List
import numpy.typing as npt
import numpy as np
def reverse(array: npt.NDArray[np.float64]) -&gt; npt.NDArray[np.float64]:
"""Reverse a numpy array"""
return array[::-1]
def cumprod_noninclusive(array: npt.NDArray[np.float64]) -&gt; npt.NDArray[np.float64]:
ret = np.ones_like(array)
for i in range(1, len(array)):
ret[i] = ret[i-1] * array[i-1]
def cumsum(
array: npt.NDArray[np.float64],
inclusive: bool = True,
one_larger_output: bool = False,
) -&gt; npt.NDArray[np.float64]:
if one_larger_output:
assert not inclusive
ret = np.zeros(len(array) + int(one_larger_output))
for i in range(int(not inclusive), len(ret)):
ret[i] = ret[i-1] + array[i - int(not inclusive)]
return ret
def cumprod(
array: npt.NDArray[np.float64],
inclusive: bool = True,
one_larger_output: bool = False,
) -&gt; npt.NDArray[np.float64]:
if one_larger_output:
assert not inclusive
ret = np.ones(len(array) + int(one_larger_output))
for i in range(int(not inclusive), len(ret)):
ret[i] = ret[i-1] * array[i - int(not inclusive)]
return ret
def next_best(nt: int, d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]):
def next_best(d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]) -&gt; List[int]:
"""
Returns the index of the best test to run.
nt, d, and p should already include a fictitious initial and fictitious final test.
"""
# Compute D and P for the unmodified order.
# np.cumsum is inclusive, D[k] = d[0] + d[1] + ... + d[k]
# cumprod_noninclusive is not, P[k] = p[0] * p[1] * ... * p[k-1] * (1 - p[k])
D = np.cumsum(d)
P = cumprod_noninclusive(p) * (1 - p)
# Precompute intermediate values.
A = np.cumsum(D * P) - d[0] * p[0]
B = np.cumsum(P) - p[0]
if len(d) == 1:
return 0
assert len(d) == len(p)
# D[ell] := d[0] + d[1] + ... + d[ell]
D = cumsum(d)
# P[ell] := p[0] * p[1] * ... * p[ell-1] * (1 - p[ell])
P = cumprod(p, inclusive=False) * (1 - p)
# A[ell] := D[1] * P[1] + D[2] * P[2] + ... + D[ell - 1] * P[ell - 1] for ell &gt; 0
A = cumsum(D * P, inclusive=False) - D[0] * P[0]
# B[ell] := P[1] + P[2] + ... + P[ell - 1] for ell &gt; 0
B = cumsum(P, inclusive=False) - P[0]
# C[ell] := D[ell+1] * P[ell+1] + D[ell+2] * P[ell+2] + ... + D[len(D) - 1] * P[len(D) - 1]
# TODO: this
C = np.cumsum(reverse(D * P))
# Compute E of the swapped order 0 &lt;-&gt; k for all k.
......@@ -102,42 +125,33 @@ def next_best(nt: int, d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]):
+ C
)
# Find min, excluding the trailing fictitious test, which would otherwise be optimal.
return np.argmin(E_of_swap[:-1])
return np.argmin(E_of_swap)
def best_order(
nt: int,
d: npt.NDArray[np.float64],
p: npt.NDArray[np.float64],
) -&gt; npt.NDArray[np.float64]:
) -&gt; npt.NDArray[int]:
"""
Returns a list of indices which consist the best order to run.
nt, d, and p should not contain fictitious tests.
"""
order = np.zeros(nt)
assert len(d) == len(p)
# Add the fictitious test to these.
nt = nt + 2
d = np.concatenate([0], d, [0])
p = np.concatenate([1], p, [0])
tests = np.arange(len(d), dtype=int)
# Loop through each test and find the best order.
for k in range(nt - 1):
for k in range(len(d)):
# The first k have already been assigned.
best = next_best(nt - k, d[k:], p[k:])
order[k] = best
best = next_best(d[k:], p[k:])
# Swap best with k + 1, the first real test in the next round.
np.swap(p, k + 1, best)
np.swap(d, k + 1, best)
# Swap best with k, so that arr[k+1:] are the remaining tests
if k != best:
p[[k, best]] = p[[best, k]]
d[[k, best]] = d[[best, k]]
tests[[k, best]] = tests[[best, k]]
# Set the first test in the next round to be fictitious.
# This deletes the test we select for this round.
d[k + 1] = 0
p[k + 1] = 1
assert list(sorted(tests)) == list(range(len(d)))
return order
return tests
</code></pre></div></p>
</body></html>
\ No newline at end of file
◊begin[
(require
txexpr
"lib/txexpr.rkt")
(define-values (rest template-data) (splitf-txexpr* doc 'template-data))
(define title (get-child template-data 'title #:default '("Working Title")))
(define head (get-child template-data 'head #:default '()))
(define body (get-elements rest))
]
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8" />
<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
<script
id="MathJax-script"
async
src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
></script>
</head>
<body>
◊(->html doc #:splice? #t)
</body>
</html>
◊(->html
`(html ((lang "en"))
(head
(meta ((charset "utf-8")))
(title ,@title)
,@head)
(body
(h1 ,@title)
,@body)))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment