Line data    Source code

       1             : //# ProfileFit1D.h: Class to fit profiles to vectors
       2             : //# Copyright (C) 2004
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //#   $Id: ProfileFit1D.h 20299 2008-04-03 05:56:44Z gervandiepen $
      27             : 
      28             : #ifndef COMPONENTS_PROFILEFIT1D_H
      29             : #define COMPONENTS_PROFILEFIT1D_H
      30             : 
      31             : //# Includes
      32             : #include <casacore/casa/aips.h>
      33             : #include <casacore/casa/Arrays/Vector.h>
      34             : 
      35             : #include <components/SpectralComponents/SpectralList.h>
      36             : #include <components/SpectralComponents/SpectralFit.h>
      37             : 
      38             : #include <set>
      39             : 
      40             : namespace casa {
      41             : 
      42             : class SpectralElement;
      43             : 
      44             : // <summary>
      45             : // Fit spectral components to a casacore::Vector of data
      46             : // </summary>
      47             : 
      48             : // <use visibility=export>
      49             : 
      50             : // <reviewed reviewer="" date="" tests="tProfileFit1D.cc">
      51             : // </reviewed>
      52             : 
      53             : // <prerequisite>
      54             : //   <li> <linkto class="SpectralElement">SpectralElement</linkto> 
      55             : //   <li> <linkto class="SpectralList">SpectralList</linkto> 
      56             : //   <li> <linkto class="SpectralFit">SpectralFit</linkto> 
      57             : // </prerequisite>
      58             : 
      59             : // <synopsis> 
      60             : // Fit lists (held in class SpectralList) of SpectralElements to a 
      61             : // casacore::Vector of data.  Each SpectralElement can  be one from a variety
      62             : // of types.  The values of the parameters for each SpectralElement 
      63             : // provide the initial starting guesses for the fitting process.  
      64             : // Also, a SpectralElement object holds a mask indicating whether 
      65             : // a parameter should be held fixed or solved for.   After the 
      66             : // fitting is done, a new SpectralList holding SpectralElements with 
      67             : // the fitted parameters is created.  
      68             : //
      69             : // For all the functions that return a status casacore::Bool, true is good. If
      70             : // false is returned, an error message can be recovered with function
      71             : // <src>errorMessage</src>,  You should not proceed if false is returned.
      72             : // </synopsis> 
      73             : 
      74             : // <example>
      75             : // <srcblock>
      76             : // const casacore::uInt n = 512;
      77             : // casacore::Vector<casacore::Double> x(n);
      78             : // casacore::Vector<casacore::Double> y(n);
      79             : // casacore::Vector<casacore::Bool> m(n);
      80             : // 
      81             : // // Code to fill data vectors x,y,m
      82             : //
      83             : // ProfileFit1D<casacore::Float> fitter; 
      84             : // casacore::Bool ok = fitter.setData (x, y, m);
      85             : // ok = fitter.setGaussianElements (2);              
      86             : // ok = fitter.fit();       
      87             : // const SpectralList& fitList = fitter.getList(true);
      88             : // </srcblock>
      89             : // </example>
      90             : 
      91             : // <todo asof="2004/07/10">
      92             : //   <li> Add constraints
      93             : // </todo>
      94             : //
      95             : //
      96             : //  <note role=caution>
      97             : //  At the moment, because of templating limitations of the underlying
      98             : //  fitting classes, you must use template type Double.  This restriction
      99             : //  will be lifted in the future.
     100             : //  </note>
     101             : //
     102             : 
     103             : template <class T> class ProfileFit1D 
     104             : {
     105             : public:
     106             :     // Constructor
     107             :     ProfileFit1D();
     108             : 
     109             :     // Destructor
     110             :     ~ProfileFit1D();
     111             : 
     112             :     // Copy constructor.  Uses copy semantics.
     113             :     ProfileFit1D(const ProfileFit1D& other);
     114             : 
     115             :     // Assignment operator. Uses copy semantics.
     116             :     ProfileFit1D& operator=(const ProfileFit1D& other);
     117             : 
     118             :     // Set abcissa, ordinate, mask and weights. A true mask value means the 
     119             :     // data are good. If you don't specify the weights vector, all weights 
     120             :     // are assumed to be unity.  If you don't specify a mask it will be
     121             :     // created as all good.
     122             :     // Status is returned, if false, error message can be recovered with <src>errorMessage</src>
     123             :     // <group>
     124             :     casacore::Bool setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask,
     125             :                   const casacore::Vector<casacore::Double>& weight);
     126             :     casacore::Bool setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y, const casacore::Vector<casacore::Bool>& mask);
     127             :     casacore::Bool setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y);
     128             :     // </group>
     129             : 
     130             :     // Set a SpectralList of SpectralElements to fit for.
     131             :     // The SpectralElements in the list hold the
     132             :     // initial estimates and must reflect the abcissa and ordinate units.  
     133             :     // They also contain the information about whether
     134             :     // specific parameters are to be held fixed or allowed to vary in
     135             :     // the fitting process.
     136             :     // You can recover the list of elements with function getList.
     137             :     void setElements (const SpectralList& list);
     138             : 
     139             :     // Set a SpectralList of Gaussian SpectralElements to fit for.  
     140             :     // The initial estimates for the Gaussians will be automatically determined.
     141             :     // All of the parameters created by this function will be solved for
     142             :     // by default. You can recover the list of elements with function getList.
     143             :     // Status is returned, if false, error message can be 
     144             :     // recovered with <src>errorMessage</src>
     145             :     casacore::Bool setGaussianElements (casacore::uInt nGauss);
     146             : 
     147             :     // Add new SpectralElement(s) to the SpectralList (can be empty)
     148             :     // of SpectralElements to be fit for.  
     149             :     //<group>
     150             :     void addElement (const SpectralElement& el);
     151             :     void addElements (const SpectralList& list);
     152             :     // </group>
     153             : 
     154             :     // Clear the SpectralList of elements to be fit for
     155             :     void clearList ();
     156             : 
     157             :     // Set abscissa range mask.  You can specify a number of ranges
     158             :     // via a vector of start indices (or X values) and a vector of end
     159             :     // indices (or X values).   When argument insideIsGood is true,
     160             :     // a mask will be created which
     161             :     // will be true (good) inside the ranges and bad (false)
     162             :     // outside of those ranges.  When argument insideIsGood is false,
     163             :     // the mask will be false (bad) inside the ranges and 
     164             :     // true (good) outside of those ranges. When the data are fit, a total
     165             :     // mask is formed combining (via a logical AND) the 
     166             :     // data mask (setData) and this range mask.
     167             :     // Status is returned, if false, error message can be recovered with <src>errorMessage</src>
     168             :     // In the single set version, the values in the set indicate the pixels to set the mask for,
     169             :     // ie no ranges, just specific pixels are to be provided. In this case, specified values
     170             :     // which are greater than or equal to the number of pixels are tacitly ignored.
     171             :     // <group>
     172             :     casacore::Bool setXRangeMask (const casacore::Vector<casacore::uInt>& startIndex,
     173             :                        const casacore::Vector<casacore::uInt>& endIndex,
     174             :                        casacore::Bool insideIsGood=true);
     175             :     casacore::Bool setXRangeMask (const casacore::Vector<T>& startIndex,
     176             :                        const casacore::Vector<T>& endIndex,
     177             :                        casacore::Bool insideIsGood=true);
     178             : 
     179             :     casacore::Bool setXMask(const std::set<casacore::uInt>& indices, casacore::Bool specifiedPixelsAreGood);
     180             :     // </group>
     181             : 
     182             :     // Recover masks.  These are the data mask (setData) the range
     183             :     // mask (setRangeMask may be length zero) and the total 
     184             :     // mask combining the two.
     185             :     // <group>
     186      529098 :     casacore::Vector<casacore::Bool> getDataMask() const {return itsDataMask;};
     187             :     casacore::Vector<casacore::Bool> getRangeMask() const {return itsRangeMask;}
     188             :     casacore::Vector<casacore::Bool> getTotalMask() const {return makeTotalMask();};
     189             :     // </group>
     190             : 
     191             :     // Do the fit and return status.  Returns convergence status.
     192             :     // Error conditions in the solution process will generate
     193             :     // an casacore::AipsError exception and you should catch these yourself.
     194             :     casacore::Bool fit ();
     195             : 
     196             :     // Get Chi Squared of fit
     197           0 :     casacore::Double getChiSquared () const {return itsFitter.chiSq();}
     198             : 
     199             :     // Get number of iterations for last fit
     200       16755 :     casacore::Double getNumberIterations () const {return itsFitter.nIterations();}
     201             : 
     202             :     // Recover the list of elements.  You can get the elements
     203             :     // as initially estimated (fit=false), or after fitting 
     204             :     // (fit=true).  In the latter case, the SpectralElements
     205             :     // hold the parameters and errors of the fit.
     206             :     const SpectralList& getList (casacore::Bool fit=true) const;
     207             : 
     208             :     // Recover vectors for the estimate, fit and residual.
     209             :     // If you don't specify which element, all elements are included
     210             :     // If the Vectors are returned with zero length, it means an error
     211             :     // condition exists (e.g. asking for fit before you do one). In this
     212             :     // case an error message can be recovered with function <src>errorMessage</src>.
     213             :     //<group>
     214             :     casacore::Vector<T> getEstimate (casacore::Int which=-1) const;
     215             :     casacore::Vector<T> getFit (casacore::Int which=-1) const;
     216             :     casacore::Vector<T> getResidual (casacore::Int which=-1, casacore::Bool fit=true)  const;
     217             :     //</group>
     218             : 
     219             :     // Recover the error message
     220           0 :     casacore::String errorMessage () const {return itsError;};
     221             : 
     222             : private:
     223             :    casacore::Vector<casacore::Double> itsX;                    // Abcissa (really should not be type T)
     224             :    casacore::Vector<T> itsY;                         // Ordinate 
     225             :    casacore::Vector<casacore::Double> itsWeight;               // Weights 
     226             :    casacore::Vector<casacore::Bool> itsDataMask;               // casacore::Data mask
     227             :    casacore::Vector<casacore::Bool> itsRangeMask;              // Mask associated with ranges
     228             : //
     229             :    SpectralList itsList;                   // casacore::List of elements to fit for
     230             : //
     231             :    SpectralFit itsFitter;                  // Fitter
     232             :    mutable casacore::String itsError;                // Error message
     233             : 
     234             : // Functions
     235             :    casacore::Vector<casacore::Bool> makeTotalMask() const;
     236             :    SpectralList getSubsetList (const SpectralList& list, casacore::Int which)  const;
     237             :    void checkType() const;
     238             :    void copy(const ProfileFit1D<T>& other);
     239             : };
     240             : 
     241             : } //#End casa namespace
     242             : #ifndef CASACORE_NO_AUTO_TEMPLATES
     243             : #include <components/SpectralComponents/ProfileFit1D.tcc>
     244             : #endif //# CASACORE_NO_AUTO_TEMPLATES
     245             : #endif