toolhead: Fixed junction deviation calculation for straight segments

[dmbutyugin]

As per suggestion in #6723, opening a separate PR for the toolhead junction deviation fix. As for

I'm concerned in may trade one unusual "corner case" for another one. Specifically, it'll close a loop-hole where we may pessimize speeds if acceleration changes midway through a straight line movement, but I fear it opens a loop-hole where a slicer could emit an arc with a series of infinitesimal moves such that each of those moves appears to be "close to straight" which then results in the arc proceeding with no slow-down at all.

I believe right now the constant +/-0.999999 is somewhat arbitrary, and it is easy to add one-two more 9s to it. And for example currently we have tan_theta_d2 = 1414.21.. for junction_cos_theta = -0.999999 and .5 * tan_theta_d2 * self.accel = 707106.6... (with accel = 1000), and for a move of 0.01 mm it would give the cornering velocity of 84 mm/sec (from centripetal acceleration). Adding one more 9 to the constant would give ~150 mm/sec, and one more - 266 mm/sec, if I did not mess up the calculations. So, if we run into issues with slicers, we could adjust the thresholds to +/-0.99999999 and that would likely address them.

[dmbutyugin]

FWIW, I'm unsure about the checks failure, but I suspect this is unrelated to my change. scripts/check_whitespace.sh gives no errors, and running ~/klippy-env/bin/python scripts/test_klippy.py -d dict/ test/klippy/*.test with the latest available klipper-dict-20240313.tar.gz gives me an error msgproto.error: mcu 'mcu': Unknown command: query_ads1220 on the test test/klippy/load_cell.test.

[MRX8024]

Hi Dmitry, I also noticed this problem once (with stops), but it was not fundamental from the point of view of normal use, and manifests itself only at very low acceleration. You probably fixed this for your new resonance measurement algorithm, for which I also want to thank you.

I didn't really understand what cases Kevin was talking about with "close to straight" moves. I've plotted the old and new code, and I haven't seen any changes. I wanted to open a topic on discourse, but I caught sight of this PR today, and I outlined the part of the material here.

For example, this is how Klipper slices the movements generated by gcode from SuperSlicer, at maximum resolution, without using arcs. On the left is the old, on the right is the new one, the toolhead parameters (speed, etc) are shown below the graph -

I am also stupefied by the current way of determining the speed on circles. The worse your gcode file is, the larger length of polygons it have, the slower the circles\turns that are so fond of adding to all 3d models now will be printed. But I do not know the right alternatives, unfortunately.

I include arcs to simulate more detail of polygons.
Fact: current slicers are very bad at generating arcs from stl files, and they can't do it directly from step. For example, they can place at least 20 small arches on one circle, instead of one, which will not connect well with each other, and cause toolhead slows.

The same code, but with arcs enabled, the resolution of arcs in the slicer is 0.01, he divided this circle into two arcs. The resolution of arcs for the Klipper ArcParser is also signed at the bottom of the graph. By default, the arcs in the Klipper have a section resolution of 1 mm, which simply does not make sense in conventional printing. Let's start with him -

To put it mildly, it does not look like a circle, and the average speed has noticeably decreased, due to the fact that the angle between the sections is large.

Next, the arc resolution is 0.1 mm, this is the most popular value among users-

This already looks like a circle), please also note that the average printing speed of the circle has risen significantly. But, again, there are no changes regarding the previous code. And you can still see the curved line after the arc, that is, the arc ended in the wrong place for some reason.

And also let's see the resolution of 0.01mm, for example, rpi5 can safely withstand such large-scale calculations during printing.

I noticed this bug a couple of months ago, the more accurate the resolution (< 0.1mm), the more the speed decreases, although the angle between the sections should be closer to straight. So far, I have not figured out what is the matter, perhaps the ArcParser somehow crookedly convert the arcs..

I will also attach the GCode files that I used -
Body1_.txt
Body1_arc.txt

In conclusion, I have not found any practical places where Dmitry's new code would change anything in a negative way. But I think that the klipper should not slow down when printing, for example, a 100mm cylinder, just because it was not given a sufficiently accurate GCode. I will continue to experiment and look for solutions to this problem. And I'll probably duplicate it into a Klipper discourse soon, too. Correct me if I don't understand something. Thanks!

