Enable boundary check on proposed value in MC. Refs #6747.

mantidproject · Mar 25, 2013 · 0eb712f · 0eb712f
1 parent 0362755
commit 0eb712f
Show file tree

Hide file tree

Showing 2 changed files with 96 additions and 8 deletions.
diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/LeBailFit.h
@@ -234,6 +234,9 @@ namespace CurveFitting
     bool proposeNewValues(vector<string> mcgroup, double m_totRwp,
                           map<string, Parameter> &curparammap, map<string, Parameter> &newparammap, bool prevBetterRwp);
 
+    ///  Limit proposed value in the specified boundary
+    double limitProposedValueInBound(Parameter param, double newvalue, double direction, int choice);
+
     /// Book keep the (sopposed) best MC result
     void bookKeepBestMCResult(map<string, Parameter> parammap,
                               vector<double> &bkgddata, double rwp, size_t istep);

diff --git a/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp b/Code/Mantid/Framework/CurveFitting/src/LeBailFit.cpp
@@ -27,6 +27,7 @@
 #define INPUTBACKGROUNDINDEX  6   // Input background
 #define SMOOTHEDBACKGROUND    8   // Smoothed background
 const double NEG_DBL_MAX(-1.*DBL_MAX);
+const double NOBOUNDARYLIMIT(1.0E10);
 
 using namespace Mantid;
 using namespace Mantid::CurveFitting;
@@ -3157,6 +3158,20 @@ namespace CurveFitting
         newvalue = fabs(newvalue);
       }
 
+      // keep the new value in the boundary
+      if (newvalue < param.minvalue)
+      {
+        int toss = rand()%2;
+        double direction = -1.0;
+        newvalue = limitProposedValueInBound(param, newvalue, direction, toss);
+      }
+      else if (newvalue > param.maxvalue)
+      {
+        int toss = rand()%2;
+        double direction = 1.0;
+        newvalue = limitProposedValueInBound(param, newvalue, direction, toss);
+      }
+
       // apply to new parameter map
       newparammap[paramname].value = newvalue;
 
@@ -3189,6 +3204,76 @@ namespace CurveFitting
     return anyparameterrefined;
   }
 
+  //-----------------------------------------------------------------------------------------------
+  /** Limit proposed value in the specified boundary
+    * @param param     :: Parameter
+    * @param newvalue  :: proposed new value that is out of boundary
+    * @param direction :: direction of parameter moved.  -1 for lower.  1 for upper
+    * @param option    :: option for various method  0: half distance.  1: periodic / reflection
+    *                     based on boundary
+    *
+    * @return :: new value in boundary
+    */
+  double LeBailFit::limitProposedValueInBound(Parameter param, double newvalue, double direction, int choice)
+  {
+    if (choice == 0)
+    {
+      // Half distance
+      if (direction > 0)
+      {
+        newvalue = (param.maxvalue - param.value)*0.5 + param.value;
+      }
+      else
+      {
+        newvalue = param.minvalue + 0.5 * (param.value - param.minvalue);
+      }
+    }
+    else
+    {
+      double deltaX = param.maxvalue-param.minvalue;
+
+      if (deltaX < NOBOUNDARYLIMIT)
+      {
+        choice = 1;  // periodic
+      }
+      else
+      {
+        choice = 2;  // reflection
+      }
+
+      if (choice == 1)
+      {
+        // Periodic boundary
+        if (direction > 0)
+        {
+          // newvalue = param.minvalue + (newvalue - param.maxvalue) % deltaX;
+          double dval = (newvalue - param.maxvalue)/deltaX;
+          newvalue = param.minvalue + deltaX * (dval - floor(dval));
+        }
+        else
+        {
+          // newvalue = param.maxvalue - (param.minvalue - newvalue) % deltaX;
+          double dval = (param.minvalue - newvalue)/deltaX;
+          newvalue = param.maxvalue - deltaX * (dval - floor(dval));
+        }
+      }
+      else
+      {
+        // Reflective boundary
+        if (direction > 0)
+        {
+          newvalue = param.maxvalue - (newvalue - param.maxvalue);
+        }
+        else
+        {
+          newvalue = param.minvalue + (param.maxvalue - newvalue);
+        }
+      }
+    }
+
+    return newvalue;
+  }
+
   //-----------------------------------------------------------------------------------------------
   /** Determine whether the proposed value should be accepted or denied
     * @param currwp:  current Rwp
@@ -3241,14 +3326,14 @@ namespace CurveFitting
       m_bestMCStep = istep;
 
       if (m_bestParameters.size() == 0)
-    {
-      // If not be initialized, initialize it!
-      m_bestParameters = parammap;
-      // copyParameterMap(parammap, m_bestParameters);
-    }
-    else
-    {
-      applyParameterValues(parammap, m_bestParameters);
+      {
+        // If not be initialized, initialize it!
+        m_bestParameters = parammap;
+        // copyParameterMap(parammap, m_bestParameters);
+      }
+      else
+      {
+        applyParameterValues(parammap, m_bestParameters);
       }
 
       m_bestBackgroundData = bkgddata;