Category Archives: Julia

Changes for DynamicHMC.jl and related packages

By: Tamás K. Papp

Re-posted from: https://tamaspapp.eu/post/dynamichmc-reorganization/

TL;DR: new wrapper types for AD (with support for Flux.jl), new transformation interface, see the examples.

DynamicHMC.jl is Julia package for Hamiltonian MCMC using the NUTS algorithm. Its interface and the related packages underwent some reorganization which involves API changes. This post is a brief summary, as I prefer to maintain the actual examples in a repository I can test. You are encouraged to read these, and the documentation for each package (some of which is WIP).

A summary of the changes:

  1. TransformVariables.jl replaces the previous DSL for transforming parameters (which was ContinuousTransformations.jl). It has a different syntax, and should make it easier to specify custom transformations in the future (this is work in progress).

    For example,

    as((μ = asℝ, σ = asℝ₊))

    will specify a transformation of 2 real numbers to a NamedTuple with fields μ (unconstrained) and σ (constrained to be positive, eg a standard deviation).

  2. LogDensityProblems.jl contains the wrapper types for log density functions \(\mathbb{R}^n \to \mathbb{R}\). The most common workflow is making your function map a tuple of parameters to a ::Real (which may be -Inf), define a transformation, and then obtain the gradient via automatic differentiation. Currently ForwardDiff.jl (suboptimal for gradients, but very robust) and Flux.jl (faster) are supported. Support for Zygote.jl and ReverseDiff.jl should be easy to add, just open an issue.

    For example, once you have specified a callable p that takes the transformed parameters which result from transformation t,

    P = TransformedLogDensity(t, p)
    ∇P = FluxGradientLogDensity(P)

    will wrap them so that Flux.jl is used for obtaining the gradient via AD.

  3. DynamicHMC.jl expects a callable that supports the dimension query function, takes a vector of reals, and returns a LogDensityProblems.ValueGradient type. The recommended way of obtaining this is using the wrapper libraries above, but you are of course free to code your own.

    To continue the example, obtaining samples is now as simple as

    chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

    because ∇P now contains the dimension information.

All of the above packages are registered. The new version for DynamicHMC.jl is v1.0.0, as the API change is breaking. All of the packages require at least Julia v0.7, but v1.0 or later is recommended.

DynamicHMCExamples.jl has some examples on suggested usage. I am happy to provide support with coding specific models, just please open an issue there. As usual, bug reports, feature requests and PRs are welcome.

Finally, I would like to thank users for their patience and feedback. I hope that you will find new interface more convenient, and the improved speed useful.

PS.: I have removed the Disqus comments from the blog; they broke the style and hardly anyone was using them anyway.

Infant Formula, Cow’s Milk and Nutrional Blending

We compare the macronutrional content of baby formula with that of cow’s milk and then try to recreate the formula by blending milk with water and adding supplements as needed.

Human milk has a vastly different macronutrient composition compared to that of other mammals, in particular cows. While the energy content is about the same, it is much richer in carbs and a lot lower in protein. The reasoning is that human babies already have a large and active brain for learning and need the sugar to fuel it. At the same time, their body does not need to grow as rapidly as that of other mammals who are already much more mature when born.

For parents who do not (or no longer) breastfeed their baby, a substitute can be used made from a powder, which mimics the nutritional content of human milk. As the baby turns older, and its digestive system is already used to different kinds of foods, cow milk (and other dairy products) may be introduced in the diet.

From other parents, I have heard that they will not feed cow’s milk as is, but dilute it with water during the night. I wondered whether this approach is justified with respect to the macronutrients which led to the following analysis.

This analysis is done purely out of curiosity. Please consult a health care professional prior to changing your infant’s diet!

Data

As a reference for baby formula, we will use Hipp PRE BIO COMBIOTIK. They provide the following data on macronutrients:

(Energy: 66 kcal / 100 ml)
Fat: 3.5 g / 100 ml
Carbs: 7.3 g / 100 ml
Protein: 1.25 g / 100 ml

For cow’s milk, we will use both whole milk as well as reduced-fat milk. The macronutrients as given for milk at LIDL are, for whole milk:

(Energy: 68 kcal / 100 ml)
Fat: 3.9 g / 100 ml
Carbs: 4.9 g / 100 ml
Protein: 3.4 g / 100 ml

And for reduced-fat milk:

(Energy: 48 kcal / 100 ml)
Fat: 1.6 g / 100 ml
Carbs: 4.9 g / 100 ml
Protein: 3.5 g / 100 ml

Model

We setup an optimization model inspired by the classical diet problem, where we want to use four ingredients (water, whole milk, reduced-fat milk, lactose) to match the nutrional profile of baby formula.

My thinking was that a blend of whole milk and reduced-fat milk should get the amount of fat right, while a dilution with water will get the protein level down to the target. Finally, we can add sugar in the form of lactose powder. This is widely available in drug stores, since it is also used as a mild laxative for infants.

In [1]:
using JuMP
using GLPKMathProgInterface
In [2]:
m = Model(solver=GLPKSolverLP())

