update readme

README.Rmd

@@ -37,28 +37,73 @@ After that, you can install the development version of `dyingforacup` like so:
 remotes::install_github('andrewGhazi/dyingforacup', type = "source")
 ```
 
-## Example
+## 1D Example animation
+
+The point of this package is to suggest coffee brewing configurations in brew parameter space that balance A) improving the expected rating and B) exploring the space. There are a **lot** of dials to turn when brewing coffee, and it's practically impossible to try every combination of grind size, temperature, bloom time, filter thickness, etc.
+
+Say you only had one brew parameter: the coarseness dial on the grinder. Imagine the true, unknown relationship between grinder setting and coffee quality looks like this:
+
+![](man/figures/bo_anim0-fs8.png)
+
+You have one starting observation at 4.5 (too fine, so the observed rating is low). There's a bit of noise about the true function. What setting should you try next?
+
+If you use Gaussian processes and [Bayesian Optimization](https://www.youtube.com/watch?v=wZODGJzKmD0), you can get suggestions that are, in a sense, optimal. Let's see how automated suggestions work out:
+
+![](man/figures/vid.mp4)
+
+There's a lot going on in that video. This explains each panel:
+
+* Top panel: 
+  * dark grey line: same true (unknown) relationship between grinder setting and coffee quality as above
+  * points: noisy observed ratings about that function
+  * light grey lines: posterior draws for the relationship from the GP 
+  * grey ribbon: 90% predictive interval of a new coffee brewed at the given setting
+  * red arrow & dashed line: the previous best observation shifted down by a set amount (`offset`, the length of the red arrow). 
+* 2nd panel, blue function: fit sd, the sd of the light grey lines
+* 3rd panel, red function: expected improvement: the expectation value of a new coffee at the given setting falling above the dashed red line.
+* bottom pane, green function: the acquisition function, a weighted mixture of the two curves above. `lambda * blue + (1-lambda) * red = green`
+  * orange diamond: maximum point on the acquisition function = next suggested point
+
+You can see that first it suggests a very high setting because that's where there's the most uncertainty given the first point. After that turns out badly as well, it tries in the middle. That does much better, after which it hones in on the global maximum (it gets somewhat lucky and finds a near optimal point at only the fourth suggestion). After that it tries elsewhere in the space, collapsing uncertainty wherever it's high to see if there's some other hidden peak. 
+
+## Usage
+
+As the last frame of the video suggests, this process can be extended to an arbitrary number of brew parameters. Bear in mind that this isn't magic, and finding optima in higher dimensional spaces will require many more observations. This is especially true if the ratings are noisy, so try hard to give each cup a fair, normally-distributed rating. Speaking of, integer ratings of 0-10 aren't allowed, the ratings have to be normally distributed. That might change if I feel like implementing it.
 
 Give the `suggest_next()` function a data frame of brew parameters with ratings and it will suggest a point to try next that has high predicted probability of improving the rating.
 
 ```{r eval=FALSE}
 library(dyingforacup)
-options(mc.cores = 4)
+
+options(mc.cores = 4, digits = 3)
 
 dat = data.frame(grinder_setting = c(  8,    7,   9), 
                  temp            = c(193,  195, 179),
                  bloom_time      = c( 25,   20,  45),
                  rating          = c(1.1, -0.7,  -1))
 
 suggest_next(dat,
-             iter_sampling = 4000, 
+             iter_sampling = 1000, 
              refresh = 0, 
+             offset = .33, 
+             lambda = .1,
              show_exceptions = FALSE, 
-             adapt_delta = .95, 
              parallel_chains = 4)
 ```
 
-## 1D Example animation
+```
+$draws_df
+...
+
+$acq_df
+...
+
+$suggested
+   post_sd exp_imp   acq grinder_setting  temp bloom_time
+     <num>   <num> <num>           <num> <num>      <num>
+1:    1.19   0.464 0.536               9   195         30
+```
+This returns a list of MCMC draws, the acquisition function values over a grid of brew parameters, and a suggestion on where to go next. `offset` and `lambda` can be tweaked to control exploration vs exploitation, but expect to be suggested some combinations that result in really bad coffee sometimes.  See `?suggest_next` for more detail on these and more function arguments.
 
 ## TODO list
 
@@ -70,7 +115,7 @@ Easy:
 Medium:
 
 * Non-normal outcome
-* Fast GP approximations for 1D/2D datasets with [`gptools`](https://github.com/onnela-lab/gptools/tree/main)
+* Fast GP approximations for 1D/2D datasets with [`gptools`](https://github.com/onnela-lab/gptools/tree/main) or [Hilbert spaces](https://arxiv.org/abs/2004.11408)
 
 Hard:
 

README.md

@@ -30,6 +30,26 @@ like so:
 remotes::install_github('andrewGhazi/dyingforacup', type = "source")
 ```
 
+## 1D Example animation
+
+The point of this package is to suggest coffee brewing configurations in
+brew parameter space that balance A) improving the expected rating and
+B) exploring the space. Say you only had one brew parameter: the
+coarseness dial on the grinder. Imagine the true, unknown relationship
+between grinder setting and coffee quality looks like this:
+
+![](man/figures/bo_anim0-fs8.png)
+
+You have one starting observation at 4.5 (too fine). There’s a bit of
+noise about the true function. What setting should you try next?
+
+If you use [Bayesian
+Optimization](https://www.youtube.com/watch?v=wZODGJzKmD0), you can
+balance exploration and exploitation. Let’s see how automated
+suggestions work out:
+
+![](man/figures/vid.mp4)
+
 ## Example
 
 Give the `suggest_next()` function a data frame of brew parameters with
@@ -53,8 +73,6 @@ suggest_next(dat,
              parallel_chains = 4)
 ```
 
-## 1D Example animation
-
 ## TODO list
 
 Easy:
@@ -65,7 +83,8 @@ Easy:
 Medium:
 
 - Non-normal outcome
-- Fast GP approximations for 1D/2D datasets
+- Fast GP approximations for 1D/2D datasets with
+  [`gptools`](https://github.com/onnela-lab/gptools/tree/main)
 
 Hard:
 
@@ -76,3 +95,4 @@ Nightmare:
 
 - Fast GP approximations for 3D+
   - I think this would require writing my own ND FFT function?
+- Refactor to use INLA (preferably from scratch over `R-INLA`)