Mini-lecture 7-12 - Intro to optimization.ipynb

Kernel: SageMath 9.8

from itertools import product
import numpy as np
import plotly.graph_objects
from ipywidgets import HBox, Layout, interact as _original_interact

%%html
<!-- To disable the modebar IN ALL PLOT.LY PLOTS -->
<style>
.modebar {
    display: none !important;
}
</style>

_original_interact = interact
def interact(_function_to_wrap=None, _layout="horizontal", **kwargs):
    """interact, but with widgets laid out in a horizontal flexbox layout

    This function works exactly like 'interact' (from SageMath or ipywidgets), 
    except that instead of putting all of the widgets into a vertical box 
    (VBox), it uses a horizontal box (HBox) by default. The HBox uses a flexbox 
    layout, so that if there are many widgets, they'll wrap onto a second row. 
    
    Options:
        '_layout' - 'horizontal' by default. Anything else, and it will revert 
        back to using the default layout of 'interact' (a VBox). 
    """
    def decorator(f):
        retval = _original_interact(f, **kwargs)
        if _layout == "horizontal":
            widgets = retval.widget.children[:-1]
            output = retval.widget.children[-1]
            hbox = HBox(widgets, layout=Layout(flex_flow="row wrap"))
            retval.widget.children = (hbox, output)
        return retval
    if _function_to_wrap is None:
        # No function passed in, so this function must *return* a decorator
        return decorator
    # This function was called directly, *or* was used as a decorator directly
    return decorator(_function_to_wrap)

