From 2615da3d2030df3b6c068f83009daea7ea217730 Mon Sep 17 00:00:00 2001
From: Gleb Belov <gleb@ampl.com>
Date: Wed, 1 Mar 2023 22:22:22 +1100
Subject: [PATCH] Improve bounds for Pow() and quadratics

---
 include/mp/flat/constr_prepro.h | 31 ++++++++++++++++++-------------
 include/mp/flat/converter.h     |  6 ++++--
 include/mp/flat/expr_bounds.h   | 13 ++++++++++---
 3 files changed, 32 insertions(+), 18 deletions(-)

diff --git a/include/mp/flat/constr_prepro.h b/include/mp/flat/constr_prepro.h
index cc874c847..76d3fcfdb 100644
--- a/include/mp/flat/constr_prepro.h
+++ b/include/mp/flat/constr_prepro.h
@@ -55,19 +55,24 @@ class ConstraintPreprocessors {
     auto& m = MP_DISPATCH( GetModel() );
     auto lb = std::pow(m.lb(arg), pwr),
         ub = std::pow(m.ub(arg), pwr);
-    if (MPD( is_integer_value(pwr) ) && pwr>=0.0) {
-      // result integer if arg is and integer, >=0 exponent
-      prepro.set_result_type( m.var_type(arg) );
-      if (MPD( is_integer_value(pwr / 2.0) )) {  // exponent is even, >=0
-        bool lb_neg = m.lb(arg)<0.0;
-        bool ub_pos = m.ub(arg)>0.0;
-        if (lb_neg && ub_pos) {
-          ub = std::max(lb, ub); lb = 0.0;
-        } else if (lb_neg) {
-          std::swap(lb, ub);
-        }
-      }
-    }
+		if (MPD( is_integer_value(pwr) )) {
+			if (pwr>=0.0) {
+				// result integer if arg is and integer, >=0 exponent
+				prepro.set_result_type( m.var_type(arg) );
+				if (MPD( is_integer_value(pwr / 2.0) )) {  // exponent is even, >=0
+					bool lb_neg = m.lb(arg)<0.0;
+					bool ub_pos = m.ub(arg)>0.0;
+					if (lb_neg && ub_pos) {
+						ub = std::max(lb, ub); lb = 0.0;
+					} else if (lb_neg) {
+						std::swap(lb, ub);
+					}
+				}
+			}
+		} else {                      // fractional power
+			if (lb<0.0)
+				lb = 0.0;
+		}
     prepro.narrow_result_bounds( std::min(lb, ub),
                           std::max(lb, ub) );
   }
diff --git a/include/mp/flat/converter.h b/include/mp/flat/converter.h
index 0fef9a2c6..2ee525bd5 100644
--- a/include/mp/flat/converter.h
+++ b/include/mp/flat/converter.h
@@ -194,8 +194,10 @@ class FlatConverter :
 	void FixUnusedDefinedVars() {
 		for (auto i=num_vars(); i--; ) {
 			if (HasInitExpression(i) &&
-					! VarUsageRef(i))
-				set_var_ub(i, lb(i));
+					! VarUsageRef(i)) {
+				set_var_lb(i, 0.0);      // fix to 0
+				set_var_ub(i, 0.0);
+			}
 		}
 	}
 
diff --git a/include/mp/flat/expr_bounds.h b/include/mp/flat/expr_bounds.h
index b9d502743..e17535332 100644
--- a/include/mp/flat/expr_bounds.h
+++ b/include/mp/flat/expr_bounds.h
@@ -91,9 +91,16 @@ class BoundComputations {
   std::pair<double, double> ProductBounds(Var x, Var y) const {
     const auto& m = MPCD(GetModel());
     auto lx=m.lb(x), ly=m.lb(y), ux=m.ub(x), uy=m.ub(y);
-    std::array<double, 4> pb{lx*ly, lx*uy, ux*ly, ux*uy};
-    return { *std::min_element(pb.begin(), pb.end()),
-          *std::max_element(pb.begin(), pb.end()) };
+		if (x!=y) {                        // different vars
+			std::array<double, 4> pb{lx*ly, lx*uy, ux*ly, ux*uy};
+			return { *std::min_element(pb.begin(), pb.end()),
+						*std::max_element(pb.begin(), pb.end()) };
+		}
+		return {
+			lx<=0.0 && ux>=0.0 ?   // lb_result: 0 included?
+						0.0 : std::min(lx*lx, ux*ux),
+			std::max(lx*lx, ux*ux)
+		};
   }
 
   /// Add / merge bounds and type