Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
10094 lines (9163 sloc) 410 KB
title author output editor_options
Food for Thought
name email twitter
Amanda Dobbyn
amanda.e.dobbyn@gmail.com
dobbleobble
html_document github_document pdf_document
keep_md toc theme
true
true
yeti
toc
true
keep_tex toc
true
false
chunk_output_type
inline

dirs <- c("prep", "build", "score", "scrape", "solve", "simulate")
paths <- stringr::str_c("./scripts/", dirs)

# Import all .R scripts from all the dirs above 
for (p in paths) {
  suppressMessages(suppressPackageStartupMessages(dobtools::import_scripts(p)))
}
## Warning in read_fun(path = path, sheet = sheet, limits = limits, shim =
## shim, : partial argument match of 'sheet' to 'sheet_i'

About

This is an ongoing project on menu optimizing. It's mainly an excuse for me to use several data science techniques in various proportions: along the way I query an API, generate menus, solve them algorithmically, simulate solving them, scrape the web for real menus, and touch on some natural language processing techniques. Don't worry about getting too hungry: this project has been fairly nicknamed slandered "Eat, Pray, Barf."

The meat of the project surrounds building menus and changing them until they are in compliance with daily nutritional guidelines. We'll simulate the curve of the proportion of these that are solvable as we increase the minimum portion size that each food must meet. Finally, I start about trying to improve the quality of the menus (i.e., decrease barf factor) by taking a cue from actual recipes scraped from Allrecipes.com.

Getting from A to Beef

The data we'll be using here is conveniently located in an Excel file called ABBREV.xlsx on the USDA website. As the name suggests, this is an abbreviated version of all the foods in their database.

If you do want the full list, they provide a Microsoft Access SQL dump as well (which requires that you have Access). The USDA also does have an open API so you can create an API key and grab foods from them with requests along the lines of a quick example I'll go through. The API documentation walks through the format for requesting data in more detail. I'll walk through an example of how to get some foods and a few of their associated nutrient values.

The base URL we'll want is http://api.nal.usda.gov/ndb/.

The default number of results per request is 50 so we specify 1500 as our max. In this example I set subset to 1 in order to grab the most common foods. (Otherwise 1:1500 query only gets you from a to beef 😆.) If you do want to grab all foods, you can send requests of 1500 iteratively specifying offset, which refers to the number of the first row you want, and then glue them together. We've specified just 4 nutrient values we want here: calories, sugar, lipids, and carbohydrates.

After attaching those parameters to the end of our base URL, we'd have:

http://api.nal.usda.gov/ndb/nutrients/?format=json&api_key="<YOUR_KEY_HERE>&subset=1&max=1500&nutrients=205&nutrients=204&nutrients=208&nutrients=269

In the browser, you could paste that same thing in to see:

We'll use the jsonlite package to turn that fromJSON into an R object.

foods_raw <- jsonlite::fromJSON(paste0("http://api.nal.usda.gov/ndb/nutrients/?format=json&api_key=", 
                       key, "&subset=1&max=1500&nutrients=205&nutrients=204&nutrients=208&nutrients=269"), flatten = FALSE)

foods <- as_tibble(foods_raw$report$foods)
head(foods) %>% kable(format = "html")
ndbno name weight measure nutrients
14007 Alcoholic beverage, beer, light, BUD LIGHT 29.5 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 9, --, 0.00, 0.38, 29, --, 0, 1.3
14009 Alcoholic beverage, daiquiri, canned 30.5 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 38, --, 0.00, 4.79, 125, --, 0, 15.7
14534 Alcoholic beverage, liqueur, coffee, 63 proof 34.8 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 107, --, 0.10, 11.21, 308, --, 0.3, 32.2
14015 Alcoholic beverage, pina colada, canned 32.6 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 77, --, 2.48, 9.00, 237, --, 7.6, 27.6
14019 Alcoholic beverage, tequila sunrise, canned 31.1 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 34, --, 0.03, 3.51, 110, --, 0.1, 11.3
14027 Alcoholic beverage, whiskey sour, canned 30.8 1.0 fl oz 208, 269, 204, 205, Energy, Sugars, total, Total lipid (fat), Carbohydrate, by difference, kcal, g, g, g, 37, --, 0.00, 4.13, 119, --, 0, 13.4

We've got one row per food and a nested list-col of nutrients.

str(foods$nutrients[1:3])
## List of 3
##  $ :'data.frame':	4 obs. of  5 variables:
##   ..$ nutrient_id: chr [1:4] "208" "269" "204" "205"
##   ..$ nutrient   : chr [1:4] "Energy" "Sugars, total" "Total lipid (fat)" "Carbohydrate, by difference"
##   ..$ unit       : chr [1:4] "kcal" "g" "g" "g"
##   ..$ value      : chr [1:4] "9" "--" "0.00" "0.38"
##   ..$ gm         : chr [1:4] "29" "--" "0" "1.3"
##  $ :'data.frame':	4 obs. of  5 variables:
##   ..$ nutrient_id: chr [1:4] "208" "269" "204" "205"
##   ..$ nutrient   : chr [1:4] "Energy" "Sugars, total" "Total lipid (fat)" "Carbohydrate, by difference"
##   ..$ unit       : chr [1:4] "kcal" "g" "g" "g"
##   ..$ value      : chr [1:4] "38" "--" "0.00" "4.79"
##   ..$ gm         : chr [1:4] "125" "--" "0" "15.7"
##  $ :'data.frame':	4 obs. of  5 variables:
##   ..$ nutrient_id: chr [1:4] "208" "269" "204" "205"
##   ..$ nutrient   : chr [1:4] "Energy" "Sugars, total" "Total lipid (fat)" "Carbohydrate, by difference"
##   ..$ unit       : chr [1:4] "kcal" "g" "g" "g"
##   ..$ value      : chr [1:4] "107" "--" "0.10" "11.21"
##   ..$ gm         : chr [1:4] "308" "--" "0.3" "32.2"

If we tried to unnest this right now we'd get an error.

foods %>% unnest()   # error :(

That's because missing values are coded as --.

foods$nutrients[[100]] %>% as_tibble() 
## # A tibble: 4 x 5
##   nutrient_id                    nutrient  unit value    gm
## *       <chr>                       <chr> <chr> <chr> <chr>
## 1         208                      Energy  kcal    --    --
## 2         269               Sugars, total     g    --    --
## 3         204           Total lipid (fat)     g  0.00     0
## 4         205 Carbohydrate, by difference     g    --    --

This becomes an issue for two of these columns, gm and value because gm gets coded as type numeric if there are no musing values and character otherwise. Consider the case where we have no missing values: here we see that gm is numeric.

foods$nutrients[[200]] %>% as_tibble()
## # A tibble: 4 x 5
##   nutrient_id                    nutrient  unit value     gm
## *       <chr>                       <chr> <chr> <chr>  <dbl>
## 1         208                      Energy  kcal   216 540.00
## 2         269               Sugars, total     g 17.00  42.50
## 3         204           Total lipid (fat)     g 12.00  30.00
## 4         205 Carbohydrate, by difference     g 23.98  59.95

We can't unnest yet because a single column in a dataframe can only have values of one type; without changing the types of the various gm columns to a lowest common denominator, we won't be able to combine them.

We can replace our --s with NAs no problem

foods$nutrients <- foods$nutrients %>% 
  map(na_if, "--") 

but unnesting is the challenge.

foods %>% unnest()
## Error in bind_rows_(x, .id): Column `gm` can't be converted from character to numeric

The naive approach of mapping character over all the nutrients columns doesn't give us the output we expect

foods$nutrients[1] %>% map(as.character)
## [[1]]
## [1] "c(\"208\", \"269\", \"204\", \"205\")"                                                   
## [2] "c(\"Energy\", \"Sugars, total\", \"Total lipid (fat)\", \"Carbohydrate, by difference\")"
## [3] "c(\"kcal\", \"g\", \"g\", \"g\")"                                                        
## [4] "c(\"9\", NA, \"0.00\", \"0.38\")"                                                        
## [5] "c(\"29\", NA, \"0\", \"1.3\")"

and we're stuck with the usual sense of not quite being able to reach the part of the data we want. (All credit here to the fantastic Jenny Bryan.)

So instead we'll dive into the second level of our nested list, take everything in there to character, and then unnest.

foods$nutrients <- foods$nutrients %>% modify_depth(2, as.character)

Which is a nicer way of saying

for (i in 1:length(foods$nutrients)) {
  for (j in 1:nrow(foods$nutrients[[1]])) {
    foods$nutrients[[i]]$nutrient_id[j] <- as.character(foods$nutrients[[i]]$nutrient_id[j])
    foods$nutrients[[i]]$nutrient[j] <- as.character(foods$nutrients[[i]]$nutrient[j])
    foods$nutrients[[i]]$unit[j] <- as.character(foods$nutrients[[i]]$unit[j])
    foods$nutrients[[i]]$gm[j] <- as.character(foods$nutrients[[i]]$gm[j])
    foods$nutrients[[i]]$value[j] <- as.character(foods$nutrients[[i]]$value[j])
  }
}

Now we can unnest the whole thing.

foods <- foods %>% unnest()

Finally, let's set value and gm to numeric.

foods$value <- as.numeric(foods$value)
foods$gm <- as.numeric(foods$gm)
foods[1:20, ] %>% kable(format = "html")
ndbno name weight measure nutrient_id nutrient unit value gm
14007 Alcoholic beverage, beer, light, BUD LIGHT 29.5 1.0 fl oz 208 Energy kcal 9.00 29.0
14007 Alcoholic beverage, beer, light, BUD LIGHT 29.5 1.0 fl oz 269 Sugars, total g NA NA
14007 Alcoholic beverage, beer, light, BUD LIGHT 29.5 1.0 fl oz 204 Total lipid (fat) g 0.00 0.0
14007 Alcoholic beverage, beer, light, BUD LIGHT 29.5 1.0 fl oz 205 Carbohydrate, by difference g 0.38 1.3
14009 Alcoholic beverage, daiquiri, canned 30.5 1.0 fl oz 208 Energy kcal 38.00 125.0
14009 Alcoholic beverage, daiquiri, canned 30.5 1.0 fl oz 269 Sugars, total g NA NA
14009 Alcoholic beverage, daiquiri, canned 30.5 1.0 fl oz 204 Total lipid (fat) g 0.00 0.0
14009 Alcoholic beverage, daiquiri, canned 30.5 1.0 fl oz 205 Carbohydrate, by difference g 4.79 15.7
14534 Alcoholic beverage, liqueur, coffee, 63 proof 34.8 1.0 fl oz 208 Energy kcal 107.00 308.0
14534 Alcoholic beverage, liqueur, coffee, 63 proof 34.8 1.0 fl oz 269 Sugars, total g NA NA
14534 Alcoholic beverage, liqueur, coffee, 63 proof 34.8 1.0 fl oz 204 Total lipid (fat) g 0.10 0.3
14534 Alcoholic beverage, liqueur, coffee, 63 proof 34.8 1.0 fl oz 205 Carbohydrate, by difference g 11.21 32.2
14015 Alcoholic beverage, pina colada, canned 32.6 1.0 fl oz 208 Energy kcal 77.00 237.0
14015 Alcoholic beverage, pina colada, canned 32.6 1.0 fl oz 269 Sugars, total g NA NA
14015 Alcoholic beverage, pina colada, canned 32.6 1.0 fl oz 204 Total lipid (fat) g 2.48 7.6
14015 Alcoholic beverage, pina colada, canned 32.6 1.0 fl oz 205 Carbohydrate, by difference g 9.00 27.6
14019 Alcoholic beverage, tequila sunrise, canned 31.1 1.0 fl oz 208 Energy kcal 34.00 110.0
14019 Alcoholic beverage, tequila sunrise, canned 31.1 1.0 fl oz 269 Sugars, total g NA NA
14019 Alcoholic beverage, tequila sunrise, canned 31.1 1.0 fl oz 204 Total lipid (fat) g 0.03 0.1
14019 Alcoholic beverage, tequila sunrise, canned 31.1 1.0 fl oz 205 Carbohydrate, by difference g 3.51 11.3

Prep Time: 20mins

Great, we've successfully unnested. As I mentioned before, we'll use our nice ABBREV.xlsx rather than using data pulled from the API. So:

