File indexing completed on 2024-12-22 04:18:13
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_unweighted.h" 0017 #include "objectstore.h" 0018 #include "ui_fitlorentzian_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& VECTOR_OUT_Y_FITTED = "Fit"; 0029 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals"; 0030 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector"; 0031 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance"; 0032 static const QString& SCALAR_OUT = "chi^2/nu"; 0033 0034 class ConfigWidgetFitLorentzianUnweightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitLorentzian_UnweightedConfig { 0035 public: 0036 ConfigWidgetFitLorentzianUnweightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitLorentzian_UnweightedConfig() { 0037 _store = 0; 0038 setupUi(this); 0039 } 0040 0041 ~ConfigWidgetFitLorentzianUnweightedPlugin() {} 0042 0043 void setObjectStore(Kst::ObjectStore* store) { 0044 _store = store; 0045 _vectorX->setObjectStore(store); 0046 _vectorY->setObjectStore(store); 0047 } 0048 0049 void setupSlots(QWidget* dialog) { 0050 if (dialog) { 0051 connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified())); 0052 connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified())); 0053 } 0054 } 0055 0056 0057 void setVectorX(Kst::VectorPtr vector) { 0058 setSelectedVectorX(vector); 0059 } 0060 0061 void setVectorY(Kst::VectorPtr vector) { 0062 setSelectedVectorY(vector); 0063 } 0064 0065 void setVectorsLocked(bool locked = true) { 0066 _vectorX->setEnabled(!locked); 0067 _vectorY->setEnabled(!locked); 0068 } 0069 0070 Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); }; 0071 void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); }; 0072 0073 Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); }; 0074 void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); }; 0075 0076 virtual void setupFromObject(Kst::Object* dataObject) { 0077 if (FitLorentzianUnweightedSource* source = static_cast<FitLorentzianUnweightedSource*>(dataObject)) { 0078 setSelectedVectorX(source->vectorX()); 0079 setSelectedVectorY(source->vectorY()); 0080 } 0081 } 0082 0083 virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) { 0084 Q_UNUSED(store); 0085 Q_UNUSED(attrs); 0086 0087 bool validTag = true; 0088 0089 // QStringRef av; 0090 // av = attrs.value("value"); 0091 // if (!av.isNull()) { 0092 // _configValue = QVariant(av.toString()).toBool(); 0093 // } 0094 0095 return validTag; 0096 } 0097 0098 public slots: 0099 virtual void save() { 0100 if (_cfg) { 0101 _cfg->beginGroup("Fit Lorentzian Plugin"); 0102 _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name()); 0103 _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name()); 0104 _cfg->endGroup(); 0105 } 0106 } 0107 0108 virtual void load() { 0109 if (_cfg && _store) { 0110 _cfg->beginGroup("Fit Lorentzian Plugin"); 0111 QString vectorName = _cfg->value("Input Vector X").toString(); 0112 Kst::Object* object = _store->retrieveObject(vectorName); 0113 Kst::Vector* vectorx = static_cast<Kst::Vector*>(object); 0114 if (vectorx) { 0115 setSelectedVectorX(vectorx); 0116 } 0117 vectorName = _cfg->value("Input Vector Y").toString(); 0118 object = _store->retrieveObject(vectorName); 0119 Kst::Vector* vectory = static_cast<Kst::Vector*>(object); 0120 if (vectory) { 0121 setSelectedVectorX(vectory); 0122 } 0123 _cfg->endGroup(); 0124 } 0125 } 0126 0127 private: 0128 Kst::ObjectStore *_store; 0129 0130 }; 0131 0132 0133 FitLorentzianUnweightedSource::FitLorentzianUnweightedSource(Kst::ObjectStore *store) 0134 : Kst::BasicPlugin(store) { 0135 } 0136 0137 0138 FitLorentzianUnweightedSource::~FitLorentzianUnweightedSource() { 0139 } 0140 0141 0142 QString FitLorentzianUnweightedSource::_automaticDescriptiveName() const { 0143 return tr("%1 Lorentzian").arg(vectorY()->descriptiveName());; 0144 } 0145 0146 0147 void FitLorentzianUnweightedSource::change(Kst::DataObjectConfigWidget *configWidget) { 0148 if (ConfigWidgetFitLorentzianUnweightedPlugin* config = static_cast<ConfigWidgetFitLorentzianUnweightedPlugin*>(configWidget)) { 0149 setInputVector(VECTOR_IN_X, config->selectedVectorX()); 0150 setInputVector(VECTOR_IN_Y, config->selectedVectorY()); 0151 } 0152 } 0153 0154 0155 void FitLorentzianUnweightedSource::setupOutputs() { 0156 setOutputVector(VECTOR_OUT_Y_FITTED, ""); 0157 setOutputVector(VECTOR_OUT_Y_RESIDUALS, ""); 0158 setOutputVector(VECTOR_OUT_Y_PARAMETERS, ""); 0159 setOutputVector(VECTOR_OUT_Y_COVARIANCE, ""); 0160 setOutputScalar(SCALAR_OUT, ""); 0161 } 0162 0163 0164 void function_initial_estimate(const double X[], const double Y[], int npts, double P[]) { 0165 double min_y = 1E300; 0166 double max_y = -1E300; 0167 double min_x = 1E300; 0168 double max_x = -1E300; 0169 double mean_y = 0.0; 0170 double x_at_min_y=0; 0171 double x_at_max_y=0; 0172 0173 double A, B, D; 0174 0175 // find peak, vally, and mean 0176 for (int i = 0; i<npts; i++) { 0177 if (min_y > Y[i]) { 0178 min_y = Y[i]; 0179 x_at_min_y = X[i]; 0180 } 0181 if (max_y < Y[i]) { 0182 max_y = Y[i]; 0183 x_at_max_y = X[i]; 0184 } 0185 mean_y += Y[i]; 0186 0187 if (min_x > X[i]) { 0188 min_x = X[i]; 0189 } 0190 if (max_x < X[i]) { 0191 max_x = X[i]; 0192 } 0193 } 0194 if (npts>0) { 0195 mean_y /= double(npts); 0196 } 0197 0198 // Heuristic for finding the sign : less time is spent in the peak than 0199 // in background if the range covers more than ~+- 2 sigma. 0200 // It will fail if you are 0201 // really zoomed in. Not sure what happens then :-( 0202 if (max_y - mean_y > mean_y - min_y) { // positive going gaussian 0203 A = max_y - min_y; 0204 D = min_y; 0205 B = x_at_max_y; 0206 } else { // negative going gaussian 0207 A = min_y - mean_y; 0208 D = max_y; 0209 B = x_at_min_y; 0210 } 0211 0212 P[0] = A; // amplitude 0213 P[1] = B; // x0 0214 P[2] = (max_x - min_x)*0.1; // Half Width: guess 1/10 the fit range 0215 P[3] = D; // offset 0216 } 0217 0218 0219 double function_calculate(double x, double P[]) { 0220 double A = P[0]; // amplitude 0221 double B = P[1]; // x0 0222 double C = P[2]; // Half Width 0223 double D = P[3]; // offset 0224 0225 x -= B; 0226 0227 double Y = A/(1.0 + x*x/(C*C)) + D; 0228 0229 return Y; 0230 0231 } 0232 0233 void function_derivative(double x, double P[], double dFdP[]) { 0234 0235 double A = P[0]; // amplitude 0236 double B = P[1]; // x0 0237 double C = P[2]; // Half Width 0238 0239 double C2 = C*C; 0240 x -= B; 0241 double x2 = x*x; 0242 double m = (x2 + C2); 0243 0244 0245 dFdP[0] = 1.0/(1.0 + x2/C2); // dF/dA 0246 dFdP[1] = 2.0*A*C2*x/(m*m); 0247 dFdP[2] = 2.0*A*C*x2/(m*m); 0248 dFdP[3] = 1.0; 0249 } 0250 0251 bool FitLorentzianUnweightedSource::algorithm() { 0252 Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X]; 0253 Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y]; 0254 0255 Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED]; 0256 Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS]; 0257 Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS]; 0258 Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE]; 0259 Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT]; 0260 0261 0262 Kst::LabelInfo label_info = inputVectorY->labelInfo(); 0263 label_info.name = tr("Lorentzian Fit to %1").arg(label_info.name); 0264 outputVectorYFitted->setLabelInfo(label_info); 0265 0266 label_info.name = tr("Lorentzian Fit Residuals"); 0267 outputVectorYResiduals->setLabelInfo(label_info); 0268 0269 bool bReturn = false; 0270 0271 bReturn = kstfit_nonlinear( inputVectorX, inputVectorY, 0272 outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters, 0273 outputVectorYCovariance, outputScalar ); 0274 return bReturn; 0275 } 0276 0277 0278 Kst::VectorPtr FitLorentzianUnweightedSource::vectorX() const { 0279 return _inputVectors[VECTOR_IN_X]; 0280 } 0281 0282 0283 Kst::VectorPtr FitLorentzianUnweightedSource::vectorY() const { 0284 return _inputVectors[VECTOR_IN_Y]; 0285 } 0286 0287 0288 QStringList FitLorentzianUnweightedSource::inputVectorList() const { 0289 QStringList vectors(VECTOR_IN_X); 0290 vectors += VECTOR_IN_Y; 0291 return vectors; 0292 } 0293 0294 0295 QStringList FitLorentzianUnweightedSource::inputScalarList() const { 0296 return QStringList(); 0297 } 0298 0299 0300 QStringList FitLorentzianUnweightedSource::inputStringList() const { 0301 return QStringList( /*STRING_IN*/ ); 0302 } 0303 0304 0305 QStringList FitLorentzianUnweightedSource::outputVectorList() const { 0306 QStringList vectors(VECTOR_OUT_Y_FITTED); 0307 vectors += VECTOR_OUT_Y_RESIDUALS; 0308 vectors += VECTOR_OUT_Y_PARAMETERS; 0309 vectors += VECTOR_OUT_Y_COVARIANCE; 0310 vectors += VECTOR_OUT_Y_PARAMETERS; 0311 return vectors; 0312 } 0313 0314 0315 QStringList FitLorentzianUnweightedSource::outputScalarList() const { 0316 return QStringList( SCALAR_OUT ); 0317 } 0318 0319 0320 QStringList FitLorentzianUnweightedSource::outputStringList() const { 0321 return QStringList( /*STRING_OUT*/ ); 0322 } 0323 0324 0325 void FitLorentzianUnweightedSource::saveProperties(QXmlStreamWriter &s) { 0326 Q_UNUSED(s); 0327 // s.writeAttribute("value", _configValue); 0328 } 0329 0330 0331 QString FitLorentzianUnweightedSource::parameterName(int index) const { 0332 QString parameter; 0333 switch (index) { 0334 case 0: 0335 parameter = "Amplitide"; 0336 break; 0337 case 1: 0338 parameter = "x_o"; 0339 break; 0340 case 2: 0341 parameter = "Half Width"; 0342 break; 0343 case 3: 0344 parameter = "Offset"; 0345 break; 0346 } 0347 0348 return parameter; 0349 } 0350 0351 0352 // Name used to identify the plugin. Used when loading the plugin. 0353 QString FitLorentzianUnweightedPlugin::pluginName() const { return tr("Lorentzian Fit"); } 0354 QString FitLorentzianUnweightedPlugin::pluginDescription() const { return tr("Generates a lorentzian fit for a set of data."); } 0355 0356 0357 Kst::DataObject *FitLorentzianUnweightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const { 0358 0359 if (ConfigWidgetFitLorentzianUnweightedPlugin* config = static_cast<ConfigWidgetFitLorentzianUnweightedPlugin*>(configWidget)) { 0360 0361 FitLorentzianUnweightedSource* object = store->createObject<FitLorentzianUnweightedSource>(); 0362 0363 if (setupInputsOutputs) { 0364 object->setupOutputs(); 0365 object->setInputVector(VECTOR_IN_X, config->selectedVectorX()); 0366 object->setInputVector(VECTOR_IN_Y, config->selectedVectorY()); 0367 } 0368 0369 object->setPluginName(pluginName()); 0370 0371 object->writeLock(); 0372 object->registerChange(); 0373 object->unlock(); 0374 0375 return object; 0376 } 0377 return 0; 0378 } 0379 0380 0381 Kst::DataObjectConfigWidget *FitLorentzianUnweightedPlugin::configWidget(QSettings *settingsObject) const { 0382 ConfigWidgetFitLorentzianUnweightedPlugin *widget = new ConfigWidgetFitLorentzianUnweightedPlugin(settingsObject); 0383 return widget; 0384 } 0385 0386 #ifndef QT5 0387 Q_EXPORT_PLUGIN2(kstplugin_FitLorentzianUnweightedPlugin, FitLorentzianUnweightedPlugin) 0388 #endif 0389 0390 // vim: ts=2 sw=2 et