Maple bug in which inequalities in {con/dis}junction with 'chilled' subexpressions are discarded.

[yuriy0]

Consider the following program fragment, which occurs as part of simplifying the program in #81:

(sum(piecewise(Or(n1 = wordUpdate, kkf2 <> Hakaru:-idx(z, n1), Hakaru:-idx(w, wordUpdate) <> Hakaru:-
idx(w, n1)), 0, And(kkf2 = Hakaru:-idx(z, n1), Hakaru:-idx(w, wordUpdate) = Hakaru:-idx(w, n1), n1 <> 
wordUpdate), 1), n1 = 0 .. Hakaru:-size(w)-1)+1)*(sum(piecewise(Or(n1 = wordUpdate, kkf2 <> Hakaru:-
idx(z, n1), Hakaru:-idx(doc, wordUpdate) <> Hakaru:-idx(doc, n1)), 0, And(kkf2 = Hakaru:-idx(z, n1), Hakaru:-
idx(doc, wordUpdate) = Hakaru:-idx(doc, n1), n1 <> wordUpdate), 1), n1 = 0 .. Hakaru:-
size(w)-1)+1)/(sum(piecewise(Or(n1 = wordUpdate, kkf2 <> Hakaru:-idx(z, n1)), 0, And(kkf2 = Hakaru:-idx(z, 
n1), n1 <> wordUpdate), 1), n1 = 0 .. Hakaru:-size(w)-1)+Hakaru:-size(word_prior)), KB:-build_unsafely( 
KB(Introduce(kkf2, HInt(Bound(`>=`, 0), Bound(`<=`, Hakaru:-size(topic_prior)-1))), Introduce(wordUpdate, 
HInt(Bound(`>=`, 0))), Introduce(z, HArray(HInt(Bound(`>=`, 0)))), Introduce(doc, HArray(HInt(Bound(`>=`, 
0)))), Introduce(w, HArray(HInt(Bound(`>=`, 0)))), Introduce(numDocs, HInt(Bound(`>=`, 0))), 
Introduce(word_prior, HArray(HReal(Bound(`>=`, 0)))), Introduce(topic_prior, HArray(HReal(Bound(`>=`, 0))))) 
);

NB:

• this is an expression sequence of length 2
• it contains three occurences of sum
• KB:-build_unsafely was inserted manually so avid readers of this issue tracker could copy-paste from above

Calling simplify_factor_assuming on this expression (sequence) produces completely bogus output. The error occurs inside simplify_assuming, inside the call to simplify. The mistake is that each sum is erroneously simplified to 1.

This is almost certain the same problem as this, since the program mentioned here contains the condition:

And(kkf2 = Hakaru:-idx(z, n1), Hakaru:-idx(doc, wordUpdate) = Hakaru:-idx(doc, n1), n1 <> wordUpdate)

which is the same condition discussed on that commit; and the expression contains additional conditions similar to this one, i.e. expressions which have an "And with an expression containing a 'chilled' [sub-]expression". The discussion there is about solve, but I wouldn't be surprised if solve and simplify are sharing a good deal of machinery (and every Maple function shares Maple evaluation semantics, which seem to be at the core of the issue here).

The particular example in this issue is worked-around by 8773d7011, but it doesn't seem feasible to try to deal with this bug in such a manner in perpetuity, since having access to Maple evaluation semantics, and solve and simplify, all seems very important to the task at hand. More importantly, this fix applies at a certain stage in the simplification pipeline (right before fromLO) so a call to simplify or solve somewhere before could still break.

