df = pd.read_csv('vitamin_costs.csv')
vitamins = df.vitamin.values
vitamin_cost = df.set_index('vitamin').to_dict()['cost']
df.head()

"""
Explanation: A new high performance and optimally cheap multivitamin. Saving the world one mixed integer program at a time
You have been tasked with developing a new superior multivitamin. You have been given free reign to select the ingredients and their relative amounts in the vitamin but you have been asked to keep the cost of the raw materials as low as possible.
Further complicating things, you must abide by some restrictions in the formula.

The formula should not have more than 20% of any one vitamin
The formula must have at least 10% Iron, Zinc, or Magnesium
The formula must have at least 20% Vitamin A, Vitamin C, or Vitamin D
Each vitamin that is used must account for at least 5% of the total
The formula may contain as few as 5 vitamins but no more than 10 vitamins
If the formula contains Magnesium it must also contain Calcium and Zinc
The formula must have one of the B vitamins either B6 or B12 but not both

The possible Ingredients and their per mg cost
{{HTML(df.to_html(index=False)) }}
You might see a simple greedy strategy to solve this, but the vitamin B constrains and the Magnesium/Zinc/Calcium constraints make things a bit more complicated. Instead of trial and error we can try and create this vitamin by writing a fairly simple mixed integer linear program. First we need to come up with an expression that captures our ultimate goal, in this case, to minimize the cost of raw materials.
We can calculate the total cost of the formula by adding up the individual cost of each vitamin in it.
$$\text{Total Cost} = \text{ sum over all the vitamins (cost of vitamin) * (percent of vitamin) } $$
First lets load the data we have and organize it so we can easily grab the cost of particular vitamin
End of explanation
"""


u = LpVariable.dicts('percent', vitamins, 0, 1, LpContinuous)

"""
Explanation: Lets model the percent each vitamin is included as by a variable u. so if Vitamin C is included at 15% then u[Vitamin C] = .15. We can use pulp.LpVariable.dicts to return a dictionary of variables that are indexed by the vitamin names. This will make referring to variables later on very easy.
End of explanation
"""


prob = LpProblem('super-awesome-vitamin', LpMinimize)

"""
Explanation: The next thing we need to do is create an instance of the pulp.LpProblem class. This creates a problem variable that will hold the cost and constraints and tells PuLP that we want to minimize our cost. We pass in a name and either LpMinimize or LpMaximize for either minimizing our cost or maximizing it.
End of explanation
"""


cost = lpSum([ u[v]*vitamin_cost[v] for v in vitamins])
prob += cost

"""
Explanation: Now we can define our cost. We will use pulp.lpSum instead of the regular python sum function for efficiency. In this problem it won't make a difference so feel free to experiment. This returns a pulp.LpAffineExpression which we will discuss in detail later on. To add our cost expression to the problem we literally just add it to the problem.
End of explanation
"""


#no more than 20% of any one vitamin
for v in vitamins:
    prob += u[v] <= .2

"""
Explanation: We can now start adding our constraints. The first few are very straightforward. Like with the cost, each constraint just gets added to the problem variable. The key difference is that the constraint will be an inequality.
1: The formula should not have more than 20% of any one vitamin
End of explanation
"""


#The formula must have at least 10% Iron, Zinc, or Magnesium
prob += u['Iron'] + u['Zinc'] + u['Magnesium'] >= .1

"""
Explanation: 2: The formula must have at least 10% Iron, Zinc, or Magnesium
End of explanation
"""


#The formula must have at least 20% Vitamin A, Vitamin C, or Vitamin D
prob += u['Vitamin A'] + u['Vitamin C'] + u['Vitamin D'] >= .2

"""
Explanation: 3: The formula must have at least 20% Vitamin A, Vitamin C, or Vitamin D
End of explanation
"""


#binary variable that captures if this vitamin will be used in the formula
b = LpVariable.dicts('use',vitamins,0,1, LpBinary)

"""
Explanation: These were fairly straightforward, the code for the constraints and the description of the constraint are almost identical. However the next constraint is a bit stranger.
"4. Each vitamin that is used must account for at least 5% of the total"
We basically need the u variables to be at least 5% or 0%. To model this we need to introduce some new variables that will track in simple yes/no manner if the vitamin is included in the final formula. Once we have these variables we will link them with the u variables somehow and satisfy the rest of the constraints. For now lets just look at the code.
End of explanation
"""


#Each vitamin that is used must account for at least 5% of the total
for v in vitamins:
    #if we don't use this vitamin then the percent must be zero
    prob += u[v] <= b[v] 
    #likewise if we do use this vitamin, then the percent must not be zero
    prob += u[v] >= .05 -100*(1-b[v]) # > .05 or > .05 -100

"""
Explanation: 4: Each vitamin that is used must account for at least 5% of the total
End of explanation
"""


#The formula may contain as few as 5 vitamins but no more than 10 vitamins
prob += lpSum([ b[v] for v in vitamins]) >= 5
prob += lpSum([ b[v] for v in vitamins]) <= 10

"""
Explanation: That looks a bit confusing but its actually a very common modeling technique we will explore in detail later on. These new b variables make the rest of the constraints really easy to model.
5: The formula may contain as few as 5 vitamins but no more than 10 vitamins
End of explanation
"""


#If the formula contains Magnesium it must also contain Calcium and Zinc
prob += 2*b['Magnesium'] <= b['Calcium'] + b['Zinc'] 

"""
Explanation: 6: If the formula contains Magnesium it must also contain Calcium and Zinc
End of explanation
"""


#The formula must have one of the B vitamins either B6 or B12 but not both
prob += b['Vitamin B12'] + b['Vitamin B6'] == 1

"""
Explanation: 7: The formula must have one of the B vitamins either B6 or B12 but not both
End of explanation
"""


#the percentages must add up to 100
prob += lpSum([ u[v] for v in vitamins]) == 1.0

"""
Explanation: Finally, while it wasn't stated as a constraint, our u variables are supposed to be percents, so they must add up to 100%
End of explanation
"""


LpStatus[prob.solve()]

print('total cost: $%.2f'%prob.objective.value())

for v in vitamins:
    if value(u[v]) >0:
        print( '%s %.0f%% at unit cost of: $%.2f' %(v, 100*value(u[v]), vitamin_cost[v]))

"""
Explanation: We can now solve the problem and relax knowing we have made the best possible multivitamin (with our highly customized and formalized definition of "best")
End of explanation
"""


