From a79d975f2370695afe1efe64aefac206f56ecdd7 Mon Sep 17 00:00:00 2001
From: Julien Schueller <schueller@phimeca.com>
Date: Thu, 7 Nov 2024 11:56:51 +0100
Subject: [PATCH] PointConditionalD: Revert ratio/unif search bounds

---
 lib/etc/openturns.conf.in                           |  1 +
 lib/src/Base/Common/ResourceMap.cxx                 |  1 +
 .../Distribution/PointConditionalDistribution.cxx   | 13 ++++++++++---
 3 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/lib/etc/openturns.conf.in b/lib/etc/openturns.conf.in
index 3177b27642..c2ba791092 100644
--- a/lib/etc/openturns.conf.in
+++ b/lib/etc/openturns.conf.in
@@ -494,6 +494,7 @@
   <DeconditionedDistribution-ContinuousDiscretizationMethod value_str="GaussProduct" />
 
   <!-- OT::PointConditionalDistribution parameters -->
+  <PointConditionalDistribution-RatioUniformBoundsFactor value_float="3.0" />
   <PointConditionalDistribution-CubaAlgorithm value_str="cuhre" />
   <PointConditionalDistribution-OptimizationAlgorithm value_str="Cobyla" />
   <PointConditionalDistribution-SmallDimension value_int="6" />
diff --git a/lib/src/Base/Common/ResourceMap.cxx b/lib/src/Base/Common/ResourceMap.cxx
index 6126a11487..29402efd26 100644
--- a/lib/src/Base/Common/ResourceMap.cxx
+++ b/lib/src/Base/Common/ResourceMap.cxx
@@ -1114,6 +1114,7 @@ void ResourceMap::loadDefaultConfiguration()
   addAsString("DeconditionedDistribution-ContinuousDiscretizationMethod", "GaussProduct");
 
   // PointConditionalDistribution parameters //
+  addAsScalar("PointConditionalDistribution-RatioUniformBoundsFactor", 3.0);
   addAsString("PointConditionalDistribution-CubaAlgorithm", "cuhre");
   addAsString("PointConditionalDistribution-OptimizationAlgorithm", "Cobyla");
   addAsUnsignedInteger("PointConditionalDistribution-SmallDimension", 6);
diff --git a/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx b/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx
index 070343b937..ec9f1fba9e 100644
--- a/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx
+++ b/lib/src/Uncertainty/Distribution/PointConditionalDistribution.cxx
@@ -235,6 +235,7 @@ void PointConditionalDistribution::update()
 
   // enable simplified path
   useSimplifiedVersion_ = hasSimplifiedVersion(simplifiedVersion_);
+  LOGDEBUG(OSS() << "useSimplifiedVersion_=" << useSimplifiedVersion_);
 
   // cache marginal for reuse
   if (!useSimplifiedVersion_)
@@ -257,10 +258,12 @@ void PointConditionalDistribution::update()
   const UnsignedInteger dimension = getDimension();
   if (!useSimplifiedVersion_ && isContinuous() && (dimension <= ResourceMap::GetAsUnsignedInteger("PointConditionalDistribution-SmallDimension")))
   {
-    const Interval bounds(getRange());
+    const Point stddev(getStandardDeviation());
+    const Point mean(getMean());
+    const Scalar factor = ResourceMap::GetAsScalar("PointConditionalDistribution-RatioUniformBoundsFactor");
+    const Interval bounds(getRange().intersect(Interval(mean - factor * stddev, mean + factor * stddev)));
     const Point lb(bounds.getLowerBound());
     const Point ub(bounds.getUpperBound());
-    const Point middle(0.5 * (lb + ub));
 
     // First, the upper bound on U
     const Function objectiveU(new PointConditionalDistributionUBoundEvaluation(*this, r_));
@@ -269,9 +272,10 @@ void PointConditionalDistribution::update()
     problemU.setBounds(bounds);
     OptimizationAlgorithm algo(OptimizationAlgorithm::GetByName(ResourceMap::GetAsString("PointConditionalDistribution-OptimizationAlgorithm")));
     algo.setProblem(problemU);
-    algo.setStartingPoint(middle);
+    algo.setStartingPoint(mean);
     algo.run();
     supU_ = std::exp(algo.getResult().getOptimalValue()[0]);
+    LOGDEBUG(OSS() << "supU_=" << supU_);
 
     // Second, the lower and upper bounds on V
     const Function objectiveV(new PointConditionalDistributionVBoundEvaluation(*this, r_));
@@ -290,6 +294,7 @@ void PointConditionalDistribution::update()
         algo.setStartingPoint(ub * 0.5);
         algo.run();
         supV_[i] = std::exp(algo.getResult().getOptimalValue()[0]);
+        LOGDEBUG(OSS() << "supV_[" << i << "]=" << supV_[i]);
       }
       if (lb[i] < 0.0)
       {
@@ -298,6 +303,7 @@ void PointConditionalDistribution::update()
         algo.setStartingPoint(lb * 0.5);
         algo.run();
         infV_[i] = -std::exp(algo.getResult().getOptimalValue()[0]);
+        LOGDEBUG(OSS() << "infV_[" << i << "]=" << infV_[i]);
       }
     }
   }
@@ -312,6 +318,7 @@ void PointConditionalDistribution::update()
   x.add(getRange().getLowerBound());
   conditioningCDF_ = reorderedDistribution_.computeSequentialConditionalCDF(x);
   conditioningCDF_.resize(conditioningIndices_.getSize());
+  LOGDEBUG(OSS() << "conditioningCDF_=" << conditioningCDF_);
 }