Note that changing any one of the sums to Sums causes an additional Sum to remain in the output; these Sums don't simplify and so should return unevaluated from simplify. Changing all three to Sum causes simplify (and by extension, simplify_assuming and simplify_factor_assuming) to produce the correct answer. It seems the only way to convince simplify and friends to not evaluate things 'incorrectly' is to use the inert form, but that doesn't work everywhere (see https://github.com/hakaru-dev/hakaru/issues/84#issuecomment-308809937).

In light of this, it seems that rewriting as much of the code as possible to not use sum or product (etc.) anywhere where these expressions may be passed to solve or simplify while also chilled, is the way to go. However, I haven't been able to pinpoint all the places where this perfect storm of circumstances occurs. It seems there a great many places where sum or product are used by themselves, and I suppose that potentially any of these could contain a condition which would trigger this bug.

[JacquesCarette]

I think this is probably one of the issues that @ccshan doesn't understand. I get the gist of it (as you showed me the solve bug), but the big block of code above is too big to be helpful. Although it does seem like it arises from LDA?

Somehow seeing the Maple-level symptoms would help me. And by this I mean stripping away as much of Hakaru as possible, and show a Maple call on Maple expressions that goes astray.

[JacquesCarette]

I got the following reply from Maplesoft:

I am not entirely certain whether a fix is going to be easy isolate in a patch. An easy work around on the user end is to separate out the <> types involving variable disjoint from the other equations (that is what solve needs to do, but it's not clear where to do that step without breaking other stuff - this particular issue has existed through the entire history of solve as near as I can tell).

[JacquesCarette]

And a follow-up. The listing comes out not-so-nice, but should be sufficient to allow a patch.

This doesn’t completely fix the problem (since it might cause down stream side effects in routines that don’t expect <> expressions to appear in solutions), but it might work for you. The if condition above line 19 and line 25 have changed vs. Maple 18+.

SolveTools:-Transformers:-NonEquality:-Apply := proc(syst::SubSystem) local eqs, noneqs, linnoneqs, news, vars, i, badnoneqs, neq, sol, sols; 1 noneqs, eqs := selectremove(type,[op(SolveTools:-Utilities:-GetEquations(syst)), op(SolveTools:-Utilities:-GetInequations(syst))],<>); 2 noneqs := map(SolveTools:-Utilities:-SimpleCond,noneqs); 3 noneqs := remove(z -> ormap(y -> evalb(type(y,<) and (lhs(y) = lhs(z) and rhs(y) = rhs(z) or lhs(y) = rhs(z) and rhs(y) = lhs(z))) = true,eqs),noneqs); 4 for i to nops(eqs) do 5 if type(eqs[i],<=) then 6 badnoneqs, noneqs := selectremove(z -> evalb(lhs(z) = lhs(eqs[i]) and rhs(z) = rhs(eqs[i]) or lhs(z) = rhs(eqs[i]) and rhs(z) = lhs(eqs[i])) = true,noneqs); 7 if nops(badnoneqs) <> 0 then 8 eqs := [op(1 .. i-1,eqs), <(op(eqs[i])), op(i+1 .. -1,eqs)] end if end if end do; 9 noneqs := remove(proc (z) local y, subst; try for y in eqs while true do if has(y,lhs(z)) then subst := eval(if(type(y,relation),y,y = 0),lhs(z) = rhs(z)); if SolveTools:-Complexity(subst) < 400 and evalb(coulditbe(subst) = false) then return true end if end if end do; return false catch: return false end try end proc,noneqs); 10 if hastype(eqs,{<, <=}) and 0 < nops(noneqs) then 11 vars := SolveTools:-Utilities:-GetVariables(syst); 12 linnoneqs, noneqs := selectremove(z -> (not has(lhs(z),vars) or type(lhs(z),('linear')(vars))) and (not has(rhs(z),vars) or type(rhs(z),('linear')(vars))),noneqs); 13 linnoneqs := SolveTools:-SortByComplexity(linnoneqs); 14 noneqs := map(z -> lhs(z)-rhs(z) <> 0,noneqs); 15 news := [SolveTools:-Utilities:-SetInequations(SolveTools:-Utilities:-SetEquations(syst,{op(eqs)}),{op(noneqs)})]; 16 for i in linnoneqs do 17 news := map(z -> op([SolveTools:-Utilities:-SetEquations(z,{op(SolveTools:-Utilities:-GetEquations(z)), lhs(i) < rhs(i)}), SolveTools:-Utilities:-SetEquations(z,{op(SolveTools:-Utilities:-GetEquations(z)), rhs(i) < lhs(i)})]),news) end do; 18 return op(news) elif nops(eqs) = 0 and 0 < nops(noneqs) then 19 vars := SolveTools:-Utilities:-GetVariables(syst); 20 sols := {}; 21 for neq in noneqs do 22 sol := {op(Engine:-Recurse({lhs(neq) = rhs(neq)},{},vars))}; 23 sol := map(x -> op(map(y -> lhs(y) <> rhs(y),x)),sol); 24 sols := union(sols,sol) end do; 25 return SolveTools:-Utilities:-SetFinished(SolveTools:-Utilities:-SetSolutions(syst,map(SolveTools:-Utilities:-Union,SolveTools:-Utilities:-GetOriginalSolutions(syst),sols)),true) else 26 noneqs := map(z -> lhs(z)-rhs(z) <> 0,noneqs); 27 return SolveTools:-Utilities:-SetInequations(SolveTools:-Utilities:-SetEquations(syst,{op(eqs)}),{op(noneqs)}) end if end proc

[JacquesCarette]

So it seems that this patch is not necessarily safe. But there is a 'cute trick' that one can do:

1. for every inequation a <> b, create a new variable (call it 's')
2. replace the inequation with (a-b)*s = 1

Basically s stands for 1/(a-b). If a-b=0, i.e. a = b, then (a-b)*s=1 has no solutions.

You can invert this process upon return from solve.

[yuriy0]

Several updates:

• I'm not sure how to get the 'trick' to work - it gives a solution which contains a = (s*b + 1)/s, but I'm not sure how to tell if that should become the original inequality, or if it has genuinely been solved to something better. Moreover, I can't actually find a case where we are solving an inequality on integers (or a condition containing one) where we expect anything but the inequality to return unchanged.

• The patch does seem to fix things for solve, and the test suite seems to run the same. It doesn't affect the behaviour inside simplify_assuming though. I tracked down the cause - it is inside assuming, before the call to simplify even happens. The broken function is eval, on line 212 of assuming (a very innocuous looking line: n := eval(n);.). By sheer coincidence, I found a workaround for this one as well (424aead11c16). I am not certain why it works. The execution paths inside assuming with and without this change diverge fairly quickly, and with this change, assuming sees the unchilled expression and starts trying to solve the integer {in}equalities in terms of RootOf. It seems to eventually decide that the RootOf solution is no good and return the original expression.

• This change was introduced because there used to be a call to value and I found (by bisection) that that change caused LDA to fail to simplify. I've tried very hard to find out why that call to value is needed (94240919b763ea646c96f17c923c95ad5e217909 and ba95d0f309da979d86e6cd37b1209dcdfccbe541), but as far as I can tell, there is absolutely no difference between calling value or calling eval_for_Simplify or even just subs(Sum=sum, ..). And yet, removing the call to value_for_Sum causes LDA to simplify to a different program! This is only tangentially related to this issue (due to the fact that the solve and assuming bugs don't happen when we have a Sum, but they do happen when we have a sum). It seems the correct way to fix this is to track down whichever code is relying on receiving sum and not Sum (it is probably a pattern match on sum) and to fix that code. This was my original plan for fixing the regression in simplifying LDA, but I didn't figure out what effect value actually had. I think I need another pair of eyes on this - either I've missed something obvious, or value actually has a global, stateful effect which causes the simplifier to work the way it should.

[JacquesCarette]

• To get the 'trick' to work, you might want to substitute s = 1/(a+b) in the solution, simplify, and remove all tautologies.

• what is n when that line 212 is executed? There are times where things could get weird - especially if that n has somehow captured n itself as part of its value. Details matter quite a bit. [You might want to ask John some questions related to this, if you can boil it down to Maple-only (aka Hakaru-less) calls].

• Feels like you should split off this last point into a precise issue of its own.

[yuriy0]

• That's what I've tried (I think):

  
  > with(Hakaru): with(KB):
  > expr := And(idx(w, n2) = idx(w, n1), a <> b):
  > expr := And( op(1, expr), (a-b)*S=1 );
                        expr := And(Hakaru:-idx(w, n2) = Hakaru:-idx(w, n1), (a - b) S = 1)

expr := solve( KB:-chill(expr) ); expr := remove(x -> x::= and lhs(x)::name and evalb(x), expr); S b + 1 expr := {S = S, a = -------, b = b, Hakaru:-idx[w, n1] = Hakaru:-idx[w, n2], Hakaru:-idx[w, n2] = Hakaru:-idx[w, n2]} S

                                     S b + 1
                        expr := {a = -------, Hakaru:-idx[w, n1] = Hakaru:-idx[w, n2]}
                                        S

expr := simplify(subs(S=(1/(a+b)), expr)); expr := {a = 2 b + a, Hakaru:-idx[w, n1] = Hakaru:-idx[w, n2]}

• It is [ simplify@(e->subsindets(e,Partition,Partition:-PartitionToPW))( .. ) ], where .. is the first expression in the expression sequence at the top of the first post of this issue (strangely, I can't copy it from Maple). That simplify never sees that big expression.

• Done: #88

[JacquesCarette]

• typo: that should have been 1/(a-b) but you did 1/(a+b). (my typo too).

• I can't debug that now. Tomorrow or Thursday, I hope.

[yuriy0]

I tried that too:

> expr := simplify(subs(S=(1/(a-b)), expr));
expr := {a=a,Hakaru:-idx[w,n1]=Hakaru:-idx[w,n2]}

Maybe Maple is actually using your trick to solve inequalities and that is what leads to the bug in the first place?

[JacquesCarette]

You replaced the inequality with an equality -- now you need to re-add it back after. So every tautology corresponds to an inequality that you had replaced.

[yuriy0]

Can you give me an example of how this looks where the result should be something other than returning the input? I don't think I understand the logic of replacing the tautology with the original inequality; for example, for

expr := And(idx[w, n2] = idx[w, n1], 2*a <> 3*b)

applying the process produces

expr := {a=b+1/2 a,`Hakaru:-idx`[w,n1]=`Hakaru:-idx`[w,n2]}

there is no tautology to replace, and the original expression can't be simplified, so the process should yield exactly the input. I don't see how the final step could be applied here. Of course, the example is a bit contrived, but I see no reason it shouldn't work.

[JacquesCarette]

The process should introduce s = 1/(2a - 3b). The results should contain a tautology. It feels like there's a bug if the result contains {a = b + 1/2*a}. Maybe due to a typo?

[yuriy0]

Maybe due to a typo?

Indeed, sorry about that. The tautology does indeed occur here again.

But:

I can't actually find a case where we are solving an inequality on integers (or a condition containing one) where we expect anything but the inequality to return unchanged.

still holds. I think I need to find and work through such a case (one with non-trivial solution) to be confident that I understand how to write a program to apply this trick in general.

But it still seems that #88 is the better solution - I will take Ken's advice on that issue and see if I can pinpoint where execution begins to diverge in terms of sum vs. Sum. And at least for now, the patch and the workaround in simplify_assuming seem to avoid the bug. According to Maplesoft, the patch might not be perfect, but I haven't encountered it breaking things (with one caveat: when solve is passed an un-chilled expression, it errors out with weird errors. I believe this is an improvement over the old behaviour, which was to return RootOfs. We shouldn't be passing unchilled things to solve anyways!). All this to say, I think this issue should be closed, at least for now.

[JacquesCarette]

For the integer case, it is a little tricky. We can't rely on isolve (Maple's routine dedicated to solving over the integers) as it is extremely weak. We are better off using solve which solves the 'relaxed' problem (over the complex, in general, but we do try hard to say it is to be real), and then map things down to integers.

In any case, if this isn't a 'live' issue, in that it is causing test failures, we can close it, and come back to it if it pops up again.

[yuriy0]

It turns out my previous fix wasn't exactly ideal. It causes the following call to happen (it's quite long but it is of the form simplify(chill(..)) assuming ..):