def find_zeros(field, *ranges, **options):
    """Numerically approximate all zeros of a vector field within a box

    Each 'range' should have the form (x, xmin, xmax), where 'x' is one of the
    variables appearing in the vector field, and 'xmin' and 'xmax' are the
    bounds on that variables.

    Options:
        'intervals' - How many subintervals to use in each dimension (default
            20). If there are zeros that are very close together, you may need
            to increase this to find them. But be warned that the running time
            is roughly this to the n power, where n is the number of variables.
        'tolerance' - How close to approximate roots, roughly (default 1e-4)
        'maxiter' - Maximum number of iterations of Newton's Method to run for
            any one (potential) solution. This defaults to -2*log(tolerance).
            There usually isn't much harm in increasing this, unless there are
            many false hits.
        'round' - Round the components of the solutions to this many decimal
            places (default None, meaning do not round them at all)
    """

    # Initialization
    intervals = options.get("intervals", 20)
    tolerance = options.get("tolerance", 1e-5)
    maxiter = options.get("maxiter", int(round(-2*log(tolerance))))
    roundto = options.get("round", None)
    n = len(ranges)
    if len(field) != n:
        raise ValueError("Dimension of vector field is {}, but {} ranges "
                         "given".format(len(field), n))

    mins  = [xmin - (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    maxes = [xmax + (xmax - xmin)/intervals*tolerance for x, xmin, xmax in ranges]
    deltas = [(xmax - xmin) / intervals for xmin, xmax in zip(mins, maxes)]
    powers = [1 << i for i in range(n)]
    J = jacobian(field, [x for x, xmin, xmax in ranges])
    def dist(v, w):
        return sqrt(sum(((a - b) / d)**2 for a, b, d in zip(v, w, deltas)))

    # Set up the array of positive/negative signs of the vector field
    signs = np.zeros((intervals + 1,) * n, dtype=int)
    sranges = [srange(m, m + d*(intervals + 0.5), d) for m, d in zip(mins, deltas)]
    for vertex, index in zip(product(*sranges), product(range(intervals + 1), repeat=n)):
        v = field(*vertex)
        signs[index] = sum(powers[i] for i in range(n) if v[i] > 0)

    # Now search through that array for potential solutions
    solutions = []
    mask = int((1 << n) - 1)
    for index in product(range(intervals), repeat=n):
        indexpowers = list(zip(index, powers))
        all0s = 0
        all1s = mask
        for k in range(mask + 1):
            newindex = [i + (1 if k & p else 0) for i, p in indexpowers]
            vertex = signs[tuple(newindex)]
            all0s |= vertex
            all1s &= vertex
        if all0s & ~all1s == mask:
            # Now do Newton's method!
            v = vector(m + d*(i + 0.5) for i, m, d in zip(index, mins, deltas))
            for i in range(maxiter):
                previous_v = v
                v = v - J(*v).solve_right(field(*v))
                if dist(v, previous_v) < tolerance:
                    break
            else:
                warn("{} iterations reached without convergence for solution "
                     "{}".format(maxiter, v), RuntimeWarning)
            if solutions and min(dist(v, w) for w in solutions) < 2*tolerance:
                continue
            if not all(xmin <= x <= xmax for x, xmin, xmax in zip(v, mins, maxes)):
                continue
            solutions.append(vector(RDF, v))

    # A convenience: round to some number of decimal places, if requested
    if roundto is not None:
        solutions = [vector(round(x, roundto) for x in v) for v in solutions]
    return solutions

def find_zeros1d(f, x_range, **options):
    x, xmin, xmax = x_range
    intervals = options.get("intervals", 20)
    tolerance = options.get("tolerance", 1e-5)
    f = fast_float(f, x)

    x_min = xmin - (xmax - xmin)/intervals*tolerance
    x_max = xmax + (xmax - xmin)/intervals*tolerance
    delta = (xmax - xmin) / intervals
    x_values = np.linspace(x_min, x_max, intervals + 1)
    f_values = np.array([int(f(x0) > 0) for x0 in x_values])
    sign_changes = (f_values[:-1] - f_values[1:]).nonzero()[0]
    solutions = []
    for i in sign_changes:
        x0 = find_root(f, x_values[i], x_values[i+1])
        if solutions and min([abs(x0 - a) for a in solutions]) < 2*tolerance:
            continue
        solutions.append(x0)
    return solutions

# Our mixin class for Figure and FigureWidget. Should never be instantiated! 
class MyFigure(object):
    def add(self, *items, subplot=None):
        if subplot is not None:
            try:
                subplot[0]
            except:
                subplot = (1, subplot)
            row = int(subplot[0])
            col = int(subplot[1])
        retval = []
        text_indices = []
        text3d_indices = []
        for item in items:
            if isinstance(item, plotly.graph_objects.layout.Annotation):
                if subplot is None:
                    self.add_annotation(item)
                else:
                    self.add_annotation(item, row=row, col=col)
                text_indices.append(len(retval))
                retval.append(None)
            elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):
                self.layout.scene.annotations += (item,)
                text3d_indices.append(len(retval))
                retval.append(None)
            else:
                if subplot is None:
                    self.add_trace(item)
                else:
                    self.add_trace(item, row=row, col=col)
                retval.append(self.data[-1])
        for i, pos in enumerate(text_indices, start=-len(text_indices)):
            retval[pos] = self.layout.annotations[i]
        for i, pos in enumerate(text3d_indices, start=-len(text3d_indices)):
            retval[pos] = self.layout.scene.annotations[i]
        if len(retval) == 1:
            return retval[0]
        return retval

    def __iadd__(self, item):
        if isinstance(item, plotly.graph_objects.layout.Annotation):
            self.add_annotation(item)
        elif isinstance(item, plotly.graph_objects.layout.scene.Annotation):
            self.layout.scene.annotations += (item,)
        else:
            self.add_trace(item)
        return self

    def axes_labels(self, *labels):
        if len(labels) == 2:
            self.layout.xaxis.title.text = labels[0]
            self.layout.yaxis.title.text = labels[1]
        elif len(labels) == 3:
            self.layout.scene.xaxis.title.text = labels[0]
            self.layout.scene.yaxis.title.text = labels[1]
            self.layout.scene.zaxis.title.text = labels[2]
        else:
            raise ValueError("You must specify labels for either 2 or 3 axes.")

    def axes_ranges(self, *ranges, scale=None):
        if len(ranges) == 2:
            (xmin, xmax), (ymin, ymax) = ranges
            self.layout.xaxis.range = (xmin, xmax)
            self.layout.yaxis.range = (ymin, ymax)
            if scale is not None:
                x, y = scale
                self.layout.xaxis.constrain = "domain"
                self.layout.yaxis.constrain = "domain"
                self.layout.yaxis.scaleanchor = "x"
                self.layout.yaxis.scaleratio = y / x
        elif len(ranges) == 3:
            (xmin, xmax), (ymin, ymax), (zmin, zmax) = ranges
            self.layout.scene.xaxis.range = (xmin, xmax)
            self.layout.scene.yaxis.range = (ymin, ymax)
            self.layout.scene.zaxis.range = (zmin, zmax)
            if isinstance(scale, str):
                self.layout.scene.aspectmode = scale
            elif scale is not None:
                x, y, z = scale
                x *= xmax - xmin
                y *= ymax - ymin
                z *= zmax - zmin
                c = sorted((x, y, z))[1]
                self.layout.scene.aspectmode = "manual"
                self.layout.scene.aspectratio.update(x=x/c, y=y/c, z=z/c)
        else:
            raise ValueError("You must specify ranges for either 2 or 3 axes.")


class Figure(plotly.graph_objects.Figure, MyFigure):
    def __init__(self, *args, **kwargs):
        specs = kwargs.pop("subplots", None)
        if specs is None:
            super().__init__(*args, **kwargs)
        else:
            try:
                specs[0][0]
            except:
                specs = [specs]
            rows = len(specs)
            cols = len(specs[0])
            fig = plotly.subplots.make_subplots(rows=rows, cols=cols, 
                                                specs=specs, **kwargs)
            super().__init__(fig)


class FigureWidget(plotly.graph_objects.FigureWidget, MyFigure):
    def __init__(self, *args, **kwargs):
        self._auto_items = []
        specs = kwargs.pop("subplots", None)
        if specs is None:
            super().__init__(*args, **kwargs)
        else:
            try:
                specs[0][0]
            except:
                specs = [specs]
            rows = len(specs)
            cols = len(specs[0])
            fig = plotly.subplots.make_subplots(rows=rows, cols=cols, 
                                                specs=specs, **kwargs)
            super().__init__(fig)

    def auto_update(self, *items):
        if self._auto_items:
            for item, newitem in zip(self._auto_items, items):
                if isinstance(newitem, tuple):
                    newitem = newitem[0]
                item.update(newitem)
        else:
            for item in items:
                if isinstance(item, tuple):
                    item, subplot = item
                    self._auto_items.append(self.add(item, subplot=subplot))
                else:
                    self._auto_items.append(self.add(item))

    def initialized(self):
        return len(self._auto_items) > 0


# Below are all the actual plotting methods. First, the 2D graphics:
def plotly_text(text, location, **options):
    options = options.copy()
    x, y = np.array(location, dtype=float)
    if options.pop("update", False):
        return dict(text=text, x=x, y=y)
    options.setdefault("font_color", options.pop("color", "black"))
    size = options.pop("size", None)
    if size is not None:
        options.setdefault("font_size", size)
    arrow = options.pop("arrow", None)
    if arrow:
        options.setdefault("ax", float(arrow[0]))
        options.setdefault("ay", float(arrow[1]))
        options.setdefault("showarrow", True)
    else:
        options.setdefault("showarrow", False)
    if options.pop("paper", False):
        options.update(xref="paper", yref="paper")
        options.setdefault("xanchor", "left")
        options.setdefault("yanchor", "bottom")
    return plotly.graph_objects.layout.Annotation(
            text=text, x=x, y=y, **options)


def plotly_points(points, **options):
    options = options.copy()
    x, y = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    options.setdefault("marker_color", options.pop("color", "black"))
    options.setdefault("marker_size", options.pop("size", 8))
    options.setdefault("mode", "markers")
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_lines(points, **options):
    options = options.copy()
    x, y = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    color = options.pop("color", "blue")
    options.setdefault("line_color", color)
    options.setdefault("mode", "lines")
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_function(f, x_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 81)
    color = options.pop("color", "blue")
    options.setdefault("line_color", color)
    options.setdefault("mode", "lines")
    options.setdefault("line_shape", "spline")
    options.setdefault("line_smoothing", 1.3)
    x, xmin, xmax = x_range
    f = fast_float(f, x)
    x = np.linspace(float(xmin), float(xmax), plotpoints)
    y = np.array([f(x0) for x0 in x])
    if options.pop("update", False):
        return dict(x=x, y=y)
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_parametric(f, t_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 101)
    color = options.pop("color", "blue")
    options.setdefault("line_color", color)
    options.setdefault("mode", "lines")
    options.setdefault("line_shape", "spline")
    options.setdefault("line_smoothing", 1.3)
    t, tmin, tmax = t_range
    f1, f2 = [fast_float(f_, t) for f_ in f]
    t = np.linspace(float(tmin), float(tmax), plotpoints)
    x, y = np.array([(f1(t0), f2(t0)) for t0 in t]).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y)
    return plotly.graph_objects.Scatter(x=x, y=y, **options)


