finish spacebook

aaaa/spacebook/README.md

@@ -13,4 +13,69 @@
 ## Writeup
 by [erin (`barzamin`)](https://imer.in).
 
-![Distribution of star magnitudes for catalog and unknown sets](magnitude-distribution.png)
+SpaceBook is a very similar problem to Attitude Adjustment; we still need to solve an orthogonal Procrustes problem to align a set of stars to a catalog. However, we don't know which stars in the unknown set match to which stars in the catalog, unlike in Attitude Adjustment.
+
+We're given a boresight vector $\vec{v}\in\mathbb{R}^3$ for each star, and a magnitude $m$. The vectors in the catalog set are off by some rotation from the vectors in the unknown set; the magnitudes match. Looking at the distribution of magnitudes, we see that there is a significant outlier population in both the catalog and reference sets, and that these outlier populations match:
+
+![Distribution of star magnitudes for catalog and unknown sets](spacebook-magnitude-distribution.png)
+
+*Note that code from here on in this report is reconstructed, since I accidentally overwrote the original copy. I took the chance to use scipy for Kabsch this time instead of the original package I used, [`rmsd`](https://github.com/charnley/rmsd).*
+
+My immediate suspicion was that we could take the outliers (simply by thresholding for stars with $m>500$) and, for each, directly match it to the catalog by finding the star with closest magnitude in the catalog:
+
+```{.python}
+def match_brightest_stars(unknown, catalog):
+    brightest = []
+    catalog_matches = []
+    for star in unknown:
+        if star['m'] < 500:
+            continue
+
+        # find closest magnitude in catalog
+        match = np.argmin(np.abs(np.array([s['m'] for s in catalog]) - star['m']))
+        print(f'[+] matched {star} to {catalog[match]}')
+        brightest.append(star)
+        catalog_matches.append(catalog[match])
+
+    return brightest, catalog_matches
+```
+
+We can then align those brightest stars to their closest magnitude matches using the Kabsch algorithm:
+```{.python}
+def orient(ref, catalog):
+    P = [s['v'] for s in catalog]
+    Q = [s['v'] for s in ref]
+
+    rot, rmsd = Rotation.align_vectors(P, Q)
+    print(f'[+] aligned; rmsd = {rmsd}')
+    return rot
+```
+
+Then for each unknown star with vector $\vec{v}$, we can find its catalog index by evaluating
+$$\arg\min_{\vec{u}\in\text{catalog}} \norm{\vec{u}-\vec{v}}_2$$
+for every challenge:
+```{.python}
+for _ in range(5):
+    unknown_stars = read_starfile(r.recvuntil('\n\n').decode())
+
+    # find pairs
+    unknown_ref, catalog_ref = match_brightest_stars(unknown_stars, catalog)
+    # get attitude
+    attitude = orient(unknown_ref, catalog_ref)
+    # rotate each star to catalog-referenced coordinates and match by L2 norm
+    index_guesses = []
+    for star in unknown_stars:
+        v = attitude.apply(star['v'])
+        index_guesses.append(np.argmin([np.linalg.norm(catalog_star['v'] - v) for catalog_star in catalog]))
+
+    r.send(','.join(map(str, index_guesses))+'\n')
+    r.recvuntil('Left...\n')
+print(r.clean())
+```
+
+### Full code
+```{.python include=spacebook.py}
+```
+
+## Resources
+- https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.align_vectors.html