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 : #include <mstransform/TVI/DenoisingLib.h>
24 :
25 : namespace casa { //# NAMESPACE CASA - BEGIN
26 :
27 : namespace denoising { //# NAMESPACE DENOISING - BEGIN
28 :
29 : using std::vector;
30 :
31 : // -----------------------------------------------------------------------
32 : // Wrap CASA Vector with a gsl_vector structure
33 : // -----------------------------------------------------------------------
34 139558 : void GslVectorWrap(Vector<Double> casa_vector, gsl_vector &wrapper)
35 : {
36 139558 : wrapper.size = casa_vector.size();
37 139558 : wrapper.stride = casa_vector.steps()(0);
38 139558 : wrapper.data = casa_vector.data();
39 139558 : wrapper.owner = false;
40 :
41 139558 : return;
42 : }
43 :
44 : // -----------------------------------------------------------------------
45 : //
46 : // Wrap CASA Matrix with a gsl_matrix structure
47 : //
48 : // GSL Matrices are stored in row-major order, meaning that
49 : // each row of elements forms a contiguous block in memory.
50 : // This is the standard “C-language ordering” of two-dimensional arrays.
51 : //
52 : // CASA Matrices are however stored in column-major order
53 : // so the elements of each column forms a contiguous block in memory.
54 : //
55 : // Therefore we have to swap rows-cols in order to match.
56 : //
57 : // Note that FORTRAN stores arrays in column-major order.
58 : // -----------------------------------------------------------------------
59 6954 : void GslMatrixWrap(Matrix<Double> &casa_matrix, gsl_matrix &wrapper, size_t ncols)
60 : {
61 6954 : ThrowIf (not casa_matrix.contiguousStorage(),
62 : "Cannot map a non-contiguous CASA matrix to GSL matrix");
63 :
64 6954 : wrapper.size1 = casa_matrix.ncolumn(); // Number of rows
65 6954 : wrapper.size2 = ncols > 0? ncols : casa_matrix.nrow(); // Number of columns
66 6954 : wrapper.tda = casa_matrix.nrow();
67 6954 : wrapper.block = nullptr;
68 6954 : wrapper.data = casa_matrix.data();
69 6954 : wrapper.owner = false;
70 :
71 6954 : return;
72 : }
73 :
74 : //////////////////////////////////////////////////////////////////////////
75 : // GslMultifitLinearBase class
76 : //////////////////////////////////////////////////////////////////////////
77 :
78 : // -----------------------------------------------------------------------
79 : //
80 : // -----------------------------------------------------------------------
81 3354 : GslMultifitLinearBase::GslMultifitLinearBase()
82 : {
83 3354 : model_p = nullptr;
84 3354 : ndata_p = 0;
85 3354 : ncomponents_p = 0;
86 3354 : max_ncomponents_p = 0;
87 :
88 3354 : gsl_coeff_real_p = nullptr;
89 3354 : gsl_coeff_imag_p = nullptr;
90 3354 : gsl_covariance_p = nullptr;
91 3354 : gsl_workspace_p = nullptr;
92 :
93 3354 : errno_p = 0;
94 3354 : chisq_p = 0;
95 :
96 3354 : debug_p = false;
97 3354 : }
98 :
99 : // -----------------------------------------------------------------------
100 : //
101 : // -----------------------------------------------------------------------
102 0 : GslMultifitLinearBase::GslMultifitLinearBase(GslLinearModelBase<Double> &model)
103 : {
104 0 : setModel(model);
105 :
106 0 : errno_p = 0;
107 0 : chisq_p = 0;
108 0 : }
109 :
110 : // -----------------------------------------------------------------------
111 : //
112 : // -----------------------------------------------------------------------
113 3354 : GslMultifitLinearBase::~GslMultifitLinearBase()
114 : {
115 3354 : freeGslResources();
116 3354 : }
117 :
118 : // -----------------------------------------------------------------------
119 : //
120 : // -----------------------------------------------------------------------
121 6954 : void GslMultifitLinearBase::allocGslResources()
122 : {
123 6954 : gsl_covariance_p = gsl_matrix_alloc (ncomponents_p, ncomponents_p);
124 6954 : gsl_workspace_p = gsl_multifit_linear_alloc (ndata_p, ncomponents_p);
125 6954 : gsl_coeff_real_p = gsl_vector_alloc(ncomponents_p);
126 6954 : gsl_coeff_imag_p = gsl_vector_alloc(ncomponents_p);
127 6954 : }
128 :
129 : // -----------------------------------------------------------------------
130 : //
131 : // -----------------------------------------------------------------------
132 10308 : void GslMultifitLinearBase::freeGslResources()
133 : {
134 10308 : if (gsl_covariance_p != nullptr) gsl_matrix_free (gsl_covariance_p);
135 10308 : if (gsl_workspace_p != nullptr) gsl_multifit_linear_free (gsl_workspace_p);
136 10308 : if (gsl_coeff_real_p != nullptr) gsl_vector_free (gsl_coeff_real_p);
137 10308 : if (gsl_coeff_imag_p != nullptr) gsl_vector_free (gsl_coeff_imag_p);
138 10308 : }
139 :
140 : // -----------------------------------------------------------------------
141 : //
142 : // -----------------------------------------------------------------------
143 3354 : void GslMultifitLinearBase::setModel(GslLinearModelBase<Double> &model)
144 : {
145 3354 : model_p = &model;
146 3354 : ndata_p = model_p->ndata();
147 3354 : max_ncomponents_p = model_p->ncomponents();
148 :
149 3354 : ncomponents_p = max_ncomponents_p;
150 3354 : GslMatrixWrap(model_p->getModelMatrix(),gsl_model_p);
151 :
152 3354 : data_p.resize(ndata_p, 1, false);
153 :
154 3354 : allocGslResources();
155 3354 : }
156 :
157 : // -----------------------------------------------------------------------
158 : //
159 : // -----------------------------------------------------------------------
160 3600 : void GslMultifitLinearBase::resetNComponents(size_t ncomponents)
161 : {
162 3600 : ThrowIf (ncomponents > max_ncomponents_p,
163 : "Maximum number of components is " + String::toString(max_ncomponents_p));
164 :
165 3600 : ncomponents_p = ncomponents;
166 3600 : GslMatrixWrap(model_p->getModelMatrix(),gsl_model_p,ncomponents_p);
167 :
168 3600 : freeGslResources();
169 3600 : allocGslResources();
170 3600 : }
171 :
172 : // -----------------------------------------------------------------------
173 : //
174 : // -----------------------------------------------------------------------
175 3354 : void GslMultifitLinearBase::resetModel(GslLinearModelBase<Double> &model)
176 : {
177 3354 : freeGslResources();
178 3354 : setModel(model);
179 3354 : }
180 :
181 : // -----------------------------------------------------------------------
182 : //
183 : // -----------------------------------------------------------------------
184 0 : void GslMultifitLinearBase::setData(Vector<Double> &data)
185 : {
186 0 : ThrowIf (data.size() != ndata_p,"Size of data does not match model");
187 :
188 0 : data_p.column(0).reference(data);
189 0 : }
190 :
191 : // -----------------------------------------------------------------------
192 : //
193 : // -----------------------------------------------------------------------
194 0 : void GslMultifitLinearBase::setData(Vector<Float> &data)
195 : {
196 0 : ThrowIf (data.size() != ndata_p,"Size of data does not match model");
197 :
198 0 : for (size_t idx=0;idx<ndata_p;idx++)
199 : {
200 0 : data_p(idx,0) = data(idx);
201 : }
202 0 : }
203 :
204 : // -----------------------------------------------------------------------
205 : //
206 : // -----------------------------------------------------------------------
207 68102 : void GslMultifitLinearBase::setData(Vector<Complex> &data)
208 : {
209 68102 : ThrowIf (data.size() != ndata_p,"Size of data does not match model");
210 :
211 68102 : if (data_p.ncolumn() != 2) data_p.resize(ndata_p, 2, false);
212 :
213 19574290 : for (size_t idx=0;idx<ndata_p;idx++)
214 : {
215 19506188 : data_p(idx,0) = real(data(idx));
216 19506188 : data_p(idx,1) = imag(data(idx));
217 : }
218 68102 : }
219 :
220 : // -----------------------------------------------------------------------
221 : //
222 : // Perform least-squares fits to a general linear model, y = X c where
223 : // y is a vector of n observations, X is an n by p matrix of predictor variables,
224 : // and the elements of the vector c are the p unknown best-fit parameters which
225 : // are to be estimated.
226 : //
227 : // NOTE: gsl_multifit_linear expects that the model matrix is organized as follows
228 : // X = [ 1 , x_1 , x_1^2 , ..., x_1^order;
229 : // 1 , x_2 , x_2^2 , ..., x_2^order;
230 : // 1 , x_3 , x_3^2 , ..., x_3^order;
231 : // ... , ... , ... , ..., ...
232 : // 1 , x_n , x_n^2 , ..., x_n^order]
233 : //
234 : // -----------------------------------------------------------------------
235 0 : double GslMultifitLinearBase::calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff)
236 : {
237 : // Wrap data vector in a gsl_vector
238 : gsl_vector data_gsl;
239 0 : GslVectorWrap(data,data_gsl);
240 :
241 : // Perform coeff calculation
242 0 : errno_p = gsl_multifit_linear (&gsl_model_p, &data_gsl,
243 : coeff, gsl_covariance_p, &chisq_p, gsl_workspace_p);
244 :
245 0 : return chisq_p;
246 : }
247 :
248 : // -----------------------------------------------------------------------
249 : //
250 : // -----------------------------------------------------------------------
251 : std::pair<vector<Complex>, Complex>
252 68102 : GslMultifitLinearBase::calcFitCoeff(Vector<Complex> &data)
253 : {
254 : // Set data
255 68102 : setData(data);
256 :
257 : // Call fit method to calculate real/imag coefficients
258 68102 : const auto chisqReal = calcFitCoeffCore(data_p.column(0), gsl_coeff_real_p);
259 68102 : const auto chisqImag = calcFitCoeffCore(data_p.column(1), gsl_coeff_imag_p);
260 :
261 : // Get imag coefficients
262 68102 : vector<Complex> fitCoeff(ncomponents_p);
263 155888 : for (size_t coeff_idx=0;coeff_idx<ncomponents_p;coeff_idx++)
264 : {
265 87786 : fitCoeff[coeff_idx] = Complex( gsl_vector_get(gsl_coeff_real_p,coeff_idx),
266 87786 : gsl_vector_get(gsl_coeff_imag_p,coeff_idx));
267 : }
268 :
269 136204 : return make_pair(fitCoeff, Complex(chisqReal, chisqImag));
270 68102 : }
271 :
272 : // -----------------------------------------------------------------------
273 : //
274 : // -----------------------------------------------------------------------
275 0 : void GslMultifitLinearBase::calcFitModelStd(Vector<Complex> &model,Vector<Complex> &std)
276 : {
277 : // Get imag coefficients
278 : gsl_vector xGSL;
279 : double y_real, y_imag, yerr_real, yerr_imag;
280 0 : for (size_t data_idx=0;data_idx<ndata_p;data_idx++)
281 : {
282 0 : Vector<Double> xCASA = model_p->getModelAt(data_idx);
283 0 : if (xCASA.size() != ncomponents_p) xCASA.resize(ncomponents_p,True);
284 0 : GslVectorWrap(xCASA,xGSL);
285 :
286 0 : y_real = 0;
287 0 : yerr_real = 0;
288 0 : errno_p = gsl_multifit_linear_est (&xGSL, gsl_coeff_real_p, gsl_covariance_p, &y_real, &yerr_real);
289 :
290 0 : y_imag = 0;
291 0 : yerr_imag = 0;
292 0 : errno_p = gsl_multifit_linear_est (&xGSL, gsl_coeff_imag_p, gsl_covariance_p, &y_imag, &yerr_imag);
293 :
294 0 : if (model.size() > 0) model(data_idx) = Complex(y_real,y_imag);
295 0 : if (std.size() > 0 ) std(data_idx) = Complex(yerr_real,yerr_imag);
296 0 : }
297 :
298 0 : return;
299 : }
300 :
301 : // -----------------------------------------------------------------------
302 : //
303 : // -----------------------------------------------------------------------
304 68102 : void GslMultifitLinearBase::calcFitModel(Vector<Complex> &model)
305 : {
306 68102 : Complex coeff_i;
307 68102 : Matrix<Double> &modelMatrix = model_p->getModelMatrix();
308 :
309 68102 : coeff_i = Complex(gsl_vector_get(gsl_coeff_real_p,0),
310 68102 : gsl_vector_get(gsl_coeff_imag_p,0));
311 19574290 : for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
312 : {
313 19506188 : model(data_idx) = coeff_i*modelMatrix(0,data_idx);
314 : }
315 :
316 87786 : for (size_t model_idx=1;model_idx<ncomponents_p;model_idx++)
317 : {
318 19684 : coeff_i = Complex(gsl_vector_get(gsl_coeff_real_p,model_idx),
319 19684 : gsl_vector_get(gsl_coeff_imag_p,model_idx));
320 3946340 : for (size_t data_idx=0; data_idx<ndata_p; data_idx++)
321 : {
322 3926656 : model(data_idx) += coeff_i*modelMatrix(model_idx,data_idx);
323 : }
324 : }
325 :
326 136204 : return;
327 : }
328 :
329 : //////////////////////////////////////////////////////////////////////////
330 : // GslMultifitWeightedLinear class
331 : //////////////////////////////////////////////////////////////////////////
332 :
333 : // -----------------------------------------------------------------------
334 : //
335 : // -----------------------------------------------------------------------
336 3354 : GslMultifitWeightedLinear::GslMultifitWeightedLinear() :
337 3354 : GslMultifitLinearBase()
338 : {
339 :
340 3354 : }
341 :
342 : // -----------------------------------------------------------------------
343 : //
344 : // -----------------------------------------------------------------------
345 0 : GslMultifitWeightedLinear::GslMultifitWeightedLinear(GslLinearModelBase<Double> &model) :
346 0 : GslMultifitLinearBase()
347 : {
348 0 : setModel(model);
349 0 : }
350 :
351 : // -----------------------------------------------------------------------
352 : //
353 : // -----------------------------------------------------------------------
354 3354 : GslMultifitWeightedLinear::~GslMultifitWeightedLinear()
355 : {
356 :
357 3354 : }
358 :
359 : // -----------------------------------------------------------------------
360 : //
361 : // -----------------------------------------------------------------------
362 3354 : void GslMultifitWeightedLinear::setModel(GslLinearModelBase<Double> &model)
363 : {
364 3354 : GslMultifitLinearBase::setModel(model);
365 3354 : weights_p.resize(ndata_p, false);
366 3354 : GslVectorWrap(weights_p,gls_weights_p);
367 3354 : }
368 :
369 : // -----------------------------------------------------------------------
370 : //
371 : // -----------------------------------------------------------------------
372 0 : void GslMultifitWeightedLinear::setWeights(Vector<Float> &weights)
373 : {
374 : // Dim check
375 0 : ThrowIf (weights.size() != ndata_p,"Size of weights does not match model");
376 :
377 : // Fill in
378 0 : for (size_t idx=0;idx<ndata_p;idx++)
379 : {
380 0 : weights_p(idx) = weights(idx);
381 : }
382 0 : }
383 :
384 : // -----------------------------------------------------------------------
385 : //
386 : // -----------------------------------------------------------------------
387 0 : void GslMultifitWeightedLinear::setFlags(Vector<Bool> &flags, Bool goodIsTrue)
388 : {
389 : // Dim check
390 0 : ThrowIf (flags.size() != ndata_p,"Size of flags does not match model");
391 :
392 : // Fill in
393 0 : if (goodIsTrue)
394 : {
395 0 : for (size_t idx=0;idx<ndata_p;idx++)
396 : {
397 0 : weights_p(idx) = flags(idx);
398 : }
399 : }
400 : else
401 : {
402 0 : for (size_t idx=0;idx<ndata_p;idx++)
403 : {
404 0 : weights_p(idx) = !flags(idx);
405 : }
406 : }
407 0 : }
408 :
409 : // -----------------------------------------------------------------------
410 : //
411 : // -----------------------------------------------------------------------
412 68102 : void GslMultifitWeightedLinear::setWeightsAndFlags(Vector<Float> &weights, Vector<Bool> &flags, Bool goodIsTrue)
413 : {
414 : // Dim check
415 68102 : ThrowIf (weights.size() != ndata_p,"Size of weights does not match model");
416 68102 : ThrowIf (flags.size() != ndata_p,"Size of flags does not match model");
417 :
418 : // Fill in
419 68102 : if (goodIsTrue)
420 : {
421 19574290 : for (size_t idx=0;idx<ndata_p;idx++)
422 : {
423 19506188 : weights_p(idx) = weights(idx)*flags(idx);
424 : }
425 : }
426 : else
427 : {
428 0 : for (size_t idx=0;idx<ndata_p;idx++)
429 : {
430 0 : weights_p(idx) = weights(idx)*!flags(idx);
431 : }
432 : }
433 68102 : }
434 :
435 : // -----------------------------------------------------------------------
436 : //
437 : // Perform least-squares fits to a general linear weighted model, y = X c where
438 : // y is a vector of n observations, with weights w, X is an n by p matrix of
439 : // predictor variables, and the elements of the vector c are the p unknown best-fit
440 : // parameters which are to be estimated.
441 : //
442 : // NOTE: gsl_multifit_linear expects that the model matrix is organized as follows
443 : // X = [ 1 , x_1 , x_1^2 , ..., x_1^order;
444 : // 1 , x_2 , x_2^2 , ..., x_2^order;
445 : // 1 , x_3 , x_3^2 , ..., x_3^order;
446 : // ... , ... , ... , ..., ...
447 : // 1 , x_n , x_n^2 , ..., x_n^order]
448 : //
449 : // NOTE: More than one data series can use the same weights / workspace
450 : // Therefore input data is a matrix where each row represents a data series
451 : //
452 : // -----------------------------------------------------------------------
453 0 : double GslMultifitWeightedLinear::calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff)
454 : {
455 : // Wrap data vector in a gsl_vector
456 : gsl_vector data_gsl;
457 0 : GslVectorWrap(data,data_gsl);
458 :
459 : // Perform coeff calculation
460 0 : errno_p = gsl_multifit_wlinear (&gsl_model_p, &gls_weights_p, &data_gsl,
461 : coeff, gsl_covariance_p, &chisq_p, gsl_workspace_p);
462 :
463 0 : return chisq_p;
464 : }
465 :
466 : //////////////////////////////////////////////////////////////////////////
467 : // IterativelyReweightedLeastSquares class
468 : //////////////////////////////////////////////////////////////////////////
469 :
470 : // -----------------------------------------------------------------------
471 : //
472 : // -----------------------------------------------------------------------
473 3354 : IterativelyReweightedLeastSquares::IterativelyReweightedLeastSquares() :
474 3354 : GslMultifitWeightedLinear()
475 : {
476 3354 : nIter_p = 1;
477 3354 : }
478 :
479 : // -----------------------------------------------------------------------
480 : //
481 : // -----------------------------------------------------------------------
482 0 : IterativelyReweightedLeastSquares::IterativelyReweightedLeastSquares(GslLinearModelBase<Double> &model,size_t nIter) :
483 0 : GslMultifitWeightedLinear(model)
484 : {
485 0 : nIter_p = nIter;
486 0 : }
487 :
488 : // -----------------------------------------------------------------------
489 : //
490 : // -----------------------------------------------------------------------
491 3354 : IterativelyReweightedLeastSquares::~IterativelyReweightedLeastSquares()
492 : {
493 :
494 3354 : }
495 :
496 : // -----------------------------------------------------------------------
497 : //
498 : // -----------------------------------------------------------------------
499 136204 : double IterativelyReweightedLeastSquares::calcFitCoeffCore(Vector<Double> data, gsl_vector* coeff)
500 : {
501 : // Wrap data vector in a gsl_vector
502 : gsl_vector data_gsl;
503 136204 : GslVectorWrap(data,data_gsl);
504 :
505 136204 : if (nIter_p == 1)
506 : {
507 136204 : errno_p = gsl_multifit_wlinear (&gsl_model_p, &gls_weights_p, &data_gsl,
508 : coeff, gsl_covariance_p, &chisq_p, gsl_workspace_p);
509 : }
510 : else
511 : {
512 : // Create a vector to store iterative weights and wrap it in a gsl_vector
513 0 : Vector<Double> reweights(ndata_p);
514 0 : reweights = weights_p; // Deep copy
515 : gsl_vector reweights_gsl;
516 0 : GslVectorWrap(reweights,reweights_gsl);
517 :
518 : // Create vectors to store model
519 0 : Vector<Double> model(ndata_p);
520 :
521 : // Iterative process
522 0 : for (size_t iter=0; iter<nIter_p; iter++)
523 : {
524 : // Calculate coefficients for this iteration
525 0 : errno_p = gsl_multifit_wlinear (&gsl_model_p, &reweights_gsl, &data_gsl,
526 : coeff, gsl_covariance_p, &chisq_p, gsl_workspace_p);
527 :
528 0 : if (iter<nIter_p-1)
529 : {
530 : // Calculate output std
531 0 : calcFitModelCore(model,coeff);
532 :
533 : // Update weights
534 0 : updateWeights(data,model,reweights);
535 : }
536 : }
537 0 : }
538 :
539 136204 : return chisq_p;
540 : }
541 :
542 : // -----------------------------------------------------------------------
543 : //
544 : // -----------------------------------------------------------------------
545 0 : void IterativelyReweightedLeastSquares::updateWeights(Vector<Double> &data, Vector<Double> &model, Vector<Double> &weights)
546 : {
547 : double currentResidual, maxResidual;
548 :
549 : // Find max residual
550 0 : maxResidual = 0;
551 0 : for (size_t idx=0; idx<ndata_p; idx++)
552 : {
553 0 : currentResidual = 0;
554 0 : if (weights(idx) > 0)
555 : {
556 0 : currentResidual = abs(data(idx)-model(idx));
557 0 : if (currentResidual > maxResidual) maxResidual = currentResidual;
558 : }
559 0 : weights(idx) = currentResidual;
560 : }
561 :
562 : // Normalize
563 0 : for (size_t idx=0; idx<ndata_p; idx++)
564 : {
565 0 : if (weights(idx) > 0) weights(idx) = (maxResidual - weights(idx))/maxResidual;
566 : }
567 :
568 0 : return;
569 : }
570 :
571 : } //# NAMESPACE DENOISING - END
572 :
573 : } //# NAMESPACE CASA - END
574 :
|