# our ingredients (without the implicit water)
@variable(m, whole_milk  0)   # in 100 ml
@variable(m, lowfat_milk  0)  # in 100 ml
@variable(m, lactose  0)      # in g

# slack in our macro targets 
@variable(m, slack_fat  0)
@variable(m, slack_carb  0)
@variable(m, slack_protein  0)

# fit liquid volume (100 ml)
@constraint(m, whole_milk + lowfat_milk  1)

# match fat content (g)
@constraint(m, 3.9 * whole_milk + 1.6 * lowfat_milk == 3.5 + slack_fat)

# match carb content (g)
@constraint(m, 4.9 * whole_milk + 4.9 * lowfat_milk + lactose == 7.3  + slack_carb )

# match protein content (g)
@constraint(m, 3.4 * whole_milk + 3.5 * lowfat_milk == 1.25  + slack_protein)

# minimize the mismatch of target
@objective(m, :Min, slack_carb + slack_fat + slack_protein);

status = solve(m)
Out[2]:
:Optimal
In [3]:
@show getvalue(whole_milk) getvalue(lowfat_milk) getvalue(lactose)
@show getvalue(slack_fat) getvalue(slack_carb) getvalue(slack_protein);
getvalue(whole_milk) = 0.8974358974358976
getvalue(lowfat_milk) = 0.0
getvalue(lactose) = 2.902564102564101
getvalue(slack_fat) = 0.0
getvalue(slack_carb) = 0.0
getvalue(slack_protein) = 1.8012820512820515

It seems that we can not get a perfect match since the slack variables are not zero.

In particular, we will overshoot the protein content, or else not meet the fat content. So, with these ingredients, we will not be able to reproduce the baby formula as planned.

Adding Canola Oil

In order to reduce the protein content, we need another non-dairy fat source. Let us use rapeseed oil which is often recommended for preparation of (non-dairy) infant meals because of its neutral taste and rich content of essential fatty acids (omega-3).

In [4]:
m = Model(solver=GLPKSolverLP())

# our ingredients
@variable(m, whole_milk  0)   # in 100 ml
@variable(m, lowfat_milk  0)  # in 100 ml
@variable(m, lactose  0)      # in g
@variable(m, oil  0)  # in g

# no slack variables this time, we are confident in feasibility!

# fit liquid volume (100ml)
@constraint(m, whole_milk + lowfat_milk  1)

# match fat content (g)
@constraint(m, 3.9 * whole_milk + 1.6 * lowfat_milk + oil == 3.5)

# match carb content (g)
@constraint(m, 4.9 * whole_milk + 4.9 * lowfat_milk + lactose == 7.3)

# match protein content (g)
@constraint(m, 3.4 * whole_milk + 3.5 * lowfat_milk == 1.25)

# minimize supplements
@objective(m, :Min, lactose + oil);

status = solve(m)
Out[4]:
:Optimal
In [5]:
@show getvalue(whole_milk) getvalue(lowfat_milk) getvalue(lactose) getvalue(oil);
getvalue(whole_milk) = 0.36764705882352944
getvalue(lowfat_milk) = 0.0
getvalue(lactose) = 5.498529411764705
getvalue(oil) = 2.0661764705882355

The interpretation of this solution is as follows: To reproduce 100 ml of baby formula, we should mix 37 ml of whole milk with 5.5 g of lactose, 2 ml of oil and about 60 ml of water. Reduced-fat milk is not useful but both lactose and oil supplements are necessary.

Of course, this only takes the macronutrients into account, while the sophisticated baby formulas also strive to meet the micronutrient needs of the baby.

Diluting the Drink

Once the infant has reached a certain age, it is not necessary to supply its body with food every 3-4 hours anymore. In particular, they do not really need to drink milk during the night and might only do so out of habit. In fact, it seems that at one year old, the infant does not need even need to get any liquid and will be fine anyways.

At the same time, frequent food intake (in particular drinking) will increase the risk of tooth decay. So it makes sense to avoid night-time feeding.

In order to ensure some peaceful sleep nonetheless, I would like to scale down the energy content of our night-time drink over time. We will end up with just water, but go there in small steps. We will use our ingredients (with supplements) from above but change the targets: Rather then fixing the macronutrient levels, we will only ask for a specific energy content and put a limit on the protein level, since the infant’s kidney might not be ready to handle larger amounts yet.

In [6]:
function dilution(rel_energy=1.0)
    m = Model(solver=GLPKSolverLP())

    # our ingredients
    @variable(m, whole_milk  0)   # in 100 ml
    @variable(m, lowfat_milk  0)  # in 100 ml
    @variable(m, lactose  0)      # in g
    @variable(m, oil  0)  # in g

    # fit liquid volume (100ml)
    @constraint(m, whole_milk + lowfat_milk  1)

    # match energy content (kcal)
    @constraint(m, 68 * whole_milk + 48 * lowfat_milk + 4 * lactose + 9 * oil == rel_energy * 66)
    
    # limit protein content (g)
    @constraint(m, 3.4 * whole_milk + 3.5 * lowfat_milk <= 1.25)

    # minimize supplements
    @objective(m, :Min, lactose + oil);

    status = solve(m)
    sol = getvalue(whole_milk), getvalue(lowfat_milk), getvalue(lactose), getvalue(oil)
    
    return status, sol
