fits/lorentzian_weighted/fitlorentzian_weighted.cpp

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 "fitlorentzian_weighted.h"
0017 #include "objectstore.h"
0018 #include "ui_fitlorentzian_weightedconfig.h"
0019
0020 #define NUM_PARAMS 3
0021 #define MAX_NUM_ITERATIONS 500
0022
0023 #include <gsl/gsl_fit.h>
0024 #include "../non_linear_weighted.h"
0025
0026 static const QString& VECTOR_IN_X = "X Vector";
0027 static const QString& VECTOR_IN_Y = "Y Vector";
0028 static const QString& VECTOR_IN_WEIGHTS = "Weights Vector";
0029 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
0030 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
0031 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
0032 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance";
0033 static const QString& SCALAR_OUT = "chi^2/nu";
0034
0035 class ConfigWidgetFitLorentzianWeightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitLorentzian_WeightedConfig {
0036   public:
0037     ConfigWidgetFitLorentzianWeightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitLorentzian_WeightedConfig() {
0038       _store = 0;
0039       setupUi(this);
0040     }
0041
0042     ~ConfigWidgetFitLorentzianWeightedPlugin() {}
0043
0044     void setObjectStore(Kst::ObjectStore* store) {
0045       _store = store;
0046       _vectorX->setObjectStore(store);
0047       _vectorY->setObjectStore(store);
0048       _vectorWeights->setObjectStore(store);
0049     }
0050
0051     void setupSlots(QWidget* dialog) {
0052       if (dialog) {
0053         connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0054         connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0055         connect(_vectorWeights, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0056       }
0057     }
0058
0059     void setVectorX(Kst::VectorPtr vector) {
0060       setSelectedVectorX(vector);
0061     }
0062
0063     void setVectorY(Kst::VectorPtr vector) {
0064       setSelectedVectorY(vector);
0065     }
0066
0067     void setVectorsLocked(bool locked = true) {
0068       _vectorX->setEnabled(!locked);
0069       _vectorY->setEnabled(!locked);
0070     }
0071
0072     Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
0073     void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
0074
0075     Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
0076     void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
0077
0078     Kst::VectorPtr selectedVectorWeights() { return _vectorWeights->selectedVector(); };
0079     void setSelectedVectorWeights(Kst::VectorPtr vector) { return _vectorWeights->setSelectedVector(vector); };
0080
0081     virtual void setupFromObject(Kst::Object* dataObject) {
0082       if (FitLorentzianWeightedSource* source = static_cast<FitLorentzianWeightedSource*>(dataObject)) {
0083         setSelectedVectorX(source->vectorX());
0084         setSelectedVectorY(source->vectorY());
0085         setSelectedVectorWeights(source->vectorWeights());
0086       }
0087     }
0088
0089     virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
0090       Q_UNUSED(store);
0091       Q_UNUSED(attrs);
0092
0093       bool validTag = true;
0094
0095 //       QStringRef av;
0096 //       av = attrs.value("value");
0097 //       if (!av.isNull()) {
0098 //         _configValue = QVariant(av.toString()).toBool();
0099 //       }
0100
0101       return validTag;
0102     }
0103
0104   public slots:
0105     virtual void save() {
0106       if (_cfg) {
0107         _cfg->beginGroup("Fit Lorentzian Weighted Plugin");
0108         _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
0109         _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
0110         _cfg->setValue("Input Vector Weights", _vectorWeights->selectedVector()->Name());
0111         _cfg->endGroup();
0112       }
0113     }
0114
0115     virtual void load() {
0116       if (_cfg && _store) {
0117         _cfg->beginGroup("Fit Lorentzian Weighted Plugin");
0118         QString vectorName = _cfg->value("Input Vector X").toString();
0119         Kst::Object* object = _store->retrieveObject(vectorName);
0120         Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
0121         if (vectorx) {
0122           setSelectedVectorX(vectorx);
0123         }
0124         vectorName = _cfg->value("Input Vector Y").toString();
0125         object = _store->retrieveObject(vectorName);
0126         Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
0127         if (vectory) {
0128           setSelectedVectorX(vectory);
0129         }
0130         vectorName = _cfg->value("Input Vector Weights").toString();
0131         object = _store->retrieveObject(vectorName);
0132         Kst::Vector* vectorweights = static_cast<Kst::Vector*>(object);
0133         if (vectorweights) {
0134           setSelectedVectorX(vectorweights);
0135         }
0136         _cfg->endGroup();
0137       }
0138     }
0139
0140   private:
0141     Kst::ObjectStore *_store;
0142
0143 };
0144
0145
0146 FitLorentzianWeightedSource::FitLorentzianWeightedSource(Kst::ObjectStore *store)
0147 : Kst::BasicPlugin(store) {
0148 }
0149
0150
0151 FitLorentzianWeightedSource::~FitLorentzianWeightedSource() {
0152 }
0153
0154
0155 QString FitLorentzianWeightedSource::_automaticDescriptiveName() const {
0156   return tr("%1 Weighted Lorentzian").arg(vectorY()->descriptiveName());;
0157 }
0158
0159
0160 void FitLorentzianWeightedSource::change(Kst::DataObjectConfigWidget *configWidget) {
0161   if (ConfigWidgetFitLorentzianWeightedPlugin* config = static_cast<ConfigWidgetFitLorentzianWeightedPlugin*>(configWidget)) {
0162     setInputVector(VECTOR_IN_X, config->selectedVectorX());
0163     setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0164     setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
0165   }
0166 }
0167
0168
0169 void FitLorentzianWeightedSource::setupOutputs() {
0170   setOutputVector(VECTOR_OUT_Y_FITTED, "");
0171   setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
0172   setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
0173   setOutputVector(VECTOR_OUT_Y_COVARIANCE, "");
0174   setOutputScalar(SCALAR_OUT, "");
0175 }
0176
0177
0178 void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates ) {
0179   double dMin;
0180   double dMax;
0181
0182   gsl_stats_minmax( &dMin, &dMax, pdX, 1, iLength );
0183
0184   pdParameterEstimates[0] = gsl_stats_mean( pdX, 1, iLength );
0185   pdParameterEstimates[1] = ( dMax - dMin ) / 2.0;
0186   pdParameterEstimates[2] = gsl_stats_max( pdY, 1, iLength );
0187 }
0188
0189
0190 double function_calculate( double dX, double* pdParameters ) {
0191   double dMean  = pdParameters[0];
0192   double dHW    = pdParameters[1];
0193   double dScale = pdParameters[2];
0194   double dY;
0195
0196   dY  = ( dScale / M_PI ) * ( dHW / 2.0 );
0197   dY /= ( ( dX - dMean ) * ( dX - dMean ) )+( ( dHW / 2.0 ) * ( dHW / 2.0 ) );
0198
0199   return dY;
0200 }
0201
0202
0203 void function_derivative( double dX, double* pdParameters, double* pdDerivatives ) {
0204   double dMean  = pdParameters[0];
0205   double dHW    = pdParameters[1];
0206   double dScale = pdParameters[2];
0207   double dDenom;
0208   double ddMean;
0209   double ddHW;
0210   double ddScale;
0211
0212   dDenom  = ( ( dX - dMean ) * ( dX - dMean ) ) + ( ( dHW / 2.0 ) * ( dHW / 2.0 ) );
0213   ddMean  = ( dScale / M_PI ) * dHW * ( dMean - dX ) / ( dDenom * dDenom );
0214   ddHW    = ( dScale / ( 2.0 * M_PI ) ) / ( dDenom * dDenom );
0215   ddHW   *= dDenom - ( dHW * dHW / 2.0 );
0216   ddScale = ( 1.0 / M_PI ) * ( dHW / 2.0 ) / dDenom;
0217
0218   pdDerivatives[0] = ddMean;
0219   pdDerivatives[1] = ddHW;
0220   pdDerivatives[2] = ddScale;
0221 }
0222
0223
0224 bool FitLorentzianWeightedSource::algorithm() {
0225   Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
0226   Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
0227   Kst::VectorPtr inputVectorWeights = _inputVectors[VECTOR_IN_WEIGHTS];
0228
0229   Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
0230   Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
0231   Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
0232   Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
0233   Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];
0234
0235   Kst::LabelInfo label_info = inputVectorY->labelInfo();
0236   label_info.name = tr("Lorentzian Fit to %1").arg(label_info.name);
0237   outputVectorYFitted->setLabelInfo(label_info);
0238
0239   label_info.name = tr("Lorentzian Fit Residuals");
0240   outputVectorYResiduals->setLabelInfo(label_info);
0241
0242
0243   bool bReturn = false;
0244
0245   bReturn = kstfit_nonlinear_weighted( inputVectorX, inputVectorY, inputVectorWeights,
0246                               outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters,
0247                               outputVectorYCovariance, outputScalar );
0248   return bReturn;
0249 }
0250
0251
0252 Kst::VectorPtr FitLorentzianWeightedSource::vectorX() const {
0253   return _inputVectors[VECTOR_IN_X];
0254 }
0255
0256
0257 Kst::VectorPtr FitLorentzianWeightedSource::vectorY() const {
0258   return _inputVectors[VECTOR_IN_Y];
0259 }
0260
0261
0262 Kst::VectorPtr FitLorentzianWeightedSource::vectorWeights() const {
0263   return _inputVectors[VECTOR_IN_WEIGHTS];
0264 }
0265
0266
0267 QStringList FitLorentzianWeightedSource::inputVectorList() const {
0268   QStringList vectors(VECTOR_IN_X);
0269   vectors += VECTOR_IN_Y;
0270   vectors += VECTOR_IN_WEIGHTS;
0271   return vectors;
0272 }
0273
0274
0275 QStringList FitLorentzianWeightedSource::inputScalarList() const {
0276   return QStringList();
0277 }
0278
0279
0280 QStringList FitLorentzianWeightedSource::inputStringList() const {
0281   return QStringList( /*STRING_IN*/ );
0282 }
0283
0284
0285 QStringList FitLorentzianWeightedSource::outputVectorList() const {
0286   QStringList vectors(VECTOR_OUT_Y_FITTED);
0287   vectors += VECTOR_OUT_Y_RESIDUALS;
0288   vectors += VECTOR_OUT_Y_PARAMETERS;
0289   vectors += VECTOR_OUT_Y_COVARIANCE;
0290   vectors += VECTOR_OUT_Y_PARAMETERS;
0291   return vectors;
0292 }
0293
0294
0295 QStringList FitLorentzianWeightedSource::outputScalarList() const {
0296   return QStringList( SCALAR_OUT );
0297 }
0298
0299
0300 QStringList FitLorentzianWeightedSource::outputStringList() const {
0301   return QStringList( /*STRING_OUT*/ );
0302 }
0303
0304
0305 void FitLorentzianWeightedSource::saveProperties(QXmlStreamWriter &s) {
0306   Q_UNUSED(s);
0307 //   s.writeAttribute("value", _configValue);
0308 }
0309
0310
0311 QString FitLorentzianWeightedSource::parameterName(int index) const {
0312   QString parameter;
0313   switch (index) {
0314     case 0:
0315       parameter = "Mean";
0316       break;
0317     case 1:
0318       parameter = "Half-width";
0319       break;
0320     case 2:
0321       parameter = "Scale";
0322       break;
0323   }
0324
0325   return parameter;
0326 }
0327
0328
0329 // Name used to identify the plugin.  Used when loading the plugin.
0330 QString FitLorentzianWeightedPlugin::pluginName() const { return tr("Lorentzian Weighted Fit"); }
0331 QString FitLorentzianWeightedPlugin::pluginDescription() const { return tr("Generates a lorentzian weighted fit for a set of data."); }
0332
0333
0334 Kst::DataObject *FitLorentzianWeightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
0335
0336   if (ConfigWidgetFitLorentzianWeightedPlugin* config = static_cast<ConfigWidgetFitLorentzianWeightedPlugin*>(configWidget)) {
0337
0338     FitLorentzianWeightedSource* object = store->createObject<FitLorentzianWeightedSource>();
0339
0340     if (setupInputsOutputs) {
0341       object->setupOutputs();
0342       object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
0343       object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0344       object->setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
0345     }
0346
0347     object->setPluginName(pluginName());
0348
0349     object->writeLock();
0350     object->registerChange();
0351     object->unlock();
0352
0353     return object;
0354   }
0355   return 0;
0356 }
0357
0358
0359 Kst::DataObjectConfigWidget *FitLorentzianWeightedPlugin::configWidget(QSettings *settingsObject) const {
0360   ConfigWidgetFitLorentzianWeightedPlugin *widget = new ConfigWidgetFitLorentzianWeightedPlugin(settingsObject);
0361   return widget;
0362 }
0363
0364 #ifndef QT5
0365 Q_EXPORT_PLUGIN2(kstplugin_FitLorentzianWeightedPlugin, FitLorentzianWeightedPlugin)
0366 #endif
0367
0368 // vim: ts=2 sw=2 et