To replicate, here is an example of a bad surface. It just has a single twist of 0.6.

import superficial._
val egBad = SkewPantsSurface(
  4,
  Set(
    SkewCurve(PantsBoundary(0, Z3(0)), PantsBoundary(2, Z3(1)), 0.0, 1.0),
    SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.0, 1.0),
    SkewCurve(PantsBoundary(0, Z3(2)), PantsBoundary(3, Z3(0)), 0.0, 1.0),
    SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(3, Z3(1)), 0.0, 1.0),
    SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0),
    SkewCurve(PantsBoundary(2, Z3(2)), PantsBoundary(3, Z3(2)), 0.0, 1.0)
  )
)

We see the extra curves, i.e., those that are not in the boundary of a face (show is ammonite repl syntax):

@ show(egBad.edges.filter(e => egBad.edgeOccurences(e) == 0)) 
HashSet(
  PantsSeam(
    1,
    SkewCurveVertex(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 1.1),
    SkewCurveVertex(SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.0, 1.0), 0.5),
    false
  ),
  PantsSeam(
    1,
    SkewCurveVertex(SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.0, 1.0), 0.5),
    SkewCurveVertex(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.10000000000000009),
    true
  ),
  SkewCurveEdge(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.10000000000000009, false),
  PantsSeam(
    1,
    SkewCurveVertex(SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(3, Z3(1)), 0.0, 1.0), 0.0),
    SkewCurveVertex(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.10000000000000009),
    false
  ),
  PantsSeam(
    1,
    SkewCurveVertex(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.6),
    SkewCurveVertex(SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(3, Z3(1)), 0.0, 1.0), 0.0),
    true
  ),
  SkewCurveEdge(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.0, false),
  SkewCurveEdge(SkewCurve(PantsBoundary(0, Z3(1)), PantsBoundary(1, Z3(0)), 0.6, 1.0), 0.6000000000000001, false)
)

Warning: I am likely to work on this later this week.

To construct surfaces efficiently, one should be able to give data of the form:

• vertices by name
• positively oriented edges by name and initial and terminal vertices, with the negative ones automatic.
• faces by name and boundary.

We can even use the sourcecode for some syntactic sugar.

All enumerations should be done with views or lazy lists, or sometimes even iterators.

When enumerating, we get an error because a PLArc does not satisfy its initialization requirement. This is a new error presumably introduced by recent corrections.

• Code to generate error:

import superficial._
import PantsSurface._
import scala.language.postfixOps
val surf = PantsSurface.allClosed(2)(1)
val skewsurf = SkewPantsSurface.enumerate(surf, Vector(0.2), Vector(1, 2))(1)
val pathsuptolen5 = NormalPath.enumerate[SkewPantsHexagon](skewsurf, Some(5))
val clpathsuptolen5 = pathsuptolen5.filter(p => p.isClosed)
val uniqclpathsuptolen5 = NormalPath.uniqueUptoFlipAndCyclicPerm(clpathsuptolen5)
val plpathsuptolen5 = PLPath.enumMinimalClosedFamily(uniqclpathsuptolen5.values.toVector, 0.05, 1.3)

The error seen:

