One more simplification for Domain to do

[JacquesCarette]

Praveen and I were talking briefly about slice sampling, and this nugget that Rob and I wanted to simplify came up:

restart: with(Hakaru): with(NewSLO): y, kb := KB:-genType(y, HReal(Bound(>,0), Bound(<,densityGaussian(0))), KB:-empty):

m := Bind(Gaussian(0,1),x, piecewise(And(0<y,y<densityGaussian(x)), Weight(densityUniform(y), Ret(x)), Msum())):

fromLO(improve(toLO(m),_ctx=kb),_ctx=kb); Weight(2 sqrt(-ln(2)-ln(Pi)-2ln(y)), Uniform(-sqrt(-ln(2)-ln(Pi)-2ln(y)), sqrt(-ln(2)-ln(Pi)-2ln(y))))

It's great that simplification produces such a nice term, but it's annoying to have to provide the bounds of y in the input kb. Shouldn't the new domain geometry facility figure out those bounds? I.e., I want non-urgently for this to happen:

restart: with(Hakaru): with(NewSLO):

y, kb := KB:-genType(y, HReal(), KB:-empty):

m := Bind(Gaussian(0,1),x, piecewise(And(0<y,y<densityGaussian(x)), Weight(densityUniform(y), Ret(x)), Msum())):

fromLO(improve(toLO(m),_ctx=kb),_ctx=kb);

piecewise(And(0<y, y<densityGaussian(0)), Weight(2 sqrt(-ln(2)-ln(Pi)-2ln(y)), Uniform(-sqrt(-ln(2)-ln(Pi)-2ln(y)), sqrt(-ln(2)-ln(Pi)-2ln(y)))), Msum())

[yuriy0]

This seemed like a particularly interesting problem; and something we should definitely be able to do! So I was compelled to try to fix it.

Two things needed to change for this simplification to occur:

• Domain:-Apply will now place constraints (that is, piecewise/Partition) outside of a subterm if that constraint is free of variables enclosed in that subterm. This was a fairly straightforward change, but could potentially lead to the simplification of lots of new things! Other parts of reduce will occasionally do something like this as a byproduct of other simplification, but this is something that should fire much more reliable now. (026a821ae6)

• The much less trivial change is that the condition y<densityGaussian(x) must become And(y<densityGaussian(x),y<densityGaussian(0)). This is only necessary because solve doesn't get the solution to the former constraint without the latter constraint in the context. I've implemented this with minimize, which could potentially be very costly (and simultaneously do nothing, or even make the condition much worse); so this new simplification is applied very conservatively. (5c64091e3, 0c46ba3c, 1c97d48d0)

The new test case added by 09084c3aff now passes; but I leave the issue open because I feel like the 2nd point needs a much better solution.

[JacquesCarette]

Yes, this is rather a non-trivial issue.

@ccshan I completely understand why this is desirable, but did you think this would be easy/hard?

[ccshan]

I only vaguely thought that the new domain geometry facility would make this "less trivial change" easier, because after all the entailment from 0<y<densityGaussian(x) to 0<y<densityGaussian(0) is a matter of domain geometry. In particular, note that the relationship between ln(y) and x is a quadratic inequality. But I am ok with this change being on the back burner. It would make a nice demo.

[JacquesCarette]

I don't see how this is "domain geometry"? If it were, then minimize/maximize would not be needed. If you can explain your reasoning, then maybe we can "get" this particular case much more easily (i.e. we might have missed a way to "reason out" this case).

It is worth noting that although it is easy for humans to see ln(y) as a (fresh) variable, this is significantly harder to automate reliably.

[ccshan]

In the old example from disintegration, domain geometry is reasoning from 0<x<1, 0<y<1, t=x-y to discover which values of t are possible (i.e., yield a non-empty set of x,y,t triples). Here we'd like to reason from -infinity<x<infinity, 0<y<exp(-x^2/2)/sqrt(2*Pi) to discover which values of y are possible (i.e., yield a non-empty set of x,y pairs).

[JacquesCarette]