def plotly_points3d(points, **options):
    options = options.copy()
    options.setdefault("marker_color", options.pop("color", "black"))
    options.setdefault("marker_size", options.pop("size", 2.5))
    options.setdefault("mode", "markers")
    x, y, z = np.array(points, dtype=float).transpose()
    if options.pop("update", False):
        return dict(x=x, y=y, z=z)
    return plotly.graph_objects.Scatter3d(x=x, y=y, z=z, **options)


def plotly_function3d(f, x_range, y_range, **options):
    options = options.copy()
    plotpoints = options.pop("plotpoints", 81)
    try:
        plotpointsx, plotpointsy = plotpoints
    except:
        plotpointsx = plotpointsy = plotpoints
    color = options.pop("color", "lightblue")
    options.setdefault("colorscale", (color, color))
    x, xmin, xmax = x_range
    y, ymin, ymax = y_range
    f = fast_float(f, x, y)
    x = np.linspace(float(xmin), float(xmax), plotpointsx)
    y = np.linspace(float(ymin), float(ymax), plotpointsy)
    x, y = np.meshgrid(x, y)
    xy = np.array([x.flatten(), y.flatten()]).transpose()
    z = np.array([f(x0, y0) for x0, y0 in xy]).reshape(plotpointsy, plotpointsx)
    if options.pop("update", False):
        return dict(x=x, y=y, z=z)
    return plotly.graph_objects.Surface(x=x, y=y, z=z, **options)