java.lang.IllegalArgumentException: requirement failed
  scala.Predef$.require(Predef.scala:325)
  superficial.PLPath.<init>(PLArc.scala:261)
  superficial.PLPath$.$anonfun$appendPLArc$6(PLArc.scala:392)
  scala.collection.Iterator$$anon$9.next(Iterator.scala:575)
  scala.collection.mutable.Growable.addAll(Growable.scala:65)
  scala.collection.mutable.Growable.addAll$(Growable.scala:60)
  scala.collection.immutable.VectorBuilder.addAll(Vector.scala:1565)
  scala.collection.immutable.Vector$.from(Vector.scala:1349)
  scala.collection.immutable.Vector$.from(Vector.scala:34)
  scala.collection.SeqFactory$Delegate.from(Factory.scala:306)
  scala.collection.immutable.IndexedSeq$.from(Seq.scala:115)
  scala.collection.immutable.IndexedSeq$.from(Seq.scala:112)
  scala.collection.IterableOps$WithFilter.map(Iterable.scala:884)
  superficial.PLPath$.$anonfun$appendPLArc$4(PLArc.scala:380)
  scala.collection.parallel.AugmentedIterableIterator.flatmap2combiner(RemainsIterator.scala:123)
  scala.collection.parallel.AugmentedIterableIterator.flatmap2combiner$(RemainsIterator.scala:121)
  scala.collection.parallel.immutable.ParHashSet$ParHashSetIterator.flatmap2combiner(ParHashSet.scala:77)
  scala.collection.parallel.ParIterableLike$FlatMap.leaf(ParIterableLike.scala:1024)
  scala.collection.parallel.Task.$anonfun$tryLeaf$1(Tasks.scala:52)
  scala.runtime.java8.JFunction0$mcV$sp.apply(JFunction0$mcV$sp.scala:18)
  scala.util.control.Breaks$$anon$1.catchBreak(Breaks.scala:97)
  scala.collection.parallel.Task.tryLeaf(Tasks.scala:55)
  scala.collection.parallel.Task.tryLeaf$(Tasks.scala:49)
  scala.collection.parallel.ParIterableLike$FlatMap.tryLeaf(ParIterableLike.scala:1020)
  scala.collection.parallel.AdaptiveWorkStealingTasks$AWSTWrappedTask.internal(Tasks.scala:159)
  scala.collection.parallel.AdaptiveWorkStealingTasks$AWSTWrappedTask.internal$(Tasks.scala:156)
  scala.collection.parallel.AdaptiveWorkStealingForkJoinTasks$AWSFJTWrappedTask.internal(Tasks.scala:304)
  scala.collection.parallel.AdaptiveWorkStealingTasks$AWSTWrappedTask.compute(Tasks.scala:149)
  scala.collection.parallel.AdaptiveWorkStealingTasks$AWSTWrappedTask.compute$(Tasks.scala:148)
  scala.collection.parallel.AdaptiveWorkStealingForkJoinTasks$AWSFJTWrappedTask.compute(Tasks.scala:304)
  java.util.concurrent.RecursiveAction.exec(RecursiveAction.java:189)
  java.util.concurrent.ForkJoinTask.doExec(ForkJoinTask.java:290)
  java.util.concurrent.ForkJoinPool$WorkQueue.topLevelExec(ForkJoinPool.java:1020)
  java.util.concurrent.ForkJoinPool.scan(ForkJoinPool.java:1656)
  java.util.concurrent.ForkJoinPool.runWorker(ForkJoinPool.java:1594)
  java.util.concurrent.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:183)

• Hyperbolic surfaces are given by pairs of pants glued with twists, with pairs of pants given by right-angled hexagons.
• For each of these, we determine if there is an approximate systole with a large number of curves, where the systole is the set of geodesics of minimal length.
• Goal is to try to determine, by compact enumeration, the largest systole for a surface of genus 3.

