From e93e42300aa8561f6566e7f60a1f751c94fda427 Mon Sep 17 00:00:00 2001
From: aml <aml@opencascade.com>
Date: Thu, 9 Jul 2015 14:12:21 +0300
Subject: [PATCH] 0026356: Wrong result done by projection algorithm

Changed internal one dimension search algorithm in case of fast changing curve.

Test-case for issue #26356
---
 src/math/math_FunctionSetRoot.cxx | 68 +++++++++++++++++++++++--------
 tests/bugs/modalg_6/bug26356      | 37 +++++++++++++++++
 2 files changed, 88 insertions(+), 17 deletions(-)
 create mode 100644 tests/bugs/modalg_6/bug26356

diff --git a/src/math/math_FunctionSetRoot.cxx b/src/math/math_FunctionSetRoot.cxx
index a57b5177b2..03510af9e4 100644
--- a/src/math/math_FunctionSetRoot.cxx
+++ b/src/math/math_FunctionSetRoot.cxx
@@ -125,18 +125,16 @@ Standard_Boolean MyDirFunction::Value(const Standard_Real x,
     p = Dir->Value(i);
     P->Value(i) = p * x + P0->Value(i);
   }
-  if( F->Value(*P, *FV) ) {
+  if( F->Value(*P, *FV) )
+  {
 
-    Standard_Real aVal = 0.;
+    Standard_Real aVal = 0.0;
 
-    for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++) {
+    for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++)
+    {
       aVal = FV->Value(i);
-      if(aVal < 0.) {
-        if(aVal <= -1.e+100) // Precision::HalfInfinite() later
-          //       if(Precision::IsInfinite(Abs(FV->Value(i)))) {
-          //	fval = Precision::Infinite();
-          return Standard_False;
-      }
+      if(aVal <= -1.e+100) // Precision::HalfInfinite() later
+        return Standard_False;
       else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
         return Standard_False;
     }
@@ -317,7 +315,8 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
 
   // (2) Si l'on a pas assez progresser on realise une recherche 
   //     en bonne et due forme, a partir des inits precedents
-  if ((fsol > 0.2*PValue) && (tol1d < 0.5)) {
+  if ((fsol > 0.2*PValue) && (tol1d < 0.5)) 
+  {
 
     if (tsol <0) {
       ax = tsol; bx = 0.0; cx = 1.0;
@@ -328,17 +327,52 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
     FSR_DEBUG(" minimisation dans la direction");
 
     math_BrentMinimum Sol(tol1d, 100, tol1d);
-    Sol.Perform(F, ax, bx, cx);
 
-    if(Sol.IsDone()) {
-      if (Sol.Minimum() <= Result) {
+    // Base invocation.
+    Sol.Perform(F, ax, bx, cx);
+    if(Sol.IsDone())
+    {
+      if (Sol.Minimum() <= Result)
+      {
         tsol = Sol.Location();
         good = Standard_True;
-        FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
-      }
-    }
+        Result = Sol.Minimum();
+
+        // Objective function changes too fast ->
+        // perform additional computations.
+        if (Gradient.Norm2() > 1.0 / Precision::SquareConfusion() &&
+            tsol > ax && 
+            tsol < cx) // Solution inside of (ax, cx) interval.
+        {
+          // First and second part invocation.
+          Sol.Perform(F, ax, (ax + tsol) / 2.0, tsol);
+          if(Sol.IsDone())
+          {
+            if (Sol.Minimum() <= Result)
+            {
+              tsol = Sol.Location();
+              good = Standard_True;
+              Result = Sol.Minimum();
+            }
+          }
+
+          Sol.Perform(F, tsol, (cx + tsol) / 2.0, cx);
+          if(Sol.IsDone())
+          {
+            if (Sol.Minimum() <= Result)
+            {
+              tsol = Sol.Location();
+              good = Standard_True;
+              Result = Sol.Minimum();
+            }
+          }
+        }
+      } // if (Sol.Minimum() <= Result)
+    } // if(Sol.IsDone())
   }
-  if (good) {
+
+  if (good)
+  {
     // mise a jour du Delta
     Dir.Multiply(tsol);
   }
diff --git a/tests/bugs/modalg_6/bug26356 b/tests/bugs/modalg_6/bug26356
new file mode 100644
index 0000000000..0049b2f430
--- /dev/null
+++ b/tests/bugs/modalg_6/bug26356
@@ -0,0 +1,37 @@
+puts "========"
+puts "OCC26356"
+puts "========"
+puts ""
+#############################################
+# Wrong result done by projection algorithm
+#############################################
+
+restore [locate_data_file OCC26356-f.brep] b1
+restore [locate_data_file OCC26356-w.brep] b2
+
+explode b2 v
+
+point p1 31350.009765625 7100 -2.17374844144233e-013
+set bug_info_1 [projponf b1 p1 -min -g]
+set bug_info_1 [string trim [string range $bug_info_1 [expr {[string first "=" $bug_info_1] + 1}] [expr {[string length $bug_info_1] - 1}]]]
+set bug_info_1 [string trim [string range $bug_info_1 0 [expr {[string first " " $bug_info_1] - 1}]]]
+set bug_info_2 [projponf b1 p1 -minmax -g]
+set bug_info_2 [string trim [string range $bug_info_2 [expr {[string first "=" $bug_info_2] + 1}] [expr {[string length $bug_info_2] - 1}]]]
+set bug_info_2 [string trim [string range $bug_info_2 0 [expr {[string first " " $bug_info_2] - 1}]]]
+
+point p2 29200.099609375 7100 -2.17374753743702e-013
+set bug_info_3 [projponf b1 p2 -min -g]
+set bug_info_3 [string trim [string range $bug_info_3 [expr {[string first "=" $bug_info_3] + 1}] [expr {[string length $bug_info_3] - 1}]]]
+set bug_info_3 [string trim [string range $bug_info_3 0 [expr {[string first " " $bug_info_3] - 1}]]]
+set bug_info_4 [projponf b1 p2 -minmax -g]
+set bug_info_4 [string trim [string range $bug_info_4 [expr {[string first "=" $bug_info_4] + 1}] [expr {[string length $bug_info_4] - 1}]]]
+set bug_info_4 [string trim [string range $bug_info_4 0 [expr {[string first " " $bug_info_4] - 1}]]]
+
+if {$bug_info_1 != $bug_info_2} {
+  puts "ERROR: OCC26356 is reproduced."
+  puts "For point #1: distance min is: ${bug_info_1}, distance minmax is: ${bug_info_2}."
+}
+if {$bug_info_3 != $bug_info_4} {
+  puts "ERROR: OCC26356 is reproduced."
+  puts "For point #2: distance min is: ${bug_info_3}, distance minmax is: ${bug_info_4}."
+}