Changeset 517 for trunk/src


Ignore:
Timestamp:
02/28/05 19:10:07 (20 years ago)
Author:
mar637
Message:
  • updated to reflect Wim's changes to Functionals
Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDFitter.cc

    r483 r517  
    5353SDFitter::~SDFitter()
    5454{
    55     reset();
     55  reset();
    5656}
    5757
    5858void SDFitter::clear()
    5959{
    60     for (uInt i=0;i< funcs_.nelements();++i) {
    61         delete funcs_[i]; funcs_[i] = 0;
    62     }
    63     funcs_.resize(0, True);
    64     parameters_.resize();
    65     error_.resize();
    66     thefit_.resize();
    67     estimate_.resize();
    68     chisquared_ = 0.0;
    69 }
     60  for (uInt i=0;i< funcs_.nelements();++i) {
     61    delete funcs_[i]; funcs_[i] = 0;
     62  }
     63  parameters_.resize();
     64  error_.resize();
     65  thefit_.resize();
     66  estimate_.resize();
     67  chisquared_ = 0.0;
     68}
     69
    7070void SDFitter::reset()
    7171{
    72     clear();
    73     x_.resize();
    74     y_.resize();
    75     m_.resize();
     72  clear();
     73  x_.resize();
     74  y_.resize();
     75  m_.resize();
    7676}
    7777
    7878
    7979bool SDFitter::computeEstimate() {
    80     if (x_.nelements() == 0 || y_.nelements() == 0)
    81         throw (AipsError("No x/y data specified."));
    82 
    83     if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) == 0)
    84         return false;
    85     uInt n = funcs_.nelements();
    86     SpectralEstimate estimator(n);
    87     estimator.setQ(5);
    88     Int mn,mx;
    89     mn = 0;
    90     mx = m_.nelements()-1;
    91     for (uInt i=0; i<m_.nelements();++i) {
    92       if (m_[i]) {
    93         mn = i;
    94         break;
    95       }
    96     }
    97     for (uInt j=m_.nelements()-1; j>=0;--j) {
    98       if (m_[j]) {
    99         mx = j;
    100         break;
    101       }
    102     }
    103     //mn = 0+x_.nelements()/10;
    104     //mx = x_.nelements()-x_.nelements()/10;
    105     estimator.setRegion(mn,mx);
    106     //estimator.setWindowing(True);
    107     SpectralList listGauss = estimator.estimate(x_, y_);
    108     Gaussian1D<Float>* g;
    109     parameters_.resize(n*3);
    110     uInt count = 0;
    111     for (uInt i=0; i<n;i++) {
    112         g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
    113         if (g) {
    114             (*g)[0] = listGauss[i].getAmpl();
    115             (*g)[1] = listGauss[i].getCenter();
    116             (*g)[2] = listGauss[i].getFWHM();
    117             ++count;
    118         }
    119     }
    120     estimate_.resize();
    121     listGauss.evaluate(estimate_,x_);
    122     return true;
     80  if (x_.nelements() == 0 || y_.nelements() == 0)
     81    throw (AipsError("No x/y data specified."));
     82
     83  if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) == 0)
     84    return false;
     85  uInt n = funcs_.nelements();
     86  SpectralEstimate estimator(n);
     87  estimator.setQ(5);
     88  Int mn,mx;
     89  mn = 0;
     90  mx = m_.nelements()-1;
     91  for (uInt i=0; i<m_.nelements();++i) {
     92    if (m_[i]) {
     93      mn = i;
     94      break;
     95    }
     96  }
     97  for (uInt j=m_.nelements()-1; j>=0;--j) {
     98    if (m_[j]) {
     99      mx = j;
     100      break;
     101    }
     102  }
     103  mn = 0+x_.nelements()/10;
     104  mx = x_.nelements()-x_.nelements()/10;
     105  estimator.setRegion(mn,mx);
     106  //estimator.setWindowing(True);
     107  SpectralList listGauss = estimator.estimate(x_, y_);
     108  parameters_.resize(n*3);
     109  Gaussian1D<Float>* g = 0;
     110  for (uInt i=0; i<n;i++) {
     111    g = dynamic_cast<Gaussian1D<Float>* >(funcs_[i]);
     112    if (g) {
     113      (*g)[0] = listGauss[i].getAmpl();
     114      (*g)[1] = listGauss[i].getCenter();
     115      (*g)[2] = listGauss[i].getFWHM();
     116    }
     117  }
     118  estimate_.resize();
     119  listGauss.evaluate(estimate_,x_);
     120  return true;
    123121}
    124122
    125123std::vector<float> SDFitter::getEstimate() const
    126124{
    127     if (estimate_.nelements() == 0)
    128         throw (AipsError("No estimate set."));
    129     std::vector<float> stlout;
    130     estimate_.tovector(stlout);
    131     return stlout;
     125  if (estimate_.nelements() == 0)
     126    throw (AipsError("No estimate set."));
     127  std::vector<float> stlout;
     128  estimate_.tovector(stlout);
     129  return stlout;
    132130}
    133131
    … …  
    135133bool SDFitter::setExpression(const std::string& expr, int ncomp)
    136134{
    137     clear();
    138     if (expr == "gauss") {
    139         if (ncomp < 1) throw (AipsError("Need at least one gaussian to fit."));
    140         funcs_.resize(ncomp);
    141         for (Int k=0; k<ncomp; ++k) {
    142             funcs_[k] = new Gaussian1D<Float>();
    143         }
    144     } else if (expr == "poly") {
    145         funcs_.resize(1);
    146         funcs_[0] = new Polynomial<Float>(ncomp);
    147     } else {
    148         //cerr << " compiled functions not yet implemented" << endl;
    149         //funcs_.resize(1);
    150         //funcs_[0] = new CompiledFunction<Float>();
    151         //funcs_[0]->setFunction(String(expr));
    152         return false;
    153     };
    154     return true;
     135  clear();
     136  if (expr == "gauss") {
     137    if (ncomp < 1) throw (AipsError("Need at least one gaussian to fit."));
     138    funcs_.resize(ncomp);
     139    for (Int k=0; k<ncomp; ++k) {
     140      funcs_[k] = new Gaussian1D<Float>();
     141    }
     142  } else if (expr == "poly") {
     143    funcs_.resize(1);
     144    funcs_[0] = new Polynomial<Float>(ncomp);
     145  } else {
     146    cerr << " compiled functions not yet implemented" << endl;
     147    //funcs_.resize(1);
     148    //funcs_[0] = new CompiledFunction<Float>();
     149    //funcs_[0]->setFunction(String(expr));
     150    return false;
     151  }
     152  return true;
    155153}
    156154
    … …  
    278276
    279277bool SDFitter::fit() {
    280     NonLinearFitLM<Float> fitter;
    281     //CompiledFunction<AutoDiff<Float> > comp;
    282     //Polynomial<AutoDiff<Float> > poly;
    283     CompoundFunction<AutoDiff<Float> > func;
    284     if (dynamic_cast<Gaussian1D<Float>* >(funcs_[0]) != 0) {
    285         //computeEstimates();
    286         for (uInt i=0; i<funcs_.nelements(); i++) {
    287             Gaussian1D<AutoDiff<Float> > gauss;//(*funcs_[i]);
    288            
    289             for (uInt j=0; j<funcs_[i]->nparameters(); j++) {
    290                 gauss[j] = AutoDiff<Float>((*funcs_[i])[j],
    291                                            gauss.nparameters(), j);
    292                 gauss.mask(j) = funcs_[i]->mask(j);
    293             }
    294            
    295             func.addFunction(gauss);
    296         }
    297     } else if (dynamic_cast<Polynomial<Float>* >(funcs_[0]) != 0) {
    298         Polynomial<AutoDiff<Float> > poly(funcs_[0]->nparameters()-1);
    299         //Polynomial<AutoDiff<Float> > poly(*funcs_[0]);
    300         for (uInt j=0; j<funcs_[0]->nparameters(); j++) {
    301             poly[j] = AutoDiff<Float>(0, poly.nparameters(), j);
    302             poly.mask(j) = funcs_[0]->mask(j);
    303         }
    304         func.addFunction(poly);
    305     } else if (dynamic_cast<CompiledFunction<Float>* >(funcs_[0]) != 0) {
    306        
    307         //         CompiledFunction<AutoDiff<Float> > comp;
    308         //         for (uInt j=0; j<funcs_[0]->nparameters(); j++) {
    309 //             comp[j] = AutoDiff<Float>(0, comp.nparameters(), j);
    310 //             comp.mask(j) = funcs_[0]->mask(j);
    311 //         }
    312 //         func.addFunction(comp);
    313        
    314         cout << "NYI." << endl;
    315     } else {
    316         throw(AipsError("Fitter not set up correctly."));
    317     }
    318     fitter.setFunction(func);
    319     fitter.setMaxIter(50+funcs_.nelements()*10);
    320     // Convergence criterium
    321     fitter.setCriteria(0.001);
    322 
    323     // Fit
    324     Vector<Float> sigma(x_.nelements());
    325     sigma = 1.0;
    326 
    327     parameters_.resize();
    328     parameters_ = fitter.fit(x_, y_, sigma, &m_);
    329 
    330     error_.resize();
    331     error_ = fitter.errors();
    332 
    333     chisquared_ = fitter.getChi2();
    334 
    335     residual_.resize();
    336     residual_ =  y_;
    337     fitter.residual(residual_,x_);
    338 
    339     // use fitter.residual(model=True) to get the model
    340     thefit_.resize(x_.nelements());
    341     fitter.residual(thefit_,x_,True);
    342     return true;
    343 }
     278  NonLinearFitLM<Float> fitter;
     279  CompoundFunction<Float> func;
     280  const uInt n = funcs_.nelements();
     281  for (uInt i=0; i<n; ++i) {
     282    func.addFunction(*funcs_[i]);
     283  }
     284  fitter.setFunction(func);
     285  fitter.setMaxIter(50+n*10);
     286  // Convergence criterium
     287  fitter.setCriteria(0.001);
     288 
     289  // Fit
     290  Vector<Float> sigma(x_.nelements());
     291  sigma = 1.0;
     292 
     293  parameters_.resize();
     294  parameters_ = fitter.fit(x_, y_, sigma, &m_);
     295  std::vector<float> ps;
     296  parameters_.tovector(ps);
     297  setParameters(ps);
     298  error_.resize();
     299  error_ = fitter.errors();
     300 
     301  chisquared_ = fitter.getChi2();
     302 
     303  residual_.resize();
     304  residual_ =  y_;
     305  fitter.residual(residual_,x_);
     306 
     307  // use fitter.residual(model=True) to get the model
     308  thefit_.resize(x_.nelements());
     309  fitter.residual(thefit_,x_,True);
     310  return true;
     311}
     312
     313
     314std::vector<float> SDFitter::evaluate(int whichComp) const
     315{ 
     316  std::vector<float> stlout;
     317  uInt idx = uInt(whichComp);
     318  Float y;
     319  if ( idx < funcs_.nelements() ) {
     320    for (uInt i=0; i<x_.nelements(); ++i) {
     321      y = (*funcs_[idx])(x_[i]);
     322      stlout.push_back(float(y));
     323    }
     324  }
     325  return stlout;
     326}
     327

  • trunk/src/SDFitter.h

    r125 r517  
    4444class SDFitter {
    4545public:
    46     SDFitter();
    47     virtual ~SDFitter();
    48     // allowed "gauss" and "poly". ncomp is either numvber of gaussions
    49     // or order of the polynomial
    50     bool setExpression(const std::string& expr, int ncomp=1);
    51     bool setData(std::vector<float> absc, std::vector<float> spec,
    52                  std::vector<bool> mask);
    53     bool setParameters(std::vector<float> params);
    54     bool setFixedParameters(std::vector<bool> fixed);
    55 
    56     std::vector<float> getResidual() const;
    57     std::vector<float> getFit() const;
    58     std::vector<float> getParameters() const;
    59     std::vector<bool> getFixedParameters() const;
    60 
    61     std::vector<float> getEstimate() const;
    62     std::vector<float> getErrors() const;
    63     float getChisquared() const;
    64     void reset();
    65     bool fit();
    66     bool computeEstimate();
    67     //std::vector<float> getEstimate() const;
     46  SDFitter();
     47  virtual ~SDFitter();
     48  // allowed "gauss" and "poly". ncomp is either numvber of gaussions
     49  // or order of the polynomial
     50  bool setExpression(const std::string& expr, int ncomp=1);
     51  bool setData(std::vector<float> absc, std::vector<float> spec,
     52               std::vector<bool> mask);
     53  bool setParameters(std::vector<float> params);
     54  bool setFixedParameters(std::vector<bool> fixed);
     55 
     56  std::vector<float> getResidual() const;
     57  std::vector<float> getFit() const;
     58  std::vector<float> getParameters() const;
     59  std::vector<bool> getFixedParameters() const;
     60 
     61  std::vector<float> getEstimate() const;
     62  std::vector<float> getErrors() const;
     63  float getChisquared() const;
     64  void reset();
     65  bool fit();
     66  bool computeEstimate();
     67 
     68  std::vector<float> evaluate(int whichComp) const;
    6869private:
    69     void clear();
    70     casa::Vector<casa::Float> x_;
    71     casa::Vector<casa::Float> y_;
    72     casa::Vector<casa::Bool> m_;
    73     casa::PtrBlock<casa::Function<casa::Float>* > funcs_;
    74     casa::CompoundFunction<casa::Float> cfunc_;
    75     //Bool estimateSet_;
    76     casa::Float chisquared_;
    77     casa::Vector<casa::Float> parameters_;
    78     casa::Vector<casa::Bool> fixedpar_;
    79 
    80     casa::Vector<casa::Float> error_;
    81     casa::Vector<casa::Float> thefit_;
    82     casa::Vector<casa::Float> residual_;
    83     casa::Vector<casa::Float> estimate_;
     70  void clear();
     71  casa::Vector<casa::Float> x_;
     72  casa::Vector<casa::Float> y_;
     73  casa::Vector<casa::Bool> m_;
     74  casa::PtrBlock<casa::Function<casa::Float>* > funcs_;
     75  //Bool estimateSet_;
     76  casa::Float chisquared_;
     77  casa::Vector<casa::Float> parameters_;
     78  casa::Vector<casa::Bool> fixedpar_;
     79 
     80  casa::Vector<casa::Float> error_;
     81  casa::Vector<casa::Float> thefit_;
     82  casa::Vector<casa::Float> residual_;
     83  casa::Vector<casa::Float> estimate_;
    8484};
    8585

  • trunk/src/python_SDFitter.cc

    r125 r517  
    5555        .def("reset", &SDFitter::reset)
    5656        .def("fit", &SDFitter::fit)
     57        .def("evaluate", &SDFitter::evaluate)
    5758      ;
    5859    };
Note: See TracChangeset for help on using the changeset viewer.