[Kst] branches/work/kst/portto4/kst/src

[Barth Netterfield]

SVN commit 1329652 by netterfield:

Fix a crash case in plot items (uninitialized variable)
Fix unweighted gaussian fits so that they actually fit.


 M  +2 -1      libkstapp/plotitem.cpp  
 M  +79 -35    plugins/fits/gaussian_unweighted/fitgaussian_unweighted.cpp  
 M  +2 -3      plugins/fits/non_linear.h  


--- branches/work/kst/portto4/kst/src/libkstapp/plotitem.cpp #1329651:1329652
@@ -93,6 +93,7 @@
   _fitMenu(0),
   _editMenu(0),
   _sharedAxisBoxMenu(0),
+  _copyMenu(0),
   _sharedBox(0),
   _axisLabelsDirty(true),
 
@@ -4322,9 +4323,9 @@
   if (shareBox) {
     if (applyX) {
       shareBox->zoomXOut(0);
+    }
     }
   }
-}
 
 
 /*
--- branches/work/kst/portto4/kst/src/plugins/fits/gaussian_unweighted/fitgaussian_unweighted.cpp #1329651:1329652
@@ -17,7 +17,7 @@
 #include "objectstore.h"
 #include "ui_fitgaussian_unweightedconfig.h"
 
-#define NUM_PARAMS 3
+#define NUM_PARAMS 4
 #define MAX_NUM_ITERATIONS 500
 
 #include <gsl/gsl_fit.h>
@@ -160,52 +160,93 @@
 }
 
 
-void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates ) {
-  double dMin;
-  double dMax;
+void function_initial_estimate( const double X[], const double Y[], int npts, double P[] ) {
+  double min_y = 1E300;
+  double max_y = -1E300;
+  double min_x = 1E300;
+  double max_x = -1E300;
+  double mean_y = 0.0;
+  double x_at_min_y;
+  double x_at_max_y;
 
-  gsl_stats_minmax( &dMin, &dMax, pdX, 1, iLength );
+  double A, B, C, D;
 
-  pdParameterEstimates[0] = gsl_stats_mean( pdX, 1, iLength );
-  pdParameterEstimates[1] = ( dMax - dMin ) / 2.0;
-  pdParameterEstimates[2] = gsl_stats_max( pdY, 1, iLength );
+  // find peak, vally, and mean
+  for (int i = 0; i<npts; i++) {
+    if (min_y > Y[i]) {
+      min_y = Y[i];
+      x_at_min_y = X[i];
 }
+    if (max_y < Y[i]) {
+      max_y = Y[i];
+      x_at_max_y = X[i];
+    }
+    mean_y += Y[i];
 
+    if (min_x > X[i]) {
+      min_x = X[i];
+    }
+    if (max_x < X[i]) {
+      max_x = X[i];
+    }
+  }
+  if (npts>0) {
+    mean_y /= double(npts);
+  }
 
-double function_calculate( double dX, double* pdParameters ) {
-  double dMean  = pdParameters[0];
-  double dSD    = pdParameters[1];
-  double dScale = pdParameters[2];
-  double dY;
+  // Heuristic for finding the sign of the : less time is spent in the peak than
+  // in background if the range covers more than ~+- 2 sigma.
+  // It will fail if you are
+  // really zoomed into the gaussian.  Not sure what happens then :-(
+  if (max_y - mean_y > mean_y - min_y) { // positive going gaussian
+    A = max_y - min_y;
+    D = min_y;
+    C = x_at_max_y;
+  } else { // negative going gaussian
+    A = min_y - mean_y;
+    D = max_y;
+    C = x_at_min_y;
+  }
+  // guess that the width of the gaussian is around 1/4 of the x range (?)
+  B = (max_x - min_x)*0.25;
 
-  dY  = ( dScale / ( dSD * M_SQRT2 * M_SQRTPI ) );
-  dY *= exp( -( ( dX - dMean ) * ( dX - dMean ) ) / ( 2.0 * dSD * dSD ) );
+  P[0] = A;
+  P[1] = 1.0/(sqrt(2.0*B));
+  P[2] = C;
+  P[3] = D;
 
-  return dY;
 }
 
 
-void function_derivative( double dX, double* pdParameters, double* pdDerivatives ) {
-  double dMean  = pdParameters[0];
-  double dSD    = pdParameters[1];
-  double dScale = pdParameters[2];
-  double dExp;  
-  double ddMean;
-  double ddSD;
-  double ddScale;
+double function_calculate(double x, double P[]) {
+  double A = P[0];
+  double B = 0.5/(P[1]*P[1]);
+  double C = P[2];
+  double D = P[3];
 
-  dExp    = exp( -( ( dX - dMean ) * ( dX - dMean ) ) / ( 2.0 * dSD * dSD ) );
-  ddMean  = ( dX - dMean ) * dScale * dExp / ( dSD * dSD * dSD * M_SQRT2 * M_SQRTPI );
-  ddSD    = dScale * dExp / ( dSD * dSD * M_SQRT2 * M_SQRTPI );
-  ddSD   *= ( ( dX - dMean ) * ( dX - dMean ) / ( dSD * dSD ) ) - 1.0;
-  ddScale = dExp / ( dSD * M_SQRT2 * M_SQRTPI );
+  x -= C;
 
-  pdDerivatives[0] = ddMean;
-  pdDerivatives[1] = ddSD;
-  pdDerivatives[2] = ddScale;
+  return A*exp(-B*x*x) + D;
 }
 
+void function_derivative( double x, double P[], double dFdP[] ) {
+  double A = P[0];
+  double s = P[1];
+  double B = 0.5/(s*s);
+  double C = P[2];
+  //double D = P[3];
 
+  x -= C;
+
+  double E = exp(-B*x*x);
+
+  dFdP[0] = E;
+  dFdP[1] = A*x*x*E/(s*s*s);
+  dFdP[2] = 2*A*B*x*E;
+  dFdP[3] = 1.0;
+
+}
+
 bool FitGaussianUnweightedSource::algorithm() {
   Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
   Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
@@ -290,14 +331,17 @@
   QString parameter;
   switch (index) {
     case 0:
-      parameter = "Mean";
+    parameter = "Amplitude";
       break;
     case 1:
-      parameter = "SD";
+    parameter = "\\sigma";
       break;
     case 2:
-      parameter = "Scale";
+    parameter = "x_o";
       break;
+  case 3:
+    parameter = "Offset";
+    break;
   }
 
   return parameter;
--- branches/work/kst/portto4/kst/src/plugins/fits/non_linear.h #1329651:1329652
@@ -173,11 +173,10 @@
           do {
             iStatus = gsl_multifit_fdfsolver_iterate( pSolver );
             if( iStatus == GSL_SUCCESS ) {
-              iStatus = gsl_multifit_test_delta( pSolver->dx, pSolver->x, 1.0e-4, 1.0e-4 );
+              iStatus = gsl_multifit_test_delta( pSolver->dx, pSolver->x, 1.0e-6, 1.0e-6 );
             }
             iIterations++;
-          }
-          while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
+          } while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
 
           gsl_multifit_covar( pSolver->J, 0.0, pMatrixCovariance );