Add MOC vs. other geometries casts

.gitignore

@@ -8,3 +8,4 @@
 regression.out
 regression.diffs
 tags
+buildpod

HACKING

@@ -0,0 +1,50 @@
+SQL definitions
+===============
+
+Long version:
+<https://www.postgresql.org/docs/current/extend-extensions.html>.
+
+If you're writing new features that require SQL support, pick some
+descriptive name; let's say my_new_op.  
+
+Put your new code into a file called pgs_my_new_op.sql.in.  The .in
+extension here usually indicates "it's for copying stuff together";
+usally, not much processing is done on such files.
+
+Then edit the Makefile.  The PGS_SQL variable contains a list of the
+SQL files eventually copied together, without the .in.  Add your new
+file there.
+
+You will also need to create an upgrade file.  In order to tell postgres
+to execute it, increase PGSPHERE_VERSION as appropriate.  As a
+consequence, you will have to::
+
+	git mv pg_sphere--<old version>.sql.in pg_sphere--<new version>.sql.in
+
+and also to update the version in pg_sphere.control.
+
+Then create a make rule:: 
+
+	pg_sphere--<old version>--<new version>.sql: pgs_my_new_op.sql.in
+		cat $^ > $@
+
+(of course, this will extend to having multiple sql.in files).
+
+Finally, add the target of that rule to the DATA_built variable.
+
+
+Regression tests
+================
+
+Regressions tests are as per
+<https://www.postgresql.org/docs/current/extend-pgxs.html>.
+
+In short, write queries executing your new features into a file
+sql/my_new_op.sql, and add "my_new_op" (without the extension or the
+directory name) to both REGRESS and TESTS in the Makefile.
+
+Then touch expected/my_new_op.out, run make test.  This will of course
+fail, because your tests hopefully will output something.  But then you
+can pick out the diff from
+/var/lib/postgresql/pgsphere/regression.diffs, have another critical
+look at it and generatoe your .out file from it.

Makefile

@@ -1,4 +1,4 @@
-PGSPHERE_VERSION = 1.2.0
+PGSPHERE_VERSION = 1.2.1
 
 # the base dir name may be changed depending on git clone command
 SRC_DIR = $(shell basename $(shell pwd))
@@ -17,18 +17,19 @@ DATA_built  = $(RELEASE_SQL) \
 			  pg_sphere--1.0_gavo--1.1.5beta0gavo.sql \
 			  pg_sphere--1.1.5beta0gavo--1.1.5beta2gavo.sql \
 			  pg_sphere--1.1.5beta2gavo--1.1.5beta4gavo.sql \
-			  pg_sphere--1.1.5beta4gavo--1.2.0.sql
+			  pg_sphere--1.1.5beta4gavo--1.2.0.sql \
+			  pg_sphere--1.2.0--1.2.1.sql
 
 DOCS        = README.pg_sphere COPYRIGHT.pg_sphere
 REGRESS     = init tables points euler circle line ellipse poly path box index \
 			  contains_ops contains_ops_compat bounding_box_gist gnomo healpix \
-			  moc
+			  moc mocautocast
 
 REGRESS_9_5 = index_9.5 # experimental for spoint3
 
 TESTS       = init_test tables points euler circle line ellipse poly path box index \
 			  contains_ops contains_ops_compat bounding_box_gist gnomo healpix \
-			  moc
+			  moc mocautocast
 
 ifndef CXXFLAGS
 # no support for CXXFLAGS in PGXS before v11
@@ -45,7 +46,8 @@ PGS_SQL     = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \
    pgs_line.sql pgs_ellipse.sql pgs_polygon.sql pgs_path.sql \
    pgs_box.sql pgs_contains_ops.sql pgs_contains_ops_compat.sql \
    pgs_gist.sql gnomo.sql \
-   healpix.sql pgs_gist_spoint3.sql pgs_moc_type.sql pgs_moc_compat.sql pgs_moc_ops.sql
+   healpix.sql pgs_gist_spoint3.sql pgs_moc_type.sql pgs_moc_compat.sql pgs_moc_ops.sql \
+   pgs_moc_geo_casts.sql
 PGS_SQL_9_5 = pgs_9.5.sql # experimental for spoint3
 
 USE_PGXS = 1
@@ -155,7 +157,7 @@ ifeq ($(pg_version_9_5_plus),y)
 else
 endif
 
-# local stuff follows here, next will be "beta2"
+# local stuff follows here
 
 AUGMENT_GAVO_111 = $(AUGMENT_UNP_111) healpix.sql # for vanilla 1.1.1 users
 UPGRADE_GAVO_111 = $(UPGRADE_UNP_COMMON)
@@ -197,6 +199,12 @@ ifeq ($(has_parallel), n)
 	sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL
 endif
 
+pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in
+	cat $^ > $@
+ifeq ($(has_parallel), n)
+	sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL
+endif
+
 # end of local stuff
 
 sscan.o : sparse.c

doc/functions.sgm

@@ -638,6 +638,34 @@
             </programlisting>
           </example>
 
+          <funcsynopsis>
+            <funcprototype>
+              <funcdef><function>max_order</function></funcdef>
+                <paramdef>smoc</paramdef>
+            </funcprototype>
+          </funcsynopsis>
+          <para>
+            Returns the maximal order of an smoc.
+          </para>
+          <para>
+          	The maximal order of a MOC can be higher than than the highest
+          	order appearing in the serialisation.  For instance, in the example
+          	below the full-sky coverage given at order 6 means that no patch
+          	larger than about 1 degree is not covered.  In the ASCII
+          	serialisation, give the maximum order with an with an empty cell
+          	list.
+          </para>
+
+          <example>
+            <title>Obtaining a MOC order</title>
+            <programlisting>
+<![CDATA[SELECT max_order(smoc('0/0-11 6/')) as order;]]>
+<![CDATA[    order]]>
+<![CDATA[-----------]]>
+<![CDATA[ 6 ]]>
+            </programlisting>
+          </example>
+
         </sect2>
 
     </sect1>

doc/gen_moccast.py

@@ -0,0 +1,156 @@
+# A script to create the automatic casts for overlaps and intersects
+# between MOCs and spolys/scircles.
+#
+# This has originally been used to create pg_sphere--1.2.0--1.2.1.sql.
+# Before 1.2.1 is release, this can be fixed to improve that SQL.
+# After the 1.2.1 release, this is just documentation on how MOC
+# casts were generated that is perhaps just a bit more readable than
+# that bunch of SQL.
+
+import datetime
+import re
+import sys
+
+
+OVERLAP_DEFS = [
+    # func_stem, operator, commutator
+    ('subset', '<@', '@>'),
+    ('not_subset', '!<@', '!@>'),
+    ('superset', '@>', '<@'),
+    ('not_superset', '!@>', '!<@'),
+]
+
+INTERSECT_DEFS = [
+    # func_stem, operator, commutator
+    ('intersect', '&&', '&&'),
+    ('not_intersect', '!&&', '!&&'),
+]
+
+
+GEO_TYPES = ["scircle", "spoly"]
+
+OP_DEFS = OVERLAP_DEFS
+
+
+class Accum:
+    """an accumulator for our output.
+    """
+    def __init__(self):
+        self.parts = []
+
+    @property
+    def content(self):
+        return "".join(self.parts)
+
+    def write(self, s):
+        self.parts.append(s)
+
+    def writeln(self, *strings):
+        self.parts.append("\n".join(strings)+"\n")
+
+    def replace_last(self, subs):
+        """replaces the last non-whitespace char with the string subs.
+        """
+        for index, part in enumerate(reversed(self.parts)):
+            if part.strip():
+                break
+        else:
+            # nothing to replace
+            return
+
+        index = -1-index
+        self.parts[index] = re.sub("[^\s](\s*)$", 
+            lambda mat: subs+mat.group(1),
+            self.parts[index])
+
+    def introduce_section(self, sec_name):
+        self.writeln()
+        self.writeln("-- #################################")
+        self.writeln(f"-- {sec_name}")
+
+
+def emit_drop_code(accum):
+    accum.introduce_section("Cleanup")
+
+    accum.writeln("DROP OPERATOR IF EXISTS")
+    for _, op, _ in OP_DEFS:
+        for geo_type in GEO_TYPES:
+            accum.writeln(f"    {op} (smoc, {geo_type}),")
+            accum.writeln(f"    {op} ({geo_type}, smoc),")
+    accum.replace_last(";")
+
+
+def make_negator(op):
+    if op.startswith("!"):
+        return op[1:]
+    else:
+        return "!"+op
+
+
+def emit_op_def(accum, operator, leftarg, rightarg, procedure, commutator):
+    accum.writeln(
+        f"CREATE OPERATOR {operator} (",
+        f"  LEFTARG    = {leftarg},",
+        f"  RIGHTARG   = {rightarg},",
+        f"  PROCEDURE  = {procedure},",
+        f"  COMMUTATOR = '{commutator}',",
+        f"  NEGATOR    = '{make_negator(operator)}',",
+        f"  RESTRICT   = contsel,",
+        f"  JOIN       = contjoinsel",
+        f");")
+
+
+def emit_op_and_func(accum, op_def):
+    func_stem, operator, commutator = op_def
+    for geo_type in GEO_TYPES:
+        func_name = f"{geo_type}_{func_stem}_smoc"
+        accum.writeln(
+            f"CREATE OR REPLACE FUNCTION {func_name}(",
+            f"  geo_arg {geo_type}, a_moc smoc) RETURNS BOOL AS $body$",
+            f"    SELECT smoc(max_order(a_moc), geo_arg) {operator} a_moc",
+            f"  $body$ LANGUAGE SQL STABLE;")
+        emit_op_def(accum, operator,
+            geo_type, "smoc",
+            func_name,
+            commutator)
+
+        accum.writeln()
+
+        func_name = f"smoc_{func_stem}_{geo_type}"
+        accum.writeln(
+            f"CREATE OR REPLACE FUNCTION {func_name}(",
+            f"  a_moc smoc, geo_arg {geo_type}) RETURNS BOOL AS $body$",
+            f"    SELECT a_moc {operator} smoc(max_order(a_moc), geo_arg)",
+            f"  $body$ LANGUAGE SQL STABLE;")
+        emit_op_def(accum, operator,
+            "smoc", geo_type,
+            func_name,
+            commutator)
+
+        accum.writeln()
+
+
+def main():
+    accum = Accum()
+
+    accum.writeln("-- MOC/geometry automatic casts.")
+    accum.writeln(f"-- Generated {datetime.date.today()} by {sys.argv[0]}.")
+    accum.writeln(f"-- Re-generation needs to be triggered manually.")
+    accum.writeln()
+    emit_drop_code(accum)
+
+    accum.introduce_section(" smoc/geo OVERLAPS")
+    for op_def in OVERLAP_DEFS:
+        emit_op_and_func(accum, op_def)
+        accum.writeln()
+
+    accum.introduce_section(" smoc/geo INTERSECTS")
+    for op_def in INTERSECT_DEFS:
+        emit_op_and_func(accum, op_def)
+        accum.writeln()
+
+    print(accum.content)
+
+
+if __name__=="__main__":
+    main()