simplify(KB:-chill(product(product(Hakaru:-idx(xsf,j),j = 0 .. piecewise(i = docUpdate,zNewb,Hakaru:-idx(z,i))-1),i = 0 .. Hakaru:-size(z)-1)*product(piecewise(piecewise(i = docUpdate,zNewb,Hakaru:-idx(z,i))+1 = Hakaru:-size(as),1,1-Hakaru:-idx(xsf,piecewise(i = docUpdate,zNewb,Hakaru:-idx(z,i)))),i = 0 .. Hakaru:-size(z)-1)*sum(product(Hakaru:-idx(xsf,i),i = 0 .. k-1)*piecewise(k+1 = Hakaru:-size(as),1,1-Hakaru:-idx(xsf,k)),k = 0 .. Hakaru:-size(as)-1)^(-Hakaru:-size(z))*2^(-Hakaru:-size(z)-Hakaru:-size(as))*exp(-1/2*sum(Hakaru:-idx(t,k)^2,k = 0 .. Hakaru:-size(z)-1))*(2^Hakaru:-size(z))^(1/2)*Pi^(-Hakaru:-size(z))*(2^Hakaru:-size(as))^(1/2)/(Pi^(Hakaru:-size(as)-Hakaru:-size(z)))^(1/2)*product(Hakaru:-idx(xsf,i)^(sum(Hakaru:-idx(as,k),k = i+1 .. Hakaru:-size(as)-1)-1),i = 0 .. Hakaru:-size(as)-2)*product((1-Hakaru:-idx(xsf,i))^(Hakaru:-idx(as,i)-1),i = 0 .. Hakaru:-size(as)-2)/product(Beta(Hakaru:-idx(as,i),sum(Hakaru:-idx(as,k),k = i+1 .. Hakaru:-size(as)-1)),i = 0 .. Hakaru:-size(as)-2)*product(piecewise(csgn(1/2*sum(piecewise(k = piecewise(kj = docUpdate,zNewb,Hakaru:-idx(z,kj)),1,0),kj = 0 .. Hakaru:-size(z)-1)+1/2) = 1,2*exp(1/2*sum(piecewise(k = piecewise(kh = docUpdate,zNewb,Hakaru:-idx(z,kh)),Hakaru:-idx(t,kh),0),kh = 0 .. Hakaru:-size(z)-1)^2/(sum(piecewise(k = piecewise(kj = docUpdate,zNewb,Hakaru:-idx(z,kj)),1,0),kj = 0 .. Hakaru:-size(z)-1)+1))/(2*sum(piecewise(k = piecewise(kj = docUpdate,zNewb,Hakaru:-idx(z,kj)),1,0),kj = 0 .. Hakaru:-size(z)-1)+2)^(1/2)*Pi^(1/2),infinity),k = 0 .. Hakaru:-size(as)-1))) assuming op([zNewb::integer, 0 <= zNewb, zNewb <= Hakaru:-size[as]-1, Hakaru:-size[xsf] = Hakaru:-size[as]-1, xsf::TopProp, Hakaru:-size[t] = Hakaru:-size[z], docUpdate::integer, 0 <= docUpdate, docUpdate <= Hakaru:-size[z]-1, t::TopProp, z::TopProp, as::TopProp, Hakaru:-size[as]::nonnegint, Hakaru:-size[t]::nonnegint, Hakaru:-size[xsf]::nonnegint, Hakaru:-size[z]::nonnegint, Hakaru:-idx[as,i]::real, Hakaru:-idx[as,k]::real, Hakaru:-idx[t,k]::real, Hakaru:-idx[t,kh]::real, Hakaru:-idx[xsf,i]::real, Hakaru:-idx[xsf,j]::real, Hakaru:-idx[xsf,k]::real, Hakaru:-idx[xsf,piecewise(i = docUpdate,zNewb,Hakaru:-idx[z,i])]::real, Hakaru:-idx[z,i]::integer, Hakaru:-idx[z,kh]::integer, Hakaru:-idx[z,kj]::integer, 0 <= Hakaru:-idx[as,i], 0 <= Hakaru:-idx[as,k], 0 <= Hakaru:-idx[z,i], 0 <= Hakaru:-idx[z,kh], 0 <= Hakaru:-idx[z,kj], 0 < Hakaru:-idx[xsf,i], 0 < Hakaru:-idx[xsf,j], 0 < Hakaru:-idx[xsf,k], 0 < Hakaru:-idx[xsf,piecewise(i = docUpdate,zNewb,Hakaru:-idx[z,i])], Hakaru:-idx[xsf,i] < 1, Hakaru:-idx[xsf,j] < 1, Hakaru:-idx[xsf,k] < 1, Hakaru:-idx[xsf,piecewise(i = docUpdate,zNewb,Hakaru:-idx[z,i])] < 1])

and this either takes an absurdly long time (at least 90 minutes), or goes into an infinite loop. In either case, clearly assuming should not be allowed to see unchilled idx. When chill is move out of assuming, it takes about 10 seconds.

I did not figure out what exactly is doing all this extra work. There are lots of calls to sum/hypergeom/poss and SolveTools:-CancelInverses, which both make lots of calls to simplify({RootOf(size(_Z))}, RootOf).

This arises from simplifying gmm_gibbs.

This changes achieves the same thing as the one it replaces, and doesn't cause this call to occur.

[JacquesCarette]

idx and size cause solve a lot of headaches. chill is definitely needed, and it needs to be done "for sure" -- which assuming kind of subverts in weird ways. This is definitely the better way to do things.