Implement solid angle measure for two-dimensional cones

[c13a4060-9513-4188-9d47-dacb5011e673]

Implement the solid angle measure (the plane angle) of a two-dimensional cone, using the dot-product (trigonometry) formula.

We do not assume nor compute the minimal description of the cone. When the input contains more than two vectors as generators of the cone, we decompose the cone into simplicial cones. The helper function simplicial_subcones_decomposition will be re-used in higher dimension cases.

For example, the solid angle of the upper half plane (generated by three rays, hence not simplicial) solid_angle_2d([[1,0],[0,1],[-1,0]]) is 1/2.

Depends on #33656

CC: @yuan-zhou @mkoeppe

Component: geometry

Keywords: solid angle

Author: Allison Fitisone, Yuan Zhou

Branch/Commit: u/gh-allisonfitisone/implement_solid_angle_measure_for_two_dimensional_cones @ 7a7c253

Issue created by migration from https://trac.sagemath.org/ticket/33752

[mkoeppe]

comment:1

The only non-simplicial cones in 2d are the non-pointed cones, i.e., cones that contain a nontrivial subspace. So the answer is either 0, 1/2, or 1; no triangulation is needed.

[yuan-zhou]

Description changed:

--- 
+++ 
@@ -1,3 +1,6 @@
-Implement the solid angle measure (the plane angle) of a two-dimensional cone, using the dot-product (trigonometry) formula. When the cone is not simplicial, it is decomposed into simplicial cones.
+Implement the solid angle measure (the plane angle) of a two-dimensional cone, using the dot-product (trigonometry) formula. 
+
+We do not assume nor compute the minimal description of the cone.
+When the input contains more than two vectors as generators of the cone, we decompose the cone into simplicial cones. The helper function `simplicial_subcones_decomposition` will be re-used in higher dimension cases.

 For example, the solid angle of the upper half plane (generated by three rays, hence not simplicial) `solid_angle_2d([[1,0],[0,1],[-1,0]])` is 1/2.

[yuan-zhou]

Branch: u/yzh/implement_solid_angle_measure_for_two_dimensional_cones

[yuan-zhou]

comment:4

Merged in #33656.

[yuan-zhou]

Commit: 29aa49e

[c13a4060-9513-4188-9d47-dacb5011e673]

Changed branch from u/yzh/implement_solid_angle_measure_for_two_dimensional_cones to u/gh-allisonfitisone/implement_solid_angle_measure_for_two_dimensional_cones

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

d4b7c16

fixed typos/formatting]

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 29aa49e to d4b7c16

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d4b7c16 to 1e62f8b

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

f31325b	wip
fe77320	change output type
72df52e	Merge branch 't/33656/implement_solid_angle_measure_for_two_dimensional_simplicial_cones' into t/33752/implement_solid_angle_measure_for_two_dimensional_cones
1e62f8b	updates to logging and output type

[yuan-zhou]

comment:8

When P.is_exact() or isinstance(P, sage.rings.abc.SymbolicRing), the return type should be SymbolicConstantsSubring element

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 1e62f8b to d3ad665

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

3c9aa58	modified output type and added tests to docstring
d3ad665	wip

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d3ad665 to 800d88d

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

800d88d

altered logging format

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 800d88d to a8e1152

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

9d5ab98	change output type for input in symbolic ring
165d13f	Merge branch 't/33656/implement_solid_angle_measure_for_two_dimensional_simplicial_cones' into t/33752/implement_solid_angle_measure_for_two_dimensional_cones
a8e1152	change output type

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a8e1152 to d23872e

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

1381c23	remove unnecessary line
4673707	wip
c48941b	output in Symbolic Constants Subring when input is exact
6535715	regroup the doctests for A and B with exact and inexact base rings.
d8a5f5a	fix doctests typo
76cb13b	Merge branch 't/33656/implement_solid_angle_measure_for_two_dimensional_simplicial_cones' into t/33752/implement_solid_angle_measure_for_two_dimensional_cones
d23872e	wip

[mkoeppe]

comment:14

