162 changes: 119 additions & 43 deletions src/plugins/fits/gaussian_weighted/fitgaussian_weighted.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
#include "objectstore.h"
#include "ui_fitgaussian_weightedconfig.h"

#define NUM_PARAMS 3
#define NUM_PARAMS 4
#define MAX_NUM_ITERATIONS 500

#include <gsl/gsl_fit.h>
Expand Down Expand Up @@ -46,13 +46,22 @@ class ConfigWidgetFitGaussianWeightedPlugin : public Kst::DataObjectConfigWidget
_vectorX->setObjectStore(store);
_vectorY->setObjectStore(store);
_vectorWeights->setObjectStore(store);

_scalarOffset->setObjectStore(store);
_scalarOffset->setDefaultValue(0.0);
_forceOffset->setChecked(false);
_scalarOffset->setEnabled(false);

}

void setupSlots(QWidget* dialog) {
if (dialog) {
connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
connect(_vectorWeights, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
connect(_scalarOffset, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
connect(_forceOffset, SIGNAL(toggled(bool)), dialog, SIGNAL(modified()));
connect(_forceOffset, SIGNAL(toggled(bool)), _scalarOffset, SLOT(setEnabled(bool)));
}
}

Expand Down Expand Up @@ -83,6 +92,8 @@ class ConfigWidgetFitGaussianWeightedPlugin : public Kst::DataObjectConfigWidget
setSelectedVectorX(source->vectorX());
setSelectedVectorY(source->vectorY());
setSelectedVectorWeights(source->vectorWeights());
_forceOffset->setChecked(source->_forceOffset);
_scalarOffset->setSelectedScalar(source->_offset);
}
}

Expand Down Expand Up @@ -145,6 +156,7 @@ class ConfigWidgetFitGaussianWeightedPlugin : public Kst::DataObjectConfigWidget

FitGaussianWeightedSource::FitGaussianWeightedSource(Kst::ObjectStore *store)
: Kst::BasicPlugin(store) {
_forceOffset = false;
}


Expand All @@ -162,6 +174,8 @@ void FitGaussianWeightedSource::change(Kst::DataObjectConfigWidget *configWidget
setInputVector(VECTOR_IN_X, config->selectedVectorX());
setInputVector(VECTOR_IN_Y, config->selectedVectorY());
setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
_forceOffset = config->_forceOffset->isChecked();
_offset = config->_scalarOffset->selectedScalar();
}
}

Expand All @@ -174,50 +188,93 @@ void FitGaussianWeightedSource::setupOutputs() {
setOutputScalar(SCALAR_OUT, "");
}

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=0;
double x_at_max_y=0;

double A, C, D;

// find peak, valley, 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];

void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates ) {
double dMin;
double dMax;
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);
}

gsl_stats_minmax( &dMin, &dMax, pdX, 1, iLength );
// 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/10 of the x range (?)

P[0] = A;
P[1] = (max_x - min_x)*0.1;
P[2] = C;
P[3] = D;

pdParameterEstimates[0] = gsl_stats_mean( pdX, 1, iLength );
pdParameterEstimates[1] = ( dMax - dMin ) / 2.0;
pdParameterEstimates[2] = gsl_stats_max( pdY, 1, iLength );
}


double function_calculate( double dX, double* pdParameters ) {
double dMean = pdParameters[0];
double dSD = pdParameters[1];
double dScale = pdParameters[2];
double dY;
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 = offset_;

dY = ( dScale / ( dSD * M_SQRT2 * M_SQRTPI ) );
dY *= exp( -( ( dX - dMean ) * ( dX - dMean ) ) / ( 2.0 * dSD * dSD ) );
if (n_params == 4) {
D = P[3];
}
x -= C;

return dY;
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;

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;

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 );

pdDerivatives[0] = ddMean;
pdDerivatives[1] = ddSD;
pdDerivatives[2] = ddScale;
}


Expand All @@ -232,6 +289,17 @@ bool FitGaussianWeightedSource::algorithm() {
Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];

if (_forceOffset) {
if (_offset) {
offset_ = _offset->value();
} else {
offset_ = 0;
}
n_params = 3;
} else {
n_params = 4;
}


