Vector Constraints with Shared Vertices

[nickmccleery]

Hey @adamchalmers, dumping content here for discussion as I think this is quite insightful.

I think your formula for the derivative of the perpendicular constraint (see here)) only works if the two lines are between 4 separate points, and doesn't work if two of the lines share a common point.

Residual: R = (x1-x0)*(x3-x2) + (y1-y0)*(y3-y2)

If you have this situation where the lines L0 = (p0, p1) and L2 = (p0, p2) then the partial derivatives don't quite work out anymore, e.g. partial derivative of x0 is -x1-x3+2 or something.

Here's the general residual R of the perpendicular constraint: Then if we assume the two lines share a common point such that instead of p0, p1, p2 and p3 we have p0, p1 and p2, we get

[nickmccleery]

So, when I first thought about this yesterday, I was in agreement.

In the case that we go from four points to three points, our residual would absolutely change. We'd move from Residual: R = (x2-x1)*(x4-x3) + (y2-y1)*(y4-y3) to say R = (x2 - x1)(x3 - x1) + (y2 - y1)(y3 - y1), assuming P1 is the shared point. This would also come with our Jacobian moving from a set of 8 partial derivatives to a set of 6.

This scared me on the grounds that one possible path through the code I wrote was to 'rewrite' constraints where we had coincident constraints, for example, eliminating one of the two coincident points from our consideration.

However, having gone through it today, I now think:

• A) We don't need to implement any special handling of the common vertex case, as in your example - with one important proviso.
• B) I would scrap the constraint 'rewriting' and simply do variable level substitution and subsequent fan-out/distribution of solved values.

Example

Ignoring B, let me walk through the common vertex example. So, in my case, we're going to imagine your example just flipped horizontally. We have three points:

• $P_1$,
• $P_2$,
• $P_3$

Then we have two lines assembled from those:

• $L1 = \vec{P_1, P_2}$
• $L2 = \vec{P_2, P_3}$

Then we set up the system like this:

def constrain_perpendicular_with_shared_vertex():
    """
    Build a simple polyline p1-p2-p3 with two segments:
    L1 = (p1, p2), L2 = (p2, p3)

    Then constrain L1 perp L2 at the shared vertex p2.
    """
    p1 = Point(x=0.0, y=0.0, id="P1")
    p2 = Point(x=3.2, y=0.7, id="P2")
    p3 = Point(x=4.8, y=2.1, id="P3")

    l1 = Line(p1, p2, id="L1")
    l2 = Line(p2, p3, id="L2")

    # Ensure segments share the middle vertex
    assert l1.p2 is l2.p1, "Segments must share the middle vertex"

    primitives = [p1, p2, p3, l1, l2]

    constraints = [
        PointFixed(point=p1),
        LineHorizontal(l1),
        PointPointXDistance(p1, p2, distance=4.0),
        PointPointYDistance(p2, p3, distance=3.0),
        LinesPerpendicular(l1, l2),
    ]

In my case, the relevant constraint looks like this:

@dataclass
class LinesPerpendicular(BaseConstraint):
    line1: Line
    line2: Line

    def get_residual(self, variable_values: Mapping[str, float]) -> nb.Vector:
        p1 = self.line1.p1.get_state(variable_values)
        p2 = self.line1.p2.get_state(variable_values)

        p3 = self.line2.p1.get_state(variable_values)
        p4 = self.line2.p2.get_state(variable_values)

        # Calculate the vectors.
        v1 = p2 - p1
        v2 = p4 - p3

        # The residual is the dot product of the two vectors.
        # If they are perpendicular, this should be zero.
        dot = v1[0] * v2[0] + v1[1] * v2[1]
        return nb.np.array([dot])

    def get_jacobian_row_values(
        self, variable_values: Mapping[str, float]
    ) -> List[Tuple[str, float, int]]:
        # Residual: R = (x2-x1)*(x4-x3) + (y2-y1)*(y4-y3)
        # ∂R/∂x1 = x3 - x4
        # ∂R/∂y1 = y3 - y4
        # ∂R/∂x2 = -x3 + x4
        # ∂R/∂y2 = -y3 + y4
        # ∂R/∂x3 = x1 - x2
        # ∂R/∂y3 = y1 - y2
        # ∂R/∂x4 = -x1 + x2
        # ∂R/∂y4 = -y1 + y2

        # Get points.
        p1 = self.line1.p1.get_state(variable_values)
        p2 = self.line1.p2.get_state(variable_values)
        p3 = self.line2.p1.get_state(variable_values)
        p4 = self.line2.p2.get_state(variable_values)

        # Get their components.
        x1, y1 = p1[0], p1[1]
        x2, y2 = p2[0], p2[1]
        x3, y3 = p3[0], p3[1]
        x4, y4 = p4[0], p4[1]

        # Calculate derivatives.
        dr_dx1 = x3 - x4
        dr_dy1 = y3 - y4
        dr_dx2 = -x3 + x4
        dr_dy2 = -y3 + y4
        dr_dx3 = x1 - x2
        dr_dy3 = y1 - y2
        dr_dx4 = -x1 + x2
        dr_dy4 = -y1 + y2

        # Make floats.
        dr_dx1 = float(dr_dx1)
        dr_dy1 = float(dr_dy1)
        dr_dx2 = float(dr_dx2)
        dr_dy2 = float(dr_dy2)
        dr_dx3 = float(dr_dx3)
        dr_dy3 = float(dr_dy3)
        dr_dx4 = float(dr_dx4)
        dr_dy4 = float(dr_dy4)

        # This constraint has a scalar residual.
        i_residual = 0

        # Get the variable IDs for the points involved.
        p1_x_var = self.line1.p1.get_variable_ids()[0]
        p1_y_var = self.line1.p1.get_variable_ids()[1]
        p2_x_var = self.line1.p2.get_variable_ids()[0]
        p2_y_var = self.line1.p2.get_variable_ids()[1]
        p3_x_var = self.line2.p1.get_variable_ids()[0]
        p3_y_var = self.line2.p1.get_variable_ids()[1]
        p4_x_var = self.line2.p2.get_variable_ids()[0]
        p4_y_var = self.line2.p2.get_variable_ids()[1]

        return [
            (p1_x_var, dr_dx1, i_residual),
            (p1_y_var, dr_dy1, i_residual),
            (p2_x_var, dr_dx2, i_residual),
            (p2_y_var, dr_dy2, i_residual),
            (p3_x_var, dr_dx3, i_residual),
            (p3_y_var, dr_dy3, i_residual),
            (p4_x_var, dr_dx4, i_residual),
            (p4_y_var, dr_dy4, i_residual),
        ]

Now, when the solver runs and we come through the residual function, line1.p2 is equal to line2.p1:

But if you look at our actual residual function, there's nothing about our dot product that is really affected by this setup compared to a four point equivalent. Under the hood, we're just pulling X and Y coordinates and doing R = (x2-x1)*(x4-x3) + (y2-y1)*(y4-y3); it doesn't affect anything that we're using a shared point.

Then, as we come through the derivative calc building the partials to slot into the Jacobian, we pull out the eight values we need, assuming four points: x1 - x4 and y1 - y4.

Again, nothing's really affected at this point, but we can see that the we're pulling through the P2 point twice. What is important to note here is anything that involves these terms, i.e., anything that involves x2 or x3 or y2 or y3. In reality, all of our partials reference these.

The important bit is that we flag that these refer to the same variable in our output:

In effect, what we're saying is that dr_dx2, dr_dy2, dr_dx3, and dr_dy3 aren't really unique/independent partial derivatives...

• dr_dx2 tells us about our residual's sensitivity to a change in P2_x
• dr_dx3 also tells us about our residual's sensitivity to a change in P2_x
• dr_dy2 tells us about our residual's sensitivity to a change in P2_y
• dr_dy3 also tells us about our residual's sensitivity to a change in P2_y.

Basically, each of these values tells us part of the residual's sensitivity to that variable; all we have to do is sum them to find the overall $\frac{\partial R}{\partial P2_x}$ and $\frac{\partial R}{\partial P2_y}$ terms.

If you look at build_sparse_jacobian as an example, we actually already take the response of get_jacobian_row_values, then iterate across it summing any partials which refer to the same variable:

...
                entries = constraint.get_jacobian_row_values(variable_values)

                for var_id, value, i_row_local in entries:
                    root_var = find(var_id, substitution_map)
                    if root_var in var_map:
                        i_col = var_map[root_var]
                        jacobian[i_row + i_row_local, i_col] += value # <----------------- here
...

In the concrete case of the example above, you can see that dr_dx2 was 1.5999999999999996. As we come through our Jacobian builder, where the overall $\frac{\partial R}{\partial P2_x}$ element already has a nonzero value, we've set it to dr_dx2, and we're about to update it by summing it with the value of dr_dx3, which is -3.2.

So that's basically it; the proviso is that we can continue with our existing approach so long as we recognise that the terms falling out of the get_jacobian_row_values may actually be only part of the partial derivative term for a given variable, and our overall Jacobian assembly routine needs to be able to identify this and accumulate terms for variables which appear multiple times.

TL;DR

We can just stick with the existing method, referencing a common vertex more than once, provided we recognise that that the partial derivatives returned by get_jacobian_row_values are not the whole story. Where more than one of these terms relates to a single variable in our system, e.g. P2_x, we need to identify that and sum the constituents together.

Anecdotally, have tested this with a three point shared vertex example and test and seems to work as expected.

[nickmccleery]

I punted my waffle at resident numerical methods oracle @davreev and he has some nice terminology recs here: local vs. global Jacobian.

It's quite a nice way to think about it: inside our get_jacobian_row_values methods we can treat $P1 - P4$ as all being unique things and our eight variables being locally unique.

However, jump up one level to our global Jacobian assembly method in build_sparse_jacobian, and there we can actually see that there are some aliases involved; two of our points which we thought were unique actually point at the same thing, and we collapse from eight local variables to six global ones.

[nickmccleery]

Python PR with some relevant changes and a test case: KittyCAD/geometric-constraint-playground#8

[nickmccleery]

I actually think this should be fine now since I snuck in partial accumulation under #80

Maybe we can add a test case?