From 8d2ba20f6ea44156c4e149aeeea75c48c6946f14 Mon Sep 17 00:00:00 2001 From: Paul Tristan Wagner Date: Wed, 6 Mar 2024 02:08:54 +0100 Subject: [PATCH] Non-linear optimization, bug fixes in CAD --- .../parse/LinearConstraintParser.java | 20 ---- ...MultivariatePolynomialConstraintLexer.java | 19 +++ ...ultivariatePolynomialConstraintParser.java | 39 +++++- .../satchecking/parse/PolynomialParser.java | 13 -- .../satchecking/smt/VariableAssignment.java | 14 +++ .../satchecking/theory/nonlinear/CAD.java | 105 +++++------------ .../theory/nonlinear/Interval.java | 91 ++++++++++---- .../nonlinear/MultivariatePolynomial.java | 20 ---- .../MultivariatePolynomialConstraint.java | 53 +++++++++ .../theory/nonlinear/Polynomial.java | 46 +++++--- .../theory/nonlinear/RealAlgebraicNumber.java | 24 ++-- .../solver/NonLinearRealArithmeticSolver.java | 111 ++++++++++++++---- 12 files changed, 358 insertions(+), 197 deletions(-) create mode 100644 src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintLexer.java diff --git a/src/main/java/me/paultristanwagner/satchecking/parse/LinearConstraintParser.java b/src/main/java/me/paultristanwagner/satchecking/parse/LinearConstraintParser.java index a19c68e..4422dba 100644 --- a/src/main/java/me/paultristanwagner/satchecking/parse/LinearConstraintParser.java +++ b/src/main/java/me/paultristanwagner/satchecking/parse/LinearConstraintParser.java @@ -14,26 +14,6 @@ public class LinearConstraintParser implements Parser { - public static void main(String[] args) { - LinearConstraintParser parser = new LinearConstraintParser(); - - Scanner scanner = new Scanner(System.in); - String line; - while ((line = scanner.nextLine()) != null) { - try { - LinearConstraint constraint = parser.parse(line); - System.out.println(constraint); - System.out.println("lhs = " + constraint.getLeftHandSide()); - System.out.println("rhs = " + constraint.getRightHandSide()); - System.out.println("bound = " + constraint.getBound()); - System.out.println("lhs - rhs = " + constraint.getDifference()); - } catch (SyntaxError e) { - e.printWithContext(); - e.printStackTrace(); - } - } - } - /* * Grammar for Linear constraints: * ::= '=' diff --git a/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintLexer.java b/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintLexer.java new file mode 100644 index 0000000..6a70b78 --- /dev/null +++ b/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintLexer.java @@ -0,0 +1,19 @@ +package me.paultristanwagner.satchecking.parse; + +import static me.paultristanwagner.satchecking.parse.TokenType.*; + +public class MultivariatePolynomialConstraintLexer extends Lexer { + + public MultivariatePolynomialConstraintLexer(String input) { + super(input); + + registerTokenTypes( + MIN, + MAX, + LPAREN, + RPAREN + ); + + this.initialize(input); + } +} diff --git a/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintParser.java b/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintParser.java index 8024350..f8c7db6 100644 --- a/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintParser.java +++ b/src/main/java/me/paultristanwagner/satchecking/parse/MultivariatePolynomialConstraintParser.java @@ -1,6 +1,5 @@ package me.paultristanwagner.satchecking.parse; -import me.paultristanwagner.satchecking.sat.Result; import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomial; import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint; @@ -10,10 +9,22 @@ import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.Comparison.LESS_THAN; import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.Comparison.LESS_THAN_OR_EQUALS; import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.Comparison.NOT_EQUALS; +import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.MultivariateMaximizationConstraint.maximize; +import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.MultivariateMinimizationConstraint.minimize; import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.multivariatePolynomialConstraint; public class MultivariatePolynomialConstraintParser implements Parser { + /* + * Grammar: + * ::= + * | 'min' '(' ')' + * | 'max' '(' ')' + * + * ::= "=" | "!=" | "<" | ">" | "<=" | ">=" + * + */ + @Override public MultivariatePolynomialConstraint parse(String string) { ParseResult result = parseWithRemaining(string); @@ -26,12 +37,36 @@ public MultivariatePolynomialConstraint parse(String string) { @Override public ParseResult parseWithRemaining(String string) { + Lexer lexer = new MultivariatePolynomialConstraintLexer(string); + if(lexer.canConsumeEither(MIN, MAX)) { + boolean isMin = lexer.canConsume(MIN); + lexer.consumeEither(MIN, MAX); + lexer.consume(LPAREN); + + Parser parser = new PolynomialParser(); + ParseResult result = parser.parseWithRemaining(lexer.getRemaining()); + MultivariatePolynomial p = result.result(); + + lexer.skip(result.charsRead()); + + lexer.consume(RPAREN); + + MultivariatePolynomialConstraint constraint; + if(isMin) { + constraint = minimize(p); + } else { + constraint = maximize(p); + } + + return new ParseResult<>(constraint, lexer.getCursor(), true); + } + Parser parser = new PolynomialParser(); ParseResult pResult = parser.parseWithRemaining(string); MultivariatePolynomial p = pResult.result(); - Lexer lexer = new ComparisonLexer(string.substring(pResult.charsRead())); + lexer = new ComparisonLexer(string.substring(pResult.charsRead())); MultivariatePolynomialConstraint.Comparison comparison = parseComparison(lexer); ParseResult qResult = parser.parseWithRemaining(lexer.getRemaining()); diff --git a/src/main/java/me/paultristanwagner/satchecking/parse/PolynomialParser.java b/src/main/java/me/paultristanwagner/satchecking/parse/PolynomialParser.java index b23c289..1318b44 100644 --- a/src/main/java/me/paultristanwagner/satchecking/parse/PolynomialParser.java +++ b/src/main/java/me/paultristanwagner/satchecking/parse/PolynomialParser.java @@ -3,8 +3,6 @@ import me.paultristanwagner.satchecking.theory.arithmetic.Number; import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomial; -import java.util.Scanner; - import static me.paultristanwagner.satchecking.parse.TokenType.*; import static me.paultristanwagner.satchecking.theory.arithmetic.Number.number; import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomial.constant; @@ -12,17 +10,6 @@ public class PolynomialParser implements Parser { - public static void main(String[] args){ - Scanner scanner = new Scanner(System.in); - PolynomialParser parser = new PolynomialParser(); - - String line; - while ((line = scanner.nextLine()) != null) { - MultivariatePolynomial polynomial = parser.parse(line); - System.out.println(polynomial); - } - } - /* * Grammar for polynomial constraints: * ::= EQUALS 0 diff --git a/src/main/java/me/paultristanwagner/satchecking/smt/VariableAssignment.java b/src/main/java/me/paultristanwagner/satchecking/smt/VariableAssignment.java index ce18f8b..44b2c41 100644 --- a/src/main/java/me/paultristanwagner/satchecking/smt/VariableAssignment.java +++ b/src/main/java/me/paultristanwagner/satchecking/smt/VariableAssignment.java @@ -36,4 +36,18 @@ public String toString() { return builder.toString(); } + + @Override + public boolean equals(Object o) { + if (this == o) return true; + if (o == null || getClass() != o.getClass()) return false; + VariableAssignment that = (VariableAssignment) o; + return Objects.equals(this.keySet(), that.keySet()) && + this.keySet().stream().allMatch(key -> Objects.equals(this.get(key), that.get(key))); + } + + @Override + public int hashCode() { + return Objects.hash(this.keySet().stream().map(key -> Objects.hash(key, this.get(key))).toArray()); + } } diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/CAD.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/CAD.java index bb6c0b9..44a93f0 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/CAD.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/CAD.java @@ -1,92 +1,59 @@ package me.paultristanwagner.satchecking.theory.nonlinear; -import me.paultristanwagner.satchecking.parse.Parser; -import me.paultristanwagner.satchecking.parse.PolynomialParser; -import me.paultristanwagner.satchecking.smt.VariableAssignment; -import me.paultristanwagner.satchecking.theory.arithmetic.Rational; - import java.util.*; import static me.paultristanwagner.satchecking.theory.nonlinear.Cell.emptyCell; import static me.paultristanwagner.satchecking.theory.nonlinear.Interval.IntervalBoundType.OPEN; import static me.paultristanwagner.satchecking.theory.nonlinear.Interval.*; -import static me.paultristanwagner.satchecking.theory.nonlinear.RealAlgebraicNumber.realAlgebraicNumber; public class CAD { - public static void main(String[] args) { - Parser parser = new PolynomialParser(); - MultivariatePolynomial p = parser.parse("x^2 + y^2 - 1"); - MultivariatePolynomial q = parser.parse("x^2 + y^3 - 1/2"); - - RealAlgebraicNumber x = realAlgebraicNumber(parser.parse("x^6-2x^4+2x^2-3/4").toUnivariatePolynomial(), Rational.parse("105/128"), Rational.parse("27/32")); - RealAlgebraicNumber y = realAlgebraicNumber(parser.parse("x^6-1x^4+x^2-1/4").toUnivariatePolynomial(), Rational.parse("1/2"), Rational.parse("3/4")); - - System.out.println(p); - System.out.println(q); - System.out.println(x); - System.out.println(y); - - - - } - - private List variables; - private Set polynomials; - - public CAD(Set polynomials) { - this.polynomials = polynomials; - - Set variablesSet = new HashSet<>(); - for (MultivariatePolynomial polynomial : polynomials) { - variablesSet.addAll(polynomial.variables); - } - this.variables = new ArrayList<>(variablesSet); - - - } - - public Set> compute( + public Set compute( Set constraints ) { - return compute(constraints, false); + return compute(constraints, null); } - public Set> compute( + public Set compute( Set constraints, - boolean onlyEqualities + List variableOrdering ) { - this.polynomials = new HashSet<>(); + Set polynomials = new HashSet<>(); for (MultivariatePolynomialConstraint constraint : constraints) { - this.polynomials.add(constraint.getPolynomial()); + polynomials.add(constraint.getPolynomial()); } - Set variablesSet = new HashSet<>(); - for (MultivariatePolynomial polynomial : polynomials) { - variablesSet.addAll(polynomial.variables); + if(variableOrdering == null) { + Set variablesSet = new HashSet<>(); + for (MultivariatePolynomial polynomial : polynomials) { + variablesSet.addAll(polynomial.variables); + } + variableOrdering = new ArrayList<>(variablesSet); } - this.variables = new ArrayList<>(variablesSet); // phase 1: projection Map> p = new HashMap<>(); - p.put(variables.size(), polynomials); + p.put(variableOrdering.size(), polynomials); + + for (int r = variableOrdering.size() - 1; r >= 1; r--) { + String variable = variableOrdering.get(r); - for (int r = variables.size() - 1; r >= 1; r--) { - String variable = variables.get(r); Set proj = mcCallumProjection(p.get(r + 1), variable); p.put(r, proj); - String previousVariable = variables.get(r); - p.get(r + 1).stream().filter(poly -> !poly.highestVariable().equals(previousVariable)); + String previousVariable = variableOrdering.get(r); + + // todo: highestVariable could cause errors + p.get(r + 1).removeIf(poly -> !poly.highestVariable().equals(previousVariable)); } // phase 2: lifting List> D = new ArrayList<>(); D.add(List.of(emptyCell())); - for (int i = 1; i <= variables.size(); i++) { + for (int i = 1; i <= variableOrdering.size(); i++) { List D_i = new ArrayList<>(); - String variable = variables.get(i - 1); + String variable = variableOrdering.get(i - 1); for (Cell R : D.get(i - 1)) { Map s = R.chooseSamplePoint(); @@ -114,40 +81,26 @@ public Set> compute( } if (sortedRoots.isEmpty()) { - if (!onlyEqualities) { - D_i.add(R.extend(variable, unboundedInterval())); - } + D_i.add(R.extend(variable, unboundedInterval())); } else { - if (!onlyEqualities) { - D_i.add(R.extend(variable, intervalLowerUnbounded(sortedRoots.get(0), OPEN))); - D_i.add(R.extend(variable, intervalUpperUnbounded(sortedRoots.get(sortedRoots.size() - 1), OPEN))); - - } - D_i.add(R.extend(variable, pointInterval(sortedRoots.get(sortedRoots.size() - 1)))); + D_i.add(R.extend(variable, intervalLowerUnbounded(sortedRoots.get(0).copy(), OPEN))); + D_i.add(R.extend(variable, intervalUpperUnbounded(sortedRoots.get(sortedRoots.size() - 1).copy(), OPEN))); + D_i.add(R.extend(variable, pointInterval(sortedRoots.get(sortedRoots.size() - 1).copy()))); } for (int j = 0; j < sortedRoots.size() - 1; j++) { RealAlgebraicNumber a = sortedRoots.get(j); RealAlgebraicNumber b = sortedRoots.get(j + 1); - D_i.add(R.extend(variable, pointInterval(a))); - if (!onlyEqualities) { - D_i.add(R.extend(variable, interval(a, b, OPEN, OPEN))); - } + D_i.add(R.extend(variable, pointInterval(a.copy()))); + D_i.add(R.extend(variable, interval(a.copy(), b.copy(), OPEN, OPEN))); } } D.add(D_i); } - List result = D.get(variables.size()); - Set> assignments = new HashSet<>(); - for (Cell cell : result) { - VariableAssignment assignment = new VariableAssignment<>(cell.chooseSamplePoint()); - assignments.add(assignment); - } - - return assignments; + return new HashSet<>(D.get(variableOrdering.size())); } public Set mcCallumProjection( diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Interval.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Interval.java index b68cbd4..fdc7ced 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Interval.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Interval.java @@ -1,8 +1,9 @@ package me.paultristanwagner.satchecking.theory.nonlinear; -import me.paultristanwagner.satchecking.parse.PolynomialParser; import me.paultristanwagner.satchecking.theory.arithmetic.Number; +import java.util.Comparator; + import static me.paultristanwagner.satchecking.theory.arithmetic.Number.ZERO; import static me.paultristanwagner.satchecking.theory.arithmetic.Number.number; import static me.paultristanwagner.satchecking.theory.nonlinear.Interval.IntervalBoundType.*; @@ -10,13 +11,6 @@ public class Interval { - public static void main(String[] args) { - PolynomialParser parser = new PolynomialParser(); - Polynomial p = parser.parse("x^2").toUnivariatePolynomial(); - Interval interval = interval(number(-1), number(5), CLOSED, CLOSED); - System.out.println(p.evaluate(interval)); - } - public enum IntervalBoundType { UNBOUNDED, OPEN, @@ -93,15 +87,15 @@ public static Interval pointInterval(Number number) { // todo: improve return values when rational numbers are possible public RealAlgebraicNumber chooseSample() { - if (lowerBoundType == UNBOUNDED && upperBoundType == UNBOUNDED) { + if (lowerBoundType == UNBOUNDED && upperBoundType == UNBOUNDED) { // (-oo, oo) return realAlgebraicNumber(ZERO()); - } else if (lowerBoundType == CLOSED) { + } else if (lowerBoundType == CLOSED) { // [a, ?? return lowerBound; // todo: here we might be able to return a rational number - } else if (upperBoundType == CLOSED) { + } else if (upperBoundType == CLOSED) { // ??, b] return upperBound; // todo: here we might be able to return a rational number } - if (lowerBoundType == UNBOUNDED && upperBoundType == OPEN) { + if (lowerBoundType == UNBOUNDED && upperBoundType == OPEN) { // (-oo, b) if (upperBound.isNumeric()) { return realAlgebraicNumber(upperBound.numericValue().subtract(number(1))); } else { @@ -109,7 +103,7 @@ public RealAlgebraicNumber chooseSample() { } } - if (lowerBoundType == OPEN && upperBoundType == UNBOUNDED) { + if (lowerBoundType == OPEN && upperBoundType == UNBOUNDED) { // (a, oo) if (lowerBound.isNumeric()) { return realAlgebraicNumber(lowerBound.numericValue().add(number(1))); } else { @@ -117,26 +111,29 @@ public RealAlgebraicNumber chooseSample() { } } - if (lowerBound.isNumeric() && upperBound.isNumeric()) { - Number rationalMidpoint = - lowerBound.numericValue().add(upperBound.numericValue()).divide(number(2)); + if (lowerBound.isNumeric() && upperBound.isNumeric()) { // (q, r) + Number rationalMidpoint = lowerBound.numericValue().midpoint(upperBound.numericValue()); return realAlgebraicNumber(rationalMidpoint); } - if (!upperBound.isNumeric() && lowerBound.isNumeric()) { - if (upperBound.getLowerBound().equals(lowerBound.numericValue())) { + if (!upperBound.isNumeric() && lowerBound.isNumeric()) { // (q, b) + while (upperBound.getLowerBound().lessThanOrEqual(lowerBound.numericValue())) { upperBound.refine(); } return realAlgebraicNumber(upperBound.getLowerBound()); - } else if (!lowerBound.isNumeric() && upperBound.isNumeric()) { - if (lowerBound.getUpperBound().equals(upperBound.numericValue())) { + } else if (!lowerBound.isNumeric() && upperBound.isNumeric()) { // (a, r) + while (lowerBound.getUpperBound().greaterThanOrEqual(upperBound.numericValue())) { lowerBound.refine(); } return realAlgebraicNumber(lowerBound.getUpperBound()); } + while (lowerBound.getUpperBound().greaterThanOrEqual(upperBound.getLowerBound())) { // (a, b) + lowerBound.refine(); + upperBound.refine(); + } return realAlgebraicNumber(lowerBound.getUpperBound()); } @@ -333,4 +330,58 @@ public String toString() { return sb.toString(); } + + public static class LowerBoundIntervalComparator implements Comparator { + + @Override + public int compare(Interval a, Interval b) { + if(a.lowerBoundType == UNBOUNDED) { // (-oo vs. ?? + if(b.lowerBoundType == UNBOUNDED) { // (-oo vs. (-oo + return 0; + } + + // (-oo vs. ?a + return -1; + } + + if(b.lowerBoundType == UNBOUNDED) { // ?a vs. (-oo + return 1; + } + + if(a.lowerBoundType == CLOSED) { // [a vs. ?? + if(b.lowerBoundType == CLOSED) { // [a vs. [a + if(a.lowerBound.equals(b.lowerBound)) { + return 0; + } else if(a.lowerBound.lessThan(b.lowerBound)) { + return -1; + } else { + return 1; + } + } else { // [a vs. (b + if(a.lowerBound.lessThanOrEqual(b.lowerBound)) { + return -1; + } else { + return 1; + } + } + } + + if(b.lowerBoundType == CLOSED) { // (a vs. [b + if(a.lowerBound.lessThan(b.lowerBound)) { + return -1; + } else { + return 1; + } + } + + // (a vs. (b + if(a.lowerBound.lessThan(b.lowerBound)) { + return -1; + } else if (a.lowerBound.equals(b.lowerBound)) { + return 0; + } else { + return 1; + } + } + } } diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomial.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomial.java index 5529786..56c4f17 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomial.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomial.java @@ -12,30 +12,10 @@ import java.util.*; import java.util.stream.Collectors; -import me.paultristanwagner.satchecking.parse.Parser; -import me.paultristanwagner.satchecking.parse.PolynomialParser; import me.paultristanwagner.satchecking.theory.arithmetic.Number; public class MultivariatePolynomial { - public static void main(String[] args) { - Parser parser = new PolynomialParser(); - Scanner scanner = new Scanner(System.in); - System.out.print("p: "); - MultivariatePolynomial p = - parser.parse("t - y^5 - x^5 + x"); // parser.parse(scanner.nextLine()); - System.out.print("q: "); - MultivariatePolynomial q = - parser.parse("x^6-1/2*x^4+2*x^2-1/2"); // parser.parse(scanner.nextLine()); - System.out.print("var: "); - String var = scanner.nextLine(); - - System.out.println("p = " + p); - System.out.println("q = " + q); - - System.out.println("Resultant[p, q, " + var + "] = " + p.resultant(q, var)); - } - public Map coefficients; public List variables; diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomialConstraint.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomialConstraint.java index 9820c2c..bc63192 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomialConstraint.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/MultivariatePolynomialConstraint.java @@ -102,4 +102,57 @@ public MultivariatePolynomial getPolynomial() { public Comparison getComparison() { return comparison; } + + @Override + public String toString() { + return polynomial + " " + comparison + " 0"; + } + + public static abstract class MultivariateOptimizationConstraint extends MultivariatePolynomialConstraint { + + protected final MultivariatePolynomial objective; + + private MultivariateOptimizationConstraint(MultivariatePolynomial objective) { + super(objective, null); + + this.objective = objective; + } + + public MultivariatePolynomial getObjective() { + return objective; + } + } + + public static class MultivariateMinimizationConstraint extends MultivariateOptimizationConstraint { + + private MultivariateMinimizationConstraint(MultivariatePolynomial objective) { + super(objective); + } + + public static MultivariateMinimizationConstraint minimize(MultivariatePolynomial objective) { + return new MultivariateMinimizationConstraint(objective); + } + + @Override + public String toString() { + return "min(" + objective + ")"; + } + } + + public static class MultivariateMaximizationConstraint extends MultivariateOptimizationConstraint { + + + private MultivariateMaximizationConstraint(MultivariatePolynomial objective) { + super(objective); + } + + public static MultivariateMaximizationConstraint maximize(MultivariatePolynomial objective) { + return new MultivariateMaximizationConstraint(objective); + } + + @Override + public String toString() { + return "max(" + objective + ")"; + } + } } diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Polynomial.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Polynomial.java index cc3f6eb..6eeb283 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Polynomial.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/Polynomial.java @@ -1,7 +1,5 @@ package me.paultristanwagner.satchecking.theory.nonlinear; -import me.paultristanwagner.satchecking.parse.Parser; -import me.paultristanwagner.satchecking.parse.PolynomialParser; import me.paultristanwagner.satchecking.theory.arithmetic.Number; import me.paultristanwagner.satchecking.theory.arithmetic.Rational; @@ -16,16 +14,6 @@ public class Polynomial { - public static void main(String[] args) { - Parser parser = new PolynomialParser(); - Polynomial p = parser.parse("x^5+x^4+x^2+x+2").toUnivariatePolynomial(); - Polynomial q = parser.parse("3x^5-7x^3+3x^2").toUnivariatePolynomial(); - - System.out.println(p); - System.out.println(q); - System.out.println(p.pow(100).squareFreeFactorization()); - } - private final int degree; private final Number[] coefficients; // smallest degree first @@ -205,7 +193,7 @@ public Number content() { } } - return content; + return content.abs(); } public List pseudoDivision(Polynomial other) { @@ -218,6 +206,17 @@ public List pseudoDivision(Polynomial other) { return polynomial(pow).multiply(this).divide(other); } + // pseudo division with absolute leading coefficient + public List pseudoDivision2(Polynomial other) { + int a = degree; + int b = other.degree; + + Number bLC = other.getLeadingCoefficient().abs(); // critical difference + Number pow = bLC.pow(a - b + 1); + + return polynomial(pow).multiply(this).divide(other); + } + public List divide(Polynomial other) { if (other.isZero()) { throw new ArithmeticException("Cannot divide by zero"); @@ -342,17 +341,19 @@ public List sturmSequence() { } List sturmSequence = new ArrayList<>(); - sturmSequence.add(this); - sturmSequence.add(this.getDerivative()); + sturmSequence.add(this.toIntegerPolynomial()); + sturmSequence.add(this.getDerivative().toIntegerPolynomial()); while (true) { Polynomial p = sturmSequence.get(sturmSequence.size() - 2); Polynomial q = sturmSequence.get(sturmSequence.size() - 1); - Polynomial negRem = p.divide(q).get(1).negate(); + Polynomial negRem = p.pseudoDivision2(q).get(1).negate(); if (negRem.isZero()) { break; } + negRem = negRem.divide(constant(negRem.content())).get(0); + sturmSequence.add(negRem); } @@ -405,6 +406,10 @@ public int numberOfRealRoots(Number a, Number b, List sturmSequence) roots--; } + if(roots < 0) { + throw new IllegalStateException("Number of roots cannot be negative"); + } + return roots; } @@ -468,13 +473,18 @@ public Set isolateRoots(Number lowerBound, Number upperBoun } Set leftRoots = new HashSet<>(isolateRoots(lowerBound, split)); - Set rightRoots = isolateRoots(split, upperBound); - leftRoots.addAll(rightRoots); if (hasRealRootAt(split)) { leftRoots.add(realAlgebraicNumber(split)); } + if(leftRoots.size() == numberOfRealRoots) { + return leftRoots; + } + + Set rightRoots = isolateRoots(split, upperBound); + leftRoots.addAll(rightRoots); + return leftRoots; } diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/RealAlgebraicNumber.java b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/RealAlgebraicNumber.java index 552805f..4120f94 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/RealAlgebraicNumber.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/nonlinear/RealAlgebraicNumber.java @@ -3,9 +3,7 @@ import me.paultristanwagner.satchecking.theory.arithmetic.Number; import java.util.List; -import java.util.Map; import java.util.Objects; -import java.util.Set; import static me.paultristanwagner.satchecking.theory.arithmetic.Number.ZERO; import static me.paultristanwagner.satchecking.theory.arithmetic.Number.number; @@ -17,8 +15,7 @@ public class RealAlgebraicNumber { private Number lowerBound; private Number upperBound; - private RealAlgebraicNumber( - Number value, Polynomial polynomial, Number lowerBound, Number upperBound) { + private RealAlgebraicNumber(Number value, Polynomial polynomial, Number lowerBound, Number upperBound, boolean minimal) { this.value = value; this.lowerBound = lowerBound; this.upperBound = upperBound; @@ -27,6 +24,11 @@ private RealAlgebraicNumber( return; } + if(minimal) { + this.polynomial = polynomial; + return; + } + // finding minimal polynomial // todo: improve efficiency List squareFreeFactors = polynomial.squareFreeFactorization(); @@ -59,12 +61,16 @@ private RealAlgebraicNumber( } public static RealAlgebraicNumber realAlgebraicNumber(Number value) { - return new RealAlgebraicNumber(value, null, value, value); + return new RealAlgebraicNumber(value, null, value, value, true); } public static RealAlgebraicNumber realAlgebraicNumber( Polynomial polynomial, Number lowerBound, Number upperBound) { - return new RealAlgebraicNumber(null, polynomial, lowerBound, upperBound); + return new RealAlgebraicNumber(null, polynomial, lowerBound, upperBound, false); + } + + public RealAlgebraicNumber copy() { + return new RealAlgebraicNumber(value, polynomial, lowerBound, upperBound, true); } public boolean isNumeric() { @@ -130,8 +136,10 @@ public Number approximate(Number epsilon) { return numericValue(); } - refine(epsilon); - return lowerBound.add(upperBound).divide(number(2)); + RealAlgebraicNumber copy = copy(); + + copy.refine(epsilon); + return copy.lowerBound.add(copy.upperBound).divide(number(2)); } public boolean isZero() { diff --git a/src/main/java/me/paultristanwagner/satchecking/theory/solver/NonLinearRealArithmeticSolver.java b/src/main/java/me/paultristanwagner/satchecking/theory/solver/NonLinearRealArithmeticSolver.java index 6e469c5..e175f27 100644 --- a/src/main/java/me/paultristanwagner/satchecking/theory/solver/NonLinearRealArithmeticSolver.java +++ b/src/main/java/me/paultristanwagner/satchecking/theory/solver/NonLinearRealArithmeticSolver.java @@ -2,21 +2,24 @@ import me.paultristanwagner.satchecking.smt.VariableAssignment; import me.paultristanwagner.satchecking.theory.TheoryResult; -import me.paultristanwagner.satchecking.theory.nonlinear.CAD; -import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomial; -import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint; -import me.paultristanwagner.satchecking.theory.nonlinear.RealAlgebraicNumber; +import me.paultristanwagner.satchecking.theory.nonlinear.*; +import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.MultivariateMaximizationConstraint; +import me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.MultivariateOptimizationConstraint; +import org.apache.commons.lang3.tuple.Pair; -import java.util.HashSet; -import java.util.Set; +import java.util.*; import static me.paultristanwagner.satchecking.theory.TheoryResult.satisfiable; import static me.paultristanwagner.satchecking.theory.TheoryResult.unsatisfiable; -import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomialConstraint.Comparison.EQUALS; +import static me.paultristanwagner.satchecking.theory.nonlinear.Interval.IntervalBoundType.CLOSED; +import static me.paultristanwagner.satchecking.theory.nonlinear.MultivariatePolynomial.variable; public class NonLinearRealArithmeticSolver implements TheorySolver { + private static final String TARGET_VARIABLE_NAME_BASE = "target_"; + private final Set constraints = new HashSet<>(); + private MultivariateOptimizationConstraint objective; @Override public void clear() { @@ -26,43 +29,111 @@ public void clear() { @Override public void load(Set constraints) { clear(); - this.constraints.addAll(constraints); + + constraints.forEach(this::addConstraint); } @Override public void addConstraint(MultivariatePolynomialConstraint constraint) { + if (constraint instanceof MultivariateOptimizationConstraint optimizationConstraint) { + if (objective != null) { + throw new IllegalArgumentException("More than one optimization objective provided"); + } + + objective = optimizationConstraint; + return; + } + constraints.add(constraint); } @Override public TheoryResult solve() { - boolean onlyEqualities = true; - Set polynomials = new HashSet<>(); + // collect all variables + Set variablesSet = new HashSet<>(); for (MultivariatePolynomialConstraint constraint : constraints) { - polynomials.add(constraint.getPolynomial()); - if (!constraint.getComparison().equals(EQUALS)) { - onlyEqualities = false; + MultivariatePolynomial polynomial = constraint.getPolynomial(); + variablesSet.addAll(polynomial.variables); + } + + // if we have an optimization objective, we need to add a fresh variable and fix the variable ordering + List variableOrdering = null; + String freshVariableName = null; + if(objective != null) { + freshVariableName = freshVariableName(variablesSet); + MultivariatePolynomial freshVariable = variable(freshVariableName); + MultivariatePolynomialConstraint helper = MultivariatePolynomialConstraint.equals(freshVariable, objective.getObjective()); + + constraints.add(helper); + + variableOrdering = new ArrayList<>(); + variableOrdering.add(freshVariableName); + variableOrdering.addAll(variablesSet); + } + + CAD cad = new CAD(); + Set result = cad.compute(constraints, variableOrdering); + + List>> pairs = new ArrayList<>(); + for (Cell cell : result) { + Map samplePoint = cell.chooseSamplePoint(); + VariableAssignment variableAssignment = new VariableAssignment<>(samplePoint); + + pairs.add(Pair.of(cell, variableAssignment)); + } + + if (objective != null) { + Comparator comparator; + + if(objective instanceof MultivariateMaximizationConstraint) { + comparator = new Interval.LowerBoundIntervalComparator().reversed(); + } else { + comparator = new Interval.LowerBoundIntervalComparator(); } + + pairs.sort((a, b) -> comparator.compare(a.getLeft().getIntervals().get(0), b.getLeft().getIntervals().get(0))); } - CAD cad = new CAD(polynomials); - Set> result = cad.compute(constraints, onlyEqualities); + for (Pair> pair : pairs) { + Cell cell = pair.getLeft(); + VariableAssignment variableAssignment = pair.getRight(); - for (VariableAssignment realAlgebraicNumberVariableAssignment : result) { boolean satisfied = true; for (MultivariatePolynomialConstraint constraint : constraints) { - int sign = constraint.getPolynomial().evaluateSign(realAlgebraicNumberVariableAssignment); - if(!constraint.getComparison().evaluateSign(sign)) { + int sign = constraint.getPolynomial().evaluateSign(variableAssignment); + + if (!constraint.getComparison().evaluateSign(sign)) { satisfied = false; break; } } - if(satisfied) { - return satisfiable(realAlgebraicNumberVariableAssignment); + if (satisfied) { + if (objective != null) { + variableAssignment.remove(freshVariableName); + + Interval targetInterval = cell.getIntervals().get(0); + if (targetInterval.getLowerBoundType() != CLOSED) { + return unsatisfiable(constraints); + } else { + return satisfiable(new VariableAssignment(variableAssignment)); + } + } + + return satisfiable(new VariableAssignment(variableAssignment)); } } return unsatisfiable(constraints); } + + private String freshVariableName(Set variables) { + int i = 0; + String freshVariableName = TARGET_VARIABLE_NAME_BASE + i; + while (variables.contains(freshVariableName)) { + freshVariableName = TARGET_VARIABLE_NAME_BASE + i; + i++; + } + return freshVariableName; + } }