Add erfc implementation extracted from boost.

Author: wcjohns

InterSpec/PeakDists.h

@@ -43,16 +43,19 @@ namespace PeakDists
    20191230: wcjohns extracted the boost::math::erf() function implementation
    from boost 1.65.1 for double precision (53 bit mantissa) into this function,
    boost_erf_imp(). Removing some of the supporting code structure, and
-   explicitly writing out the polynomial equation evaluation seems to speed
-   things up by about a factor of ~3 over calling boost::math::erf().
+   explicitly writing out the polynomial equation evaluation, and removing some
+   very minor corrections, seems to speed things up by about a factor of ~3 over
+   calling boost::math::erf() (and agrees with boosts implementation to within
+   1E-10%, across the entire range).
    
    In the implementation file, there is a commented out erf_approx() function looks to be about 25% faster than
    this boost version, but I havent carefully checked out the precision implications
    so not switching to it yet.
    */
   double boost_erf_imp( double z );
 
-  /** Returns `1-boost_erf_imp(z)` */
+  /** Returns `1-boost_erf_imp(z)`, however, due to rounding, is a separate implementation than erf.
+   */
   double boost_erfc_imp( double z );
 
 

src/D3SpectrumDisplayDiv.cpp

@@ -1850,7 +1850,7 @@ void D3SpectrumDisplayDiv::performExistingRoiEdgeDragWork(
          && !m_last_being_drug_peaks.empty()
          && (dt < std::chrono::seconds(15)) )
       {
-        cout << "Using previously drug-n-fit peaks" << endl;
+        //cout << "Using previously drug-n-fit peaks" << endl;
 
         const auto old_continuum = m_last_being_drug_peaks.front()->continuum();
         assert( old_continuum );

src/PeakDists.cpp

@@ -172,6 +172,121 @@ namespace PeakDists
   }//double boost_erf_imp( double z )
 
 
+  /** 20231123: The JavaScript needed a better `erfc` implementation than just `1-erf`, or otherwise there were huge
+   artifacts, and regions that just wouldnt draw, so wcjohns extracted this `erfc` function from boost; boost implements `erf`
+   and `erfc` in the same function, but since the JS just needs `erfc`, wcjohns separated the erf and erf implementation, even
+   though they are largely duplicate - saves some small amount of code/overhead.
+   */
+  double boost_erfc_imp( double z )
+  {
+    /* Since this function is a (slight) modification of boost source code, it
+     is subject to the Boost Version 1 Software License, whose text is:
+     
+     Boost Software License - Version 1.0 - August 17th, 2003
+     
+     Permission is hereby granted, free of charge, to any person or organization
+     obtaining a copy of the software and accompanying documentation covered by
+     this license (the "Software") to use, reproduce, display, distribute,
+     execute, and transmit the Software, and to prepare derivative works of the
+     Software, and to permit third-parties to whom the Software is furnished to
+     do so, all subject to the following:
+     
+     The copyright notices in the Software and this entire statement, including
+     the above license grant, this restriction and the following disclaimer,
+     must be included in all copies of the Software, in whole or in part, and
+     all derivative works of the Software, unless such copies or derivative
+     works are solely in the form of machine-executable object code generated by
+     a source language processor.
+     
+     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+     FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+     SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+     FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+     ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+     DEALINGS IN THE SOFTWARE.
+     */
+
+    if( z < 0 )
+      return (z < -0.5) ? (2 - boost_erfc_imp( -z )) : (1 + boost_erf_imp( -z ));
+    
+    if( z < 0.5 )
+      return 1 - boost_erf_imp( z );
+    
+    if( z >= 28 )
+      return 0;
+    
+    if(z < 1.5f)
+    {
+      static const double P[] = { -0.098090592216281240205, 0.178114665841120341155,
+        0.191003695796775433986, 0.0888900368967884466578, 0.0195049001251218801359,
+        0.00180424538297014223957
+      };
+      static const double Q[] = { 1.0, 1.84759070983002217845,
+        1.42628004845511324508, 0.578052804889902404909, 0.12385097467900864233,
+        0.0113385233577001411017, 0.337511472483094676155e-5
+      };
+      
+      const double zarg = z - 0.5;
+      const double P_eval = ((((zarg*P[5] + P[4])*zarg + P[3])*zarg + P[2])*zarg + P[1])*zarg + P[0];
+      const double Q_eval = (((((Q[6]*zarg + Q[5])*zarg + Q[4])*zarg + Q[3])*zarg + Q[2])*zarg + Q[1])*zarg + Q[0];
+      
+      return (0.405935764312744140625f + P_eval / Q_eval) * (exp(-z * z) / z);
+    }//if(z < 1.5f)
+    
+    if( z < 2.5f )
+    {
+      static const double P[] = { -0.0243500476207698441272, 0.0386540375035707201728,
+        0.04394818964209516296, 0.0175679436311802092299, 0.00323962406290842133584,
+        0.000235839115596880717416
+      };
+      static const double Q[] = { 1.0, 1.53991494948552447182, 0.982403709157920235114,
+        0.325732924782444448493, 0.0563921837420478160373, 0.00410369723978904575884
+      };
+      
+      const double zarg = z - 1.5;
+      const double P_eval = ((((zarg*P[5] + P[4])*zarg + P[3])*zarg + P[2])*zarg + P[1])*zarg + P[0];
+      const double Q_eval = ((((zarg*Q[5] + Q[4])*zarg + Q[3])*zarg + Q[2])*zarg + Q[1])*zarg + Q[0];
+      
+      return (0.50672817230224609375f + P_eval / Q_eval) * (exp(-z * z) / z);
+      // Boost implementation has an additional minor error correction here
+    }//if( z < 2.5f
+    
+    if( z < 4.5f )
+    {
+      static const double P[] = { 0.00295276716530971662634, 0.0137384425896355332126,
+        0.00840807615555585383007, 0.00212825620914618649141, 0.000250269961544794627958,
+        0.113212406648847561139e-4
+      };
+      static const double Q[] = { 1.0, 1.04217814166938418171,
+        0.442597659481563127003, 0.0958492726301061423444, 0.0105982906484876531489,
+        0.000479411269521714493907
+      };
+      
+      const double zarg = z - 3.5;
+      const double P_eval = ((((zarg*P[5] + P[4])*zarg + P[3])*zarg + P[2])*zarg + P[1])*zarg + P[0];
+      const double Q_eval = ((((zarg*Q[5] + Q[4])*zarg + Q[3])*zarg + Q[2])*zarg + Q[1])*zarg + Q[0];
+      return (0.5405750274658203125f + P_eval / Q_eval) * (exp(-z * z) / z);
+      // Boost implementation has an additional minor error correction here
+    }//if( z < 4.5f )
+    
+    static const double P[] = { 0.00628057170626964891937, 0.0175389834052493308818,
+      -0.212652252872804219852, -0.687717681153649930619, -2.5518551727311523996,
+      -3.22729451764143718517, -2.8175401114513378771
+    };
+    static const double Q[] = { 1.0, 2.79257750980575282228, 11.0567237927800161565,
+      15.930646027911794143, 22.9367376522880577224, 13.5064170191802889145,
+      5.48409182238641741584
+    };
+    
+    const double zarg = 1.0 / z;
+    const double P_eval = (((((P[6]*zarg + P[5])*zarg + P[4])*zarg + P[3])*zarg + P[2])*zarg + P[1])*zarg + P[0];
+    const double Q_eval = (((((Q[6]*zarg + Q[5])*zarg + Q[4])*zarg + Q[3])*zarg + Q[2])*zarg + Q[1])*zarg + Q[0];
+    
+    return (0.5579090118408203125f + P_eval / Q_eval) * (exp(-z * z) / z);
+    // Boost implementation has an additional minor error correction here
+}//double boost_erfc_imp( double z )
+  
   /*
   double erf_approx( double x )
   {
@@ -186,11 +301,7 @@ namespace PeakDists
   }
    */
 
