<jupyter_start><jupyter_text>IntroductionThe Ugrainium recipe uses groupings of ingredients to try to cover nutritional requirements with enough overlap that a single ingredient with quality issues does not cause a failure for the whole recipe. A opimizer is used to find the right amount of each ingredient to fulfill the nutritional and practicle requirements for a recipe. To Do* Nutrients without an upper limit should have the upper limit constraint removed* Add constraints for the NIH essential protein combinations as a limit* Add a radar graph for vitamins showing the boundry between RDI and UL* Add a radar graph for vitamins without an upper limit but showing the RDI* Add a radar graph for essential proteins showing the range between RDI and UL* Add a radar graph for essential proteins without an upper limit, but showing the RDI as the lower limit* Add a radar graph pair for non-essential proteins with the above UL and no UL pairing* Add equality constraints for at least energy, and macro nutrients if possible<jupyter_code># Import all of the helper libraries 

from scipy.optimize import minimize
from scipy.optimize import Bounds
from scipy.optimize import least_squares, lsq_linear
import pandas as pd
import numpy as np
import os
import json
import requests
from math import e
import matplotlib.pyplot as plt
import seaborn as sns
from math import log10
from ipysheet import from_dataframe, to_dataframe
#%matplotlib inline
#!pip install seaborn
#!pip install ipysheet
#!pip install ipywidgets
# Setup the notebook context