Learning goals:

• Be able to explain the significance of optimization in biology, and give several examples.

• Be able to describe the main biological process that underlies all optimization problems in biology.

• Be able to distinguish local maxima and local minima of a function from global extrema.

• Know the significance of local maxima in evolution.

Example 1: Optimal foraging

Consider a bird or small mammal that has a nest, and forages for food in some area around its nest.

• The larger that area, the more food it can gather, so the more energy it gains, or the more energy it is able to provide for its offspring.

• But the larger the area, the more energy it must spend defending its territory from competitors.

What should be optimized here?

def foraging_interactive():
    energy_in(x) = 0.0001 * pi*x^2
    energy_out(x) = 0.000001 * x^3

    figure = FigureWidget(subplots=[[{"rowspan": 2}, {}], [None, {}]])
    figure.axes_ranges((-322, 322), (-322, 350), scale=(1, 1))
    figure.layout.xaxis1.visible = False
    figure.layout.yaxis1.visible = False
    figure.layout.yaxis2.range = [-2, 33]
    figure.layout.xaxis3.range = [0, 320]
    figure.layout.yaxis3.range = [-1, 5]
    figure.layout.yaxis2.title.text = r"Energy (kcal/day)"
    figure.layout.xaxis3.title.text = r"Foraging radius (meters)"
    figure.layout.yaxis3.title.text = r"Net energy gain<br>(kcal/day)"
    figure.layout.showlegend = False
    figure.add(plotly_text("Foraging area", (0,350), font_size=18), subplot=(1,1))
    figure.add(plotly_text("Nest", (0,-2), font_size=12, color="saddlebrown", yanchor="top"), subplot=(1,1))
    color1 = plotly.colors.sequential.Viridis[6]
    color2 = plotly.colors.sequential.Viridis[0]
    color3 = plotly.colors.sequential.Viridis[9]
    labels = ["Energy from food<br>(∝ area)", "Energy spent<br>defending territory", "Net energy gain"]
    energy_chart = plotly.graph_objects.Bar(x=labels, orientation="v", 
            marker_color=[color1, color3, color2], marker_colorscale="Viridis")
    energy_chart = figure.add(energy_chart, subplot=(1,2))
    netenergy_graph = plotly_function(energy_in - energy_out, (x, 0, 1), color=color2)
    netenergy_graph = figure.add(netenergy_graph, subplot=(2,2))
    maximum = find_root((energy_in - energy_out).derivative(), 10, 300)
    maxline = plotly_lines([(maximum, -0.5), (maximum, 6)], color=color1, line_dash="dash")
    maxline = figure.add(maxline, subplot=(2, 2))
    maxlabel = plotly_text(fr"$r = {round(maximum, 1)}\,\text{{m}}$", 
            (maximum - 0.5, 2), font_size=14, xanchor="right", xref="x3", yref="y3")
    maxlabel = figure.add(maxlabel)

    @interact(r=slider(1, 320, 1, default=100, label="Radius"), 
              show_graph=checkbox(False, label="Show graph"), 
              show_max=checkbox(False, label="Show maximum"))
    def update(r, show_graph, show_max):
        initialized = figure.initialized()
        f(t) = (r*cos(t), r*sin(t))
        energy = np.array([energy_in(r), energy_out(r), energy_in(r) - energy_out(r)], dtype=float)
        forage_area = plotly_parametric(f, (t, 0, 2*pi), plotpoints=41, 
                fill="toself", color=color1, fillcolor=color1, update=initialized)
        nest = plotly_points([(0,0)], color="saddlebrown")
        netenergy_update = plotly_function(energy_in - energy_out, (x, 0, r), update=True)
        with figure.batch_update():
            figure.auto_update(forage_area, nest)
            energy_chart.update(y=energy)
            netenergy_graph.update(netenergy_update)
            figure.layout.xaxis3.visible = show_graph
            figure.layout.yaxis3.visible = show_graph
            netenergy_graph.visible = show_graph
            maxline.visible = show_graph and show_max and r > maximum
            maxlabel.visible = show_graph and show_max and r > maximum
            if show_graph:
                figure.layout.yaxis2.domain = [0.42, 1]
                figure.layout.yaxis3.domain = [0, 0.28]
            else:
                figure.layout.yaxis2.domain = [0.1, 1]
                figure.layout.yaxis3.domain = [0, 0.1]

    return figure

foraging_interactive()

Interactive function <function foraging_interactive.<locals>.update at 0x7f01cfd71550> with 3 widgets
  r: Tra…

