Part 5: 10 foods

Select 10 foods: bread, cheese, eggs, fish, fruit and berries, milk and yoghurt, potatoes, red meat, sugar and sweets, vegetables.

library(data.table)
fd <- read.csv('data/foods_0914.csv', sep = ',')
setDT(fd) # use data.table format
fd

              food intake energy protein   fat carbs sugar alcohol  ghge
            <char>  <num>  <num>   <num> <num> <num> <num>   <int> <num>
 1:          Bread  175.4 10.696   0.091 0.030 0.441 0.002       0 0.001
 2:         Cheese   43.4 13.502   0.217 0.242 0.048 0.002       0 0.010
 3:           Eggs   24.6  6.179   0.130 0.106 0.004 0.000       0 0.002
 4:           Fish   69.5  6.086   0.170 0.075 0.024 0.006       0 0.003
 5: Fruit, berries  171.5  2.729   0.008 0.004 0.134 0.029       0 0.001
 6:  Milk, yoghurt  306.1  1.980   0.036 0.011 0.056 0.010       0 0.001
 7:       Potatoes   67.8  3.791   0.021 0.007 0.178 0.000       0 0.000
 8:       Red meat  117.6  8.342   0.173 0.139 0.014 0.000       0 0.013
 9:  Sugar, sweets   16.9 17.988   0.053 0.178 0.609 0.444       0 0.004
10:     Vegetables  154.6  1.565   0.015 0.008 0.050 0.005       0 0.001
    intake_mean intake_lwr intake_upr
          <num>      <num>      <num>
 1:   188.31866  18.831866  407.77324
 2:    46.59652   4.659652  122.39639
 3:    26.41185   2.641185  101.99699
 4:    74.61885   7.461885  267.87632
 5:   184.13142  18.413142  480.13743
 6:   328.64505  32.864505  966.28731
 7:    72.79364   7.279364  218.48830
 8:   122.50376  12.250376  328.75242
 9:    18.14473   1.814473   69.78742
10:   165.98669  16.598669  365.04187

contrib_pergram <- fd[, .(energy, protein, fat, carbs, sugar, alcohol, ghge)]
contrib_pergram

    energy protein   fat carbs sugar alcohol  ghge
     <num>   <num> <num> <num> <num>   <int> <num>
 1: 10.696   0.091 0.030 0.441 0.002       0 0.001
 2: 13.502   0.217 0.242 0.048 0.002       0 0.010
 3:  6.179   0.130 0.106 0.004 0.000       0 0.002
 4:  6.086   0.170 0.075 0.024 0.006       0 0.003
 5:  2.729   0.008 0.004 0.134 0.029       0 0.001
 6:  1.980   0.036 0.011 0.056 0.010       0 0.001
 7:  3.791   0.021 0.007 0.178 0.000       0 0.000
 8:  8.342   0.173 0.139 0.014 0.000       0 0.013
 9: 17.988   0.053 0.178 0.609 0.444       0 0.004
10:  1.565   0.015 0.008 0.050 0.005       0 0.001

# note: intake from the food.csv is slightly different from the new file
current_diet <- fd$intake
current_diet

 [1] 175.4  43.4  24.6  69.5 171.5 306.1  67.8 117.6  16.9 154.6

# 10 foods contribution (maximum)
const_max_10foods <- t(as.matrix(current_diet)) %*% as.matrix(contrib_pergram)
const_max_10foods

       energy protein    fat    carbs   sugar alcohol   ghge
[1,] 5895.142 77.7671 48.704 153.0605 17.1657       0 3.0957

# exclude sugar, alcohol
const_max_10foods <- const_max_10foods[, c('energy', 'protein', 'fat', 'carbs', 'ghge')]


# target constraint on energy and nutrients
# set lower to be 0.9; upper remain the current max
const_lwrupr <- rbind(const_max_10foods*0.9, const_max_10foods*1)
rownames(const_lwrupr) <- c('lwr', 'upr')
const_lwrupr <- data.table(const_lwrupr)
const_lwrupr

     energy  protein     fat    carbs    ghge
      <num>    <num>   <num>    <num>   <num>
