File indexing completed on 2024-12-22 04:18:11

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 "fitexponential_weighted.h"
0017 #include "objectstore.h"
0018 #include "ui_fitexponential_weightedconfig.h"
0019 
0020 #define NUM_PARAMS 4
0021 #define MAX_NUM_ITERATIONS 500
0022 
0023 #include <gsl/gsl_fit.h>
0024 #include "../non_linear_weighted.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_IN_WEIGHTS = "Weights Vector";
0029 static const QString& VECTOR_OUT_Y_FITTED = "Fit";
0030 static const QString& VECTOR_OUT_Y_RESIDUALS = "Residuals";
0031 static const QString& VECTOR_OUT_Y_PARAMETERS = "Parameters Vector";
0032 static const QString& VECTOR_OUT_Y_COVARIANCE = "Covariance";
0033 static const QString& SCALAR_OUT = "chi^2/nu";
0034 
0035 class ConfigWidgetFitExponentialWeightedPlugin : public Kst::DataObjectConfigWidget, public Ui_FitExponential_WeightedConfig {
0036   public:
0037     ConfigWidgetFitExponentialWeightedPlugin(QSettings* cfg) : DataObjectConfigWidget(cfg), Ui_FitExponential_WeightedConfig() {
0038       _store = 0;
0039       setupUi(this);
0040     }
0041 
0042     ~ConfigWidgetFitExponentialWeightedPlugin() {}
0043 
0044     void setObjectStore(Kst::ObjectStore* store) { 
0045       _store = store; 
0046       _vectorX->setObjectStore(store);
0047       _vectorY->setObjectStore(store);
0048       _vectorWeights->setObjectStore(store);
0049     }
0050 
0051     void setupSlots(QWidget* dialog) {
0052       if (dialog) {
0053         connect(_vectorX, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0054         connect(_vectorY, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0055         connect(_vectorWeights, SIGNAL(selectionChanged(QString)), dialog, SIGNAL(modified()));
0056       }
0057     }
0058 
0059     void setVectorX(Kst::VectorPtr vector) {
0060       setSelectedVectorX(vector);
0061     }
0062 
0063     void setVectorY(Kst::VectorPtr vector) {
0064       setSelectedVectorY(vector);
0065     }
0066 
0067     void setVectorsLocked(bool locked = true) {
0068       _vectorX->setEnabled(!locked);
0069       _vectorY->setEnabled(!locked);
0070     }
0071 
0072     Kst::VectorPtr selectedVectorX() { return _vectorX->selectedVector(); };
0073     void setSelectedVectorX(Kst::VectorPtr vector) { return _vectorX->setSelectedVector(vector); };
0074 
0075     Kst::VectorPtr selectedVectorY() { return _vectorY->selectedVector(); };
0076     void setSelectedVectorY(Kst::VectorPtr vector) { return _vectorY->setSelectedVector(vector); };
0077 
0078     Kst::VectorPtr selectedVectorWeights() { return _vectorWeights->selectedVector(); };
0079     void setSelectedVectorWeights(Kst::VectorPtr vector) { return _vectorWeights->setSelectedVector(vector); };
0080 
0081     virtual void setupFromObject(Kst::Object* dataObject) {
0082       if (FitExponentialWeightedSource* source = static_cast<FitExponentialWeightedSource*>(dataObject)) {
0083         setSelectedVectorX(source->vectorX());
0084         setSelectedVectorY(source->vectorY());
0085         setSelectedVectorWeights(source->vectorWeights());
0086       }
0087     }
0088 
0089     virtual bool configurePropertiesFromXml(Kst::ObjectStore *store, QXmlStreamAttributes& attrs) {
0090       Q_UNUSED(store);
0091       Q_UNUSED(attrs);
0092 
0093       bool validTag = true;
0094 
0095 //       QStringRef av;
0096 //       av = attrs.value("value");
0097 //       if (!av.isNull()) {
0098 //         _configValue = QVariant(av.toString()).toBool();
0099 //       }
0100 
0101       return validTag;
0102     }
0103 
0104   public slots:
0105     virtual void save() {
0106       if (_cfg) {
0107         _cfg->beginGroup("Fit Exponential Weighted Plugin");
0108         _cfg->setValue("Input Vector X", _vectorX->selectedVector()->Name());
0109         _cfg->setValue("Input Vector Y", _vectorY->selectedVector()->Name());
0110         _cfg->setValue("Input Vector Weights", _vectorWeights->selectedVector()->Name());
0111         _cfg->endGroup();
0112       }
0113     }
0114 
0115     virtual void load() {
0116       if (_cfg && _store) {
0117         _cfg->beginGroup("Fit Exponential Weighted Plugin");
0118         QString vectorName = _cfg->value("Input Vector X").toString();
0119         Kst::Object* object = _store->retrieveObject(vectorName);
0120         Kst::Vector* vectorx = static_cast<Kst::Vector*>(object);
0121         if (vectorx) {
0122           setSelectedVectorX(vectorx);
0123         }
0124         vectorName = _cfg->value("Input Vector Y").toString();
0125         object = _store->retrieveObject(vectorName);
0126         Kst::Vector* vectory = static_cast<Kst::Vector*>(object);
0127         if (vectory) {
0128           setSelectedVectorX(vectory);
0129         }
0130         vectorName = _cfg->value("Input Vector Weights").toString();
0131         object = _store->retrieveObject(vectorName);
0132         Kst::Vector* vectorweights = static_cast<Kst::Vector*>(object);
0133         if (vectorweights) {
0134           setSelectedVectorX(vectorweights);
0135         }
0136         _cfg->endGroup();
0137       }
0138     }
0139 
0140   private:
0141     Kst::ObjectStore *_store;
0142 
0143 };
0144 
0145 
0146 FitExponentialWeightedSource::FitExponentialWeightedSource(Kst::ObjectStore *store)
0147 : Kst::BasicPlugin(store) {
0148 }
0149 
0150 
0151 FitExponentialWeightedSource::~FitExponentialWeightedSource() {
0152 }
0153 
0154 
0155 QString FitExponentialWeightedSource::_automaticDescriptiveName() const {
0156   return tr("%1 Weighted Exponential").arg(vectorY()->descriptiveName());
0157 }
0158 
0159 
0160 void FitExponentialWeightedSource::change(Kst::DataObjectConfigWidget *configWidget) {
0161   if (ConfigWidgetFitExponentialWeightedPlugin* config = static_cast<ConfigWidgetFitExponentialWeightedPlugin*>(configWidget)) {
0162     setInputVector(VECTOR_IN_X, config->selectedVectorX());
0163     setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0164     setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
0165   }
0166 }
0167 
0168 
0169 void FitExponentialWeightedSource::setupOutputs() {
0170   setOutputVector(VECTOR_OUT_Y_FITTED, "");
0171   setOutputVector(VECTOR_OUT_Y_RESIDUALS, "");
0172   setOutputVector(VECTOR_OUT_Y_PARAMETERS, "");
0173   setOutputVector(VECTOR_OUT_Y_COVARIANCE, "");
0174   setOutputScalar(SCALAR_OUT, "");
0175 }
0176 
0177 double _X0 = 0; // this use of a global is a hack to inject the first sample into the function.
0178 
0179 void swapDouble(double *A, double *B) {
0180   double C;
0181 
0182   C = *A;
0183   *A = *B;
0184   *B = C;
0185 }
0186 
0187 
0188 void function_initial_estimate( const double x[], const double y[], int npts, double P0[] ) {
0189   Q_UNUSED( x )
0190   Q_UNUSED( y )
0191   Q_UNUSED( npts )
0192 
0193   _X0 = x[0];
0194 
0195   // determine the signs of the terms.
0196   // get the average of the first 5%, last 5% and middle 5% of points
0197   int n_ave = npts/20;
0198   if (n_ave < 1) n_ave = 1;
0199   if (n_ave > 100) n_ave = 100;
0200 
0201   double y0 = 0;
0202   double x0 = 0;
0203   double x1 = 0;
0204   double y1 = 0;
0205   double x2 = 0;
0206   double y2 = 0;
0207   int d1 = npts/2 - n_ave/2;
0208   int d2 = npts-n_ave;
0209 
0210   if ((d1 + n_ave > npts) || (d2 + n_ave > npts)) { // bail if not enough points.
0211     P0[0] =  1.0;
0212     P0[1] =  0.0;
0213     P0[2] =  0.0;
0214     return;
0215   }
0216 
0217   for (int i=0; i<n_ave; i++) {
0218     x0+=x[i];
0219     y0+=y[i];
0220     x1+=x[i+d1];
0221     y1+=y[i+d1];
0222     x2+=x[i+d2];
0223     y2+=y[i+d2];
0224   }
0225   x0 /= (double)n_ave;
0226   y0 /= (double)n_ave;
0227   x1 /= (double)n_ave;
0228   y1 /= (double)n_ave;
0229   x2 /= (double)n_ave;
0230   y2 /= (double)n_ave;
0231 
0232   // Make sure x0, x1, x2 are monotonic.
0233   if (x2 > x0) {
0234     if (x1 > x2) {
0235       swapDouble(&x1, &x2);
0236       swapDouble(&y1, &y2);
0237     }
0238     if (x1 < x0) {
0239       swapDouble(&x1, &x0);
0240       swapDouble(&y1, &y0);
0241     }
0242   } else {
0243     if (x1 < x2) {
0244       swapDouble(&x1, &x2);
0245       swapDouble(&y1, &y2);
0246     }
0247     if (x1 > x0) {
0248       swapDouble(&x1, &x0);
0249       swapDouble(&y1, &y0);
0250     }
0251   }
0252 
0253   if ((x1 == x0) || (x2 == x0) || (x1 == x2)) { // bail if no x range
0254     P0[0] =  1.0;
0255     P0[1] =  0.0;
0256     P0[2] =  0.0;
0257     return;
0258   }
0259 
0260   P0[0] =  fabs(y2 - y0)/M_E;
0261   P0[1] =  M_E/fabs(x2 - x0);
0262   P0[2] =  y0;
0263 
0264   double m = (y2 - y0)/(x2 - x0);
0265   if (m > 0) { // rising
0266     if ((x1-x0)*m + y0 > y1) { // neg curvature +A, +B
0267       //P0[0] *= -1;
0268       //P0[1] *= -1.0;
0269     } else {   // -A, -B
0270       P0[0] *= -1;
0271       P0[1] *= -1.0;
0272     }
0273   } else { // falling
0274     if ((x1-x0)*m + y0 > y1) { // Curving down +A, -B
0275       //P0[0] *= -1;
0276       P0[1] *= -1.0;
0277     } else { // -A, +B
0278       P0[0] *= -1;
0279       //P0[1] *= -1.0;
0280     }
0281   }
0282 
0283   fflush(stdout);
0284 }
0285 
0286 // It would be nice to fit exp(B*(x-_X0)) + C, and not have A,
0287 // but A is required to at least hold the sign.  We we just fix
0288 // _X0 to the first X value in the data set, and fit for A.
0289 double function_calculate( double x, double* P ) {
0290   double A = P[0];
0291   double B = P[1];
0292   double C = P[2];
0293   double y;
0294 
0295   y  = ( A*exp( B*(x - _X0) ) ) + C;
0296 
0297   return y;
0298 }
0299 
0300 
0301 void function_derivative( double x, double* P, double* dFdP ) {
0302   double A  = P[0];
0303   double B = P[1];
0304   double exp_BxXo;
0305 
0306   // dFdA = exp(b*(x-_X0)
0307   // dFdB = A*(x-_X0)*exp(b*(x-_X0))
0308 
0309   exp_BxXo = exp( B * (x-_X0) );
0310 
0311   dFdP[0]   = exp_BxXo;
0312   dFdP[1]   = (x-_X0) * A * exp_BxXo;
0313   dFdP[2]   = 1.0;
0314 
0315 }
0316 
0317 
0318 
0319 
0320 bool FitExponentialWeightedSource::algorithm() {
0321   Kst::VectorPtr inputVectorX = _inputVectors[VECTOR_IN_X];
0322   Kst::VectorPtr inputVectorY = _inputVectors[VECTOR_IN_Y];
0323   Kst::VectorPtr inputVectorWeights = _inputVectors[VECTOR_IN_WEIGHTS];
0324 
0325   Kst::VectorPtr outputVectorYFitted = _outputVectors[VECTOR_OUT_Y_FITTED];
0326   Kst::VectorPtr outputVectorYResiduals = _outputVectors[VECTOR_OUT_Y_RESIDUALS];
0327   Kst::VectorPtr outputVectorYParameters = _outputVectors[VECTOR_OUT_Y_PARAMETERS];
0328   Kst::VectorPtr outputVectorYCovariance = _outputVectors[VECTOR_OUT_Y_COVARIANCE];
0329   Kst::ScalarPtr outputScalar = _outputScalars[SCALAR_OUT];
0330 
0331   Kst::LabelInfo label_info = inputVectorY->labelInfo();
0332 
0333   _X0 = inputVectorX->noNanValue()[0];
0334   n_params = 3;
0335 
0336   label_info.name = tr("A\\Cdotexp((x-x_o)/\\tau) + C fit to %1").arg(label_info.name);
0337   outputVectorYFitted->setLabelInfo(label_info);
0338 
0339   label_info.name = tr("Exponential Fit Residuals");
0340   outputVectorYResiduals->setLabelInfo(label_info);
0341 
0342 
0343   bool bReturn = false;
0344 
0345   bReturn = kstfit_nonlinear_weighted( inputVectorX, inputVectorY, inputVectorWeights,
0346                               outputVectorYFitted, outputVectorYResiduals, outputVectorYParameters,
0347                               outputVectorYCovariance, outputScalar );
0348 
0349   outputVectorYParameters->raw_V_ptr()[1] = 1.0/outputVectorYParameters->raw_V_ptr()[1];
0350   outputVectorYParameters->raw_V_ptr()[3] = _X0;
0351 
0352   return bReturn;
0353 }
0354 
0355 
0356 Kst::VectorPtr FitExponentialWeightedSource::vectorX() const {
0357   return _inputVectors[VECTOR_IN_X];
0358 }
0359 
0360 
0361 Kst::VectorPtr FitExponentialWeightedSource::vectorY() const {
0362   return _inputVectors[VECTOR_IN_Y];
0363 }
0364 
0365 
0366 Kst::VectorPtr FitExponentialWeightedSource::vectorWeights() const {
0367   return _inputVectors[VECTOR_IN_WEIGHTS];
0368 }
0369 
0370 
0371 QStringList FitExponentialWeightedSource::inputVectorList() const {
0372   QStringList vectors(VECTOR_IN_X);
0373   vectors += VECTOR_IN_Y;
0374   vectors += VECTOR_IN_WEIGHTS;
0375   return vectors;
0376 }
0377 
0378 
0379 QStringList FitExponentialWeightedSource::inputScalarList() const {
0380   return QStringList();
0381 }
0382 
0383 
0384 QStringList FitExponentialWeightedSource::inputStringList() const {
0385   return QStringList( /*STRING_IN*/ );
0386 }
0387 
0388 
0389 QStringList FitExponentialWeightedSource::outputVectorList() const {
0390   QStringList vectors(VECTOR_OUT_Y_FITTED);
0391   vectors += VECTOR_OUT_Y_RESIDUALS;
0392   vectors += VECTOR_OUT_Y_PARAMETERS;
0393   vectors += VECTOR_OUT_Y_COVARIANCE;
0394   vectors += VECTOR_OUT_Y_PARAMETERS;
0395   return vectors;
0396 }
0397 
0398 
0399 QStringList FitExponentialWeightedSource::outputScalarList() const {
0400   return QStringList( SCALAR_OUT );
0401 }
0402 
0403 
0404 QStringList FitExponentialWeightedSource::outputStringList() const {
0405   return QStringList( /*STRING_OUT*/ );
0406 }
0407 
0408 
0409 void FitExponentialWeightedSource::saveProperties(QXmlStreamWriter &s) {
0410   Q_UNUSED(s);
0411 //   s.writeAttribute("value", _configValue);
0412 }
0413 
0414 
0415 QString FitExponentialWeightedSource::parameterName(int index) const {
0416   QString parameter;
0417   switch (index) {
0418     case 0:
0419       parameter = "A";
0420       break;
0421     case 1:
0422       parameter = "\\tau";
0423       break;
0424     case 2:
0425       parameter = "C";
0426       break;
0427     case 3:
0428       parameter = "X_o";
0429       break;
0430     default:
0431       parameter = "";
0432       break;
0433   }
0434 
0435   return parameter;
0436 }
0437 
0438 
0439 // Name used to identify the plugin.  Used when loading the plugin.
0440 QString FitExponentialWeightedPlugin::pluginName() const { return tr("Exponential Weighted Fit"); }
0441 QString FitExponentialWeightedPlugin::pluginDescription() const { return tr("Generates an exponential weighted fit for a set of data."); }
0442 
0443 
0444 Kst::DataObject *FitExponentialWeightedPlugin::create(Kst::ObjectStore *store, Kst::DataObjectConfigWidget *configWidget, bool setupInputsOutputs) const {
0445 
0446   if (ConfigWidgetFitExponentialWeightedPlugin* config = static_cast<ConfigWidgetFitExponentialWeightedPlugin*>(configWidget)) {
0447 
0448     FitExponentialWeightedSource* object = store->createObject<FitExponentialWeightedSource>();
0449 
0450     if (setupInputsOutputs) {
0451       object->setupOutputs();
0452       object->setInputVector(VECTOR_IN_X, config->selectedVectorX());
0453       object->setInputVector(VECTOR_IN_Y, config->selectedVectorY());
0454       object->setInputVector(VECTOR_IN_WEIGHTS, config->selectedVectorWeights());
0455     }
0456 
0457     object->setPluginName(pluginName());
0458 
0459     object->writeLock();
0460     object->registerChange();
0461     object->unlock();
0462 
0463     return object;
0464   }
0465   return 0;
0466 }
0467 
0468 
0469 Kst::DataObjectConfigWidget *FitExponentialWeightedPlugin::configWidget(QSettings *settingsObject) const {
0470   ConfigWidgetFitExponentialWeightedPlugin *widget = new ConfigWidgetFitExponentialWeightedPlugin(settingsObject);
0471   return widget;
0472 }
0473 
0474 #ifndef QT5
0475 Q_EXPORT_PLUGIN2(kstplugin_FitExponentialWeightedPlugin, FitExponentialWeightedPlugin)
0476 #endif
0477 
0478 // vim: ts=2 sw=2 et