Correcting some error estimates: preventing negative errors and stopping...

Correcting some error estimates: preventing negative errors and stopping errors from being set to zero when the Y-value is zero. refs #841

Code/Mantid/Algorithms/src/Divide.cpp

@@ -23,24 +23,15 @@ namespace Mantid
         const double& leftY = lhsY[j];
         const double& rightY = rhsY[j];
 
-        // Calculate result and store in local variable to avoid overwriting original data if
-        // output workspace is same as one of the input ones
-        const double Y = leftY/rightY;
+        //  error dividing two uncorrelated numbers, re-arrange so that you don't get infinity if leftY==0 (when rightY=0 the Y value and the result will both be infinity)
+        // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
+        // (Sa c/a)2 + (Sb c/b)2 = (Sc)2
+        // = (Sa 1/b)2 + (Sb (a/b2))2
+        // (Sc)2 = (1/b)2( (Sa)2 + (Sb a/b)2 )
+        EOut[j] = sqrt( pow(lhsE[j], 2)+pow( leftY*rhsE[j]/rightY, 2) )/rightY;
 
-        if (fabs(rightY)>1.0e-12 && fabs(Y)>1.0e-12)
-        {
-          //  gaussian errors
-          // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
-          //  So after taking proportions, squaring, summing,
-          //  and taking the square root, you get a proportional error to the product c.
-          //  Multiply that proportional error by c to get the actual standard deviation Sc.
-          const double lhsFactor = (lhsE[j]<1.0e-12|| fabs(leftY)<1.0e-12) ? 0.0 : pow((lhsE[j]/leftY),2);
-          const double rhsFactor = rhsE[j]<1.0e-12 ? 0.0 : pow((rhsE[j]/rightY),2);
-          EOut[j] = std::abs(Y) * sqrt(lhsFactor+rhsFactor);
-        }
-
-        // Now store the result
-        YOut[j] = Y;
+        // Copy the result last in case one of the input workspaces is also any output
+        YOut[j] = leftY/rightY;;
       }
     }
 
@@ -48,23 +39,17 @@ namespace Mantid
                                        const double& rhsY, const double& rhsE, MantidVec& YOut, MantidVec& EOut)
     {
       // Do the right-hand part of the error calculation just once
-      const double rhsFactor = (rhsE<1.0e-12) ? 0.0 : pow((rhsE/rhsY),2);
+      const double rhsFactor = pow(rhsE/rhsY,2);
       const int bins = lhsE.size();
       for (int j=0; j<bins; ++j)
       {
         // Get reference to input Y
         const double& leftY = lhsY[j];
-        // Calculate result into local variable
-        const double Y = leftY/rhsY;
-
-        if (fabs(Y)>1.0e-12)
-        {
-          const double lhsFactor = (lhsE[j]<1.0e-12 || fabs(leftY)<1.0e-12) ? 0.0 : pow((lhsE[j]/leftY),2);
-          EOut[j] = std::abs(Y) * sqrt(lhsFactor+rhsFactor);
-        }
 
-        // Copy result in
-        YOut[j] = Y;
+        // see comment in the function above for the error formula
+        EOut[j] = sqrt( pow(lhsE[j], 2)+pow( leftY, 2)*rhsFactor )/rhsY;
+        // Copy the result last in case one of the input workspaces is also any output
+        YOut[j] = leftY/rhsY;
       }
     }
 

Code/Mantid/Algorithms/src/Minus.cpp

