From 5933b00394b566dfed6a67825c2c3090bf5aeafe Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Fri, 19 Jan 2024 11:35:15 +0100 Subject: [PATCH 1/2] add a decent quickstart tutorial --- docs/src/0-quickstart.jl | 274 +++++++++++++++++++++++ docs/src/1-metabolic-modeling.jl | 6 +- docs/src/2-quadratic-optimization.jl | 6 +- docs/src/3-mixed-integer-optimization.jl | 7 +- 4 files changed, 277 insertions(+), 16 deletions(-) create mode 100644 docs/src/0-quickstart.jl diff --git a/docs/src/0-quickstart.jl b/docs/src/0-quickstart.jl new file mode 100644 index 0000000..f66baa9 --- /dev/null +++ b/docs/src/0-quickstart.jl @@ -0,0 +1,274 @@ + +# Copyright (c) 2023, University of Luxembourg #src +# Copyright (c) 2023, Heinrich-Heine University Duesseldorf #src +# #src +# Licensed under the Apache License, Version 2.0 (the "License"); #src +# you may not use this file except in compliance with the License. #src +# You may obtain a copy of the License at #src +# #src +# http://www.apache.org/licenses/LICENSE-2.0 #src +# #src +# Unless required by applicable law or agreed to in writing, software #src +# distributed under the License is distributed on an "AS IS" BASIS, #src +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. #src +# See the License for the specific language governing permissions and #src +# limitations under the License. #src + +# # Quick start +# +# The primary purpose of ConstraintTrees.jl is to make the representation of +# constraint systems neat, and thus make their manipulation easy and +# high-level. In short, the package abstracts the users from keeping track of +# variable and constraint indexes in matrix form, and gives a nice data +# structure that describes the system, while keeping all variable allocation +# and constraint organization completely implicit. +# +# Here we demonstrate the absolutely basic concepts on the "field allocation" +# problem. +# +# ## The problem: Field area allocation +# +# Suppose we have 100 square kilometers of field, 500 kilos of fertilizer and +# 300 kilos of insecticide. We also have a practically infinite supply of wheat +# and barley seeds. If we decide to sow barley, we can make 550🪙 per square +# kilometer of harvest; if we decide to go with wheat instead, we can make +# 350🪙. Unfortunately each square kilometer of wheat requires 6 kilos of +# fertilizer, and 1 kilo of insecticide, whereas each square kilometer of +# barley requires 2 kilos of fertilizer and 4 kilos of insecticide, because +# insects love barley. +# +# How much of our fields should we allocate to wheat and barley to maximize +# our profit? +# +# ## Field area allocation with ConstraintTrees.jl +# +# Let's import the package and start constructing the problem: + +import ConstraintTrees as C + +# Let's name our system `s`. We first need a few variables: + +s = C.variables(keys = [:wheat, :barley]) + +# With ConstraintTrees.jl, we can (and want to!) label everything very nicely +# -- the constraint trees are essentially directory structures, so one can +# prefix everything with symbols to put it into nice directories. + +s = :area^s + +# To be absolutely realistic, we also want to make sure that all areas are +# non-negative. To demonstrate how to do it nicely from the start, we here +# rather re-do the constraints with an interval bound: + +s = :area^C.variables(keys = [:wheat, :barley], bounds = C.Between(0, Inf)) + +# Constraint trees can be browsed using dot notation, or just like +# dictionaries: + +s.area + +# + +s[:area].barley + +# Now let's start rewriting the problem into the constraint-tree-ish +# description. First, we only have 100 square kilometers of area: + +total_area = s.area.wheat.value + s.area.barley.value + +total_area_constraint = C.Constraint(total_area, (0, 100)) + +# We can add any kind of constraint into the existing constraint trees by +# "merging" multiple trees with operator `*`: + +s *= :total_area^total_area_constraint + +# Now let's add constraints for resources. We can create whole +# [`ConstraintTree`](@ref ConstraintTrees.ConstraintTree) structures like +# dictionaries in place, as follows: + +s *= + :resources^C.ConstraintTree( + :fertilizer => + C.Constraint(s.area.wheat.value * 6 + s.area.barley.value * 2, (0, 500)), + :insecticide => + C.Constraint(s.area.wheat.value * 1 + s.area.barley.value * 4, (0, 300)), + ) + +# We can also represent the expected profit as a constraint (although we do not +# need to actually put a constraint bound there): + +s *= :profit^C.Constraint(s.area.wheat.value * 350 + s.area.barley.value * 550) + +# ## Solving the system with JuMP +# +# We can now take the structure of the constraint tree, translate it to any +# suitable linear optimizer interface, and have it solved. For popular reasons +# we choose JuMP with GLPK -- the code is left uncommented here; see the other +# examples for a slightly more detailed explanation: + +import JuMP +function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimizer) + model = JuMP.Model(optimizer) + JuMP.@variable(model, x[1:C.var_count(cs)]) + JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) + C.traverse(cs) do c + b = c.bound + if b isa C.EqualTo + JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to) + elseif b isa C.Between + val = C.substitute(c.value, x) + isinf(b.lower) || JuMP.@constraint(model, val >= b.lower) + isinf(b.upper) || JuMP.@constraint(model, val <= b.upper) + end + end + JuMP.optimize!(model) + JuMP.value.(model[:x]) +end + +import GLPK +optimal_variable_assignment = optimized_vars(s, s.profit.value, GLPK.Optimizer) + +# That gives us the optimized variable values! If we cared to remember what +# they stand for, we might already see how much barley to use. On the other +# hand, the main point of ConstraintTrees is that one should not be forced to +# remember things like variable ordering and indexes, or be forced to manually +# calculate how much money we actually make or how much fertilizer is going to +# be left -- we can simply feed the variable values back to the constraint +# tree, and get a really good overview of all values in our constrained system: + +optimal_s = C.substitute_values(s, optimal_variable_assignment) + +# ## Browsing the result +# +# `optimal_s` is now like the original constraint tree, just the contents are +# "plain old values" instead of the constraints as above. Thus we can easily +# see our profit: + +optimal_s.profit + +@test isapprox(optimal_s.profit, 48333.33333333334) #src + +# The occupied area for each crop: + +optimal_s.area + +# The consumed resources: + +optimal_s.resources + +# ## Increasing the complexity +# +# A crucial property of constraint trees is that the users do not need to care +# about what kind of value they are manipulating -- no matter if something is a +# variable or a derived value, the code that works with it is the same. For +# example, we can use the actual prices for our resources (30🪙 and 110🪙 for a +# kilo of fertilizer and insecticide, respectively) to make a corrected profit: + +s *= + :actual_profit^C.Constraint( + s.profit.value - 30 * s.resources.fertilizer.value - + 110 * s.resources.insecticide.value, + ) + +# Is the result going to change if we optimize for the corrected profit? + +realistically_optimal_s = + C.substitute_values(s, optimized_vars(s, s.actual_profit.value, GLPK.Optimizer)) + +# + +realistically_optimal_s.area + +# ## Combining constraint systems: Let's have a factory! +# +# The second crucial property of constraint trees is the ability to easily +# combine different constraint systems into one. Let's pretend we also somehow +# obtained a food factory that produces malty sweet bread and wheaty +# weizen-style beer, with various extra consumptions of water and heat for each +# of the products. For simplicity, let's just create the corresponding +# constraint system (`f` as a factory) here: + +f = :products^C.variables(keys = [:bread, :weizen], bounds = C.Between(0, Inf)) +f *= :profit^C.Constraint(25 * f.products.weizen.value + 35 * f.products.bread.value) + +# We can make the constraint systems more complex by adding additional +# variables. To make sure the variables do not "conflict", one must use the `+` +# operator. While constraint systems combined with `*` always share variables, +# constraint systems combined with `+` are independent. +f += :materials^C.variables(keys = [:wheat, :barley], bounds = C.Between(0, Inf)) + +# How much resources are consumed by each product, with a limit on each: + +f *= + :resources^C.ConstraintTree( + :heat => C.Constraint( + 5 * f.products.bread.value + 3 * f.products.weizen.value, + (0, 1000), + ), + :water => C.Constraint( + 2 * f.products.bread.value + 10 * f.products.weizen.value, + (0, 3000), + ), + ) + +# How much raw materials are required for each product: + +f *= + :material_allocation^C.ConstraintTree( + :wheat => C.Constraint( + 8 * f.products.bread.value + 2 * f.products.weizen.value - + f.materials.wheat.value, + 0, + ), + :barley => C.Constraint( + 0.5 * f.products.bread.value + 10 * f.products.weizen.value - + f.materials.barley.value, + 0, + ), + ) + +# Having the 2 systems at hand, we can connect the factory "system" `f` to the +# field "system" `s` (into a compound system `c` here): + +c = :factory^f + :fields^s + +#md # !!! warning "Operators for combining constraint trees" +#md # Always remember to use `+` instead of `*` when combining _independent_ constraint systems. If we use `*`, the variables in both systems will become implicitly shared, which is rarely what one wants in the first place. Use `*` only if adding additional constraints to an existing system. As a rule of thumb, one can remember the boolean interpretation of `*` as "and" and of `+` as "or". +#md # +#md # On a side note, the operator `^` was chosen mainly to match the algebraic view of the tree combination, and nicely fit into Julia's operator priority structure. + +# To connect the systems, let's add a transport -- the barley and wheat +# produced on the fields is going to be the only barley and wheat consumed by +# the factory, thus their production and consumption must sum to net zero: + +c *= :transport^C.zip(c.fields.area, c.factory.materials) do area, material + C.Constraint(area.value - material.value, 0) +end + +#md # !!! info "High-level constraint tree manipulation" +#md # There is also a [dedicated example](4-functional-tree-processing.md) with many more useful functions like `zip` above. + +# Finally, let's see how much money can we make from having the factory +# supported by our fields in total! + +optimal_c = + C.substitute_values(c, optimized_vars(c, c.factory.profit.value, GLPK.Optimizer)) + +# How much field area did we allocate? + +optimal_c.fields.area + +# How much of each of the products does the factory make in the end? + +optimal_c.factory.products + +# How much extra resources is consumed by the factory? + +optimal_c.factory.resources + +# And what is the factory profit in the end? + +optimal_c.factory.profit + +@test isapprox(optimal_c.factory.profit, 361.5506329113926) #src diff --git a/docs/src/1-metabolic-modeling.jl b/docs/src/1-metabolic-modeling.jl index 76bb96c..a31ff86 100644 --- a/docs/src/1-metabolic-modeling.jl +++ b/docs/src/1-metabolic-modeling.jl @@ -239,7 +239,7 @@ function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimize model = JuMP.Model(optimizer) JuMP.@variable(model, x[1:C.var_count(cs)]) JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) - function add_constraint(c::C.Constraint) + C.traverse(cs) do c b = c.bound if b isa C.EqualTo JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to) @@ -249,10 +249,6 @@ function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimize isinf(b.upper) || JuMP.@constraint(model, val <= b.upper) end end - function add_constraint(c::C.ConstraintTree) - add_constraint.(values(c)) - end - add_constraint(cs) JuMP.optimize!(model) JuMP.value.(model[:x]) end diff --git a/docs/src/2-quadratic-optimization.jl b/docs/src/2-quadratic-optimization.jl index c216a82..281498c 100644 --- a/docs/src/2-quadratic-optimization.jl +++ b/docs/src/2-quadratic-optimization.jl @@ -113,7 +113,7 @@ function quad_optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer model = JuMP.Model(optimizer) JuMP.@variable(model, x[1:C.var_count(cs)]) JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) - function add_constraint(c::C.Constraint) + C.traverse(cs) do c b = c.bound if b isa C.EqualTo JuMP.@constraint(model, C.substitute(c.value, x) == b.equal_to) @@ -123,10 +123,6 @@ function quad_optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer isinf(b.upper) || JuMP.@constraint(model, val <= b.upper) end end - function add_constraint(c::C.ConstraintTree) - add_constraint.(values(c)) - end - add_constraint(cs) JuMP.set_silent(model) JuMP.optimize!(model) JuMP.value.(model[:x]) diff --git a/docs/src/3-mixed-integer-optimization.jl b/docs/src/3-mixed-integer-optimization.jl index 2801955..2c2d42c 100644 --- a/docs/src/3-mixed-integer-optimization.jl +++ b/docs/src/3-mixed-integer-optimization.jl @@ -14,7 +14,6 @@ # See the License for the specific language governing permissions and #src # limitations under the License. #src - # # Example: Mixed integer optimization (MILP) # # This example demonstrates the extension of `ConstraintTree` bounds structures @@ -78,13 +77,9 @@ function milp_optimized_vars(cs::C.ConstraintTree, objective::C.Value, optimizer model = JuMP.Model(optimizer) JuMP.@variable(model, x[1:C.var_count(cs)]) JuMP.@objective(model, JuMP.MAX_SENSE, C.substitute(objective, x)) - function add_constraint(c::C.Constraint) + C.traverse(cs) do c isnothing(c.bound) || jump_constraint(model, x, c.value, c.bound) end - function add_constraint(c::C.ConstraintTree) - add_constraint.(values(c)) - end - add_constraint(cs) JuMP.set_silent(model) JuMP.optimize!(model) JuMP.value.(model[:x]) From 7d85ceb9e18f86dd20d0cb2d788e985f8922bcc8 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Fri, 19 Jan 2024 14:58:45 +0100 Subject: [PATCH 2/2] fix reviews, add a few links --- docs/src/0-quickstart.jl | 45 ++++++++++++++++++++++------------------ 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/docs/src/0-quickstart.jl b/docs/src/0-quickstart.jl index f66baa9..635ddea 100644 --- a/docs/src/0-quickstart.jl +++ b/docs/src/0-quickstart.jl @@ -46,19 +46,21 @@ import ConstraintTrees as C -# Let's name our system `s`. We first need a few variables: +# Let's name our system `s`. We first need a few [`variables:`](@ref +# ConstraintTrees.variables) s = C.variables(keys = [:wheat, :barley]) # With ConstraintTrees.jl, we can (and want to!) label everything very nicely # -- the constraint trees are essentially directory structures, so one can -# prefix everything with symbols to put it into nice directories. +# prefix everything with symbols to put it into nice directories, e.g. as such: -s = :area^s +:area^s # To be absolutely realistic, we also want to make sure that all areas are -# non-negative. To demonstrate how to do it nicely from the start, we here -# rather re-do the constraints with an interval bound: +# non-negative. To demonstrate how to do that nicely from the start, we rather +# re-do the constraints with an appropriate [interval bound](@ref +# ConstraintTrees.Between): s = :area^C.variables(keys = [:wheat, :barley], bounds = C.Between(0, Inf)) @@ -78,8 +80,8 @@ total_area = s.area.wheat.value + s.area.barley.value total_area_constraint = C.Constraint(total_area, (0, 100)) -# We can add any kind of constraint into the existing constraint trees by -# "merging" multiple trees with operator `*`: +# We can add any kind of [constraint](@ref ConstraintTrees.Constraint) into +# the existing constraint trees by "merging" multiple trees with operator `*`: s *= :total_area^total_area_constraint @@ -104,8 +106,9 @@ s *= :profit^C.Constraint(s.area.wheat.value * 350 + s.area.barley.value * 550) # # We can now take the structure of the constraint tree, translate it to any # suitable linear optimizer interface, and have it solved. For popular reasons -# we choose JuMP with GLPK -- the code is left uncommented here; see the other -# examples for a slightly more detailed explanation: +# we choose [JuMP](https://jump.dev/) with +# [GLPK](https://www.gnu.org/software/glpk/) -- the code is left uncommented +# here as-is; see the other examples for a slightly more detailed explanation: import JuMP function optimized_vars(cs::C.ConstraintTree, objective::C.LinearValue, optimizer) @@ -129,13 +132,14 @@ end import GLPK optimal_variable_assignment = optimized_vars(s, s.profit.value, GLPK.Optimizer) -# That gives us the optimized variable values! If we cared to remember what -# they stand for, we might already see how much barley to use. On the other -# hand, the main point of ConstraintTrees is that one should not be forced to +# This gives us the optimized variable values! If we cared to remember what +# they stand for, we might already know how much barley to sow. On the other +# hand, the main point of ConstraintTree.jl is that one should not be forced to # remember things like variable ordering and indexes, or be forced to manually # calculate how much money we actually make or how much fertilizer is going to -# be left -- we can simply feed the variable values back to the constraint -# tree, and get a really good overview of all values in our constrained system: +# be left -- instead, we can simply [feed the variable values back to the +# constraint tree](@ref ConstraintTrees.substitute_values), and get a really +# good overview of all values in our constrained system: optimal_s = C.substitute_values(s, optimal_variable_assignment) @@ -228,8 +232,8 @@ f *= ), ) -# Having the 2 systems at hand, we can connect the factory "system" `f` to the -# field "system" `s` (into a compound system `c` here): +# Having the two systems at hand, we can connect the factory "system" `f` to +# the field "system" `s`, making a compound system `c` as such: c = :factory^f + :fields^s @@ -238,16 +242,17 @@ c = :factory^f + :fields^s #md # #md # On a side note, the operator `^` was chosen mainly to match the algebraic view of the tree combination, and nicely fit into Julia's operator priority structure. -# To connect the systems, let's add a transport -- the barley and wheat -# produced on the fields is going to be the only barley and wheat consumed by -# the factory, thus their production and consumption must sum to net zero: +# To actually connect the systems (which now exist as completely independent +# parts of `s`), let's add a transport -- the barley and wheat produced on the +# fields is going to be the only barley and wheat consumed by the factory, thus +# their production and consumption must sum to net zero: c *= :transport^C.zip(c.fields.area, c.factory.materials) do area, material C.Constraint(area.value - material.value, 0) end #md # !!! info "High-level constraint tree manipulation" -#md # There is also a [dedicated example](4-functional-tree-processing.md) with many more useful functions like `zip` above. +#md # There is also a [dedicated example](4-functional-tree-processing.md) with many more useful functions like [`zip`](@ref ConstraintTrees.zip) above. # Finally, let's see how much money can we make from having the factory # supported by our fields in total!