File indexing completed on 2024-12-22 04:18:11

0001 /***************************************************************************
0002  *                                                                         *
0003  *   copyright : (C) 2007 The University of Toronto                        *
0004  *                   netterfield@astro.utoronto.ca                         *
0005  *   copyright : (C) 2005  University of British Columbia                  *
0006  *                   dscott@phas.ubc.ca                                    *
0007  *                                                                         *
0008  *   This program is free software; you can redistribute it and/or modify  *
0009  *   it under the terms of the GNU General Public License as published by  *
0010  *   the Free Software Foundation; either version 2 of the License, or     *
0011  *   (at your option) any later version.                                   *
0012  *                                                                         *
0013  ***************************************************************************/
0014 
0015 
0016 #include "fitgaussian_unweighted.h"
0017 #include "objectstore.h"
0018 #include "ui_fitgaussian_unweightedconfig.h"
0019 
0020 #define NUM_PARAMS 4
0021 #define MAX_NUM_ITERATIONS 500
0022 
0023 #include <gsl/gsl_fit.h>
0024 #include "../non_linear.h"
0025 
0026 static const QString& VECTOR_IN_X = "X Vector";
0027 static const QString& VECTOR_IN_Y = "Y Vector";
0028 static const QString& SCALAR_IN_OFFSET = "Offset";
0029 
0030 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
0031 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
0032 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
0033 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance";
0034 static const QString& SCALAR_OUT = "chi^2/nu";
0035 
0036 class ConfigWidgetFitGaussianUnweightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitGaussian_UnweightedConfig {
0037   public:
0038     ConfigWidgetFitGaussianUnweightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitGaussian_UnweightedConfig() {
0039       _store = 0;
0040       setupUi(this);
0041     }
0042 
0043     ~ConfigWidgetFitGaussianUnweightedPlugin() {}
0044 
0045     void setObjectStore(Kst::ObjectStore* store) { 
0046       _store = store; 
0047       _vectorX->setObjectStore(store);
0048       _vectorY->setObjectStore(store);
0049       _scalarOffset->setObjectStore(store);
0050       _forceOffset->setChecked(false);
0051       _scalarOffset->setEnabled(false);
0052     }
0053 
0054     void setupSlots(QWidget* dialog) {
0055       if (dialog) {
0056         connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0057         connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0058         connect(_scalarOffset, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0059         _scalarOffset->setDefaultValue(0.0);
0060         connect(_forceOffset, SIGNAL(toggled(bool)), dialog, SIGNAL(modified()));
0061         connect(_forceOffset, SIGNAL(toggled(bool)), _scalarOffset, SLOT(setEnabled(bool)));
0062       }
0063     }
0064 
0065 
0066     void setVectorX(Kst::VectorPtr vector) {
0067       setSelectedVectorX(vector);
0068     }
0069 
0070     void setVectorY(Kst::VectorPtr vector) {
0071       setSelectedVectorY(vector);
0072     }
0073 
0074     void setVectorsLocked(bool locked = true) {
0075       _vectorX->setEnabled(!locked);
0076       _vectorY->setEnabled(!locked);
0077     }
0078 
0079     Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
0080     void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
0081 
0082     Kst::ScalarPtr scalarOffset() {return _scalarOffset->selectedScalar(); };
0083     void setScalarOffset(Kst::ScalarPtr scalar) {_scalarOffset->setSelectedScalar(scalar);};
0084 
0085     Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
0086     void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
0087 
0088     virtual void setupFromObject(Kst::Object* dataObject) {
0089       if (FitGaussianUnweightedSource* source = static_cast<FitGaussianUnweightedSource*>(dataObject)) {
0090         setSelectedVectorX(source->vectorX());
0091         setSelectedVectorY(source->vectorY());
0092         _forceOffset->setChecked(source->_forceOffset);
0093         setScalarOffset(source->scalarOffset());
0094       }
0095     }
0096 
0097     virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
0098       bool validTag = true;
0099 
0100       setObjectStore(store);
0101 
0102       QStringRef av;
0103       bool force_offset = false;
0104       av = attrs.value("ForceOffset");
0105       if (!av.isNull()) {
0106         force_offset = QVariant(av.toString()).toBool();
0107       }
0108       _forceOffset->setChecked(force_offset);
0109 //      if (force_offset) {
0110 //   The tickets will also be taken care of.     av = attrs.value("Offset");
0111 //        if (!av.isNull()) {
0112 //          QString name = av.toString();
0113 //          Kst::ObjectPtr object = store->retrieveObject(name);
0114 //          Kst::ScalarPtr scalar = Kst::kst_cast<Kst::Scalar>(object);
0115 //          setScalarOffset(scalar);
0116 //        }
0117 
0118 //      }
0119 
0120       return validTag;
0121     }
0122 
0123   public slots:
0124     virtual void save() {
0125       if (_cfg) {
0126         _cfg->beginGroup("Fit Gaussian Plugin");
0127         _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
0128         _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
0129         _cfg->setValue("Force Offset", _forceOffset->isChecked());
0130         if (_forceOffset->isChecked()) {
0131           _cfg->setValue("Offset", _scalarOffset->selectedScalar()->Name());
0132         }
0133         _cfg->endGroup();
0134       }
0135     }
0136 
0137     virtual void load() {
0138       if (_cfg && _store) {
0139         _cfg->beginGroup("Fit Gaussian Plugin");
0140         QString vectorName = _cfg->value("Input Vector X").toString();
0141         Kst::Object* object = _store->retrieveObject(vectorName);
0142         Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
0143         if (vectorx) {
0144           setSelectedVectorX(vectorx);
0145         }
0146         vectorName = _cfg->value("Input Vector Y").toString();
0147         object = _store->retrieveObject(vectorName);
0148         Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
0149         if (vectory) {
0150           setSelectedVectorX(vectory);
0151         }
0152         bool force_offset = _cfg->value("Force Offset").toBool();
0153         _forceOffset->setChecked(force_offset);
0154         if (force_offset) {
0155           QString scalarName = _cfg->value("Offset").toString();
0156           object = _store->retrieveObject(scalarName);
0157           Kst::Scalar* scalar = static_cast<Kst::Scalar*>(object);
0158           if (scalar) {
0159             setScalarOffset(scalar);
0160           }
0161         }
0162         _cfg->endGroup();
0163       }
0164     }
0165 
0166   private:
0167     Kst::ObjectStore *_store;
0168 
0169 };
0170 
0171 
0172 FitGaussianUnweightedSource::FitGaussianUnweightedSource(Kst::ObjectStore *store)
0173 : Kst::BasicPlugin(store) {
0174   _forceOffset = false;
0175 }
0176 
0177 
0178 FitGaussianUnweightedSource::~FitGaussianUnweightedSource() {
0179 }
0180 
0181 
0182 QString FitGaussianUnweightedSource::_automaticDescriptiveName() const {
0183     return tr("%1 Unweighted Gaussian").arg(vectorY()->descriptiveName());
0184 }
0185 
0186 
0187 void FitGaussianUnweightedSource::change(Kst::DataObjectConfigWidget *configWidget) {
0188   if (ConfigWidgetFitGaussianUnweightedPlugin* config = static_cast<ConfigWidgetFitGaussianUnweightedPlugin*>(configWidget)) {
0189     setInputVector(VECTOR_IN_X, config->selectedVectorX());
0190     setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0191     setInputScalar(SCALAR_IN_OFFSET, config->scalarOffset());
0192     _forceOffset = config->_forceOffset->isChecked();
0193 //    _offset = config->_scalarOffset->selectedScalar();
0194   }
0195 }
0196 
0197 
0198 void FitGaussianUnweightedSource::setupOutputs() {
0199   setOutputVector(VECTOR_OUT_Y_FITTED, "");
0200   setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
0201   setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
0202   setOutputVector(VECTOR_OUT_Y_COVARIANCE, "");
0203   setOutputScalar(SCALAR_OUT, "");
0204 
0205   int i=0;
0206   for (QString paramName = parameterName(i); !paramName.isEmpty(); paramName = parameterName(++i)) {
0207     setOutputScalar(paramName, "");
0208   }
0209 }
0210 
0211 
0212 void function_initial_estimate( const double X[], const double Y[], int npts, double P[] ) {
0213   double min_y = 1E300;
0214   double max_y = -1E300;
0215   double min_x = 1E300;
0216   double max_x = -1E300;
0217   double mean_y = 0.0;
0218   double x_at_min_y=0;
0219   double x_at_max_y=0;
0220 
0221   double A, C, D;
0222 
0223   // find peak, valley, and mean
0224   for (int i = 0; i<npts; i++) {
0225     if (min_y > Y[i]) {
0226       min_y = Y[i];
0227       x_at_min_y = X[i];
0228     }
0229     if (max_y < Y[i]) {
0230       max_y = Y[i];
0231       x_at_max_y = X[i];
0232     }
0233     mean_y += Y[i];
0234 
0235     if (min_x > X[i]) {
0236       min_x = X[i];
0237     }
0238     if (max_x < X[i]) {
0239       max_x = X[i];
0240     }
0241   }
0242   if (npts>0) {
0243     mean_y /= double(npts);
0244   }
0245 
0246   // Heuristic for finding the sign of the : less time is spent in the peak than
0247   // in background if the range covers more than ~+- 2 sigma.
0248   // It will fail if you are
0249   // really zoomed into the gaussian.  Not sure what happens then :-(
0250   if (max_y - mean_y > mean_y - min_y) { // positive going gaussian
0251     A = max_y - min_y;
0252     D = min_y;
0253     C = x_at_max_y;
0254   } else { // negative going gaussian
0255     A = min_y - mean_y;
0256     D = max_y;
0257     C = x_at_min_y;
0258   }
0259   // guess that the width of the gaussian is around 1/10 of the x range (?)
0260 
0261   P[0] = A;
0262   P[1] = (max_x - min_x)*0.1;
0263   P[2] = C;
0264   P[3] = D;
0265 
0266 }
0267 
0268 
0269 double function_calculate(double x, double P[]) {
0270   double A = P[0];
0271   double B = 0.5/(P[1]*P[1]);
0272   double C = P[2];
0273   double D = offset_;
0274 
0275   if (n_params == 4) {
0276     D = P[3];
0277   }
0278   x -= C;
0279 
0280   return A*exp(-B*x*x) + D;
0281 }
0282 
0283 void function_derivative( double x, double P[], double dFdP[] ) {
0284   double A = P[0];
0285   double s = P[1];
0286   double B = 0.5/(s*s);
0287   double C = P[2];
0288   //double D = P[3];
0289 
0290   x -= C;
0291 
0292   double E = exp(-B*x*x);
0293 
0294   dFdP[0] = E;
0295   dFdP[1] = A*x*x*E/(s*s*s);
0296   dFdP[2] = 2*A*B*x*E;
0297   dFdP[3] = 1.0;
0298 
0299 }
0300 
0301 bool FitGaussianUnweightedSource::algorithm() {
0302   Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
0303   Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
0304   Kst::ScalarPtr inputScalarOffset = _inputScalars[SCALAR_IN_OFFSET];
0305 
0306   Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
0307   Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
0308   Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
0309   Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
0310   Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];
0311 
0312 
0313   if (_forceOffset) {
0314     if (inputScalarOffset) {
0315       offset_ = inputScalarOffset->value();
0316     } else {
0317       offset_ = 0;
0318     }
0319     n_params = 3;
0320   } else {
0321     n_params = 4;
0322   }
0323 
0324   Kst::LabelInfo label_info = inputVectorY->labelInfo();
0325   label_info.name = tr("A\\cdotexp((x-x_o)^2/2\\sigma^2) + C fit to %1").arg(label_info.name);
0326   outputVectorYFitted->setLabelInfo(label_info);
0327 
0328   label_info.name = tr("Gaussian Fit Residuals");
0329   outputVectorYResiduals->setLabelInfo(label_info);
0330 
0331   bool bReturn = false;
0332 
0333   bReturn = kstfit_nonlinear( inputVectorX, inputVectorY,
0334                               outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters,
0335                               outputVectorYCovariance, outputScalar );
0336   return bReturn;
0337 }
0338 
0339 
0340 Kst::VectorPtr FitGaussianUnweightedSource::vectorX() const {
0341   return _inputVectors[VECTOR_IN_X];
0342 }
0343 
0344 
0345 Kst::VectorPtr FitGaussianUnweightedSource::vectorY() const {
0346   return _inputVectors[VECTOR_IN_Y];
0347 }
0348 
0349 
0350 Kst::ScalarPtr FitGaussianUnweightedSource::scalarOffset() const {
0351   return _inputScalars[SCALAR_IN_OFFSET];
0352 }
0353 
0354 
0355 QStringList FitGaussianUnweightedSource::inputVectorList() const {
0356   QStringList vectors(VECTOR_IN_X);
0357   vectors += VECTOR_IN_Y;
0358   return vectors;
0359 }
0360 
0361 
0362 QStringList FitGaussianUnweightedSource::inputScalarList() const {
0363   QStringList scalars(SCALAR_IN_OFFSET);
0364   return scalars;
0365 }
0366 
0367 
0368 QStringList FitGaussianUnweightedSource::inputStringList() const {
0369   return QStringList( /*STRING_IN*/ );
0370 }
0371 
0372 
0373 QStringList FitGaussianUnweightedSource::outputVectorList() const {
0374   QStringList vectors(VECTOR_OUT_Y_FITTED);
0375   vectors += VECTOR_OUT_Y_RESIDUALS;
0376   vectors += VECTOR_OUT_Y_PARAMETERS;
0377   vectors += VECTOR_OUT_Y_COVARIANCE;
0378   vectors += VECTOR_OUT_Y_PARAMETERS;
0379   return vectors;
0380 }
0381 
0382 
0383 QStringList FitGaussianUnweightedSource::outputScalarList() const {
0384   return QStringList( SCALAR_OUT );
0385 }
0386 
0387 
0388 QStringList FitGaussianUnweightedSource::outputStringList() const {
0389   return QStringList( /*STRING_OUT*/ );
0390 }
0391 
0392 
0393 void FitGaussianUnweightedSource::saveProperties(QXmlStreamWriter &s) {
0394   QString force_offset;
0395   force_offset.setNum(_forceOffset);
0396   s.writeAttribute("ForceOffset", force_offset);
0397 }
0398 
0399 
0400 QString FitGaussianUnweightedSource::parameterName(int index) const {
0401   QString parameter;
0402   switch (index) {
0403   case 0:
0404     parameter = "A";
0405     break;
0406   case 1:
0407     parameter = "\\sigma";
0408     break;
0409   case 2:
0410     parameter = "x_o";
0411     break;
0412   case 3:
0413     parameter = "C";
0414     break;
0415   default:
0416     parameter = "";
0417   }
0418 
0419   return parameter;
0420 }
0421 
0422 
0423 // Name used to identify the plugin.  Used when loading the plugin.
0424 QString FitGaussianUnweightedPlugin::pluginName() const { return tr("Gaussian Fit"); }
0425 QString FitGaussianUnweightedPlugin::pluginDescription() const { return tr("Generates a gaussian fit for a set of data."); }
0426 
0427 
0428 Kst::DataObject *FitGaussianUnweightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
0429 
0430   if (ConfigWidgetFitGaussianUnweightedPlugin* config = static_cast<ConfigWidgetFitGaussianUnweightedPlugin*>(configWidget)) {
0431 
0432     FitGaussianUnweightedSource* object = store->createObject<FitGaussianUnweightedSource>();
0433 
0434     object->_forceOffset = config->_forceOffset->isChecked();
0435 
0436     if (setupInputsOutputs) {
0437       object->setupOutputs();
0438       object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
0439       object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0440       object->setInputScalar(SCALAR_IN_OFFSET, config->scalarOffset());
0441     }
0442 
0443     object->setPluginName(pluginName());
0444 
0445     object->writeLock();
0446     object->registerChange();
0447     object->unlock();
0448 
0449     return object;
0450   }
0451   return 0;
0452 }
0453 
0454 
0455 Kst::DataObjectConfigWidget *FitGaussianUnweightedPlugin::configWidget(QSettings *settingsObject) const {
0456   ConfigWidgetFitGaussianUnweightedPlugin *widget = new ConfigWidgetFitGaussianUnweightedPlugin(settingsObject);
0457   return widget;
0458 }
0459 
0460 #ifndef QT5
0461 Q_EXPORT_PLUGIN2(kstplugin_FitGaussianUnweightedPlugin, FitGaussianUnweightedPlugin)
0462 #endif
0463 
0464 // vim: ts=2 sw=2 et