@@ -26,7 +26,7 @@ namespace Mantid
     {
       std::transform(lhsY.begin(),lhsY.end(),YOut.begin(),std::bind2nd(std::minus<double>(),rhsY));
       // Only do E if non-zero, otherwise just copy
-      if (rhsE > 1.0e-12)
+      if (rhsE != 0)
         std::transform(lhsE.begin(),lhsE.end(),EOut.begin(),std::bind2nd(VectorHelper::SumGaussError<double>(),rhsE));
       else
         EOut = lhsE;

Code/Mantid/Algorithms/src/Multiply.cpp

@@ -23,48 +23,31 @@ namespace Mantid
         const double& leftY = lhsY[j];
         const double& rightY = rhsY[j];
 
-        // Calculate result and store in local variable to avoid overwriting original data if
-        // output workspace is same as one of the input ones
-        const double Y = leftY*rightY;
-
-        if (fabs(Y)>1.0e-12)
-        {
-          //  gaussian errors
-          // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
-          //  So after taking proportions, squaring, summing,
-          //  and taking the square root, you get a proportional error to the product c.
-          //  Multiply that proportional error by c to get the actual standard deviation Sc.
-          const double lhsFactor = (lhsE[j]<1.0e-12 || fabs(leftY)<1.0e-12) ? 0.0 : pow((lhsE[j]/leftY),2);
-          const double rhsFactor = (rhsE[j]<1.0e-12 || fabs(rightY)<1.0e-12) ? 0.0 : pow((rhsE[j]/rightY),2);
-          EOut[j] = std::abs(Y) * sqrt(lhsFactor+rhsFactor);
-        }
-
-        // Now store the result
-        YOut[j] = Y;
+         // error multiplying two uncorrelated numbers, re-arrange so that you don't get infinity if leftY or rightY == 0
+        // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
+        // (Sc)2 = (Sa c/a)2 + (Sb c/b)2 
+        //       = (Sa b)2 + (Sb a)2 
+        EOut[j] = sqrt(pow(lhsE[j]*rightY, 2) + pow(rhsE[j]*leftY, 2));
+
+        // Copy the result last in case one of the input workspaces is also any output
+        YOut[j] = leftY*rightY;
       }
     }
 
     void Multiply::performBinaryOperation(const MantidVec& lhsX, const MantidVec& lhsY, const MantidVec& lhsE,
                                           const double& rhsY, const double& rhsE, MantidVec& YOut, MantidVec& EOut)
     {
-      // Do the right-hand part of the error calculation just once
-      const double rhsFactor = (rhsE<1.0e-12) ? 0.0 : pow((rhsE/rhsY),2);
       const int bins = lhsE.size();
       for (int j=0; j<bins; ++j)
       {
         // Get reference to input Y
         const double& leftY = lhsY[j];
-        // Calculate result into local variable
-        const double Y = leftY*rhsY;
 
-        if (fabs(Y)>1.0e-12)
-        {
-          const double lhsFactor = (lhsE[j]<1.0e-12 || fabs(leftY)<1.0e-12) ? 0.0 : pow((lhsE[j]/leftY),2);
-          EOut[j] = std::abs(Y) * sqrt(lhsFactor+rhsFactor);
-        }
+        // see comment in the function above for the error formula
+        EOut[j] = sqrt(pow(lhsE[j]*rhsY, 2) + pow(rhsE*leftY, 2));
 
-        // Copy result in
-        YOut[j] = Y;
+        // Copy the result last in case one of the input workspaces is also any output
+        YOut[j] = leftY*rhsY;
       }
     }
 

Code/Mantid/Algorithms/src/Plus.cpp

@@ -26,7 +26,7 @@ namespace Mantid
     {
       std::transform(lhsY.begin(),lhsY.end(),YOut.begin(),std::bind2nd(std::plus<double>(),rhsY));
       // Only do E if non-zero, otherwise just copy
-      if (rhsE > 1.0e-12)
+      if (rhsE != 0)
         std::transform(lhsE.begin(),lhsE.end(),EOut.begin(),std::bind2nd(VectorHelper::SumGaussError<double>(),rhsE));
       else
         EOut = lhsE;

Code/Mantid/Algorithms/src/PolynomialCorrection.cpp

@@ -3,6 +3,7 @@
 //----------------------------------------------------------------------
 #include "MantidAlgorithms/PolynomialCorrection.h"
 #include "MantidKernel/ArrayProperty.h"
+#include <cmath>
 
 using namespace Mantid::API;
 using namespace Mantid::Kernel;
@@ -41,7 +42,7 @@ namespace Algorithms
     
     // Multiply the data and error by the correction factor
     YOut = YIn*factor;
-    EOut = EIn*factor;
+    EOut = EIn*std::abs(factor);
   }
 
 } // namespace Algorithms

Code/Mantid/Algorithms/test/QxyTest.h

@@ -72,7 +72,7 @@ public:
     TS_ASSERT_EQUALS( (*(result->getAxis(1)))(0), -0.1 )
     TS_ASSERT_DELTA( (*(result->getAxis(1)))(31), -0.038, 0.001 )
     TS_ASSERT_EQUALS( (*(result->getAxis(1)))(100), 0.1 )
-    
+
     TS_ASSERT_EQUALS( result->readX(0).size(), 101 )
     TS_ASSERT_EQUALS( result->readX(0).front(), -0.1 )
     TS_ASSERT_DELTA( result->readX(0)[64], 0.028, 0.01 )
@@ -80,9 +80,9 @@ public:
     TS_ASSERT_DIFFERS( result->readY(0).front(), result->readY(0).front() )  // NaN
     TS_ASSERT_DELTA( result->readY(26)[73], 4438798, 1 )
     TS_ASSERT_DELTA( result->readY(18)[36], 174005, 1 )
-    TS_ASSERT_EQUALS( result->readE(0).front(), 0 )
-    TS_ASSERT_EQUALS( result->readE(55)[50], 0 )
-    TS_ASSERT_EQUALS( result->readE(97).back(), 0 )
+    TS_ASSERT_DELTA( result->readE(20)[67], 0.0, 1e-10 )
+    TS_ASSERT_DELTA( result->readE(27)[70], 0, 1e-10 )
+    TS_ASSERT_DELTA( result->readE(23)[34], 0.0, 1e-10 )
     
     Mantid::API::AnalysisDataService::Instance().remove(inputWS);
     Mantid::API::AnalysisDataService::Instance().remove(outputWS);