end
Out[6]:
dilution (generic function with 2 methods)
In [7]:
status, sol = dilution(1)
Out[7]:
(:Optimal, (0.36764705882352944, 0.0, 0.0, 4.555555555555555))

By using a parameter value of 1, we will get a drink with the same energy content as our baby formula, but with a more relaxed relation of macronutrients.

Again, reduced-fat milk is not used and we find that only an oil supplement is needed but no sugar supplement. This makes sense, since the energy content of fat is higher than that of sugar and our objective here was to minimize the total mass of supplements.

Let us now compute a sequence of solutions with increasingly low energy content, over the course of two weeks.

In [8]:
# iterate over parameter range and collect solutions
r = range(1.0, stop=1e-6, length=14)
sols = []
for rel_nrg in r
    status, sol = dilution(rel_nrg)
    sol = collect(sol)' # tuple to row vector
    push!(sols, sol)
end

# collect all solutions values into one 2d-array
sols = vcat(sols...)
Out[8]:
14×4 Array{Float64,2}:
 0.367647    0.0  0.0  4.55556  
 0.367647    0.0  0.0  3.99145  
 0.367647    0.0  0.0  3.42735  
 0.367647    0.0  0.0  2.86325  
 0.367647    0.0  0.0  2.29915  
 0.367647    0.0  0.0  1.73505  
 0.367647    0.0  0.0  1.17094  
 0.367647    0.0  0.0  0.606842 
 0.367647    0.0  0.0  0.0427396
 0.298643    0.0  0.0  0.0      
 0.223983    0.0  0.0  0.0      
 0.149322    0.0  0.0  0.0      
 0.0746615   0.0  0.0  0.0      
 9.70588e-7  0.0  0.0  0.0      

We can see that we should never use more than 37 ml of whole milk per 100 ml of our drink.
In the beginning, this amount of milk stays the same and the oil supplementation will go down to account for the reduced energy content. Only when there is no more oil to be added, does the amount of milk go down. At this point, we will only mix whole milk and water.

Let us also analyze the macronutrient composition of our dilutions:

In [9]:
whole_milk, oil = sols[:, 1], sols[:, 4]

fat = 3.9 * whole_milk + oil
carbs = 4.9 * whole_milk
protein = 3.4 * whole_milk

energy = 9 * fat + 4 * carbs + 4 * protein

macros = hcat(fat, carbs, protein, energy)
Out[9]:
14×4 Array{Float64,2}:
 5.98938     1.80147     1.25      66.1103    
 5.42528     1.80147     1.25      61.0334    
 4.86118     1.80147     1.25      55.9565    
 4.29707     1.80147     1.25      50.8795    
 3.73297     1.80147     1.25      45.8026    
 3.16887     1.80147     1.25      40.7257    
 2.60477     1.80147     1.25      35.6488    
 2.04067     1.80147     1.25      30.5719    
 1.47656     1.80147     1.25      25.495     
 1.16471     1.46335     1.01539   20.3973    
 0.873532    1.09751     0.761541  15.298     
 0.582356    0.731678    0.507695  10.1987    
 0.29118     0.365841    0.253849   5.09938   
 3.78529e-6  4.75588e-6  3.3e-6     6.62912e-5

Not surprisingly, our drink is relatively high in fat, more so than either whole milk or baby formula.

Similarly, we can compute the relative energy from macros in the above sequence:

In [10]:
macro_nrg = hcat(9 * fat, 4 * carbs, 4 * protein) ./ energy
Out[10]:
14×3 Array{Float64,2}:
 0.815371  0.108998  0.0756312
 0.800013  0.118065  0.0819224
 0.781868  0.128777  0.0893552
 0.760102  0.141626  0.0982713
 0.733511  0.157325  0.109164 
 0.70029   0.176937  0.122773 
 0.657607  0.202135  0.140257 
 0.600748  0.235703  0.163549 
 0.521243  0.28264   0.196117 
 0.513909  0.286969  0.199122 
 0.513909  0.286969  0.199122 
 0.513909  0.286969  0.199122 
 0.513909  0.286969  0.199122 
 0.513909  0.286969  0.199122 

Energy from fat gos down, while the contribution from carbs and protein goes up at the same rate:

In [11]:
macro_nrg[:, 2] ./ macro_nrg[:, 3]
Out[11]:
14-element Array{Float64,1}:
 1.4411764705882357
 1.4411764705882355
 1.4411764705882355
 1.4411764705882357
 1.4411764705882353
 1.4411764705882355
 1.4411764705882355
 1.4411764705882355
 1.4411764705882357
 1.4411764705882353
 1.4411764705882355
 1.4411764705882353
 1.4411764705882353
 1.4411764705882353

That’s it for this analysis. I’m sorry for the lack of visualization!