And here we can do it because exp is monotone increasing (or, if you prefer, exp(-z) is monotone decreasing). That is very special. And tricky, because exp(-x^2) is not so well-behaved. In other words, back to 'discovering' that ln(y) is a better "coordinate" than y.

[ccshan]

I hope I've managed to answer your question 'how this is "domain geometry"?'!

As for using ln(y) as coordinate, I tried doing it manually as below, and it didn't work just now:

> restart: with(Hakaru): with(NewSLO):
> ly, kb := KB:-genType(ly, HReal(), KB:-empty);
> m := Bind(Gaussian(0,1),x,
       piecewise(And(exp(ly)<density[Gaussian](0,1)(x)),
                 Weight(density[Uniform](0,density[Gaussian](0,1)(x))(exp(ly)),
                        Ret(x)),
                 Msum())):
> fromLO(improve(toLO(m),_ctx=kb),_ctx=kb);
Bind(Lebesgue(-infinity, infinity), x, PARTITION(...))

(The Weight is missing a Jacobian factor, but it doesn't really matter.)

[yuriy0]

the entailment from 0<y<densityGaussian(x) to 0<y<densityGaussian(0) is a matter of domain geometry

Yes! But it just isn't the type of domain geometry we have been solving up until now.

In particular, we have been 'solving' things of the form A = B, where A,B are domains; and by 'solving', I mean replacing A with B (or vice versa) based on which is "simpler" or "better". The most basic explanation of how that is done is to manipulate A into a form which a Maple solver (e.g. solve, SemiAlgebraic, previously LinearMultivariateSystem) understands, call that solver to produce B (which we simply trust will be "simpler"), and just replace A with B.

This particular case has to be solved with a more complicated technique. Here we have A' => A, where A,A' are domains, and A is the input domain. Of course, we cannot simply replace A with A', or v.v.! In the previous case, the reasoning we use to arrive to the answer is obvious; here we need more complicated reasoning: A becomes A' and A; and then A is solved in a context containing A' (which is solved and which is placed in the context could be decided by the variables occurring in each; right now, since {max/min}imize are used to determine A', we know A' must be free of the domain variables, so it is always in the context; this is just an implementation detail) to produce A' and B, where (A = B) assuming A'. The solved domain (B) may actually be weaker than the input (A), since we have assumed A' - so while A' and A = A, it may not be the case that A' and B = A. We want to eventually see the expression B(e) (where e is the body of the expression from which the domain was extracted; and B(.) is Domain:-Apply(B)) because that expression gets recognized as a "nice" Hakaru term. By this point, the information of B will have been moved into the bounds of a binder like Int, because that is the form which gets recognized by fromLO. This simplification was already being done, but is also an essential part of seeing the right output in this case (indeed, it is essentially to basically every domain improvement of which I can conceive). As Ken's last example shows, this is currently quite fragile - the new simplification only happens when a relation is a bound on a variable (one side of the relation is precisely a name), and exp(ly)<density[Gaussian](0,1)(x)) is not a bound on a variable. This relation might (should?) be turned into a bound by KB:-try_improve_exp - I'm not certain this isn't just a very mechanical issue (I only have one test case for this new logic, and one test case does not define a semantics).

To summarize, the call to {min/max}imize is needed only to discover that A' => A given the domain A; everything else is genuine "domain geometry" manipulations; it is just that these new manipulations are almost entirely disjoint, i.e. using different laws of domain geometry (which are really just laws of propositional logic), from the old manipulations which were being done.

[yuriy0]

As for using ln(y) as coordinate, I tried doing it manually as below, and it didn't work just now:

Fixed by 4753a28, and added as a test case by 4e965064.

Please reopen if there are more cases like these which don't simplify.

[ccshan]

I'm not sure why, but the tests "clamp condition to move it out of integral" and "clamp condition to move it out of integral (ln coordinate)" in NewSLOTests.mpl fail nowdays...

[yuriy0]

Thanks for finding this. Should be fixed by a063fc58356bd2b11619974f36aa729ec5daf030.