Task #1892: variable step size for dnds

CEP/DP3/AOFlagger/include/AOFlagger/gui/quality/histogrampage.h

@@ -107,6 +107,7 @@ class HistogramPage : public Gtk::HBox {
 		Gtk::Frame _functionFrame;
 		Gtk::VBox _functionBox;
 		Gtk::RadioButton _nsButton, _dndsButton;
+		Gtk::Entry _deltaSEntry;
 		
 		Gtk::Button _plotPropertiesButton, _dataExportButton;
 		

CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h

@@ -291,6 +291,25 @@ class LogHistogram : public Serializable
 			return (double) sqrtl(sqErrorSum/(xsqErrorSum * (long double) (n-2)));
 		}
 		
+		double NormalizedSlopeStdDevBySampling(double startAmplitude, double endAmplitude, double slope, double stepFactor) const
+		{
+			long double sum = 0.0;
+			unsigned long n = 0;
+			if(stepFactor <= 1.0001) stepFactor = 1.0001;
+			while(startAmplitude < endAmplitude)
+			{
+				const double stepEnd = startAmplitude * stepFactor;
+				double sampledSlope = NormalizedSlope(startAmplitude, stepEnd);
+				double sampleError = sampledSlope - slope;
+				sum += sampleError * sampleError;
+				++n;
+				
+				startAmplitude = stepEnd;
+			}
+			
+			return (double) sqrtl(sum / ((long double) n*n - n));
+		}
+		
 		double PowerLawUpperLimit(double constrainingAmplitude, double exponent, double factor) const
 		{
 			const double count = NormalizedCountAbove(constrainingAmplitude);

CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp

@@ -49,6 +49,7 @@ HistogramPage::HistogramPage() :
 	_functionFrame("Function"),
 	_nsButton("N(S)"),
 	_dndsButton("dN(S)/dS"),
+	_deltaSEntry(),
 	_plotPropertiesButton("Properties"),
 	_dataExportButton("Data"),
 	_slopeFrame("Slope"),
@@ -123,6 +124,9 @@ HistogramPage::HistogramPage() :
 	_dndsButton.signal_clicked().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot));
 	_dndsButton.set_group(group);
 	_nsButton.set_active(true);
+	_functionBox.pack_start(_deltaSEntry, Gtk::PACK_SHRINK);
+	_deltaSEntry.set_text("2");
+	_deltaSEntry.signal_activate().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot));
 	
 	_functionFrame.add(_functionBox);
 	_sideBox.pack_start(_functionFrame, Gtk::PACK_SHRINK);
@@ -341,13 +345,15 @@ void HistogramPage::plotFit(const LogHistogram &histogram, const std::string &ti
 void HistogramPage::addHistogramToPlot(const LogHistogram &histogram)
 {
 	const bool derivative = _dndsButton.get_active();
+	double deltaS = atof(_deltaSEntry.get_text().c_str());
+	if(deltaS <= 1.0001) deltaS = 1.0001;
 	for(LogHistogram::iterator i=histogram.begin();i!=histogram.end();++i)
 	{
 		if(derivative)
 		{
 			const double x = i.value();
 			const double logx = log10(x);
-			const double cslope = histogram.NormalizedSlope(x*0.5, x*2.0);
+			const double cslope = histogram.NormalizedSlope(x/deltaS, x*deltaS);
 			if(std::isfinite(logx) && std::isfinite(cslope))
 				_plot.PushDataPoint(logx, cslope);
 		} else {
@@ -463,8 +469,6 @@ void HistogramPage::onDataExportClicked()
 void HistogramPage::updateSlopeFrame(const LogHistogram &histogram)
 {
 	std::stringstream str;
-	str << "Slopes:";
-	
 	addSlopeText(str, histogram, true);
 	
 	_slopeTextView.get_buffer()->set_text(str.str());
@@ -472,6 +476,8 @@ void HistogramPage::updateSlopeFrame(const LogHistogram &histogram)
 
 void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &histogram, bool updateRange)
 {
+	double deltaS = atof(_deltaSEntry.get_text().c_str());
+	if(deltaS <= 1.0001) deltaS = 1.0001;
 	double minRange, maxRange;
 	if(_slopeAutoRangeButton.get_active())
 	{
@@ -494,11 +500,12 @@ void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &his
 		slope = histogram.NormalizedSlope(minRange, maxRange),
 		offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope),
 		error = histogram.NormalizedSlopeStdDev(minRange, maxRange, slope, offset),
+		errorB = histogram.NormalizedSlopeStdDevBySampling(minRange, maxRange, slope, deltaS),
 		upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset)),
 		lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio),
 		lowerError = fabs(lowerLimit - histogram.PowerLawLowerLimit(minRange, slope - error, pow10(offset), rfiRatio)),
 		lowerLimit2 = histogram.PowerLawLowerLimit2(minRange, slope, pow10(offset), rfiRatio);
-	str << '\n' << slope << "±" << error << "\n["
+	str << '\n' << slope << "±" << error << "\n/±" << errorB << "\n["
 		<< log10(lowerLimit) << "±" << lowerError << ';' << log10(upperLimit) << ']' << '\n'
 		<< log10(lowerLimit2);
 }