计算两个向量之间的有符号角度的数值稳定方法

huangapple go评论68阅读模式
英文:

Numerically stable approach to computing a signed angle between vectors

问题

以下是您提供的文本的翻译:

我的第一个Stack Exchange问题涉及几何中的数值稳定性。我编写了下面的这个函数,用于计算从一个向量到另一个向量的有符号角度,灵感来自此处提出的想法,我对此表示感激。

在这个函数中,使用三重积来从第一部分得到原始锐角,然后得到从第一个向量到第二个向量的有向角度。

import numpy as np
from numpy import ndarray, sign
from math import pi

Vec3D = ndarray

def directed_angle(
        from_vec: Vec3D,
        to_vec: Vec3D,
        normal: Vec3D,
        debug=False):

    """返回从第一个向量到第二个向量开放的[0, 2π)有向角度,
    需要第一个两个向量之前线性独立的第三个向量,该向量用于提供
    首先计算的锐角的方向性(通常是未定向的)"""

    assert np.any(from_vec) and np.any(to_vec), \
           f'from_vec、to_vec和normal都不得为零向量,但from_vec: {from_vec},to_vec: {to_vec}'

    magnitudes = np.linalg.norm(from_vec) * np.linalg.norm(to_vec)
    dot_product = np.matmul(from_vec, to_vec)
    angle = np.arccos(np.clip(dot_product / magnitudes, -1, 1))

    # 三重积反映了从第一个向量到第二个向量开放的角度的符号,
    # 基于输入的normal值提供方向性,否则完全不存在方向性
    triple_product = np.linalg.det(np.array([from_vec, to_vec, normal]))

    if triple_product < 0:
        result_angle = angle
    else:
        result_angle = 2*pi - angle

    # 将输出范围从(0, 2π]翻转到[0, 2π),这是我们语义的首选
    if result_angle == 2*pi:
        result_angle = 0

    if debug:
        print(f'\nfrom_vec: {from_vec}, '
              f'\nto_vec: {to_vec}'
              f'\nnormal: {normal}, '
              f'\n未定向角度: {angle},'
              f'\n未定向角度/pi: {angle/pi}, '
              f'\n三重积: {triple_product}'
              f'\n三重积符号: {sign(triple_product)}'
              f'\n有向角度/pi: {result_angle/pi}')

    return result_angle


def test():

    v0 = np.array([1, 0, 0])
    normal = np.array([0, 0, 1])

    vecs = dict(
        same_direction = np.array([1, 0, 0]),
        opposite_direction = np.array([-1, 0, 0]),
        close_to_same_direction_from_side_1 = np.array([1, 0.00000001, 0]),
        close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))

    for desc, v in vecs.items():
        print(f'\n{desc}:')
        directed_angle(v0, v, normal, debug=True)

当像上面的示例测试代码一样给定的第二个向量的角度非常接近于第一个给定向量的角度无论从哪一侧看而不是得到结果接近于0和接近于2π的结果函数的结果实际上是零因此这两个在语义上不同的向量方向产生相同的输出角度

您如何处理使该函数在这个意义上在数值上稳定呢

假设有向角从0增长到接近2π的数字您将如何避免函数的输出角度在它无限接近完整圆周的时候跳至0只要它的输出在数值上不精确我们都没问题但如果它从接近2π跳至0这种不连续的值跳跃将对我的应用程序中的输出解释以及其他许多情况产生破坏性影响

显然我们可以经验性地找到函数输出接近于2π的角度时候崩溃到0的上界边界 - 实际上下面附带了一个用于此的函数 - 但当输入是产生该角度的向量时它完全没有帮助我们无法从向量中得知它将跳转到0因此当函数的结果是0时我们无法确定实际角度是否为0还是接近于2π的值)。

我猜您可以将术语有向角度替换为顺时针或逆时针角度之一

诊断函数如下