FigureWidget({
    'data': [{'marker': {'color': [#35b779, #fde725, #440154],
                         'colors…

Example 2: Optimal clutch size

Birds reproduce annually at a certain time of year. The group of eggs a bird lays at one time is called its clutch. So the number of eggs laid is called its clutch size.

• The survival rate is the probability that any one egg will survive to adulthood, until it too can reproduce.

• We can also think of the survival rate as the fraction of eggs from a clutch that will survive, on average. (expected value)

• The more eggs in a clutch, the more difficult it is for the parents to provide for them, protect them, etc. So as the clutch size increases, the survival rate decreases.

What should be optimized here?

def clutchsize_interactive():
    survival(x) = (1 - x/9)^2

    figure = FigureWidget(subplots=[[{"type": "pie"}, {}], [{}, {}]])
    figure.layout.yaxis1.range = [0, 9.5]
    figure.layout.xaxis2.range = [0, 9.5]
    figure.layout.yaxis2.range = [0, 1]
    figure.layout.yaxis2.rangemode = "tozero"
    figure.layout.yaxis2.domain = [0, 0.28]
    figure.layout.xaxis3.range = [0, 9.5]
    figure.layout.yaxis3.range = [0, 1.5]
    figure.layout.yaxis1.title.text = r"# of offspring"
    figure.layout.xaxis2.title.text = r"Clutch size"
    figure.layout.yaxis2.title.text = r"Survival rate"
    figure.layout.xaxis3.title.text = r"Clutch size"
    figure.layout.yaxis3.title.text = r"Surviving offspring"
    figure.layout.showlegend = False
    figure.add(plotly_text("Survival rate<br>(fraction of offspring that survive)", (0.22,0.42), 
            font_size=18, paper=True, xanchor="center", yanchor="top"))
    color1 = plotly.colors.sequential.Viridis[6]
    color2 = plotly.colors.sequential.Viridis[0]
    color3 = plotly.colors.sequential.Viridis[9]
    labels = ["Don't survive", "Survive"]
    survival_chart = plotly.graph_objects.Pie(labels=labels, sort=False, 
            marker_colors=[color2, color1], textinfo="label+percent")
    survival_chart.domain = {"x": [0, 0.45], "y": [0.44, 0.95]}
    survival_chart = figure.add(survival_chart)
    survival_graph = plotly_function(survival, (x, 0, 9), color=color1)
    survival_graph = figure.add(survival_graph, subplot=(2,1))
    labels = ["Clutch size", "Offspring<br>that survive"]
    offspring_chart = plotly.graph_objects.Bar(x=labels, orientation="v", 
            marker_color=[color3, color1])
    offspring_chart = figure.add(offspring_chart, subplot=(1,2))
    netoffspring_graph = plotly_function(x, (x, 0, 1), color=color1)
    netoffspring_graph = figure.add(netoffspring_graph, subplot=(2,2))
    maximum = find_root(diff(x * survival(x), x), 1, 8)
    maxline = plotly_lines([(maximum, -1), (maximum, 2)], color=color2, line_dash="dash")
    maxline = figure.add(maxline, subplot=(2, 2))
    maxlabel = plotly_text(fr"{int(maximum)} eggs", 
            (maximum - 0.1, 0.7), font_size=14, xanchor="right", xref="x3", yref="y3")
    maxlabel = figure.add(maxlabel)

    @interact(C=slider(0, 9, default=0, label="Clutch size"), 
              show_graph=checkbox(False, label="Show graph"), 
              show_max=checkbox(False, label="Show maximum"))
    def update(C, show_graph, show_max):
        survival_pie = np.array([1 - survival(C), survival(C)], dtype=float)
        offspring = np.array([C, C*survival(C)], dtype=float)
        netoffspring_update = plotly_function(x*survival(x), (x, 0, C), update=True)
        with figure.batch_update():
            survival_chart.update(values=survival_pie)
            offspring_chart.update(y=offspring)
            netoffspring_graph.update(netoffspring_update)
            figure.layout.xaxis2.visible = show_graph
            figure.layout.yaxis2.visible = show_graph
            survival_graph.visible = show_graph
            figure.layout.xaxis3.visible = show_graph
            figure.layout.yaxis3.visible = show_graph
            netoffspring_graph.visible = show_graph
            maxline.visible = show_graph and show_max and C > maximum
            maxlabel.visible = show_graph and show_max and C > maximum
            if show_graph:
                figure.layout.yaxis1.domain = [0.42, 1]
                figure.layout.yaxis3.domain = [0, 0.28]
            else:
                figure.layout.yaxis1.domain = [0.1, 1]
                figure.layout.yaxis3.domain = [0, 0.1]

    return figure

clutchsize_interactive()

Interactive function <function clutchsize_interactive.<locals>.update at 0x7f01cfc3a790> with 3 widgets
  C: T…

FigureWidget({
    'data': [{'domain': {'x': [0, 0.45], 'y': [0.44, 0.95]},
              'labels': [Don't sur…

Example 3: Optimal vascular branching

Major arteries branch off into smaller arteries, which branch into arterioles, then smaller arterioles, etc, on down to the level of capillaries, the thinnest blood vessels.

• In each artery, there is some vascular resistance: essentially friction that pushes against the blood flowing through the vessels.

• In thinner arteries, the resistance is much higher per unit of length, so it's advantageous to keep these thinner arteries shorter.

• But blood still needs to reach tissues that are not directly along the wider arteries. Making the thinner arteries as short as possible requires making the wider arteries longer, resulting in more vascular resistance in the wider artery.

What should be optimized here?

def vascular_interactive():
    k = 1e7
    r1 = 60
    r2 = 48
    d = 1.6
    p = 0.3
    distance(x) = d - p * cot(x*pi/180)
    resistance1(x) = k/r1^4 * distance(x)
    resistance2(x) = k/r2^4 * p * csc(x*pi/180)
    resistance = resistance1 + resistance2

    figure = FigureWidget(subplots=[[{"colspan": 2}, None], [{}, {}]])
    figure.layout.margin = dict(t=5, l=60)
    figure.axes_ranges((0, 1.4*d), (-0.3*p, 1.35*p), scale=(1,1))
    figure.layout.xaxis1.visible = False
    figure.layout.yaxis1.visible = False
    figure.layout.yaxis2.range = [0, 2.5]
    figure.layout.xaxis3.range = [15, 90]
    figure.layout.yaxis3.range = [1.6, 2.5]
    figure.layout.yaxis2.title.text = r"Resistance (HRU)"
    figure.layout.xaxis3.title.text = r"Branching angle (degrees)"
    figure.layout.yaxis3.title.text = r"Resistance (HRU)"
    figure.layout.showlegend = False
    figure.add(plotly_text("The smaller blood vessel is 80% as wide as the main artery.", 
            (0, 0.90), paper=True, font_size=14))
    figure.add(plotly_text("Vascular resistance", (0.22, 0.44), paper=True, 
            xanchor="center", font_size=18))
    title = plotly_text("Vascular resistance", (0.72, 0.44), paper=True, 
            xanchor="center", font_size=18)
    color1 = plotly.colors.sequential.Viridis[6]
    color2 = plotly.colors.sequential.Viridis[9]
    color3 = plotly.colors.sequential.Viridis[0]
    figure.add(plotly_lines([(0,0), (2*d,0)], line_width=r1, color="darkred"))
    figure.add(plotly_text("Main artery", (0.15*d, 0.15*p), color="white", font_size=16))
    labels = ["Main artery<br>to branch point", "Smaller blood<br>vessel", "Total<br>resistance"]
    resistance_chart = plotly.graph_objects.Bar(x=labels, orientation="v", 
            marker_color=[color1, color2, color3])
    resistance_chart = figure.add(resistance_chart, subplot=(2,1))
    resistance_graph = plotly_function(x, (x, 0, 1), color=color3)
    resistance_graph = figure.add(resistance_graph, subplot=(2,2))
    tissue(t) = (d*(1 + 0.06*cos(t) + 0.01*cos(11*t)), p*(1.05 + 0.24*sin(t) + 0.04*sin(11*t)))
    tissue = plotly_parametric(tissue, (t, 0, 2*pi), fill="toself", color=color3, fillcolor=color3)
    figure.add(plotly_text("Tissue", (d, 1.05*p), color="white", font_size=16))
    label = plotly_text("Smaller<br>artery", (0,0), color="white", font_size=16)
    minimum = arccos((r2 / r1)^4) * 180 / pi
    minline = plotly_lines([(minimum, 0), (minimum, 3)], color=color1, line_dash="dash")
    minline = figure.add(minline, subplot=(2, 2))
    minlabel = plotly_text(fr"$\theta = {round(minimum, 2)}^\circ$", 
            (minimum - 0.5, 2.3), font_size=14, xanchor="right", xref="x3", yref="y3")
    title, label, minlabel = figure.add(title, label, minlabel)

    @interact(theta=slider(15, 89, default=15, label="Branch angle"), 
              show_graph=checkbox(False, label="Show graph"), 
              show_min=checkbox(False, label="Show minimum"))
    def update(theta, show_graph, show_min):
        initialized = figure.initialized()
        values = np.array([resistance1(theta), resistance2(theta), resistance(theta)], dtype=float)
        graph_update = plotly_function(resistance1(x) + resistance2(x), (x, 10, theta), update=True)
        branch = plotly_lines([(distance(theta),0), (d,p)], line_width=r2, 
                color="darkred", update=initialized)
        tobranchpoint = plotly_lines([(0,0), (distance(theta),0)], line_width=2, 
                color=color1, update=initialized)
        frombranchpoint = plotly_lines([(distance(theta),0), (d,p)], line_width=2, 
                color=color2, update=initialized)
        with figure.batch_update():
            figure.auto_update(branch, tobranchpoint, frombranchpoint, tissue)
            label.update(x=float(0.5*distance(theta) + 0.5*d), y=float(0.5*p), textangle=-theta)
            resistance_chart.update(y=values)
            resistance_graph.update(graph_update)
            figure.layout.xaxis3.visible = show_graph
            figure.layout.yaxis3.visible = show_graph
            resistance_graph.visible = show_graph
            title.visible = show_graph
            minline.visible = show_graph and show_min
            minlabel.visible = show_graph and show_min

    return figure

vascular_interactive()

Interactive function <function vascular_interactive.<locals>.update at 0x7f01cfc64ca0> with 3 widgets
  theta:…

FigureWidget({
    'data': [{'line': {'color': 'darkred', 'width': 60},
              'mode': 'lines',
       …

What's the biological mechanism underlying all of these?

The bird doesn't do a calculus problem to decide how large its foraging radius should be, or how many eggs it should lay.

The cells don't solve math problems when arranging themselves to form arteries, in order to minimize the total vascular resistance.

So what's actually behind all this?

Evolution!

More precisely, natural selection: survival of the fittest.

def fitness_interactive():
    fitness(x) = 0.79 / (1 + ((x - 60)/30)^2) * (1 + (x - 60)/100)

    figure = FigureWidget(subplots=[{}, {}])
    figure.axes_ranges((0, 100), (0, 60), scale=(1,1))
    figure.axes_labels("Body length (cm)", "")
    figure.layout.margin.r = 100
    figure.layout.xaxis1.zeroline = False
    figure.layout.yaxis1.visible = False
    figure.layout.xaxis2.range = [0, 100]
    figure.layout.yaxis2.range = [0, 1]
    figure.layout.xaxis2.title.text = "Body length (cm)"
    figure.layout.yaxis2.title.text = "Fitness"
    figure.layout.showlegend = False
    figure.add(plotly_text("Fitness graph", (0.78,0.95), font_size=18, paper=True, xanchor="center"))
    figure.add_layout_image(source="images/FishImage.png", x=0, y=0, 
            xref="x", yref="y", xanchor="left", yanchor="bottom")
    fish = figure.layout.images[0]
    color1 = plotly.colors.sequential.Viridis[6]
    color2 = plotly.colors.sequential.Viridis[0]
    figure.add(plotly_function(fitness, (x, 5, 100), color=color1), subplot=2)
    point_on_graph = figure.add(plotly_points([(0, 0)]), subplot=2)
    maximum = find_root(fitness.derivative(), 10, 100)
    maxline = plotly_lines([(maximum, -0.1), (maximum, 0.9)], color=color2, line_dash="dash")
    maxline = figure.add(maxline, subplot=2)
    maxlabel = plotly_text(fr"{round(maximum, 1)} cm", 
            (maximum - 0.5, 0.4), font_size=14, xanchor="right", xref="x2", yref="y2")
    maxlabel = figure.add(maxlabel)

    @interact(l=slider(5, 100, 1, default=20, label="Body length"), 
              show_max=checkbox(False, label="Show maximum"))
    def update(l, show_max):
        with figure.batch_update():
            fish.update(sizex=l, sizey=l)
            point_on_graph.update(x=[float(l)], y=[float(fitness(l))])
            maxline.visible = show_max
            maxlabel.visible = show_max

    return figure

fitness_interactive()

Interactive function <function fitness_interactive.<locals>.update at 0x7f01cfc610d0> with 2 widgets
  l: Tran…

FigureWidget({
    'data': [{'line': {'color': '#35b779', 'shape': 'spline', 'smoothing': 1.30000000000000},
 …

Fitness is a function of genotype

(or phenotype)

Definition: The fitness of a particular genotype is the average number (expected value) of offspring that an individual with that exact genotype will have.

In other words, for any combination of genes, fitness measures how much an individual with those genes is likely to contribute to the gene pool of the next generation.

Example 4: Darwin's finches

finches = {
    (0.85, 0.15): "images/DarwinsFinches1.png", 
    (0.40, 0.25): "images/DarwinsFinches2.png", 
    (0.15, 0.65): "images/DarwinsFinches3.png", 
    (0.50, 0.85): "images/DarwinsFinches4.png", 
}
var("x, y")
bump(x0, y0, h) = h^2/(h^2 + (x - x0)^2 + (y - y0)^2)
fitness(x, y) = 0
genotypes = np.array(list(finches.keys()), dtype=float)
for (x0, y0), a in zip(genotypes, (0.7, 0.3, 0.35, 0.6)):
    fitness += a * bump(x0, y0, 0.15 + 0.08*random())
gradient = fitness.gradient()
hessian = jacobian(gradient, (x, y))
crit_points = find_zeros(gradient, (x, 0, 1), (y, 0, 1))
crit_points = [pt for pt in crit_points if all(l < 0 for l in hessian(*pt).eigenvalues())]
closest = lambda pt: np.linalg.norm(genotypes - pt, axis=1).argmin()
crit_points.sort(key=closest)
for oldpt, newpt in zip(genotypes, crit_points):
    finches[tuple(newpt)] = finches.pop(tuple(oldpt))

def finches_interactive():
    figure = FigureWidget(subplots=[{}, {"type": "scene"}])
    figure.layout.yaxis.domain = [0, 0.95]
    figure.layout.scene.domain.y = [0, 0.95]
    figure.axes_ranges((0, 1), (0, 1), scale=(1,1))
    figure.axes_labels("Beak length", "Beak pointiness")
    figure.axes_ranges((0, 1), (0, 1), (0, 1), scale=(1,1,1))
    figure.axes_labels("Beak length", "Beak pointiness", "Fitness")
    figure.layout.showlegend = False
    figure.add(plotly_text("Genotype space", (0.22,0.95), font_size=18, 
            paper=True, xanchor="center"))
    title = figure.add(plotly_text("Fitness landscape", (0.78,0.95), font_size=18, 
            paper=True, xanchor="center"))
    for (x0, y0), source in finches.items():
        figure.add_layout_image(source=source, x=x0, y=y0, sizex=0.2, sizey=0.2, 
                xref="x", yref="y", xanchor="center", yanchor="middle")
        figure.add_layout_image(source=source, x=0.5, y=0.5, sizex=0.5, sizey=1, 
                xref="paper", yref="paper", xanchor="left", yanchor="middle")
    fitness_landscape = figure.add(plotly_function3d(fitness, (x, 0, 1), (y, 0, 1), opacity=0.8))
    maxima = [(x0, y0, fitness(x0, y0)) for x0, y0 in finches]
    maxima = figure.add(plotly_points3d(maxima, size=1.5, color="darkgreen"))

    @interact(finch=slider(["None", 1, 2, 3, 4, "All"], label="Finches"), 
              show_graph=checkbox(False, label="Show fitness landscape"), 
              show_maxima=checkbox(False, label="Show species"))
    def update(finch, show_graph, show_maxima):
        finchnum = -1 if finch == "None" else 99 if finch == "All" else finch - 1
        with figure.batch_update():
            for i in range(4):
                figure.layout.images[2*i].visible = (i < finchnum)
                figure.layout.images[2*i + 1].visible = (i == finchnum)
            title.visible = show_graph
            fitness_landscape.visible = show_graph
            maxima.visible = show_maxima

    return figure

finches_interactive()

interactive(children=(SelectionSlider(description='Finches', options=('None', 1, 2, 3, 4, 'All'), value='None'…

FigureWidget({
    'data': [{'colorscale': [[0.0, 'lightblue'], [1.0, 'lightblue']],
              'opacity': …

Local maxima and local minima

Definition:

• A function has a local maximum at a point $x$ if there is some region around $x$ such that, within that region, the maximum value of the function occurs at $x$.

• A function has a local minimum at a point $x$ if there is some region around $x$ such that, within that region, the minimum value of the function occurs at $x$.

Also, maxima and minima collectively are called extrema, as in the extreme values of a function.

f(x) = (16*(x-3)^4 - 10*(x-3)^2 + 3*(x-3) + 4) * exp(-(x-3)^2)
df = f.derivative()
d2f = df.derivative()
critpts = find_zeros1d(df, (x, 0, 6))
colors = []
for x0 in critpts:
    evalue = d2f(x0)
    colors.append("darkgreen" if evalue < 0 else "darkred")
critpts = [(x0, f(x0)) for x0 in critpts]

figure = Figure()
figure.axes_ranges((0, 6), (0, 8))
figure.axes_labels("$x$", "$f(x)$")
figure.layout.margin.b = 80
figure.layout.margin.r = 200
figure.add(plotly_function(f, (x, 0, 6), name="Graph of function", plotpoints=161))
figure.add(plotly_points(critpts, color=colors, name="Extrema", visible="legendonly"))
figure

Conclusions:

1. Many many features (anatomical, behavioral, even biochemical) in biology seem to be the result of an optimization process. Any time one can say a species is “adapted to do ... very well”, that's probably describing such a result.

2. The biological mechanism underlying all of these optimizations is evolution, or more specifically, natural selection: survival of the fittest. The fitness of a certain genotype (or phenotype) is defined roughly as the average number of offspring an individual with that genotype will produce.

3. A function has a local maximum at a point $x$ if there is some region around $x$ such that, within that region, the maximum value of the function occurs at $x$. A local minimum is defined similarly.

4. One mathematical model for speciation is that species can form at any one of the local maxima of the fitness function, or fitness landscape.