1: 5305.628 69.99039 43.8336 137.7544 2.78613
2: 5895.142 77.76710 48.7040 153.0605 3.09570

Constraints for 3 foods

We compute the constraints for 3 foods together. Note that the consumption for each of the 28 food groups are different, hence we can not use 3/28 times the total energy, protein…; we need to use a weighted average.

For example, the total energy for the 10 foods together should be within [5305, 5895].

Standardize food contribution per gram

Since the range of constraints for 5 categories differ hugely, it could affect the numeric evaluation. We want them to be on comparative scales.

The current solution is to standardize the contribution in each category (e.g. energy) by its original value divided by the standard deviation.

Similarly, the upper and lower limit of the constraints also need to be re-scaled. After rescaling, the target will be on a range of hundreds, rather than 3000 vs 1.8.

contrib_pergram <- contrib_pergram[, c('energy', 'protein', 'fat', 'carbs', 'ghge')]

sd_coef <- apply(contrib_pergram, MARGIN = 2, sd)

contrib_pergram_std <- sweep(contrib_pergram, MARGIN = 2, 1/sd_coef, FUN = '*')
contrib_pergram_std

      energy   protein        fat      carbs      ghge
1  1.9807223 1.1926698 0.35692934 2.14338903 0.2284823
2  2.5003470 2.8440587 2.87923000 0.23329404 2.2848233
3  1.1442486 1.7038140 1.26115033 0.01944117 0.4569647
4  1.1270265 2.2280644 0.89232335 0.11664702 0.6854470
5  0.5053657 0.1048501 0.04759058 0.65127921 0.2284823
6  0.3666632 0.4718254 0.13087409 0.27217638 0.2284823
7  0.7020305 0.2752315 0.08328351 0.86513208 0.0000000
8  1.5448004 2.2673832 1.65377260 0.06804410 2.9702703
9  3.3310800 0.6946318 2.11778074 2.95991819 0.9139293
10 0.2898121 0.1965939 0.09518116 0.24301463 0.2284823

# standardize constraint
# test the previous constraint
const_lwrupr

     energy  protein     fat    carbs    ghge
      <num>    <num>   <num>    <num>   <num>
1: 5305.628 69.99039 43.8336 137.7544 2.78613
2: 5895.142 77.76710 48.7040 153.0605 3.09570

const_lwrupr_std <- sweep(const_lwrupr, MARGIN = 2, 1/sd_coef, FUN = '*')
const_lwrupr_std

     energy   protein      fat    carbs     ghge
1  982.5146  917.3123 521.5166 669.5269 636.5815
2 1091.6828 1019.2359 579.4629 743.9188 707.3128

In the future, different scaling factor can be applied; but it should be a positive number after scaling. Could try dividing the difference between max and min of this variable.

Run optimization

Now we are going to solve the optimization problem using nloptr. The details of the algorithm please refer to part 3.

The objective function remains the same for 10 foods. We want to minimize the difference between the current and optimized diet.

current_diet

 [1] 175.4  43.4  24.6  69.5 171.5 306.1  67.8 117.6  16.9 154.6

# 10 foods, therefore 10 elements
# minimize the deviation from the current intake, on 10 foods
objective <- function(x)
{
  return ( (x[1]- current_diet[1])^2 + 
             (x[2]- current_diet[2])^2 + 
             (x[3]- current_diet[3])^2 + 
             (x[4]- current_diet[4])^2 + 
             (x[5]- current_diet[5])^2 + 
             (x[6]- current_diet[6])^2 + 
             (x[7]- current_diet[7])^2 + 
             (x[8]- current_diet[8])^2 + 
             (x[9]- current_diet[9])^2 + 
             (x[10]- current_diet[10])^2)
}

We optimize using 5 categories (10 inequality constraints):

• energy
• protein
• fat
• carbs
• ghge

We show two examples: one reduces ghge to 90% what the current diet produces; the other reduces ghge to 85%.

Example 1: reduce ghge to 90%

First we reduce the GHGE to 90%.

