Update random(QQ) to use standard uniform distribution

[d-torrance]

This is split off from #3478. I think it makes sense to keep them separate since that deals with the Probability package and this deals with how rationals are randomly generated in the engine. The big change to the commits that were removed from that PR is the proposed probability distribution (no longer normal).

Rather than using the quotient of two discrete random variables with the uniform distribution on {1, ..., h}, we use the continuous uniform distribution on [0, 1], but rounded to the nearest rational number with denominator bounded by the number given by the Height option (default 10). For example:

i1 : random(QQ^3, QQ^3, Height => 3)

o1 = | 2/3 2/3 1/3 |
     | 1   1   2/3 |
     | 1   2/3 1/2 |

              3       3
o1 : Matrix QQ  <-- QQ

i2 : random(QQ^3, QQ^3)

o2 = | 1/8 1   7/8  |
     | 1/3 5/6 9/10 |
     | 1   3/4 1/7  |

              3       3
o2 : Matrix QQ  <-- QQ

i3 : random(QQ^3, QQ^3, Height => 100)

o3 = | 61/65 33/59 3/16  |
     | 15/26 7/26  31/39 |
     | 3/26  7/48  62/95 |

              3       3
o3 : Matrix QQ  <-- QQ

[d-torrance]

I think the build failures are coming from examples that use random in some way and that have some nonzero probability of failing. This PR is changing things so that we're only calling the pseudorandom number generator once instead of twice for each call to rawRandomQQ, we're getting different random numbers and getting some unexpected results.

In particular, CoincidentRootLoci::realrank sometimes calls QEPCAD, which isn't available in the macOS builds, and GraphicslModelsMLE::solverMLE sometimes results in division by zero.

I was able to reproduce both of these pretty quickly using Macaulay2 1.24.05 on Debian unstable, i.e., before the proposed changes:

i1 : needsPackage "CoincidentRootLoci";

i2 : realrank randomBinaryForm(6, 4, 4)

o2 = 4

i3 : realrank randomBinaryForm(6, 4, 4)

o3 = 4

i4 : realrank randomBinaryForm(6, 4, 4)
stdio:4:10:(3):[1]: error: the ideal is not apolar

i5 : realrank randomBinaryForm(6, 4, 4)

o5 = 4

i6 : realrank randomBinaryForm(6, 4, 4)

o6 = 4

i7 : realrank randomBinaryForm(6, 4, 4)

o7 = 4

i8 : realrank randomBinaryForm(6, 4, 4)
/bin/sh: 1: qepcad: not found
stdio:8:1:(3): error: error occurred while executing QEPCAD. Please make sure that it is installed and configured correctly.

i1 : needsPackage "GraphicalModelsMLE";
-- storing configuration for package Graphs in /root/.Macaulay2/init-Graphs.m2
-- storing configuration for package NumericalAlgebraicGeometry in /root/.Macaulay2/init-NumericalAlgebraicGeometry.m2
-- storing configuration for package Bertini in /root/.Macaulay2/init-Bertini.m2

i2 : G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});