```python
def test_monotonicity(debug=False, interpolation_interval=1e-4):

    """测试函数是否在其域的区间内([0, 2π))是单调的(根据使用的插值间隔),
    并且正确的"""

    v0 = np.array([1, 0, 0])
    normal = np.array([0, 0, 1])

    # 在0到2π的范围内以0.1间隔
    angles = np.arange(0, 2*pi, step=interpolation_interval)
    prev_angle = None
    for angle in angles:

        # 创建一个表示当前角度离v0的向量,使其随着“angle”的增加而顺时针移动
        # 离开v0(假设您从正常向量的角度看两个向量,否则方向性将无效)
        vec = np.array([np.cos(angle), -np.sin(angle), 0])
        result_angle = directed_angle(v0, vec, normal)

        if debug:  print(f'angle/pi: {angle

<details>
<summary>英文:</summary>

My first Stack Exchange question regarding numeric stability in geometry. I wrote this function (below) for computing a signed angle from one vector to the other, inspired by [the idea presented here](https://math.stackexchange.com/a/1139228/66486), to which I am indebted. 

In the function, [the triple product](https://en.wikipedia.org/wiki/Triple_product) is used for getting from the raw acute angle that you get on the first part, to the directed angle opening up from the first vector to the second.

&lt;!-- language: lang-py --&gt;
    import numpy as np
    from numpy import ndarray, sign
    from math import pi
    
    Vec3D = ndarray
    
    def directed_angle(
            from_vec: Vec3D,
            to_vec: Vec3D,
            normal: Vec3D,
            debug=False):
    
        &quot;&quot;&quot; returns the [0, 2&#120587;) directed angle opening up from the first to the second vector,
        given a 3rd vector linearly independent of the first two vectors, which is used to provide
        the directionality of the acute (regular, i.e. undirected) angle firstly computed &quot;&quot;&quot;
    
        assert np.any(from_vec) and np.any(to_vec), \
               f&#39;from_vec, to_vec, and normal must not be zero vectors, but from_vec: {from_vec}, to_vec: {to_vec}&#39;
    
        magnitudes = np.linalg.norm(from_vec) * np.linalg.norm(to_vec)
        dot_product = np.matmul(from_vec, to_vec)
        angle = np.arccos(np.clip(dot_product / magnitudes, -1, 1))
    
        # the triplet product reflects the sign of the angle opening up from the first vector to the second,
        # based on the input normal value providing the directionality which is otherwise totally abscent
        triple_product = np.linalg.det(np.array([from_vec, to_vec, normal]))
    
        if triple_product &lt; 0:
            result_angle = angle
        else:
            result_angle = 2*pi - angle
    
        # flip the output range from (0, 2&#120587;] to [0, 2&#120587;) which is the prefernce for our semantics
        if result_angle == 2*pi:
            result_angle = 0
    
        if debug:
            print(f&#39;\nfrom_vec: {from_vec}, &#39;
                  f&#39;\nto_vec: {to_vec}&#39;
                  f&#39;\nnormal: {normal}, &#39;
                  f&#39;\nundirected angle: {angle},&#39;
                  f&#39;\nundirected angle/pi: {angle/pi}, &#39;
                  f&#39;\ntriple product: {triple_product}&#39;
                  f&#39;\ntriple product sign: {sign(triple_product)}&#39;
                  f&#39;\nresult_angle/pi: {result_angle/pi}&#39;)
    
        return result_angle
    
    
    
    def test():
    
        v0 = np.array([1, 0, 0])
        normal = np.array([0, 0, 1])
    
        vecs = dict(
            same_direction = np.array([1, 0, 0]),
            opposite_direction = np.array([-1, 0, 0]),
            close_to_same_direction_from_side_1 = np.array([1, 0.00000001, 0]),
            close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))
    
        for desc, v in vecs.items():
            print(f&#39;\n{desc}:&#39;)
            directed_angle(v0, v, normal, debug=True)

When as in the example test code above, the given second vector&#39;s angle is extremely close to the first given vector&#39;s angle from either side of it, instead of getting the resulting angle  respectively very close to 0 and very close to 2&#120587;, the function&#39;s result is actually zero, so these two semantically different vector directions yield the same output angle. 

How would you approach making this function numerically stable in that sense? 

Assume the directed angle grows from 0 up to numbers approaching 2&#120587;, how would you avoid the output angle of the function jumping to 0 as it infinitesimally approaches the full circle directed angle? as long as its output is infinitesimally imprecise we are fine, but if it goes from close to 2&#120587; to 0 as it infinitesimally approaches 2&#120587; ― this kind of discontinuous value jump will be devastating for the interpretation of its output in my application as well as I guess in many other cases. 

Obviously we can empirically find upper-bound boundaries for where the function output for angles approaching 2&#120587; collapses to 0 ― a function for that is actually attached below ― but this doesn&#39;t help at all when an input is the vector which yields that angle (we cannot tell from the vector that it is going to yield a flip to 0, so when we have 0 as the resulting angle of the function ― we cannot tell if the real angle was 0 or rather a value approaching 2&#120587;).

I guess you can replace the term _directed angle_ with any of the terms _clockwise (or counter-clockwise) angle_. 

diagnostic functions follow.

&lt;!-- language: lang-py --&gt;
    def test_monotonicity(debug=False, interpolation_interval=1e-4):
    
        &quot;&quot;&quot; test that the function is monotonic (up to the used interpolation interval),
        and correct, across the interval of its domain ([0, 2&#120587;)) &quot;&quot;&quot;
    
        v0 = np.array([1, 0, 0])
        normal = np.array([0, 0, 1])
    
        # range from 0 to 2&#120587; at 0.1 intervals
        angles = np.arange(0, 2*pi, step=interpolation_interval)
        prev_angle = None
        for angle in angles:
    
            # make a vector representing the current angle away from v0, such that it goes clockwise
            # away from v0 as `angle` increases (clockwise assuming you look at the two vectors from
            # the perspective of the normal vector, otherwise directionality is void)
            vec = np.array([np.cos(angle), -np.sin(angle), 0])
            result_angle = directed_angle(v0, vec, normal)
    
            if debug:  print(f&#39;angle/pi: {angle/pi}, computed angle/pi: {result_angle/pi}&#39;)
    
            if prev_angle is None:
                assert angle == 0
            else:
                assert angle &gt; prev_angle
                drift = result_angle - angle
                if angle == result_angle:
                    pass
                else:
                    print(f&#39;angle/pi: {angle/pi}, result_angle/pi: {result_angle/pi}, drift/pi: {drift/pi}&#39;)
                    assert isclose(result_angle, angle, rel_tol=1e-6)
    
            prev_angle = angle
    
    
    def test_demonstrating_the_concern():
    
        &quot;&quot;&quot; demonstration of sufficiently small angles from both sides of the reference vector (v0)
        collapsing to zero, which is a problem only when the angle should be close to 2 pi &quot;&quot;&quot;
    
        v0 = np.array([1, 0, 0])
        normal = np.array([0, 0, 1])
    
        vecs = dict(
            same_direction = np.array([1, 0, 0]),
            opposite_direction = np.array([-1, 0, 0]),
            close_to_same_direction_from_side_1 = np.array([1, +0.00000001, 0]),
            close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))
    
        expected_results = [0, 1*pi, None, None]
    
        for (desc, v), expected_result in zip(vecs.items(), expected_results):
            print(f&#39;\n{desc}:&#39;)
            v = v / np.linalg.norm(v)
            result = directed_angle(v0, v, normal, debug=True)
            if expected_result is not None:
                assert(result == expected_result)
    
    
    def test_angle_approaching_2pi(epsilon=1e-7, interpolation_interval=1e-10, verbose=True):
    
        &quot;&quot;&quot; stress test the angle values approaching 2pi where the result flips to zero &quot;&quot;&quot;
    
        v0 = np.array([1, 0, 0])
        normal = np.array([0, 0, 1])
    
        # range from 2&#120587;-epsilon to close to 2&#120587;
        angles = np.arange(2*pi-epsilon, 2*pi, step=interpolation_interval)
        prev_angle = None
        for angle in angles:
    
            # vector representing the current angle away from v0, clockwise
            vec = np.array([np.cos(angle), -np.sin(angle), 0])
            result_angle = directed_angle(v0, vec, normal)
    
            if prev_angle is None:
                pass
            else:
                assert angle &gt; prev_angle
                drift = result_angle - angle
                if angle == result_angle:
                    pass
                else:
                    if verbose: print(f&#39;angle/pi: {angle / pi}, result_angle/pi: {result_angle / pi}, drift/pi: {drift / pi}&#39;)
                    if result_angle == 0:
                        print(f&#39;function angle output hit 0 at angle: {angle};&#39;
                              f&#39;\nangle/pi: {angle / pi}&#39;
                              f&#39;\nangle distance from 2&#120587;: {2*pi - angle}&#39;)
                        break
                    else:
                        assert isclose(result_angle, angle, rel_tol=1e-6)
    
            prev_angle = angle


And all that said, I&#39;m happy to also have insights about possible other numerical stability issues which may arise in this function arising in the calculation steps that it is making. I have seen comments about preferring arctan to arccos and variations, but wouldn&#39;t know whether they fully apply to numpy without introducing other stability drawbacks. 

Also thankful for comments to any general methodology for thinking numerical stability or guidelines specific to numpy for safely writing numerically stable computation in the realm of 3D geometry over arbitrary inputs; although I come to realize that different applications will worry about different ranges where inaccuracies become detrimental to their specific logic. 

As my inputs to this function come in from machine learning predictions, I can&#39;t avoid extreme vector values even though their probability is low at any single invocation of the function.

Thanks!




</details>


# 答案1
**得分**: 1

将此计算方法调整为NumPy格式后我所有的数值稳定性问题似乎都得到了很好的满足

```python
angle = 2 * np.arctan2(
                norm(from_vec / norm(from_vec) - to_vec / norm(to_vec)),
                norm(from_vec / norm(from_vec) + to_vec / norm(to_vec)))

