Solver for minimum number of food items per day

LICENSE.md

Copyright 2026 Dyllon Gagnier

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Nutrition Per Dollar-v1.5.csv

solve_nutrition.py

# Data courtesy of https://efficiencyiseverything.com/food-nutrition-per-dollar/
import numpy as np
import pandas as pd
from scipy.optimize import milp, LinearConstraint, Bounds

food_data = pd.read_excel(r"C:\gdrive\code\python\nutrition-optim\Nutrition Per Dollar-v1.5.xlsx")
food_data = food_data[~food_data["Food"].isin(["Cilantro"])].reset_index(drop=True)

calories = 2000
bw = 62
vita_requirements = {
	"Vit_A_RAE": 900,          # µg RAE
	"Thiamin_(mg)": 1.2,       # mg
	"Riboflavin_(mg)": 1.3,    # mg
	"Niacin_(mg)": 16,         # mg
	"Panto_Acid_mg)": 5,       # mg  (AI)
	"Vit_B6_(mg)": 1.3,        # mg
	"Folate_Tot_(µg)": 400,    # µg DFE (RDA/AI target)
	"Vit_B12_(µg)": 2.4,       # µg
	"Vit_C_(mg)": 90,          # mg
	"Vit_E_(mg)": 15,          # mg
	"Vit_D_µg": 15,            # µg
	"Vit_K_(µg)": 120,         # µg  (AI)
	"Calcium_(mg)": 1000,      # mg
	"Iron_(mg)": 8,            # mg
	"Magnesium_(mg)": 420,     # mg
	"Potassium_(mg)": 3400,    # mg  (AI)
	"Fiber_TD_(g)": 38,        # g   (AI)
	"Phosphorus_(mg)": 700,    # mg
	"Zinc_(mg)": 11,           # mg
	"Copper_mg)": 0.9,         # mg (900 µg)
	"Selenium_(µg)": 55,       # µg
	"Choline_Tot_ (mg)": 550,   # mg (AI)
	"Energ_Kcal": calories,
	"Sodium_(mg)": 1500,
	"Protein_(g)": bw * 1.2
}

# Nutrient constraint types
lower_bound_vitas = {"Sodium_(mg)", "FA_Sat_(g)"}   # nutrient <= bound
equality_vitas    = {"Energ_Kcal"}                   # nutrient == bound
upper_bound_vitas = {"Protein_(g)"}                  # additional upper bound (separate row)
fat_cols = ["FA_Sat_(g)", "FA_Mono_(g)", "FA_Poly_(g)"]

foods = list(food_data["Food"])
n = len(foods)

# Total variables: n food quantities + n on/off binaries
n_vars = 2 * n

# Offset into variable vector
food_offset  = 0   # x[food_offset  + i] = grams of food i
onoff_offset = n   # x[onoff_offset + i] = binary on/off for food i

food_index = {food: i for i, food in enumerate(foods)}

INF = np.inf

A_rows = []
lb_rows = []
ub_rows = []

# Fat sum: FA_Sat + FA_Mono + FA_Poly in [0.2*calories/9, 0.35*calories/9] grams
# Fat provides 9 kcal/g, so grams = kcal_fraction * calories / 9
fat_vec = np.zeros(n_vars)
for col in fat_cols:
    if col not in food_data.columns:
        continue
    for i, (idx, _) in enumerate(food_data["Food"].items()):
        val = food_data.loc[idx, col]
        if not pd.isna(val):
            fat_vec[food_offset + i] += val / 100.0

A_rows.append(fat_vec.copy())
lb_rows.append(0.2 * calories / 9)   # grams: 20% of kcal ÷ 9 kcal/g
ub_rows.append(0.35 * calories / 9)  # grams: 35% of kcal ÷ 9 kcal/g

protein_row = None   # we'll save it for the protein upper-bound row
calcium_row = None   # we'll save it for the calcium upper-bound row

for vita, rda in vita_requirements.items():
    if vita not in food_data.columns:
        print(f"  WARNING: '{vita}' not found in food_data columns, skipping")
        continue

    row_vec = np.zeros(n_vars)
    for i, (idx, _) in enumerate(food_data["Food"].items()):
        val = food_data.loc[idx, vita]
        if pd.isna(val):
            continue
        row_vec[food_offset + i] = val / 100.0

    if vita in equality_vitas:
        A_rows.append(row_vec); lb_rows.append(rda);  ub_rows.append(rda)
    elif vita in lower_bound_vitas:
        A_rows.append(row_vec); lb_rows.append(-INF); ub_rows.append(rda)
    else:
        A_rows.append(row_vec); lb_rows.append(rda);  ub_rows.append(INF)

    if vita == "Protein_(g)":
        protein_row = row_vec.copy()
    if vita == "Calcium_(mg)":
        calcium_row = row_vec.copy()

# Protein upper bound: protein <= 2 * bw
if protein_row is not None:
    A_rows.append(protein_row); lb_rows.append(-INF); ub_rows.append(2 * bw)

# Calcium upper bound: calcium <= 2500 mg (tolerable upper intake level)
if calcium_row is not None:
    A_rows.append(calcium_row); lb_rows.append(-INF); ub_rows.append(2500)

print(f"Constraint matrix: {len(A_rows)} rows × {n_vars} columns")

