LCOV - code coverage report
Current view: top level - mstransform/TVI - DenoisingLib.h (source / functions) Hit Total Coverage
Test: casa_coverage.info Lines: 29 50 58.0 %
Date: 2023-10-25 08:47:59 Functions: 8 13 61.5 %

          Line data    Source code
       1             : //# DenoisingLib.h: This file contains the interface definition of the MSTransformManager class.
       2             : //#
       3             : //#  CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
       4             : //#  Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
       5             : //#  Copyright (C) European Southern Observatory, 2011, All rights reserved.
       6             : //#
       7             : //#  This library is free software; you can redistribute it and/or
       8             : //#  modify it under the terms of the GNU Lesser General Public
       9             : //#  License as published by the Free software Foundation; either
      10             : //#  version 2.1 of the License, or (at your option) any later version.
      11             : //#
      12             : //#  This library is distributed in the hope that it will be useful,
      13             : //#  but WITHOUT ANY WARRANTY, without even the implied warranty of
      14             : //#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      15             : //#  Lesser General Public License for more details.
      16             : //#
      17             : //#  You should have received a copy of the GNU Lesser General Public
      18             : //#  License along with this library; if not, write to the Free Software
      19             : //#  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
      20             : //#  MA 02111-1307  USA
      21             : //# $Id: $
      22             : 
      23             : #ifndef DenoisingLib_H_
      24             : #define DenoisingLib_H_
      25             : 
      26             : // casacore containers
      27             : #include <casacore/casa/Arrays/Matrix.h>
      28             : #include <casacore/casa/Arrays/Vector.h>
      29             : 
      30             : // logger
      31             : #include <casacore/casa/Logging/LogIO.h>
      32             : 
      33             : // GSL includes
      34             : #include <gsl/gsl_multifit.h>
      35             : 
      36             : using namespace casacore;
      37             : 
      38             : namespace casa { //# NAMESPACE CASA - BEGIN
      39             : 
      40             : namespace denoising { //# NAMESPACE DENOISING - BEGIN
      41             : 
      42             : void GslVectorWrap(Vector<Double> casa_vector, gsl_vector &wrapper);
      43             : void GslMatrixWrap(Matrix<Double> &casa_matrix, gsl_matrix &wrapper, size_t ncols=0);
      44             : 
      45             : //////////////////////////////////////////////////////////////////////////
      46             : // GslLinearModelBase class
      47             : //////////////////////////////////////////////////////////////////////////
      48             : 
      49             : template<class T> class GslLinearModelBase
      50             : {
      51             : 
      52             : public:
      53             : 
      54        1674 :         GslLinearModelBase(size_t ndata, size_t ncomponents)
      55        1674 :         {
      56        1674 :                 ndata_p = ndata;
      57        1674 :                 ncomponents_p = ncomponents;
      58        1674 :                 model_p.resize(ncomponents_p,ndata_p, false);
      59        1674 :         }
      60             : 
      61        3354 :         size_t ndata() {return ndata_p;}
      62        6776 :         size_t ncomponents() {return ncomponents_p;}
      63       78478 :         Matrix<T>& getModelMatrix(){return model_p;}
      64           0 :         Vector<T> getModelAt(size_t pos){return model_p.column(pos);}
      65             : 
      66             : protected:
      67             : 
      68             :         Matrix<T> model_p;
      69             :         size_t ndata_p;
      70             :         size_t ncomponents_p;
      71             : };
      72             : 
      73             : //////////////////////////////////////////////////////////////////////////
      74             : // GslPolynomialModel class
      75             : //////////////////////////////////////////////////////////////////////////
      76             : 
      77             : template<class T> class GslPolynomialModel : public GslLinearModelBase<T>
      78             : {
      79             :         using GslLinearModelBase<T>::model_p;
      80             :         using GslLinearModelBase<T>::ndata_p;
      81             :         using GslLinearModelBase<T>::ncomponents_p;
      82             : 
      83             : public:
      84             : 
      85             :         GslPolynomialModel(size_t ndata, size_t order) :
      86             :                 GslLinearModelBase<T>(ndata,order+1)
      87             :         {
      88             :                 // Initialize model
      89             :                 model_p.row(0) = 1.0;
      90             : 
      91             :                 // Populate linear components
      92             :                 if (order > 0)
      93             :                 {
      94             :                         Vector<T> linearComponent;
      95             :                         linearComponent.reference(model_p.row(1));
      96             :                         indgen(linearComponent,-1.0,2.0/(ndata_p-1));
      97             :                 }
      98             : 
      99             :                 // Populate higher orders
     100             :                 for (size_t model_idx=2;model_idx<ncomponents_p;model_idx++)
     101             :                 {
     102             :                         for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
     103             :                         {
     104             :                                 model_p(model_idx,data_idx) = pow(model_p(1,data_idx),model_idx);
     105             :                         }
     106             :                 }
     107             :         }
     108             : 
     109        1674 :         GslPolynomialModel(const Vector<T> &data_x, size_t order) :
     110        1674 :                 GslLinearModelBase<T>(data_x.size(),order+1)
     111             :         {
     112             :                 // Calculate scale
     113             :                 T min,max,mid,scale;
     114        1674 :                 minMax(min,max,data_x);
     115        1674 :                 mid = 0.5*(min+max);
     116        1674 :                 scale = 1.0 / (min-mid);
     117             : 
     118             :                 // Populate linear components
     119        1674 :                 model_p.row(0) = 1.0;
     120        1674 :                 if (order > 0) model_p.row(1) = scale*(data_x-mid);
     121             : 
     122             :                 // Populate higher orders
     123        1828 :                 for (size_t model_idx=2;model_idx<ncomponents_p;model_idx++)
     124             :                 {
     125       25285 :                         for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
     126             :                         {
     127       25131 :                                 model_p(model_idx,data_idx) = pow(model_p(1,data_idx),model_idx);
     128             :                         }
     129             :                 }
     130        1674 :         }
     131             : 
     132        3422 :         Vector<Float>& getLinearComponentFloat()
     133             :         {
     134        3422 :                 if (linear_component_float_p.size() != ndata_p) initLinearComponentFloat();
     135        3422 :                 return linear_component_float_p;
     136             :         }
     137             : 
     138             : private:
     139             : 
     140        1674 :         void initLinearComponentFloat()
     141             :         {
     142        1674 :                 linear_component_float_p.resize(ndata_p);
     143     1245338 :                 for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
     144             :                 {
     145     1243664 :                         linear_component_float_p(data_idx) = model_p(1,data_idx);
     146             :                 }
     147        1674 :         }
     148             : 
     149             :         Vector<Float> linear_component_float_p; // Float-compatibility
     150             : 
     151             : };
     152             : 
     153             : //////////////////////////////////////////////////////////////////////////
     154             : // GslMultifitLinearBase class
     155             : //////////////////////////////////////////////////////////////////////////
     156             : 
     157             : class GslMultifitLinearBase
     158             : {
     159             : 
     160             : public:
     161             : 
     162             :         GslMultifitLinearBase();
     163             :         GslMultifitLinearBase(GslLinearModelBase<Double> &model);
     164             : 
     165             :         ~GslMultifitLinearBase();
     166             : 
     167             :         void resetModel(GslLinearModelBase<Double> &model);
     168             : 
     169             :         void resetNComponents(size_t ncomponents);
     170             : 
     171             :         void setDebug(Bool debug) {debug_p = debug;};
     172             : 
     173             :         std::pair<std::vector<Complex>, Complex> calcFitCoeff(Vector<Complex> &data);
     174           0 :         template<class T> std::pair<std::vector<T>, double> calcFitCoeff(Vector<T> &data)
     175             :         {
     176             :                 // Set data
     177           0 :                 setData(data);
     178             : 
     179             :                 // Call fit method to calculate coefficients
     180           0 :                 double chisq = calcFitCoeffCore(data_p.column(0),gsl_coeff_real_p);
     181             : 
     182             :                 // Convert GSL vector into vector
     183           0 :                 std::vector<T> coeffCASA(ncomponents_p);
     184           0 :                 for (size_t coeff_idx=0;coeff_idx<ncomponents_p;coeff_idx++)
     185             :                 {
     186           0 :                         coeffCASA[coeff_idx] = gsl_vector_get(gsl_coeff_real_p, coeff_idx);
     187             :                 }
     188             : 
     189           0 :                 return std::make_pair(coeffCASA, chisq);
     190             :         }
     191             : 
     192             :         void getFitCoeff(Vector<Complex> &coeff);
     193             :         template<class T> void getFitCoeff(Vector<T> &coeff)
     194             :         {
     195             :                 coeff.resize(ncomponents_p, false);
     196             :                 for (size_t model_idx=0;model_idx<ncomponents_p;model_idx++)
     197             :                 {
     198             :                         coeff(model_idx) = gsl_vector_get(gsl_coeff_real_p,model_idx);
     199             :                 }
     200             : 
     201             :                 return;
     202             :         }
     203             : 
     204             :         void calcFitModelStd(Vector<Complex> &model,Vector<Complex> &std);
     205             :         template<class T> void calcFitModelStd(   Vector<T> &model, Vector<T> &std)
     206             :         {
     207             :                 calcFitModelStdCore(model,std,gsl_coeff_real_p);
     208             : 
     209             :                 return;
     210             :         }
     211             : 
     212             :         void calcFitModel(Vector<Complex> &model);
     213           0 :         template<class T> void calcFitModel(Vector<T> &model)
     214             :         {
     215           0 :                 calcFitModelCore(model,gsl_coeff_real_p);
     216             : 
     217           0 :                 return;
     218             :         }
     219             : 
     220             : 
     221             : protected:
     222             : 
     223             :         void allocGslResources();
     224             :         void freeGslResources();
     225             : 
     226             :         virtual void setModel(GslLinearModelBase<Double> &model);
     227             : 
     228             :         void setData(Vector<Float> &data);
     229             :         void setData(Vector<Double> &data);
     230             :         void setData(Vector<Complex> &data);
     231             : 
     232             :         virtual double calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff);
     233             : 
     234             :         template<class T> void calcFitModelStdCore(       Vector<T> &model, Vector<T> &std, gsl_vector *coeff)
     235             :         {
     236             :                 // Get imag coefficients
     237             :                 gsl_vector xGSL;
     238             :                 double y, yerr;
     239             :                 for (size_t data_idx=0;data_idx<ndata_p;data_idx++)
     240             :                 {
     241             :                         Vector<Double> xCASA = model_p->getModelAt(data_idx);
     242             :                         if (xCASA.size() != ncomponents_p) xCASA.resize(ncomponents_p,True);
     243             :                         GslVectorWrap(xCASA,xGSL);
     244             : 
     245             :                         y = 0;
     246             :                         yerr = 0;
     247             :                         errno_p = gsl_multifit_linear_est (&xGSL, coeff, gsl_covariance_p, &y, &yerr);
     248             : 
     249             :                         if (model.size() > 0) model(data_idx) = y;
     250             :                         if (std.size() > 0 ) std(data_idx) = yerr;
     251             :                 }
     252             : 
     253             :                 return;
     254             :         }
     255             : 
     256           0 :         template<class T> void calcFitModelCore(Vector<T> &model, gsl_vector *coeff)
     257             :         {
     258             :                 Double coeff_i;
     259           0 :                 Matrix<Double> &modelMatrix = model_p->getModelMatrix();
     260             : 
     261           0 :                 coeff_i = gsl_vector_get(coeff,0);
     262           0 :                 for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
     263             :                 {
     264           0 :                         model(data_idx) = coeff_i*modelMatrix(0,data_idx);
     265             :                 }
     266             : 
     267           0 :                 for (size_t model_idx=0;model_idx<ncomponents_p;model_idx++)
     268             :                 {
     269           0 :                         coeff_i = gsl_vector_get(coeff,model_idx);
     270           0 :                         for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
     271             :                         {
     272           0 :                                 model(data_idx) = coeff_i*modelMatrix(model_idx,data_idx);
     273             :                         }
     274             :                 }
     275             : 
     276           0 :                 return;
     277             :         }
     278             : 
     279             :         // Model
     280             :         size_t ndata_p;
     281             :         size_t ncomponents_p;
     282             :         size_t max_ncomponents_p;
     283             :         gsl_matrix gsl_model_p;
     284             :         GslLinearModelBase<Double> *model_p;
     285             : 
     286             :         // GSL Resources
     287             :         gsl_vector *gsl_coeff_real_p;
     288             :         gsl_vector *gsl_coeff_imag_p;
     289             :         gsl_matrix *gsl_covariance_p;
     290             :         gsl_multifit_linear_workspace *gsl_workspace_p;
     291             : 
     292             :         // Data
     293             :         Matrix<Double> data_p;
     294             : 
     295             :         // Fit result
     296             :         int errno_p;
     297             :         double chisq_p;
     298             : 
     299             :         // Misc
     300             :         Bool debug_p;
     301             : };
     302             : 
     303             : //////////////////////////////////////////////////////////////////////////
     304             : // GslMultifitWeightedLinear class
     305             : //////////////////////////////////////////////////////////////////////////
     306             : 
     307             : class GslMultifitWeightedLinear : public GslMultifitLinearBase
     308             : {
     309             : 
     310             : public:
     311             : 
     312             :         GslMultifitWeightedLinear();
     313             :         GslMultifitWeightedLinear(GslLinearModelBase<Double> &model);
     314             :         ~GslMultifitWeightedLinear();
     315             : 
     316             :         void setWeights(Vector<Float> &weights);
     317             :         void setFlags(Vector<Bool> &flags,Bool goodIsTrue=True);
     318             :         void setWeightsAndFlags(Vector<Float> &weights, Vector<Bool> &flags, Bool goodIsTrue=True);
     319             : 
     320             : protected:
     321             : 
     322             :         void setModel(GslLinearModelBase<Double> &model);
     323             :         virtual double calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff);
     324             : 
     325             :         // Weights
     326             :         Vector<Double> weights_p;
     327             :         gsl_vector gls_weights_p;
     328             : };
     329             : 
     330             : //////////////////////////////////////////////////////////////////////////
     331             : // Iteratively Reweighted Least Squares class
     332             : //////////////////////////////////////////////////////////////////////////
     333             : 
     334             : class IterativelyReweightedLeastSquares : public GslMultifitWeightedLinear
     335             : {
     336             : 
     337             : public:
     338             : 
     339             :         IterativelyReweightedLeastSquares();
     340             :         IterativelyReweightedLeastSquares(GslLinearModelBase<Double> &model,size_t nIter);
     341             :         ~IterativelyReweightedLeastSquares();
     342             : 
     343        3354 :         void setNIter(size_t nIter) {nIter_p = nIter;};
     344             : 
     345             :         virtual double calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff);
     346             :         virtual void updateWeights(Vector<Double> &data, Vector<Double> &model,  Vector<Double> &weights);
     347             : 
     348             : protected:
     349             : 
     350             :         size_t nIter_p;
     351             : 
     352             : };
     353             : 
     354             : } //# NAMESPACE DENOISING - END
     355             : 
     356             : } //# NAMESPACE CASA - END
     357             : 
     358             : 
     359             : #endif /* DenoisingLib_H_ */
     360             : 

Generated by: LCOV version 1.16