cstr <- copy(const_lwrupr_std)
cstr$ghge <- cstr$ghge * 0.9 # reduce to 0.9
cstr

     energy   protein      fat    carbs     ghge
1  982.5146  917.3123 521.5166 669.5269 572.9233
2 1091.6828 1019.2359 579.4629 743.9188 636.5815

# use the standardized food contribution per gram 
# contrib_pergram_std

# define the inequality constraints
inequalconstr <- function (x) {
  
  cps <- contrib_pergram_std
  
  constr <- c(
    # energy
    - x[1]*cps$energy[1] - x[2]*cps$energy[2] - x[3]*cps$energy[3] - x[4]*cps$energy[4] - x[5]*cps$energy[5] - x[6]*cps$energy[6] -  x[7]*cps$energy[7] - x[8]*cps$energy[8] - x[9]*cps$energy[9]- x[10]*cps$energy[10] + cstr$energy[1], # lower
    x[1]*cps$energy[1] + x[2]*cps$energy[2] + x[3]*cps$energy[3] + x[4]*cps$energy[4] + x[5]*cps$energy[5] + x[6]*cps$energy[6] + x[7]*cps$energy[7] + x[8]*cps$energy[8] + x[9]*cps$energy[9] + x[10]*cps$energy[10] - cstr$energy[2], # upper
    
    # protein

    - x[1]*cps$protein[1] - x[2]*cps$protein[2] - x[3]*cps$protein[3] - x[4]*cps$protein[4] - x[5]*cps$protein[5] - x[6]*cps$protein[6] -  x[7]*cps$protein[7] - x[8]*cps$protein[8] - x[9]*cps$protein[9]- x[10]*cps$protein[10] + cstr$protein[1], # lower
    x[1]*cps$protein[1] + x[2]*cps$protein[2] + x[3]*cps$protein[3] + x[4]*cps$protein[4] + x[5]*cps$protein[5] + x[6]*cps$protein[6] + x[7]*cps$protein[7] + x[8]*cps$protein[8] + x[9]*cps$protein[9] + x[10]*cps$protein[10] - cstr$protein[2], # upper
    
    # fat

    - x[1]*cps$fat[1] - x[2]*cps$fat[2] - x[3]*cps$fat[3] - x[4]*cps$fat[4] - x[5]*cps$fat[5] - x[6]*cps$fat[6] -  x[7]*cps$fat[7] - x[8]*cps$fat[8] - x[9]*cps$fat[9]- x[10]*cps$fat[10] + cstr$fat[1], # lower
    x[1]*cps$fat[1] + x[2]*cps$fat[2] + x[3]*cps$fat[3] + x[4]*cps$fat[4] + x[5]*cps$fat[5] + x[6]*cps$fat[6] + x[7]*cps$fat[7] + x[8]*cps$fat[8] + x[9]*cps$fat[9] + x[10]*cps$fat[10] - cstr$fat[2], # upper
    
    # carbs

    - x[1]*cps$carbs[1] - x[2]*cps$carbs[2] - x[3]*cps$carbs[3] - x[4]*cps$carbs[4] - x[5]*cps$carbs[5] - x[6]*cps$carbs[6] -  x[7]*cps$carbs[7] - x[8]*cps$carbs[8] - x[9]*cps$carbs[9]- x[10]*cps$carbs[10] + cstr$carbs[1], # lower
    x[1]*cps$carbs[1] + x[2]*cps$carbs[2] + x[3]*cps$carbs[3] + x[4]*cps$carbs[4] + x[5]*cps$carbs[5] + x[6]*cps$carbs[6] + x[7]*cps$carbs[7] + x[8]*cps$carbs[8] + x[9]*cps$carbs[9] + x[10]*cps$carbs[10] - cstr$carbs[2], # upper
    
    # ghge
    - x[1]*cps$ghge[1] - x[2]*cps$ghge[2] - x[3]*cps$ghge[3] - x[4]*cps$ghge[4] - x[5]*cps$ghge[5] - x[6]*cps$ghge[6] -  x[7]*cps$ghge[7] - x[8]*cps$ghge[8] - x[9]*cps$ghge[9]- x[10]*cps$ghge[10] + cstr$ghge[1], # lower
    x[1]*cps$ghge[1] + x[2]*cps$ghge[2] + x[3]*cps$ghge[3] + x[4]*cps$ghge[4] + x[5]*cps$ghge[5] + x[6]*cps$ghge[6] + x[7]*cps$ghge[7] + x[8]*cps$ghge[8] + x[9]*cps$ghge[9] + x[10]*cps$ghge[10] - cstr$ghge[2] # upper
    
  )
  return (constr)
}