+                    var('a b')
+                    eqn = (-a * v[0][0] - b * v[k][0] - v[j][0] == 0,
+                           -a * v[0][1] - b * v[k][1] - v[j][1] == 0,
+                           a >= 0,
+                           b >= 0)
+                    if solve(eqn, (a, b)):

Best not to use the general symbolics facility (var, solve) for linear algebra

[mkoeppe]

comment:15

Functions that looks like this:

+        matrices = []
+        for simplex in triangulation:
+            matrices.append(matrix(A[i] for i in simplex if i != origin))
+        return matrices

... are best rewritten as generators. https://realpython.com/introduction-to-python-generators/#understanding-generators

[mkoeppe]

comment:16

+    This example shows that the cone in `\RR^4` spanned by the rows of ``A``
+    (which is input as a list of lists) is actually a halfspace of affine
+    dimension `2`. The triangulation dissects it into three 2-d subcones::
+
+        sage: A_in = [[-3,0,5,0],[0,0,1,0],[-4,0,0,0],[-1,0,0,0],[0,0,-4,0]]
+        sage: simplicial_subcones_decomposition(A_in)
+        [
+        [-3  0  5  0]  [-3  0  5  0]  [-4  0  0  0]
+        [ 0  0  1  0], [-4  0  0  0], [ 0  0 -4  0]
+        ]

This decomposition looks correct, but I don't think this is the best way of computing the solid angle - as I explained in comment:1.

[mkoeppe]

comment:17

All of this can be replaced by just using Polyhedron(rays=A.rows()). To see whether it is at least a halfplane, check lines()

[c13a4060-9513-4188-9d47-dacb5011e673]

comment:18

When using Polyhedron(rays=A.rows()), we can't do computations over SR (see below). Any suggestions for a workaround?

sage: A = matrix([[sqrt(2),0],[0,1]])
sage: solid_angle_2d(A)
ValueError: no default backend for computations with Symbolic Ring

Replying to @mkoeppe:

All of this can be replaced by just using Polyhedron(rays=A.rows()). To see whether it is at least a halfplane, check lines()

[mkoeppe]

comment:19

If you only need algebraic numbers, you can use A = matrix(AA, [[sqrt(2),0],[0,1]])

[yuan-zhou]

comment:20

Right, if the input is of type AA, then Polyhedron can deal with it, and the output is also of type AA. The question is, what if the input is just A = matrix([[sqrt(2),0],[0,1]])? Would it be better to avoid Polyhedron and just use more elementary operations like rank? Replying to @mkoeppe:

If you only need algebraic numbers, you can use A = matrix(AA, [[sqrt(2),0],[0,1]])

[mkoeppe]

comment:21

Replying to @yuan-zhou:

Right, if the input is of type AA, then Polyhedron can deal with it, and the output is also of type AA. The question is, what if the input is just A = matrix([[sqrt(2),0],[0,1]])?

Well, you can try to convert it to a matrix space over AA, and if that fails, just raise an error.

Would it be better to avoid Polyhedron and just use more elementary operations like rank?

Do you have an application for non-algebraic irrational input?

If yes, then yes, one would have to either avoid Polyhedron or extend Polyhedron so it can handle well-behaved symbolic elements. After all, this error:

sage: Polyhedron(vertices=[[1], [sqrt(2)]], base_ring=SR, backend='field')
ValueError: the 'field' backend for polyhedron cannot be used with non-exact fields

is artificial - it's there to protect users from weird stuff that can happen.

If no, then maybe this generality is not worth working on?

[yuan-zhou]

comment:22

You are right, I don't have an application where the input is non-algebraic irrationals.

I agree that symbolics facility (var and solve) should be avoid.

However, I thought that checking the special cases of half space and whole space is straightforward enough, which make the use of Polyhedron and the conversion to AA seem overkill. It's just some if else statements, and with that, the input A = matrix([[sqrt(2),0],[0,1]]) can work. In my option, it's unusual to input matrix(AA, [[sqrt(2),0],[0,1]]).

