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 "fitgaussian_weighted.h"
17 #include "objectstore.h"
18 #include "ui_fitgaussian_weightedconfig.h"
19
20 #define NUM_PARAMS 4
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& SCALAR_IN_OFFSET = "Offset";
30
31 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
32 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
33 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
34 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance";
35 static const QString& SCALAR_OUT = "chi^2/nu";
36
37 class ConfigWidgetFitGaussianWeightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitGaussian_WeightedConfig {
38 public:
ConfigWidgetFitGaussianWeightedPlugin(QSettings * cfg)39 ConfigWidgetFitGaussianWeightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitGaussian_WeightedConfig() {
40 _store = 0;
41 setupUi(this);
42 }
43
~ConfigWidgetFitGaussianWeightedPlugin()44 ~ConfigWidgetFitGaussianWeightedPlugin() {}
45
setObjectStore(Kst::ObjectStore * store)46 void setObjectStore(Kst::ObjectStore* store) {
47 _store = store;
48 _vectorX->setObjectStore(store);
49 _vectorY->setObjectStore(store);
50 _vectorWeights->setObjectStore(store);
51
52 _scalarOffset->setObjectStore(store);
53 _forceOffset->setChecked(false);
54 _scalarOffset->setEnabled(false);
55
56 }
57
setupSlots(QWidget * dialog)58 void setupSlots(QWidget* dialog) {
59 if (dialog) {
60 connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
61 connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
62 connect(_vectorWeights, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
63 connect(_scalarOffset, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
64 _scalarOffset->setDefaultValue(0.0);
65 connect(_forceOffset, SIGNAL(toggled(bool)), dialog, SIGNAL(modified()));
66 connect(_forceOffset, SIGNAL(toggled(bool)), _scalarOffset, SLOT(setEnabled(bool)));
67 }
68 }
69
setVectorX(Kst::VectorPtr vector)70 void setVectorX(Kst::VectorPtr vector) {
71 setSelectedVectorX(vector);
72 }
73
setVectorY(Kst::VectorPtr vector)74 void setVectorY(Kst::VectorPtr vector) {
75 setSelectedVectorY(vector);
76 }
77
scalarOffset()78 Kst::ScalarPtr scalarOffset() {return _scalarOffset->selectedScalar(); };
setScalarOffset(Kst::ScalarPtr scalar)79 void setScalarOffset(Kst::ScalarPtr scalar) {_scalarOffset->setSelectedScalar(scalar);};
80
setVectorsLocked(bool locked=true)81 void setVectorsLocked(bool locked = true) {
82 _vectorX->setEnabled(!locked);
83 _vectorY->setEnabled(!locked);
84 }
85
selectedVectorX()86 Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
setSelectedVectorX(Kst::VectorPtr vector)87 void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
88
selectedVectorY()89 Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
setSelectedVectorY(Kst::VectorPtr vector)90 void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
91
selectedVectorWeights()92 Kst::VectorPtr selectedVectorWeights() { return _vectorWeights->selectedVector(); };
setSelectedVectorWeights(Kst::VectorPtr vector)93 void setSelectedVectorWeights(Kst::VectorPtr vector) { return _vectorWeights->setSelectedVector(vector); };
94
setupFromObject(Kst::Object * dataObject)95 virtual void setupFromObject(Kst::Object* dataObject) {
96 if (FitGaussianWeightedSource* source = static_cast<FitGaussianWeightedSource*>(dataObject)) {
97 setSelectedVectorX(source->vectorX());
98 setSelectedVectorY(source->vectorY());
99 setSelectedVectorWeights(source->vectorWeights());
100 _forceOffset->setChecked(source->_forceOffset);
101 setScalarOffset(source->scalarOffset());
102 }
103 }
104
configurePropertiesFromXml(Kst::ObjectStore * store,QXmlStreamAttributes & attrs)105 virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
106
107 bool validTag = true;
108
109 setObjectStore(store);
110
111 QStringRef av;
112 bool force_offset = false;
113 av = attrs.value("ForceOffset");
114 if (!av.isNull()) {
115 force_offset = QVariant(av.toString()).toBool();
116 }
117 _forceOffset->setChecked(force_offset);
118 // if (force_offset) {
119 // av = attrs.value("Offset");
120 // if (!av.isNull()) {
121 // QString name = av.toString();
122 // Kst::ObjectPtr object = store->retrieveObject(name);
123 // Kst::ScalarPtr scalar = Kst::kst_cast<Kst::Scalar>(object);
124 // setScalarOffset(scalar);
125 // }
126 // }
127
128 return validTag;
129 }
130
131 public slots:
save()132 virtual void save() {
133 if (_cfg) {
134 _cfg->beginGroup("Fit Gaussian Weighted Plugin");
135 _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
136 _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
137 _cfg->setValue("Input Vector Weights", _vectorWeights->selectedVector()->Name());
138 _cfg->setValue("Force Offset", _forceOffset->isChecked());
139 if (_forceOffset->isChecked()) {
140 _cfg->setValue("Offset", _scalarOffset->selectedScalar()->Name());
141 }
142 _cfg->endGroup();
143 }
144 }
145
load()146 virtual void load() {
147 if (_cfg && _store) {
148 _cfg->beginGroup("Fit Gaussian Weighted Plugin");
149 QString vectorName = _cfg->value("Input Vector X").toString();
150 Kst::Object* object = _store->retrieveObject(vectorName);
151 Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
152 if (vectorx) {
153 setSelectedVectorX(vectorx);
154 }
155 vectorName = _cfg->value("Input Vector Y").toString();
156 object = _store->retrieveObject(vectorName);
157 Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
158 if (vectory) {
159 setSelectedVectorX(vectory);
160 }
161 vectorName = _cfg->value("Input Vector Weights").toString();
162 object = _store->retrieveObject(vectorName);
163 Kst::Vector* vectorweights = static_cast<Kst::Vector*>(object);
164 if (vectorweights) {
165 setSelectedVectorX(vectorweights);
166 }
167 bool force_offset = _cfg->value("Force Offset").toBool();
168 _forceOffset->setChecked(force_offset);
169 if (force_offset) {
170 QString scalarName = _cfg->value("Offset").toString();
171 object = _store->retrieveObject(scalarName);
172 Kst::Scalar* scalar = static_cast<Kst::Scalar*>(object);
173 if (scalar) {
174 setScalarOffset(scalar);
175 }
176 }
177
178 _cfg->endGroup();
179 }
180 }
181
182 private:
183 Kst::ObjectStore *_store;
184
185 };
186
187
FitGaussianWeightedSource(Kst::ObjectStore * store)188 FitGaussianWeightedSource::FitGaussianWeightedSource(Kst::ObjectStore *store)
189 : Kst::BasicPlugin(store) {
190 _forceOffset = false;
191 }
192
193
~FitGaussianWeightedSource()194 FitGaussianWeightedSource::~FitGaussianWeightedSource() {
195 }
196
197
_automaticDescriptiveName() const198 QString FitGaussianWeightedSource::_automaticDescriptiveName() const {
199 return tr("%1 Weighted Gaussian").arg(vectorY()->descriptiveName());
200 }
201
202
change(Kst::DataObjectConfigWidget * configWidget)203 void FitGaussianWeightedSource::change(Kst::DataObjectConfigWidget *configWidget) {
204 if (ConfigWidgetFitGaussianWeightedPlugin* config = static_cast<ConfigWidgetFitGaussianWeightedPlugin*>(configWidget)) {
205 setInputVector(VECTOR_IN_X, config->selectedVectorX());
206 setInputVector(VECTOR_IN_Y, config->selectedVectorY());
207 setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
208 setInputScalar(SCALAR_IN_OFFSET, config->scalarOffset());
209 _forceOffset = config->_forceOffset->isChecked();
210 }
211 }
212
213
setupOutputs()214 void FitGaussianWeightedSource::setupOutputs() {
215 setOutputVector(VECTOR_OUT_Y_FITTED, "");
216 setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
217 setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
218 setOutputVector(VECTOR_OUT_Y_COVARIANCE, "");
219 setOutputScalar(SCALAR_OUT, "");
220
221 int i=0;
222 for (QString paramName = parameterName(i); !paramName.isEmpty(); paramName = parameterName(++i)) {
223 setOutputScalar(paramName, "");
224 }
225 }
226
function_initial_estimate(const double X[],const double Y[],int npts,double P[])227 void function_initial_estimate( const double X[], const double Y[], int npts, double P[] ) {
228 double min_y = 1E300;
229 double max_y = -1E300;
230 double min_x = 1E300;
231 double max_x = -1E300;
232 double mean_y = 0.0;
233 double x_at_min_y=0;
234 double x_at_max_y=0;
235
236 double A, C, D;
237
238 // find peak, valley, and mean
239 for (int i = 0; i<npts; i++) {
240 if (min_y > Y[i]) {
241 min_y = Y[i];
242 x_at_min_y = X[i];
243 }
244 if (max_y < Y[i]) {
245 max_y = Y[i];
246 x_at_max_y = X[i];
247 }
248 mean_y += Y[i];
249
250 if (min_x > X[i]) {
251 min_x = X[i];
252 }
253 if (max_x < X[i]) {
254 max_x = X[i];
255 }
256 }
257 if (npts>0) {
258 mean_y /= double(npts);
259 }
260
261 // Heuristic for finding the sign of the : less time is spent in the peak than
262 // in background if the range covers more than ~+- 2 sigma.
263 // It will fail if you are
264 // really zoomed into the gaussian. Not sure what happens then :-(
265 if (max_y - mean_y > mean_y - min_y) { // positive going gaussian
266 A = max_y - min_y;
267 D = min_y;
268 C = x_at_max_y;
269 } else { // negative going gaussian
270 A = min_y - mean_y;
271 D = max_y;
272 C = x_at_min_y;
273 }
274 // guess that the width of the gaussian is around 1/10 of the x range (?)
275
276 P[0] = A;
277 P[1] = (max_x - min_x)*0.1;
278 P[2] = C;
279 P[3] = D;
280
281 }
282
283
function_calculate(double x,double P[])284 double function_calculate(double x, double P[]) {
285 double A = P[0];
286 double B = 0.5/(P[1]*P[1]);
287 double C = P[2];
288 double D = offset_;
289
290 if (n_params == 4) {
291 D = P[3];
292 }
293 x -= C;
294
295 return A*exp(-B*x*x) + D;
296 }
297
function_derivative(double x,double P[],double dFdP[])298 void function_derivative( double x, double P[], double dFdP[] ) {
299 double A = P[0];
300 double s = P[1];
301 double B = 0.5/(s*s);
302 double C = P[2];
303 //double D = P[3];
304
305 x -= C;
306
307 double E = exp(-B*x*x);
308
309 dFdP[0] = E;
310 dFdP[1] = A*x*x*E/(s*s*s);
311 dFdP[2] = 2*A*B*x*E;
312 dFdP[3] = 1.0;
313
314 }
315
316
algorithm()317 bool FitGaussianWeightedSource::algorithm() {
318 Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
319 Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
320 Kst::VectorPtr inputVectorWeights = _inputVectors[VECTOR_IN_WEIGHTS];
321 Kst::ScalarPtr inputScalarOffset = _inputScalars[SCALAR_IN_OFFSET];
322
323 Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
324 Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
325 Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
326 Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
327 Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];
328
329 if (_forceOffset) {
330 if (inputScalarOffset) {
331 offset_ = inputScalarOffset->value();
332 } else {
333 offset_ = 0;
334 }
335 n_params = 3;
336 } else {
337 n_params = 4;
338 }
339
340
341 Kst::LabelInfo label_info = inputVectorY->labelInfo();
342 label_info.name = tr("A\\cdotexp((x-x_o)^2/2\\sigma^2) + C fit to %1").arg(label_info.name);
343 outputVectorYFitted->setLabelInfo(label_info);
344
345 label_info.name = tr("Gaussian Fit Residuals");
346 outputVectorYResiduals->setLabelInfo(label_info);
347
348
349 bool bReturn = false;
350
351 bReturn = kstfit_nonlinear_weighted( inputVectorX, inputVectorY, inputVectorWeights,
352 outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters,
353 outputVectorYCovariance, outputScalar );
354 return bReturn;
355 }
356
357
vectorX() const358 Kst::VectorPtr FitGaussianWeightedSource::vectorX() const {
359 return _inputVectors[VECTOR_IN_X];
360 }
361
362
vectorY() const363 Kst::VectorPtr FitGaussianWeightedSource::vectorY() const {
364 return _inputVectors[VECTOR_IN_Y];
365 }
366
367
vectorWeights() const368 Kst::VectorPtr FitGaussianWeightedSource::vectorWeights() const {
369 return _inputVectors[VECTOR_IN_WEIGHTS];
370 }
371
scalarOffset() const372 Kst::ScalarPtr FitGaussianWeightedSource::scalarOffset() const {
373 return _inputScalars[SCALAR_IN_OFFSET];
374 }
375
inputVectorList() const376 QStringList FitGaussianWeightedSource::inputVectorList() const {
377 QStringList vectors(VECTOR_IN_X);
378 vectors += VECTOR_IN_Y;
379 vectors += VECTOR_IN_WEIGHTS;
380 return vectors;
381 }
382
383
inputScalarList() const384 QStringList FitGaussianWeightedSource::inputScalarList() const {
385 QStringList scalars(SCALAR_IN_OFFSET);
386 return scalars;
387 }
388
389
inputStringList() const390 QStringList FitGaussianWeightedSource::inputStringList() const {
391 return QStringList( /*STRING_IN*/ );
392 }
393
394
outputVectorList() const395 QStringList FitGaussianWeightedSource::outputVectorList() const {
396 QStringList vectors(VECTOR_OUT_Y_FITTED);
397 vectors += VECTOR_OUT_Y_RESIDUALS;
398 vectors += VECTOR_OUT_Y_PARAMETERS;
399 vectors += VECTOR_OUT_Y_COVARIANCE;
400 vectors += VECTOR_OUT_Y_PARAMETERS;
401 return vectors;
402 }
403
404
outputScalarList() const405 QStringList FitGaussianWeightedSource::outputScalarList() const {
406 return QStringList( SCALAR_OUT );
407 }
408
409
outputStringList() const410 QStringList FitGaussianWeightedSource::outputStringList() const {
411 return QStringList( /*STRING_OUT*/ );
412 }
413
414
saveProperties(QXmlStreamWriter & s)415 void FitGaussianWeightedSource::saveProperties(QXmlStreamWriter &s) {
416 QString force_offset;
417 force_offset.setNum(_forceOffset);
418 s.writeAttribute("ForceOffset", force_offset);
419 }
420
421
parameterName(int index) const422 QString FitGaussianWeightedSource::parameterName(int index) const {
423 QString parameter;
424 switch (index) {
425 case 0:
426 parameter = "A";
427 break;
428 case 1:
429 parameter = "\\sigma";
430 break;
431 case 2:
432 parameter = "x_o";
433 break;
434 case 3:
435 parameter = "C";
436 break;
437 default:
438 parameter = "";
439 break;
440 }
441
442 return parameter;
443 }
444
445
446 // Name used to identify the plugin. Used when loading the plugin.
pluginName() const447 QString FitGaussianWeightedPlugin::pluginName() const { return tr("Gaussian Weighted Fit"); }
pluginDescription() const448 QString FitGaussianWeightedPlugin::pluginDescription() const { return tr("Generates a gaussian weighted fit for a set of data."); }
449
450
create(Kst::ObjectStore * store,Kst::DataObjectConfigWidget * configWidget,bool setupInputsOutputs) const451 Kst::DataObject *FitGaussianWeightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
452
453 if (ConfigWidgetFitGaussianWeightedPlugin* config = static_cast<ConfigWidgetFitGaussianWeightedPlugin*>(configWidget)) {
454
455 FitGaussianWeightedSource* object = store->createObject<FitGaussianWeightedSource>();
456 object->_forceOffset = config->_forceOffset->isChecked();
457
458 if (setupInputsOutputs) {
459 object->setupOutputs();
460 object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
461 object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
462 object->setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
463 object->setInputScalar(SCALAR_IN_OFFSET, config->scalarOffset());
464 }
465
466 object->setPluginName(pluginName());
467
468 object->writeLock();
469 object->registerChange();
470 object->unlock();
471
472 return object;
473 }
474 return 0;
475 }
476
477
configWidget(QSettings * settingsObject) const478 Kst::DataObjectConfigWidget *FitGaussianWeightedPlugin::configWidget(QSettings *settingsObject) const {
479 ConfigWidgetFitGaussianWeightedPlugin *widget = new ConfigWidgetFitGaussianWeightedPlugin(settingsObject);
480 return widget;
481 }
482
483 #ifndef QT5
484 Q_EXPORT_PLUGIN2(kstplugin_FitGaussianWeightedPlugin, FitGaussianWeightedPlugin)
485 #endif
486
487 // vim: ts=2 sw=2 et
488