abbrev_raw <- readxl::read_excel("./data/raw/ABBREV.xlsx") %>% as_tibble()
## Warning in read_fun(path = path, sheet = sheet, limits = limits, shim =
## shim, : partial argument match of 'sheet' to 'sheet_i'
abbrev_raw %>% sample_n(20) %>% kable(format = "html")
NDB_No Shrt_Desc Water_(g) Energ_Kcal Protein_(g) Lipid_Tot_(g) Ash_(g) Carbohydrt_(g) Fiber_TD_(g) Sugar_Tot_(g) Calcium_(mg) Iron_(mg) Magnesium_(mg) Phosphorus_(mg) Potassium_(mg) Sodium_(mg) Zinc_(mg) Copper_mg) Manganese_(mg) Selenium_(µg) Vit_C_(mg) Thiamin_(mg) Riboflavin_(mg) Niacin_(mg) Panto_Acid_mg) Vit_B6_(mg) Folate_Tot_(µg) Folic_Acid_(µg) Food_Folate_(µg) Folate_DFE_(µg) Choline_Tot_ (mg) Vit_B12_(µg) Vit_A_IU Vit_A_RAE Retinol_(µg) Alpha_Carot_(µg) Beta_Carot_(µg) Beta_Crypt_(µg) Lycopene_(µg) Lut+Zea_ (µg) Vit_E_(mg) Vit_D_µg Vit_D_IU Vit_K_(µg) FA_Sat_(g) FA_Mono_(g) FA_Poly_(g) Cholestrl_(mg) GmWt_1 GmWt_Desc1 GmWt_2 GmWt_Desc2 Refuse_Pct
10998 CANADIAN BACON,CKD,PAN-FRIED 62.50 146 28.31 2.78 4.60 1.80 0.0 1.20 7 0.56 27 309 999 993 1.73 0.063 0.016 50.4 0.0 0.669 0.185 9.988 0.720 0.280 4 0 4 4 104.8 0.43 0 0 0 0 0 0 0 0 0.41 0.2 9 0.2 1.039 1.255 0.485 67 13.80 1 slice NA NA 0
21382 MCDONALD'S,FILET-O-FISH (WITHOUT TARTAR SAUCE) 46.90 243 12.47 7.62 1.93 31.08 1.0 4.27 128 1.64 23 131 194 464 0.54 0.065 0.230 NA 0.1 0.283 0.201 2.748 NA NA 57 NA NA NA NA 0.81 93 NA NA NA NA NA NA NA NA NA NA NA 1.773 2.768 2.323 25 124.00 1 item NA NA 0
16596 MORNINGSTAR FARMS GRILLERS QRTR PND VEGGIE BRGR,FRZ,UNPRP 55.60 219 22.80 10.50 2.20 8.90 2.5 0.60 79 4.90 NA 124 235 429 1.10 NA NA NA 28.0 1.150 0.300 10.800 NA 0.780 2 NA 2 NA NA 7.70 1 NA NA NA NA NA NA NA NA NA 0 NA 1.900 2.000 4.900 1 114.00 1 burger NA NA 0
21018 FAST FOODS,EGG,SCRAMBLED 66.70 212 13.84 16.18 1.20 2.08 0.0 1.64 57 2.59 14 242 147 187 1.66 0.067 0.043 22.5 3.3 0.080 0.520 0.210 0.940 0.190 29 0 29 29 180.6 1.01 679 176 171 0 62 6 0 233 0.96 1.1 46 9.0 6.153 5.889 1.969 426 96.00 2 eggs NA NA 0
01070 DESSERT TOPPING,POWDERED 1.47 577 4.90 39.92 1.17 52.54 0.0 52.54 17 0.03 7 74 166 122 0.08 0.118 0.225 0.6 0.0 0.000 0.000 0.000 0.000 0.000 0 0 0 0 0.1 0.00 0 0 0 0 0 0 0 0 1.52 0.0 0 9.9 36.723 0.600 0.447 0 43.00 1.5 oz 1.3 1 portion, amount to make 1 tbsp 0
03009 BABYFOOD,MEAT,HAM,JUNIOR 80.50 97 11.30 3.80 0.70 3.70 0.0 0.00 5 1.01 11 89 210 44 1.70 0.068 NA 15.0 2.1 0.142 0.194 2.840 0.531 0.200 2 0 2 2 45.2 0.10 0 0 0 0 0 0 0 0 0.40 0.4 18 0.0 1.271 1.804 0.516 29 28.35 1 oz 71.0 1 jar 0
23584 BEEF,TOP SIRLOIN,STEAK,LN,1/8" FAT,SEL,RAW 73.31 127 22.27 3.54 1.19 0.00 0.0 0.00 22 1.61 23 211 357 56 4.00 0.077 0.011 30.8 0.0 0.075 0.120 6.469 0.654 0.628 13 0 13 13 93.0 0.94 0 0 0 0 0 0 0 0 0.28 0.0 2 1.1 1.307 1.422 0.167 59 28.35 1 oz 453.6 1 lb 14
06740 SOUP,CHICK VEG,CHUNKY,RED FAT,RED NA,RTS,SINGLE BRAND 89.50 40 2.70 0.50 1.00 6.30 NA NA NA NA NA NA NA 192 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 1283 NA NA NA 770 NA NA NA NA NA NA NA 0.121 0.168 0.107 4 240.00 1 serving 454.0 1 package, yields 0
21234 MCDONALD'S,QUARTER POUNDER 50.37 244 14.10 11.55 1.81 22.17 1.6 5.13 84 2.41 22 124 227 427 2.69 0.107 0.199 NA 0.9 0.183 0.344 4.452 NA NA 56 NA NA NA NA 1.28 56 NA NA NA NA NA NA NA NA NA NA NA 4.008 4.202 0.283 39 171.00 1 item NA NA 0
11296 ONION RINGS,BREADED,PAR FR,FRZ,PREP,HTD IN OVEN 46.37 276 4.14 14.30 1.40 33.79 2.2 5.10 31 1.25 17 71 123 370 0.42 0.073 0.348 5.6 1.6 0.185 0.116 1.349 0.294 0.117 33 0 33 33 10.7 0.00 0 0 0 0 0 0 0 0 0.46 0.0 0 34.1 2.137 3.000 7.633 0 48.00 1 cup 71.0 10 rings, large (3-4" dia) 0
23449 BEEF,NZ,IMP,BRISKET NAVEL END,LN & FAT,RAW 53.33 345 15.81 31.27 0.61 0.00 0.0 0.00 9 1.11 13 118 223 54 2.66 0.037 0.003 2.5 0.0 0.056 0.076 2.683 0.255 0.177 NA NA NA NA NA 1.38 29 9 9 0 0 0 0 0 0.29 0.2 8 NA 12.915 10.892 0.766 62 114.00 4 oz NA NA 0
10956 PORK,LOIN,LEG CAP STEAK,BNLESS,LN & FAT,CKD,BRLD 68.72 158 27.57 4.41 1.07 0.00 0.0 0.00 7 0.97 23 221 366 76 4.11 0.073 0.010 32.1 0.0 0.481 0.387 8.205 0.772 0.430 1 0 1 1 88.4 0.69 0 0 0 0 0 0 0 0 0.10 NA NA 0.0 1.358 1.928 0.578 81 85.00 3 oz 194.0 1 piece 0
08575 CEREALS,CRM OF WHT,2 1/2 MIN COOK TIME,CKD W/H2O,MW,WO/ SALT 87.11 52 1.88 0.37 0.54 10.10 0.7 0.34 131 4.98 8 57 26 45 0.26 0.052 0.226 NA 0.0 0.087 0.045 0.715 0.312 0.033 32 20 12 46 NA NA NA NA NA NA NA NA NA NA 0.06 NA NA 0.0 0.099 0.058 0.176 NA 231.00 1 cup NA NA 0
20022 CORNMEAL,DEGERMED,ENR,YEL 11.18 370 7.11 1.75 0.51 79.45 3.9 1.61 3 4.36 32 99 142 7 0.66 0.076 0.174 10.5 0.0 0.551 0.382 4.968 0.240 0.182 209 180 30 335 8.6 0.00 214 11 0 63 97 0 0 1628 0.12 0.0 0 0.0 0.220 0.390 0.828 0 157.00 1 cup NA NA 0
10019 PORK,FRSH,LEG (HAM),SHANK HALF,LN,CKD,RSTD 65.28 175 28.69 5.83 1.20 0.00 0.0 0.00 14 0.92 24 261 376 84 2.67 0.123 0.019 28.6 0.0 0.440 0.373 8.082 0.877 0.478 0 0 0 0 98.6 0.50 4 1 1 0 0 0 0 0 0.26 0.4 15 0.0 1.837 2.482 1.130 93 85.00 3 oz 2900.0 1 roast 27
23321 BEEF,AUS,IMP,WGU,SMLEDRIB STK/RST,BNLES,LN&FAT,MRBSCR4/5,RAW 54.63 317 17.07 27.64 0.78 0.00 0.0 0.00 4 1.68 NA NA NA 52 NA NA NA NA 0.0 NA NA NA NA NA NA NA NA NA NA NA 17 5 5 0 0 0 0 0 0.00 NA NA NA 10.329 13.685 0.763 76 114.00 4 oz 342.0 1 roast 4
12061 ALMONDS 4.41 579 21.15 49.93 2.97 21.55 12.5 4.35 269 3.71 270 481 733 1 3.12 1.031 2.179 4.1 0.0 0.205 1.138 3.618 0.471 0.137 44 0 44 44 52.1 0.00 2 0 0 0 1 0 0 1 25.63 0.0 0 0.0 3.802 31.551 12.329 0 143.00 1 cup, whole 92.0 1 cup, sliced 60
10868 PORK,CURED,HAM -- H2O ADDED,SLICE,BONE-IN,LN,HTD,PAN-BROIL 68.58 131 22.04 4.30 4.06 1.48 0.0 1.48 11 1.12 20 260 289 1374 2.29 0.118 0.022 30.0 0.0 0.384 0.182 5.436 0.726 0.461 2 0 2 2 90.2 0.58 38 12 12 0 0 0 0 0 0.26 NA NA 0.0 1.439 1.989 0.687 65 85.00 1 serving, (3 oz) 436.0 1 slice 13
17398 LAMB,NZ,IMP,LOIN CHOP,LN,CKD,FAST FRIED 60.21 208 27.43 10.70 1.25 0.41 0.0 0.00 42 1.88 27 233 367 84 3.48 0.152 0.009 7.3 0.0 0.097 0.222 6.651 0.806 0.189 NA NA NA NA NA 1.84 20 6 6 0 0 0 0 0 0.37 0.1 2 NA 4.039 2.745 0.539 85 85.00 3 oz NA NA 44
17413 LAMB,NZ,IMP,NETTED SHLDR,ROLLED,BNLESS,L & F,CKD,SLOW RSTD 56.45 287 21.45 22.33 0.96 0.05 0.0 0.00 7 1.33 21 178 305 61 4.28 0.101 0.009 4.2 0.0 0.116 0.167 3.739 0.673 0.107 NA NA NA NA NA 1.94 46 14 14 0 0 0 0 0 0.47 0.1 3 NA 8.431 5.677 0.851 78 85.00 3 oz NA NA 2

How much data are we working with here?

dim(abbrev_raw)
## [1] 8790   53

You can read in depth the prep I did on this file in /scripts/prep. Mainly this involved a bit of cleaning like stripping out parentheses from column names, e.g., Vit_C_(mg) becomes Vit_C_mg.

In there you'll also find a dataframe called all_nut_and_mr_df where I define the nutritional constraints on menus. If a nutrient is among the "must restricts," that is, it's one of Lipid_Tot_g, Sodium_mg, Cholestrl_mg, FA_Sat_g, then its corresponding value is a daily upper bound. Otherwise, the nutrient is a "positive nutrient" and its value is a lower bound. For example, you're supposed to have at least 18mg of Iron and no more than 2400mg of Sodium per day. (As someone who puts salt on everything indiscriminately I'd be shocked if I've ever been under that threshold.)

all_nut_and_mr_df %>% kable(format = "html")
nutrient value
Lipid_Tot_g 65
Sodium_mg 2400
Cholestrl_mg 300
FA_Sat_g 20
Protein_g 56
Calcium_mg 1000
Iron_mg 18
Magnesium_mg 400
Phosphorus_mg 1000
Potassium_mg 3500
Zinc_mg 15
Copper_mg 2
Manganese_mg 2
Selenium_µg 70
Vit_C_mg 60
Thiamin_mg 2
Riboflavin_mg 2
Niacin_mg 20
Panto_Acid_mg 10
Vit_B6_mg 2

In /scripts/prep we also create a z-scored version of abbrev with:

scaled <- abbrev %>% 
  drop_na_(all_nut_and_mr_df$nutrient) %>% filter(!(is.na(Energ_Kcal)) & !(is.na(GmWt_1))) %>% 
  mutate_at(
    vars(nutrient_names, "Energ_Kcal"), dobtools::z_score   # <-- equivalent to scale(), but simpler
  )

I usually shunt Shrt_Desc all the way off screen to the right but here I'll put it next to its shorter sibling, shorter_desc so you can see how the truncation of name looks. Here's a random sample of 20 foods.

scaled %>% sample_n(20) %>% 
  select(shorter_desc, Shrt_Desc, everything()) %>% kable(format = "html")
shorter_desc Shrt_Desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal NDB_No score scaled_score
TUNA TUNA,LT,CND IN OIL,WO/SALT,DRND SOL 1 85.00 85.00 1.49 -0.1329359 -0.2243601 -0.2369700 -0.3085528 1.4434110 -0.2594454 -0.1957240 -0.0542244 0.6145169 -0.2197257 -0.3739253 -0.1985987 -0.0803815 1.8350398 -0.1318713 -0.3130598 -0.2579559 1.8612349 -0.2079963 -0.3881159 -0.0955685 15183 -2859.920 0.6169279
BROCCOLI BROCCOLI,RAW 1 91.00 91.00 3.90 -0.6926310 -0.2394511 -0.3613150 -0.5482143 -0.9547258 -0.0979257 -0.3496006 -0.2454015 -0.4803418 0.0496502 -0.5073419 -0.2319270 -0.0560104 -0.4358511 1.5091280 -0.2446894 -0.2644760 -0.6419115 -0.0610370 -0.2347818 -1.1698735 11090 -2927.355 0.4880110
EGG CUSTARDS EGG CUSTARDS,DRY MIX,PREP W/ WHL MILK 1 141.00 141.00 3.15 -0.4334865 -0.1941781 -0.0090043 -0.2287191 -0.8480812 0.3391276 -0.4405278 -0.3409900 -0.1943379 -0.2197257 -0.4801140 -0.2607106 -0.0813813 -0.3277134 -0.1300316 -0.2654077 -0.0406215 -0.7498186 0.0301791 -0.4966292 -0.5934172 19170 -2861.999 0.6129527
CRAYFISH CRAYFISH,MXD SP,WILD,RAW 1 85.00 85.00 7.39 -0.6512250 -0.2172585 0.4262031 -0.5289773 0.2438868 -0.1929373 -0.3239545 -0.1306952 0.3687323 0.0150514 -0.2650138 0.3285947 -0.0540107 0.4632363 -0.1097951 -0.2467613 -0.4492102 -0.3079742 -0.0805833 -0.3928339 -0.8881960 15145 -2954.602 0.4359224
GARLIC BREAD GARLIC BREAD,FRZ 1 43.00 43.00 7.23 0.4667376 0.2141670 -0.3613150 0.2921239 -0.4497589 -0.1929373 0.1912992 -0.2071661 -0.3864967 -0.4767449 -0.3820936 -0.1425466 -0.0271400 -0.0002108 -0.1281919 0.5653959 -0.1188619 0.1036477 0.0019455 -0.4518085 0.9001288 18963 -3499.014 -0.6048485
CHEESE CHEESE,ROQUEFORT 1 28.35 28.35 5.23 1.4683349 1.3371159 0.2604098 2.5335603 0.7515882 2.8236809 -0.3892355 -0.0733421 0.9764906 -0.5064010 -0.0526363 -0.2546509 -0.0785068 -0.0650934 -0.1318713 -0.3089162 0.7548225 -0.6216922 0.7772824 -0.3550901 1.0245909 01039 -3581.506 -0.7625507
ANISE SEED ANISE SEED 1 2.10 2.10 6.73 0.4160509 -0.2545421 -0.3613150 -0.4605254 0.3924601 2.7476716 8.0972960 2.6031363 1.1909935 2.8299059 0.8241013 1.0724222 0.2051977 -0.3586099 0.2544626 0.3126327 0.1115126 -0.1266392 0.1011249 0.8857366 0.8149705 02002 -3316.067 -0.2551022
SCHOOL LUNCH SCHOOL LUNCH,PIZZA,SAUSAGE TOP,THICK CRUST,WHL GRN,FRZ,CKD 1 129.00 129.00 9.06 -0.0715407 0.1156315 -0.1817056 0.1218761 -0.0031276 0.4768945 -0.0045439 -0.2071661 0.4715149 -0.0343753 -0.1833301 -0.1198227 -0.0353887 0.2593196 -0.1245126 0.2836271 0.1571528 0.0947086 -0.1602164 -0.3598080 0.2909192 21605 -2950.914 0.4429727
PICKLE&PIMIENTO LOAF PICKLE&PIMIENTO LOAF,PORK 1 38.00 38.00 9.32 0.4196204 0.6544695 0.0393521 0.2930858 -0.1881606 0.1966102 -0.2097128 0.0031287 -0.0915552 0.1855738 -0.1615478 -0.1470913 -0.0683834 -0.2690101 0.0116241 0.4203679 -0.2688226 -0.2488062 -0.1551489 0.3384519 0.0812988 07058 -3532.926 -0.6696790
CANDIES CANDIES,CAROB,UNSWTND 1 28.35 28.35 2.25 1.5197355 -0.1737608 -0.3544070 4.0973718 -0.4689003 1.1182227 -0.2190386 0.0413641 -0.2122131 0.8330645 0.3421679 -0.0289273 -0.0647590 -0.3524306 -0.1226729 -0.1846064 -0.1319019 -0.5565649 0.0670999 -0.3409362 2.1447504 19071 -3104.445 0.1494628
ONIONS ONIONS,CND,SOL&LIQUIDS 1 63.00 63.00 7.65 -0.7126202 0.0605937 -0.3613150 -0.5519014 -1.1342898 -0.1074269 -0.4894885 -0.5321670 -0.6501566 -0.4569742 -0.5400153 -0.2228375 -0.0695082 -0.5038233 -0.0527648 -0.3254908 -0.5057172 -0.7649298 -0.4056312 -0.3244233 -1.2681331 11285 -3484.090 -0.5763172
PIZZA HUT 12" PEPPERONI PIZZA PIZZA HUT 12" PEPPERONI PIZZA,HAND-TOSSED CRUST 1 96.00 96.00 3.96 0.0933695 0.4423076 -0.1817056 0.2728869 -0.0395873 0.4151369 -0.0208641 -0.2071661 0.1989175 -0.2197257 -0.1615478 -0.1486062 -0.0366385 0.1450027 -0.1318713 0.2711961 0.0702191 0.0793846 -0.1377744 -0.2890385 0.4415840 21274 -3562.980 -0.7271334
MUNG BNS MUNG BNS,MATURE SEEDS,SPROUTED,RAW 1 104.00 104.00 1.67 -0.7061951 -0.2634192 -0.3613150 -0.5470922 -0.9346729 -0.2594454 -0.3076343 -0.2454015 -0.5339675 -0.3630633 -0.5073419 -0.0577108 -0.0587599 -0.4945544 0.1109671 -0.2177557 -0.2492625 -0.6184997 -0.2007569 -0.4400136 -1.1960761 11043 -3113.261 0.1326089
GRAPES GRAPES,CND,THOMPSON SEEDLESS,H2O PK,SOL&LIQUIDS 1 245.00 245.00 7.22 -0.7111924 -0.2634192 -0.3613150 -0.5488555 -1.1661920 -0.2736972 -0.2913140 -0.5321670 -0.6948447 -0.4668595 -0.6053622 -0.2213226 -0.0773819 -0.5100026 -0.1134745 -0.3275626 -0.4687703 -0.7502443 -0.4461717 -0.4942703 -1.1305697 09133 -3036.218 0.2798926
PANCAKES PANCAKES,PLN,DRY MIX,INCOMPLETE (INCL BTTRMLK) 1 28.35 28.35 8.24 -0.5976828 1.0592637 -0.3613150 -0.5143892 -0.3002742 1.3034953 0.1912992 -0.0542244 2.0266611 -0.2592671 -0.3902620 -0.1319421 -0.0328891 -0.1145277 -0.1318713 0.8720267 0.2875535 0.0219193 -0.1674558 -0.2159099 0.9328820 18291 -3451.765 -0.5145198
PORK PORK,FRSH,LOIN,SIRLOIN (CHOPS OR ROASTS),BNLESS,LN&FAT,RAW 1 85.00 85.00 4.88 -0.4299170 -0.2128199 0.0738923 -0.3311563 0.8381800 -0.2784477 -0.3915670 -0.2262838 0.3285130 0.1435610 -0.1561023 -0.2016286 -0.0808814 0.6208083 -0.1318713 0.8513084 0.0289255 0.7627961 0.0461057 0.7677873 -0.5213601 10210 -2888.234 0.5627985
FAST FOODS FAST FOODS,CHSEBURGER; DBLE,REG,PATTY & BN; W/ CONDMNT & VEG 1 228.00 228.00 5.10 0.3853533 0.0898880 -0.0780848 0.3432624 -0.0231805 0.0303400 -0.0371844 -0.3409900 -0.0915552 -0.3086939 -0.1261516 -0.1985987 -0.0672586 0.0214168 -0.1097951 0.1261680 -0.1058218 0.0010615 -0.2731507 -0.3645260 0.4743372 21095 -3401.267 -0.4179810
RUTABAGAS RUTABAGAS,CKD,BLD,DRND,W/SALT 1 120.00 120.00 7.44 -0.7061951 -0.0432680 -0.3613150 -0.5498174 -1.1269979 -0.2356925 -0.4778312 -0.4556962 -0.5920620 -0.1974837 -0.5863027 -0.2622255 -0.0701331 -0.4914647 0.2139895 -0.2218993 -0.4296501 -0.6257361 -0.3636428 -0.4069878 -1.1960761 11851 -3310.710 -0.2448605
PORK PORK,SHLDR PETITE TENDER,BNLESS,LN & FAT,RAW 1 105.00 105.00 2.62 -0.4399115 -0.2243601 0.0946165 -0.3922339 0.7616146 -0.2926995 -0.2703308 -0.1880484 0.2480744 0.2350006 -0.0308540 -0.1698152 -0.0808814 0.1511819 -0.1318713 1.3299010 0.7852494 0.3526640 0.2350534 0.8857366 -0.5541133 10961 -2760.869 0.8062862
BEANS BEANS,KIDNEY,RED,MATURE SEEDS,CND,SOL & LIQ,LO NA 1 256.00 256.00 7.85 -0.6933449 -0.1648838 -0.3613150 -0.5344278 -0.7359676 -0.1834361 -0.2283645 -0.0733421 -0.3015893 -0.0887448 -0.4501633 -0.0834645 -0.0458870 -0.4791061 -0.1171538 -0.1721754 -0.3731432 -0.6727725 -0.3817413 -0.4588855 -0.8619934 16337 -2560.381 1.1895663

Then we do a few mutates to abbrev using the function below. This is a function we can use on any menu dataframe, not just abbrev, which is why it's called do_menu_mutates(). Turns out that the short descriptions of foods in the Shrt_Desc column actually aren't so short so we'll create a shorter_desc column by taking only the values in Shrt_Desc up to the first comma. That turns "BUTTER,WHIPPED,W/ SALT" into just "BUTTER".

Since we'll need a cost associated with each row in order to optimize something, for now each item gets a random cost between $1 and $10.

What we'll do when we eventually "solve" these menus is change the amount we have of each item, i.e. its GmWt_1. We'll vary that by multiplying it by some solution_amount. In order to keep a record of what the gram weight of a single serving of a food is, we'll save that in serving_gmwt. Since we know that all foods in abbrev are exactly one serving, for now GmWt_1 and serving_gmwt are the same thing, and solution_amounts is 1.

Finally, we rearrange columns a bit.

cols_to_keep <- c(nutrient_names, "Shrt_Desc", "GmWt_1", "NDB_No")

do_menu_mutates <- function(menu, to_keep = cols_to_keep) {

  quo_to_keep <- quo(to_keep)
  
  menu <- menu %>% 
    mutate(
      shorter_desc = map_chr(Shrt_Desc, grab_first_word, splitter = ","), # Take only the fist word
      cost = runif(nrow(.), min = 1, max = 10) %>% round(digits = 2) # Add a cost column
    ) 
  
  if (!("serving_gmwt" %in% names(menu))) {
    menu <- menu %>% mutate(
      serving_gmwt = GmWt_1   # Single serving gram weight
    )
  }
  
  if (!("solution_amounts" %in% names(menu))) {
    menu <- menu %>% mutate(
      solution_amounts = 1   # Single serving gram weight
    )
  }
  
  menu <- menu %>%
    select(shorter_desc, solution_amounts, GmWt_1, serving_gmwt, cost, !!quo_to_keep,  Shrt_Desc, NDB_No)
  
  return(menu)
}

We'll do these mutates and score each item (see /scripts/score/rank_foods.R for add_ranked_foods(); also more on scoring in the Scoring section below).

abbrev <- abbrev %>% do_menu_mutates() %>% add_ranked_foods() 
## score column doesn't exist; creating it
## scaled_score doesn't exist; creating it

And now we've got our main bucket of foods to pull from, each with a single portion size for now.

abbrev[1:20, ] %>% kable(format = "html")
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
BUTTER 1 5.00 5.00 6.26 81.11 643 215 51.368 0.85 24 0.02 2 24 24 0.09 0.000 0.000 1.0 0.0 0.005 0.034 0.042 0.110 0.003 717 BUTTER,WITH SALT 01001 -3419.716 -0.4532518
BUTTER 1 3.80 3.80 3.88 78.30 583 225 45.390 0.49 23 0.05 1 24 41 0.05 0.010 0.001 0.0 0.0 0.007 0.064 0.022 0.097 0.008 718 BUTTER,WHIPPED,W/ SALT 01002 -3405.992 -0.4270146
BUTTER OIL 1 12.80 12.80 8.05 99.48 2 256 61.924 0.28 4 0.00 0 3 5 0.01 0.001 0.000 0.0 0.0 0.001 0.005 0.003 0.010 0.001 876 BUTTER OIL,ANHYDROUS 01003 -3426.108 -0.4654711
CHEESE 1 28.35 28.35 6.20 28.74 1146 75 18.669 21.40 528 0.31 23 387 256 2.66 0.040 0.009 14.5 0.0 0.029 0.382 1.016 1.729 0.166 353 CHEESE,BLUE 01004 -3383.120 -0.3832890
CHEESE 1 132.00 132.00 7.83 29.68 560 94 18.764 23.24 674 0.43 24 451 136 2.60 0.024 0.012 14.5 0.0 0.014 0.351 0.118 0.288 0.065 371 CHEESE,BRICK 01005 -2550.059 1.2092995
CHEESE 1 28.35 28.35 1.06 27.68 629 100 17.410 20.75 184 0.50 20 188 152 2.38 0.019 0.034 14.5 0.0 0.070 0.520 0.380 0.690 0.235 334 CHEESE,BRIE 01006 -3427.868 -0.4688367
CHEESE 1 28.35 28.35 1.33 24.26 842 72 15.259 19.80 388 0.33 20 347 187 2.38 0.021 0.038 14.5 0.0 0.028 0.488 0.630 1.364 0.227 300 CHEESE,CAMEMBERT 01007 -3365.981 -0.3505239
CHEESE 1 28.35 28.35 8.82 29.20 690 93 18.584 25.18 673 0.64 22 490 93 2.94 0.024 0.021 14.5 0.0 0.031 0.450 0.180 0.190 0.074 376 CHEESE,CARAWAY 01008 -3234.675 -0.0995030
CHEESE 1 132.00 132.00 2.61 33.31 653 99 18.867 22.87 710 0.14 27 455 76 3.64 0.030 0.027 28.5 0.0 0.029 0.428 0.059 0.410 0.066 404 CHEESE,CHEDDAR 01009 -2687.571 0.9464129
CHEESE 1 28.35 28.35 7.93 30.60 700 103 19.475 23.37 643 0.21 21 464 95 2.79 0.042 0.012 14.5 0.0 0.046 0.293 0.080 0.413 0.074 387 CHEESE,CHESHIRE 01010 -3257.267 -0.1426935
CHEESE 1 132.00 132.00 4.12 32.11 604 95 20.218 23.76 685 0.76 26 457 127 3.07 0.042 0.012 14.5 0.0 0.015 0.375 0.093 0.210 0.079 394 CHEESE,COLBY 01011 -2599.704 1.1143912
CHEESE 1 113.00 113.00 8.72 4.30 364 17 1.718 11.12 83 0.07 8 159 104 0.40 0.029 0.002 9.7 0.0 0.027 0.163 0.099 0.557 0.046 98 CHEESE,COTTAGE,CRMD,LRG OR SML CURD 01012 -3386.210 -0.3891963
CHEESE 1 113.00 113.00 3.85 3.85 344 13 2.311 10.69 53 0.16 7 113 90 0.33 0.040 0.003 7.7 1.4 0.033 0.142 0.150 0.181 0.068 97 CHEESE,COTTAGE,CRMD,W/FRUIT 01013 -3463.568 -0.5370853
CHEESE 1 145.00 145.00 3.07 0.29 372 7 0.169 10.34 86 0.15 11 190 137 0.47 0.030 0.022 9.4 0.0 0.023 0.226 0.144 0.446 0.016 72 CHEESE,COTTAGE,NONFAT,UNCRMD,DRY,LRG OR SML CURD 01014 -3278.578 -0.1834343
CHEESE 1 113.00 113.00 7.66 2.27 308 12 1.235 10.45 111 0.13 9 150 125 0.51 0.033 0.015 11.9 0.0 0.020 0.251 0.103 0.524 0.057 81 CHEESE,COTTAGE,LOWFAT,2% MILKFAT 01015 -3266.099 -0.1595762
CHEESE 1 113.00 113.00 3.56 1.02 406 4 0.645 12.39 61 0.14 5 134 86 0.38 0.028 0.003 9.0 0.0 0.021 0.165 0.128 0.215 0.068 72 CHEESE,COTTAGE,LOWFAT,1% MILKFAT 01016 -3490.534 -0.5886355
CHEESE 1 14.50 14.50 8.59 34.44 314 101 20.213 6.15 97 0.11 9 107 132 0.50 0.018 0.011 8.6 0.0 0.023 0.230 0.091 0.517 0.056 350 CHEESE,CREAM 01017 -3389.710 -0.3958887
CHEESE 1 28.35 28.35 4.73 27.80 812 89 17.572 24.99 731 0.44 30 536 188 3.75 0.036 0.011 14.5 0.0 0.037 0.389 0.082 0.281 0.076 357 CHEESE,EDAM 01018 -3208.657 -0.0497637
CHEESE 1 150.00 150.00 3.77 21.28 917 89 14.946 14.21 493 0.65 19 337 62 2.88 0.032 0.028 15.0 0.0 0.154 0.844 0.991 0.967 0.424 264 CHEESE,FETA 01019 -3516.569 -0.6384083
CHEESE 1 132.00 132.00 6.64 31.14 800 116 19.196 25.60 550 0.23 14 346 64 3.50 0.025 0.014 14.5 0.0 0.021 0.204 0.150 0.429 0.083 389 CHEESE,FONTINA 01020 -3304.806 -0.2335737

Per 100g vs. Raw

Note how our column titles have always end with _g or _mg. That's because this column is giving us the value of each nutrient, per 100g of this food. The value we've got in that column isn't the raw value. Our constraints, though, are in raw terms. We'll need a way to know whether we've gotten our 1000mg of Calcium from the foods in our menu, each of which list how much Calcium they provide per 100g of that food.

In order to get to the raw value of a nutrient, for each food in our menu we'll multiply the 100g value of that nutrient by the weight of the food in grams, or its GmWt_1:

$TotalNutrientVal = \sum_{i=1}^{k} Per100gVal_{i} * GmWt_{i}$

where k the total number of foods in our menu.

Two helper functions get_per_g_vals() and get_raw_vals() in /scripts/solve allow us to go back and forth between raw and per 100g values. We'll try to keep everything in per 100g whenever possible, as that's the format our raw data is in. Our main solving function does accept both formats, however.

abbrev %>% sample_n(10) %>% get_raw_vals() %>% kable(format = "html")
shorter_desc GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal solution_amounts Shrt_Desc NDB_No score scaled_score
SAUCE 30.00 30.00 9.83 5.01000 200.1000 2.1000 0.9999000 0.300000 7.8000 0.07500 1.8000 5.1000 20.400 0.03600 0.0069000 0.03210 0.27000 0.69000 0.00480 0.008700 0.02790 0.0219000 0.013200 63.300 1 SAUCE,TARTAR,RTS 27049 -3545.6234 -0.6939525
PASTA 117.00 117.00 9.78 2.00070 4.6800 0.0000 0.2843100 7.008300 15.2100 2.01240 63.1800 148.5900 112.320 1.56780 0.2632500 1.54557 42.47100 0.00000 0.18252 0.115830 3.65742 0.3135600 0.108810 174.330 1 PASTA,WHOLE-WHEAT,CKD 20125 -2982.4185 0.3827436
CHEESE 132.00 132.00 7.85 23.23200 811.8000 72.6000 14.9160000 32.604000 997.9200 0.68640 36.9600 654.7200 182.160 4.26360 0.0343200 0.01320 19.14000 0.00000 0.02508 0.423720 0.20592 0.6283200 0.096360 361.680 1 CHEESE,PROVOLONE,RED FAT 01208 -2366.6671 1.5598948
BABYFOOD 15.00 15.00 3.47 0.07800 2.1000 0.1500 0.0555000 0.165000 4.5000 0.02100 1.5000 4.2000 15.000 0.03900 0.0030000 NA 0.13500 2.08500 0.00150 0.006000 0.02850 NA 0.012000 11.700 1 BABYFOOD,DSSRT,BANANA YOGURT,STR 43539 -3348.6875 -0.3174641
MACKEREL 28.35 28.35 6.69 1.78605 107.4465 22.3965 0.5264595 6.574365 68.3235 0.57834 10.4895 85.3335 54.999 0.28917 0.0416745 0.01134 10.68795 0.25515 0.01134 0.060102 1.75203 0.0864675 0.059535 44.226 1 MACKEREL,JACK,CND,DRND SOL 15048 -3266.6025 -0.1605397
NUTS 157.00 157.00 4.44 86.61690 224.5100 0.0000 6.6049900 33.331100 456.8700 5.77760 430.1800 731.6200 1097.430 4.81990 1.4993500 3.86220 6.43700 0.00000 0.14444 1.226170 5.75405 0.3595300 0.185260 952.990 1 NUTS,ALMONDS,OIL RSTD,LIGHTLY SALTED 12665 -944.2775 4.2791211
PORK 85.00 85.00 2.25 11.57700 40.8000 68.0000 4.3435000 23.145500 17.8500 0.90950 16.1500 153.8500 317.900 2.02300 0.0654500 0.01020 38.50500 0.51000 0.53720 0.215900 3.75615 0.5508000 0.311100 203.150 1 PORK,FRSH,LOIN,WHL,LN&FAT,CKD,BRSD 10021 -2922.4307 0.4974243
YOGURT 170.00 170.00 5.29 0.66300 61.2000 8.5000 0.1989000 17.323000 187.0000 0.11900 18.7000 229.5000 239.700 0.88400 0.0289000 0.01530 16.49000 0.00000 0.03910 0.472600 0.35360 0.5627000 0.107100 100.300 1 YOGURT,GREEK,PLN,NONFAT 01256 -2733.2666 0.8590551
RUFFED GROUSE 113.00 113.00 4.41 0.99440 56.5000 45.2000 0.1469000 29.312200 5.6500 0.65540 36.1600 258.7700 351.430 0.57630 0.0655400 0.01808 NA 0.00000 0.04746 0.316400 13.10800 NA 1.440750 126.560 1 RUFFED GROUSE,BREAST MEAT,SKINLESS,RAW 05363 -2779.2912 0.7710685
PATE 13.00 13.00 7.78 3.64000 90.6100 33.1500 1.2441000 1.846000 9.1000 0.71500 1.6900 26.0000 17.940 0.37050 0.0520000 0.01560 5.40800 0.26000 0.00390 0.078000 0.42900 0.1560000 0.007800 41.470 1 PATE,LIVER,NOT SPECIFIED,CND 07055 -3438.5723 -0.4892996

Creating and Solving Menus

Building

Now that we've got our data, on to building a menu. The only constraint we'll worry about for now is that menus have to contain at least 2300 calories. Our strategy is simple; pick one serving of a food at random from our dataset and, if it doesn't yet exist in our menu, add it. We do this until we're no longer under 2300 calories.

That's implemented in add_calories() below, which we'll as a helper inside the main building function, build_menu(). The reason I've spun add_calories() out into its own function is so that we can easily add more foods to existing menus. Unlike build_menu() which takes a dataframe of possible foods to choose from as its first argument, add_calories() takes menu as its first argument, making it convenient to pipe a menu in that needs more calories.

add_calories <- function(menu = NULL, df = abbrev, seed = NULL, ...) {

  # If we're starting from an existing menu
  if (! is.null(menu)) {
    menu <- menu %>% drop_na_(all_nut_and_mr_df$nutrient) %>%  filter(!(is.na(Energ_Kcal)) & !(is.na(GmWt_1)))
    cals <- sum((menu$Energ_Kcal * menu$GmWt_1), na.rm = TRUE)/100   # Set calories to our current number of calories
  # If we're starting from scratch
  } else {
    cals <- 0   
    menu <- NULL
  }

  while (cals < 2300) {
    df <- df %>% filter(!NDB_No %in% menu$NDB_No)   # Only add foods we don't already have

    if (nrow(df) == 0) {
      message("No more elligible foods to sample from. Returning menu too low in calories.")
      return(menu)
    } else {
      food_i <- df %>%
        sample_n(1)   # Sample a new index from a food that doesn't already exist in our menu
    }

    this_food_cal <- (food_i$Energ_Kcal * food_i$GmWt_1)/100   
    cals <- cals + this_food_cal    

    menu <- bind_rows(menu, food_i)   
  }
  return(menu)   
}

Okay now for build_menu(). We'll make sure we don't have missing values in any of our nutrient columns, calories, or the food weight. The from_better_cutoff argument allows us to specify that we only want to pull foods that have at least a certain z-score on our scaled_score dimension. More on scoring in a bit.

The default, though, will just be to pull foods from our main abbrev dataframe.

build_menu <- function(df = abbrev, menu = NULL, seed = NULL, from_better_cutoff = NULL, do_mutates = TRUE) {
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  df <- df %>% drop_na_(all_nut_and_mr_df$nutrient) %>%
    filter(!(is.na(Energ_Kcal)) & !(is.na(GmWt_1)))    # Filter out rows that have NAs in columns that we need
  
  # Optionally choose a floor for what the z-score of each food to build from should be
  if (!is.null(from_better_cutoff)) {
    assert_that(is.numeric(from_better_cutoff), msg = "from_better_cutoff must be numeric or NULL")
    if (! "scaled_score" %in% names(df)) {
      df <- df %>% 
        add_ranked_foods()
    }
    df <- df %>% 
      filter(scaled_score > from_better_cutoff)
  }
  
  if (nrow(df) == 0) {
    stop("No foods to speak of; you might try a lower cutoff.")
  }
  
  # Add one serving of food until we hit 2300
  menu <- add_calories(menu = menu, df = df)
  
  return(menu)
}
our_random_menu <- build_menu()
our_random_menu %>% kable(format = "html")
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
CANDIES 1 40.00 40.00 5.43 34.50 39 11 20.600 5.58 109 1.33 38 129 306 0.79 0.180 0.000 0.0 1.0 0.040 0.160 0.330 0.250 0.060 563 CANDIES,HERSHEY'S,ALMOND JOY BITES 19248 -3179.352 0.0062599
SALMON 1 85.00 85.00 3.75 7.31 75 44 1.644 20.47 239 1.06 29 326 377 1.02 0.084 0.030 35.4 0.0 0.016 0.193 5.480 0.550 0.300 153 SALMON,SOCKEYE,CND,WO/SALT,DRND SOL W/BONE 15182 -2602.498 1.1090489
BROADBEANS 1 109.00 109.00 1.47 0.60 50 0 0.138 5.60 22 1.90 38 95 250 0.58 0.074 0.320 1.2 33.0 0.170 0.110 1.500 0.086 0.038 72 BROADBEANS,IMMAT SEEDS,RAW 11088 -2939.264 0.4652428
GELATIN DSSRT 1 85.00 85.00 4.71 0.00 466 0 0.000 7.80 3 0.13 2 141 7 0.01 0.118 0.011 6.7 0.0 0.003 0.041 0.009 0.014 0.001 381 GELATIN DSSRT,DRY MIX 19172 -3627.439 -0.8503611
CRAYFISH 1 85.00 85.00 5.25 1.20 94 133 0.181 16.77 60 0.83 33 270 296 1.76 0.685 0.522 36.7 0.9 0.050 0.085 2.280 0.580 0.076 82 CRAYFISH,MXD SP,WILD,CKD,MOIST HEAT 15146 -2955.922 0.4333988
CRACKERS 1 30.00 30.00 3.35 11.67 1167 0 3.333 10.00 67 4.80 22 162 141 0.90 0.118 0.517 26.2 0.0 1.087 0.750 7.170 0.528 0.044 418 CRACKERS,CHS,RED FAT 18965 -3595.367 -0.7890484
CHEESE 1 28.35 28.35 7.10 30.64 1809 90 19.263 21.54 662 0.56 30 392 91 2.08 0.034 0.030 14.5 0.0 0.040 0.586 0.734 1.731 0.124 369 CHEESE,ROQUEFORT 01039 -3581.506 -0.7625507
SPICES 1 0.70 0.70 3.60 4.07 76 0 2.157 22.98 2240 89.80 711 274 2630 7.10 2.100 9.800 3.0 0.8 0.080 1.200 4.900 0.838 1.340 233 SPICES,BASIL,DRIED 02003 -3332.583 -0.2866766
RAVIOLI 1 242.00 242.00 9.97 1.45 306 3 0.723 2.48 33 0.74 15 50 232 0.36 0.142 0.176 3.5 0.0 0.074 0.080 1.060 0.272 0.102 77 RAVIOLI,CHEESE-FILLED,CND 22899 -3306.693 -0.2371810
LUXURY LOAF 1 28.00 28.00 5.14 4.80 1225 36 1.580 18.40 36 1.05 20 185 377 3.05 0.100 0.041 21.5 0.0 0.707 0.297 3.482 0.515 0.310 141 LUXURY LOAF,PORK 07060 -3541.980 -0.6869870
CANDIES 1 14.50 14.50 9.75 30.00 11 0 17.750 4.20 32 3.13 115 132 365 1.62 0.700 0.800 4.2 0.0 0.055 0.090 0.427 0.105 0.035 480 CANDIES,SEMISWEET CHOC 19080 -3286.911 -0.1993645
ONIONS 1 210.00 210.00 2.31 0.05 8 0 0.009 0.71 27 0.34 8 2 101 0.09 0.024 0.040 0.4 5.1 0.016 0.018 0.132 0.078 0.070 28 ONIONS,FRZ,WHL,CKD,BLD,DRND,WO/SALT 11290 -3086.386 0.1839856
PIZZA HUT 14" PEPPERONI PIZZA 1 113.00 113.00 6.17 13.07 676 23 4.823 11.47 147 2.57 22 193 187 1.36 0.104 0.425 15.5 1.0 0.420 0.210 3.750 0.323 0.090 291 PIZZA HUT 14" PEPPERONI PIZZA,PAN CRUST 21297 -3521.658 -0.6481376
BEANS 1 104.00 104.00 3.07 0.70 13 0 0.085 6.15 15 1.93 101 100 307 0.89 0.356 0.408 0.6 18.8 0.390 0.215 1.220 0.825 0.191 67 BEANS,NAVY,MATURE SEEDS,SPROUTED,RAW 11046 -2811.162 0.7101393
GRAPE JUC 1 253.00 253.00 8.25 0.13 5 0 0.025 0.37 11 0.25 10 14 104 0.07 0.018 0.239 0.0 0.1 0.017 0.015 0.133 0.048 0.032 60 GRAPE JUC,CND OR BTLD,UNSWTND,WO/ ADDED VIT C 09135 -3032.103 0.2877596
PIE 1 28.35 28.35 4.42 12.50 211 0 3.050 2.40 7 1.12 7 28 79 0.19 0.053 0.185 7.8 1.7 0.148 0.107 1.230 0.093 0.032 265 PIE,APPL,PREP FROM RECIPE 18302 -3399.654 -0.4148992
PORK 1 85.00 85.00 8.42 2.59 63 63 0.906 22.81 9 0.56 23 251 354 1.72 0.069 0.011 37.4 0.0 0.610 0.256 7.348 0.728 0.611 121 PORK,FRSH,LOIN,SIRLOIN (CHOPS OR ROASTS),BNLESS,LN,RAW 10214 -2881.317 0.5760225
FAST FOODS 1 226.00 226.00 9.28 11.75 350 54 4.654 15.17 45 2.59 22 139 252 2.51 0.097 0.110 11.3 0.5 0.160 0.170 3.350 0.240 0.240 239 FAST FOODS,HAMBURGER; DOUBLE,LRG PATTY; W/ CONDMNT & VEG 21114 -3206.685 -0.0459943

Alright nice -- we've got a random menu that's at least compliant on calories. Is it compliant on nutrients and must restricts?

Testing Compliance

A few quick functions for testing whether we're compliant on the other dimensions. Nothing fancy here; all we're doing is going through positives and must restricts and figuring out how much of a given nutrient we've got and comparing that to the requirement. If we're below the minimum on any positives, above the maximum on any must restricts, or we're below 2300 calories we're out of compliance. To make it easier to see where we're out of compliance, we'll return a dataframe of the nutrients we're uncompliant on.

For must restricts:

test_mr_compliance <- function(orig_menu, capitalize_colname = TRUE) {
  
  compliance_df <- list(must_restricts_uncompliant_on = vector(), 
                        `difference_(g)` = vector()) %>% as_tibble()
  
  for (m in seq_along(mr_df$must_restrict)) {    
    nut_to_restrict <- mr_df$must_restrict[m]    # Grab the name of the nutrient we're restricting
    to_restrict <- (sum(orig_menu[[nut_to_restrict]] * orig_menu$GmWt_1, na.rm = TRUE))/100   # Get the amount of that must restrict nutrient in our original menu
    
    if ((to_restrict - mr_df$value[m]) > 0.01) {    # Account for rounding error
      this_compliance <- list(must_restricts_uncompliant_on = nut_to_restrict,
                              `difference_(g)` = (to_restrict - mr_df$value[m]) %>% round(digits = 2)) %>% as_tibble()
      compliance_df <- bind_rows(compliance_df, this_compliance)
    }
  }
  if (capitalize_colname == TRUE) {
    compliance_df <- compliance_df %>% cap_df()
  }
  return(compliance_df)
}

Same idea for positives. Then to test whether we're compliant overall, we'll see whether we pass all of these tests. If not, we're not compliant.

test_all_compliance <- function(orig_menu) {
  combined_compliance <- "Undetermined"
  
  if (nrow(test_mr_compliance(orig_menu)) + nrow(test_pos_compliance(orig_menu)) > 0 |
      test_calories(orig_menu) == "Calories too low") {
    combined_compliance <- "Not Compliant"
    
  } else if (nrow(test_mr_compliance(orig_menu)) + nrow(test_pos_compliance(orig_menu)) == 0 &
             test_calories(orig_menu) == "Calorie compliant") {
    combined_compliance <- "Compliant"
    
  } else {
    combined_compliance <- "Undetermined"
  }
  return(combined_compliance)
}

Let's see where we are with our random menu.

our_random_menu %>% test_all_compliance
## [1] "Not Compliant"

😔 We've got some work to do!

Scoring

Now I want an objective and preferably single scalar metric by which to judge menus. We want a metric that takes into account the following things:

  • You don't get extra credit for going above the daily minimum on positive nutrients
    • This reflects the fact that your body can only absorb up to a certain amount of a vitamin
  • You do, however, keep getting penalized for going farther and farther above the minimum on must_restricts
    • There's no cap on how bad an increase in bad stuff will keep

For simplicity and because I'm not a doctor, we'll assume a linear relationship between increasing and decreasing nutrients and their effect on our score. Though they're really two different dimensions, I want to be able to combine a must restrict score with a positive score to get a single number out. The directionality of the scores will also have to be the same if we want to combine them; so in both cases, more positive scores mean worse.

Similar to how we tested compliance, I'll do is go through a given menu and multiply the nutrient value per 100g by GmWt_1, the amount of the food we have. That will give us the raw amount of this nutrient. Then I'll see how much that raw amount differs from the minimum or maximum daily value of that nutrient we're supposed to have and give it a score accordingly. Then I'll add it up.

First, the must restricts. For each must restrict, we find the difference between our maximum allowed value and the value of that must restrict and add those all up. (Perhaps percent above the maximum would be a better metric.)

$\sum_{i=1}^{k} MaxAllowedAmount_{i} - AmountWeHave_{i}$

where k is the total number of must restricts.

So the farther above our max we are on must restricts, the higher our score will be.

mr_score <- function(orig_menu) {
  total_mr_score <- 0
  
  for (m in seq_along(mr_df$must_restrict)) {    
    mr_considering <- mr_df$must_restrict[m]    
    val_mr_considering <- (sum((orig_menu[[mr_considering]] * orig_menu$GmWt_1), 
                                na.rm = TRUE))/100   
    
    mr_score <- mr_df$value[m] - val_mr_considering  # max amount it's supposed to be - amount it is

    total_mr_score <- total_mr_score + mr_score
  }
  return(total_mr_score)
}

Similar story for the positives: we'll take the difference between our minimum required amount and the amount of the nutrient we've got in our menu and multiply that by -1:

$\sum_{i=1}^{k} (-1) * (MaxAllowedAmount_{i} - AmountWeHave_{i})$

where k is the total number of positive nutrients in our constraints.

That means that if we've got less than the amount we need, this value will be negative; if we've got more than we need it'll be positive. Next, to make the best score 0 I'll turn everything greater than 0 into 0. This takes away the extra brownie points for going above and beyond. Same as with must restricts, lower scores mean "better."

pos_score <- function(orig_menu) {
  total_nut_score <- 0
  
  for (p in seq_along(pos_df$positive_nut)) {    
    nut_considering <- pos_df$positive_nut[p]    
    val_nut_considering <- (sum((orig_menu[[nut_considering]] * orig_menu$GmWt_1), 
                                na.rm = TRUE))/100   
    
    nut_score <- (-1)*(pos_df$value[p] - val_nut_considering)    # (-1)*(min amount it's supposed to be - amount it is here)
    
    if (nut_score > 0) {
      nut_score <- 0
    } else if (is.na(nut_score)) {
      message("Nutrient has no score")
      break
    }
    total_nut_score <- total_nut_score + nut_score
  }
  return(total_nut_score)
}

Last step is just a simple sum:

score_menu <- function(orig_menu) {
  healthiness_score <- pos_score(orig_menu) + mr_score(orig_menu)
  return(healthiness_score)
}

Let's see what our menu's score is.

our_random_menu %>% score_menu()
## [1] -2018.912

Solving

Getting a Solution

Solving our menus is the next step. We've got a fixes set of constraints and an objective function: to minimize cost.

Given these conditions, it makes sense to use a simple linear programming algorithm. The implementation we use for solving is the GNU linear program solver which has an R interface via the Rglpk package. (You'll need to brew install glpk or otherwise install the GLPK in order to be able to install the Rglpk package.)

The Rglpk_solve_LP() function is going to do the work for us. What solve_it() below will do is grab the elements of a menu that we need for this function, pass them to the solver in the format it needs, and return a solution that is a list of a few things we're interested in: the cost of our final menu, the original menu, and the multiplier on each food's portion size.

Kind of a lot going on in solve_it(), which I'll walk through below. If you're only interested in what we get out of it, feel free to skip this section 😁.

Into the bowels of solve_it()

What we first have to do is get the raw values of every nutrient, if our nutrients are in per 100g form. (If they're already in raw form, we're all set.) We know they're in raw form already if the df_is_per_100g flag is FALSE. Whichever form we get our menu data in, we'll transform it to the other form in order to return that in our list at the end.

Rglpk_solve_LP() needs something to optimize for, which for us will be df[["cost"]]. We'll tell it we want to minimize that.

Next we need to set up a series of constraint inequalities. On the left hand side of each inequality will be the raw values of each nutrient we've got in our menu. That will be followed by a directionality, either ">" if the value of that nutrient is a positive or a "<" if it is a must restrict. Last we'll supply the upper or lower bound for that particular nutrient, which we supply in bounds. If we're thinking about Riboflavin in our menu and we've got n items in our menu each with some amount of Riboflavin, that would look like:

$\sum_{i=1}^{n} OurRawRiboflavin_{i} > MinRequiredDailyRiboflavin$

Now to construct the constraint matrix which I'm cleverly calling constraint_matrix for all the nutritional constraints that need to be met. We'll make this by essentially transposing our menu's nutrient columns; whereas in a typical menu dataframe we have one row per food and one column per nutrient, we'll turn this into a matrix with one row per constraint and one column per food. (In practice we do this by taking our vector of nutrient constraint values, and, in the matrix call of construct_matrix(), creating byrow = TRUE matrix from them.) We can print out the constraint matrix by turning v_v_verbose on.

Cool, so now we can read a given row in this matrix pertaining to a certain nutrient left to right as adding up the value of that nutrient contained in all of our menu foods. That gives us the sum total of that nutrient in the menu.

What we'd need next if we keep reading from left to right is the directionality which the solver accepts as a vector of ">"s and "<"s. Farthest to the right, we'll need the min or max value of that nutrient, which we'll supply in the vector rhs for "right hand side" of the equation. We get rhs from the user-supplied nut_df, or dataframe of nutrients and their daily upper or lower bounds.

We'll specify minimum and serving sizes per food by creating bounds from min_food_amount and max_food_amount. This acts as the other half of the constraint on the solver; not only do we need a certain amount of each nutrient, we also need a certain amount of each food.

Finally, we can specify that we want only full serving sizes by setting only_full_servings to TRUE. If we do that, we'll tell the solver that the types must be integer, "I" rather than continuous, "C".

If we turn verbose on we'll know whether a solution could be found or not, and what the cost of the solved menu is.
Return

The native return value of the call to Rglpk_solve_LP() is a list including: vector of solution amounts, the cost, and the status (0 for solved, 1 for not solvable). solve_it() will take that list and append a few more things to it, so we can check that it's working correctly and pipe it into other things to distill menus and nutritional information out . We'll append the nutritional constraints we supplied in nut_df, the constraint matrix we constructed, and the nutrient values in our original menu in both raw and per 100g form.

solve_it <- function(df, nut_df = nutrient_df, df_is_per_100g = TRUE, only_full_servings = FALSE, 
                     min_food_amount = 1, max_food_amount = 100, 
                     verbose = TRUE, v_v_verbose = FALSE, maximize = FALSE) {
  
  # If our nutrient values are per 100g (i.e., straight from menu_builder)
  if (df_is_per_100g == TRUE) {
    df_per_100g <- df        # Save our original df in df_per_100g
    df <- get_raw_vals(df)   # Get the raw values
  } else {
    df_per_100g <- get_per_g_vals(df)
    df <- df
  }
  
  n_foods <- length(df$shorter_desc)
  nut_quo <- quo(nut_df$nutrient)
  
  dir_mr <- rep("<", nut_df %>% filter(is_must_restrict == TRUE) %>% ungroup() %>% count() %>% as_vector())       # And less than on all the must_restricts
  dir_pos <- rep(">", nut_df %>% filter(is_must_restrict == FALSE) %>% ungroup() %>% count() %>% as_vector())     # Final menu must be greater than on all the positives
  
  dir <- c(dir_mr, dir_pos)
  rhs <- nut_df[["value"]]      # The right-hand side of the equation is all of the min or max nutrient values
  obj_fn <- df[["cost"]]             # Objective function will be to minimize total cost
  
  bounds <- list(lower = list(ind = seq(n_foods), 
                              val = rep(min_food_amount, n_foods)),
                 upper = list(ind = seq(n_foods), 
                              val = rep(max_food_amount, n_foods)))
  
  construct_matrix <- function(df, nut_df) {       # Set up matrix constraints
    mat_base <- df %>% select(!!nut_quo) %>% as_vector()    # Get a vector of all our nutrients
    mat <- matrix(mat_base, nrow = nrow(nut_df), byrow = TRUE)       # One row per constraint, one column per food (variable)
    return(mat)
  }
  
  const_mat_names <- str_c(df$shorter_desc,  # Use combo of shorter_desc and NDB_No
        df$NDB_No, sep = ", ")  # so that names are interpretable but also unique
  
  mat <- construct_matrix(df, nut_df)
  constraint_matrix <- mat %>% dplyr::as_data_frame() 
  names(constraint_matrix) <- const_mat_names
  
  constraint_matrix <- constraint_matrix %>% 
    mutate(
      dir = dir,
      rhs = rhs
    ) %>% left_join(nut_df, by = c("rhs" = "value")) %>% 
    select(nutrient, everything())
  
  if(only_full_servings == TRUE) {    # Integer values of coefficients if only full servings
    types <- rep("I", n_foods)
  } else {
    types <- rep("C", n_foods)
  }
  
  if(v_v_verbose == TRUE) {
    v_v_verbose <- TRUE
    message("Constraint matrix below:")
    print(constraint_matrix)
  } else {
    v_v_verbose <- FALSE
  }
  
  out <- Rglpk_solve_LP(obj_fn, mat, dir, rhs,                    # Do the solving; we get a list back
                        bounds = bounds, types = types, 
                        max = maximize, verbose = v_v_verbose)   
  
  out <- append(append(append(                                           # Append the dataframe of all min/max nutrient values
    out, list(necessary_nutrients = nut_df)),
    list(constraint_matrix = constraint_matrix)),                        # our constraint matrix
    list(original_menu_raw = df))                                            # and our original menu
  
  if (!is.null(df_per_100g)) {
    out <- append(out, list(original_menu_per_g = df_per_100g))
  }
  
  if (verbose == TRUE) {
    message(paste0("Cost is $", round(out$optimum, digits = 2), ".")) 
    if (out$status == 0) {
      message("Optimal solution found :)")
    } else {
      message("No optimal solution found :'(")
    }
  }
  
  return(out)
}

Let's try it out.

our_menu_solution <- our_random_menu %>% solve_it()
## Cost is $101.44.
## No optimal solution found :'(
our_menu_solution
## $optimum
## [1] 101.44
## 
## $solution
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $status
## [1] 1
## 
## $solution_dual
##  [1] 5.43 3.75 1.47 4.71 5.25 3.35 7.10 3.60 9.97 5.14 9.75 2.31 6.17 3.07
## [15] 8.25 4.42 8.42 9.28
## 
## $auxiliary
## $auxiliary$primal
##  [1]   91.337680 4269.667000  399.285000   38.956894  143.787350
##  [6] 1019.891500   21.990730  432.931500 2032.328000 3677.960000
## [11]   16.206145    2.316085    3.516593  173.897050   70.397550
## [16]    2.861833    2.313175   34.651609    5.334605    2.381441
## [21] 2681.570000
## 
## $auxiliary$dual
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 
## 
## $necessary_nutrients
## # A tibble: 21 x 3
##         nutrient value is_must_restrict
##            <chr> <dbl>            <lgl>
##  1   Lipid_Tot_g    65             TRUE
##  2     Sodium_mg  2400             TRUE
##  3  Cholestrl_mg   300             TRUE
##  4      FA_Sat_g    20             TRUE
##  5     Protein_g    56            FALSE
##  6    Calcium_mg  1000            FALSE
##  7       Iron_mg    18            FALSE
##  8  Magnesium_mg   400            FALSE
##  9 Phosphorus_mg  1000            FALSE
## 10  Potassium_mg  3500            FALSE
## # ... with 11 more rows
## 
## $constraint_matrix
## # A tibble: 45 x 22
##         nutrient `CANDIES, 19248` `SALMON, 15182` `BROADBEANS, 11088`
##            <chr>            <dbl>           <dbl>               <dbl>
##  1   Lipid_Tot_g           13.800          6.2135             0.65400
##  2     Sodium_mg           15.600         63.7500            54.50000
##  3  Cholestrl_mg            4.400         37.4000             0.00000
##  4      FA_Sat_g            8.240          1.3974             0.15042
##  5     Niacin_mg            8.240          1.3974             0.15042
##  6     Protein_g            2.232         17.3995             6.10400
##  7    Calcium_mg           43.600        203.1500            23.98000
##  8 Phosphorus_mg           43.600        203.1500            23.98000
##  9       Iron_mg            0.532          0.9010             2.07100
## 10  Magnesium_mg           15.200         24.6500            41.42000
## # ... with 35 more rows, and 18 more variables: `GELATIN DSSRT,
## #   19172` <dbl>, `CRAYFISH, 15146` <dbl>, `CRACKERS, 18965` <dbl>,
## #   `CHEESE, 01039` <dbl>, `SPICES, 02003` <dbl>, `RAVIOLI, 22899` <dbl>,
## #   `LUXURY LOAF, 07060` <dbl>, `CANDIES, 19080` <dbl>, `ONIONS,
## #   11290` <dbl>, `PIZZA HUT 14" PEPPERONI PIZZA, 21297` <dbl>, `BEANS,
## #   11046` <dbl>, `GRAPE JUC, 09135` <dbl>, `PIE, 18302` <dbl>, `PORK,
## #   10214` <dbl>, `FAST FOODS, 21114` <dbl>, dir <chr>, rhs <dbl>,
## #   is_must_restrict <lgl>
## 
## $original_menu_raw
## # A tibble: 18 x 30
##                        shorter_desc GmWt_1 serving_gmwt  cost Lipid_Tot_g
##                               <chr>  <dbl>        <dbl> <dbl>       <dbl>
##  1                          CANDIES  40.00        40.00  5.43    13.80000
##  2                           SALMON  85.00        85.00  3.75     6.21350
##  3                       BROADBEANS 109.00       109.00  1.47     0.65400
##  4                    GELATIN DSSRT  85.00        85.00  4.71     0.00000
##  5                         CRAYFISH  85.00        85.00  5.25     1.02000
##  6                         CRACKERS  30.00        30.00  3.35     3.50100
##  7                           CHEESE  28.35        28.35  7.10     8.68644
##  8                           SPICES   0.70         0.70  3.60     0.02849
##  9                          RAVIOLI 242.00       242.00  9.97     3.50900
## 10                      LUXURY LOAF  28.00        28.00  5.14     1.34400
## 11                          CANDIES  14.50        14.50  9.75     4.35000
## 12                           ONIONS 210.00       210.00  2.31     0.10500
## 13 "PIZZA HUT 14\" PEPPERONI PIZZA" 113.00       113.00  6.17    14.76910
## 14                            BEANS 104.00       104.00  3.07     0.72800
## 15                        GRAPE JUC 253.00       253.00  8.25     0.32890
## 16                              PIE  28.35        28.35  4.42     3.54375
## 17                             PORK  85.00        85.00  8.42     2.20150
## 18                       FAST FOODS 226.00       226.00  9.28    26.55500
## # ... with 25 more variables: Sodium_mg <dbl>, Cholestrl_mg <dbl>,
## #   FA_Sat_g <dbl>, Protein_g <dbl>, Calcium_mg <dbl>, Iron_mg <dbl>,
## #   Magnesium_mg <dbl>, Phosphorus_mg <dbl>, Potassium_mg <dbl>,
## #   Zinc_mg <dbl>, Copper_mg <dbl>, Manganese_mg <dbl>, Selenium_µg <dbl>,
## #   Vit_C_mg <dbl>, Thiamin_mg <dbl>, Riboflavin_mg <dbl>,
## #   Niacin_mg <dbl>, Panto_Acid_mg <dbl>, Vit_B6_mg <dbl>,
## #   Energ_Kcal <dbl>, solution_amounts <dbl>, Shrt_Desc <chr>,
## #   NDB_No <chr>, score <dbl>, scaled_score <dbl>
## 
## $original_menu_per_g
## # A tibble: 18 x 30
##                        shorter_desc solution_amounts GmWt_1 serving_gmwt
##                               <chr>            <dbl>  <dbl>        <dbl>
##  1                          CANDIES                1  40.00        40.00
##  2                           SALMON                1  85.00        85.00
##  3                       BROADBEANS                1 109.00       109.00
##  4                    GELATIN DSSRT                1  85.00        85.00
##  5                         CRAYFISH                1  85.00        85.00
##  6                         CRACKERS                1  30.00        30.00
##  7                           CHEESE                1  28.35        28.35
##  8                           SPICES                1   0.70         0.70
##  9                          RAVIOLI                1 242.00       242.00
## 10                      LUXURY LOAF                1  28.00        28.00
## 11                          CANDIES                1  14.50        14.50
## 12                           ONIONS                1 210.00       210.00
## 13 "PIZZA HUT 14\" PEPPERONI PIZZA"                1 113.00       113.00
## 14                            BEANS                1 104.00       104.00
## 15                        GRAPE JUC                1 253.00       253.00
## 16                              PIE                1  28.35        28.35
## 17                             PORK                1  85.00        85.00
## 18                       FAST FOODS                1 226.00       226.00
## # ... with 26 more variables: cost <dbl>, Lipid_Tot_g <dbl>,
## #   Sodium_mg <dbl>, Cholestrl_mg <dbl>, FA_Sat_g <dbl>, Protein_g <dbl>,
## #   Calcium_mg <dbl>, Iron_mg <dbl>, Magnesium_mg <dbl>,
## #   Phosphorus_mg <dbl>, Potassium_mg <dbl>, Zinc_mg <dbl>,
## #   Copper_mg <dbl>, Manganese_mg <dbl>, Selenium_µg <dbl>,
## #   Vit_C_mg <dbl>, Thiamin_mg <dbl>, Riboflavin_mg <dbl>,
## #   Niacin_mg <dbl>, Panto_Acid_mg <dbl>, Vit_B6_mg <dbl>,
## #   Energ_Kcal <dbl>, Shrt_Desc <chr>, NDB_No <chr>, score <dbl>,
## #   scaled_score <dbl>

How long did that take?

system.time(our_random_menu %>% solve_it(verbose = FALSE))
##    user  system elapsed 
##   0.023   0.000   0.023

Not long. Thanks for being written in C, GLPK!

Solve menu

Okay so our output of solve_it() is an informative but long list. It has all the building blocks we need to create a solved menu; now we just need to extract those parts and glue them together in the right ways. Here's where solve_menu() comes in.

solve_menu() takes one main argument: the result of a call to solve_it(). Since we've written the return value of solve_it() to contain the original menu and a vector of solution amounts -- that is, the amount we're multiplying each portion size by in order to arrive at our solution -- we can combine these to get our solved menu.

We also return a message, if verbose is TRUE, telling us which food we've got the most servings of, as this might be something we'd want to decrease. (Now that I'm thinking about it, maybe a more helpful message would take a threshold portion size and only alert us if we've exceeded that threshold.)

solve_menu <- function(sol, verbose = TRUE) {
  
  solved_col <-  tibble(solution_amounts = sol$solution)    # Grab the vector of solution amounts
  
  if (! "solution_amounts" %in% names(sol$original_menu_per_g)) {   # If we don't yet have a solution amounts column add it
    df_solved <- sol$original_menu_per_g %>% 
      bind_cols(solved_col)            # cbind that to the original menu
  } else {
    df_solved <- sol$original_menu_per_g %>% 
      mutate(
        solution_amounts = solved_col %>% as_vector()    # If we've already got a solution amounts column, replace the old one with the new
      ) 
  }
  
  df_solved <- df_solved %>% 
    mutate(
      GmWt_1 = GmWt_1 * solution_amounts,
      cost = cost * solution_amounts
    ) %>% 
    select(shorter_desc, solution_amounts, GmWt_1, serving_gmwt, everything()) 
  
  max_food <- df_solved %>%                                   # Find what the most of any one food we've got is
    filter(solution_amounts == max(df_solved$solution_amounts)) %>% 
    slice(1:1)                                           # If we've got multiple maxes, take only the first
  
  if (verbose == TRUE) {
    message(paste0("We've got a lot of ", max_food$shorter_desc %>% as_vector()), ". ",
            max_food$solution_amounts %>% round(digits = 2), " servings of ",
            max_food$shorter_desc %>% as_vector() %>% is_plural(return_bool = FALSE), ".")  
  }
  
  return(df_solved)
}

Let's see what our tidied output looks like.

our_solved_menu <- our_menu_solution %>% solve_menu()
## We've got a lot of CANDIES. 1 servings of them.
our_solved_menu %>% kable(format = "html")
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
CANDIES 1 40.00 40.00 2.24 34.50 39 11 20.600 5.58 109 1.33 38 129 306 0.79 0.180 0.000 0.0 1.0 0.040 0.160 0.330 0.250 0.060 563 CANDIES,HERSHEY'S,ALMOND JOY BITES 19248 -4490.352 0.0062599
SALMON 1 85.00 85.00 9.02 7.31 75 44 1.644 20.47 239 1.06 29 326 377 1.02 0.084 0.030 35.4 0.0 0.016 0.193 5.480 0.550 0.300 153 SALMON,SOCKEYE,CND,WO/SALT,DRND SOL W/BONE 15182 -3913.498 1.1090489
BROADBEANS 1 109.00 109.00 3.93 0.60 50 0 0.138 5.60 22 1.90 38 95 250 0.58 0.074 0.320 1.2 33.0 0.170 0.110 1.500 0.086 0.038 72 BROADBEANS,IMMAT SEEDS,RAW 11088 -4250.264 0.4652428
GELATIN DSSRT 1 85.00 85.00 1.96 0.00 466 0 0.000 7.80 3 0.13 2 141 7 0.01 0.118 0.011 6.7 0.0 0.003 0.041 0.009 0.014 0.001 381 GELATIN DSSRT,DRY MIX 19172 -4938.439 -0.8503611
CRAYFISH 1 85.00 85.00 6.64 1.20 94 133 0.181 16.77 60 0.83 33 270 296 1.76 0.685 0.522 36.7 0.9 0.050 0.085 2.280 0.580 0.076 82 CRAYFISH,MXD SP,WILD,CKD,MOIST HEAT 15146 -4266.922 0.4333988
CRACKERS 1 30.00 30.00 3.41 11.67 1167 0 3.333 10.00 67 4.80 22 162 141 0.90 0.118 0.517 26.2 0.0 1.087 0.750 7.170 0.528 0.044 418 CRACKERS,CHS,RED FAT 18965 -4906.367 -0.7890484
CHEESE 1 28.35 28.35 7.25 30.64 1809 90 19.263 21.54 662 0.56 30 392 91 2.08 0.034 0.030 14.5 0.0 0.040 0.586 0.734 1.731 0.124 369 CHEESE,ROQUEFORT 01039 -4892.506 -0.7625507
SPICES 1 0.70 0.70 9.03 4.07 76 0 2.157 22.98 2240 89.80 711 274 2630 7.10 2.100 9.800 3.0 0.8 0.080 1.200 4.900 0.838 1.340 233 SPICES,BASIL,DRIED 02003 -4643.583 -0.2866766
RAVIOLI 1 242.00 242.00 7.63 1.45 306 3 0.723 2.48 33 0.74 15 50 232 0.36 0.142 0.176 3.5 0.0 0.074 0.080 1.060 0.272 0.102 77 RAVIOLI,CHEESE-FILLED,CND 22899 -4617.693 -0.2371810
LUXURY LOAF 1 28.00 28.00 4.64 4.80 1225 36 1.580 18.40 36 1.05 20 185 377 3.05 0.100 0.041 21.5 0.0 0.707 0.297 3.482 0.515 0.310 141 LUXURY LOAF,PORK 07060 -4852.980 -0.6869870
CANDIES 1 14.50 14.50 9.99 30.00 11 0 17.750 4.20 32 3.13 115 132 365 1.62 0.700 0.800 4.2 0.0 0.055 0.090 0.427 0.105 0.035 480 CANDIES,SEMISWEET CHOC 19080 -4597.911 -0.1993645
ONIONS 1 210.00 210.00 4.83 0.05 8 0 0.009 0.71 27 0.34 8 2 101 0.09 0.024 0.040 0.4 5.1 0.016 0.018 0.132 0.078 0.070 28 ONIONS,FRZ,WHL,CKD,BLD,DRND,WO/SALT 11290 -4397.386 0.1839856
PIZZA HUT 14" PEPPERONI PIZZA 1 113.00 113.00 5.69 13.07 676 23 4.823 11.47 147 2.57 22 193 187 1.36 0.104 0.425 15.5 1.0 0.420 0.210 3.750 0.323 0.090 291 PIZZA HUT 14" PEPPERONI PIZZA,PAN CRUST 21297 -4832.658 -0.6481376
BEANS 1 104.00 104.00 8.77 0.70 13 0 0.085 6.15 15 1.93 101 100 307 0.89 0.356 0.408 0.6 18.8 0.390 0.215 1.220 0.825 0.191 67 BEANS,NAVY,MATURE SEEDS,SPROUTED,RAW 11046 -4122.162 0.7101393
GRAPE JUC 1 253.00 253.00 4.53 0.13 5 0 0.025 0.37 11 0.25 10 14 104 0.07 0.018 0.239 0.0 0.1 0.017 0.015 0.133 0.048 0.032 60 GRAPE JUC,CND OR BTLD,UNSWTND,WO/ ADDED VIT C 09135 -4343.103 0.2877596
PIE 1 28.35 28.35 1.75 12.50 211 0 3.050 2.40 7 1.12 7 28 79 0.19 0.053 0.185 7.8 1.7 0.148 0.107 1.230 0.093 0.032 265 PIE,APPL,PREP FROM RECIPE 18302 -4710.654 -0.4148992
PORK 1 85.00 85.00 2.46 2.59 63 63 0.906 22.81 9 0.56 23 251 354 1.72 0.069 0.011 37.4 0.0 0.610 0.256 7.348 0.728 0.611 121 PORK,FRSH,LOIN,SIRLOIN (CHOPS OR ROASTS),BNLESS,LN,RAW 10214 -4192.317 0.5760225
FAST FOODS 1 226.00 226.00 5.02 11.75 350 54 4.654 15.17 45 2.59 22 139 252 2.51 0.097 0.110 11.3 0.5 0.160 0.170 3.350 0.240 0.240 239 FAST FOODS,HAMBURGER; DOUBLE,LRG PATTY; W/ CONDMNT & VEG 21114 -4517.685 -0.0459943

Solve nutrients

We'll want to do something with nutrients that's analogous to what we're doing in solve_menu(). This function will let us find what the raw nutrient amounts in our solved menu are, and let us know which nutrient we've overshot the lower bound on the most. Like solve_menu(), a result from solve_it() can be piped nicely in here.

One part of the solution returned by the solver is a vector of the values of the constraints -- that is, our nutrients -- at solution. That lives in $auxiliary$primal and becomes our solved_nutrient_value in the function below.

Recall also that we took nut_df, the dataframe of nutritional requirements handed to us by the user, and appended it to the solution so that it's also returned as a result of our call to solve_it(). This means the outcome of solve_it() will let us compare the required_value for each nutrient to its solved_nutrient_value. We calculate the ratio of these two for every nutrient, and if verbose is TRUE, let the user know which nutrient they've overshot the daily minimum on the most.

solve_nutrients <- function(sol, verbose = TRUE) {
  
  solved_nutrient_value <- list(solution_nutrient_value =       # Grab the vector of nutrient values in the solution
                              sol$auxiliary$primal) %>% as_tibble()
  
  nut_df_small_solved <- sol$necessary_nutrients %>%       # cbind it to the nutrient requirements
    bind_cols(solved_nutrient_value)  %>% 
    rename(
      required_value = value
    ) %>% 
    select(nutrient, is_must_restrict, required_value, solution_nutrient_value)
  
  ratios <- nut_df_small_solved %>%                # Find the solution:required ratios for each nutrient
    mutate(
      ratio = solution_nutrient_value/required_value
    )
  
  max_pos_overshot <- ratios %>%             # Find where we've overshot our positives the most
    filter(is_must_restrict == FALSE) %>% 
    filter(ratio == max(.$ratio))
  
  if (verbose == TRUE) {
    message(paste0("We've overshot the most on ", max_pos_overshot$nutrient %>% as_vector()), 
            ". It's ", 
        max_pos_overshot$ratio %>% round(digits = 2), " times what is needed.")
  }
  
  return(nut_df_small_solved)
}

Remember that we saved the result of solve_it() in our_menu_solution. Let's see what those ratios look like in our solution.

our_menu_solution %>% solve_nutrients() %>% kable(format = "html")
## We've overshot the most on Protein_g. It's 2.57 times what is needed.
nutrient is_must_restrict required_value solution_nutrient_value
Lipid_Tot_g TRUE 65 91.337680
Sodium_mg TRUE 2400 4269.667000
Cholestrl_mg TRUE 300 399.285000
FA_Sat_g TRUE 20 38.956894
Protein_g FALSE 56 143.787350
Calcium_mg FALSE 1000 1019.891500
Iron_mg FALSE 18 21.990730
Magnesium_mg FALSE 400 432.931500
Phosphorus_mg FALSE 1000 2032.328000
Potassium_mg FALSE 3500 3677.960000
Zinc_mg FALSE 15 16.206145
Copper_mg FALSE 2 2.316085
Manganese_mg FALSE 2 3.516592
Selenium_µg FALSE 70 173.897050
Vit_C_mg FALSE 60 70.397550
Thiamin_mg FALSE 2 2.861833
Riboflavin_mg FALSE 2 2.313175
Niacin_mg FALSE 20 34.651609
Panto_Acid_mg FALSE 10 5.334605
Vit_B6_mg FALSE 2 2.381441
Energ_Kcal FALSE 2300 2681.570000

Swapping

Our menus often aren't solvable. That is, at the minimum portion size we set, there's no way we can change portion sizes in such a way that we stay under all the maximum values for each must restrict and under the minimums for all positive nutrients as well as have enough calories.

In these cases, we'll need to change up our lineup.

Single Swap

I only use single swapping for the cases where we're above the max threshold on must restricts, but you could imagine implementing the same functions to deal with positives.

The idea with a single swap is to see which must restricts are not satisfied, find the food that is the max_offender on that must restrict (i.e., contributes the most in absolute terms to the value of the must restrict) and then swap it out. We try to replace_food_w_better(), that is, swap it out for a food from a pool of better foods on that dimension. We define better as foods that score above a user-specified z-score cutoff on that must_restrict. If there are no foods that satisfy that cutoff, we choose a food at random from the pool of all possible foods.

smart_swap_single <- function(menu, max_offender, cutoff = 0.5, df = abbrev, verbose = FALSE) {
  
  swap_count <- 0

    for (m in seq_along(mr_df$must_restrict)) {   
      nut_to_restrict <- mr_df$must_restrict[m]    # grab the name of the nutrient we're restricting
      message(paste0("------- The nutrient we're restricting is ", nut_to_restrict, ". It has to be below ", mr_df$value[m])) 
      to_restrict <- (sum(menu[[nut_to_restrict]] * menu$GmWt_1, na.rm = TRUE))/100   # get the amount of that must restrict nutrient in our original menu
      message(paste0("The original total value of that nutrient in our menu is ", to_restrict)) 
      
      if (to_restrict > mr_df$value[m]) {     # if the amount of the must restrict in our current menu is above the max value it should be according to mr_df
        swap_count <- swap_count + 1
        
        max_offender <- which(menu[[nut_to_restrict]] == max(menu[[nut_to_restrict]]))   # Find the food that's the worst offender in this respect
        
        message(paste0("The worst offender in this respect is ", menu[max_offender, ]$Shrt_Desc))
        
        menu[max_offender, ] <- replace_food_w_better(menu, max_offender, 
                                                           nutrient_to_restrict = nut_to_restrict, cutoff = cutoff)
        
        to_restrict <- (sum(menu[[nut_to_restrict]] * menu$GmWt_1, na.rm = TRUE))/100   # recalculate the must restrict nutrient content
        message(paste0("Our new value of this must restrict is ", to_restrict)) 
      } else {
        message("We're all good on this nutrient.") 
      }
    }
  
  if (verbose == TRUE) {
    print(paste0(swap_count, " swaps were completed."))
  }
  
  return(menu)
}

do_single_swap <- function(menu, solve_if_unsolved = TRUE, verbose = FALSE,
                        new_solution_amount = 1){  # What should the solution amount of the newly swapped in foods be?
  
  if (verbose == FALSE) {
    out <- suppressWarnings(suppressMessages(menu %>% 
      smart_swap_single(menu))) 
  } else {
    out <- menu
      smart_swap_single(menu) 
  }

  return(out)
}

In practice, that looks like:

our_random_menu %>% do_single_swap() %>% kable(format = "html")
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
RUTABAGAS 1 170.00 170.00 7.72 0.18 5 0 0.029 0.93 18 0.18 10 41 216 0.12 0.029 0.097 0.7 18.8 0.082 0.041 0.715 0.155 0.102 30 RUTABAGAS,CKD,BLD,DRND,WO/SALT 11436 -2861.039 0.6147894
SALMON 1 85.00 85.00 3.75 7.31 75 44 1.644 20.47 239 1.06 29 326 377 1.02 0.084 0.030 35.4 0.0 0.016 0.193 5.480 0.550 0.300 153 SALMON,SOCKEYE,CND,WO/SALT,DRND SOL W/BONE 15182 -2602.498 1.1090489
BROADBEANS 1 109.00 109.00 1.47 0.60 50 0 0.138 5.60 22 1.90 38 95 250 0.58 0.074 0.320 1.2 33.0 0.170 0.110 1.500 0.086 0.038 72 BROADBEANS,IMMAT SEEDS,RAW 11088 -2939.264 0.4652428
GELATIN DSSRT 1 85.00 85.00 4.71 0.00 466 0 0.000 7.80 3 0.13 2 141 7 0.01 0.118 0.011 6.7 0.0 0.003 0.041 0.009 0.014 0.001 381 GELATIN DSSRT,DRY MIX 19172 -3627.439 -0.8503611
COWPEAS 1 171.00 171.00 6.26 0.71 255 0 0.185 8.13 26 3.05 96 142 375 1.87 0.271 0.473 2.5 0.4 0.162 0.046 0.714 0.386 0.092 117 COWPEAS,CATJANG,MATURE SEEDS,CKD,BLD,W/SALT 16361 -2687.950 0.9456888
CRACKERS 1 30.00 30.00 3.35 11.67 1167 0 3.333 10.00 67 4.80 22 162 141 0.90 0.118 0.517 26.2 0.0 1.087 0.750 7.170 0.528 0.044 418 CRACKERS,CHS,RED FAT 18965 -3595.367 -0.7890484
PORK 1 113.00 113.00 8.59 32.93 94 100 11.311 22.83 20 1.15 19 192 280 2.58 0.038 0.013 38.0 0.0 0.341 0.488 7.522 0.984 0.508 393 PORK,GROUND,72% LN / 28% FAT,CKD,CRUMBLES 10974 -2981.649 0.3842142
SPICES 1 0.70 0.70 3.60 4.07 76 0 2.157 22.98 2240 89.80 711 274 2630 7.10 2.100 9.800 3.0 0.8 0.080 1.200 4.900 0.838 1.340 233 SPICES,BASIL,DRIED 02003 -3332.583 -0.2866766
RAVIOLI 1 242.00 242.00 9.97 1.45 306 3 0.723 2.48 33 0.74 15 50 232 0.36 0.142 0.176 3.5 0.0 0.074 0.080 1.060 0.272 0.102 77 RAVIOLI,CHEESE-FILLED,CND 22899 -3306.693 -0.2371810
LUXURY LOAF 1 28.00 28.00 5.14 4.80 1225 36 1.580 18.40 36 1.05 20 185 377 3.05 0.100 0.041 21.5 0.0 0.707 0.297 3.482 0.515 0.310 141 LUXURY LOAF,PORK 07060 -3541.980 -0.6869870
BEVERAGES 1 240.00 240.00 3.77 1.04 63 0 0.000 0.42 188 0.30 7 8 50 0.63 0.017 0.033 0.1 0.0 0.015 0.177 0.075 0.009 0.003 38 BEVERAGES,ALMOND MILK,SWTND,VANILLA FLAVOR,RTD 14016 -2916.226 0.5092852
ONIONS 1 210.00 210.00 2.31 0.05 8 0 0.009 0.71 27 0.34 8 2 101 0.09 0.024 0.040 0.4 5.1 0.016 0.018 0.132 0.078 0.070 28 ONIONS,FRZ,WHL,CKD,BLD,DRND,WO/SALT 11290 -3086.386 0.1839856
PIZZA HUT 14" PEPPERONI PIZZA 1 113.00 113.00 6.17 13.07 676 23 4.823 11.47 147 2.57 22 193 187 1.36 0.104 0.425 15.5 1.0 0.420 0.210 3.750 0.323 0.090 291 PIZZA HUT 14" PEPPERONI PIZZA,PAN CRUST 21297 -3521.658 -0.6481376
BEANS 1 104.00 104.00 3.07 0.70 13 0 0.085 6.15 15 1.93 101 100 307 0.89 0.356 0.408 0.6 18.8 0.390 0.215 1.220 0.825 0.191 67 BEANS,NAVY,MATURE SEEDS,SPROUTED,RAW 11046 -2811.162 0.7101393
GRAPE JUC 1 253.00 253.00 8.25 0.13 5 0 0.025 0.37 11 0.25 10 14 104 0.07 0.018 0.239 0.0 0.1 0.017 0.015 0.133 0.048 0.032 60 GRAPE JUC,CND OR BTLD,UNSWTND,WO/ ADDED VIT C 09135 -3032.103 0.2877596
PIE 1 28.35 28.35 4.42 12.50 211 0 3.050 2.40 7 1.12 7 28 79 0.19 0.053 0.185 7.8 1.7 0.148 0.107 1.230 0.093 0.032 265 PIE,APPL,PREP FROM RECIPE 18302 -3399.654 -0.4148992
PORK 1 85.00 85.00 8.42 2.59 63 63 0.906 22.81 9 0.56 23 251 354 1.72 0.069 0.011 37.4 0.0 0.610 0.256 7.348 0.728 0.611 121 PORK,FRSH,LOIN,SIRLOIN (CHOPS OR ROASTS),BNLESS,LN,RAW 10214 -2881.317 0.5760225
FAST FOODS 1 226.00 226.00 9.28 11.75 350 54 4.654 15.17 45 2.59 22 139 252 2.51 0.097 0.110 11.3 0.5 0.160 0.170 3.350 0.240 0.240 239 FAST FOODS,HAMBURGER; DOUBLE,LRG PATTY; W/ CONDMNT & VEG 21114 -3206.685 -0.0459943

Wholesale Swap

The wholesale swap takes a different approach. It uses knowledge we've gained from solving to keep the foods that the solver wanted more of and offer the rest up for swapping. The intuition here is that foods that the solver increased the portion sizes of are more valuable to the menu as a whole. We sample a percent_to_swap of the foods that the solver assigned the lowest portion size to, shamefully dubbed the worst_foods.

wholesale_swap <- function(menu, df = abbrev, percent_to_swap = 0.5) {
  
  # Get foods with the lowest solution amounts
  min_solution_amount <- min(menu$solution_amounts)
  worst_foods <- menu %>% 
    filter(solution_amounts == min_solution_amount)
  
  if (nrow(worst_foods) >= 2) {
    to_swap_out <- worst_foods %>% sample_frac(percent_to_swap)
    message(paste0("Swapping out a random ", percent_to_swap*100, "% of foods: ", 
                   str_c(to_swap_out$shorter_desc, collapse = ", ")))
    
  } else if (nrow(worst_foods) == 1)  {
    message("Only one worst food. Swapping this guy out.")
    to_swap_out <- worst_foods
    
  } else {
    message("No worst foods")
  }
  
  get_swap_candidates <- function(df, to_swap_out) {
    candidate <- df %>% 
      filter(! (NDB_No %in% menu)) %>%    # We can't swap in a food that already exists in our menu
      sample_n(., size = nrow(to_swap_out)) %>% 
      mutate(solution_amounts = 1)    # Give us one serving of each of these new foods
    return(candidate)
  }
  swap_candidate <- get_swap_candidates(df = df, to_swap_out = to_swap_out)
  
  if (score_menu(swap_candidate) < score_menu(to_swap_out)) {
    message("Swap candidate not good enough; reswapping.")
    swap_candidate <- get_swap_candidates(df = df, to_swap_out = to_swap_out)
    
  } else {
      message("Swap candidate is good enough. Doing the wholesale swap.")
      return(swap_candidate)
  }
  
  newly_swapped_in <- get_swap_candidates(df, to_swap_out)
  
  message(paste0("Replacing with: ", 
                 str_c(newly_swapped_in$shorter_desc, collapse = ", ")))
  
  out <- menu %>% 
    filter(!NDB_No %in% worst_foods) %>% 
    bind_rows(newly_swapped_in)
  
  return(out)
}

Let's do a wholesale swap.

our_random_menu %>% wholesale_swap() %>% kable(format = "html")
## Swapping out a random 50% of foods: LUXURY LOAF, CRACKERS, CANDIES, SPICES, PORK, GRAPE JUC, PIE, PIZZA HUT 14" PEPPERONI PIZZA, FAST FOODS
## Swap candidate is good enough. Doing the wholesale swap.
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
POTATOES 1 245.00 245.00 4.01 3.68 335 12 2.255 2.87 57 0.57 19 63 378 0.40 0.163 0.166 1.6 10.6 0.069 0.092 1.053 0.514 0.178 88 POTATOES,SCALLPD,HOME-PREPARED W/BUTTER 11372 -2927.267 0.4881786
TURKEY BREAST 1 85.00 85.00 7.74 3.46 397 42 0.980 22.16 9 0.66 21 214 248 1.53 0.041 0.015 25.7 0.0 0.053 0.133 9.067 0.489 0.320 126 TURKEY BREAST,PRE-BASTED,MEAT&SKN,CKD,RSTD 05293 -3281.581 -0.1891749
SNACKS 1 28.35 28.35 6.90 49.60 1531 133 20.800 21.50 68 3.40 21 180 257 2.42 0.130 0.086 NA 6.8 0.141 0.436 4.540 0.328 0.205 550 SNACKS,BF STKS,SMOKED 19407 -3705.245 -0.9991068
BEEF 1 85.00 85.00 8.03 7.67 81 70 3.419 20.87 15 2.31 18 195 330 7.86 0.092 0.014 21.9 0.0 0.070 0.190 4.097 0.680 0.355 152 BEEF,CHUCK EYE COUNTRY-STYLE RIBS,BNLESS,LN,0" FAT,CHOIC,RAW 23072 -2987.803 0.3724493
WHEAT FLR 1 NA NA 9.73 1.45 2 NA 0.268 11.50 20 5.06 30 112 138 0.84 0.161 0.679 27.5 0.0 0.736 0.445 5.953 0.405 0.032 363 WHEAT FLR,WHITE (INDUSTRIAL),11.5% PROT,UNBLEACHED,ENR 20636 -3374.000 -0.3658548
OKARA 1 122.00 122.00 9.38 1.73 9 0 0.193 3.52 80 1.30 26 60 213 0.56 0.200 0.404 10.6 0.0 0.020 0.020 0.100 0.088 0.115 76 OKARA 16130 -2904.295 0.5320946
CHICKEN 1 140.00 140.00 1.78 13.60 82 88 3.790 27.30 15 1.26 23 182 223 1.94 0.066 0.020 23.9 0.0 0.063 0.168 8.487 1.030 0.400 239 CHICKEN,BROILERS OR FRYERS,MEAT&SKN,CKD,RSTD 05009 -2925.658 0.4912538
FISH 1 NA NA 6.40 12.95 870 67 2.440 23.19 55 0.55 29 270 390 0.77 0.148 0.037 30.5 0.0 0.043 0.201 8.610 0.822 0.378 209 FISH,SALMON,KING,W/ SKN,KIPPERED,(ALASKA NATIVE) 35168 -3374.000 -0.3658548
KRAFT 1 28.00 28.00 8.07 4.10 1532 4 0.800 12.60 63 4.32 NA 129 267 NA NA NA NA 3.3 0.390 0.290 3.840 NA NA 381 KRAFT,STOVE TOP STUFFING MIX CHICKEN FLAVOR 18567 -3670.005 -0.9317363

Full Solving

fully_solved <- build_menu() %>% solve_full(verbose = FALSE)
fully_solved %>% kable(format = "html")
shorter_desc solution_amounts GmWt_1 serving_gmwt cost Lipid_Tot_g Sodium_mg Cholestrl_mg FA_Sat_g Protein_g Calcium_mg Iron_mg Magnesium_mg Phosphorus_mg Potassium_mg Zinc_mg Copper_mg Manganese_mg Selenium_µg Vit_C_mg Thiamin_mg Riboflavin_mg Niacin_mg Panto_Acid_mg Vit_B6_mg Energ_Kcal Shrt_Desc NDB_No score scaled_score
LAMB 1 114.00000 114.0 1.540000 4.73 73 76 1.857 23.42 2 2.27 26 210 311 2.32 0.128 0.006 8.1 0.0 0.147 0.380 8.490 0.880 0.556 136 LAMB,AUSTRALIAN,IMP,FRSH,TENDERLOIN,BNLESS,LN,1/8" FAT,RAW 17443 -2872.275 0.5933092
CHI FORMU 1 31.00000 31.0 8.040000 4.70 36 2 1.256 2.80 92 1.32 19 80 124 0.56 0.106 0.144 3.0 9.6 0.256 0.200 0.960 0.960 0.248 99 CHI FORMU,ABBT NUTR,PEDIASU,RTF,W/ IRON & FIB (FORMER ROSS) 03870 -3283.729 -0.1932802
PUMPKIN LEAVES 1 360.87918 39.0 76.987558 0.40 11 0 0.207 3.15 39 2.22 38 104 436 0.20 0.133 0.355 0.9 11.0 0.094 0.128 0.920 0.042 0.207 19 PUMPKIN LEAVES,RAW 11418 -3130.351 0.0999373
TANGERINES 1 2958.09865 252.0 29.111447 0.10 6 0 0.012 0.45 7 0.37 8 10 78 0.24 0.044 0.032 0.4 19.8 0.053 0.044 0.445 0.125 0.042 61 TANGERINES,(MANDARIN ORANGES),CND,LT SYRUP PK 09220 -3074.289 0.2071124
CEREALS RTE 1 39.30006 27.0 7.874567 4.50 639 0 1.000 8.80 370 16.70 89 222 230 13.90 0.253 2.320 18.5 22.2 1.400 1.600 18.500 0.749 1.852 378 CEREALS RTE,GENERAL MILLS,BERRY BURST CHEERIOS,TRIPLE BERRY 08239 -3273.216 -0.1731829
EGG CUSTARDS 1 141.00000 141.0 6.590000 2.83 87 49 1.475 4.13 146 0.33 17 137 214 0.61 0.012 0.016 4.9 0.2 0.056 0.235 0.135 0.683 0.066 112 EGG CUSTARDS,DRY MIX,PREP W/ 2% MILK 19205 -2831.054 0.6721117
PUDDINGS 1 140.00000 140.0 4.490000 2.90 156 9 1.643 2.80 99 0.04 9 74 119 0.33 0.025 0.005 3.3 0.0 0.036 0.149 0.078 0.326 0.028 113 PUDDINGS,VANILLA,DRY MIX,REG,PREP W/ WHL MILK 19207 -3179.996 0.0050279
INF FORMULA 1 9.60000 9.6 2.430000 27.65 154 0 11.430 15.36 998 10.20 46 666 768 3.84 0.461 0.026 9.2 61.0 0.512 0.768 5.376 2.304 0.307 512 INF FORMULA, ABB NUTR, SIMIL, GO & GR, PDR, W/ ARA & DHA 33871 -3144.150 0.0735572
SEA BASS 1 129.00000 129.0 3.960000 2.00 68 41 0.511 18.43 10 0.29 41 194 256 0.40 0.019 0.015 36.5 0.0 0.110 0.120 1.600 0.750 0.400 97 SEA BASS,MXD SP,RAW 15091 -2795.921 0.7392762
WHEAT FLR 1 125.00000 125.0 2.420000 0.98 2 0 0.155 10.33 15 1.17 22 108 107 0.70 0.144 0.682 33.9 0.0 0.120 0.040 1.250 0.438 0.044 364 WHEAT FLR,WHITE,ALL-PURPOSE,UNENR 20481 -3001.896 0.3455075
BEANS 1 179.00000 179.0 4.840000 0.64 238 0 0.166 8.97 73 2.84 68 169 463 1.09 0.149 0.510 1.3 0.0 0.236 0.059 0.272 0.251 0.127 142 BEANS,SML WHITE,MATURE SEEDS,CKD,BLD,W/SALT 16346 -2389.504 1.5162376
MILK 1 245.00000 245.0 7.340000 1.98 59 8 1.232 3.95 143 0.06 15 112 182 0.41 0.011 0.002 2.6 1.1 0.045 0.194 0.101 0.339 0.046 56 MILK,RED FAT,FLUID,2% MILKFAT,W/ NONFAT MILK SOL,WO/ VIT A 01152 -2416.917 1.4638299
BEANS 1 169.00000 169.0 5.520000 0.49 2 0 0.126 9.06 52 2.30 65 165 508 0.96 0.271 0.548 1.4 0.0 0.257 0.063 0.570 0.299 0.175 149 BEANS,PINK,MATURE SEEDS,CKD,BLD,WO/SALT 16041 -2016.445 2.2294253
BEEF 1 85.00000 85.0 4.560000 9.04 64 64 3.544 21.22 9 1.96 21 152 275 4.89 0.073 0.073 21.1 0.0 0.100 0.269 5.080 0.540 0.488 166 BEEF,RIB EYE STK/RST,BONE-IN,LIP-ON,LN,1/8" FAT,ALL GRDS,RAW 23150 -3057.622 0.2389742
DESSERTS 1 141.00000 141.0 8.460000 3.43 351 0 0.685 1.75 35 0.82 8 28 78 0.18 0.071 0.130 3.5 2.2 0.083 0.081 0.846 0.092 0.040 161 DESSERTS,APPL CRISP,PREPARED-FROM-RECIPE 19186 -3650.814 -0.8950487


Simulating Solving

Cool, so we've got a mechanism for creating and solving menus. But what portion of our menus are even solvable from the get-go? We'll stipulate that solvable means solvable at a minimum portion size of 1 without doing any swapping. To answer that, I set about making a way to run a some simulations.

First, a helper function for just plucking the status portion of our solve_it() response telling us whether we solved the menu or not. The result of get_status() should always be either a 1 for unsolvable or 0 for solved.

get_status <- function(seed = NULL, min_food_amount = 0.5, verbose = TRUE) {  
  this_menu <- build_menu(seed = seed) %>% 
    solve_it(min_food_amount = min_food_amount, verbose = verbose, only_full_servings = FALSE) %>% 
    purrr::pluck("status")
}

Now to the simulations: or a given minimum portion size, what proportion of a random number of simulated menus can we solve?

We'll use map_dbl to get the status of each solution in our simulation. Then all we need to do is specify a minimum portion size for all menus to have and the number of simulations to run. We'll shuffle the seed at which random menus are built for each simulation and then return a vector of whether each was solved or not.

simulate_menus <- function(n_sims = 10, min_food_amount = 0.5, verbose = FALSE) {
  
  # Choose as many random seeds as we have simulations
  seeds <- sample(1:n_sims, size = n_sims, replace = FALSE)
  
  out <- seeds %>% map2_dbl(.y = min_food_amount, .f = get_status)
  return(out)
}
simulate_menus(verbose = FALSE)
## Cost is $454.63.
## No optimal solution found :'(
## Cost is $56.04.
## No optimal solution found :'(
## Cost is $56.06.
## No optimal solution found :'(
## Cost is $170.53.
## No optimal solution found :'(
## Cost is $244.35.
## No optimal solution found :'(
## Cost is $98.86.
## Optimal solution found :)
## Cost is $133.71.
## No optimal solution found :'(
## Cost is $99.43.
## Optimal solution found :)
## Cost is $53.44.
## No optimal solution found :'(
## Cost is $149.91.
## Optimal solution found :)
##  [1] 1 1 1 1 1 0 1 0 1 0

Alright so that's kinda useful for a single portion size, but what if we wanted to see how solvability varies as we change up portion size? Presumably as we decrease the lower bound for each food's portion size we'll give ourselves more flexibility and be able to solve a higher proportion of menus. But will the percent of menus that are solvable increase linearly as we decrease portion size?

Simulate Spectrum

I named this next function simulate_spectrum() because it allows us to take a lower and an upper bound of minimum portion sizes and see what happens at each point between those two intervals.

We specify the lower bound for the min portion size spectrum with from and the upper bound with to. How spaced out those points are and how many of them there are are set with n_intervals and n_sims; in other words, n_intervals is the number of chunks we want to split the spectrum of from to to into and n_sims is the number of times we want to repeat the simulation at each point.

Instead of a vector, this time we'll return a dataframe in order to be able to match up the minimum portion size (min_amount, which we're varying) with whether or not the menu was solvable.

simulate_spectrum <- function(n_intervals = 10, n_sims = 2, from = -1, to = 1,
                              min_food_amount = NULL, verbose = FALSE, ...) {

  interval <- (to - from) / n_intervals
  spectrum <- seq(from = from, to = to, by = interval) %>% rep(n_sims) %>% sort()
  
  seeds <- sample(1:length(spectrum), size = length(spectrum), replace = FALSE)
  
  out_status <- vector(length = length(spectrum))
  
  for (i in seq_along(spectrum)) {
    this_status <- get_status(seed = seeds[i], min_food_amount = spectrum[i], verbose = verbose)
    if (!is.integer(this_status)) {
      this_status <- integer(0)     # If we don't get an integer value back, make it NA
    }
    out_status[i] <- this_status
  }
  
  out <- tibble(min_amount = spectrum, status = out_status)
  
  return(out)
}
status_spectrum <- simulate_spectrum()
status_spectrum %>% kable(format = "html")
min_amount status
-1.0 0
-1.0 0
-0.8 0
-0.8 0
-0.6 0
-0.6 0
-0.4 1
-0.4 0
-0.2 1
-0.2 0
0.0 0
0.0 0
0.2 0
0.2 0
0.4 1
0.4 0
0.6 1
0.6 1
0.8 1
0.8 1
1.0 1
1.0 1

Simulate Spectrum with Swapping

The next obvious question is, how many menus are solvable within a certain number of swaps?

get_swap_status() is analogous to get_status() from our vanilla simulator above. We specify a maximum number of allowed swaps. We build a random menu and count how many swaps it takes to solve it. If we can't solve it within max_n_swaps swaps, we'll give up. At the end, we return a tibble of the status and the number of swaps done for each food.

get_swap_status <- function(seed = NULL, min_food_amount = 0.5, max_n_swaps = 3, return_status = TRUE,
                           verbose = TRUE, ...) {  
  counter <- 0
  this_solution <- build_menu(seed = seed) %>% 
    solve_it(min_food_amount = min_food_amount, verbose = verbose, only_full_servings = FALSE) 
  
  this_status <- this_solution %>% purrr::pluck("status")
  
  this_menu <- this_solution %>% solve_menu()
  
  while (counter < max_n_swaps & this_status == 1) {
    this_solution <- this_menu %>% do_single_swap() %>% 
      solve_it(min_food_amount = min_food_amount, verbose = verbose, only_full_servings = FALSE)
    this_status <- this_solution %>% purrr::pluck("status")
    
    if (this_status == 0) {
      message(paste0("Solution found in ", counter, " steps"))
      if (return_status == TRUE) {
        out <- list(status = this_status, n_swaps_done = counter) %>% as_tibble()
        return(out)
      } else {
        this_menu <- this_solution %>% solve_menu()
        return(this_menu)
      }
    }
    counter <- counter + 1
  }
  
  message(paste0("No solution found in ", counter, " steps :/"))
  out <- tibble(status = this_status, n_swaps_done = counter)
  return(out)
}

Let's test it.

get_swap_status(seed = 12345)
## Cost is $315.88.
## No optimal solution found :'(
## We've got a lot of CARROTS. 100 servings of them.
## Cost is $311.39.
## No optimal solution found :'(
## Cost is $292.1.
## No optimal solution found :'(
## Cost is $216.46.
## No optimal solution found :'(
## No solution found in 3 steps :/
## # A tibble: 1 x 2
##   status n_swaps_done
##    <int>        <dbl>
## 1      1            3

So for this particular random menu that was created, we can see whether a solution was found and how many swaps we did to get there.

Now we can do the same for a spectrum of minimum portion sizes with a fixed max number of swaps we're willing to do. Like we did with simulate_spectrum(), we split a minimum portion size (from to to) into n_intervals and do n_sims at each interval, recording the number of swaps it took to solve it. Our return tibble this time includes the minimum portion size we were allowed, the swap status (0 for good, 1 for bad), and the number of swaps we had to do

simulate_swap_spectrum <- function(n_intervals = 10, n_sims = 2, max_n_swaps = 3, from = -1, to = 1,
                                   seed = NULL, verbose = FALSE, ...) {
  
  interval <- (to - from) / n_intervals
  spectrum <- seq(from = from, to = to, by = interval) %>% rep(n_sims) %>% sort()
  
  if (!is.null(seed)) { set.seed(seed) }
  seeds <- sample(1:length(spectrum), size = length(spectrum), replace = FALSE)
  
  out_spectrum <- tibble(min_amount = spectrum)
  out_status <- tibble(status = vector(length = length(spectrum)), 
                       n_swaps_done = vector(length = length(spectrum)))
  
  for (i in seq_along(spectrum)) {
    this_status_df <- get_swap_status(seed = seeds[i], min_food_amount = spectrum[i], max_n_swaps = max_n_swaps, verbose = verbose)
    if (!is.integer(this_status_df$status)) {
      this_status_df$status <- integer(0)     # If we don't get an integer value back, make it NA
    }
    out_status[i, ] <- this_status_df
  }
  
  out <- bind_cols(out_spectrum, out_status)
  
  return(out)
}
simmed_swaps <- simulate_swap_spectrum(n_intervals = 20, n_sims = 5, max_n_swaps = 4, seed = 2018,
                                       from = -1, to = 2)
simmed_swaps %>% kable(format = "html")

Summarising a Spectrum

Let's get some summary values out of our spectra. summarise_status_spectrum() allows us to summarise either a vanilla status spectrum or a swap spectrum. For a given portion size, we'll find the proportion of those menus that are we were able to solve.

If we've allowed swapping, we'll find the average number of swaps we did at a given portion size.

summarise_status_spectrum <- function(spec) {
  
  # If this was a product of simulate_spectrum()
  if (!"n_swaps_done" %in% names(spec)){
    spec_summary <- spec %>% 
      group_by(min_amount) %>% 
      summarise(
        sol_prop = mean(status)
      )
    
    # If this was a product of simulate_swap_spectrum()
  } else {
    spec_summary <- spec %>% 
      group_by(min_amount) %>% 
      summarise(
        sol_prop = mean(status),
        mean_n_swaps_done = mean(n_swaps_done)
      )
  }
  
  return(spec_summary)
}

Let's summarise our vanilla spectrum. sol_prop refers to the proportion of the recipes that we were able to solve.

(status_spectrum_summary <- summarise_status_spectrum(status_spectrum))
## # A tibble: 11 x 2
##    min_amount sol_prop
##         <dbl>    <dbl>
##  1       -1.0      0.0
##  2       -0.8      0.0
##  3       -0.6      0.0
##  4       -0.4      0.5
##  5       -0.2      0.5
##  6        0.0      0.0
##  7        0.2      0.0
##  8        0.4      0.5
##  9        0.6      1.0
## 10        0.8      1.0
## 11        1.0      1.0

Now the fun part: visualizing the curve of minimum allowed portion size per food compared to the proportion of menus that were solvable at that portion size.

ggplot() +
  geom_smooth(data = status_spectrum, aes(min_amount, 1 - status),
              se = FALSE) +
  geom_point(data = status_spectrum_summary, aes(min_amount, 1 - sol_prop)) +
  theme_minimal() +
  ggtitle("Curve of portion size vs. solvability") +
  labs(x = "Minimum portion size", y = "Proportion of solutions") +
  ylim(0, 1) 

So we don't have a linear relationship between portion size and the proportion of menus that are solvable at that portion size.

We can do the same for our swap spectrum.

simmed_swaps_summary <- summarise_status_spectrum(simmed_swaps)
# Plot min portion size vs. whether we solved it or not
ggplot() +
  geom_smooth(data = simmed_swaps, aes(min_amount, 1 - status),
              se = FALSE) +
  geom_point(data = simmed_swaps_summary, aes(min_amount, 1 - sol_prop, colour = factor(mean_n_swaps_done))) +
  # facet_wrap( ~ n_swaps_done) +
  theme_minimal() +
  ggtitle("Curve of portion size vs. solvability") +
  labs(x = "Minimum portion size", y = "Proportion of solutions") +
  ylim(0, 1) 
## Warning: Removed 31 rows containing missing values (geom_smooth).




Scraping

I joked with my co-data scientist at Earlybird, Boaz Reisman, that this project so far could fairly be called "Eat, Pray, Barf." The menus we generate start off random and that's bad enough -- then, once we change up portion sizes, the menus only get less appetizing.

I figured the best way to decrease the barf factor was to look through how real menus are structured and try to suss out general patterns or rules in them. For instance, maybe we could learn that usually more than 1/3 of a dish should be dairy, or pork and apples tend to go well together.

I thought allrecipes.com would be likely to live up to its name and provide a good amount of data to work with. After a bit of poking a few recipes to try to discern if there was a pattern in how Allrecipes structures its URLs, I found that that all the recipe URLs followed this basic structure: http://allrecipes.com/recipe/<ID>/<NAME-OF-RECIPE>/. Omitting the <NAME-OF-RECIPE> parameter seemed to be fine in all cases; http://allrecipes.com/recipe/<ID> would redirect you to http://allrecipes.com/recipe/<ID>/<NAME-OF-RECIPE>/.

I couldn't figure out much of a pattern behind IDs save that they were always all digits and appeared to usually be between 10000 and 200000. (There's probably some pattern I'm missing here but this was good enough to start off with.)

So we know our base URL is going to be "http://allrecipes.com/recipe/".

base_url <- "http://allrecipes.com/recipe/"

Then we need to attach IDs to it, so for instance

grab_urls <- function(base_url, id) {
  id <- as.character(id)
  recipe_url <- str_c(base_url, id)
  return(recipe_url)
}

(urls <- grab_urls(base_url, 244940:244950))
##  [1] "http://allrecipes.com/recipe/244940"
##  [2] "http://allrecipes.com/recipe/244941"
##  [3] "http://allrecipes.com/recipe/244942"
##  [4] "http://allrecipes.com/recipe/244943"
##  [5] "http://allrecipes.com/recipe/244944"
##  [6] "http://allrecipes.com/recipe/244945"
##  [7] "http://allrecipes.com/recipe/244946"
##  [8] "http://allrecipes.com/recipe/244947"
##  [9] "http://allrecipes.com/recipe/244948"
## [10] "http://allrecipes.com/recipe/244949"
## [11] "http://allrecipes.com/recipe/244950"

Now that we've got URLs to scrape, we'll need to do the actual scraping.

Since we're appending some random numbers to the end of our base URL, there's a good chance some of those pages won't exist. We want a helper function that can try to read HTML on a page if it exists, and if the page doesn't exist, tell us without erroring out and exiting our loop. purrr::possibly() will let us do that. It provides a sort of try-catch set up where we try to read_url() but if we can't, return "Bad URL" and go on to the next URL.

read_url <- function(url) {
  page <- read_html(url)
}
try_read <- possibly(read_url, otherwise = "Bad URL", quiet = TRUE)

For example,

try_read("foo")
## [1] "Bad URL"

read_html() from the xml2 package will return us the raw HTML for a given page. We're only interested in the recipe portion of that, so using the Chrome inspector or the SelectorGadget Chrome extension we can figure out what the CSS tag is of the content itself.

The recipe's name gets the CSS class .recipe-summary__h1 and the content gets .checkList__line. So, we'll pluck everything tagged with those two classes using html_nodes() and return text we can use with html_text().

get_recipe_name <- function(page) {
  recipe_name <- page %>% 
    html_nodes(".recipe-summary__h1") %>% 
    html_text() 
  return(recipe_name)
}

Let's test that out on our fourth URL.

urls[4] %>% try_read() %>% get_recipe_name()
## [1] "Banana, Orange, and Ginger Smoothie"

We'll need an extra couple steps when it comes to recipe content to pare out all the stray garbage left over like \n new lines etc.

get_recipe_content <- function(page) {
  recipe <- page %>% 
    html_nodes(".checkList__line") %>% 
    html_text() %>% 
    str_replace_all("ADVERTISEMENT", "") %>% 
    str_replace_all("\n", "") %>% 
    str_replace_all("\r", "") %>% 
    str_replace_all("Add all ingredients to list", "")
  return(recipe)
}

And the content:

urls[4] %>% try_read() %>% get_recipe_content()
## [1] "                                            1 orange, peeled                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "                                  
## [2] "                                            1/2 banana                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "                                        
## [3] "                                            3 ice cubes                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "                                       
## [4] "                                            2 teaspoons honey                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "                                 
## [5] "                                            1/2 teaspoon grated fresh ginger root, or to taste                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "
## [6] "                                            1/2 cup plain yogurt                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    "                              
## [7] "                                                                                "                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  
## [8] "                                                                                                                                                                                                                                                                                                                                                                                                                                                        "                                                                                                                                                                                                          
## [9] "                                                "

Cool, so we've got three functions now, one for reading the content from a URL and turning it into a page and two for taking that page and grabbing the parts of it that we want. We'll use those functions in get_recipes() which will take a vector of URLs and return us a list of recipes. We also include parameters for how long to wait in between requests (sleep) so as to avoid getting booted from allrecipes.com and whether we want the "Bad URL"s included in our results list or not. If verbose is TRUE we'll get a message of count of the number of 404s we had and the number of duped recipes.

Note on dupes

Dupes come up because multiple IDs can point to the same recipe which means that two different URLs could resolve to the same page. I figured there were two routes we could go to see whether a recipe is a dupe or not; one, just go off of the recipe name or two, go off of the recipe name and content. By going off of the name, we don't go through the trouble of pulling in duped recipe content if we think we've got a dupe; we just skip it. Going off of content and checking whether the recipe content exists in our list so far would be safer (we'd only skip the recipes that we definitely already have), but slower because we have to both get_recipe_name() and get_recipe_content(). I went with the faster way; in get_recipes() we just check the recipe name we're on against all the recipe names in our list with if (!recipe_name %in% names(out)).

get_recipes <- function(urls, sleep = 5, verbose = TRUE, append_bad_URLs = TRUE) {
  bad_url_counter <- 0
  duped_recipe_counter <- 0
  
  out <- NULL       # In this case we don't know how long our list will be 

  
  for (url in urls) {
    Sys.sleep(sleep)    # Sleep in between requests to avoid 429 (too many requests)
    recipe_page <- try_read(url)
  
    if (recipe_page == "Bad URL" ||
       (!class(recipe_page) %in% c("xml_document", "xml_node"))) { 
      recipe_list <- recipe_page    # If we've got a bad URL, recipe_df will be "Bad URL" because of the otherwise clause
      bad_url_counter <- bad_url_counter + 1
      
      if (append_bad_URLs == TRUE) { out <- append(out, recipe_list) }

    } else {
      recipe_name <- get_recipe_name(recipe_page)
      
      if (!recipe_name %in% names(out)) {
        
        if (verbose == TRUE) { message(recipe_name) }
      
        recipe <- recipe_page %>% 
          get_recipe_content() %>% 
          map(remove_whitespace) %>% as_vector()
        
        recipe_list <- list(tmp_name = recipe) %>% as_tibble()  
        names(recipe_list) <- recipe_name
        
        out <- append(out, recipe_list)
        
      } else {
        duped_recipe_counter <- duped_recipe_counter + 1
        if (verbose == TRUE) {
          message("Skipping recipe we already have")
        }
      }
    }
  }
  if (verbose == TRUE) { 
    message(paste0("Number bad URLs: ", bad_url_counter))
    message(paste0("Number duped recipes: ", duped_recipe_counter))
  }
  
  return(out)
}

Let's give it a shot with a couple URLs.

(a_couple_recipes <- get_recipes(urls[4:5]))
## Banana, Orange, and Ginger Smoothie
## Alabama-Style White Barbecue Sauce
## Number bad URLs: 0
## Number duped recipes: 0
## [[1]]
## [1] FALSE
## 
## [[2]]
## [1] FALSE
## 
## $`Banana, Orange, and Ginger Smoothie`
## [1] "1 orange, peeled"                                  
## [2] "1/2 banana"                                        
## [3] "3 ice cubes"                                       
## [4] "2 teaspoons honey"                                 
## [5] "1/2 teaspoon grated fresh ginger root, or to taste"
## [6] "1/2 cup plain yogurt"                              
## 
## $`Alabama-Style White Barbecue Sauce`
## [1] "2 cups mayonnaise"                          
## [2] "1/2 cup apple cider vinegar"                
## [3] "1/4 cup prepared extra-hot horseradish"     
## [4] "2 tablespoons fresh lemon juice"            
## [5] "1 1/2 teaspoons freshly ground black pepper"
## [6] "2 teaspoons prepared yellow mustard"        
## [7] "1 teaspoon kosher salt"                     
## [8] "1/2 teaspoon cayenne pepper"                
## [9] "1/4 teaspoon garlic powder"

Now we've got a list of named recipes with one row per ingredient. Next step is tidying. We want to put this list of recipes into dataframe format with one observation per row and one variable per column. Our rows will contain items in the recipe content, each of which we'll associate with the recipe's name.

dfize <- function(lst, remove_bad_urls = TRUE) {

  df <- NULL
  if (remove_bad_urls == TRUE) {
    lst <- lst[!lst == "Bad URL"]
  }

  for (i in seq_along(lst)) {
    this_df <- lst[i] %>% as_tibble()
    recipe_name <- names(lst[i])
    names(this_df) <- "ingredients"
    this_df <- this_df %>% 
      mutate(recipe_name = recipe_name)
    df <- df %>% bind_rows(this_df)
  }
  return(df)
}
a_couple_recipes_df <- dfize(a_couple_recipes)
a_couple_recipes_df %>% kable(format = "html")
ingredients recipe_name
1 (12 inch) pre-baked pizza crust Johnsonville® Three Cheese Italian Style Chicken Sausage Skillet Pizza
1 1/2 cups shredded mozzarella cheese Johnsonville® Three Cheese Italian Style Chicken Sausage Skillet Pizza
1 (14 ounce) jar pizza sauce Johnsonville® Three Cheese Italian Style Chicken Sausage Skillet Pizza
1 (12 ounce) package Johnsonville® Three Cheese Italian Style Chicken Sausage, sliced Johnsonville® Three Cheese Italian Style Chicken Sausage Skillet Pizza
1 (3.5 ounce) package sliced pepperoni Johnsonville® Three Cheese Italian Style Chicken Sausage Skillet Pizza
3 teaspoons peanut oil, divided Ginger Fried Rice
4 large eggs, beaten Ginger Fried Rice
1 bunch scallions, chopped Ginger Fried Rice
2 tablespoons minced fresh ginger Ginger Fried Rice
3 cups cold cooked long-grain brown rice Ginger Fried Rice
1 cup frozen peas Ginger Fried Rice
1 cup mung bean sprouts (see Note) Ginger Fried Rice
3 tablespoons prepared stir-fry or oyster sauce Ginger Fried Rice

Great, so we've got a tidy dataframe that we can start to get some useful data out of.

One of the goals here is to see what portion of a menu tends to be devoted to, say, meat or spices or a word that appears in the recipe name etc. In order to answer that, we'll need to extract portion names and portion sizes from the text. That wouldn't be pretty simple with a fixed list of portion names ("gram", "lb") if portion sizes were always just a single number.

But, as it happens, portion sizes don't usually consist of just one number. There are a few hurdles:

  1. Complex fractions
  • 2 1/3 cups of flour should become: 2.3333 cups of flour
  1. Multiple items of the same item
  • 4 (12oz) bottles of beer should become: 48 oz of beer
  1. Ranges
  • 6-7 tomatoes should become: 6.5 tomatoes

Here is a fake recipe to illustrate some of those cases. (Certainly falls into Eat, Pray, Barf territory.)

some_recipes_tester %>% kable(format = "html")
ingredients
1.2 ounces or maybe pounds of something with a decimal
3 (14 ounce) cans o' beef broth
around 4 or 5 eels
5-6 cans spam
11 - 46 tbsp of sugar
1/3 to 1/2 of a ham
5 1/2 pounds of apples
4g cinnamon
about 17 fluid ounces of wine
4-5 cans of 1/2 caf coffee
3 7oz figs with 1/3 rind

Rather than start doing something conditional random field-level smart to get around these problems, to start off I started writing a few rules of thumb.

We'll worry first about how we find and extract numbers and next about how we'll add, multiply, or average them as necessary.

Extracting Numbers

We'll need a few regexes to extract our numbers.

portions_reg will match any digit even if it contains a decimal or a slash in it, which will be important for capturing complex fractions.

multiplier_reg covers all cases of numbers that might need to be multiplied in the Allrecipes data, because these are always separated by " (", whereas multiplier_reg_looser is a more loosely-defined case matching numbers separated just by " ".

# Match any number, even if it has a decimal or slash in it
portions_reg <- "[[:digit:]]+\\.*[[:digit:]]*+\\/*[[:digit:]]*"

# Match numbers separated by " (" as in "3 (5 ounce) cans of broth" for multiplying
multiplier_reg <- "[[:digit:]]+ \\(+[[:digit:]]"   

# Match numbers separated by " "
multiplier_reg_looser <- "[0-9]+\ +[0-9]"

Now the multiplier_reg regexes will allow us to detect that we've got something that needs to be multiplied, like "4 (12 oz) hams" or a fraction like "1 2/3 pound of butter". If we do, then we'll multiply or add those numbers as appropriate. The only_mult_after_paren parameter is something I put in that is specific to Allrecipes. On Allrecipes, it seems that if we do have multiples, they'll always be of the form "number_of_things (quantity_of_single_thing)". There are always parentheses around quantity_of_single_thing. If we're only using Allrecipes data, that gives us some more security that we're only multiplying quantities that actually should be multiplied. If we want to make this extensible in the future we'd want to set only_mult_after_paren to FALSE to account for cases like "7 4oz cans of broth".

We use str_extract() to check that our regexes are grabbing the parts of a string that we'll need to do computation on.

str_extract_all("3 1/4 lb patties", portions_reg)
## [[1]]
## [1] "3"   "1/4"

And check that multiplier_reg

str_extract_all("3 (4 pound patties) for grilling", multiplier_reg)
## [[1]]
## [1] "3 (4"

We'll clean that up by passing it to portions_reg to just grab the numbers:

str_extract_all("3 (4 pound patties) for grilling", multiplier_reg) %>% str_extract_all(portions_reg)
## [[1]]
## [1] "3" "4"

Finally, let's make sure that our stricter multiplier regex doesn't want to multiply something shouldn't be multiplied.

str_extract_all("3 or 4 lb patties", multiplier_reg)
## [[1]]
## character(0)

Okay, now to the multiplying and adding.

First, let's consider complex fractions. Off the bat, we know we'll need a way to turn a single fraction into a decimal form. We keep them as.character for now and turn them into numeric later down the pipe.

frac_to_dec <- function(e) {
  if (length(e) == 0) {    # If NA because there are no numbers, make the portion 0
    out <- 0
  } else {
    out <- parse(text = e) %>% eval() %>% as.character()
  }
  return(out)
}

eval(), which is what does the work inside frac_to_dec) only only evaluates the last string in a vector, not multiple, so as a workaround I put it into a helper that will turn all fractions into decimal strings:

map_frac_to_dec <- function(e) {
  out <- NULL
  for (i in e) {
    out <- e %>% map_chr(frac_to_dec)
  }
  return(out)
}

For example:

map_frac_to_dec(c("1/2", "1/8", "1/3"))
## [1] "0.5"               "0.125"             "0.333333333333333"

Cool, so for a given ingredient we'll need to look for numbers that are occur next to other numbers, and then and add complex fractions and multiply multiples.

If we've got two numbers next to each other and the second number evaluates to a decimal less than 1, we've got a complex fraction. For example, if we're extracting digits and turning all fractions among them into decimals if we consider "4 1/2 loaves of bread" we'd end up with "4" and "0.5". We know 0.5 is less than 1, so we've got a complex fraction on our hands. We need to add 4 + 0.5 to end up with 4.5 loaves of bread.

It's true that this function doesn't address the issue of having both a complex fraction and multiples in a recipe. That would look like "3 (2 1/4 inch) blocks of cheese." I haven't run into that issue too much but it certainly could use a workaround.

multiply_or_add_portions <- function(e) {
  if (length(e) == 0) {
    e <- 0    
  } else if (length(e) > 1) {
    if (e[2] < 1) {  # If our second element is a fraction, we know this is a complex fraction so we add the two
      e <- e[1:2] %>% reduce(`+`)
    } else {   # Otherwise, we multiply them
      e <- e[1:2] %>% reduce(`*`)
    }   
  }
  return(e)
}
multiply_or_add_portions(c(4, 0.5))
## [1] 4.5
multiply_or_add_portions(c(4, 5))
## [1] 20

This function will allow us to add a new column to our dataframe called mult_add_portion. If we've done any multiplying or adding of numbers, we'll have a value greater than 0 there, and 0 otherwise.

get_mult_add_portion <- function(e, only_mult_after_paren = FALSE) {
  if ((str_detect(e, multiplier_reg) == TRUE | str_detect(e, multiplier_reg_looser) == TRUE)
      & only_mult_after_paren == FALSE) {  # If either matches and we don't care about where there's a parenthesis there or not
      if (str_detect(e, multiplier_reg) == TRUE) {
        out <- e %>% str_extract_all(portions_reg) %>% 
          map(map_frac_to_dec) %>%   
          map(as.numeric) %>% 
          map_dbl(multiply_or_add_portions) %>%   
          round(digits = 2)
    } else {    # If we do care, and we have a parenthesis
      out <- e %>% str_extract_all(portions_reg) %>% 
        map(map_frac_to_dec) %>%   
        map(as.numeric) %>% 
        map_dbl(multiply_or_add_portions) %>%   
        round(digits = 2)
    }
  } else {
    out <- 0
  }
  return(out)
}
get_mult_add_portion("4 1/2 steaks") 
## [1] 4.5
get_mult_add_portion("4 (5 lb melons)") 
## [1] 20

Ranges

Finally, let's deal with ranges. If two numbers are separated by an "or" or a "-" like "4-5 teaspoons of sugar" we know that this is a range. We'll take the average of those two numbers.

We'll add a new column to our dataframe called range_portion for the result of any range calculations. If we don't have a range, just like mult_add_portion, we set this value to 0.

to_reg <- "([0-9])(( to ))(([0-9]))"
or_reg <- "([0-9])(( or ))(([0-9]))"
dash_reg_1 <- "([0-9])((-))(([0-9]))"
dash_reg_2 <- "([0-9])(( - ))(([0-9]))"

First, a couple helpers. If two numbers are separated by an "or" or a "-" we know that this is a range, e.g., 4-5 teaspoons of sugar.

determine_if_range <- function(ingredients) {
  if (str_detect(ingredients, pattern = to_reg) | 
      str_detect(ingredients, pattern = or_reg) |
      str_detect(ingredients, pattern = dash_reg_1) |
      str_detect(ingredients, pattern = dash_reg_2)) {
    contains_range <- TRUE
  } else {
    contains_range <- FALSE
  }
  return(contains_range)
}

And, we'll want to be able to get the mean of the first two elements in a numeric vector.

get_portion_means <- function(e) {
  if (length(e) == 0) {
    e <- 0    # NA to 0
  } else if (length(e) > 1) {
      e <- mean(e[1:2])
  }
  return(e)
}
get_ranges <- function(e) {
  
  if (determine_if_range(e) == TRUE) {
    out <- str_extract_all(e, portions_reg) %>%  
      
      map(str_split, pattern = " to ", simplify = FALSE) %>%  
      map(str_split, pattern = " - ", simplify = FALSE) %>%  
      map(str_split, pattern = "-", simplify = FALSE) %>%
      
      map(map_frac_to_dec) %>%
      map(as.numeric) %>% 
      map_dbl(get_portion_means) %>% round(digits = 2)
    
  } else {
    out <- 0
  }
  return(out)
}

Let's make sure we get the average.

get_ranges("7 to 21 peaches")
## [1] 14

At the end of the day, we want to end up with a single number describing how much of our recipe item we want. So, let's put all that together into one function. Either range_portion or mult_add_portion will always be 0, so we add them together to get our final portion size. If we neither need to get a range nor multiply or add numbers, we'll just take whatever the first number is in there.

get_portion_values <- function(df, only_mult_after_paren = FALSE) {
  df <- df %>% 
    mutate(
      range_portion = map_dbl(ingredients, get_ranges),
      mult_add_portion = map_dbl(ingredients, get_mult_add_portion, only_mult_after_paren = only_mult_after_paren),
      portion = ifelse(range_portion == 0 & mult_add_portion == 0,
                       str_extract_all(ingredients, portions_reg) %>%
                         map(map_frac_to_dec) %>%
                         map(as.numeric) %>%
                         map_dbl(first),
                       range_portion + mult_add_portion)   # Otherwise, take either the range or the multiplied value
    )
  return(df)
}

Let's see what that looks like in practice.

some_recipes_tester %>% get_portion_values() %>% kable(format = "html")
ingredients range_portion mult_add_portion portion
1.2 ounces or maybe pounds of something with a decimal 0.00 0.0 1.20
3 (14 ounce) cans o' beef broth 0.00 42.0 42.00
around 4 or 5 eels 4.50 0.0 4.50
5-6 cans spam 5.50 0.0 5.50
11 - 46 tbsp of sugar 28.50 0.0 28.50
1/3 to 1/2 of a ham 0.42 0.0 0.42
5 1/2 pounds of apples 0.00 5.5 5.50
4g cinnamon 0.00 0.0 4.00
about 17 fluid ounces of wine 0.00 0.0 17.00
4-5 cans of 1/2 caf coffee 4.50 0.0 4.50
3 7oz figs with 1/3 rind 0.00 21.0 21.00

Looks pretty solid.

Extracting Measurement Units

Now onto easier waters: portion names. You can check out /scripts/scrape/get_measurement_types.R if you're interested in the steps I took to find some usual portion names and create an abbreviation dictionary, abbrev_dict. What we also do there is create measures_collapsed which is a single vector of all portion names separated by "|" so we can find all the portion names that might occur in a given item.

measures_collapsed
## [1] "[[:digit:]]oz |[[:digit:]]pt |[[:digit:]]lb |[[:digit:]]kg |[[:digit:]]g |[[:digit:]]l |[[:digit:]]dl |[[:digit:]]ml |[[:digit:]] oz |[[:digit:]] pt |[[:digit:]] lb |[[:digit:]] kg |[[:digit:]] g |[[:digit:]] l |[[:digit:]] dl |[[:digit:]] ml |ounce|pint|pound|kilogram|gram|liter|deciliter|milliliter"

Then if there are multiple portions that match, we'll take the last one.

We'll also add approximate to our dataframe which is just a boolean value indicating whether this item is exact or approximate. If the item contains one of approximate (about ,around ,as desired ,as needed ,optional ,or so ,to taste) then we give it a TRUE.

str_detect("8 or so cloves of garlic", approximate)
## [1] TRUE
str_detect("8 cloves of garlic", approximate)
## [1] FALSE
get_portion_text <- function(df) {
  
  df <- df %>% 
    mutate(
      raw_portion_num = str_extract_all(ingredients, portions_reg, simplify = FALSE) %>%   # Extract the raw portion numbers,
        map_chr(str_c, collapse = ", ", default = ""),   # separating by comma if multiple
      
      portion_name = str_extract_all(ingredients, measures_collapsed) %>%
        map(nix_nas) %>%  
        str_extract_all("[a-z]+") %>% 
        map(nix_nas) %>%   # Get rid of numbers
        map_chr(last),       # If there are multiple arguments that match, grab the last one

      approximate = str_detect(ingredients, approximate)
    )
  return(df)
}

Last thing for us for now on this subject (though there's a lot more to do here!) will be to add abbreviations. This will let us standardize things like "ounces" and "oz" which actually refer to the same thing.

All add_abbrevs() will do is let us mutate our dataframe with a new column for the abbreviation of our portion size, if we've got a recognized portion size.

add_abbrevs <- function(df) {

  out <- vector(length = nrow(df))
  for (i in seq_along(out)) {
    if (df$portion_name[i] %in% abbrev_dict$name) {
      out[i] <- abbrev_dict[which(abbrev_dict$name == df$portion_name[i]), ]$key
      
    } else {
      out[i] <- df$portion_name[i]
    }
  }
  
  out <- df %>% bind_cols(list(portion_abbrev = out) %>% as_tibble())
  return(out)
}
tibble(ingredients = "10 pounds salt, or to taste") %>% 
  get_portion_text() %>% add_abbrevs() %>% kable(format = "html")
ingredients raw_portion_num portion_name approximate portion_abbrev
10 pounds salt, or to taste 10 pound TRUE lb

All together now. Get the portion text and values. If we only want our best guess as to the portion size, that is, final_portion_size, we'll chuck range_portion and mult_add_portion.

get_portions <- function(df, add_abbrevs = FALSE, pare_portion_info = FALSE) {
  df %<>% get_portion_text() 
  if (add_abbrevs == TRUE) {
    df %<>% add_abbrevs()
  }
  df %<>% get_portion_values()
  if (pare_portion_info == TRUE) {
    df %<>% select(-range_portion, -mult_add_portion)
  }
  return(df)
}
some_recipes_tester %>%
  get_portions(pare_portion_info = TRUE) %>% 
  add_abbrevs() %>% kable(format = "html")
ingredients raw_portion_num portion_name approximate portion portion_abbrev
1.2 ounces or maybe pounds of something with a decimal 1.2 pound FALSE 1.20 lb
3 (14 ounce) cans o' beef broth 3, 14 ounce FALSE 42.00 oz
around 4 or 5 eels 4, 5 TRUE 4.50
5-6 cans spam 5, 6 FALSE 5.50
11 - 46 tbsp of sugar 11, 46 tbsp FALSE 28.50 tbsp
1/3 to 1/2 of a ham 1/3, 1/2 FALSE 0.42
5 1/2 pounds of apples 5, 1/2 pound FALSE 5.50 lb
4g cinnamon 4 g FALSE 4.00 g
about 17 fluid ounces of wine 17 ounce TRUE 17.00 oz
4-5 cans of 1/2 caf coffee 4, 5, 1/2 FALSE 4.50
3 7oz figs with 1/3 rind 3, 7, 1/3 oz FALSE 21.00 oz

We've got some units! Next step will be to convert all units into grams, so that we have them all in a standardized format.

Converting to Grams

Rather than rolling our own conversion dictionary, let's turn to the measurements package that sports the conv_unit() function for going from one unit to another. For example, converting 12 inches to centimeters, we get:

conv_unit(12, "inch", "cm")
## [1] 30.48

Let's see how that'll work with our data. Let's take our 14 oz to grams.

conv_unit(more_recipes_df[3, ]$portion, more_recipes_df[3, ]$portion_abbrev, "g")
## [1] 396.8933

Let's see which of our units conv_unit() can successfully convert out of the box.

We'll set up exception handling so that conv_unit() gives us an NA rather than an error if it encounters a value it can't convert properly.

try_conv <- possibly(conv_unit, otherwise = NA)

We'll mutate our abbreviation dictionary, adding a new column to convert to either grams in the case that our unit is a solid or milliliters if it's a liquid. These have a 1-to-1 conversion (1g = 1ml) so we'll take whichever one of these is not a missing value and put that in our converted column.

We'll use a sample value of 10 for everything.

test_abbrev_dict_conv <- function(dict, key_col, val = 10) {
  
  quo_col <- enquo(key_col)
  
  out <- dict %>% 
    rowwise() %>% 
    mutate(
      converted_g = try_conv(val, !!quo_col, "g"),
      converted_ml = try_conv(val, !!quo_col, "ml"),
      converted = case_when(
        !is.na(converted_g) ~ converted_g,
        !is.na(converted_ml) ~ converted_ml
      )
    )
  
  return(out)
}
test_abbrev_dict_conv(abbrev_dict, key)
## Source: local data frame [12 x 5]
## Groups: <by row>
## 
## # A tibble: 12 x 5
##           name      key converted_g converted_ml  converted
##          <chr>    <chr>       <dbl>        <dbl>      <dbl>
##  1       ounce       oz    283.4952           NA   283.4952
##  2        pint       pt          NA           NA         NA
##  3       pound       lb          NA           NA         NA
##  4    kilogram       kg  10000.0000           NA 10000.0000
##  5        gram        g     10.0000           NA    10.0000
##  6       liter        l          NA        10000 10000.0000
##  7   deciliter       dl          NA         1000  1000.0000
##  8  milliliter       ml          NA           10    10.0000
##  9  tablespoon     tbsp          NA           NA         NA
## 10    teaspoon      tsp          NA           NA         NA
## 11         cup      cup          NA           NA         NA
## 12 fluid ounce fluid oz          NA           NA         NA

What proportion of the portion abbreviations are we able to to convert to grams off the bat?

converted_units <- test_abbrev_dict_conv(abbrev_dict, key)
length(converted_units$converted[!is.na(converted_units$converted)]) / length(converted_units$converted)
## [1] 0.5

We can take a look at the units that measurements provides conversions for to see if we'll need to go elsewhere to do the conversion math ourselves.

conv_unit_options$volume
##  [1] "ul"        "ml"        "dl"        "l"         "cm3"      
##  [6] "dm3"       "m3"        "km3"       "us_tsp"    "us_tbsp"  
## [11] "us_oz"     "us_cup"    "us_pint"   "us_quart"  "us_gal"   
## [16] "inch3"     "ft3"       "mi3"       "imp_tsp"   "imp_tbsp" 
## [21] "imp_oz"    "imp_cup"   "imp_pint"  "imp_quart" "imp_gal"

This explains why pint, cup, etc. weren't convertable. It looks like we need to put the prefix "us_" before some of our units. We'll create a new accepted column of abbrev_units that provides the convertable

to_usize <- c("tsp", "tbsp", "cup", "pint")  

accepted <- c("oz", "pint", "lbs", "kg", "g", "l", "dl", "ml", "tbsp", "tsp", "cup", "oz")
accepted[which(accepted %in% to_usize)] <- 
  stringr::str_c("us_", accepted[which(accepted %in% to_usize)])

# cbind this to our dictionary 
abbrev_dict_w_accepted <- abbrev_dict %>% bind_cols(accepted = accepted)

What percentage of units are we able to convert now?

test_abbrev_dict_conv(abbrev_dict_w_accepted, accepted)
## Source: local data frame [12 x 6]
## Groups: <by row>
## 
## # A tibble: 12 x 6
##           name      key accepted converted_g converted_ml   converted
##          <chr>    <chr>    <chr>       <dbl>        <dbl>       <dbl>
##  1       ounce       oz       oz    283.4952           NA   283.49523
##  2        pint       pt  us_pint          NA   4731.76473  4731.76473
##  3       pound       lb      lbs   4535.9243           NA  4535.92428
##  4    kilogram       kg       kg  10000.0000           NA 10000.00000
##  5        gram        g        g     10.0000           NA    10.00000
##  6       liter        l        l          NA  10000.00000 10000.00000
##  7   deciliter       dl       dl          NA   1000.00000  1000.00000
##  8  milliliter       ml       ml          NA     10.00000    10.00000
##  9  tablespoon     tbsp  us_tbsp          NA    147.86765   147.86765
## 10    teaspoon      tsp   us_tsp          NA     49.28922    49.28922
## 11         cup      cup   us_cup          NA   2365.88236  2365.88236
## 12 fluid ounce fluid oz       oz    283.4952           NA   283.49523

Looks like all of them! Good stuff.

Let's write a function to convert units for our real dataframe.

convert_units <- function(df, name_col = accepted, val_col = portion,
                          pare_down = TRUE, round_to = 2) {
  
  quo_name_col <- enquo(name_col)
  quo_val_col <- enquo(val_col)
  
  out <- df %>% 
    rowwise() %>% 
    mutate(
      converted_g = try_conv(!!quo_val_col, !!quo_name_col, "g") %>% round(digits = round_to),
      converted_ml = try_conv(!!quo_val_col, !!quo_name_col, "ml") %>% round(digits = round_to), 
      converted = case_when(
        !is.na(converted_g) ~ as.numeric(converted_g), 
        !is.na(converted_ml) ~ as.numeric(converted_ml), 
        is.na(converted_g) && is.na(converted_ml) ~ NA_real_ 
      )
    ) 
  
  if (pare_down == TRUE) {
    out <- out %>% 
      select(-converted_g, -converted_ml)
  }
  
  return(out)
}

Next let's add an accepted column onto our dataframe to get our units in the right format and run our function.

more_recipes_df %>% 
  left_join(abbrev_dict_w_accepted, by = c("portion_abbrev" = "key")) %>% 
  sample_n(30) %>%
  convert_units() %>% 
  kable(format = "html")
ingredients recipe_name raw_portion_num portion_name approximate range_portion mult_add_portion portion portion_abbrev name accepted converted
1 (14.5 ounce) can chicken broth African Chicken in Spicy Red Sauce 1, 14.5 ounce FALSE 0 14.50 14.5000000 oz ounce oz 411.06809
1 cup blueberries Fourth of July Salad 1 FALSE 0 0.00 1.0000000 NA NA NA
1/2 cup sweetened flaked coconut Baked Coconut Shrimp with Spicy Dipping Sauce 1/2 FALSE 0 0.00 0.5000000 NA NA NA
2 eggs Tuna Fish Patties 2 FALSE 0 0.00 2.0000000 NA NA NA
1 (.25 ounce) package yeast White Pizza with Mushrooms 1, 25 ounce FALSE 0 0.00 1.0000000 oz ounce oz 28.34952
1/4 teaspoon pepper Caramelized Pork Slices 1/4 FALSE 0 0.00 0.2500000 NA NA NA
1 teaspoon salt Seasoned Hamburger Mix 1 FALSE 0 0.00 1.0000000 NA NA NA
1 cup chopped baby spinach White Bean Tabbouleh 1 FALSE 0 0.00 1.0000000 NA NA NA
1/2 teaspoon sugar Baby Lettuces with Green Apple, Walnuts, and Dried Cranberries 1/2 FALSE 0 0.00 0.5000000 NA NA NA
1 cup milk Sweet Baked Onion Rings 1 FALSE 0 0.00 1.0000000 NA NA NA
2 cups water, or as needed - divided Bayrischer Schweinebraten (Bavarian Roast Pork) 2 TRUE 0 0.00 2.0000000 NA NA NA
1/2 teaspoon garlic salt Quinoa Crust Quiche 1/2 FALSE 0 0.00 0.5000000 NA NA NA
1 (15 ounce) can white or yellow corn, drained Black-Eyed Peas Salad 1, 15 ounce FALSE 0 15.00 15.0000000 oz ounce oz 425.24285
1 (1.25 ounce) package TACO BELL HOME ORIGINALS Taco Seasoning Mix Burrito Bake 1, 1.25 ounce FALSE 0 1.25 1.2500000 oz ounce oz 35.43690
2 cups cubed fully cooked ham Ham N Cheese Potato Bake 2 FALSE 0 0.00 2.0000000 NA NA NA
4 boneless lean pork loin chops Apples 'n' Onion Topped Chops 4 FALSE 0 0.00 4.0000000 NA NA NA
2/3 cup fat-free sour cream Mexican Beef and Corn Casserole from Country Crock® 2/3 FALSE 0 0.00 0.6666667 NA NA NA
1/2 cup chopped red onion Roast Beef Sandwich Roll 1/2 FALSE 0 0.00 0.5000000 NA NA NA
1 (200 g) package KRAFT Dinner Extra Creamy Macaroni and Cheese One-Pot Salsa Beef Skillet 1, 200 FALSE 0 200.00 200.0000000 NA NA NA
2 1/2 tablespoons extra-virgin olive oil Baby Lettuces with Green Apple, Walnuts, and Dried Cranberries 2, 1/2 FALSE 0 2.50 2.5000000 NA NA NA
1 cup salsa Veggie Fajitas 1 FALSE 0 0.00 1.0000000 NA NA NA
4 dried pasilla chilies Chicken Mole with Four Chiles 4 FALSE 0 0.00 4.0000000 NA NA NA
1/2 teaspoon dried oregano Extra-Virgin Olive Oil Dipping Sauce 1/2 FALSE 0 0.00 0.5000000 NA NA NA
1 tablespoon cold unsalted butter, cut into chunks Crust for Two 1 FALSE 0 0.00 1.0000000 NA NA NA
6 cloves garlic, minced Grilled Open-Face Feta and Veggie Sandwich 6 FALSE 0 0.00 6.0000000 NA NA NA
2 cups frozen corn One-Pot Salsa Beef Skillet 2 FALSE 0 0.00 2.0000000 NA NA NA
Salt and black pepper Easy Barbecued Spareribs FALSE 0 0.00 NA NA NA NA
2 eggs Carrot Pancake Cake 2 FALSE 0 0.00 2.0000000 NA NA NA
1 onion, chopped Claire's Curried Butternut Squash Soup 1 FALSE 0 0.00 1.0000000 NA NA NA
1 1/4 cups milk Broccoli Sausage Breakfast Bake 1, 1/4 FALSE 0 1.25 1.2500000 NA NA NA

All the data

Let's put it all together, scraping all of our URLs and doing all of the munging in order.

recipes_raw <- more_urls %>% get_recipes(sleep = 3)
recipes <- recipes_raw[!recipes_raw == "Bad URL"]

recipes_df <- recipes %>% 
  dfize() %>% 
  get_portions() %>% 
  add_abbrevs() %>% 
  left_join(abbrev_dict_w_accepted, by = c("portion_abbrev" = "key")) %>% 
  convert_units()
set.seed(1234)

recipes_df %>% 
  sample_n(20) %>%
  select(ingredients, recipe_name, portion, portion_abbrev, converted) %>% 
  kable(format = "html")
ingredients recipe_name portion portion_abbrev converted
1 (14.5 ounce) can diced tomatoes No Ordinary Meatloaf 14.5000000 oz 411.068085
2 eggs Blueberry Banana Nut Bread 2.0000000 NA
1/2 cup KRAFT LIGHT DONE RIGHT! Raspberry Vinaigrette Reduced Fat Dressing Chicken and Citrus Salad 0.5000000 cup 118.294118
2 tablespoons cold butter Blueberry Banana Nut Bread 2.0000000 tbsp 29.573530
1 cup superfine sugar Banana Coffee Cake with Pecans 1.0000000 cup 236.588236
1 onion, chopped Claire's Curried Butternut Squash Soup 1.0000000 NA
2/3 cup milk Cabbage Cakes 0.6666667 cup 157.725491
1/2 cup red salsa Chorizo, Potato and Green Chile Omelet 0.5000000 cup 118.294118
1/2 cup raisins Fudge Drops 0.5000000 cup 118.294118
1/4 cup brown sugar Mango Chicken with Greens 0.2500000 cup 59.147059
4 tablespoons Safeway SELECT Verdi Olive Oil Bistro Beef Salad 4.0000000 tbsp 59.147059
1/4 cup canola oil Hungarian Salad 0.2500000 cup 59.147059
4 ounces PHILADELPHIA Neufchatel Cheese, 1/3 Less Fat than Cream Cheese Creamy Two-Layer Pumpkin Pie 4.0000000 oz 113.398093
1 tablespoon milk Caramel Pear Crumble 1.0000000 tbsp 14.786765
3 tablespoons Worcestershire sauce Traveling Oven-Barbecued Baby Back Ribs 3.0000000 tbsp 44.360294
1 teaspoon crushed red pepper flakes Artichoke and Shrimp Linguine 1.0000000 tsp 4.928922
1 tablespoon freshly ground black pepper African Chicken in Spicy Red Sauce 1.0000000 tbsp 14.786765
2 cups sugar snap peas, trimmed Sugar Snap Salad 2.0000000 cup 473.176473
1 teaspoon canola oil Apples 'n' Onion Topped Chops 1.0000000 tsp 4.928922
1 (8 ounce) package KRAFT Mexican Style Shredded Four Cheese with a Touch of PHILADELPHIA, divided Chorizo, Potato and Green Chile Omelet 8.0000000 oz 226.796185

Recipe Text NLP Analysis

Now we'll do a quick text analysis of words in these freshly scraped recipes. The text we'll be pulling ngrams out of is the ingredients column, where all the action happens.

Let's take a sample of recipes and extract all the ngrams in those recipes to get into n_buckets bags of words. We want to take a sample of whole recipes and not just of ingredients because we'll want a full set of all the ingredients for a particular recipe; an ingredient isn't a meaningful unit for us outside the context of its recipe.

We'll use tidytext::unnest_tokens() to get sets of n_grams for each recipe. If n_grams is 1, that's every word in a recipe. If n_grams is 3, that's every instance of 3 words that occur in a row.

grab_words <- function(df, n_buckets = 50, n_grams = 1) {
  recipes_to_keep <- df %>% 
    distinct(recipe_name) %>% 
    sample_n(n_buckets) %>% 
    pull("recipe_name")
  
  df <- df %>% 
    filter(recipe_name %in% recipes_to_keep) %>% 
    group_by(recipe_name) %>% 
    mutate(ingredient_num = row_number()) %>% 
    ungroup() %>% 
    unnest_tokens(word, ingredients, token = "ngrams", n = n_grams) %>% 
    select(recipe_name, word, everything())
  
  return(df)
}

First, a glance at what our data look like before we slice it up in to unigrams:

recipes_df %>% 
  sample_n(10) %>% 
  select(ingredients, recipe_name, portion_abbrev, portion, converted) %>% 
  kable()

ingredients recipe_name portion_abbrev portion converted


12 slices thinly sliced country ham, halved Country Ham and Cheddar Mini Hand Pies 12.00 NA 3 teaspoons cornstarch Saucy Apricot Chicken tsp 3.00 14.786765 Cake: Ultimate Black Forest Cake NA NA 1 teaspoon lemon-pepper seasoning Seasoned Broccoli Spears tsp 1.00 4.928922 2 tablespoons chopped onion Macaroni 'n' Cheese for Two tbsp 2.00 29.573530 2 cups assorted wild mushrooms, chopped Broccoflower® Risotto with Wild Mushrooms cup 2.00 473.176473 1 avocado, peeled and pitted Avocado Lime Hummus 1.00 NA 1/2 cup golden raisins Creamy Asian Slaw cup 0.50 118.294118 1 teaspoon salt Sweet Jalapeno Cornbread tsp 1.00 4.928922 1/4 teaspoon cayenne pepper Pork Picante tsp 0.25 1.232230

Let's unspool this into one unigram (single word) per row:

unigrams <- grab_words(recipes_df)

unigrams[1:10, ] %>% kable()

recipe_name word raw_portion_num portion_name approximate range_portion mult_add_portion portion portion_abbrev name accepted converted ingredient_num


20-Minute Green Bean Chicken and Rice Dinner 1 1 tablespoon FALSE 0 0.00 1.00 tbsp tablespoon us_tbsp 14.78676 1 20-Minute Green Bean Chicken and Rice Dinner tablespoon 1 tablespoon FALSE 0 0.00 1.00 tbsp tablespoon us_tbsp 14.78676 1 20-Minute Green Bean Chicken and Rice Dinner cooking 1 tablespoon FALSE 0 0.00 1.00 tbsp tablespoon us_tbsp 14.78676 1 20-Minute Green Bean Chicken and Rice Dinner oil 1 tablespoon FALSE 0 0.00 1.00 tbsp tablespoon us_tbsp 14.78676 1 20-Minute Green Bean Chicken and Rice Dinner 1 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5 20-Minute Green Bean Chicken and Rice Dinner 10.75 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5 20-Minute Green Bean Chicken and Rice Dinner ounce 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5 20-Minute Green Bean Chicken and Rice Dinner can 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5 20-Minute Green Bean Chicken and Rice Dinner condensed 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5 20-Minute Green Bean Chicken and Rice Dinner cream 1, 10.75 ounce FALSE 0 10.75 10.75 oz ounce oz 304.75737 5

and do the same for bigrams. This replaces our ingredients column with a new word column

bigrams <- grab_words(recipes_df, n_grams = 2)

bigrams[1:10, ] %>% kable()

recipe_name word raw_portion_num portion_name approximate range_portion mult_add_portion portion portion_abbrev name accepted converted ingredient_num


Crust for Two 1 tablespoon 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two tablespoon cold 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two cold unsalted 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two unsalted butter 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two butter cut 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two cut into 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two into chunks 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 5 Crust for Two 1 tablespoon 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 6 Crust for Two tablespoon canola 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 6 Crust for Two canola or 1 tablespoon FALSE 0 0 1 tbsp tablespoon us_tbsp 14.78676 6

We'll focus on just the unigrams for now.

Cleaning unigrams

Now that we've got all our individual words, we want to do some cleaning to make sure we only keep words that are actually meaningful in our context. The tidytext package helpfully supplies a stop_words dataset that can be easily anti_joined to our unigrams. We'll do that in a minute.

Just like we don't want "the"s and "and"s muddying our recipe content waters, we also want to pare out any ngrams that are (or possibly contain) numbers. Here a quick helper for determining that by adding a column with a logical is_num. We could as easily have used a regex, but instead we'll use an approach of seeing if we can coerce a value to numeric. If we can, it's definitely a number.

find_nums <- function(df, col = word, add_contains_num = TRUE) {
  quo_col <- rlang::enquo(col)

  df <- df %>% dplyr::mutate(
    num = suppressWarnings(as.numeric(!!quo_col)),  
    is_num = dplyr::case_when(
      !is.na(num) ~ TRUE,
      is.na(num) ~ FALSE
    )
  )

  if (add_contains_num == TRUE) {
    df <- df %>%
      dplyr::mutate(
        contains_num = dplyr::case_when(
          grepl("\\d+", !!quo_col) ~ TRUE,
          !(grepl("\\d+", !!quo_col)) ~ FALSE
        )
      )
  }

  df <- df %>% dplyr::select(-num)

  return(df)
}

This works like so:

test_nums <- tibble(
  a = runif(3),
  b = janeaustenr::mansfieldpark[1:3]
)

test_nums %>% 
  find_nums(a)
## # A tibble: 3 x 4
##       a b              is_num contains_num
##   <dbl> <chr>          <lgl>  <lgl>       
## 1 0.456 MANSFIELD PARK TRUE   TRUE        
## 2 0.265 ""             TRUE   TRUE        
## 3 0.305 (1814)         TRUE   TRUE
test_nums %>% 
  find_nums(b)
## # A tibble: 3 x 4
##       a b              is_num contains_num
##   <dbl> <chr>          <lgl>  <lgl>       
## 1 0.456 MANSFIELD PARK FALSE  FALSE       
## 2 0.265 ""             FALSE  FALSE       
## 3 0.305 (1814)         FALSE  TRUE

For now we won't worry about whether a unigram contains a number (for instance, the A1 in A1 Sauce contains a number but I'd consider it fair game as an ingredient.)

But we will filter out anything in unigrams that is a straight up number.

unigrams <- unigrams %>%
  find_nums(add_contains_num = FALSE) %>%
  filter(is_num == FALSE) %>% 
  select(-is_num)

We also want to filter out anything that's a unit. We'll put all of our units, including the plurals, into a dataframe so it can be anti_joined on our words.

all_units <- c(units, abbrev_dict$name, abbrev_dict$key, "inch")
(all_units_df <- list(word = all_units) %>% as_tibble())
## # A tibble: 77 x 1
##    word       
##    <chr>      
##  1 cups       
##  2 cup        
##  3 tablespoons
##  4 tablespoon 
##  5 teaspoons  
##  6 teaspoon   
##  7 pounds     
##  8 pound      
##  9 ounces     
## 10 ounce      
## # ... with 67 more rows

Finding Patterns

Now looking at pairs of words within a recipe, which pairs tend to co-occur, i.e., have a higher frequency within the same recipe?

per_rec_freq <- unigrams %>% 
  anti_join(stop_words) %>% 
  anti_join(all_units_df) %>% 
  group_by(recipe_name) %>% 
  add_count(word) %>%    # Count of number of times this word appears in this recipe
  rename(n_this_rec = n) %>% 
  ungroup() %>% 
  add_count(word) %>%    # Count of number of times this word appears in all recipes
  rename(n_all_rec = n) %>%
  select(recipe_name, word, n_this_rec, n_all_rec)

# Get the total number of words per recipe
per_rec_totals <- per_rec_freq %>% 
  group_by(recipe_name) %>%
  summarise(total_this_recipe = sum(n_this_rec))

# Get the total number of times a word is used across all the recipes
all_rec_totals <- per_rec_freq %>% 
  ungroup() %>% 
  summarise(total_this_recipe = sum(n_this_rec))
  
# Join that on the sums we've found
per_rec_freq <- per_rec_freq %>% 
  mutate(
    total_overall = sum(n_this_rec)
  ) %>% 
  left_join(per_rec_totals) %>% 
  left_join(all_rec_totals)

per_rec_freq %>% 
  arrange(desc(n_all_rec)) %>% 
  slice(1:10) %>% 
  kable()

recipe_name word n_this_rec n_all_rec total_overall total_this_recipe


20-Minute Green Bean Chicken and Rice Dinner cream 1 4 141 27 Black Bean & Roasted Red Pepper Dip chopped 1 4 141 23 Black Bean & Roasted Red Pepper Dip cream 1 4 141 23 PHILADELPHIA Apple Crumble cream 1 4 141 18 PHILADELPHIA Apple Crumble chopped 1 4 141 18 Salmon On-Sea-of-Salt Bed cream 1 4 141 47 Salmon On-Sea-of-Salt Bed chopped 1 4 141 47 Saucy Apricot Chicken chopped 1 4 141 26 20-Minute Green Bean Chicken and Rice Dinner white 1 3 141 27 20-Minute Green Bean Chicken and Rice Dinner halves 1 3 141 27

Let's use our per_rec_freq to find the TFIDF (term frequency inverse document frequency), a measure of how frequently a word appears in this document as compared to how much it's seen in all the other documents. If a word has a high TFIDF, it's more common in this document (i.e., recipe) than it is in the general population of recipes, which is an indication that it's an important word for the recipe in question

recipe_tf_idf <- per_rec_freq %>% 
  bind_tf_idf(word, recipe_name, n_this_rec) %>% 
  arrange(desc(tf_idf))

recipe_tf_idf %>% 
  top_n(10) %>% 
  kable()
## Selecting by tf_idf

recipe_name word n_this_rec n_all_rec total_overall total_this_recipe tf idf tf_idf


PHILADELPHIA Apple Crumble tub 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble philadelphia 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble light 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble spread 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble sugar 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble crushed 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble nilla 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble vanilla 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble wafers 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble granny 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble smith 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble apples 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble peeled 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble ground 1 1 141 18 0.0555556 1.609438 0.0894132 PHILADELPHIA Apple Crumble cinnamon 1 1 141 18 0.0555556 1.609438 0.0894132

Next we'll get a measure of how frequently words co-occur with one another in these recipes. We'll find the pairwise correlation between words in each recipe using widyr::pairwise_cor(). The feature term we're grouping on is recipe_name.

pairwise_per_rec <- per_rec_freq %>% 
  pairwise_cor(word, recipe_name, sort = TRUE) 

pairwise_per_rec %>% sample_n(10) %>% kable()

item1 item2 correlation


taco cooking -0.2500000 divided milk -0.2500000 dill onion 0.6123724 drained seasonings 1.0000000 preserves wine 0.6123724 fresh apples -0.3750000 thawed philadelphia -0.2500000 finely cheese -0.4082483 cheese soup -0.4082483 roasted rinsed 1.0000000

Now that we have the correlations between each word, we can graph the relationships between words using igraph::graph_from_data_frame(), which takes a dataframe and turns it into a graph structure, and ggraph to display the network graph using ggplot2.

Since there are a lot of words, it'll make sense to choose just a few to make the graph readable. We'll also set a correlation_cutoff to keep only relationships above a certain threshold.

food_features <- c("cheese", "garlic", "onion", "sugar")
correlation_cutoff <- 0.5

Now all that remains is to graph the correlations between a few words and their highest correlated neighbors.

pairwise_per_rec %>%
  filter(item1 %in% food_features) %>% 
  filter(correlation > correlation_cutoff) %>%
  graph_from_data_frame() %>%
  ggraph(layout = "fr") +
  geom_edge_link(aes(edge_alpha = correlation), show.legend = FALSE) +
  geom_node_point(color = "lightblue", size = 5) +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_void()