Set other parameters

# Initial values
x0 <- fd$intake

# lower and upper bounds of x (10 foods)
lb <- fd$intake_lwr
ub <- fd$intake_upr


opts <- list( "algorithm" = "NLOPT_GN_ISRES",
              "xtol_rel"= 1.0e-15,
              "maxeval"= 160000,
              "tol_constraints_ineq" = rep( 1.0e-10, 10 ))

# run the algorithm
res <- nloptr::nloptr(
  x0          = x0,        # initial value for x
  eval_f      = objective, # objective function
  lb          = lb,        # lower bound for x
  ub          = ub,        # upper bound for x
  eval_g_ineq = inequalconstr, # inequality constraint
  opts        = opts       # options
)

print(res)

Call:
nloptr::nloptr(x0 = x0, eval_f = objective, lb = lb, ub = ub, 
    eval_g_ineq = inequalconstr, opts = opts)


Minimization using NLopt version 2.7.1 

NLopt solver status: 5 ( NLOPT_MAXEVAL_REACHED: Optimization stopped because 
maxeval (above) was reached. )

Number of Iterations....: 160000 
Termination conditions:  xtol_rel: 1e-15    maxeval: 160000 
Number of inequality constraints:  10 
Number of equality constraints:    0 
Current value of objective function:  333.453953531023 
Current value of controls: 174.6512 34.47073 24.18368 66.96662 170.151 304.9251 67.96771 102.1279 15.35428 
153.3486

Now we print out the results in a more readable way

res_diet <- res$solution

# print with the current diet, and change percentage
# also print the boundary, in case it hit boundary

diet_result <- data.frame(
  name = fd$food, # food names
  current = current_diet, 
  new = res_diet, 
  percent_change = round((res_diet - current_diet)/current_diet, 3),
  lower_limit = lb, 
  upper_limit = ub
)
diet_result

             name current       new percent_change lower_limit upper_limit
1           Bread   175.4 174.65125         -0.004   18.831866   407.77324
2          Cheese    43.4  34.47073         -0.206    4.659652   122.39639
3            Eggs    24.6  24.18368         -0.017    2.641185   101.99699
4            Fish    69.5  66.96662         -0.036    7.461885   267.87632
5  Fruit, berries   171.5 170.15104         -0.008   18.413142   480.13743
6   Milk, yoghurt   306.1 304.92508         -0.004   32.864505   966.28731
7        Potatoes    67.8  67.96771          0.002    7.279364   218.48830
8        Red meat   117.6 102.12788         -0.132   12.250376   328.75242
9   Sugar, sweets    16.9  15.35428         -0.091    1.814473    69.78742
10     Vegetables   154.6 153.34856         -0.008   16.598669   365.04187

# verify whether it falls within 
output_newdiet <- t(as.matrix(res_diet)) %*% as.matrix(contrib_pergram_std)
# output_newdiet

# cstr
const_result <- t(rbind(output_newdiet, cstr))
colnames(const_result) <- c('new_diet','const_lwr', 'const_upr')
const_result <- data.table(const_result)
# conditions
const_result[, is_ok := 'Yes']
const_result[new_diet < const_lwr, is_ok := 'beyond lower']
const_result[new_diet > const_upr, is_ok := 'beyond upper']

# relative difference (since we rescaled the targets)
const_result[, relative_dev := 0]
const_result[is_ok == 'beyond lower', relative_dev := round((new_diet - const_lwr)/const_lwr, 3)]
const_result[is_ok == 'beyond upper', relative_dev := round((new_diet - const_upr)/const_upr, 3)]

