Automatic planning of cheap, varied, nutritional meals

mealplanner.py

# -*- coding: utf-8 -*-
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
MEAL PLANNER - Automatic planning of cheap, varied, nutritional meals
---------------------------------------------------------------------

Summary
-------
This Python module uses optimization (Linear Programming) and stochastic methods to
generate a meal with lots of variety while still obeying user-defined nutritional
constraints. Solving the Stigler Diet problem directly yields sparse, boring solutions.
This program samples some foods, solves a small Stigler Diet problem, and then evalutes
the many solutions to come up with a meal plan which has good variety.

Created by GitHub users @tommyod and @JensWahl on a slow weekend, July 2019.

Installation
------------
This program was tested with Python 3.7.3 and ortools 7.2.6977. To install and run:
    1. Install Python 3.7
    2. Install ortools (`pip install ortools`)
    3. Run `python <this file>`

If you have markdown2 (`pip install markdown2`), you can render the output of this
program as a HTML file and print it for your fridge. Tested with markdown2 2.3.8.
    1. $ python <this file> > meal.md
    2. $ markdown2 meal.md > meal.html

"""


import operator
import dataclasses
import collections
import warnings
import statistics
import pprint
import typing
import random
import copy

from ortools.linear_solver import pywraplp

# =============================================================================
# CLASSES - Used to store data and make some computations more easily
# =============================================================================


@dataclasses.dataclass(frozen=True)
class Food:
    """
    A food consists of a name and nutritional data given in units of 100 grams.
    
    Examples
    --------
    >>> eggs = Food(name='eggs', protein=13.0, fat=10.6, carbs=0.3, kcal=149, 
    ...             price_per_product=32.9, grams_per_product=690)
    >>> eggs.price
    4.768115942028985
    >>> eggs.name
    'eggs'
    """

    name: str

    # Values are per 100 grams of food
    protein: float
    fat: float
    carbs: float
    kcal: float

    # Used to infer the price per 100 grams
    price_per_product: float
    grams_per_product: float

    @property
    def price(self):
        return self.price_per_product / self.grams_per_product * 100

    def __post_init__(self):
        """Verify the relationship between macros and kcal."""
        computed_kcal = 4 * self.protein + 4 * self.carbs + 9 * self.fat
        relative_error = abs((self.kcal - computed_kcal) / computed_kcal)
        if relative_error > 0.1:
            warnings.warn(f"Got a {relative_error:.2f} error on kcal: '{self.name}'.")


@dataclasses.dataclass(frozen=True)
class Meal:
    """
    A meal consists of several foods given at some "base" unit of grams. In the example
    below 'eggs' is the general food, an a meal consists of discrete units of eggs,
    weighing approximately 65 grams each.
    
    Examples
    --------
    >>> eggs = Food(name='eggs', protein=13.0, fat=10.6, carbs=0.3, kcal=149, 
    ...             price_per_product=32.9, grams_per_product=690)
    >>> egg = Meal(name='egg', foods={eggs:65}, discrete=True)
    >>> egg.price # price of the meal
    3.0992753623188407
    """

    # Foods are added as: foods={all_foods["lettmelk"]:100, all_foods["musli"]:100}
    # This means that a baseline
    name: str
    foods: typing.Dict[Food, float] = dataclasses.field(default_factory=dict)
    discrete: bool = True
    type: str = None

    def __getattr__(self, key):
        """Allow accessing attributes of foods, summing over them."""
        return sum(
            getattr(food, key) * quantity / 100
            for (food, quantity) in self.foods.items()
        )

    def __hash__(self):
        return hash(self.name) + hash(frozenset(self.foods.keys()))

    def __iter__(self):
        return iter(self.foods.keys())

    def __copy__(self):
        return type(self)(
            name=self.name, foods=self.foods.copy(), discrete=self.discrete
        )

    def __str__(self):
        name = type(self).__name__
        foods_names = (
            "{"
            + ", ".join(
                f"{quantity}g {food.name}" for food, quantity in self.foods.items()
            )
            + "}"
        )
        return name + f"(name='{self.name}', foods={foods_names})"


class Day(collections.UserList):
    def __getattr__(self, key):
        """Allow accessing attributes of meals"""
        return sum(getattr(meal, key) for meal in self.data)


# =============================================================================
# FUNCTIONS - Used to generate the meal plan
# =============================================================================
def sample_single_day(meals, *, num_meals=3, type_constraints=None):
    """Sample a day consisting of several meals, meeting all meal type constraints."""

    if type_constraints is None:
        type_constraints = dict()
    else:
        type_constraints = copy.deepcopy(type_constraints)

    type_constraints[None] = num_meals - sum(type_constraints.values())

    chosen_meals = Day()

    for meal_type, num in type_constraints.items():

        # Filters the meal of a given meal type
        valid_meals = [meal for meal in meals if meal.type == meal_type]

        # Choose the desired number of meals of the type and add to chosen list
        to_add = random.choices(valid_meals, k=num)
        chosen_meals += to_add

    assert len(chosen_meals) == num_meals
    return chosen_meals


def optimize_day(day, *, dietary_constraints, price_epsilon=0.001, low_range=(1, 4)):
    """Optimize the quantitiy of each meal in a day, given constraints."""

    # Create a solver and an objective function
    solver = pywraplp.Solver("food", pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
    objective_function = 0
    INF = solver.infinity()

    # Create variables representing quantities of food
    solution = []
    for i, food in enumerate(day):

        # A lower limit for the quanitity, forces variation in the solutions
        low_limit = random.choice(range(*low_range))

        # Create a variable - use str(i) to avoid duplicate names if same foods are
        # present in the solution
        if food.discrete:
            x = solver.IntVar(low_limit, INF, food.name + str(i))
        else:
            x = solver.NumVar(low_limit, INF, food.name + str(i))

        # Gently perturb the solution towards cheaper combinations
        objective_function += price_epsilon * food.price * x
        solution.append(x)

    # Set up the contraints. Here goal programming is used to create "soft" constraints
    # instead of hard contraints, since a feasible solution might not exist in general
    # with hard constraints: https://en.wikipedia.org/wiki/Goal_programming
    for macro, (low, high) in dietary_constraints.items():

        food_macros = [getattr(meal, macro) for meal in day]

        # Create the sum: sum_i food_i * macro_i
        total_macro = sum(c * x for x, c in zip(solution, food_macros))

        # Slack variables related to the lower limit. Only "undershooting" is penalized.
        low_positive = solver.NumVar(0, INF, "over_lower_limit_" + macro)
        low_negative = solver.NumVar(0, INF, "under_lower_limit_" + macro)
        solver.Add(total_macro + low_positive - low_negative == low)

        # Slack variables related to the upper limit. Only "overshooting" is penalized.
        high_positive = solver.NumVar(0, INF, "over_upper_limit_" + macro)
        high_negative = solver.NumVar(0, INF, "under_upper_limit_" + macro)
        solver.Add(total_macro + high_positive - high_negative == high)

        # Scale the influence of calorier compared to the macros
        scale = 0.2 if macro == "kcal" else 1

        objective_function += scale * (low_positive + high_negative)

    # Minimize the deviation from the goal (with a preference for low price too)
    solver.Minimize(objective_function)
    result_status = solver.Solve()

    assert result_status == pywraplp.Solver.OPTIMAL
    assert solver.VerifySolution(1e-7, True)

    answer = [(food, round(x.solution_value(), 1)) for food, x in zip(day, solution)]
    return answer, round(solver.Objective().Value(), 2)


def sort_by_penality(*, days, obj_func_values):
    """Choose the best days based on heuristics."""

    def get_price(result):
        """Preference towards cheaper days."""
        return sum(food.price * num for (food, num) in result)

    def calorie_spread(result):
        """Preference for lower standard deviation of calories."""
        calories = [food.kcal * num for food, num in result]
        return statistics.pstdev(calories)

    def get_duplicates(result):
        """Preference for many distinct meals in a day."""
        return len(result) - len(set([f for f, num in result]))

    def standarize(values):
        """Standardize the penality scores."""
        mu = statistics.mean(values)
        sigma = statistics.pstdev(values)
        return [0 if sigma == 0 else (x - mu) / sigma for x in values]

    prices = standarize([get_price(result) for result in days])
    calorie_stds = standarize([calorie_spread(result) for result in days])
    duplicates = standarize([get_duplicates(result) for result in days])
    obj_func_values = standarize(obj_func_values)

    # A user can change these weights to penalize different properties
    generator = zip(prices, calorie_stds, duplicates, obj_func_values)
    total_penalty = [p * 1 + c * 1 + d * 1 + o * 5 for (p, c, d, o) in generator]

    # Join the solutions with the scores in preparation of argsort
    joined = sorted(zip(days, total_penalty), key=operator.itemgetter(1))
    return list(map(operator.itemgetter(0), joined))


def print_results(results, *, verbose=True):
    """Print the results in markdown."""
    print(f"# Meal plan")
    for day_num, result in enumerate(results, 1):

        # Heuristics to get more carbohydrates earlier in the day
        result = sorted(result, key=lambda r: r[0].carbs, reverse=True)

        price = int(sum(meal.price * qnty for (meal, qnty) in result))
        print(f"\n## Day {day_num} (price: {price} NOK)")

        print(f"\n### Meals\n")
        for meal, qnty in result:
            print(f"- {qnty} x {str(meal)}")

        if not verbose:
            continue

        print(f"\n### Statistics\n")
        for macro in ["kcal", "protein", "fat", "carbs"]:
            macro_distribution = [
                int(getattr(meal, macro) * qnty) for (meal, qnty) in result
            ]
            print(f"- Total {macro}: {sum(macro_distribution)} {macro_distribution}")


# =============================================================================
# DATA - Stored in dataclasses for convenience
# =============================================================================
def get_data():
    """Stored in a function so as not to pollute global namespace."""

    all_foods = [Food(name='bacon', protein=14.0, fat=32.0, carbs=1.0, kcal=350, price_per_product=50.9, grams_per_product=400),
                 Food(name='burger', protein=15.0, fat=18.0, carbs=2.0, kcal=230, price_per_product=87.3, grams_per_product=800),
                 Food(name='coop bratwurst', protein=12.0, fat=20.0, carbs=5.2, kcal=253, price_per_product=25.9, grams_per_product=240),
                 Food(name='cottage cheese', protein=13.0, fat=2.0, carbs=2.1, kcal=79, price_per_product=24.4, grams_per_product=400),
                 Food(name='egg', protein=13.0, fat=10.6, carbs=0.3, kcal=149, price_per_product=32.9, grams_per_product=690),
                 Food(name='frossen kyllingfilet', protein=19.0, fat=1.8, carbs=0.3, kcal=94, price_per_product=260.0, grams_per_product=2500),
                 Food(name='grovt brød', protein=11.0, fat=4.8, carbs=36.0, kcal=245, price_per_product=39.5, grams_per_product=750),
                 Food(name='gulost', protein=27.0, fat=27.0, carbs=0.0, kcal=351, price_per_product=110.0, grams_per_product=1000),
                 Food(name='jasmin ris', protein=2.7, fat=0.1, carbs=31.1, kcal=136, price_per_product=45.8, grams_per_product=1000),
                 Food(name='kjøttdeig', protein=19.0, fat=9.0, carbs=0.0, kcal=157, price_per_product=32.5, grams_per_product=400),
                 Food(name='lettmelk', protein=3.5, fat=0.5, carbs=4.5, kcal=37, price_per_product=16.4, grams_per_product=1000),
                 Food(name='melkesjokolade', protein=8.1, fat=33.0, carbs=55.0, kcal=550, price_per_product=38.6, grams_per_product=200),
                 Food(name='musli', protein=9.0, fat=4.8, carbs=63.0, kcal=351, price_per_product=23.1, grams_per_product=750),
                 Food(name='PF whey', protein=71.8, fat=8.1, carbs=7.9, kcal=377, price_per_product=599.0, grams_per_product=3000),
                 Food(name='svinekotelett dypfryst', protein=20.0, fat=18.0, carbs=0.0, kcal=243, price_per_product=98.6, grams_per_product=2000),
                 Food(name='sweet and sour sauce', protein=0.4, fat=0.2, carbs=16.4, kcal=71, price_per_product=35.0, grams_per_product=675),
                 Food(name='tomatbønner', protein=3.8, fat=0.6, carbs=14.0, kcal=83, price_per_product=11.8, grams_per_product=420),
                 Food(name='yoghurt', protein=3.7, fat=3.1, carbs=10.5, kcal=84, price_per_product=17.0, grams_per_product=600)]
    
    

    foods = {food.name: food for food in all_foods}
    #pprint.pprint(sorted(foods.values(), key=lambda f:f.name.lower()))

    all_meals = [
        Meal(
            name="pølse i skive",
            foods={foods["grovt brød"]: 40, foods["coop bratwurst"]: 80},
        ),
        Meal(
            name="pølse, skive og 2 egg",
            foods={
                foods["grovt brød"]: 40,
                foods["egg"]: 120,
                foods["coop bratwurst"]: 80,
            },
        ),
        Meal(
            name="burger, skive og egg",
            foods={foods["grovt brød"]: 40, foods["egg"]: 60, foods["burger"]: 80},
        ),
        Meal(
            name="1 yoghurt med ct.cheese",
            foods={foods["yoghurt"]: 150, foods["cottage cheese"]: 100},
        ),
        Meal(
            name="2 yoghurt med musli",
            foods={foods["yoghurt"]: 300, foods["musli"]: 100},
        ),
        Meal(
            name="2 scoops whey", foods={foods["PF whey"]: 45, foods["lettmelk"]: 300}
        ),
        Meal(
            name="melk og musli",
            foods={foods["lettmelk"]: 150, foods["musli"]: 100},
            discrete=False,
        ),
    ]
    
    
    meals = {meal.name: meal for meal in all_meals}
    return foods, meals


def main():
    """Main function: the user defines parameters and gets a high level overview."""
    random.seed(123)
    foods, meals = get_data()

    # =============================================================================
    #     USER DEFINED HYPERPARAMETERS ARE PRIMARILY SET HERE
    # =============================================================================

    # User defined constraints for types of meals, e.g. {"breakfast":1, "dinner":1}
    # Some Meal-instances must have Meal(..., type="dinner"). This can be used to
    # enforce one dinner every day, or one liquid meal every day, etc.
    type_constraints = dict()

    num_days = 7  # Total number of days for the meal program
    num_meals = 4  # Number of meals per day
    dietary_constraints = {"kcal": (1800, 2000), "protein": (120, 200)}

    # Algorithm hyperparameters
    num_samples = 1000  # Number of completely random samples to draw
    top_samples_to_keep = 100  # Top samples to retain

    # =============================================================================
    #     SOLVE THE PROBLEM USING OPTIMIZATION AND HEURISTICS
    # =============================================================================

    # Sample different possible days
    meal_objects = list(meals.values())
    optimized_results, obj_func_values = list(), list()
    for _ in range(num_samples):

        # Sample a possible day (i.e. a collection of meals)
        day = sample_single_day(
            meal_objects, num_meals=num_meals, type_constraints=type_constraints
        )

        # Optimize how much of each meal should go into the day to meet dietary contraints
        result, obj_func_value = optimize_day(
            day,
            dietary_constraints=dietary_constraints,
            price_epsilon=0.001,
            low_range=(1, 4),
        )

        # Store values
        optimized_results.append(result)
        obj_func_values.append(obj_func_value)

    # Sort the days by their "penalty", a weighted sum of terms
    sorted_results = sort_by_penality(
        days=optimized_results, obj_func_values=obj_func_values
    )
    sorted_results = sorted_results[:top_samples_to_keep]
    weights = list(reversed(range(1, len(sorted_results) + 1)))
    sum_weights = sum(weights)
    weights = [w / sum_weights for w in weights]

    # Choose among the best days, weighted by how good they are
    results_solved = random.choices(sorted_results, weights=weights, k=num_days)
    print_results(results_solved, verbose=True)


if __name__ == "__main__":
    main()