Containers for families of geometric objects

[mkoeppe]

We define APIs for static and dynamic containers storing finite families of geometric objects, extending Container, Set, MutableSet, Mapping, MutableMapping https://docs.python.org/3/library/collections.abc.html

The family must define a collection of "measurement maps" B_i, each sending an object to an InternalRealInterval.

The product of these intervals can be thought of as a "generalized bounding box".

Ideally the maps would be inclusion-preserving.

The implementation is provided by rtree/libspatialindex (#32000).

The Sage-specific class will take care of:

• replacing an InternalRealInterval by a rescaling or inclusion-preserving overestimation, making the coordinates suitable for the underlying library.
• providing accelerated set operations such as:
  • a fast path for computing empty intersection when the result will be empty,
  • a fast path for __contains__ when the result will be False.

Geometric lookup operations to be supported:

• ...

Applications:

• RationalPolyhedralFan
• PolyhedralComplex (PolyhedralComplex #31748)
• Families/posets of (open) subsets of manifolds (Meta-ticket: Families, posets, complexes of manifold subsets #31740): A natural choice of the maps B_i comes from the coordinates of covering charts, using (-oo, oo) for coordinates for which no coordinate definition for the subset is known

Depends on #34277

CC: @DRKWang

Component: geometry

Work Issues: Rebase on #34277

Author: Binshuai Wang

Branch/Commit: public/32170 @ 05ea99d

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

[DRKWang]

comment:1

The precision of interval values in Package libspatialindex is double.

They did not claim it explicitly, but I found it from

• https://rtree.readthedocs.io/en/latest/class.html#rtree.index.Index
• https://github.com/libspatialindex/libspatialindex/blob/master/include/spatialindex/capi/DataStream.h.

[mkoeppe]

comment:2

OK, the next step would be to find out whether the library does any arithmetic with these double values (in this case we would have to check whether they apply directed rounding correctly), or whether they only use comparisons on these double values.

[DRKWang]

comment:4

The library libspatialindex uses arithmetic operations, including +, - and *, with those double values in order to support some functions, such as

• LineSegment::Collinear( https://libspatialindex.org/en/latest/doxygen/LineSegment_8cc_source.html#l00439),

• LineSegment::getArea(https://github.com/libspatialindex/libspatialindex/blob/master/src/spatialindex/Region.cc)

and etc.

Besides, the library libspatialindex also uses std::numeric_limits<double>::epsilon() to deal with machine precision issue, such as Region::touchesPoint(const Point& p) (https://github.com/libspatialindex/libspatialindex/blob/master/src/spatialindex/Region.cc)

std::numeric_limits<double>::max();(https://github.com/libspatialindex/libspatialindex/blob/master/src/spatialindex/Point.cc).

[DRKWang]

comment:6

A few test codes have been made for testing arithmetic precision. The conclusions are as follows, more details and codes can be found from this(https://github.com/DRKWang/Data-structure-for-piecewise-linear-function/blob/main/log35.ipynb).

• The dimension of the index database should always be greater than 1. (0 or 1 dimension will fail).

• The arithmetic precision is 16 decimal digits, which comes from the double type's precision.

• The minimal capturing of arithmetic precision is 324 decimal digits, which means the minimal difference to distinguish two numbers is 1*10^{-324}. otherwise, they will be viewed as the same number.

• The intersection for retrieving is a closed intersection. (-0.000) is also considered as the same number as 0, but both representations will be displayed at the same time.

• Rounding. The round processing is a little complex and untraceable, it is neither rounding up or rounding down.

[DRKWang]

comment:7

I have made some debugs for this library. I would like to list the key points below. The entire log file can be found here https://github.com/DRKWang/Data-structure-for-piecewise-linear-function/blob/main/logsfile/log36.ipynb.

1. The python wrapper part for rtree was only used for passing the values and data to 2 basic functions, core.rt.Index_InsertData() and core.rt.Index_Intersects_obj(). Those two functions are stored here https://github.com/Toblerity/rtree/blob/master/rtree/core.py. Thus, there is no arithmetic error in this part.

2. The corresponding part in libspatialindex for core.rt.Index_InsertData() and core.rt.Index_Intersects_obj() are placed in line 499 and line 678 in https://github.com/libspatialindex/libspatialindex/blob/master/src/capi/sidx_api.cc .

3. In the Index_InsertData() in C++, the author uses the concepts "points" and "boxes" to refine the input information, and then call the function index().insertData(**) to compute. As we discussed last time, the author mistakenly uses the std::numeric_limits<double>::epsilon() as the range of approximation, that is, if the bound difference is within this range, then they will be seen as a single point instead of a box. According to this note, https://stackoverflow.com/questions/48133572/what-can-stdnumeric-limitsdoubleepsilon-be-used-for, we may correct it by multiplying the magnitudes of the comparison. However, because it is for scientific computation, we should get rid of approximation as much as possible.

4. Except for the above place, the author also mistakenly has used std::numeric_limits<double>::epsilon() in many other places. I think we could inform the author to check them deeply.

[DRKWang]

comment:8

• The extensions, .idx and .dat, are two files used to serialize the rtree to disk, according to https://rtree.readthedocs.io/en/latest/tutorial.html.

• Those two files are created by #include <fstream> internally, based on this source code, https://github.com/libspatialindex/libspatialindex/blob/master/src/storagemanager/DiskStorageManager.cc.

• The .idx is used to store vital information like the page size, the next available page, a list of empty pages, and the sequence of pages associated with every entity ID.There is no note for the detail or purpose of .dat file. Thus, I made a test where putting 50 boxes into this library. The shows that the .dat is way larger than .idx, which indicates the .dat is used to store the real data, whereas .idx is used for storage the rest of it.

• The correct way to close and reopen those serialized files can be found here, https://gis.stackexchange.com/questions/254781/saving-python-rtree-spatial-index-to-file

• Except for DiskStorageManager, it also has a MemoryStorageManager. The detail of it is that "Everything is stored in main memory using a simple vector. No properties are needed to initialize a MemoryStorageManager object. When a MemoryStorageManager instance goes out of scope, all data that it contains are lost."

• The .dat_extension = xxx and .idx_extension = xxx can be used to change the names of those extensions.

[DRKWang]

comment:9

I am trying to build a subclass called rtree_based_realset(), which is inherited from the class RealSet(). But got stuck on the purpose of this subclass. I would like to list my chain of thoughts below and hope to get a verification.

The intention of building this subclass is to keep the interface of this class, (I feel basically the interface can be considered as functions), the same and use rtree data structure to be a preliminary filter for overestimating in order to reduce the computation time. But since rtree only serves as an overestimate from one side, in order to be accurate, this subclass, which we want to define, still needs the finite unit intervals structure inside. In that case, seemingly, the subclass does not have a big difference from the RealSet() class. And the difference between those two classes only comes from whether we use rtree as an auxiliary tools or not.

[mkoeppe]

comment:10

Replying to @DRKWang:

since rtree only serves as an overestimate from one side, in order to be accurate, this subclass, which we want to define, still needs the finite unit intervals structure inside.

That's right, we are not replacing the internal representation of RealSet instances.

We are only adding an index data structure to it.

[DRKWang]

comment:11

Since the rtree.count(), which will return the number of objects whose boxes have intersections with the target point, is not exactly suitable for our expected function, I try to figure out some other functions to replace it. In my opinion, the best function would be rtree.getfirst(), if we should have, which ideally can assist with both RealSet.contains() and RealSet.interesction(). However, I did not get such a function, instead, I found rtree.nearest(), which will return the item whose box has minimal distance to the target point. However, after comparing rtree.nearest() with rtree.count(), the test results shows the rtree.nearest() is slower than rtree.count().

[mkoeppe]

comment:12

As the boxes overestimate the true objects, a count result of 0 is useful.

[mkoeppe]

comment:13

For intersection, can you propose a function for rtree that would help?

[DRKWang]

comment:14

Replying to @mkoeppe:

For intersection, can you propose a function for rtree that would help?

Yes, I think rtree.count() also can help with that when there is a non-intersection between RealSet A and RealSet B, And I have implemented the code.

[DRKWang]

comment:15

I am thinking about changing the internal code of function RealSet.union(). Instead of ordinary function union(A,B), which has two equally operands, the function RealSet.union() has already self, and this format kind of indicates the self set is dominating and the other set will be added to self. In this case, I feel when we are trying to extend this function onto the rtree_bssed_Realset class, it is better to return a rtree_based_RealSet subclass instead of a RealSet class, since the programmer kind of wants to maintain the internal rtree data structure. And since this function only covers the + operation, which correspondingly is also easy to implement with rtree, thus, I feel I can try to achieve that.

[mkoeppe]

comment:16

Replying to @DRKWang:

I feel when we are trying to extend this function onto the rtree_bssed_Realset class, it is better to return a rtree_based_RealSet subclass instead of a RealSet class,

Yes, I agree. This is a general change that should be made throughout the RealSet class. For example, in RealSet.complement, it should probably use self.__class__ instead of RealSet when it constructs the result.

[DRKWang]

comment:17

Yes, I agree. This is a general change that should be made throughout the RealSet class. For example, in RealSet.complement, it should probably use self.__class__ instead of RealSet when it constructs the result.

I see. I have made the changes inside the RealSet class by replacing RealSet with self.__class__. And I feel those replacements only need to happen under 2 following conditions,

• It was after the return. Therefore, the function will finally return the corresponding class, and apart from return, it is necessary to make a replacement.
• We can not make those replacements under the function that has a @staticmethod decorator. Since function with @staticmethod does not have self parameter, which means we can not use self.__class__ to help.

Specifically, the functions I made the change are union(), intersection(), complement(), __mul__().

[mkoeppe]

comment:18

Replying to @DRKWang:

• We can not make those replacements under the function that has a @staticmethod decorator. Since function with @staticmethod does not have self parameter, which means we can not use self.__class__ to help.

Try changing these methods to @classmethod - https://docs.python.org/3/library/functions.html#classmethod

And push the branch please to the ticket

[DRKWang]

Commit: 89595a6

[DRKWang]

Branch: public/32170

[DRKWang]

New commits:

fa4e4ab	build/pkgs/cgal: New
30b487f	build/pkgs/swig:New
df9069a	build/pkgs/cgal_swig_bindings: New
82bf167	add the upstream_url to checksums
40b8a43	Merge tag '9.3.beta6' into t/31098/cgal_swig_bindings
6be25f7	build/pkgs/cgal_swig_bindings/checksums.ini: Fix upstream_url
5695105	build/pkgs/cgal_swig_bindings: Use unpatched cgal-swig-bindings, set Python_EXECUTABLE
e380217	merge th current version to public/32000
89595a6	make a whole change for replacing RealSet() objectas a return with 'self.__class__()' or 'cls()' as a return for 'RealSet' class; add 'Rtree_based_RealSet' class for rtree extension.

[DRKWang]

comment:20

Try changing these methods to @classmethod - https://docs.python.org/3/library/functions.html#classmethod

Got it. I have made the change inside the RealSet object for placing the @staticmethod with @classmethod and use cls() as a wrapper for returning.

And push the branch please to the ticket

The above changes, including whole changes for RealSet, Rtree_based_RealSet class, have been made inside the sage/src/sage/sets/real_set.py file. And the code have been submitted to the branch to the ticket.

[DRKWang]

comment:21

Since I also want to take a test to check whether it performs correctly and since it has been changed under the src folder, for which I thought we can use it explicitly without importing any more. However, sage IDE shows that name 'Rtree_based_RealSet' is not defined.

[DRKWang]

comment:22

Also, I always have a question on how to pushing the latest changes correctly for the sagemath. Assuming there are no previous branches was uploaded for a ticket, If I want to upload a new initial branch for it, either should I merge to the latest version, for example, currently it is sage-9.5, before uploading, or is it ok to keep an old version only if we can run the code without errors?

[mkoeppe]

comment:23

Yes, your current branch has a merge conflict with the current develop, 9.5.beta0, so it is best to merge it in, resolving the conflicts

[DRKWang]

comment:46

Replying to @mkoeppe:

Another way to select a set of projection directions:

• A subset of the set of normal vectors of all facets of all cones.
• Perhaps there is a way to measure how "separating" one of these candidate normal vectors is?

I agree that "the normal vectors of all facets of all cones could also be good candidates for separating these cones". In some sense, those facets can be basically categorized into 2 types, either a facet including an intersection of two maximal cones, or just a pure facet without any intersection. Comparing those two types with each other, the first one seems to be more important and "separating", since the first one focuses on neighbor separating, whereas the second one does not show any information. Thus, we could generate some projecting directions based on the normal vectors of some facets of cones.

Generally, the facets including an intersection can still be categorized into more subtypes based on what is the dimension of the intersection. And heuristically, the higher dimension of the intersection is, the more efficient this normal vector is. But, due to the complexity problem, I feel it is not necessary to categorize those facets at that level.

Even though the above method seems to be efficient, it could be weak when it comes to the case where all cones of a given do not have any intersection except the original points. In that case, there is no "first type of facets" at all. Thus, if that happens, we may need other ways to generate those directions. A good way is to use the direction coming from 2 center points of the tokens of two neighboring cones as supplied directions.

[DRKWang]

comment:48

For the rtree.intersection() function, we can make the return objects of it either a list of boxes with an rtree data structure or a set of boxes without any structure. we can build a flag, like rtree_rebuilt = true/false to indicate the intention of it. Concretely, the boxes with an rtree data structure will be rebuilt by the class of rtree that we defined.

[DRKWang]

comment:49

For the static construction, we may build the rtree by 3 ways, top-down, bottom-up or hybrid. Supposedly, there are $n$ boxes in $k$ dimension for construction.

The idea of top-down should be made by recursive operation, which means the parent tree will be built after their children's trees were built. My expectation will be something like an H tree (https://en.wikipedia.org/wiki/H_tree), which has a property of self-similar fractal, or like a $k-D tree$. Anyway, the similarity of the subsets of boxes can be displayed on the rtree. Concretely, one way to do this is as follows:

1. Separate the whole boxes into 2 subsets, equivalently in terms of numbers, based on their X-coordinate. And take those 2 subsets as the materials for building the left sub-rtree and right sub-rtree.
2. Cyclically switch the criteria from x-coordinate to y-coordinate, to z-coordinate with separating the subsets smaller and smaller until the number of subsets equal 1. For this case, the complexity of it will be $nlog(n)+2* n/2 log(n/2) + 4* n/2 log(n/4) + ... + n*1 log(1) <= n log(n) * log(n) $

For the bottom-up, we may use the greedy strategy to do the local optimization on each step. We can always merge two boxes based on some criteria, like merging the two boxes that have minimal 1-norm distance. And replace the 2 boxes with merged box, so that the size of the input will be decreased by 1, (now the number of boxes will be $n-1$). Thus, after $n-1$ steps, the rtree will be built completely.

The complexity will be at least n^2, since we will compute the distance at least.

For the hybrid one, we can build the rtree partially by top-down and partially by bottom-up in certain steps. But, how to combining them adaptively is still a hard problem for me. I would like to implement the first two before thinking about the third one.

[DRKWang]

comment:50

The following content in latex form can be checked in this link.(https://github.com/DRKWang/cutgeneratingfunctions_new/blob/main/logsfile/log48.ipynb)

To determine whether a given fan is a pointed fan or not, we may encode it into an linear optimization.

Observation:

If a fan is pointed, then there exists a containing-the-original-point hyperplane that could divide the space into 2 half-spaces, one containing the fan and the other not.

Analysis:

We may think about finding the normal vector of this hyperplane via optimzation.
Let us encode it.

Supposing that the normal vector of this hyperplane is $x$ and 1-norm vectors of the generating rays of all cones in the fan are $r_i, i\in A $,

then we may have let the constraints be:

$x \cdot r_i \ge 0$

Since we do not hope to get $x = 0$, we may add an extra constraint, $x_1 = 1$, for getting one representing element.

If we can not find any solution for this system, then the fan is not pointed.
If we can find one solution, then the fan is pointed.

In practice, in order to give a initial solution for simplex method, we may change the system into the following form:

max. : $x_1$

s.t.

$x \cdot r_i \ge 0$

$1\ge x_1 \ge 0$

As is desired, $x = 0$ always be the solution for this system.

A better hyperplane for generating better projection

Seemingly, we may have optional hyperplanes for the same system. The ideal hyperplane would be constraints are all $>$, instead of $\ge$. We may get rid of $x_1 = 1$ and add a little "balance" for each constraints to get it. i.e.

$x \cdot r_i - b_i\ge 0$

$b_i$ are all some given small positive numbers.

Similarly, in order to have a initial solution, we may change the system into the following form:

max. :$ b_0$

$x \cdot r_i - b_i\ge 0$

$b_i \ge b_0$

$1\ge b_0 \ge 0$

As is desired, $x = 0$ always be the solution for this system.

[DRKWang]

comment:51

Except for the above view, we may think differently from the view of the convex hull of the fan, but it will be much harder than the above one, in terms of implementation and theory.

Based on the following observation, we can have:

If the convex hull of the fan is the universal space, then it is not pointed; otherwise, it is pointed.

Thus, we can figure out the convex hull to determine whether it is pointed or not.

In the beginning, my imagination to compute the convex hull will be the "Graham scan method", a popular method for computing the convex hull in 2 dimensions. But when it comes to high dimensions, this algorithm will fail and be replaced by an algorithm called "Quickhull", (https://en.wikipedia.org/wiki/Quickhull), for computing the convex hull. This algorithm is a little randomized for initial selection and a little bit complicated understanding. And in terms of performance, although this algorithm may recursively get the convex hull, we still need to compute the normal vector for the hyperplane, onto which the fan will be projected, based on the above optimization.

[DRKWang]

comment:52

As is enlightened by your idea for considering the double descriptions of cones, facets, and generating rays, we may have a facets-based method to find the proper normal vector of the separating hyperplane for a pointed fan.

Observation:

For a cone, a ray is in the (interior or boundary) of its polar cone, if and only if, the ray can induce (an interior or a boundary) separating hyperplane, whose normal vector is the ray, for this cone.

Note that the polar cone of cone $C$ is defined by the negative dual cone of cone $C$, i.e.

the polar cone of cone $C$ = ${y | y\cdot x \le 0, \forall x \in C}$

Thus, to figure out the separating hyperplane for a fan, the polar cones of cones in this fan can provide help.

Observation:

For a fan, a ray is in the interior points of the common intersection of all polar cones of the fan, if and only if, the ray can induce an interior separating hyperplane for the fan.

Proof

A ray is in the interior points of the common intersection of all polar cones,

$\iff$ it must be in the interior points of each polar cone,

$\iff$ it must induce the same interior separating hyperplane for each cone in the fan,

$\iff$ it must induce an interior separating hyperplane for the fan.

Thus, we can have a better understanding of the location of the normal vector of separating hyperplane for a fan. Meanwhile, we can also use it to determine whether a fan is a pointed fan or not.

[DRKWang]

comment:53

Based on these testing records https://github.com/DRKWang/cutgeneratingfunctions_new/blob/main/test_for_facets_and_generating_rays_based_methods.ipynb, in terms of time, the generating-rays-based method has significantly outperformed the facets-based method.
However, according to test 4 in the above records, the LP-solver may have some numerical issues and it may lead to some inaccurate determination.

In my two cents, since LP-solver only takes a little time and the solution for the linear programming can be easily verified, we may use the LP-solver first to try to get a solution, if it fails or it only gives zero solution, then we use the facets-based method instead; otherwise, if the solution is non-zero, we may verify its correctness, and if the solution is correct, then we are home and we do not need the facets-based method. This mixed-method is better than the pure facets-based method especially in the case that the fan is pointed because in that case, it is really difficult to compute common intersection for all cones with keeping in analytical status according to the testing records.

[DRKWang]

comment:54

Apart from the function contains(), I just realized that the class RationalPolyhedralFan also has another function support_contains(self,*args), which is used for checking if a point is contained in the support of the fan, and the support of a fan is the union of all cones of the fan. This function is pretty good for the usage of rtree, and based on the following test results, https://github.com/DRKWang/cutgeneratingfunctions_new/blob/main/fan_rtree_test.ipynb, it shows that the support of rtree gives better performance than the original one. The only disadvantage is that the construction of rtree in the initial part may take a little longer time. But based on the test results, it only takes extra at most 1/3 time compared with the original one, which is much smaller than saving time for using support_contains(self,*args).

[sagetrac-git]

Changed commit from d160842 to ebe6d92

[sagetrac-git]

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

c9b20a0	fixed the rtree's checksum
078f823	adding self.polar() for class ConvexRationalPolyhedralCone()
ebe6d92	add subclass: RationalPolyhedralFan_rtree for function support_contains().

[sagetrac-git]

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

a2fa0ba

Fixed bugs in docstring.

[sagetrac-git]

Changed commit from ebe6d92 to a2fa0ba

[mkoeppe]

Author: Binshuai Wang

[sagetrac-git]

Changed commit from a2fa0ba to 05ea99d

[sagetrac-git]

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

05ea99d

Fixed the failed examples in Class.

[mkoeppe]

comment:60

Docbuild fails

OSError: /sage/local/var/lib/sage/venv-python3.10.2/lib/python3.10/site-packages/sage/geometry/cone.py:docstring of sage.geometry.cone.ConvexRationalPolyhedralCone.polar:6: WARNING: Bullet list ends without a blank line; unexpected unindent.

[mkoeppe]

Dependencies: #34277

[mkoeppe]

Work Issues: Rebase on #34277