M = 2_000.0  # big-M: no food exceeds 10 kg
M_overrides = {
	"Flour": 300,
	"Milk": 2000,
	"Eggs": 300,
	"Cereal": 200,
	"Romaine Lettuce": 500,
	"Kale": 400,
	"Carrots": 500,
	"Pinto Beans": 500,
	"Russet Potatoes": 800,
	"Black Beans": 500,
	"Peanuts": 120,
	"Spinach": 400,
	"Sweet Potatoes": 800,
	"Turnip Greens": 400,
	"Split Peas": 500,
	"Brown Rice": 600,
	"Butternut Squash": 600,
	"Cheese": 200,
	"Frozen Carrots": 500,
	"Green Cabbage": 500,
	"Light Red Kidney Beans": 500,
	"Canned Tomatoes": 800,
	"Lentils": 500,
	"Wheat Noodles": 600,
	"Curly Parsley": 100,
	"Banana": 400,
	"Red Potatoes": 800,
	"Jalapeno Peppers": 100,
	"Collard Greens": 400,
	"Frozen Peas": 400,
	"Red Cabbage": 500,
	"White Rice": 600,
	"Yogurt": 1000,
	"White Noodles": 600,
	"Organic Bananas": 400,
	"Pork": 400,
	"Ground Turkey": 400,
	"Monster": 500,
	"Flat Italy Parsley": 100,
	"Iceberg lettuce": 500,
	"Watermellon": 1000,
	"Roma Tomatoes": 600,
	"Green Pepper": 400,
	"Corn Oil": 60,
	"Frozen Corn": 400,
	"Oranges": 600,
	"Plantains": 400,
	"Pineapple": 600,
	"Garlic": 30,
	"Chicken Breast": 400,
	"Veg Juice": 750,
	"Mini Watermellons": 1000,
	"Cantaloupe": 800,
	"Beats": 400,
	"Cranberry Juice": 500,
	"Peanut Oil": 60,
	"Yellow Onions": 300,
	"Sweet Onions": 300,
	"Blackberries": 300,
	"Atuulfo Mangos": 500,
	"Red Onions": 300,
	"Yellow Squash": 600,
	"Maradol Papayas": 600,
	"Avacados": 300,
	"White Onions": 300,
	"Green Beans": 450,
	"Spaghetti Squash": 600,
	"Almonds": 120,
	"Grapefruit": 500,
	"Mangos": 500,
	"Golden Honey Dew": 800,
	"Red Pepper": 400,
	"Radishes": 300,
	"Leeks": 300,
	"Strawberries": 400,
	"Grapes": 500,
	"Cilantro": 100,
	"Asparagus": 500,
	"Red Delicious": 400,
	"Olive Oil": 60,
	"Bose Pears": 400,
	"Red Pears": 400,
	"Lemons": 150,
	"Anjou Pears": 400,
	"Golden Delicious": 400,
	"Celery": 400,
	"Gala": 400,
	"Fuji Apder": 400,
	"Wine": 300,
	"Grany Smith": 400,
	"Artichokes": 400,
	"Limes ": 150,
	"Rockstar": 500,
	"Raspberries": 300,
	"Pink Lady": 400,
	"Peaches": 400,
	"Honey Crisp": 400,
	"Barlet Pears": 400,
	"Budlight": 710,
	"Ginger Root": 6,
	"Blueberries": 300,
	"Necterines": 400,
	"Redbull": 500,
	"Plums": 400,
	"Coffee": 720,
	"Black Tea": 1500,
	"Green Tea": 1500,
}
for i in range(n):
    row_vec = np.zeros(n_vars)
    row_vec[food_offset  + i] =  1.0   # x[i]
    row_vec[onoff_offset + i] = -M_overrides.get(foods[i], M)     # -M * y[i]
    A_rows.append(row_vec)
    lb_rows.append(-INF)
    ub_rows.append(0.0)   # x[i] - M*y[i] <= 0

# Stack constraint matrix
A  = np.array(A_rows,  dtype=float)
lb = np.array(lb_rows, dtype=float)
ub = np.array(ub_rows, dtype=float)

constraints = LinearConstraint(A, lb, ub)

# Integrality: 0 = continuous (food quantities), 1 = integer (on/off)
integrality = np.array([0] * n + [1] * n, dtype=float)

# Variable bounds: food_qty >= 0 (no upper), onoff in [0, 1]
var_lb = np.array([0.0] * n + [0.0] * n)
var_ub = np.array([INF]  * n + [1.0]  * n)
bounds = Bounds(lb=var_lb, ub=var_ub)

c = np.array([0.0] * n + [1.0] * n)

result = milp(
    c,
    constraints=constraints,
    integrality=integrality,
    bounds=bounds,
)

print(f"Status  : {result.status}")
print(f"Message : {result.message}")
if result.success:
    print(f"Optimal foods used: {int(round(result.fun))}")

if not result.success:
    print("Solver did not find a solution:", result.message)
else:
    x_sol = result.x
    food_qtys = x_sol[food_offset : food_offset + n]   # grams per food
    onoff_vals = x_sol[onoff_offset : onoff_offset + n]  # binary on/off

    print("=== Foods used ===")
    # Use the binary on/off variable (> 0.5) as the selection indicator,
    # not the continuous quantity (which may be small but still nonzero)
    used = [(foods[i], food_qtys[i]) for i in range(n) if onoff_vals[i] > 0.5]
    for food, qty in used:
        print(f"  {food:40s}  {qty:8.1f} g")

    print()
    print("=== Nutrient totals vs. targets ===")
    header = f"{'Nutrient':<25} {'Achieved':>10} {'Target':>10} {'Type':<10}"
    print(header)
    print("-" * len(header))
    for vita, rda in vita_requirements.items():
        if vita not in food_data.columns:
            continue
        total = sum(
            food_qtys[i] * (food_data.loc[idx, vita] / 100.0)
            for i, (idx, _) in enumerate(food_data["Food"].items())
            if not pd.isna(food_data.loc[idx, vita])
        )
        if vita in equality_vitas:
            kind = "=="
        elif vita in lower_bound_vitas:
            kind = "<="
        else:
            kind = ">="
        print(f"  {vita:<23} {total:>10.2f} {rda:>10.2f}  {kind}")