# print out
const_result

Index: <is_ok>
    new_diet const_lwr const_upr  is_ok relative_dev
       <num>     <num>     <num> <char>        <num>
1: 1034.1342  982.5146 1091.6828    Yes            0
2:  949.5431  917.3123 1019.2359    Yes            0
3:  521.5166  521.5166  579.4629    Yes            0
4:  732.9418  669.5269  743.9188    Yes            0
5:  636.5815  572.9233  636.5815    Yes            0

You should always check how much it actually deviates from the target constraints.

Example 2: reduce ghge to 85%

cstr <- copy(const_lwrupr_std) # this line is important: it keeps a copy of the original values
cstr$ghge <- cstr$ghge * 0.85 # reduce to 0.85
cstr

     energy   protein      fat    carbs     ghge
1  982.5146  917.3123 521.5166 669.5269 541.0943
2 1091.6828 1019.2359 579.4629 743.9188 601.2158

# use the standardized food contribution per gram 
# contrib_pergram_std

# define the inequality constraints
inequalconstr <- function (x) {
  
  cps <- contrib_pergram_std
  
  constr <- c(
    # energy
    - x[1]*cps$energy[1] - x[2]*cps$energy[2] - x[3]*cps$energy[3] - x[4]*cps$energy[4] - x[5]*cps$energy[5] - x[6]*cps$energy[6] -  x[7]*cps$energy[7] - x[8]*cps$energy[8] - x[9]*cps$energy[9]- x[10]*cps$energy[10] + cstr$energy[1], # lower
    x[1]*cps$energy[1] + x[2]*cps$energy[2] + x[3]*cps$energy[3] + x[4]*cps$energy[4] + x[5]*cps$energy[5] + x[6]*cps$energy[6] + x[7]*cps$energy[7] + x[8]*cps$energy[8] + x[9]*cps$energy[9] + x[10]*cps$energy[10] - cstr$energy[2], # upper
    
    # protein

    - x[1]*cps$protein[1] - x[2]*cps$protein[2] - x[3]*cps$protein[3] - x[4]*cps$protein[4] - x[5]*cps$protein[5] - x[6]*cps$protein[6] -  x[7]*cps$protein[7] - x[8]*cps$protein[8] - x[9]*cps$protein[9]- x[10]*cps$protein[10] + cstr$protein[1], # lower
    x[1]*cps$protein[1] + x[2]*cps$protein[2] + x[3]*cps$protein[3] + x[4]*cps$protein[4] + x[5]*cps$protein[5] + x[6]*cps$protein[6] + x[7]*cps$protein[7] + x[8]*cps$protein[8] + x[9]*cps$protein[9] + x[10]*cps$protein[10] - cstr$protein[2], # upper
    
    # fat

    - x[1]*cps$fat[1] - x[2]*cps$fat[2] - x[3]*cps$fat[3] - x[4]*cps$fat[4] - x[5]*cps$fat[5] - x[6]*cps$fat[6] -  x[7]*cps$fat[7] - x[8]*cps$fat[8] - x[9]*cps$fat[9]- x[10]*cps$fat[10] + cstr$fat[1], # lower
    x[1]*cps$fat[1] + x[2]*cps$fat[2] + x[3]*cps$fat[3] + x[4]*cps$fat[4] + x[5]*cps$fat[5] + x[6]*cps$fat[6] + x[7]*cps$fat[7] + x[8]*cps$fat[8] + x[9]*cps$fat[9] + x[10]*cps$fat[10] - cstr$fat[2], # upper
    
    # carbs

    - x[1]*cps$carbs[1] - x[2]*cps$carbs[2] - x[3]*cps$carbs[3] - x[4]*cps$carbs[4] - x[5]*cps$carbs[5] - x[6]*cps$carbs[6] -  x[7]*cps$carbs[7] - x[8]*cps$carbs[8] - x[9]*cps$carbs[9]- x[10]*cps$carbs[10] + cstr$carbs[1], # lower
    x[1]*cps$carbs[1] + x[2]*cps$carbs[2] + x[3]*cps$carbs[3] + x[4]*cps$carbs[4] + x[5]*cps$carbs[5] + x[6]*cps$carbs[6] + x[7]*cps$carbs[7] + x[8]*cps$carbs[8] + x[9]*cps$carbs[9] + x[10]*cps$carbs[10] - cstr$carbs[2], # upper
    
    # ghge
    - x[1]*cps$ghge[1] - x[2]*cps$ghge[2] - x[3]*cps$ghge[3] - x[4]*cps$ghge[4] - x[5]*cps$ghge[5] - x[6]*cps$ghge[6] -  x[7]*cps$ghge[7] - x[8]*cps$ghge[8] - x[9]*cps$ghge[9]- x[10]*cps$ghge[10] + cstr$ghge[1], # lower
    x[1]*cps$ghge[1] + x[2]*cps$ghge[2] + x[3]*cps$ghge[3] + x[4]*cps$ghge[4] + x[5]*cps$ghge[5] + x[6]*cps$ghge[6] + x[7]*cps$ghge[7] + x[8]*cps$ghge[8] + x[9]*cps$ghge[9] + x[10]*cps$ghge[10] - cstr$ghge[2] # upper
    
  )
  return (constr)
}