data_dir = '../data'
pd.set_option('max_columns', 70)<jupyter_output><empty_output><jupyter_text>Our DataThe [tables](https://docs.google.com/spreadsheets/d/1UOx5ZkNMjcwybx9TgFS_EaQnjbokqDx6K9NW0mRDGlg/editgid=6244197120) containing our ingredients nutrition profile are held in Google Sheets.The sheet names are "Ingredients" and "Nutrition Profile"<jupyter_code># Download our nutrition profile data from Google Sheets
!(cd ..; make sync_data_from_google)

nutrition_profile_df = pd.read_csv(f'{data_dir}/raw/nutrition_profile.csv', index_col=0, verbose=True).transpose()
ingredients_df = pd.read_csv(f'{data_dir}/raw/ingredients.csv', index_col=0, verbose=True).transpose()

# convert all values to float
for col in ingredients_df.columns:
    ingredients_df[col] = ingredients_df[col].astype(float)

#sheet = from_dataframe(ingredients_df)
#sheet
ingredients_df<jupyter_output><empty_output><jupyter_text>To install ipysheet, the following commands need to be run:- conda install ipywidgets- jupyter labextension install @jupyter-widgets/jupyterlab-managersee: https://stackoverflow.com/questions/36351109/ipython-notebook-ipywidgets-does-not-show<jupyter_code>from ipysheet import sheet, cell, row, column, cell_range
from ipywidgets import interact, interactive, fixed, interact_manual
sheet1 = sheet(rows=3, columns=4)
display(sheet1)<jupyter_output><empty_output><jupyter_text>ingredient_map = {    'wheat bran': '785757',    'oat flour': '169741',    'corn flour': '169694',    'whey protein isolate': '173177',    'soy protein isolate': '174276',    'lecithin': '171426',    'canola oil': '172336',    'barley flour': '169739',    'flaxseed': '169414',    'guar gum': '169045',    'xanthan gum': '563344',}def get_fdc(fdc_id):    key = '<KEY>'    payload={'api_key': key}    url = f'https://api.nal.usda.gov/fdc/v1/food/{fdc_id}'    r = requests.get(url, params=payload)    def nutrient_decode(el):        if 'nutrient' in el:            if 'amount' in el:                return((f'{el["nutrient"]["name"]} ({el["nutrient"]["unitName"]})', el['amount']))            else:                return(None)        return(el)    print(f'result: {r}')    res = json.loads(r.text, object_hook=nutrient_decode)    return pd.Series(dict(filter(lambda el: el is not None, res['foodNutrients'])), name=res['description'])    ing_list = []for (ing_name, fdc_id) in ingredient_map.items():    ing_list.append(get_fdc(fdc_id))ing_df = pd.DataFrame(ing_list)<jupyter_code>ing_df = ingredients_df
print(str(sorted(ing_df.columns)))
d = {'a':'b'}
d.items()
l = [d]
list(d.keys())[0]
tuple(d)
tuple(d.items())<jupyter_output><empty_output><jupyter_text>Problem SetupLet's cast our data into the from $\vec{y} = A \vec{x} + \vec{b}$ where $A$ is our ingredients data, $\vec{x}$ is the quantity of each ingredient for our recipe, and $\vec{b}$ is the nutrition profile.The problem to be solved is to find the quantity of each ingredient which will optimally satisfy the nutrition profile, or in our model, to minimize: $|A \vec{x} - \vec{b}|$.There are some nutrients we only want to track, but not optimize. For example, we want to know how much cholesterol is contained in our recipe, but we don't want to constrain our result to obtain a specific amount of cholesterol as a goal. The full list of ingredients are named: A_full, and b_full. The values to optimized are named: A and b<jupyter_code>b_full = nutrition_profile_df
A_full = ingredients_df.transpose()
A = ingredients_df.transpose()[nutrition_profile_df.loc['Report Only'] == False]

b_full = nutrition_profile_df.loc['Target']
b = nutrition_profile_df.loc['Target'][nutrition_profile_df.loc['Report Only'] == False]
ul = nutrition_profile_df.loc['UL'][nutrition_profile_df.loc['Report Only'] == False]
rdi = nutrition_profile_df.loc['RDI'][nutrition_profile_df.loc['Report Only'] == False]
weight = nutrition_profile_df.loc['Weight'][nutrition_profile_df.loc['Report Only'] == False]
ul_full = nutrition_profile_df.loc['UL']
rdi_full = nutrition_profile_df.loc['RDI']
#weight<jupyter_output><empty_output><jupyter_text>Least Squares OptimizationStart with a simple solver using a least squares optimizer<jupyter_code># Using the scipy least_squares optimizer
# Define a residual function: if y = A x + b -> A x - b is the residual function
k_thres = 0.7
k_rdi = b - ((rdi - b) * k_thres)
k_ul =  b + ((ul - b) * k_thres)

def limit_scale(df):
    r = ((df['y'] - df['b']) /  df['b']) * df['weight']
    k_thres = 0.7
    if df['y'] < df['k_rdi']:
        r *= 10
        if df['y'] < df['rdi']:
            r *= 10
    elif df['y'] >= df['b']:
        if df['k_ul'] is np.NaN:
            return(0.0)
        elif df['y'] > df['k_ul']:
            r *= 10
        if df['y'] > df['ul']:
            r *= 10
    return(r)

scale_limits = pd.DataFrame({'rdi': rdi, 'k_rdi': k_rdi, 'b': b, 'k_ul': k_ul, 'ul': ul, 'weight': weight})

def residuals(x, *args, **kwargs):
    y0 = A.dot(x.transpose())
    res0 = pd.concat([scale_limits, pd.Series(y0,scale_limits.index, name='y')], axis=1)
    res1 = res0.apply(limit_scale, axis=1)
    #print(f'x: {x}, residual: {res1}')
    return(res1)
# Constrain ingredients before the optimization process. Many of the ingredients are required for non-nutritional purposes 
# or are being limited to enhance flavor
#
# The bounds units are in fractions of 100g / day, i.e.: 0.5 represents 50g / day, of the ingredient

bounds_df = pd.DataFrame(index=ingredients_df.index, data={'lower': 0, 'upper': np.inf})             
bounds_df.loc['Guar gum'] = [1.0 * .01, 1.0 * .01 + .0001]
bounds_df.loc['Xanthan Gum'] = [1.0 * .01, 1.0 * .01 + .0001]
bounds_df.loc['Alpha-galactosidase enzyme (Beano)'] = [1.0, 1.0 + .0001]
bounds_df.loc['Multivitamin'] = [1.0, 1.0 + .0001]
bounds_df.loc['Corn flour, nixtamalized'] = [0, 1.0]
bounds_df.loc['Whey protein'] = [0,0.15]
bounds_df.loc['Ascorbic acid'] = [0.01, 0.01 + .0001]
bounds_df.loc['Peanut butter'] = [0.25, 0.50]
lower = list(bounds_df.lower.values)
upper = list(bounds_df.upper.values)

# x0 is our initial guess at a solution
x0 = np.array(lower)

len(lower), len(upper), x0
solution = least_squares(residuals, x0, jac='2-point', bounds=(lower, upper), method='trf', ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, loss='linear', f_scale=1.0, diff_step=None, tr_solver=None, tr_options={}, jac_sparsity=None, max_nfev=None, verbose=1, args=(), kwargs={})
solution.success, solution.x
A.shape, b.shape
b
#solution = lsq_linear(A, b, bounds=(lower, upper), method='trf', lsq_solver='lsmr', verbose=2)
solution = lsq_linear(A, b, method='trf', lsq_solver='lsmr', verbose=2)
solution.success, solution.x
A_full.dot(solution.x)
# Scale the ingredient nutrient amounts for the given quantity of each ingredient given by the optimizer
solution_df = A_full.transpose().mul(solution.x, axis=0) # Scale each nutrient vector per ingredient by the amount of the ingredient
solution_df.insert(0, 'Quantity (g)', solution.x * 100) # Scale to 100 g since that is basis for the nutrient quantities

# Add a row showing the sum of the scaled amount of each nutrient
total = solution_df.sum()
total.name = 'Total'
solution_df = solution_df.append(total)
# Plot the macro nutrient profile
# The ratio of Calories for protein:carbohydrates:fat is 4:4:9 kcal/g
pc = solution_df['Protein (g)']['Total'] * 4.0
cc = solution_df['Carbohydrates (g)']['Total'] * 4.0
fc = solution_df['Total Fat (g)']['Total'] * 9.0
tc = pc + cc + fc
p_pct = int(round(pc / tc * 100))
c_pct = int(round(cc / tc * 100))
f_pct = int(round(fc / tc * 100))
(p_pct, c_pct, f_pct)
# create data
names=f'Protein {p_pct}%', f'Carbohydrates {c_pct}%', f'Fat {f_pct}%', 
size=[p_pct, c_pct, f_pct]
 
fig = plt.figure(figsize=(10, 5))

fig.add_subplot(1,2,1)
# Create a circle for the center of the plot
my_circle=plt.Circle( (0,0), 0.5, color='white')

# Give color names
cmap = plt.get_cmap('Spectral')
sm = plt.cm.ScalarMappable(cmap=cmap)
colors = ['yellow','orange','red']
plt.pie(size, labels=names, colors=colors)

#p=plt.gcf()
#p.gca().add_artist(my_circle)
fig.gca().add_artist(my_circle)
#plt.show()

fig.add_subplot(1,2,2)
barWidth = 1
fs = [solution_df['Soluble Fiber (g)']['Total']]
fi = [solution_df['Insoluble Fiber (g)']['Total']]
plt.bar([0], fs, color='red', edgecolor='white', width=barWidth)
plt.bar([0], fi, bottom=fs, color='yellow', edgecolor='white', width=barWidth, label=['Insoluble Fiber (g)'])
 
# Also show the Omega-3, Omega-6 ratio
# Saturated:Monounsaturated:Polyunsaturated ratios
# Prepare data as a whole for plotting by normalizing and scaling
amounts = solution_df
total = A_full.dot(solution.x) #solution_df.loc['Total']

# Normalize as a ratio beyond RDI
norm = (total) / rdi_full
norm_ul = (ul_full) / rdi_full

nuts = pd.concat([pd.Series(norm.values, name='value'), pd.Series(norm.index, name='name')], axis=1)
# Setup categories of nutrients and a common plotting function
vitamins = ['Vitamin A (IU)','Vitamin B6 (mg)','Vitamin B12 (ug)','Vitamin C (mg)','Vitamin D (IU)',
            'Vitamin E (IU)','Vitamin K (ug)','Thiamin (mg)','Riboflavin (mg)','Niacin (mg)','Folate (ug)','Pantothenic Acid (mg)','Biotin (ug)','Choline (mg)']
minerals = ['Calcium (g)','Chloride (g)','Chromium (ug)','Copper (mg)','Iodine (ug)','Iron (mg)',
            'Magnesium (mg)','Manganese (mg)','Molybdenum (ug)','Phosphorus (g)','Potassium (g)','Selenium (ug)','Sodium (g)','Sulfur (g)','Zinc (mg)']
essential_aminoacids = ['Cystine (mg)','Histidine (mg)','Isoleucine (mg)','Leucine (mg)','Lysine (mg)',
                        'Methionine (mg)','Phenylalanine (mg)','Threonine (mg)','Tryptophan (mg)','Valine (mg)']
other_aminoacids = ['Tyrosine (mg)','Arginine (mg)','Alanine (mg)','Aspartic acid (mg)','Glutamic acid (mg)','Glycine (mg)','Proline (mg)','Serine (mg)','Hydroxyproline (mg)']

def plot_group(nut_names, title):
    nut_names_short = [s.split(' (')[0] for s in nut_names] # Snip off the units from the nutrient names
    
    # Create a bar to indicate an upper limit 
    ul_bar = (norm_ul * 1.04)[nut_names]
    ul_bar[ul_full[nut_names].isnull() == True] = 0
    
    # Create a bar to mask the UL bar so just the end is exposed
    ul_mask = norm_ul[nut_names]
    ul_mask[ul_full[nut_names].isnull() == True] = 0
    
    
    n = []  # normalized values for each bar
    for x, mx in zip(norm[nut_names], ul_mask.values):
        if mx == 0: # no upper limit
            if x < 1.0:
                n.append(1.0 - (x / 2.0))
            else:
                n.append(0.50)
        else:
            n.append(1.0 - (log10(x) / log10(mx)))
    clrs = sm.to_rgba(n, norm=False)
    
    g = sns.barplot(x=ul_bar.values, y=nut_names_short, color='red')
    g.set_xscale('log')
    sns.barplot(x=ul_mask.values, y=nut_names_short, color='white')
    bax = sns.barplot(x=norm[nut_names], y=nut_names_short, label="Total", palette=clrs)
    
    # Add a legend and informative axis label
    g.set( ylabel="",xlabel="Nutrient Mass / RDI (Red Band is UL)", title=title)

    #sns.despine(left=True, bottom=True)
# Construct a group of bar charts for each nutrient group

# Setup the colormap for each bar
cmap = plt.get_cmap('Spectral')
sm = plt.cm.ScalarMappable(cmap=cmap)

#fig = plt.figure(figsize=plt.figaspect(3.))
fig = plt.figure(figsize=(20, 20))
fig.add_subplot(4, 1, 1)
plot_group(vitamins,'Vitamin amounts relative to RDI')
fig.add_subplot(4, 1, 2)
plot_group(minerals,'Mineral amounts relative to RDI')
fig.add_subplot(4, 1, 3)
plot_group(essential_aminoacids,'Essential amino acid amounts relative to RDI')
fig.add_subplot(4, 1, 4)
plot_group(other_aminoacids,'Other amino acid amounts relative to RDI')

#fig.show()
fig.tight_layout()
solu_amount = solution_df['Quantity (g)']
solu_amount.index.name = 'Ingredient'
solu_amount.reset_index()
trad_amount = pd.Series({
    'Lecithin granules': 10.0,
    'Canola oil': 19.0,
    'Xanthan Gum': 0.9,
    'Guar gum': 0.9,
    'Salt, iodized': 3.0,
    'Potassium chloride': 5.0,
    'Ascorbic acid': 1.0,
    'Vitamin K': 20.0,
    'Multivitamin': 100.0,
    'Calcium and vitamin D': 80.0,
    'Choline Bitartrate': 2.7,
    'Asprin': 0.5,
    'Alpha-galactosidase enzyme (Beano)': 2.0,
    'Flaxseed, fresh ground': 50.0,
    'Wheat bran, crude': 36.0,
    'Whey protein': 13.0,
    'Pea Protein Powder': 26.0,
    'Soy protein isolate': 29.0,
    'Oat flour': 59.0,
    'Barley flour': 104.0,
    'Corn flour, nixtamalized': 160.0
}, name='Quantity (g)')
trad_amount.index.name = 'Ingredient'
#trad_amount.reset_index()
# Percentage difference the now solution is from the traditional recipe
(((solu_amount - trad_amount) / trad_amount) * 100).dropna().astype(int)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1454877705931965, 'fraction_numerical': 0.025653640708462187, 'mean_word_length': 4.644188813962713, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '475501', 'license': 'MIT', 'n_tokens_mistral': 5465, 'n_tokens_neox': 5033, 'n_words': 1408, 'path': 'notebooks/01-dvg-solver-2020.05.25-Copy1.ipynb', 'repo_name': 'sekondusg/dunli', 'size': '14263'}