i3 : U = matrix{{6.2849049, 10.292875, 1.038475, 1.1845757}, {3.1938475, 3.2573, 1.13847, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};

                4         4
o3 : Matrix RR    <-- RR
              53        53

i4 : solverMLE(G,U,RealPrecision=>10)

o4 = (1.83311, | 4.52904   7.81123   -.0250185 .309951   |, 2)
               | 7.81123   14.3732   -.0382332 .473665   |
               | -.0250185 -.0382332 .00337164 -.0418254 |
               | .309951   .473665   -.0418254 .742627   |

o4 : Sequence

i5 : solverMLE(G,U,RealPrecision=>10)
stdio:5:1:(3): error: atttempt to divide by zero

I'm thinking of replacing the failing examples with canned ones for the time being and also opening issues for both of these examples. Does that sound reasonable?

[mahrud]

Instead of canning the examples, try setting the random seed to something that prevents failure.

[d-torrance]

For now, I've just added a setRandomSeed line to the offending examples, but I wonder if that might be confusing to a user reading the documentation?

Maybe at some point we should add some way of setting the random seed for a particular example (or test)?

[mahrud]

For now, I've just added a setRandomSeed line to the offending examples, but I wonder if that might be confusing to a user reading the documentation?

Maybe at some point we should add some way of setting the random seed for a particular example (or test)?

This kind of situation is fairly rare, so I think it's fine. Also, I think it's helpful for users to see that they can produce the results identical to the documentation.

[mahrud]

I'm happy with this PR, but given that it's changing the behavior of random rather significantly, I think it would be good to have a second opinion from @mikestillman or @antonleykin.

[d-torrance]

I'm happy with this PR, but given that it's changing the behavior of random rather significantly, I think it would be good to have a second opinion from @mikestillman or @antonleykin.

Yes, I agree!

[mahrud]

While we're here, could you also add something like random RR for QQ? e.g.

random QQ := QQ => opts -> x -> x * random(QQ, opts)

A problem with this method is that the height is then skewed (e.g. 6 * random(QQ, Height => 3) always gives integers), so maybe we want random x + random(QQ, opts) instead? I'm not sure what's the best approach.

[d-torrance]

I pushed a commit adding random(QQ) by first calling random(RR) and then rounding it. That seems to work pretty well:

i1 : apply(20, i -> random(10/1))

      19  11  59  25  7  5  17  89  4  8  33  31  31  11  53     17  64  31  53
o1 = {--, --, --, --, -, -, --, --, -, -, --, --, --, --, --, 4, --, --, --, --}
      10   5   6   8  6  2   3  10  7  5   4  10   4   3   9      6   7   6   6

o1 : List

[d-torrance]

I went ahead and added more random methods, so now pretty much any Number object or pair of Number objects should do something sensible.

For the complex numbers, if we give it $a + bi$, it will give us something in the rectangle $[0, a]\times [0, b]$ on the complex plane, and if we give it $a + bi$ and $c + di$, it will give us something in $[a,c]\times[b, d]$. (The order on the endpoints may be reversed if necessary.)

i1 : random(1, pi)

o1 = 1.47764022754653

o1 : RR (of precision 53)

i2 : random(1/2, ii)

o2 = .398515040066406+.230213567526911*ii

o2 : CC (of precision 53)

[mahrud]

I went ahead and added more random methods, so now pretty much any Number object or pair of Number objects should do something sensible.

I don't like the "pair of numbers" methods and would even advocate for deprecating random(ZZ, ZZ). They overcrowd the methods on random for very little gain, since they are simply for non-integers random(min, max) = min + max * random class max or min + random max for integers. The reason I wanted random QQ is that otherwise getting the height correct is nontrivial.

[d-torrance]

Fair enough -- I'll remove them. It was redundant with random uniformDistribution(a, b) from the Probability package anyway

[d-torrance]

I'm wondering if maybe the support of random(QQ) should be [0, opts.Height] instead of [0, 1]? That would keep it consistent with the current behavior, and might cut down on the random seed failures since I think now there's a higher probability of singular random rational matrices, etc, with the smaller support

[mahrud]

Yeah, maybe. Height of a rational number a/b is usually defined as max(a,b) so you're right that that would make more sense. I guess I'm just annoyed that Height has different effects for ZZ and QQ and seemingly no effect for RR or CC. Can we make these consistent somehow and document them all in the same page?

Also, should random CC give a random complex number on the unit circle instead of a random complex number in the unit square? Not sure which is preferable.

[d-torrance]

I've updated the support so that random QQ returns a random rational number in [0, opts.Height]. At least locally, this fixes it so that only one of the tests/examples that previously needed a different random seed still needs one.

So for RR and CC, are you suggesting that we adjust it so that the support is [0, opts.Height] and [0, opts.Height] x [0, opts.Height], respectively (or [0, opts.Height] x [0, 2$\pi$] in polar coords if we switch to a circle in the complex plane)?.

[d-torrance]

Just pushed another commit fixing a bug I found working on this:

Before

i1 : random(QQ, Height => 0)
Floating point exception (core dumped)

After

i1 : random(QQ, Height => 0)
stdio:1:6:(3): error: expected a positive height

[mahrud]

I think it might be simpler to just assert that Height is positive in the top-level, no?

[d-torrance]

My concern was that we also want to catch this when making random matrices, which is a different call to the engine that isn't specific to the coefficient ring. Is height always positive for all the various meanings of height we might want for random? Or could it be zero/negative in other contexts?

[d-torrance]

I keep re-thinking what the support of random QQ should be lol...

Considering that the height of a rational number is the max of the absolute values of its numerator and denominator, and the proposed behavior would allow larger numerators, then maybe the current behavior is basically correct, or at least has the correct support. But instead maybe we should have a uniform distribution on that support, rather than the current distribution where 1 is more likely than anything else. So in other words, instead of the distribution being $\frac{X}{Y}$ where $X,Y\in\{1,\ldots,h\}$, we have $\frac{X}{Y}$ where $X,Y\in\{1,\ldots,h\}$ and $\gcd(X,Y)=1$.

So rawRandomQQ could be something like (pseudocode):

while (
    x = random_int(1, h)
    y = random_int(1, h)
    gcd(x, y) != 1) do nothing
return x/y

[d-torrance]

So where did we land on the support? [0, 1] or [0, height]?

[mahrud]

I think new option Support with default value [-1,1]? And perhaps this?

random QQ := o -> r -> random(QQ, o, Support => [0,r])

[mahrud]

If it makes sense to you, I think we should merge this and work on the optional arguments in follow ups.

[d-torrance]

I'm not sure if the current proposed behavior is what we want, though, since the numerator can exceed the height.

[mahrud]

Related: #999 and #2089.