In the following path,
val path1 = PLPath( NormalPath( Vector( NormalArc( 7, 5, SkewPantsHexagon( 1, true, Set( SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(2, Z3(2)), PantsBoundary(3, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(2, Z3(1)), 0.2, 2), SkewCurve(PantsBoundary(0, Z3(2)), PantsBoundary(1, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(3, Z3(1)), PantsBoundary(3, Z3(2)), 0.2, 1), SkewCurve(PantsBoundary(0, Z3(0)), PantsBoundary(0, Z3(1)), 0.2, 1) ) ) ), NormalArc( 5, 3, SkewPantsHexagon( 1, false, Set( SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(2, Z3(2)), PantsBoundary(3, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(2, Z3(1)), 0.2, 2), SkewCurve(PantsBoundary(0, Z3(2)), PantsBoundary(1, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(3, Z3(1)), PantsBoundary(3, Z3(2)), 0.2, 1), SkewCurve(PantsBoundary(0, Z3(0)), PantsBoundary(0, Z3(1)), 0.2, 1) ) ) ), NormalArc( 0, 1, SkewPantsHexagon( 2, true, Set( SkewCurve(PantsBoundary(1, Z3(2)), PantsBoundary(2, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(2, Z3(2)), PantsBoundary(3, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(1, Z3(1)), PantsBoundary(2, Z3(1)), 0.2, 2), SkewCurve(PantsBoundary(0, Z3(2)), PantsBoundary(1, Z3(0)), 0.2, 1), SkewCurve(PantsBoundary(3, Z3(1)), PantsBoundary(3, Z3(2)), 0.2, 1), SkewCurve(PantsBoundary(0, Z3(0)), PantsBoundary(0, Z3(1)), 0.2, 1) ) ) ) ) ), Vector(0.0, 0.07541963176135802, 0.2), Vector(2.15, 0.0, 0.2) )
the third arc has zero length. However, while plotting the arc in blue, we get three arcs of non zero length -

the first is across edges 20 and 21, the second is nearly parallel to edge 2, and the third runs halfway along edge 14, to continue nowhere. The third arc should just be a point, so that the path starts at a vertex of edge 21 and comes back to it.

The examples should be in the main code, and tests should be added.

• Define Boolean functions on two-complexes determining whether they are surfaces, closed surfaces and punctured surfaces.
• Also decide if a surface is connected.

The problem

• While serialization as json with upickle is clean and easy, the generated files are very large - huge in the case of the map from paths to unique classes.
• Most of the file is field names.
• This means that when upack is used to binary serialize we get a smaller file, but only by a factor of 2/3.
• Another source of inefficiency is in a normal path, the base surface is stored separately for each normal arc. Actually it should be stored once for the entire collection of PL-Arcs

Solution

• For each of the case classes, have an associated compact type made of tuples (of Integers, BigDecimals, Doubles and Booleans) with two way conversions.
• A clean design is to have
  • a statement in each companion object of the form type Cpt = (Int, BigDecimal, Boolean) or equivalent.
  • a method in the companion object fromCpt(cpt: Cpt): SkewCurve or equivalent.
  • a method in the case class mapping to the corresponding compact types.
  • for the collections and maps we want to store associated to a surface, a special associated compact type avoiding duplication of the associated surface.
• In addition we should have round-trip tests to see that serializing and deserializing gives the same object.
• As these are built from standard scala types, they can be automatically serialized by upickle (with upack for efficiency) or something else (e.g. Boopickle).

This is a summary as of June 30, 2021, commit 8a4d930.

Non-normal homotopy/isotopy

The subtle part of the code, which needs careful verification, is homotopies of normal curves that change the underlying normal curve. There are two forms of errors:

1. Missed homotopy across vertices.
2. Wrong homotopy across vertices.

Missed homotopy

• This is the less serious error, as it can be checked post-facto, i.e., after we have computed the systole, to ensure that all the curves are indeed homotopically non-trivial.
• For this, we have to account for all curves in a neighbourhood - these should be part of the systole but recognized as equivalent or be too long.
• We should also check that the underlying normal curves are simple and connected.

Wrong homotopy

• This can happen anywhere along the way, so we should guard against this.
• The best way seems to be to generate objects that record details of homotopies, and checking later that:
  • these are correct.
  • there is a homotopy starting with every PL curve that is minimal in its normal class which is not part of the systole.
• One can have a single sealed trait, with the initial curve as a method (and perhaps a check function), and with case classes for the different homotopies.
• Homotopies as used are stored as side-effects in a buffer to be verified later.
• Verification of homotopies can be done with:
  • using functions that check validity (implemented separately from creating the homotopy, so easier to check their correctness).
  • plotting the start and end curves and inspecting correctness.
  • the very act of writing the code to generate these will lead to scrutinizing curves where there are homotopies.

@shabarishch Can you check:

• if the summary looks fine and if I am missing something
• what are the functions where we have a homotopy, and what is the nature of the homotopy (including null-homotopy by returning None)

• The function should work if the two-complex given is a surface.
• If it is a 1-vertex, 1-face quadrangulation of a negative Euler characteristic closed surface, then should give non-positive quadrangulation.

This is to keep track of changes as we figure they need to be made while chatting.

The loop "a1!" in the quad for the ClosedSurface(2) does not canonicise correctly.

The function shorten moves a PLPath to the local minimum for length in its isotopy class. It now takes in two additional parameters -

1. uniqrepuptoflipandcyclicper - a map which gives, for each normal path, its unique representative upto flips and cyclic permutations that has been enumerated into a PLPath
2. enumdata - a map defined on the set of normalpaths in the image of uniqrepuptoflipandcyclicper that takes a normalpath to its minimal length plpath.

For shortening a PLPath of base length n, we need enumdata for paths of base length upto (n+1) since shifting across a vertex may increase base length by 1. The attached worksheet presents an error on the second path.
len5paths.txt

Vertex indices, not vertices should be used. Here is an example giving the error.

@ import superficial._ 
import superficial._

@ val genus2 = new StandardSurface(2) 
genus2: StandardSurface = superficial.StandardSurface@3bf4b018

@ val q2 = genus2.quadrangulate 
q2: TwoComplex = superficial.TwoComplex$quad$1$@64641998

@ q2.edges 
res17: Set[Edge] = Set(superficial.EdgePair@2e07dc23.positive, superficial.EdgePair@2e07dc23.negative)

This should be fixed @anotherArka

• We can collapse edges with endpoints distinct, and merge distinct faces to get a new two-complex for the same surface.
• By doing this, if a surface is connected, then repeat this to get a surface with a single edge and a single vertex.
• We should check that the complex is a surface and connected first, using #5.

The problems

• Consider two curves C and D with indices i and j so that the vertices of the curves at these indices coincide.
• Throughout, before and after are to be viewed cyclically.
• Let c0 and c1 be the edges of C before and after the coinciding vertex and let d0 and d1 be similar.
• Assume that both the curves are reduced, and both are not positive powers of the same curve (in the latter case we replace one by its inverse for computing intersection numbers).
• We have to decide whether the curve D comes in from above C at (or before) the given intersection. Note that we are not ruling out c0 == d0.
• Thus we have to define a boolean function fromAbove(C, D, i, j) of the curves and indices.
• We will define this recursively.

The function fromAbove

1. If c0 == d0 then we replace indices by i - 1 and j - 1 (cyclically) and compute. The hypothesis guarantees that this is not an infinite loop.
2. if d0 = inverse(c1), then we compute fromAbove(C, D, i, j) = !fromAbove(C, inverse(D), i, -j). This immediately gives the above case, so we do simplify.As the curves are reduced, we can only get this case if we have a single intersection point, and once we flip this way we will not have to flip again for the same intersection.
3. Finally, if neither of the above happen then rotate counterclockwise from inverse(c1) to c0 and return true if and only if we encounter d0 along the way.

• Given a Boolean condition such as being reduced, which is inherited by subpaths, we should enumerate all edge paths satisfying that condition.
• This should recursively give lazy lists.

• Map edge paths in a non-positive quadrangulation to their canonical geodesics.
• For pairs of such canonical geodesics, determine which of the intersection intervals are essential and determine the sign of these, to calculate the Geometric intersection number.

• For this we should have a basePoint for polygons.
• This is anyway correct in terms of equality.

• Define when a two-complex is a subdivision of another.
• Have also maps between edge-paths of the subdivided complex and the original one.
• For modifications such as 1-vertex 1-face in #6 define a common subdivision.

• This is needed even before starting #9.
• It is best to allow degenerate geodesics.
• This will mean having an initial vertex and a vector/list as fields for defining, with a require statement verifying that this is a path, and deducing terminal vertex.
• May also want to check that it is in a given 2-complex.

Theory

• As we move away from the point on a geodesic l that minimizes the total distance to two points (i.e., the intersection of geodesics) we expect that the distance increases monotonically.
• This should be straightforward to prove using the first variational formula for energy, since the change with the angle for piecewise geodesic curves is of first order and independent of curvature.

Algorithms.

Varying endpoints

• Instead of checking for all points in a range, we start with an initial point, check this and its neighbours.
• Recursively, we have the data that we have moved in a specific direction. So we only have to check moving further in that direction.
• To combine uniformly with the initial step, we can have an "increment set" with 1 or 2 points - we even get 0 points if we go out of the range.
• For a recursive step, we check whether we have a lower value for any dispacement in the increment set.
• To avoid recomputing, we can optionally call with precomputed value at the point.
• The gain is large if we start close to the minimal point, which we can do by taking the same points as a nearby surface.

Varying paths

• We proceed recursively on the number of arcs as before.
• For closed paths, we should ensure that ends are matched.

Have better printing and precision.

• In the cases where we collapse edges, merge faces or generate quadrangulations, we should also define (and return from functions) functions between edge-paths of the complexes.
• Each edge-path should be taken to a homotopic one.

• Determine the action of a Dehn twist on canonical geodesics (using #9).
• Using this, determine whether a composition of Dehn twists gives the identity in the mapping class group.