hakaru-dev / hakaru

A probabilistic programming language
BSD 3-Clause "New" or "Revised" License
311 stars 30 forks source link

RoundTrip/t62 failure #100

Open yuriy0 opened 7 years ago

yuriy0 commented 7 years ago
 $ pretty tests/RoundTrip/t62.expected.hk
fn x1 real:
fn x0 real:
if x1 / x0 <= nat2real(0): reject. measure(unit)
else:
  if x1 / x0 <= nat2real(1): weight(real2prob(x1 / x0), return ())
  else:
    if x1 / x0 <= nat2real(2):
      weight(real2prob(nat2real(2) - x1 / x0), return ())
    else: reject. measure(unit)

 $ hk-maple tests/RoundTrip/t62.expected.hk
fn x1 real:
fn x0 real:
if signum(x1) * x0 < +0/1: reject. measure(unit)
else:
  if x1 / x0 <= +1/1 && +0/1 < x1 / x0:
    weight(real2prob(x1 / x0), return ())
  else:
    weight(real2prob((x1 * (-1/1) + x0 * (+2/1)) / x0), return ())

This almost roundtrips, but the two conditions whose body is reject get combined, and the resulting condition is solved:

> solve( {Or(And(0 < x1/x0, 1 < x1/x0, 2 < x1/x0), x1/x0 <= 0)}, [x0], useassumptions=true ) assuming (x0::real, x1::real)
  [[signum(x1) x0<0]]

This is clearly the wrong answer, so I think it's a bug. Interestingly, solving for [x1] instead of [x0] gives

piecewise(x0 < 0, [[x1 < 2*x0], [0 <= x1]], x0 = 0, [], 0 < x0, [[2*x0 < x1], [x1 <= 0]])

which is a correct answer (but probably not all that useful; it will be thrown away by KB). I'm guessing we want this program to roundtrip, so throwing away both of these results seems to be the correct thing to do.


 $ pretty tests/RoundTrip/t62.0.hk
fn x2 real:
fn x1 real:
x0 <~ uniform(nat2real(0), nat2real(1))
if nat2real(0) < x2 / x1 - x0 && x2 / x1 - x0 < nat2real(1):
  return ()
else: reject. measure(unit)

 $ hk-maple tests/RoundTrip/t62.0.hk
fn x2 real:
fn x1 real:
if x2 == x1 && +0/1 < x1 || x2 == x1 && x1 < +0/1: return ()
else:
  if +0/1 < x1 && +0/1 < x2 && x2 < x1 ||
     x1 < +0/1 && x1 < x2 && x2 < +0/1:
    weight(real2prob(x2 / x1), return ())
  else:
    if +0/1 < x1 && x1 < x2 && x2 < x1 * (+2/1) ||
       x1 < +0/1 && x2 < x1 && x1 * (+2/1) < x2:
      weight(real2prob(+1/1 + (x2 * (-1/1) + x1) / x1), return ())
    else: reject. measure(unit)

Here we get a pretty nice solution. It could be nicer in terms of the operational cost of computing the inequalities. The first if could be written x2 == x1 && (x /= 0). The other two conditions feel like they could be simpler, but I don't see to what they could be simplified. In either case, we don't really have a mechanism for simplifying conditions in this way; Maple will basically always gives us conditions in disjunctive normal form. Maybe that's good enough.

The test case is currently testing that t62.0 simplifies to t62.expected but there is really no way that will happen; they should be separate tests.

ccshan commented 7 years ago

Regarding the output of simplifying t62.expected.hk: I agree that the result is buggy and I agree that the signum is weird but it doesn't seem to me that the generated branch with signum is what makes the result buggy. I'm more concerned with the situation where x1/x0 = 0 or x1/x0 > 2.

It's too bad if t62.0.hk and t62.expected.hk don't simplify to the same program without random choices.

yuriy0 commented 7 years ago

I'm more concerned with the situation where x1/x0 = 0 or x1/x0 > 2.

I'm not sure I understand. Did you mean x1/x0 <= 0 or x1/x0 > 2 (that is the condition which gets solved to a signum)? In either case, could you elaborate as to what you mean by 'concerned'?

