fits/kneefrequency/fitkneefrequency.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 "fitkneefrequency.h"
0017 #include "objectstore.h"
0018 #include "ui_fitkneefrequencyconfig.h"
0019
0020 #include <gsl/gsl_fit.h>
0021 #include "../common.h"
0022
0023 #define KNEEFREQ_NUMPARAMETERS 5
0024
0025 static const QString& VECTOR_IN_X = "X Vector";
0026 static const QString& VECTOR_IN_Y = "Y Vector";
0027 static const QString& SCALAR_IN_MAX = "Max 1/f^a Freq Scalar";
0028 static const QString& SCALAR_IN_MIN = "Min. White Noise Freq Scalar";
0029 static const QString& SCALAR_IN_WHITENOISE = "White Noise C Scalar";
0030 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
0031 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
0032 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
0033
0034 class ConfigWidgetFitKneeFrequencyPlugin : public Kst::DataObjectConfigWidget, public Ui_FitKneeFrequencyConfig {
0035   public:
0036     ConfigWidgetFitKneeFrequencyPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitKneeFrequencyConfig() {
0037       _store = 0;
0038       setupUi(this);
0039     }
0040
0041     ~ConfigWidgetFitKneeFrequencyPlugin() {}
0042
0043     void setObjectStore(Kst::ObjectStore* store) {
0044       _store = store;
0045       _vectorX->setObjectStore(store);
0046       _vectorY->setObjectStore(store);
0047       _scalarMax->setObjectStore(store);
0048       _scalarMin->setObjectStore(store);
0049       _scalarWhiteNoise->setObjectStore(store);
0050     }
0051
0052     void setupSlots(QWidget* dialog) {
0053       if (dialog) {
0054         connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0055         connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0056         connect(_scalarMax, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0057         connect(_scalarMin, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0058         connect(_scalarWhiteNoise, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0059       }
0060     }
0061
0062
0063     void setVectorX(Kst::VectorPtr vector) {
0064       setSelectedVectorX(vector);
0065     }
0066
0067     void setVectorY(Kst::VectorPtr vector) {
0068       setSelectedVectorY(vector);
0069     }
0070
0071     void setVectorsLocked(bool locked = true) {
0072       _vectorX->setEnabled(!locked);
0073       _vectorY->setEnabled(!locked);
0074     }
0075
0076     Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
0077     void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
0078
0079     Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
0080     void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
0081
0082     Kst::ScalarPtr selectedScalarMax() { return _scalarMax->selectedScalar(); };
0083     void setSelectedScalarMax(Kst::ScalarPtr scalar) { return _scalarMax->setSelectedScalar(scalar); };
0084
0085     Kst::ScalarPtr selectedScalarMin() { return _scalarMin->selectedScalar(); };
0086     void setSelectedScalarMin(Kst::ScalarPtr scalar) { return _scalarMin->setSelectedScalar(scalar); };
0087
0088     Kst::ScalarPtr selectedScalarWhiteNoise() { return _scalarWhiteNoise->selectedScalar(); };
0089     void setSelectedScalarWhiteNoise(Kst::ScalarPtr scalar) { return _scalarWhiteNoise->setSelectedScalar(scalar); };
0090
0091     virtual void setupFromObject(Kst::Object* dataObject) {
0092       if (FitKneeFrequencySource* source = static_cast<FitKneeFrequencySource*>(dataObject)) {
0093         setSelectedVectorX(source->vectorX());
0094         setSelectedVectorY(source->vectorY());
0095         setSelectedScalarMax(source->scalarMax());
0096         setSelectedScalarMin(source->scalarMin());
0097         setSelectedScalarWhiteNoise(source->scalarWhiteNoise());
0098       }
0099     }
0100
0101     virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
0102       Q_UNUSED(store);
0103       Q_UNUSED(attrs);
0104
0105       bool validTag = true;
0106
0107 //       QStringRef av;
0108 //       av = attrs.value("value");
0109 //       if (!av.isNull()) {
0110 //         _configValue = QVariant(av.toString()).toBool();
0111 //       }
0112
0113       return validTag;
0114     }
0115
0116   public slots:
0117     virtual void save() {
0118       if (_cfg) {
0119         _cfg->beginGroup("Fit Knee Frequency Plugin");
0120         _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
0121         _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
0122         _cfg->setValue("Input Scalar Max 1/f^a Freq", _scalarMax->selectedScalar()->Name());
0123         _cfg->setValue("Input Scalar Min. White Noise Freq", _scalarMin->selectedScalar()->Name());
0124         _cfg->setValue("Input Scalar White Noise C", _scalarWhiteNoise->selectedScalar()->Name());
0125         _cfg->endGroup();
0126       }
0127     }
0128
0129     virtual void load() {
0130       if (_cfg && _store) {
0131         _cfg->beginGroup("Fit Knee Frequency Plugin");
0132         QString vectorName = _cfg->value("Input Vector X").toString();
0133         Kst::Object* object = _store->retrieveObject(vectorName);
0134         Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
0135         if (vectorx) {
0136           setSelectedVectorX(vectorx);
0137         }
0138         vectorName = _cfg->value("Input Vector Y").toString();
0139         object = _store->retrieveObject(vectorName);
0140         Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
0141         if (vectory) {
0142           setSelectedVectorX(vectory);
0143         }
0144         QString scalarName = _cfg->value("Input Scalar Max 1/f^a Freq").toString();
0145         object = _store->retrieveObject(scalarName);
0146         Kst::Scalar* maxScalar = static_cast<Kst::Scalar*>(object);
0147         if (maxScalar) {
0148           setSelectedScalarMax(maxScalar);
0149         }
0150         scalarName = _cfg->value("Input Scalar Min. White Noise Freq").toString();
0151         object = _store->retrieveObject(scalarName);
0152         Kst::Scalar* minScalar = static_cast<Kst::Scalar*>(object);
0153         if (minScalar) {
0154           setSelectedScalarMin(minScalar);
0155         }
0156         scalarName = _cfg->value("Input Scalar White Noise C").toString();
0157         object = _store->retrieveObject(scalarName);
0158         Kst::Scalar* whiteNoiseScalar = static_cast<Kst::Scalar*>(object);
0159         if (whiteNoiseScalar) {
0160           setSelectedScalarWhiteNoise(whiteNoiseScalar);
0161         }
0162         _cfg->endGroup();
0163       }
0164     }
0165
0166   private:
0167     Kst::ObjectStore *_store;
0168
0169 };
0170
0171
0172 FitKneeFrequencySource::FitKneeFrequencySource(Kst::ObjectStore *store)
0173 : Kst::BasicPlugin(store) {
0174 }
0175
0176
0177 FitKneeFrequencySource::~FitKneeFrequencySource() {
0178 }
0179
0180
0181 QString FitKneeFrequencySource::_automaticDescriptiveName() const {
0182   return tr("Fit Knee Frequency Plugin");
0183 }
0184
0185
0186 void FitKneeFrequencySource::change(Kst::DataObjectConfigWidget *configWidget) {
0187   if (ConfigWidgetFitKneeFrequencyPlugin* config = static_cast<ConfigWidgetFitKneeFrequencyPlugin*>(configWidget)) {
0188     setInputVector(VECTOR_IN_X, config->selectedVectorX());
0189     setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0190     setInputScalar(SCALAR_IN_MAX, config->selectedScalarMax());
0191     setInputScalar(SCALAR_IN_MIN, config->selectedScalarMin());
0192     setInputScalar(SCALAR_IN_WHITENOISE, config->selectedScalarWhiteNoise());
0193   }
0194 }
0195
0196
0197 void FitKneeFrequencySource::setupOutputs() {
0198   setOutputVector(VECTOR_OUT_Y_FITTED, "");
0199   setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
0200   setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
0201 }
0202
0203
0204 bool FitKneeFrequencySource::algorithm() {
0205   Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
0206   Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
0207   Kst::ScalarPtr inputScalarMax = _inputScalars[SCALAR_IN_MAX];
0208   Kst::ScalarPtr inputScalarMin = _inputScalars[SCALAR_IN_MIN];
0209   Kst::ScalarPtr inputScalarWhiteNoise = _inputScalars[SCALAR_IN_WHITENOISE];
0210
0211   Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
0212   Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
0213   Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
0214
0215   if (inputVectorX->length() != inputVectorY->length())  {
0216     _errorString = tr("Error:  Input Vector Sizes do not match");
0217     return false;
0218   }
0219   if (inputVectorX->length() < 1)  {
0220     _errorString = tr("Error:  Input Vector X invalid");
0221     return false;
0222   }
0223
0224   int inArraysLength = inputVectorX->length();
0225
0226   outputVectorYFitted->resize(inArraysLength);
0227   outputVectorYResiduals->resize(inArraysLength);
0228   outputVectorYParameters->resize(KNEEFREQ_NUMPARAMETERS);
0229
0230   double xi, yi;
0231   int i;
0232   double maxOneOverFFreq, minWhiteNoiseFreq, whiteNoiseC;
0233
0234   maxOneOverFFreq = inputScalarMax->value();
0235   minWhiteNoiseFreq = inputScalarMin->value();
0236   whiteNoiseC = inputScalarWhiteNoise->value();
0237
0238   int maxOneOverFIndex, minWhiteNoiseIndex;
0239
0240   //fast calculation of index for maxOneOverFFreq
0241   int i_bot = 0;
0242   int i_top = inArraysLength - 1;
0243
0244   while (i_bot + 1 < i_top) {
0245     int i0 = (i_top + i_bot)/2;
0246     if (maxOneOverFFreq < inputVectorX->value()[i0]) {
0247       i_top = i0;
0248     } else {
0249       i_bot = i0;
0250     }
0251   }
0252   maxOneOverFIndex = i_top; //top because we use i_bot+1.
0253
0254   //fast calculation of index for minWhiteNoiseFreq
0255   i_bot = 0;
0256   i_top = inArraysLength - 1;
0257
0258   while (i_bot + 1 < i_top) {
0259     int i0 = (i_top + i_bot)/2;
0260     if (minWhiteNoiseFreq < inputVectorX->value()[i0]) {
0261       i_top = i0;
0262     } else {
0263       i_bot = i0;
0264     }
0265   }
0266   minWhiteNoiseIndex = i_top;
0267
0268   //verify calculated indices.
0269   if ( !(maxOneOverFIndex>0) || !(minWhiteNoiseIndex>=maxOneOverFIndex) || !(minWhiteNoiseIndex<(inArraysLength-1)) ) {
0270     _errorString = tr("Error:  Input Frequencies are Invalid\n");
0271     return false;
0272   }
0273
0274   // calculate white noise limit
0275   double sumY, sumY2;
0276   sumY = sumY2 = 0;
0277
0278   for (i = minWhiteNoiseIndex; i < inArraysLength; ++i) {
0279     yi = inputVectorY->value()[i];
0280     sumY    +=  yi;
0281     sumY2   +=  pow(yi,2);
0282   }
0283
0284   double ybar, ysigma;
0285   ybar = sumY/(inArraysLength - minWhiteNoiseIndex);
0286   ysigma = sqrt((sumY2 - 2*ybar*sumY + pow(ybar,2)*(inArraysLength - minWhiteNoiseIndex))/(inArraysLength - minWhiteNoiseIndex));
0287   // end calculate white noise limit
0288
0289   // fit 1/f noise
0290   double sumLnXLnY, sumLnX, sumLnY, sumLnX2;
0291   sumLnXLnY = sumLnX = sumLnY = sumLnX2 = 0;
0292
0293   for (i = 0; i < maxOneOverFIndex; ++i) {
0294     xi = inputVectorX->value()[i];
0295     yi = inputVectorY->value()[i];
0296
0297     if (!(xi>0) || !((yi-ybar)>0)) {
0298       _errorString = tr("Error:  Input Data Invalid.\n");
0299       return false;
0300     }
0301
0302     sumLnXLnY += log(xi)*log(yi-ybar); //-ybar to isolate 1/f noise.
0303     sumLnX    += log(xi);
0304     sumLnY    += log(yi-ybar);
0305     sumLnX2   += pow(log(xi),2);
0306   }
0307
0308   double a, b;
0309   a = (maxOneOverFIndex*sumLnXLnY - sumLnX*sumLnY)/(maxOneOverFIndex*sumLnX2 - pow(sumLnX,2));
0310   b = exp((sumLnY - a*sumLnX)/maxOneOverFIndex);
0311   // end fit 1/f noise
0312
0313   double knee_freq = pow(ybar*whiteNoiseC/b,1.0/a); // calculate knee frequency.
0314
0315   // output fit data
0316   for (i = 0; i < maxOneOverFIndex; ++i) {
0317       outputVectorYFitted->raw_V_ptr()[i] = b * pow(inputVectorX->value()[i],a) + ybar;
0318       outputVectorYResiduals->raw_V_ptr()[i] = inputVectorY->value()[i] - outputVectorYFitted->raw_V_ptr()[i];
0319   }
0320
0321   for (i = maxOneOverFIndex; i < minWhiteNoiseIndex; ++i) { // zeros for unfitted region.
0322       outputVectorYFitted->raw_V_ptr()[i] = 0;
0323       outputVectorYResiduals->raw_V_ptr()[i] = 0;
0324   }
0325
0326   for (i = minWhiteNoiseIndex; i < inArraysLength; ++i) {
0327       outputVectorYFitted->raw_V_ptr()[i] = ybar;
0328       outputVectorYResiduals->raw_V_ptr()[i] = outputVectorYFitted->raw_V_ptr()[i] - ybar;
0329   }
0330
0331   outputVectorYParameters->raw_V_ptr()[0] = ybar;
0332   outputVectorYParameters->raw_V_ptr()[1] = ysigma;
0333   outputVectorYParameters->raw_V_ptr()[2] = b;
0334   outputVectorYParameters->raw_V_ptr()[3] = -a;
0335   outputVectorYParameters->raw_V_ptr()[4] = knee_freq;
0336
0337   return true;
0338 }
0339
0340
0341 Kst::VectorPtr FitKneeFrequencySource::vectorX() const {
0342   return _inputVectors[VECTOR_IN_X];
0343 }
0344
0345
0346 Kst::VectorPtr FitKneeFrequencySource::vectorY() const {
0347   return _inputVectors[VECTOR_IN_Y];
0348 }
0349
0350
0351 Kst::ScalarPtr FitKneeFrequencySource::scalarMax() const {
0352   return _inputScalars[SCALAR_IN_MAX];
0353 }
0354
0355
0356 Kst::ScalarPtr FitKneeFrequencySource::scalarMin() const {
0357   return _inputScalars[SCALAR_IN_MIN];
0358 }
0359
0360
0361 Kst::ScalarPtr FitKneeFrequencySource::scalarWhiteNoise() const {
0362   return _inputScalars[SCALAR_IN_WHITENOISE];
0363 }
0364
0365
0366 QStringList FitKneeFrequencySource::inputVectorList() const {
0367   QStringList vectors(VECTOR_IN_X);
0368   vectors += VECTOR_IN_Y;
0369   return vectors;
0370 }
0371
0372
0373 QStringList FitKneeFrequencySource::inputScalarList() const {
0374   QStringList scalars(SCALAR_IN_MAX);
0375   scalars += SCALAR_IN_MIN;
0376   scalars += SCALAR_IN_WHITENOISE;
0377   return scalars;
0378 }
0379
0380
0381 QStringList FitKneeFrequencySource::inputStringList() const {
0382   return QStringList( /*STRING_IN*/ );
0383 }
0384
0385
0386 QStringList FitKneeFrequencySource::outputVectorList() const {
0387   QStringList vectors(VECTOR_OUT_Y_FITTED);
0388   vectors += VECTOR_OUT_Y_RESIDUALS;
0389   vectors += VECTOR_OUT_Y_PARAMETERS;
0390   vectors += VECTOR_OUT_Y_PARAMETERS;
0391   return vectors;
0392 }
0393
0394
0395 QStringList FitKneeFrequencySource::outputScalarList() const {
0396   return QStringList( /*SCALAR_OUT */);
0397 }
0398
0399
0400 QStringList FitKneeFrequencySource::outputStringList() const {
0401   return QStringList( /*STRING_OUT*/ );
0402 }
0403
0404
0405 void FitKneeFrequencySource::saveProperties(QXmlStreamWriter &s) {
0406   Q_UNUSED(s);
0407 //   s.writeAttribute("value", _configValue);
0408 }
0409
0410
0411 QString FitKneeFrequencySource::parameterName(int index) const {
0412   QString parameter;
0413   switch (index) {
0414     case 0:
0415       parameter = "White Noise Limit";
0416       break;
0417     case 1:
0418       parameter = "White Noise Sigma";
0419       break;
0420     case 2:
0421       parameter = "1/f^a Amplitude";
0422       break;
0423     case 3:
0424       parameter = "1/f^a Power Law a";
0425       break;
0426     case 4:
0427       parameter = "Knee Frequency";
0428       break;
0429   }
0430
0431   return parameter;
0432 }
0433
0434
0435 // Name used to identify the plugin.  Used when loading the plugin.
0436 QString FitKneeFrequencyPlugin::pluginName() const { return tr("Knee Frequency Fit"); }
0437 QString FitKneeFrequencyPlugin::pluginDescription() const { return tr("Generates a knee frequency fit for a set of data."); }
0438
0439
0440 Kst::DataObject *FitKneeFrequencyPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
0441
0442   if (ConfigWidgetFitKneeFrequencyPlugin* config = static_cast<ConfigWidgetFitKneeFrequencyPlugin*>(configWidget)) {
0443
0444     Kst::ScalarPtr max;
0445     Kst::ScalarPtr min;
0446     Kst::ScalarPtr noise;
0447
0448     // access/create scalars before creating plugin
0449     // in order to preserve continuous scalar shortnames
0450     if (setupInputsOutputs) {
0451       max = config->selectedScalarMax();
0452       min = config->selectedScalarMin();
0453       noise = config->selectedScalarWhiteNoise();
0454     }
0455
0456     FitKneeFrequencySource* object = store->createObject<FitKneeFrequencySource>();
0457
0458     if (setupInputsOutputs) {
0459       object->setInputScalar(SCALAR_IN_MAX, max);
0460       object->setInputScalar(SCALAR_IN_MIN, min);
0461       object->setInputScalar(SCALAR_IN_WHITENOISE, noise);
0462       object->setupOutputs();
0463       object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
0464       object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0465     }
0466
0467     object->setPluginName(pluginName());
0468
0469     object->writeLock();
0470     object->registerChange();
0471     object->unlock();
0472
0473     return object;
0474   }
0475   return 0;
0476 }
0477
0478
0479 Kst::DataObjectConfigWidget *FitKneeFrequencyPlugin::configWidget(QSettings *settingsObject) const {
0480   ConfigWidgetFitKneeFrequencyPlugin *widget = new ConfigWidgetFitKneeFrequencyPlugin(settingsObject);
0481   return widget;
0482 }
0483
0484 #ifndef QT5
0485 Q_EXPORT_PLUGIN2(kstplugin_FitKneeFrequencyPlugin, FitKneeFrequencyPlugin)
0486 #endif
0487
0488 // vim: ts=2 sw=2 et