完整的函数代码和测试如下:

import numpy as np
from numpy import ndarray, sign
from math import pi, isclose
from numpy.linalg import norm

Vec3D = ndarray

def directed_angle_base(
        from_vec: Vec3D,
        to_vec: Vec3D,
        normal: Vec3D,
        debug=False):

    assert np.any(from_vec) and np.any(to_vec), \
           f'from_vec, to_vec, and normal must not be zero vectors, but from_vec: {from_vec}, to_vec: {to_vec}'

    angle = 2 * np.arctan2(
        norm(from_vec / norm(from_vec) - to_vec /norm(to_vec)),
        norm(from_vec / norm(from_vec) + to_vec / norm(to_vec)))

    triple_product = np.linalg.det(np.array([from_vec, to_vec, normal]))

    if triple_product < 0:
        result_angle = angle
    else:
        result_angle = 2*pi - angle

    if result_angle == 2*pi:
        if debug: print('angle output flipped from 2π to zero')
        result_angle = 0

    if debug:
        print(f'\nfrom_vec: {from_vec}, '
              f'\nto_vec: {to_vec}'
              f'\nnormal: {normal}, '
              f'\nundirected angle: {angle},'
              f'\nundirected angle/pi: {angle/pi}, '
              f'\ntriple product: {triple_product}'
              f'\ntriple product sign: {sign(triple_product)}'
              f'\nresult_angle/pi: {result_angle/pi}')

    return result_angle


