Add chi2/dof and enrichment text to Rel. Eff. chart

On the Rel. Eff. chart added a chi2/dof (and its p-value) to the chart. If it s a U-enrich, or Pu-enrich problem, or a ratio of just two nuclides, this relevant info will also appear on the chart. The chi2 and p-values only use stat uncert, and not the "add. uncert", so if you use additional uncertainty, these values will be terrible.
sandialabs · Oct 7, 2024 · 0db5a25 · 0db5a25
1 parent c741a1c
commit 0db5a25
Show file tree

Hide file tree

Showing 7 changed files with 184 additions and 17 deletions.
diff --git a/InterSpec/RelEffChart.h b/InterSpec/RelEffChart.h
@@ -30,12 +30,14 @@ class RelEffChart : public Wt::WContainerWidget
   void setData( const double live_time,
                const std::vector<PeakDef> &fit_peaks,
                const std::vector<RelActCalcAuto::NuclideRelAct> &rel_acts,
-               const std::string &jsRelEffEqn );
+               const std::string &jsRelEffEqn,
+               const Wt::WString &chi2_title_str );
 
   /** Set data from an "manual" relative efficiency fit. */
   void setData( const std::vector<RelActCalcManual::GenericPeakInfo> &peaks,
                const std::map<std::string,std::pair<double,std::string>> &relActsColors,
-               std::string relEffEqn );
+               std::string relEffEqn,
+               const Wt::WString &chi2_title_str );
 
   void setLineColor( const Wt::WColor &color );
   void setTextColor( const Wt::WColor &color );

diff --git a/InterSpec_resources/RelEffPlot.js b/InterSpec_resources/RelEffPlot.js
@@ -169,12 +169,17 @@ RelEffPlot.prototype.setYAxisTitle = function( title, dontCallResize ){
 }//RelEffPlot.prototype.setYAxisTitle
 
 
-RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
+RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn, chi2_txt) {
   const self = this;
 
   if( !Array.isArray(data_vals) || (data_vals.length === 0) )
     data_vals = null;
 
+    console.log( 'chi2_txt:', chi2_txt );
+  this.chi2Txt = chi2_txt;
+  if( (typeof this.chi2Txt !== "string") || (this.chi2Txt.length === 0) )
+    this.chi2Txt = null;
+
   const parentWidth = this.chart.clientWidth;
   const parentHeight = this.chart.clientHeight;
 
@@ -185,6 +190,26 @@ RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
   const yaxistitleBB = this.yaxistitle ? this.yaxistitle.node().getBBox() : null;
   const ytitleh = yaxistitleBB ? yaxistitleBB.height + titlePad : 0;
 
+  if( this.chi2Txt )
+  {
+    if( !this.chartInfoTitle ){
+      this.chartInfoTitle = this.svg.append("text")
+        .attr("class", "ChartInfoTitle")
+        .attr("y",0)                    // right at the top
+        .style("text-anchor", "end")     // Align text to the end (right)
+        .attr("dominant-baseline", "hanging") // Align text to the top of the text box
+        ;
+    }
+
+    this.chartInfoTitle.attr("x", parentWidth - 10)
+      .text(this.chi2Txt);
+  }else if( this.chartInfoTitle )
+  {
+    this.chartInfoTitle.remove();
+    this.chartInfoTitle = null;
+  }
+  const chi2_txt_pad = this.chartInfoTitle ? 0.5*this.chartInfoTitle.node().getBBox().height : 0;
+
   // We will perform an initial estimate of the chart are - using the axis we have.
   //  Later on we will put current values into the axis, and then get the size, and
   //  update things
@@ -195,7 +220,7 @@ RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
   let ytickw = yticks.empty() ? 37 : Math.max( 37, d3.max(yticks[0], function(t){ return d3.select(t).node().getBBox().width; }) ); //If we have no labels, we'll get like 7 px, so will require at least 37 px, which is the nominal value we should get
 
   let chartAreaWidth = parentWidth - this.options.margins.left - this.options.margins.right - ytitleh - ytickw;
-  let chartAreaHeight = parentHeight - this.options.margins.top - this.options.margins.bottom - xtitleh - xtickh;
+  let chartAreaHeight = parentHeight - this.options.margins.top - this.options.margins.bottom - xtitleh - xtickh - chi2_txt_pad;
 
   this.data_vals = data_vals;
   this.fit_eqn = fit_eqn;
@@ -260,7 +285,7 @@ RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
   ytickw = yticks.empty() ? 7 : d3.max(yticks[0], function(t){ return d3.select(t).node().getBBox().width; });
 
   chartAreaWidth = parentWidth - this.options.margins.left - this.options.margins.right - ytitleh - ytickw;
-  chartAreaHeight = parentHeight - this.options.margins.top - this.options.margins.bottom - xtitleh - xtickh;
+  chartAreaHeight = parentHeight - this.options.margins.top - this.options.margins.bottom - xtitleh - xtickh - chi2_txt_pad;
 
   this.xScale.range([0, chartAreaWidth]);
   this.yScale.range([chartAreaHeight, 0]);
@@ -279,13 +304,13 @@ RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
   if( this.yaxistitle )
     this.yaxistitle
       .attr("y", this.options.margins.left )
-      .attr("x",this.options.margins.top - 0.5*chartAreaHeight );
+      .attr("x",this.options.margins.top + chi2_txt_pad - 0.5*chartAreaHeight );
 
   this.chartArea.selectAll('.xAxis')
     .attr("transform", "translate(0," + chartAreaHeight + ")");
   this.chartArea
     .attr("transform", "translate(" + (this.options.margins.left + ytickw + ytitleh)
-                       + "," + this.options.margins.top + ")");
+                       + "," + (this.options.margins.top + chi2_txt_pad) + ")");
 
   if( fit_eqn ){
     // Define the line
@@ -405,5 +430,5 @@ RelEffPlot.prototype.setRelEffData = function (data_vals, fit_eqn) {
 
 
 RelEffPlot.prototype.handleResize = function () {
-  this.setRelEffData(this.data_vals, this.fit_eqn);
+  this.setRelEffData(this.data_vals, this.fit_eqn, this.chi2_txt );
 };//RelEffPlot.prototype.handleResize
diff --git a/InterSpec_resources/app_text/RelActManualGui.xml b/InterSpec_resources/app_text/RelActManualGui.xml
@@ -95,4 +95,7 @@
   <message id="ramg-consider-add-uncert-u">You may want to consider adding "Add. Uncert" for Uranium problems.</message>
   <message id="ramg-you-using-add-uncert-u">You are adding "Add. Uncert" - interpret computed uncertainties with care.</message>
   <message id="ramg-you-using-add-uncert-non-u">You are currently adding "Add. Uncert", which will cause computed uncertainties to not be correct.</message>
+  <message id="ramg-chart-info-title">χ²/dof = {1}/{2}{3}{4}</message>
+  <message id="ramg-pval"> (pval={1})</message>
+  <message id="ramg-1-pval"> (1-pval={1})</message>
 </messages>
diff --git a/src/PeakSearchGuiUtils.cpp b/src/PeakSearchGuiUtils.cpp
@@ -927,7 +927,8 @@ class PeakSelectorWindow : public AuxWindow
     m_chartPanel->setHidden( (max_num_peaks < 2) );
 
     string relEffEqn = "";
-    m_rel_eff_chart->setData( peaks, relActsColors, relEffEqn );
+    WString title_chi2_info;
+    m_rel_eff_chart->setData( peaks, relActsColors, relEffEqn, title_chi2_info );
   }//void refreshRelEffChart()
 #endif // USE_REL_ACT_TOOL
 

diff --git a/src/RelActAutoGui.cpp b/src/RelActAutoGui.cpp
@@ -4664,8 +4664,66 @@ void RelActAutoGui::updateFromCalc( std::shared_ptr<RelActCalcAuto::RelActAutoSo
                                                                   answer->m_rel_eff_coefficients );
 
   const double live_time = answer->m_foreground ? answer->m_foreground->live_time() : 1.0f;
+
+
+  WString chi2_title("χ²/dof = {1}/{2}{3}");
+  chi2_title.arg( SpecUtils::printCompact(answer->m_chi2, 3) )
+            .arg( static_cast<int>(answer->m_dof) );
+
+  // If we have U or Pu, we'll give the enrichment, or if we have two nuclides we'll
+  //  give their ratio
+  set<const SandiaDecay::Nuclide *> isotopes;
+  for( const auto &relact : answer->m_rel_activities )
+  {
+    if( relact.nuclide )
+      isotopes.insert( relact.nuclide );
+  }
+
+  const SandiaDecay::SandiaDecayDataBase *db = DecayDataBaseServer::database();
+  assert( db );
+  const SandiaDecay::Nuclide * const u235 = db->nuclide( "U235" );
+  const SandiaDecay::Nuclide * const u238 = db->nuclide( "U238" );
+  const SandiaDecay::Nuclide * const pu239 = db->nuclide( "Pu239" );
+  const SandiaDecay::Nuclide * const pu240 = db->nuclide( "Pu240" );
+  assert( u235 && u238 && pu239 && pu240 );
+
+  if( (isotopes.count(u235) && isotopes.count(u238) && !isotopes.count(pu239))
+     || (isotopes.count(pu239) && isotopes.count(pu240) && !isotopes.count(u235)) )
+  {
+    const SandiaDecay::Nuclide * const iso = isotopes.count(u235) ? u235 : pu239;
+    const double nominal = answer->mass_enrichment_fraction( iso );
+    // TODO: add errors
+    string enrich = ", " + SpecUtils::printCompact(100.0*nominal, 4)
+                        // + "±" + SpecUtils::printCompact(100.0*error, 4)
+                         + "% " + iso->symbol;
+    chi2_title.arg( enrich );
+  }else if( isotopes.size() == 2 )
+  {
+    const vector<RelActCalcAuto::NuclideRelAct> &rel_acts = answer->m_rel_activities;
+    const int num_index = (rel_acts[0].rel_activity > rel_acts[1].rel_activity) ? 1 : 0;
+    const int denom_index = (num_index ? 0 : 1);
+    const SandiaDecay::Nuclide * const num_nuc = rel_acts[num_index].nuclide;
+    const SandiaDecay::Nuclide * const den_nuc = rel_acts[denom_index].nuclide;
+    assert( num_nuc && den_nuc );
+    if( num_nuc && den_nuc )
+    {
+      const double ratio = answer->activity_ratio( num_nuc, den_nuc );
+      // TODO: add errors
+      string ratio_txt = ", act(" + num_nuc->symbol + "/" + den_nuc->symbol + ")="
+      + SpecUtils::printCompact(ratio, 4);
+
+      chi2_title.arg( ratio_txt );
+    }else
+    {
+      chi2_title.arg( "" );
+    }
+  }else
+  {
+    chi2_title.arg( "" );
+  }
 
-  m_rel_eff_chart->setData( live_time, answer->m_fit_peaks, answer->m_rel_activities, rel_eff_eqn_js );
+  m_rel_eff_chart->setData( live_time, answer->m_fit_peaks, answer->m_rel_activities,
+                           rel_eff_eqn_js, chi2_title );
 
 
   bool any_nucs_updated = false;

diff --git a/src/RelActManualGui.cpp b/src/RelActManualGui.cpp
@@ -25,6 +25,8 @@
 
 #include "rapidxml/rapidxml.hpp"
 
+#include <boost/math/distributions/chi_squared.hpp>
+
 #include <Wt/WMenu>
 #include <Wt/WLabel>
 #include <Wt/WPanel>
@@ -1098,7 +1100,7 @@ void RelActManualGui::render( Wt::WFlags<Wt::RenderFlag> flags )
 void RelActManualGui::calculateSolution()
 {
   m_currentSolution.reset();
-  m_chart->setData( vector<RelActCalcManual::GenericPeakInfo>{}, {}, "" );
+  m_chart->setData( vector<RelActCalcManual::GenericPeakInfo>{}, {}, "", {} );
   m_results->clear();
 
   // TODO: should do the actual computation not on the GUI thread!
@@ -1584,7 +1586,7 @@ void RelActManualGui::updateGuiWithResults( shared_ptr<RelActCalcManual::RelEffS
     case RelActCalcManual::ManualSolutionStatus::ErrorFindingSolution:
     case RelActCalcManual::ManualSolutionStatus::ErrorGettingSolution:
     {
-      m_chart->setData( vector<RelActCalcManual::GenericPeakInfo>{}, {}, "" );
+      m_chart->setData( vector<RelActCalcManual::GenericPeakInfo>{}, {}, "", {} );
       break;
     }
 
@@ -1640,8 +1642,79 @@ void RelActManualGui::updateGuiWithResults( shared_ptr<RelActCalcManual::RelEffS
     relActsColors[act.m_isotope] = make_pair( act.m_rel_activity, color );
   }//for( const auto &act : solution.m_rel_activities )
 
+  WString pval_str;
+  try
+  {
+    boost::math::chi_squared chi_squared_dist( solution.m_dof );
+    const double prob = boost::math::cdf( chi_squared_dist, solution.m_chi2 );
+    if( prob > 0.99 )
+      pval_str = WString::tr("ramg-1-pval").arg( SpecUtils::printCompact(1.0 - prob, 3) );
+    else
+      pval_str = WString::tr("ramg-pval").arg( SpecUtils::printCompact(prob, 3) );
+  }catch( std::exception &e )
+  {
+    pval_str = "";
+  }
+
+  WString chi2_title = WString::tr("ramg-chart-info-title");
+  chi2_title.arg( SpecUtils::printCompact(solution.m_chi2, 3) )
+            .arg( static_cast<int>(solution.m_dof) )
+            .arg( pval_str );
+
+  // If we have U or Pu, we'll give the enrichment, or if we have two nuclides we'll
+  //  give their ratio.  If we have U and Pu, we wont give enrichment.
+  set<string> isotopes;
+  for( const auto &relact : solution.m_rel_activities )
+    isotopes.insert( relact.m_isotope );
+  if( (isotopes.count("U235") && isotopes.count("U238") && !isotopes.count("Pu239"))
+     || (isotopes.count("Pu239") && isotopes.count("Pu240") && !isotopes.count("U235")) )
+  {
+    string enrich;
+    const string iso = isotopes.count("U235") ? "U235" : "Pu239";
+
+    try
+    {
+      const double nominal = solution.mass_fraction(iso);
+      const double plus = solution.mass_fraction( iso, 1.0 );
+      const double minus = solution.mass_fraction( iso, -1.0 );
+      const double error = 0.5*( fabs(plus - nominal) + fabs(nominal - minus) );
+      enrich = ", " + PhysicalUnits::printValueWithUncertainty(100.0*nominal, 100.0*error, 4)
+               + "% "+ iso;
+    }catch( std::exception & )
+    {
+      // We dont have the covariance matrix required to get mass fraction errors
+      //  We'll try leaving out the errors
+      try 
+      {
+        const double nominal = solution.mass_fraction(iso);
+        enrich = ", " + SpecUtils::printCompact(100.0*nominal, 4) + "% "+ iso;
+      }catch( std::exception & )
+      {
+        // I really dont think we should ever get here
+        assert( 0 );
+      }
+    }//try / catch
+
+    chi2_title.arg( enrich );
+  }else if( isotopes.size() == 2 )
+  {
+    const vector<RelActCalcManual::IsotopeRelativeActivity> &rel_acts = solution.m_rel_activities;
+    const int num_index = (rel_acts[0].m_rel_activity > rel_acts[1].m_rel_activity) ? 1 : 0;
+    const int denom_index = (num_index ? 0 : 1);
+    const string num_nuc = rel_acts[num_index].m_isotope;
+    const string den_nuc = rel_acts[denom_index].m_isotope;
+    const double ratio = solution.activity_ratio(num_nuc, den_nuc);
+    const double uncert = solution.activity_ratio_uncert(num_nuc, den_nuc);
+
+    string ratio_txt = ", act(" + num_nuc + "/" + den_nuc + ")="
+                   + PhysicalUnits::printValueWithUncertainty(ratio, uncert, 4);
+    chi2_title.arg( ratio_txt );
+  }else
+  {
+    chi2_title.arg( "" );
+  }
 
-  m_chart->setData( solution.m_input_peak, relActsColors, relEffEqn );
+  m_chart->setData( solution.m_input_peak, relActsColors, relEffEqn, chi2_title );
 
 
   // Now update the text

diff --git a/src/RelEffChart.cpp b/src/RelEffChart.cpp
@@ -136,7 +136,8 @@ void RelEffChart::render( Wt::WFlags<Wt::RenderFlag> flags )
 void RelEffChart::setData( const double live_time,
                           const vector<PeakDef> &fit_peaks,
                           const std::vector<RelActCalcAuto::NuclideRelAct> &rel_acts,
-                          const std::string &relEffEqn )
+                          const std::string &relEffEqn,
+                          const Wt::WString &chi2_title_str )
 {
   std::vector<RelActCalcManual::GenericPeakInfo> peaks;
   map<string,pair<double,string>> relActsColors;
@@ -189,13 +190,14 @@ void RelEffChart::setData( const double live_time,
     peaks.push_back( peak );
   }//for( const PeakDef &p : fit_peaks )
 
-  setData( peaks, relActsColors, relEffEqn );
+  setData( peaks, relActsColors, relEffEqn, chi2_title_str );
 }//void setData( std::vector<PeakDef> m_fit_peaks, std::string relEffEqn )
 
 
 void RelEffChart::setData( const std::vector<RelActCalcManual::GenericPeakInfo> &peaks,
                           const map<string,pair<double,string>> &relActsColors,
-                          string relEffEqn )
+                          string relEffEqn,
+                          const Wt::WString &chi2_title_str )
 {
   char buffer[512] = { '\0' };
 
@@ -288,7 +290,10 @@ void RelEffChart::setData( const std::vector<RelActCalcManual::GenericPeakInfo>
   if( njson_entries < 1 )
     rel_eff_plot_values.str("null");
 
-  const string js = m_jsgraph + ".setRelEffData(" + rel_eff_plot_values.str() + "," + relEffEqn + ");";
+  const string js = m_jsgraph + ".setRelEffData(" + rel_eff_plot_values.str() + ","
+                     + relEffEqn + "," 
+                      + (chi2_title_str.empty() ? string("null") : chi2_title_str.jsStringLiteral())
+  + ");";
   if( isRendered() )
     doJavaScript( js );
   else