It's too bad if t62.0.hk and t62.expected.hk don't simplify to the same program without random choices.

I agree they should simplify to the same program (although I'm not sure what you mean by 'random choices'). This issue is due to what Domain can and can't represent. In particular, in t62.0, Domain gets to see an expression with a single 'body'; because of that, it gets a chance to (an successfully does) solve the domain problem. After that, the integrals are made innermost and elimed, which produces different weights on the body (which is the applyintegrand).

In t62.expected, the weights are there to begin with, so they prevent Domain from doing anything at all. That is, extracting the shape fails because the values of the Partition have different weights. In order for Domain to get its shot, it would need the ability to represent weights, which requires a non-trivial change to the Domain type. With such a change, Domain improvement could actually happen - as it stands, Domain improvement (and thus the call to SemiAlgebraic) isn't possible. Perhaps a call to SemiAlgebraic is needed elsewhere, but that would duplicate much of the code of Domain.

There is also a much more minor issue. There are no integrals (i.e. no Domain bounds), so Domain won't even try to do anything. This would be fairly easy to fix, although care would have to be taken to not call Domain simplifiers which don't have a chance of succeeding without bounds, so as to avoid lots of extra work in the general case.

ccshan commented 7 years ago

To elaborate on my "concern":

By "random choices", I mean <~ or <|> basically. Happily t62.0.hk and t62.expected.hk both simplify to programs without random choices; it's just too bad that it's not the same program without random choices.

JacquesCarette commented 7 years ago

Regarding the weights that block things: this is part of the reason why we used to push the weights in as far as possible, to expose the control flow to simplification. Then, as the end, you want to do the exact opposite: extract the weights out as far as possible, to help in recognition!

yuriy0 commented 7 years ago

To elaborate on my "concern": [...]

I agree that this is the issue, but I am fairly sure that the issue is caused precisely by the buggy result of solve.

I visualize it as follows:

plots:-inequal( Or(And(0 < x1/x0, 1 < x1/x0, 2 < x1/x0), x1/x0 <= 0), x1=-10..10, x0=-10..10 )

cond0

and

plots:-inequal( signum(x1)*x0<0, x1=-10..10, x0=-10..10 )

cond1

and (solving for x1 instead of x0, and putting the result in a form appropriate for inequal)

plots:-inequal( {Or(0 <= x0, 0 <= x1, x1 < 2*x0), x0 <> 0, Or(x0 <= 0, x1 <= 0, 2*x0 < x1)}, x1=-10..10, x0=-10..10 )

cond2

it's just too bad that it's not the same program without random choices.

Indeed. But I can't see any (simpler) way to get that behaviour without giving Domain the ability to represent weights, and that isn't a simple change at all.

yuriy0 commented 7 years ago

Regarding the weights that block things: this is part of the reason why we used to push the weights in as far as possible, to expose the control flow to simplification.

In this case, pushing the weights down wouldn't help us. Indeed, there is no place we can put the weights to allow Domain a chance to solve the domain problem, except in the Domain shape itself. We end up with an expression of the form

piecewise(C0, h(x), C1, -h(x), C2, 0)

which doesn't allow Domain to extract a shape and a body, since there is no common body (or, the common body is the entire expression, and Domain doesn't get to see the conditions, so it doesn't get to solve them).

JacquesCarette commented 7 years ago

What is the current status of this one?

yuriy0 commented 7 years ago

t62.0 is unchanged; I don't think I'll be revisiting it for a while. It is somewhat of a tricky edge case, in terms of the solution from SemiAlgebraic being slightly sub-optimal, and in terms improving the conditions any more with solve isn't very fruitful (naturally, as it probably just calls SemiAlgebraic on them again).

t62.expected is also unchanged. Making it simplify to the same program as t62.0 is entirely non-trivial, so I don't think I'll revisit that for a while either. The buggy output of solve needs to be avoided somehow, but how? Can we just throw away any solution with signum in it? Is such a solution ever a good solution? I have no other examples (in the wild) of solve producing signum.

ccshan commented 7 years ago

Maybe it's worth telling Maplesoft about this solve bug?

JacquesCarette commented 7 years ago

Yes, it is - @yuriy0, could you do this please?