1 /***************************************************************************
2 * *
3 * copyright : (C) 2007 The University of Toronto *
4 * netterfield@astro.utoronto.ca *
5 * copyright : (C) 2005 University of British Columbia *
6 * dscott@phas.ubc.ca *
7 * *
8 * This program is free software; you can redistribute it and/or modify *
9 * it under the terms of the GNU General Public License as published by *
10 * the Free Software Foundation; either version 2 of the License, or *
11 * (at your option) any later version. *
12 * *
13 ***************************************************************************/
14
15
16 #include "fitlorentzian_weighted.h"
17 #include "objectstore.h"
18 #include "ui_fitlorentzian_weightedconfig.h"
19
20 #define NUM_PARAMS 3
21 #define MAX_NUM_ITERATIONS 500
22
23 #include <gsl/gsl_fit.h>
24 #include "../non_linear_weighted.h"
25
26 static const QString& VECTOR_IN_X = "X Vector";
27 static const QString& VECTOR_IN_Y = "Y Vector";
28 static const QString& VECTOR_IN_WEIGHTS = "Weights Vector";
29 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
30 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
31 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
32 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance";
33 static const QString& SCALAR_OUT = "chi^2/nu";
34
35 class ConfigWidgetFitLorentzianWeightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitLorentzian_WeightedConfig {
36 public:
ConfigWidgetFitLorentzianWeightedPlugin(QSettings * cfg)37 ConfigWidgetFitLorentzianWeightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitLorentzian_WeightedConfig() {
38 _store = 0;
39 setupUi(this);
40 }
41
~ConfigWidgetFitLorentzianWeightedPlugin()42 ~ConfigWidgetFitLorentzianWeightedPlugin() {}
43
setObjectStore(Kst::ObjectStore * store)44 void setObjectStore(Kst::ObjectStore* store) {
45 _store = store;
46 _vectorX->setObjectStore(store);
47 _vectorY->setObjectStore(store);
48 _vectorWeights->setObjectStore(store);
49 }
50
setupSlots(QWidget * dialog)51 void setupSlots(QWidget* dialog) {
52 if (dialog) {
53 connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
54 connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
55 connect(_vectorWeights, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
56 }
57 }
58
setVectorX(Kst::VectorPtr vector)59 void setVectorX(Kst::VectorPtr vector) {
60 setSelectedVectorX(vector);
61 }
62
setVectorY(Kst::VectorPtr vector)63 void setVectorY(Kst::VectorPtr vector) {
64 setSelectedVectorY(vector);
65 }
66
setVectorsLocked(bool locked=true)67 void setVectorsLocked(bool locked = true) {
68 _vectorX->setEnabled(!locked);
69 _vectorY->setEnabled(!locked);
70 }
71
selectedVectorX()72 Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
setSelectedVectorX(Kst::VectorPtr vector)73 void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
74
selectedVectorY()75 Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
setSelectedVectorY(Kst::VectorPtr vector)76 void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
77
selectedVectorWeights()78 Kst::VectorPtr selectedVectorWeights() { return _vectorWeights->selectedVector(); };
setSelectedVectorWeights(Kst::VectorPtr vector)79 void setSelectedVectorWeights(Kst::VectorPtr vector) { return _vectorWeights->setSelectedVector(vector); };
80
setupFromObject(Kst::Object * dataObject)81 virtual void setupFromObject(Kst::Object* dataObject) {
82 if (FitLorentzianWeightedSource* source = static_cast<FitLorentzianWeightedSource*>(dataObject)) {
83 setSelectedVectorX(source->vectorX());
84 setSelectedVectorY(source->vectorY());
85 setSelectedVectorWeights(source->vectorWeights());
86 }
87 }
88
configurePropertiesFromXml(Kst::ObjectStore * store,QXmlStreamAttributes & attrs)89 virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
90 Q_UNUSED(store);
91 Q_UNUSED(attrs);
92
93 bool validTag = true;
94
95 // QStringRef av;
96 // av = attrs.value("value");
97 // if (!av.isNull()) {
98 // _configValue = QVariant(av.toString()).toBool();
99 // }
100
101 return validTag;
102 }
103
104 public slots:
save()105 virtual void save() {
106 if (_cfg) {
107 _cfg->beginGroup("Fit Lorentzian Weighted Plugin");
108 _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
109 _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
110 _cfg->setValue("Input Vector Weights", _vectorWeights->selectedVector()->Name());
111 _cfg->endGroup();
112 }
113 }
114
load()115 virtual void load() {
116 if (_cfg && _store) {
117 _cfg->beginGroup("Fit Lorentzian Weighted Plugin");
118 QString vectorName = _cfg->value("Input Vector X").toString();
119 Kst::Object* object = _store->retrieveObject(vectorName);
120 Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
121 if (vectorx) {
122 setSelectedVectorX(vectorx);
123 }
124 vectorName = _cfg->value("Input Vector Y").toString();
125 object = _store->retrieveObject(vectorName);
126 Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
127 if (vectory) {
128 setSelectedVectorX(vectory);
129 }
130 vectorName = _cfg->value("Input Vector Weights").toString();
131 object = _store->retrieveObject(vectorName);
132 Kst::Vector* vectorweights = static_cast<Kst::Vector*>(object);
133 if (vectorweights) {
134 setSelectedVectorX(vectorweights);
135 }
136 _cfg->endGroup();
137 }
138 }
139
140 private:
141 Kst::ObjectStore *_store;
142
143 };
144
145
FitLorentzianWeightedSource(Kst::ObjectStore * store)146 FitLorentzianWeightedSource::FitLorentzianWeightedSource(Kst::ObjectStore *store)
147 : Kst::BasicPlugin(store) {
148 }
149
150
~FitLorentzianWeightedSource()151 FitLorentzianWeightedSource::~FitLorentzianWeightedSource() {
152 }
153
154
_automaticDescriptiveName() const155 QString FitLorentzianWeightedSource::_automaticDescriptiveName() const {
156 return tr("%1 Weighted Lorentzian").arg(vectorY()->descriptiveName());;
157 }
158
159
change(Kst::DataObjectConfigWidget * configWidget)160 void FitLorentzianWeightedSource::change(Kst::DataObjectConfigWidget *configWidget) {
161 if (ConfigWidgetFitLorentzianWeightedPlugin* config = static_cast<ConfigWidgetFitLorentzianWeightedPlugin*>(configWidget)) {
162 setInputVector(VECTOR_IN_X, config->selectedVectorX());
163 setInputVector(VECTOR_IN_Y, config->selectedVectorY());
164 setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
165 }
166 }
167
168
setupOutputs()169 void FitLorentzianWeightedSource::setupOutputs() {
170 setOutputVector(VECTOR_OUT_Y_FITTED, "");
171 setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
172 setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
173 setOutputVector(VECTOR_OUT_Y_COVARIANCE, "");
174 setOutputScalar(SCALAR_OUT, "");
175 }
176
177
function_initial_estimate(const double * pdX,const double * pdY,int iLength,double * pdParameterEstimates)178 void function_initial_estimate( const double* pdX, const double* pdY, int iLength, double* pdParameterEstimates ) {
179 double dMin;
180 double dMax;
181
182 gsl_stats_minmax( &dMin, &dMax, pdX, 1, iLength );
183
184 pdParameterEstimates[0] = gsl_stats_mean( pdX, 1, iLength );
185 pdParameterEstimates[1] = ( dMax - dMin ) / 2.0;
186 pdParameterEstimates[2] = gsl_stats_max( pdY, 1, iLength );
187 }
188
189
function_calculate(double dX,double * pdParameters)190 double function_calculate( double dX, double* pdParameters ) {
191 double dMean = pdParameters[0];
192 double dHW = pdParameters[1];
193 double dScale = pdParameters[2];
194 double dY;
195
196 dY = ( dScale / M_PI ) * ( dHW / 2.0 );
197 dY /= ( ( dX - dMean ) * ( dX - dMean ) )+( ( dHW / 2.0 ) * ( dHW / 2.0 ) );
198
199 return dY;
200 }
201
202
function_derivative(double dX,double * pdParameters,double * pdDerivatives)203 void function_derivative( double dX, double* pdParameters, double* pdDerivatives ) {
204 double dMean = pdParameters[0];
205 double dHW = pdParameters[1];
206 double dScale = pdParameters[2];
207 double dDenom;
208 double ddMean;
209 double ddHW;
210 double ddScale;
211
212 dDenom = ( ( dX - dMean ) * ( dX - dMean ) ) + ( ( dHW / 2.0 ) * ( dHW / 2.0 ) );
213 ddMean = ( dScale / M_PI ) * dHW * ( dMean - dX ) / ( dDenom * dDenom );
214 ddHW = ( dScale / ( 2.0 * M_PI ) ) / ( dDenom * dDenom );
215 ddHW *= dDenom - ( dHW * dHW / 2.0 );
216 ddScale = ( 1.0 / M_PI ) * ( dHW / 2.0 ) / dDenom;
217
218 pdDerivatives[0] = ddMean;
219 pdDerivatives[1] = ddHW;
220 pdDerivatives[2] = ddScale;
221 }
222
223
algorithm()224 bool FitLorentzianWeightedSource::algorithm() {
225 Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
226 Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
227 Kst::VectorPtr inputVectorWeights = _inputVectors[VECTOR_IN_WEIGHTS];
228
229 Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
230 Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
231 Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
232 Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
233 Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];
234
235 Kst::LabelInfo label_info = inputVectorY->labelInfo();
236 label_info.name = tr("Lorentzian Fit to %1").arg(label_info.name);
237 outputVectorYFitted->setLabelInfo(label_info);
238
239 label_info.name = tr("Lorentzian Fit Residuals");
240 outputVectorYResiduals->setLabelInfo(label_info);
241
242
243 bool bReturn = false;
244
245 bReturn = kstfit_nonlinear_weighted( inputVectorX, inputVectorY, inputVectorWeights,
246 outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters,
247 outputVectorYCovariance, outputScalar );
248 return bReturn;
249 }
250
251
vectorX() const252 Kst::VectorPtr FitLorentzianWeightedSource::vectorX() const {
253 return _inputVectors[VECTOR_IN_X];
254 }
255
256
vectorY() const257 Kst::VectorPtr FitLorentzianWeightedSource::vectorY() const {
258 return _inputVectors[VECTOR_IN_Y];
259 }
260
261
vectorWeights() const262 Kst::VectorPtr FitLorentzianWeightedSource::vectorWeights() const {
263 return _inputVectors[VECTOR_IN_WEIGHTS];
264 }
265
266
inputVectorList() const267 QStringList FitLorentzianWeightedSource::inputVectorList() const {
268 QStringList vectors(VECTOR_IN_X);
269 vectors += VECTOR_IN_Y;
270 vectors += VECTOR_IN_WEIGHTS;
271 return vectors;
272 }
273
274
inputScalarList() const275 QStringList FitLorentzianWeightedSource::inputScalarList() const {
276 return QStringList();
277 }
278
279
inputStringList() const280 QStringList FitLorentzianWeightedSource::inputStringList() const {
281 return QStringList( /*STRING_IN*/ );
282 }
283
284
outputVectorList() const285 QStringList FitLorentzianWeightedSource::outputVectorList() const {
286 QStringList vectors(VECTOR_OUT_Y_FITTED);
287 vectors += VECTOR_OUT_Y_RESIDUALS;
288 vectors += VECTOR_OUT_Y_PARAMETERS;
289 vectors += VECTOR_OUT_Y_COVARIANCE;
290 vectors += VECTOR_OUT_Y_PARAMETERS;
291 return vectors;
292 }
293
294
outputScalarList() const295 QStringList FitLorentzianWeightedSource::outputScalarList() const {
296 return QStringList( SCALAR_OUT );
297 }
298
299
outputStringList() const300 QStringList FitLorentzianWeightedSource::outputStringList() const {
301 return QStringList( /*STRING_OUT*/ );
302 }
303
304
saveProperties(QXmlStreamWriter & s)305 void FitLorentzianWeightedSource::saveProperties(QXmlStreamWriter &s) {
306 Q_UNUSED(s);
307 // s.writeAttribute("value", _configValue);
308 }
309
310
parameterName(int index) const311 QString FitLorentzianWeightedSource::parameterName(int index) const {
312 QString parameter;
313 switch (index) {
314 case 0:
315 parameter = "Mean";
316 break;
317 case 1:
318 parameter = "Half-width";
319 break;
320 case 2:
321 parameter = "Scale";
322 break;
323 }
324
325 return parameter;
326 }
327
328
329 // Name used to identify the plugin. Used when loading the plugin.
pluginName() const330 QString FitLorentzianWeightedPlugin::pluginName() const { return tr("Lorentzian Weighted Fit"); }
pluginDescription() const331 QString FitLorentzianWeightedPlugin::pluginDescription() const { return tr("Generates a lorentzian weighted fit for a set of data."); }
332
333
create(Kst::ObjectStore * store,Kst::DataObjectConfigWidget * configWidget,bool setupInputsOutputs) const334 Kst::DataObject *FitLorentzianWeightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
335
336 if (ConfigWidgetFitLorentzianWeightedPlugin* config = static_cast<ConfigWidgetFitLorentzianWeightedPlugin*>(configWidget)) {
337
338 FitLorentzianWeightedSource* object = store->createObject<FitLorentzianWeightedSource>();
339
340 if (setupInputsOutputs) {
341 object->setupOutputs();
342 object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
343 object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
344 object->setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
345 }
346
347 object->setPluginName(pluginName());
348
349 object->writeLock();
350 object->registerChange();
351 object->unlock();
352
353 return object;
354 }
355 return 0;
356 }
357
358
configWidget(QSettings * settingsObject) const359 Kst::DataObjectConfigWidget *FitLorentzianWeightedPlugin::configWidget(QSettings *settingsObject) const {
360 ConfigWidgetFitLorentzianWeightedPlugin *widget = new ConfigWidgetFitLorentzianWeightedPlugin(settingsObject);
361 return widget;
362 }
363
364 #ifndef QT5
365 Q_EXPORT_PLUGIN2(kstplugin_FitLorentzianWeightedPlugin, FitLorentzianWeightedPlugin)
366 #endif
367
368 // vim: ts=2 sw=2 et
369