-Maxim

[KevinOConnor]

I believe right now the constant +/-0.999999 is somewhat arbitrary

Yes - it's only there to prevent a divide by zero. So, what if we did this instead:

--- a/klippy/toolhead.py
+++ b/klippy/toolhead.py
@@ -72,7 +72,7 @@ class Move:
                                + axes_r[2] * prev_axes_r[2])
         if junction_cos_theta > 0.999999:
             return
-        junction_cos_theta = max(junction_cos_theta, -0.999999)
+        junction_cos_theta = max(junction_cos_theta, -0.999999999999)
         sin_theta_d2 = math.sqrt(0.5*(1.0-junction_cos_theta))
         R_jd = sin_theta_d2 / (1. - sin_theta_d2)
         # Approximated circle must contact moves no further away than mid-move

FWIW, I'm unsure about the checks failure, but I suspect this is unrelated to my change

Yeah, I'm not sure either. You could try rebasing this to the latest Klipper master so that it picks up commit f2e69a3.

I didn't really understand what cases Kevin was talking about with "close to straight" moves.

Consider the case where the toolhead is at position 0,0 and the following two moves are made: G1 X100 Y0 and G1 X0 Y0. In this case, the toolhead makes a 180 degree turn at position x=100 and we'd expect the toolhead to perform a zero velocity junction speed between the two moves. Now consider these three moves: G1 X100 Y0, G1 X100 Y0.000000000001, G1 X0, Y0. In this case, we make a 90 degree turn, followed by a 90 degree turn, followed by a 90 degree turn. The original junction deviation code (from GRBL) would perform these moves with a junction speed at square_corner_velocity (typically 5mm/s). This is wrong because the move is effectively identical to the original move above - the addition of an infinitesimal move shouldn't have allowed a faster junction speed. Now consider the following moves: G1 X100 Y0, G1 X100.000000000001 Y0.000000000001, G1 X0, Y0 - which is three 60 degree angles, which the original GRBL junction deviation code would take at much higher than 5mm/s. Keep adding infinitesimal moves and one can turn 180 degrees at full speed - causing extreme toolhead "jerkiness". That is, keep adding infinitesimal moves and one can create a series of arbitrary small angles that ultimately add up to a full 180 degree turn with no slowdown. Because these infinitesimal moves are well below a step size, they don't alter the actual final stepper move timings - they just alter the junction speeds.

To address this, the code in commit a4439b9 was added so that the code does not calculate the junction speed solely on angles - it also looks at the move length. This avoids the problem of these infinitesimal moves that some users were running into when printing some models (and/or with some slicer settings).

The PR here would bypass the code in a4439b9 based on the angle of moves, but I'm concerned that may reintroduce a corner case because one can introduce a series of infinitesimal moves that all have an arbitrarily small intra-move angle.

Cheers,
-Kevin

[dmbutyugin]

