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 0 : }
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 :
|