Also, if we convert the inputs matrix to over AA, what would be the right type of the output? Do we always have AA as output type, never the Symbolic ring or Symbolic Constants Subring?

[mkoeppe]

comment:23

Replying to @yuan-zhou:

Also, if we convert the inputs matrix to over AA, what would be the right type of the output? Do we always have AA as output type, never the Symbolic ring or Symbolic Constants Subring?

Yes, it would make sense to convert it back to the ring where it came from if AA is used internally. (That's how we do this for the Polyhedron(backend='normaliz'), by the way.)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

e3a91f7	wip
7a7c253	add QQ

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d23872e to 7a7c253

[yuan-zhou]

comment:25

Do you mean to use number_field_elements_from_algebraics from sage.rings.qqbar? I still don't see how to convert AA(sqrt(2)) (interval AA number) back to sqrt(2) as output, if the input was in the symbolic ring.

Replying to @mkoeppe:

Replying to @yuan-zhou:

Also, if we convert the inputs matrix to over AA, what would be the right type of the output? Do we always have AA as output type, never the Symbolic ring or Symbolic Constants Subring?

Yes, it would make sense to convert it back to the ring where it came from if AA is used internally. (That's how we do this for the Polyhedron(backend='normaliz'), by the way.)

[mkoeppe]

comment:26

First check if Polyhedron(..., backend='normaliz') is able to handle the kinds of numbers you'd like to handle and gives back SR elements.

[yuan-zhou]

comment:27

Sure, using backend='normaliz' solves most cases. The problem is that this is an optional backend. Do we really want to use an optional package that most users don't install to deal with something so elementary -- 2 very special cases in 2D?

Replying to @mkoeppe:

First check if Polyhedron(..., backend='normaliz') is able to handle the kinds of numbers you'd like to handle and gives back SR elements.

[mkoeppe]

comment:28

Your code would still work without normaliz for rational input.

[yuan-zhou]

comment:29

Yes, that's right. As it's quite usual to have sqrt(2), sqrt(3) in rays (45°, 60°), I think the code should work for these cases regardless of having normaliz or not.

Replying to @mkoeppe:

Your code would still work without normaliz for rational input.

[mkoeppe]

comment:30

I think some more thought on the use cases is needed. Where do your irrational inputs come from? If you already have angles, there's no need to recompute them

[mkoeppe]

comment:31

Replying to @yuan-zhou:

Sure, using backend='normaliz' solves most cases.

I mean I know that it can deal with sqrt(2). Have you checked other cases? What are "most", are there exceptions that are motivated by an application?

[yuan-zhou]

comment:32

I don't really have use cases in 2D, 2D angles are trivial...

It feels weird that if a Calculus student wants to compute the angle between ray (1, 0) and ray (0, sqrt(3)), they must have normaliz installed.

My thought is not to use Polyhedron in 2d, but just check for the two special cases.

Alternatively, we can use P = Polyhedron(rays=A.rows(), base_ring=AA) which doesn't call normaliz. But then, I don't know how to convert the AA result back to symbolic.

[mkoeppe]

comment:33

Replying to @yuan-zhou:

Alternatively, we can use P = Polyhedron(rays=A.rows(), base_ring=AA) which doesn't call normaliz. But then, I don't know how to convert the AA result back to symbolic.

You're missing that the conversion code that is in our Normaliz backend does not need Normaliz. We are using it to communicate with Normaliz.

[yuan-zhou]

comment:34

I know that "conversion code that is in our Normaliz backend does not need Normaliz." But how does it help with 2d solid angle?

Replying to @yuan-zhou:

Alternatively, we can use P = Polyhedron(rays=A.rows(), base_ring=AA) which doesn't call normaliz. But then, I don't know how to convert the AA result back to symbolic.

You're missing that the conversion code that is in our Normaliz backend does not need Normaliz. We are using it to communicate with Normaliz.

[mkoeppe]

comment:35

See #34479

[yuan-zhou]

comment:36

I agree with #34479. That ticket will become more useful when we write code for solid angle computation in higher dimension. In 2d, it is like shooting butterfly with rifles :P, but why not...