OK, thank you for the feedback, Kevin. Let me think about your examples a bit (and your alternative proposal too, which I'd need to give a test).

[Wulfsta]

@KevinOConnor, if R_jd is just used in a min statement where it is multiplied by other floats, why not just use the following rather than this arbitrary maximum?

R_jd = sin_theta_d2 / (1. - sin_theta_d2) if sin_theta_d2 != 1.0 else float('Inf')

[dmbutyugin]

@KevinOConnor

To address this, the code in commit a4439b9 was added so that the code does not calculate the junction speed solely on angles - it also looks at the move length. This avoids the problem of these infinitesimal moves that some users were running into when printing some models (and/or with some slicer settings).

FWIW, I think the current code still does it a bit backwards. That is, it has somewhat arbitrary limit for junction_cos_theta, and depending on the actual segmentation resolution for an arc and the length of the segments, the users will get different centripetal acceleration limits. My proposal to address this is in b97d6e2 commit, where I propose to check whether the radius of curvature is finite or not. This way the users will get the consistent results more or less independent on the resolution of segmentation (up to rounding errors in GCode coordinates). While I'd need a bit more testing of it overall (e.g. with my new resonance test), I'd like to get your early feedback on this approach.

[KevinOConnor]

if R_jd is just used in a min statement where it is multiplied by other floats, why not just use the following rather than this arbitrary maximum?

Well, that's basically what the code attempts to do today. Using 'Inf' is not desirable because those special floating point values can propagate in subtle ways and lead to hard to debug errors in other parts of the code. So, the goal is to use an arbitrarily large constant instead of 'Inf' . The -0.999999 was intended to bring about that arbitrarily large constant - if the resulting value isn't big enough then we can add more 9s (eg, -0.999999999999).

My proposal to address this is in b97d6e2 commit, where I propose to check whether the radius of curvature is finite or not.

I don't understand what that additional check does (max(self.move_d, prev_move.move_d) < 99999999.9 * cos_theta_d2). At first glance though, it seems like that code would "round up" to higher junction speeds where I think we really want to "round down" in unusual corner cases like this.

At a minimum, the checks on junction_cos_theta should occur prior to running math.sqrt(), as we don't want any weird rounding issues to cause a math fault (eg, sqrt on infinitesimally small negative number).

FWIW, if this change is not a prerequisite for #6723 then I think we should merge that first.

Cheers,
-Kevin

[dmbutyugin]

Thanks Kevin,

I don't understand what that additional check does (max(self.move_d, prev_move.move_d) < 99999999.9 * cos_theta_d2). At first glance though, it seems like that code would "round up" to higher junction speeds where I think we really want to "round down" in unusual corner cases like this.

The curvature radius is equal to 0.5 * move_d * tan(theta/2) = 0.5 * move_d * sin(theta/2) / cos(theta/2). What I want to judge is whether this radius evaluates to some finite (even if relatively large) number or to an infinite number (or a very large number due to a loss of precision). The latter can happen only if theta -> pi. In this case sin(theta/2) -> 1. So in principle we want to check if move_d / cos(theta/2) is finite or not. So I basically check if move_d / cos(theta/2) < 99999999.9 or not. But in order to avoid issues with division, after rearranging the inequality, I check if move_d < 99999999.9 * cos(theta/2). Note that this condition works as expected even if cos(theta/2) == 0.0 precisely (because for any move_d >= 0, move_d < 0.0 is false). And since we have two moves, I take actually a min of the two lengths (notably, it was max previously, which is in general not correct).

Edit: and just to expand on this point a bit, if we determine that the radius of curvature is infinite (or very large), we should not add additional limits on junction_v2 that is related to centripetal acceleration calculation, which is the behavior of the new code in this PR. FWIW, the regular junction deviation should also not apply limits in this case, but since that code does not depend on move_d, it is sufficient to use the regular junction_cos_theta >= -0.999999 comparison (basically, if junction_cos_theta ~= -0.999999, R_jd is already very large, which ensures the sufficient continuity of code behavior between junction_cos_theta ~= -0.999999 but above it and junction_cos_theta < -0.999999).

At a minimum, the checks on junction_cos_theta should occur prior to running math.sqrt(), as we don't want any weird rounding issues to cause a math fault (eg, sqrt on infinitesimally small negative number).

Yes, you are right, I messed that up previously. I addressed that now, though a bit differently (by taking max(1.0+junction_cos_theta, 0.) prior to computing sqrt(), and note that sin_theta_d2 is free of this issue due to the check above). I do not check junction_cos_theta value directly because the correct way is to estimate the relation between move_d and cos(theta/2), and doing that before taking a square root would just complicate the inequality and add unnecessary calculations.

FWIW, if this change is not a prerequisite for #6723 then I think we should merge that first.

No, this change is a prerequisite for that one, so this PR should go in first.

[KevinOConnor]

No, this change is a prerequisite for that one, so this PR should go in first.

Oh - okay. I misunderstood that.

The curvature radius is equal to

Okay, thanks. I think I understand now.

However, I'm not sure about this approach for two reasons: it seems to be "rounding up" to higher junction speeds where I think we should prefer "rounding down", and this adds additional complexity to the "fast path".

On the first point, it seems this code would basically say "if it seems the junction is mostly straight then assume it is straight and don't limit speed at all". But that's seems backwards to me - if we're unsure about the math we should round down the junction speeds, not round them up. When in doubt, we want to go slower through junctions, not present a possibility that we go too fast. These are all "corner cases" that should not impact real-world prints, so better to "do no harm".

We could change the first test to if 1. - sin_theta_d2 >= 0.: and change the second test to if cos_theta_d2 >= 0.:. That would more directly test for a divide by zero, and not arbitrarily "round up" junction speeds. However, we'd still have the higher code complexity on the "fast path". This code adds two additional if branches, three additional min() calls, and two additional assignments to self.max_start_v2. I'm sure that could be optimized, but I'm not sure why new code is preferred over keeping the existing code with a different range constant (ie, junction_cos_theta = max(junction_cos_theta, -0.999999999999) . Ultimately, we should get the same results for all real-world use cases, right?

Thoughts?
-Kevin

EDIT: Just for perspective, if we limit junction_cos_theta to -0.999999999999 then R_jd would be 3999644429280 (or >=6435656mm/s with square_corner_velocity=5 and accel>=100) and tan_theta_d2 would be 1414229 (or >=8409mm/s with move_d>=1 and accel>=100).

[KevinOConnor]

No, this change is a prerequisite for that one, so this PR should go in first.

On further thought, I'm a little confused. What is the new code doing that requires this change?

Based on my quick calculations, the current code limits junction_cos_theta to -0.999999, which would result in an R_jd=3999998 (>=6435mm/s with typical values) and a tan_theta_d2=1414 (>=265mm/s with move_d>=1 and accel>=100).

I think I may be missing something subtle. What is the test case doing that would be limited by the above?

Thanks,
-Kevin

[dmbutyugin]

@KevinOConnor

it seems to be "rounding up" to higher junction speeds where I think we should prefer "rounding down", and this adds additional complexity to the "fast path".

Sorry, I'm unsure I understand where 'rounding up' and 'rounding down' comes from. FWIW,

it seems this code would basically say "if it seems the junction is mostly straight then assume it is straight and don't limit speed at all".

is mostly right, however not entirely. The problem at hand is that as the segmentation for a curve becomes finer and finer, the angle between consecutive segments becomes closer and closer to 180 degrees (and cos(theta/2) approaches 0). And so, the current mainline Klipper code has an issue: it has a minimum limit for cos(theta/2), and thus technically produces incorrect results for angles theta very close to 180 degrees, both for consecutive straight segments and for very small segments forming a curve. Adjusting the limit (e.g. setting it to -0.999999999999) just pushes the problematic cases further down the line. My approach is based on the following observation: for smooth curves, as the segmentation becomes finer, both move_d and cos(theta/2) go to 0, but their ratio does not (it does not go to either 0 or infinity), as limit(move_d / cos(theta/2)) = 2 * R. Thus, by checking the ratio move_d / cos(theta/2), we can tell apart the cases "a curve with finite curvature, just very fine"and "(sufficiently) straight segments" regardless of how fine the segmentation is, and without having to play with limits too much, and apply the junction velocity limits from centripetal acceleration only to the first case. You could still say that the code contains a threshold, but the behavior of the code does not depend on its value too much, basically if the threshold is sufficiently large, it won't matter how large exactly it is. This is not the case for junction_cos_theta, as for any threshold you could choose sufficiently small segmentation that would make code produce technically incorrect junction velocity limits.

However, we'd still have the higher code complexity on the "fast path"... I'm sure that could be optimized, but I'm not sure why new code is preferred over keeping the existing code with a different range constant (ie, junction_cos_theta = max(junction_cos_theta, -0.999999999999) . Ultimately, we should get the same results for all real-world use cases, right?

I tried to explain above why I believe the approach with adjusting the threshold is not optimal: it may still cause issues for sufficiently fine gcode. And I suspect the effect of the changes on the code performance is almost negligible: there are far more complex calculations running per move, e.g. several sqrt operations, then look-ahead processing with moving the moves in and out of the queue, and such. I also did not want to invest into premature optimizations, however some optimizations can be done (e.g. removing if junction_cos_theta >= -0.999999: conditional and putting the corresponding code under the other if, using local max_start_v2 for intermediate calculations, maybe some additional math rearrangements). However, if you insist on not using this approach, I can see if just going with a different threshold proposed by you would resolve the problems that the new resonance test runs into.

No, this change is a prerequisite for that one, so this PR should go in first.

On further thought, I'm a little confused. What is the new code doing that requires this change?

Based on my quick calculations, the current code limits junction_cos_theta to -0.999999, which would result in an R_jd=3999998 (>=6435mm/s with typical values) and a tan_theta_d2=1414 (>=265mm/s with move_d>=1 and accel>=100).

I think I may be missing something subtle. What is the test case doing that would be limited by the above?

As was reported in the other PR, the problematic code snippet (presented in this comment) is very similar to what the new test executes from time to time, leading to the wrong velocities and timing of moves.

[KevinOConnor]

Sorry, I'm unsure I understand where 'rounding up' and 'rounding down' comes from. FWIW,

What I meant is, if there is a move that results in a natural calculation of junction_cos_theta=-0.999999 (eg, G1 X0 Y0 Z0, G1 X10 Y0, G1 X110 Y0.1) then the new code would not limit the junction speed. That seems odd to me - I understand that if we truly have a straight move (eg, G1 X10 Y0, G1 X110 Y0) then we don't want to impose a limit, but it seems strange to artificially go faster through junctions that aren't actually straight.

So, by "round up" I meant that the new code identifies a set of unusual corner cases and allows them to go faster than they nominally might, where as the current code can make some moves go slower but should not permit moves to go faster than they nominally should ("round down").

it may still cause issues for sufficiently fine gcode. And I suspect the effect of the changes on the code performance is almost negligible

I agree that the performance impact is negligible and I'm okay with changing the code if needed. If we do need to change the code, could we directly test for the divide by zero case so we avoid optimizing moves that are not actually straight? (That is, if 1. - sin_theta_d2 >= 0.: and if cos_theta_d2 >= 0.:)

As was reported in the other PR, the problematic code snippet (presented in #6723 (comment)) is very similar to what the new test executes from time to time, leading to the wrong velocities and timing of moves.

So the new resonance testing code effectively issues a SET_VELOCITY_LIMIT ACCEL=0.1?

Thanks again,
-Kevin

[KevinOConnor]

Thinking about this further, I think the 0.999999 constants are arbitrary and kinda confusing. How about we just remove them with something like:

--- a/klippy/toolhead.py
+++ b/klippy/toolhead.py
@@ -65,25 +65,30 @@ class Move:
         # Allow extruder to calculate its maximum junction
         extruder_v2 = self.toolhead.extruder.calc_junction(prev_move, self)
         # Find max velocity using "approximated centripetal velocity"
+        move_jd_v2 = prev_move_jd_v2 = extruder_v2
         axes_r = self.axes_r
         prev_axes_r = prev_move.axes_r
         junction_cos_theta = -(axes_r[0] * prev_axes_r[0]
                                + axes_r[1] * prev_axes_r[1]
                                + axes_r[2] * prev_axes_r[2])
-        if junction_cos_theta > 0.999999:
-            return
-        junction_cos_theta = max(junction_cos_theta, -0.999999)
-        sin_theta_d2 = math.sqrt(0.5*(1.0-junction_cos_theta))
-        R_jd = sin_theta_d2 / (1. - sin_theta_d2)
+        sin_theta_d2 = math.sqrt(max(0.5*(1.0-junction_cos_theta), 0.))
+        one_minus_sin_theta_d2 = 1. - sin_theta_d2
+        if one_minus_sin_theta_d2 > 0.:
+            R_jd = sin_theta_d2 / one_minus_sin_theta_d2
+            move_jd_v2 = R_jd * self.junction_deviation * self.accel
+            prev_move_jd_v2 = (R_jd * prev_move.junction_deviation
+                               * prev_move.accel)
         # Approximated circle must contact moves no further away than mid-move
-        tan_theta_d2 = sin_theta_d2 / math.sqrt(0.5*(1.0+junction_cos_theta))
-        move_centripetal_v2 = .5 * self.move_d * tan_theta_d2 * self.accel
-        prev_move_centripetal_v2 = (.5 * prev_move.move_d * tan_theta_d2
-                                    * prev_move.accel)
+        move_centripetal_v2 = prev_move_centripetal_v2 = extruder_v2
+        cos_theta_d2 = math.sqrt(max(0.5*(1.0+junction_cos_theta), 0.))
+        if cos_theta_d2 > 0.:
+            tan_theta_d2 = sin_theta_d2 / cos_theta_d2
+            move_centripetal_v2 = .5 * self.move_d * tan_theta_d2 * self.accel
+            prev_move_centripetal_v2 = (.5 * prev_move.move_d * tan_theta_d2
+                                        * prev_move.accel)
         # Apply limits
         self.max_start_v2 = min(
-            R_jd * self.junction_deviation * self.accel,
-            R_jd * prev_move.junction_deviation * prev_move.accel,
+            move_jd_v2, prev_move_jd_v2,
             move_centripetal_v2, prev_move_centripetal_v2,
             extruder_v2, self.max_cruise_v2, prev_move.max_cruise_v2,
             prev_move.max_start_v2 + prev_move.delta_v2)

That should reduce the weird corner cases and is maybe a little easier to read.

Thoughts?
-Kevin

EDIT: Fixed broken patch.

[dmbutyugin]

Thanks Kevin,

What I meant is, if there is a move that results in a natural calculation of junction_cos_theta=-0.999999 (eg, G1 X0 Y0 Z0, G1 X10 Y0, G1 X110 Y0.1) then the new code would not limit the junction speed. That seems odd to me - I understand that if we truly have a straight move (eg, G1 X10 Y0, G1 X110 Y0) then we don't want to impose a limit, but it seems strange to artificially go faster through junctions that aren't actually straight.

I understand that you meant R_jd-related calculations here, right? Because for centripetal acceleration this is not correct:

junction_cos_theta = -100. / math.sqrt(100.**2 + 0.1**2) # == -0.9999995000003751
cos_theta_d2 = math.sqrt(0.5 * max(1.0+junction_cos_theta, 0.)) # == 0.0004999998124627294
min(self.move_d, prev_move.move_d) < 99999999.9 * cos_theta_d2 #  == True, because 10.0 < 49999.98119627296

Incidentally, even G1 X110 Y0.001 would work correctly for centripetal acceleration:

junction_cos_theta = -100. / math.sqrt(100.**2 + 0.001**2) # == -0.99999999995
cos_theta_d2 == 5.000000206850923e-06
99999999.9 * cos_theta_d2 == 500.0000201850923

Admittedly, G1 X110 Y0.00001 would no longer work, but we need to draw a line somewhere (Edit: it would no longer work in a sense that the two moves G1 X10 Y0, G1 X110 Y0.00001 will be considered collinear by the code and no junction limits related to centripetal acceleration applied). So, my current proposal would be along these lines (note that R_jd-related calculations are moved into the common if, this is because under that condition we can guarantee that 1 - sin_theta_d2 > 0.5 * (move_d / 99999999.9)**2):

         if junction_cos_theta > 0.999999:
             return
-        junction_cos_theta = max(junction_cos_theta, -0.999999)
-        sin_theta_d2 = math.sqrt(0.5*(1.0-junction_cos_theta))
-        R_jd = sin_theta_d2 / (1. - sin_theta_d2)
-        # Approximated circle must contact moves no further away than mid-move
-        tan_theta_d2 = sin_theta_d2 / math.sqrt(0.5*(1.0+junction_cos_theta))
-        move_centripetal_v2 = .5 * self.move_d * tan_theta_d2 * self.accel
-        prev_move_centripetal_v2 = (.5 * prev_move.move_d * tan_theta_d2
-                                    * prev_move.accel)
+        sin_theta_d2 = math.sqrt(0.5 * (1.0-junction_cos_theta))
+        cos_theta_d2 = math.sqrt(0.5 * max(1.0+junction_cos_theta, 0.))
         # Apply limits
-        self.max_start_v2 = min(
-            R_jd * self.junction_deviation * self.accel,
-            R_jd * prev_move.junction_deviation * prev_move.accel,
-            move_centripetal_v2, prev_move_centripetal_v2,
+        max_start_v2 = min(
             extruder_v2, self.max_cruise_v2, prev_move.max_cruise_v2,
             prev_move.max_start_v2 + prev_move.delta_v2)
+        # Check that the radius of curvature for this junction is finite
+        if min(self.move_d, prev_move.move_d) < 99999999.9 * cos_theta_d2:
+            R_jd = sin_theta_d2 / (1. - sin_theta_d2)
+            # Approximated circle must contact moves no further than mid-move, and
+            # centripetal_v2 = 0.5 * tan_theta_d2 * move_d * accel
+            quarter_tan_theta_d2 = .25 * sin_theta_d2 / cos_theta_d2
+            move_centripetal_v2 = quarter_tan_theta_d2 * self.delta_v2
+            prev_move_centripetal_v2 = quarter_tan_theta_d2 * prev_move.delta_v2
+            max_start_v2 = min(
+                max_start_v2,
+                R_jd * self.junction_deviation * self.accel,
+                R_jd * prev_move.junction_deviation * prev_move.accel,
+                max_start_v2, move_centripetal_v2, prev_move_centripetal_v2)
+        self.max_start_v2 = max_start_v2
         self.max_smoothed_v2 = min(
-            self.max_start_v2
-            , prev_move.max_smoothed_v2 + prev_move.smooth_delta_v2)
+            max_start_v2, prev_move.max_smoothed_v2 + prev_move.smooth_delta_v2)

What do you think about it?

How about we just remove them with something like:

My only concern is that this code does not entirely eliminate a possibility to get infinity along the calculations, e.g. 10. / 1e-320 == inf. So, while not very likely, move_jd_v2, prev_move_jd_v2, move_centripetal_v2, or prev_move_centripetal_v2 can still get inf value. You seem to be opposed of that. Otherwise, we could go with either a variant of your patch or a variant of @Wulfsta's proposal (I think conditionally computing the values for R_jd and tan_theta_d2 or switching to inf is a bit more readable, at least it gives a reader an assurance that this is an intended behavior of the code).

As was reported in the other PR, the problematic code snippet (presented in #6723 (comment)) is very similar to what the new test executes from time to time, leading to the wrong velocities and timing of moves.

So the new resonance testing code effectively issues a SET_VELOCITY_LIMIT ACCEL=0.1?

Yes, more or less, as evident by this chart (note very small non-zero accelerations occasionally in the beginning, for example):

[Wulfsta]

My only concern is that this code does not entirely eliminate a possibility to get infinity along the calculations, e.g. 10. / 1e-320 == inf. So, while not very likely, move_jd_v2, prev_move_jd_v2, move_centripetal_v2, or prev_move_centripetal_v2 can still get inf value. You seem to be opposed of that. Otherwise, we could go with either a variant of your patch or a variant of @Wulfsta's proposal (I think conditionally computing the values for R_jd and tan_theta_d2 or switching to inf is a bit more readable, at least it gives a reader an assurance that this is an intended behavior of the code).

Yes, in this case the min call that is sure to include values less than Inf should be considered defensive programming against this propagation. And in my opinion, it makes more sense to plug the hole (with the point at infinity, so not “plugging the hole” in the sense you might normally see it for this sort of piecewise operation) in this function, and not an arbitrarily small interval (as determined by some constant in the code) around the hole. Overall, I prefer the approach of checking for finite curvature - though maybe I am missing some downside as I have not had a lot of time to think it over.

[KevinOConnor]

Thanks for the additional information.

I understand that you meant R_jd-related calculations here, right?

Yes, that specific example was for R_jd.

What do you think about it?

It generally seems fine - merging the two branches and using delta_v2 seems like a good idea. I guess I'm a little uncomfortable by the new 99999999.9 constant though. Basically it does if cos_theta_d2 > 0.00000001 * move_d and I think I'd be more comfortable with if cos_theta_d2 > 0. .

My only concern is that this code does not entirely eliminate a possibility to get infinity along the calculations

The min() a few lines down would make an Inf moot. As a general goal, I don't like to encourage Infs, but I don't think we can reasonably guarantee they'll never occur. That said, if @Wulfsta's proposal for using explicit Infs is preferred then I'm okay with that.

Yes, more or less, as evident by this chart

Okay, thanks. I think I understand it better now.

Cheers,
-Kevin

[KevinOConnor]

So, your patch but without the constants might look something like:

--- a/klippy/toolhead.py
+++ b/klippy/toolhead.py
@@ -64,32 +64,32 @@ class Move:
             return
         # Allow extruder to calculate its maximum junction
         extruder_v2 = self.toolhead.extruder.calc_junction(prev_move, self)
+        max_start_v2 = min(
+            extruder_v2, self.max_cruise_v2, prev_move.max_cruise_v2,
+            prev_move.max_start_v2 + prev_move.delta_v2)
         # Find max velocity using "approximated centripetal velocity"
         axes_r = self.axes_r
         prev_axes_r = prev_move.axes_r
         junction_cos_theta = -(axes_r[0] * prev_axes_r[0]
                                + axes_r[1] * prev_axes_r[1]
                                + axes_r[2] * prev_axes_r[2])
-        if junction_cos_theta > 0.999999:
-            return
-        junction_cos_theta = max(junction_cos_theta, -0.999999)
-        sin_theta_d2 = math.sqrt(0.5*(1.0-junction_cos_theta))
-        R_jd = sin_theta_d2 / (1. - sin_theta_d2)
-        # Approximated circle must contact moves no further away than mid-move
-        tan_theta_d2 = sin_theta_d2 / math.sqrt(0.5*(1.0+junction_cos_theta))
-        move_centripetal_v2 = .5 * self.move_d * tan_theta_d2 * self.accel
-        prev_move_centripetal_v2 = (.5 * prev_move.move_d * tan_theta_d2
-                                    * prev_move.accel)
+        sin_theta_d2 = math.sqrt(max(0.5*(1.0-junction_cos_theta), 0.))
+        cos_theta_d2 = math.sqrt(max(0.5*(1.0+junction_cos_theta), 0.))
+        one_minus_sin_theta_d2 = 1. - sin_theta_d2
+        if cos_theta_d2 > 0. and one_minus_sin_theta_d2 > 0.:
+            R_jd = sin_theta_d2 / one_minus_sin_theta_d2
+            move_jd_v2 = R_jd * self.junction_deviation * self.accel
+            pmove_jd_v2 = R_jd * prev_move.junction_deviation * prev_move.accel
+            # Approximated circle must contact moves no further than mid-move
+            quarter_tan_theta_d2 = .25 * sin_theta_d2 / cos_theta_d2
+            move_centripetal_v2 = quarter_tan_theta_d2 * self.delta_v2
+            pmove_centripetal_v2 = quarter_tan_theta_d2 * prev_move.delta_v2
+            max_start_v2 = min(max_start_v2, move_jd_v2, pmove_jd_v2,
+                               move_centripetal_v2, pmove_centripetal_v2)
         # Apply limits
-        self.max_start_v2 = min(
-            R_jd * self.junction_deviation * self.accel,
-            R_jd * prev_move.junction_deviation * prev_move.accel,
-            move_centripetal_v2, prev_move_centripetal_v2,
-            extruder_v2, self.max_cruise_v2, prev_move.max_cruise_v2,
-            prev_move.max_start_v2 + prev_move.delta_v2)
+        self.max_start_v2 = max_start_v2
         self.max_smoothed_v2 = min(
-            self.max_start_v2
-            , prev_move.max_smoothed_v2 + prev_move.smooth_delta_v2)
+            max_start_v2, prev_move.max_smoothed_v2 + prev_move.smooth_delta_v2)
     def set_junction(self, start_v2, cruise_v2, end_v2):
         # Determine accel, cruise, and decel portions of the move distance
         half_inv_accel = .5 / self.accel

-Kevin

[dmbutyugin]

@KevinOConnor,

I updated this PR to match your suggestion, and I also tested it on top of the latest updated resonance testing code, and it seems to be working fine: