2. Modeling¶

ConstraintHg is a model organizer. That means it can take a given model and structure it so that other models can connect to it without information or meaning being lost. However, you still have to do all the work of providing the model. It is often more convenient to do this work in a different framework. The pendulum example, for instance, can be modeled using Newtonian physics, or setting up the standard Euler-Langrange equations.

Note

If you want to see a discussion of modeling the pendulum dynamics, check out the overview here.

In the linked example, we derive the following constraint hypergraph, which we’ll now setup in the Python script.

CHG of a pendulum — *Constraint hypergraph of the pendulum.*¶

The first thing to do is to initialize the hypergraph:

hg = Hypergraph()

There are a few options you can call when initializing this, such as setting up logging. More information is available at the API. Now we can start passing in the variables.

2.1. Nodes¶

The constraint hypergraph has several variables which we’ll want to include in our model. Think of a the set of variables as all the information we might want to know about the pendulum. Each variable becomes a node in the graph.

g = hg.add_node(Node('gravity', -10))
r = hg.add_node(Node('radius', 0.5))
theta0 = hg.add_node(Node('theta0', 3.14159/6))
theta = hg.add_node(Node('theta'))
d_theta = hg.add_node(Node('delta theta'))
s_theta = hg.add_node(Node('sine theta'))
F = hg.add_node(Node('gravitational force'))
omega0 = hg.add_node(Node('omega0', 0.0))
omega = hg.add_node(Node('omega'))
d_omega = hg.add_node(Node('delta omega'))
c = hg.add_node(Node('damping coeff', 1.5))
alpha = hg.add_node(Node('alpha'))
time_step = hg.add_node(Node('time_step', .03))
time = Node('time', 0.)

The hg.add_node() method returns a Node object, which we can use later to set up the edges in the graph.

In the above we also defined some intermediate notes, these are variables that didn’t show up in the original equations, but which will show up in the graph. This is because algebra often condenses multiple operations into a single equation, though for the hypergraph its better to break these operations out so that we can interface with the hidden variables.

2.2. Edges¶

Each edge in a constraint hypergraph maps the values of its source set to the value of its target node. Mathematically this is just function mapping, though you can think of it as assigning a value of a node to every ordered pair of another set of variables. For example, given nodes \(A, B := \lbrace 1, 2 \rbrace\), and \(C := \lbrace -1, 0, 1 \rbrace\), we could define a hyperedge such that:

A	B	C
1	1	0
1	2	-1
2	1	1
2	2	0

Note that every possible pair \(A\) and \(B\) need to be mapped to a corresponding value in \(C\). We’ll see later on how to take a subset of values in the relation only. Such a method could encode this table, or we could just note that \(C = A - B\). To pass this relationship into the hypergraph, we would just need to pass along a method that took in a value of \(A\), a value of \(B\), and returned the correct value of \(C\). For the example above, we might write a method like the following:

def subtract_A_from_B(A, B):
    return A - B

Note that although we’re writing these as methods in Python, they need to behave as functions in the strictest sense. For instance, we cannot modify any variables from the system inside the method–only the variable listed as the function’s target (or codomain).

Note

Previous versions of ConstraintHg required the variable specifiers *args and **kwargs to be given as arguments to any relation. This was updated in version 0.3.0 to no longer be needed.

To simplify coding, many relationships have been specified in the relations module. These were imported into the script with import constrainthg.relations as R, and then called as R.Rsubtract for example. We’ll use mostly these provided relations in the demo, but we will need one custom function for the Eulerian integration. Let’s define it now so we can pass it later on:

def Rintegrate(step, slope, initial_val):
    """First order Eulerian integrator."""
    return step * slope + initial_val

Returning to the pendulum, let’s put in some of the more simple edges. The main method here is Hypergraph.add_edge. The required arguments for `add_edge` are the source nodes (sources), the target node (target), and the relation (rel). Nodes, whether sources or targets, can be passed as a string or a full Node object. For the add_edge method, the source nodes can be a single node, a list, or a dictionary of key : string/node pairs (this allows you to reference the nodes by the passed keywords in the constraint method). The relation is a method used for calculating the target node value. Another argument, label, is an optional string ID that helps us uniquely identify the edge.

hg.add_edge(
    sources=theta0,
    target=theta,
    rel=R.Rmean,
    label='theta0->theta',
)
hg.add_edge(
    sources=omega0,
    target=omega,
    rel=R.Rmean,
    label='omega0->omega',
)
hg.add_edge(
    sources={'s1': g, 's2': r}, #dictionary of source nodes ensures 'g'is the numerator,
    target='g/r',
    rel=R.Rdivide,
    label='(g,r)->b1',
)
hg.add_edge(
    sources=theta,
    target=s_theta,
    rel=R.Rsin,
    label='theta->sine',
)
hg.add_edge(
    sources={'s_theta': s_theta, 'g/r': 'g/r'},
    target=F,
    rel=R.Rmultiply,
    label='(sine, b1)->F',
)
hg.add_edge(
    sources={'omega': omega, 'c': c},
    target='beta2',
    rel=R.Rmultiply,
    label='(omega, c)->b2',
)
hg.add_edge(
    sources={'F': F},
    target=alpha,
    rel=R.equal('F'),
    label='F->alpha',
)

2.3. Printing¶

Now that we have our basic hypergraph, it’s a good idea to query its structure to make sure we put everything in the right place. The package has some basic functionality for printing the paths in a Hypergraph to the terminal.

Note

A path in a hypergraph is a chain of edges through the hypergraph. In a hypergraph, this chain appears as a tree.

Let’s start by printing all the possible paths for calculating the angular acceleration (\(\alpha\)) of the pendulum. To do this, add the following line to the script and run the program:

print(hg.summary(alpha))

You should get the following output in your terminal:

└──alpha, index=1, cost=5
    └──gravitational force, index=1, cost=4
        ├◯─sine theta, index=1, cost=2
        │  └──theta, index=1, cost=1
        │     └──theta0, index=1
        └●─g/r, index=1, cost=1
            ├◯─radius, index=1
            └●─gravity, index=1

Each section of the output corresponds to an node in a path that eventually leads to alpha. The source nodes for the edges leading to the higher are indented on the next lines. If there is one source node, then it’s given with a single angle pipe └──. Otherwise, for multiple sources, the nodes are specified with circles ├◯─. The last source node for an edge is distinguished with a filled circle └●─.

In our output, you can see that there is only one path for calculating alpha, with five edges in the path.

Note

The cost of a node is taken by summing the weights of all edges in the path leading up to it. The method prints the minimum summation. Since we used the default edge weight of 1, here the cost indicates how many edges form the shortest path to the given node.

We can start simulating any path in the graph immediately, but to get full functionality we also need to learn about cycles, as well as figure out what the index parameter means. Click here to go to the next step or use the navigation below.