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