-  double boost_erfc_imp( double z )
-  {
-    return 1.0 - boost_erf_imp( z );
-  }
-   
+
   void photopeak_function_integral( const double mean,
                                           const double sigma,
                                           const double amp,
@@ -263,7 +374,6 @@ namespace PeakDists
   double bortel_indefinite_integral( const double x, const double mean,
                            const double sigma, const double skew )
   {
-    const double root_two = boost::math::constants::root_two<double>();
     const double one_div_root_two = boost::math::constants::one_div_root_two<double>();
 
     const double t = (x - mean) / sigma;
@@ -493,7 +603,7 @@ namespace PeakDists
           return t_eqn*sigma + mean;
         }
 
-        const double root_two = boost::math::constants::root_two<double>();
+        const double root_two = boost::math::constants::root_two<double>(); //1.414213562373095048801688724209698078
         const double t_eqn = root_two* boost::math::erf_inv( (1.0/gauss_indef_amp) * (prob - indefinite_of_tail + gaus_indef_at_skew) );
         return t_eqn*sigma + mean;
       };//auto x_from_eqn

src/PeakFit.cpp

@@ -1979,10 +1979,10 @@ PeakShrdVec refitPeaksThatShareROI( const std::shared_ptr<const Measurement> &da
             refit_peaks.push_back( make_shared<PeakDef>(p) );
 
           const double refit_chi2Dof = chi2_for_region( refit_peaks, dataH, lower_channel, upper_channel ) / nbin;
-          cout << "refit_chi2Dof=" << refit_chi2Dof << ", prechi2Dof=" << prechi2Dof << ", postchi2Dof=" << postchi2Dof << endl;
+          //cout << "refit_chi2Dof=" << refit_chi2Dof << ", prechi2Dof=" << prechi2Dof << ", postchi2Dof=" << postchi2Dof << endl;
           if( refit_chi2Dof <= prechi2Dof )
           {
-            cout << "Using re-fit peaks!" << endl;
+            //cout << "Using re-fit peaks!" << endl;
             answer = refit_peaks;
           }
         }//if( output_peak.size() == inpeaks.size() )

target/testing/test_PeakDists.cpp

@@ -31,6 +31,8 @@
 //#include <boost/test/unit_test.hpp>
 #include <boost/test/included/unit_test.hpp>
 
+#include <boost/math/special_functions/erf.hpp>
+
 #include "SpecUtils/DateTime.h"
 
 #include "InterSpec/PeakDists.h"
@@ -42,6 +44,24 @@ using namespace PeakDists;
 using namespace boost::unit_test;
 
 
+BOOST_AUTO_TEST_CASE( Erfc )
+{
+  for( double x = -10; x <= 10; x += 0.01 )
+  {
+    const double our_val = PeakDists::boost_erf_imp(x);
+    const double boost_val = boost::math::erf(x);
+    BOOST_CHECK_CLOSE( our_val, boost_val, 1.0E-12 );
+  }
+  
+  for( double z = -30; z <= 30; z += 0.01 )
+  {
+    const double our_val = PeakDists::boost_erfc_imp(z);
+    const double boost_val = boost::math::erfc(z);
+    BOOST_CHECK_CLOSE( our_val, boost_val, 1.0E-10 );
+  }
+}//BOOST_AUTO_TEST_CASE( Erfc )
+
+
 BOOST_AUTO_TEST_CASE( GaussianDist )
 {
   // Check Gaussian has unit area