doc/operators.sgm

@@ -231,10 +231,19 @@
             <para>
                 An overlap or contain operator does not exist for all
                 combinations of data types.
-                For instance, scircle <@ <type>spoint</type> is
+                For instance, scircle &lt;@ <type>spoint</type> is
                 useless because a spherical point can never contain a
                 spherical circle.
             </para>
+            <para>
+                When one of the arguments of such an operator is a MOC
+                and the other is an scircle or an spoly, the non-MOC argument
+                is converted to a MOC of the order of the maximum order of
+                the MOC.  When comparing against a MOC-valued column, it
+                is usually much faster to explicitly convert the geometry
+                using the smoc constructor, as the conversion will then
+                only happen once.
+            </para>
             <example>
               <title>Is the left circle contained by the right circle?</title>
               <programlisting>
@@ -255,6 +264,40 @@
 <![CDATA[(1 row)]]>
               </programlisting>
             </example>
+
+            <example>
+              <title>Overlaps between a circle and a moc</title>
+              <programlisting>
+<![CDATA[sql> SELECT scircle '<(37d, 5d), 0.25d>' <@ smoc('4/1117') AS test ;]]>
+<![CDATA[ test]]>
+<![CDATA[------]]>
+<![CDATA[ f]]>
+<![CDATA[(1 row)]]>
+              </programlisting>
+            </example>
+
+            <example>
+              <title>Overlaps between a circle and a moc with explicit order</title>
+              <programlisting>
+<![CDATA[sql> SELECT scircle '<(37d, 5d), 0.25d>' <@ smoc('4/1117 5/') AS test ;]]>
+<![CDATA[ test]]>
+<![CDATA[------]]>
+<![CDATA[ t]]>
+<![CDATA[(1 row)]]>
+              </programlisting>
+            </example>
+
+            <example>
+              <title>Overlaps between a circle and a moc with explicit cast (normally faster)</title>
+              <programlisting>
+<![CDATA[sql> SELECT smoc(5, scircle '<(37d, 5d), 0.25d>') <@ smoc('4/1117 5/') AS test ;]]>
+<![CDATA[ test]]>
+<![CDATA[------]]>
+<![CDATA[ t]]>
+<![CDATA[(1 row)]]>
+              </programlisting>
+            </example>
+
         </sect2>
 
         <sect2 id="op.cross">

expected/init_test.out

@@ -33,5 +33,5 @@ psql:pg_sphere.test.sql:8569: NOTICE:  argument type pointkey is only a shell
 psql:pg_sphere.test.sql:8575: NOTICE:  argument type pointkey is only a shell
 psql:pg_sphere.test.sql:8581: NOTICE:  argument type pointkey is only a shell
 psql:pg_sphere.test.sql:8587: NOTICE:  argument type pointkey is only a shell
-psql:pg_sphere.test.sql:9152: NOTICE:  return type smoc is only a shell
-psql:pg_sphere.test.sql:9158: NOTICE:  argument type smoc is only a shell
+psql:pg_sphere.test.sql:9154: NOTICE:  return type smoc is only a shell
+psql:pg_sphere.test.sql:9160: NOTICE:  argument type smoc is only a shell