# Initial values
x0 <- fd$intake

# lower and upper bounds of x (10 foods)
lb <- fd$intake_lwr
ub <- fd$intake_upr


opts <- list( "algorithm" = "NLOPT_GN_ISRES",
              "xtol_rel"= 1.0e-15,
              "maxeval"= 160000,
              "tol_constraints_ineq" = rep( 1.0e-10, 10 ))

# run the algorithm
res <- nloptr::nloptr(
  x0          = x0,        # initial value for x
  eval_f      = objective, # objective function
  lb          = lb,        # lower bound for x
  ub          = ub,        # upper bound for x
  eval_g_ineq = inequalconstr, # inequality constraint
  opts        = opts       # options
)

print(res)

Call:
nloptr::nloptr(x0 = x0, eval_f = objective, lb = lb, ub = ub, 
    eval_g_ineq = inequalconstr, opts = opts)


Minimization using NLopt version 2.7.1 

NLopt solver status: 5 ( NLOPT_MAXEVAL_REACHED: Optimization stopped because 
maxeval (above) was reached. )

Number of Iterations....: 160000 
Termination conditions:  xtol_rel: 1e-15    maxeval: 160000 
Number of inequality constraints:  10 
Number of equality constraints:    0 
Current value of objective function:  1130.07723624611 
Current value of controls: 172.9834 36.42209 30.73556 67.67675 167.417 303.3886 67.60507 86.28209 21.40758 
151.5364

Now we check the results for the second problem.

             name current       new percent_change lower_limit upper_limit
1           Bread   175.4 172.98342         -0.014   18.831866   407.77324
2          Cheese    43.4  36.42209         -0.161    4.659652   122.39639
3            Eggs    24.6  30.73556          0.249    2.641185   101.99699
4            Fish    69.5  67.67675         -0.026    7.461885   267.87632
5  Fruit, berries   171.5 167.41696         -0.024   18.413142   480.13743
6   Milk, yoghurt   306.1 303.38857         -0.009   32.864505   966.28731
7        Potatoes    67.8  67.60507         -0.003    7.279364   218.48830
8        Red meat   117.6  86.28209         -0.266   12.250376   328.75242
9   Sugar, sweets    16.9  21.40758          0.267    1.814473    69.78742
10     Vegetables   154.6 151.53638         -0.020   16.598669   365.04187

Index: <is_ok>
    new_diet const_lwr const_upr  is_ok relative_dev
       <num>     <num>     <num> <char>        <num>
1: 1036.9677  982.5146 1091.6828    Yes            0
2:  932.6577  917.3123 1019.2359    Yes            0
3:  521.5166  521.5166  579.4629    Yes            0
4:  743.9186  669.5269  743.9188    Yes            0
5:  601.2158  541.0943  601.2158    Yes            0