Kst::LabelInfo label_info = inputVectorY->labelInfo();
label_info.name = tr("Gaussian Fit to %1").arg(label_info.name);
Expand Down Expand Up @@ -312,15 +380,21 @@ void FitGaussianWeightedSource::saveProperties(QXmlStreamWriter &s) {
QString FitGaussianWeightedSource::parameterName(int index) const {
QString parameter;
switch (index) {
case 0:
parameter = "Mean";
break;
case 1:
parameter = "SD";
break;
case 2:
parameter = "Scale";
break;
case 0:
parameter = "Amplitude";
break;
case 1:
parameter = "\\sigma";
break;
case 2:
parameter = "x_o";
break;
case 3:
parameter = "Offset";
break;
default:
parameter = "";
break;
}

return parameter;
Expand All @@ -343,6 +417,8 @@ Kst::DataObject *FitGaussianWeightedPlugin::create(Kst::ObjectStore *store, Kst:
object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
object->setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
object->_forceOffset = config->_forceOffset->isChecked();
object->_offset = config->_scalarOffset->selectedScalar();
}

object->setPluginName(pluginName());
Expand Down
2 changes: 2 additions & 0 deletions src/plugins/fits/gaussian_weighted/fitgaussian_weighted.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ class FitGaussianWeightedSource : public Kst::BasicPlugin {
virtual QStringList outputStringList() const;

virtual void saveProperties(QXmlStreamWriter &s);
bool _forceOffset;
Kst::ScalarPtr _offset;

protected:
FitGaussianWeightedSource(Kst::ObjectStore *store);
Expand Down
36 changes: 32 additions & 4 deletions src/plugins/fits/gaussian_weighted/fitgaussian_weightedconfig.ui
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
<x>0</x>
<y>0</y>
<width>426</width>
<height>114</height>
<height>106</height>
</rect>
</property>
<property name="minimumSize">
Expand Down Expand Up @@ -36,7 +36,7 @@
</property>
</widget>
</item>
<item row="0" column="1">
<item row="0" column="1" colspan="2">
<widget class="Kst::VectorSelector" name="_vectorX" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Preferred">
Expand Down Expand Up @@ -65,7 +65,7 @@
</property>
</widget>
</item>
<item row="1" column="1">
<item row="1" column="1" colspan="2">
<widget class="Kst::VectorSelector" name="_vectorY" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Preferred">
Expand All @@ -91,7 +91,7 @@
</property>
</widget>
</item>
<item row="2" column="1">
<item row="2" column="1" colspan="2">
<widget class="Kst::VectorSelector" name="_vectorWeights" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Preferred">
Expand All @@ -101,6 +101,29 @@
</property>
</widget>
</item>
<item row="3" column="0" colspan="2">
<widget class="QCheckBox" name="_forceOffset">
<property name="sizePolicy">
<sizepolicy hsizetype="Fixed" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="text">
<string>Force Offset to:</string>
</property>
</widget>
</item>
<item row="3" column="2">
<widget class="Kst::ScalarSelector" name="_scalarOffset" native="true">
<property name="minimumSize">
<size>
<width>30</width>
<height>0</height>
</size>
</property>
</widget>
</item>
</layout>
</widget>
<layoutdefault spacing="6" margin="11"/>
Expand All @@ -110,6 +133,11 @@
<extends>QWidget</extends>
<header>vectorselector.h</header>
</customwidget>
<customwidget>
<class>Kst::ScalarSelector</class>
<extends>QWidget</extends>
<header>scalarselector.h</header>
</customwidget>
</customwidgets>
<resources/>
<connections/>
Expand Down
33 changes: 22 additions & 11 deletions src/plugins/fits/non_linear.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,9 @@ struct data {
const double* pdY;
};

int n_params = NUM_PARAMS;
double offset_ = 0.0;

void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates );
double function_calculate( double dX, double* pdParameters );
void function_derivative( double dX, double* pdParameters, double* pdDerivatives );
Expand All @@ -46,7 +49,7 @@ int function_f( const gsl_vector* pVectorX, void* pParams, gsl_vector* pVectorF
data* pData = (data*)pParams;
size_t i;

for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dParameters[i] = gsl_vector_get( pVectorX, i );
}

Expand All @@ -66,14 +69,14 @@ int function_df( const gsl_vector* pVectorX, void* pParams, gsl_matrix* pMatrixJ
size_t i;
size_t j;

for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dParameters[i] = gsl_vector_get( pVectorX, i );
}

for( i=0; i<pData->n; i++ ) {
function_derivative( pData->pdX[i], dParameters, dDerivatives );

for( j=0; j<NUM_PARAMS; j++ ) {
for( j=0; j<n_params; j++ ) {
gsl_matrix_set( pMatrixJ, i, j, dDerivatives[j] );
}
}
Expand Down Expand Up @@ -153,7 +156,7 @@ bool kstfit_nonlinear(
vectorOutYCovariance->resize(NUM_PARAMS * NUM_PARAMS);

pType = gsl_multifit_fdfsolver_lmsder;
pSolver = gsl_multifit_fdfsolver_alloc( pType, iLength, NUM_PARAMS );
pSolver = gsl_multifit_fdfsolver_alloc( pType, iLength, n_params );
if( pSolver != NULL ) {
d.n = iLength;
d.pdX = pInputX;
Expand All @@ -163,13 +166,13 @@ bool kstfit_nonlinear(
function.df = function_df;
function.fdf = function_fdf;
function.n = iLength;
function.p = NUM_PARAMS;
function.p = n_params;
function.params = &d;

pMatrixCovariance = gsl_matrix_alloc( NUM_PARAMS, NUM_PARAMS );
pMatrixCovariance = gsl_matrix_alloc( n_params, n_params );
if( pMatrixCovariance != NULL ) {
function_initial_estimate( pInputX, pInputY, iLength, dXInitial );
vectorViewInitial = gsl_vector_view_array( dXInitial, NUM_PARAMS );
vectorViewInitial = gsl_vector_view_array( dXInitial, n_params );

gsl_multifit_fdfsolver_set( pSolver, &function, &vectorViewInitial.vector );

Expand All @@ -184,7 +187,7 @@ bool kstfit_nonlinear(
iIterations++;
} while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );
#if GSL_MAJOR_VERSION >= 2
pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
pMatrixJacobian = gsl_matrix_alloc( iLength, n_params );
#else
pMatrixJacobian = pSolver->J;
#endif
Expand All @@ -197,7 +200,7 @@ bool kstfit_nonlinear(
//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}

Expand All @@ -210,9 +213,17 @@ bool kstfit_nonlinear(
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
if (i<n_params) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
} else {
vectorOutYParameters->raw_V_ptr()[i] = offset_;
}
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
if ((i<n_params) && (j<n_params)) {
vectorOutYCovariance->raw_V_ptr()[(i*n_params)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
} else {
vectorOutYCovariance->raw_V_ptr()[(i*n_params)+j] = 0.0;
}
}
}

Expand Down
33 changes: 22 additions & 11 deletions src/plugins/fits/non_linear_weighted.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,9 @@ struct data {
const double* pdWeight;
};

int n_params = NUM_PARAMS;
double offset_ = 0.0;

void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates );
double function_calculate( double dX, double* pdParameters );
void function_derivative( double dX, double* pdParameters, double* pdDerivatives );
Expand All @@ -47,7 +50,7 @@ int function_f( const gsl_vector* pVectorX, void* pParams, gsl_vector* pVectorF
data* pData = (data*)pParams;
size_t i;

for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dParameters[i] = gsl_vector_get( pVectorX, i );
}

Expand All @@ -67,14 +70,14 @@ int function_df( const gsl_vector* pVectorX, void* pParams, gsl_matrix* pMatrixJ
size_t i;
size_t j;

for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dParameters[i] = gsl_vector_get( pVectorX, i );
}

for( i=0; i<pData->n; i++ ) {
function_derivative( pData->pdX[i], dParameters, dDerivatives );

for( j=0; j<NUM_PARAMS; j++ ) {
for( j=0; j<n_params; j++ ) {
gsl_matrix_set( pMatrixJ, i, j, dDerivatives[j] * pData->pdWeight[i] );
}
}
Expand Down Expand Up @@ -166,7 +169,7 @@ bool kstfit_nonlinear_weighted(
vectorOutYCovariance->resize(NUM_PARAMS * NUM_PARAMS);

pType = gsl_multifit_fdfsolver_lmsder;
pSolver = gsl_multifit_fdfsolver_alloc( pType, iLength, NUM_PARAMS );
pSolver = gsl_multifit_fdfsolver_alloc( pType, iLength, n_params );
if( pSolver != NULL ) {
d.n = iLength;
d.pdX = pInputs[XVALUES];
Expand All @@ -177,13 +180,13 @@ bool kstfit_nonlinear_weighted(
function.df = function_df;
function.fdf = function_fdf;
function.n = iLength;
function.p = NUM_PARAMS;
function.p = n_params;
function.params = &d;

pMatrixCovariance = gsl_matrix_alloc( NUM_PARAMS, NUM_PARAMS );
pMatrixCovariance = gsl_matrix_alloc( n_params, n_params );
if( pMatrixCovariance != NULL ) {
function_initial_estimate( pInputs[XVALUES], pInputs[YVALUES], iLength, dXInitial );
vectorViewInitial = gsl_vector_view_array( dXInitial, NUM_PARAMS );
vectorViewInitial = gsl_vector_view_array( dXInitial, n_params );

gsl_multifit_fdfsolver_set( pSolver, &function, &vectorViewInitial.vector );

Expand All @@ -200,7 +203,7 @@ bool kstfit_nonlinear_weighted(
while( iStatus == GSL_CONTINUE && iIterations < MAX_NUM_ITERATIONS );

#if GSL_MAJOR_VERSION >= 2
pMatrixJacobian = gsl_matrix_alloc( iLength, NUM_PARAMS );
pMatrixJacobian = gsl_matrix_alloc( iLength, n_params );
#else
pMatrixJacobian = pSolver->J;
#endif
Expand All @@ -214,7 +217,7 @@ bool kstfit_nonlinear_weighted(
//
// determine the fitted values...
//
for( i=0; i<NUM_PARAMS; i++ ) {
for( i=0; i<n_params; i++ ) {
dXInitial[i] = gsl_vector_get( pSolver->x, i );
}

Expand All @@ -227,9 +230,17 @@ bool kstfit_nonlinear_weighted(
// fill in the parameter values and covariance matrix...
//
for( i=0; i<NUM_PARAMS; i++ ) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
if (i<n_params) {
vectorOutYParameters->raw_V_ptr()[i] = gsl_vector_get( pSolver->x, i );
} else {
vectorOutYParameters->raw_V_ptr()[i] = offset_;
}
for( j=0; j<NUM_PARAMS; j++ ) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
if ((i<n_params) && (j<n_params)) {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = gsl_matrix_get( pMatrixCovariance, i, j );
} else {
vectorOutYCovariance->raw_V_ptr()[(i*NUM_PARAMS)+j] = 0.0;
}
}
}

Expand Down