def test_monotonicity(debug=False, interpolation_interval=1e-4):
    # 省略测试代码部分

def test_near_2pi_concern():
    # 省略测试代码部分

def test_angle_approaching_2pi(epsilon=1e-11, interpolation_interval=1e-15, verbose=True):
    # 省略测试代码部分

def test_smallest_angle():
    # 省略测试代码部分

# 免责声明和其他注释已省略

请注意,我已将原文中的HTML实体编码(例如&lt;&gt;)还原为正常的Python代码。如果需要其他帮助,请随时告诉我。

英文:

Adapting this computation method to numpy, all my numerical stability concerns seem impeccably satisfied:

<!-- language: lang-py -->

angle = 2 * np.arctan2(
norm(from_vec / norm(from_vec) - to_vec / norm(to_vec)),
norm(from_vec / norm(from_vec) + to_vec / norm(to_vec)))

Full function code and tests:
<!-- language: lang-py -->

import numpy as np
from numpy import ndarray, sign
from math import pi, isclose
from numpy.linalg import norm
Vec3D = ndarray
def directed_angle_base(
from_vec: Vec3D,
to_vec: Vec3D,
normal: Vec3D,
debug=False):
&quot;&quot;&quot; returns the [0, 2&#120587;) directed angle opening up from the first to the second vector,
given a 3rd vector linearly independent of the first two vectors, which is used to provide
the directionality of the acute (regular, i.e. undirected) angle firstly computed. &quot;&quot;&quot;
assert np.any(from_vec) and np.any(to_vec), \
f&#39;from_vec, to_vec, and normal must not be zero vectors, but from_vec: {from_vec}, to_vec: {to_vec}&#39;
angle = 2 * np.arctan2(
norm(from_vec / norm(from_vec) - to_vec /norm(to_vec)),
norm(from_vec / norm(from_vec) + to_vec / norm(to_vec)))
# the triplet product reflects the sign of the angle opening up from the first vector to the second,
# based on the input normal value providing the directionality which is otherwise totally abscent
triple_product = np.linalg.det(np.array([from_vec, to_vec, normal]))
if triple_product &lt; 0:
result_angle = angle
else:
result_angle = 2*pi - angle
# flip the output range from (0, 2&#120587;] to [0, 2&#120587;) which is the prefernce for our semantics
if result_angle == 2*pi:
if debug: print(f&#39;angle output flipped from 2&#120587; to zero&#39;)
result_angle = 0
if debug:
print(f&#39;\nfrom_vec: {from_vec}, &#39;
f&#39;\nto_vec: {to_vec}&#39;
f&#39;\nnormal: {normal}, &#39;
f&#39;\nundirected angle: {angle},&#39;
f&#39;\nundirected angle/pi: {angle/pi}, &#39;
f&#39;\ntriple product: {triple_product}&#39;
f&#39;\ntriple product sign: {sign(triple_product)}&#39;
f&#39;\nresult_angle/pi: {result_angle/pi}&#39;)
return result_angle
def test_monotonicity(debug=False, interpolation_interval=1e-4):
&quot;&quot;&quot; tests that the function is monotonic (up to the used interpolation interval),
and correct, across the interval of its domain from 0 up to only 2&#120587;-ɛ, ɛ being
implied by the interpolation interval (so it doesn&#39;t hit the value flipping from
2&#120587;-ɛ to zero but only checks we are correct and monotonic up to a certain ɛ
away from that value &quot;&quot;&quot;
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
# range from 0 to 2&#120587; at 0.1 intervals
angles = np.arange(0, 2*pi, step=interpolation_interval)
prev_angle = None
for angle in angles:
# make a vector representing the current angle away from v0, such that it goes clockwise
# away from v0 as `angle` increases (clockwise assuming you look at the two vectors from
# the perspective of the normal vector, otherwise directionality is void)
vec = np.array([np.cos(angle), -np.sin(angle), 0])
result_angle = directed_angle_base(v0, vec, normal)
if debug:  print(f&#39;angle/pi: {angle/pi}, computed angle/pi: {result_angle/pi}&#39;)
if prev_angle is None:
assert angle == 0
else:
assert angle &gt; prev_angle
drift = result_angle - angle
if angle == result_angle:
pass
else:
print(f&#39;angle/pi: {angle/pi}, result_angle/pi: {result_angle/pi}, drift/pi: {drift/pi}&#39;)
assert isclose(result_angle, angle, rel_tol=1e-6)
prev_angle = angle
def test_near_2pi_concern():
&quot;&quot;&quot; demonstration of sufficiently small angles from both sides of the reference vector (v0)
collapsing to zero, which was a problem only when the angle was infinitesimally close to 2&#120587;
with the naive acos/atan formula, before switching to Velvel Kahan&#39;s formula like we have now
(https://stackoverflow.com/a/76722358/1509695) &quot;&quot;&quot;
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
vecs = dict(
same_direction = np.array([1, 0, 0]),
opposite_direction = np.array([-1, 0, 0]),
close_to_same_direction_from_side_1 = np.array([1, +0.00000001, 0]),  # should be 2&#120587;-ɛ
close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))  # should be +ɛ
expected_results = [0, 1*pi, None, None]
for (desc, v), expected_result in zip(vecs.items(), expected_results):
print(f&#39;\n{desc}:&#39;)
v = v / np.linalg.norm(v)
result = directed_angle_base(v0, v, normal, debug=True)
if expected_result is not None:
assert(result == expected_result)
def test_angle_approaching_2pi(epsilon=1e-11, interpolation_interval=1e-15, verbose=True):
&quot;&quot;&quot; stress test the angle values approaching 2pi where the result flips to zero &quot;&quot;&quot;
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
# range from 2&#120587;-epsilon to close to 2&#120587;
angles = np.arange(2*pi-epsilon, 2*pi, step=interpolation_interval)
prev_angle = None
for angle in angles:
if angle &gt; 2*pi:
print(f&#39;arange output value larger than 2&#120587;: {angle}, ignoring&#39;)
continue
# vector representing the current angle away from v0, clockwise
vec = np.array([np.cos(angle), -np.sin(angle), 0])
result_angle = directed_angle_base(v0, vec, normal)
if prev_angle is None:
pass
else:
assert angle &gt; prev_angle
drift = result_angle - angle
if verbose: print(f&#39;angle/pi: {angle / pi}, result_angle/pi: {result_angle / pi}, drift/pi: {drift / pi}&#39;)
if angle == result_angle:
pass
else:
if result_angle == 0:
print(f&#39;function angle output hit 0 at angle: {angle};&#39;
f&#39;\nangle/pi: {angle / pi}&#39;
f&#39;\nangle distance from 2&#120587;: {2*pi - angle}&#39;)
return angle, vec
else:
assert isclose(result_angle, angle, rel_tol=1e-6)
prev_angle = angle
def test_smallest_angle():
v0 = np.array([1, 0, 0])
angle = 1.999999999999999*pi
vec = np.array([np.cos(angle), -np.sin(angle), 0])
normal = np.array([0, 0, 1])
result_angle = directed_angle_base(v0, vec, normal, debug=True)
print(f&#39;angle/pi: {angle/pi}, result_angle/pi: {result_angle/pi}&#39;)
assert result_angle == angle

Disclaiming that the above tests do not test the monotonicity of the function at the smallest possible floating point intervals (but I assume that's granted anyway) and that you may or may not prefer other stability features for other implementations requiring a signed angle function. For me, it seems to solve the pain point of flipping to zero at values 2𝜋-ε for infinitesimal values of ε, while passing my other tests, and I hope not to be coming across other stability challenges as I further the use of this function for my application.

答案2

得分: 0

以下是代码部分的翻译:

# 第一个非解决方案:
# 定义一个函数,计算有向角度,覆盖了基本情况下计算结果为0的情况,
# 在这种情况下,如果检查显示两个向量不指向相同方向,则返回2π-ε。

def directed_angle(
    from_vec: Vec3D,
    to_vec: Vec3D,
    normal: Vec3D,
    epsilon=1e-7,
    debug=False):

    base_angle_calculation = directed_angle_base(from_vec, to_vec, normal, debug=debug)

    if base_angle_calculation == 0:
        # 同向向量?
        divided = from_vec / to_vec
        if isclose(divided[0], divided[1]) and isclose(divided[1], divided[2]):
            return 0
        else:
            # 弹回到一个任意的、由调用者选择的角度值2π-ε;
            # 这保持了所需的语义,即我们要避免从2π-ε翻转到零;
            # 但是如果epsilon设置得太大,我们就会在某些其他角度小于2π-ε的情况下,
            # 失去函数的单调性(即使后者大于前者,我们的函数对后者的结果将小于其他角度的结果)
            return 2*pi - epsilon

这是您提供的代码部分的翻译。如果您需要进一步的翻译或有其他问题,请告诉我。

英文:

First non-solution:

<!-- language: lang-py -->

def directed_angle(
from_vec: Vec3D,
to_vec: Vec3D,
normal: Vec3D,
epsilon=1e-7,
debug=False):
&quot;&quot;&quot; returns the directed angle while overriding for the case where the base 
calculation returns 0, in which case we return 2&#120587;-ε instead if checking shows 
that the two vectors do not point in the same direction &quot;&quot;&quot;
base_angle_calculation = directed_angle_base(from_vec, to_vec, normal, debug=debug)
if base_angle_calculation == 0:
# same direction vectors?
divided = from_vec / to_vec
if isclose(divided[0], divided[1]) and isclose(divided[1], divided[2]):
return 0
else:
# bounce back to an arbitrary, caller chosen, angle value of 2&#120587;-ε; this maintains
# the desired semantics whereby we want to avoid flipping from 2&#120587;-ε to zero; but
# it does mean that if epsilon is set too large, we give away the monotonicity of
# the function in the case that some other angle that does cause our function to
# flip to 0 is smaller than 2&#120587;-ε (even though the latter is larger than the former,
# our function result for it will be smaller than the other angle&#39;s result)
return 2*pi - epsilon

This is wrong because both:

  1. It requires figuring the threshold angle triggering the value flip from 2𝜋-ε to zero in advance (ignorantly about floating point standard implementations, I am not aware of guarantees for this threshold angle being the same on different processors or setups of numpy, so this implies a calibration step being necessary upon startup, to find that threshold angle).
  2. It's hard to tell whether the various ways of checking whether two vectors have the same direction ― which this approach hinges upon as seen in its code ― are themselves very numerically stable in the same context as this function would use them.

If this approach bears any promise, I am happy to learn.

答案3

得分: 0

第二个被放弃的想法:

"承认局限性并将其减小" -

  1. 对于给定的输入,两次使用基本角度计算函数 - 一次使用原始向量(v1),第二次使用v1的一个大致旋转版本(比如,v1旋转了𝜋/4)。对于第二次调用的基本函数的输出,将结果旋转回𝜋/4。

  2. 现在我们有两个结果。这两个结果中只有一个可能接近于从2𝜋-ε到零的“翻转区域”。

  3. 因此,我们可以通过检查这两个结果是否大致相隔2𝜋-ε来完成。

然而,这再次假设在旋转向量时不会遇到数值稳定性问题,所以被放弃...

英文:

Second abandoned idea:

"Acknowledging the limitation and reducing to it" ―

  1. For a given input, feed the base angle calculation function twice ― once with the original vector (v1) and a second time with a substantially rotated version of v1 (say, v1 rotated by 𝜋/4). For the base function's output for the second call, rotate the result back by 𝜋/4.

  2. So we have two results now. Only one of the two results may have fallen close to the "flip zone" from 2𝜋-ε to zero.

  3. Hence we can finish by checking whether the two results are approximately 2𝜋-ε apart or not.

However this again, assumes you do not bump into numerical stability issues in rotating the vector, so abandoned ...

huangapple
  • 本文由 发表于 2023年7月18日 06:55:27
  • 转载请务必保留本文链接:https://go.coder-hub.com/76708551.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定