FindInstance and boolean minimization

[anandijain]

Hello, I've been using this package a lot recently and I love it.

I know fairly little about this subject, but I have a couple of applications that are really benefiting from the Symbolic polynomial improvements and I'm looking to help out.

The first application is balancing chemical reactions and the second is minimizing boolean expressions.

Lots of chemical reactions have an infinite number of possible stoichiometric coefficients that satisfy being balanced. For example, the following equations represent the reaction OH- + H+ <-> H2O in mathematica. To the human eye, it's clear that the coefficients [1, 1, 1] are the "best". I'd like to try to write my own FindInstance, but having trouble finding resources on how this function is implemented. I'm wondering if you know.

eqs = {
    u[1] == u[3],
    u[1] + u[2] == 2*u[3],
    -u[1] + u[2] == 0
}
possible = FindInstance[eqs, {u[1], u[2], u[3]}, PositiveIntegers]
xs = Map[Last, First[possible]]
sol = xs / Min[xs] # {1,1,1}

My second question is more concrete. I stumbled on https://en.wikipedia.org/wiki/Quine%E2%80%93McCluskey_algorithm, which mentions that it is an analogous algorithm to Buchberger's for boolean expressions. I'm curious if your F4 implementation might be general enough that I don't have to reimplement Quine McCluskey, or what I'd need to do to make it work.

Thanks and feel free to close since this issue is pretty off topic to the core purpose of this package

[sumiya11]

Hi @anandijain , thanks for your kind words.

For the second question, I am not familiar with Quine–McCluskey algorithm, but I think Groebner.jl can help. As I understand, one can quite easily (at least) find all prime implicants with Gröbner bases. Take, for instance, the function from example on wiki:

$$f(A,B,C,D) = A'BC'D' + AB'C'D' + AB'CD' + AB'CD + ABC'D' + ABCD.$$

For any two terms, say $AB'CD'$ and $AB'CD$ we want to encode implication $AB'CD' + AB'CD \rightarrow AB'C$ as a rule for Gröbner computation. First thing that comes to mind is to match each minterm with the corresponding monomial according to the following: $AB'CD'$ becomes $A(1 + B)C(1 + D)$. The idea is simple, if $x = A(1 + B)C(1 + D)$ and $y = A(1 + B)CD$ are members of the ideal, then $x - y = A(1 + B)C$ is as well.

Going to the code

using Groebner, AbstractAlgebra
_, (a, b, c, d) = PolynomialRing(QQ, [:a,:b,:c,:d])

# all minterms in f, x' maps to (1 + x)
F = [
           (1 + a)*b*(1 + c)*(1 + d),
           a*(1 + b)*(1 + c)*(1 + d),
           a*(1 + b)*(1 + c)*d,
           a*(1 + b)*c*(1 + d),
           a*(1 + b)*c*d,
           a*b*(1 + c)*(1 + d),
           a*b*c*(1 + d),
           a*b*c*d
       ]

groebner(F)

4-element Vector{AbstractAlgebra.Generic.MPoly{Rational{BigInt}}}:
 b*c*d + b*c + b*d + b
 a*d + a
 a*c
 a*b + a

Switching back $(1 + x) \rightarrow x'$ we obtain exactly $f(A, B, C, D) = BC'D' + AB' + AD' + AC$ !

So, as I far as I see, Groebner can do the first step from the QMC algorithm. Probably we can figure how to fully minimize the function using only Groebner bases after a closer look.

[sumiya11]

About your first question, I am actually very interested in CRNs myself. I have little knowledge about your particular application, but I believe a lot of work has been done in this direction. I kindly ask @pogudingleb for some insight here

[pogudingleb]

For the first question. It looks like you can get only linear equations from your balancing conditions (correct me if I am wrong). In the case of linear system, finding "good" solution (nonegative, sparse) etc can often will be done using the tools of convex geometry: you intersect the subspace of balanced reactions with a cone of nonegative coefficients. In Julia, one can use Polymake.jl. We did a similar thing with Xinjian Zhang also in the context of reaction networks in this paper, our code is here. Hope this helps.

[anandijain]

Thanks both!

@sumiya11 I am fairly close to implementing what you wrote for boolean minimization in code. This amount of simplification is still really good.

The issue I've run into though is I have everything in Symbolics, (I'm using a version-updated copy of your Symbolics fork). I also was calling Symbolics.groebner_basis to make use of your code to convert Symbolics to DynamicPolynomials. I'm realizing this won't work since we need to do everything in a ring, not a field since the groebner basis had expressions with - and ^, which would obviously mess up my rules to convert from the monomial form back to boolean expressions.

I'm wondering if you have any codes to translate to AA. If not I can hack it together relatively easily I think.

edit: Wait actually, some of this might be wrong, (I'm not the best at math) but the issue of how to take the a*b + a and factor it back to (1+b)a is the main challenge

[sumiya11]

I am fairly close to implementing what you wrote for boolean minimization in code.

That's great!

I am not completely sure what you mean by ring vs. field, could you elaborate?

I probably have no functions for AA <--> DP conversions.. Btw, Groebner.jl can work with DynamicPolynomials.jl directly:

using Groebner, DynamicPolynomials
@polyvar a b c d

F = [
           (1 + a)*b*(1 + c)*(1 + d),
           a*(1 + b)*(1 + c)*(1 + d),
           .....
       ]

groebner(F)

> 4-element Vector{Polynomial{true, Int64}}:
 ad + a
 ac
 ab + a
 bcd + bc + bd + b

[sumiya11]

For factoring $bcd + bc + bd + b$ and others to $b(1 + c)(1 + d)$ you can use division; for example (example with AA 😄 )

using AbstractAlgebra: divides
R, (a, b, c, d) = PolynomialRing(QQ, ["a","b","c","d"])

f = b*c*d + b*c + b*d + b

divides(f, b)  # f / b
> (true, c*d + c + d + 1)  # implying b

f = c*d + c + d + 1

divides(f, c)  # f / c
> (false, 0) # implying 1+c

divides(f, 1 + c)  # f / (1 + c)
> (true, d + 1)

f = d + 1  

#  and so on..

[anandijain]

Ahh, I was looking for a divides in symbolics, that should help. I'll rewrite my factorization code to use that

https://github.com/anandijain/Boolin.jl/blob/aj/simplify/test/simplify.jl#L178 is the code that I've been messing around with.

I think what I meant is that it's easiest to convert back to boolean if I only allowed the ring operations (+, *) in the resulting expressions from groebner, since I'm blanking on what the rule should be to convert -x1. I don't think its !x1

I think I might have been staring at this for too long, there could be just dumb bugs of mine in the rules or the factoring.

Thanks for walking me through this, shouldn't be too much more to get it relatively robust

[sumiya11]

Aah, I see, good luck with that!

Let me know if you'd need help with Groebner stuff

[sumiya11]

A thought on ^ If each of the variables is either 0 or 1, then x^2 equals x. Could this be the case in your application?