Line data Source code
1 : /*
2 : * SDDoubleCircleGainCal.cpp
3 : *
4 : * Created on: May 31, 2016
5 : * Author: nakazato
6 : */
7 :
8 : //#include "SDDoubleCircleGainCal.h"
9 : #include "SDDoubleCircleGainCalImpl.h"
10 :
11 : #include <casacore/casa/Arrays/Vector.h>
12 : #include <casacore/casa/Arrays/ArrayMath.h>
13 : #include <casacore/casa/IO/ArrayIO.h>
14 : #include <casacore/scimath/Functionals/Interpolate1D.h>
15 : #include <casacore/scimath/Functionals/ScalarSampledFunctional.h>
16 : #include <casacore/casa/BasicSL/Constants.h>
17 : #include <casacore/casa/Logging/LogIO.h>
18 : #include <casacore/casa/Quanta/Quantum.h>
19 :
20 : #include <casacore/casa/Utilities/Sort.h>
21 :
22 : #include <cassert>
23 : #include <cmath>
24 :
25 : using namespace casacore;
26 :
27 : #define LOG logger_ << LogOrigin("SDDoubleCircleGainCal", __FUNCTION__, WHERE)
28 : #define POSTLOG LogIO::POST
29 :
30 : namespace { // anonymous namespace START
31 : // primary beam size based on observing frequency and antenna diameter in radian
32 0 : inline Double getPrimaryBeamSize(Double const observing_frequency,
33 : Double const antenna_diameter) {
34 0 : Double beam_size = -1.0;
35 0 : constexpr Double kFactorALMA = 1.13; // measured factor for ALMA
36 0 : Double const speed_of_light = C::c;
37 : //Double const rad2arcsec = 180.0 / C::pi * 3600.0;
38 :
39 0 : if (observing_frequency <= 0.0 || antenna_diameter <= 0.0) {
40 0 : beam_size = 0.0;
41 : } else {
42 0 : beam_size = kFactorALMA * speed_of_light // * rad2arcsec
43 0 : / (antenna_diameter * observing_frequency);
44 : }
45 0 : return beam_size;
46 : }
47 :
48 : // default smoothing size based on the radius of central region
49 : // radius will be given in unit of radian but the following formula
50 : // seems to assume radius in unit of arcsec
51 0 : inline Int getDefaultSmoothingSize(Double const radius) {
52 0 : auto const radius_in_arcsec = Quantity(radius, "rad").getValue("arcsec");
53 0 : return 2 * static_cast<Int>(round(radius_in_arcsec / 1.5)) + 1;
54 : }
55 :
56 : // average
57 0 : inline Double average(size_t const index_from, size_t const index_to,
58 : Vector<Double> const &data) {
59 0 : Double sum = 0.0;
60 0 : assert(index_from <= index_to);
61 : //cout << "number of data to be averaged " << index_to - index_from << endl;
62 0 : for (size_t i = index_from; i < index_to; ++i) {
63 0 : sum += data[i];
64 : }
65 0 : return sum / static_cast<Double>(index_to - index_from);
66 : }
67 0 : inline Double average(size_t const index_from, size_t const index_to,
68 : Vector<Double> const &data, Vector<Bool> const &flag) {
69 0 : Double sum = 0.0;
70 0 : size_t count = 0;
71 0 : assert(index_from <= index_to);
72 : //cout << "number of data to be averaged " << index_to - index_from << endl;
73 0 : for (size_t i = index_from; i < index_to; ++i) {
74 0 : if (flag[i] == false) {
75 0 : sum += data[i];
76 0 : count++;
77 : }
78 : }
79 :
80 0 : if (count == 0) {
81 0 : return 0.0;
82 : }
83 :
84 0 : return sum / static_cast<Double>(count);
85 : }
86 :
87 : // smoothing
88 0 : inline void smooth(Int const smooth_size, Vector<Double> const &data,
89 : Vector<Double> &smoothed_data) {
90 : // TODO replace with sakura function
91 0 : assert(data.nelements() == smoothed_data.nelements());
92 0 : size_t num_data = data.nelements();
93 0 : if (smooth_size < 2 || static_cast<size_t>(smooth_size) >= num_data) {
94 : //cout << "no smoothing" << endl;
95 0 : smoothed_data = data;
96 : } else {
97 0 : size_t left_edge = (smooth_size + 1) / 2 - 1;
98 0 : size_t right_edge = smooth_size / 2 + 1;
99 0 : for (size_t i = 0; i < left_edge; ++i) {
100 0 : size_t l = 0;
101 0 : size_t r = i + right_edge;
102 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
103 0 : smoothed_data[i] = average(l, r, data);
104 : }
105 0 : for (size_t i = left_edge; i < num_data - right_edge; ++i) {
106 0 : size_t l = i - left_edge;
107 0 : size_t r = i + right_edge;
108 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
109 0 : smoothed_data[i] = average(l, r, data);
110 : }
111 0 : for (size_t i = num_data - right_edge; i < num_data; ++i) {
112 0 : size_t l = i - left_edge;
113 0 : size_t r = num_data;
114 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
115 0 : smoothed_data[i] = average(l, r, data);
116 : }
117 : }
118 0 : }
119 :
120 0 : inline void smooth(Int const smooth_size, Vector<Double> const &data,
121 : Vector<Bool> const &flag, Vector<Double> &smoothed_data) {
122 : // TODO replace with sakura function
123 0 : assert(data.nelements() == smoothed_data.nelements());
124 0 : size_t num_data = data.nelements();
125 0 : if (smooth_size < 2 || static_cast<size_t>(smooth_size) >= num_data) {
126 : //cout << "no smoothing" << endl;
127 0 : smoothed_data = data;
128 : } else {
129 0 : size_t left_edge = (smooth_size + 1) / 2 - 1;
130 0 : size_t right_edge = smooth_size / 2 + 1;
131 0 : for (size_t i = 0; i < left_edge; ++i) {
132 0 : size_t l = 0;
133 0 : size_t r = i + right_edge;
134 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
135 0 : if (flag[i] == true) {
136 0 : smoothed_data[i] = data[i];
137 : } else {
138 0 : smoothed_data[i] = average(l, r, data, flag);
139 : }
140 : }
141 0 : for (size_t i = left_edge; i < num_data - right_edge; ++i) {
142 0 : size_t l = i - left_edge;
143 0 : size_t r = i + right_edge;
144 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
145 0 : if (flag[i] == true) {
146 0 : smoothed_data[i] = data[i];
147 : } else {
148 0 : smoothed_data[i] = average(l, r, data, flag);
149 : }
150 : }
151 0 : for (size_t i = num_data - right_edge; i < num_data; ++i) {
152 0 : size_t l = i - left_edge;
153 0 : size_t r = num_data;
154 : //cout << "i = " << i << ": l, r = " << l << "," << r << endl;
155 0 : if (flag[i] == true) {
156 0 : smoothed_data[i] = data[i];
157 : } else {
158 0 : smoothed_data[i] = average(l, r, data, flag);
159 : }
160 : }
161 : }
162 0 : }
163 :
164 : // interpolation
165 : // unused
166 : /*
167 : inline void interpolateLinear(Vector<Double> const &x0,
168 : Vector<Double> const &y0, Vector<Double> const &x1, Vector<Double> &y1) {
169 : // TODO replace with sakura function (need to add double precision version of interpolation function)
170 : Interpolate1D<Double, Double> interpolator(
171 : ScalarSampledFunctional<Double>(x0), ScalarSampledFunctional<Double>(y0),
172 : true, true);
173 : interpolator.setMethod(Interpolate1D<Double, Double>::linear);
174 : for (size_t i = 0; i < x1.nelements(); ++i) {
175 : y1[i] = interpolator(x1[i]);
176 : }
177 : }
178 : */
179 : // utility
180 0 : inline size_t toUnsigned(ssize_t const v) {
181 0 : assert(v >= 0);
182 0 : return static_cast<size_t>(v);
183 : }
184 : } // anonymous namespace END
185 :
186 : namespace casa { // namespace casa START
187 0 : SDDoubleCircleGainCalImpl::SDDoubleCircleGainCalImpl() :
188 0 : central_region_(-1.0), do_smooth_(false), smooth_size_(-1), observing_frequency_(
189 0 : 0.0), antenna_diameter_(0.0), logger_() {
190 0 : }
191 :
192 0 : SDDoubleCircleGainCalImpl::~SDDoubleCircleGainCalImpl() {
193 0 : }
194 :
195 0 : Double SDDoubleCircleGainCalImpl::getPrimaryBeamSize() const {
196 : // ::getPrimaryBeamSize returns the size in radian
197 0 : return ::getPrimaryBeamSize(observing_frequency_, antenna_diameter_);
198 : }
199 :
200 0 : Int SDDoubleCircleGainCalImpl::getDefaultSmoothingSize() const {
201 0 : Int default_size = -1;
202 0 : if (central_region_ > 0.0) {
203 0 : default_size = ::getDefaultSmoothingSize(central_region_);
204 : } else {
205 0 : default_size = ::getDefaultSmoothingSize(0.5 * getPrimaryBeamSize());
206 : }
207 0 : return default_size;
208 : }
209 :
210 0 : void SDDoubleCircleGainCalImpl::calibrate(Cube<Float> const &data,
211 : Vector<Double> const &time, Matrix<Double> const &direction,
212 : Vector<Double> &gain_time, Cube<Float> &gain) {
213 :
214 : // radius of the central region
215 0 : Double radius = getRadius();
216 :
217 0 : LOG << "radius = " << radius << POSTLOG;
218 :
219 : // select data within radius
220 0 : auto const data_shape = data.shape();
221 0 : size_t const num_pol = ::toUnsigned(data_shape[0]);
222 0 : size_t const num_chan = ::toUnsigned(data_shape[1]);
223 0 : assert(time.nelements() == uInt(data_shape[2]));
224 0 : assert(direction.shape()[1] == data_shape[2]);
225 0 : assert(direction.shape()[0] == 2);
226 0 : findDataWithinRadius(radius, time, data, direction, gain_time, gain);
227 0 : size_t num_gain = gain_time.nelements();
228 0 : LOG << "num_gain = " << num_gain << POSTLOG;
229 :
230 : //LOG << "indices within radius: " << within_radius << POSTLOG;
231 :
232 0 : if (num_gain < 100) {
233 0 : LOG << LogIO::WARN << "Probably not enough points for gain calibration: "
234 0 : << num_gain << endl << "Skipping..." << POSTLOG;
235 0 : gain_time.resize();
236 0 : gain.resize();
237 0 : return;
238 : }
239 :
240 : // for each spectral data
241 0 : Vector<Double> work_data(num_gain);
242 0 : Vector<Double> smoothed_data;
243 0 : if (do_smooth_) {
244 0 : smoothed_data.resize(num_gain);
245 : }
246 0 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
247 0 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
248 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
249 0 : work_data[idata] = gain(ipol, ichan, idata);
250 : }
251 :
252 : // LOG << "work_data[" << ipol << "," << ichan << "]=" << work_data
253 : // << POSTLOG;
254 :
255 : // smoothing if necessary
256 0 : if (do_smooth_) {
257 0 : Int smooth_size = getEffectiveSmoothingSize();
258 0 : if (smooth_size < 2 || static_cast<size_t>(smooth_size) >= num_gain) {
259 0 : LOG << LogIO::WARN
260 : << "data is not smoothed since smoothing size is invalid: "
261 : << smooth_size << " (number of data " << num_gain << ")"
262 0 : << POSTLOG;
263 : }
264 0 : smooth(smooth_size, work_data, smoothed_data);
265 : } else {
266 0 : LOG << "no smoothing" << POSTLOG;
267 0 : smoothed_data.reference(work_data);
268 : }
269 :
270 0 : LOG << LogIO::DEBUGGING << "smoothed_data[" << ipol << "," << ichan
271 0 : << "]=" << smoothed_data << POSTLOG;
272 :
273 0 : LOG << LogIO::DEBUGGING << "mean value = " << mean(smoothed_data)
274 0 : << POSTLOG;
275 :
276 : // derive gain factor: mean(smoothed_data) / smoothed_data
277 0 : work_data = mean(smoothed_data) / smoothed_data;
278 :
279 : // LOG << "gfactor[" << ipol << "," << ichan << "]=" << work_data
280 : // << POSTLOG;
281 :
282 : // conversion for G type calibration
283 0 : work_data = 1.0 / sqrt(work_data);
284 :
285 : // LOG << "fparam[" << ipol << "," << ichan << "]=" << work_data
286 : // << POSTLOG;
287 :
288 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
289 0 : gain(ipol, ichan, idata) = work_data[idata];
290 : }
291 : }
292 : }
293 0 : }
294 :
295 0 : void SDDoubleCircleGainCalImpl::calibrate(Cube<Float> const &data,
296 : Cube<Bool> const &flag, Vector<Double> const &time,
297 : Matrix<Double> const &direction, Vector<Double> &gain_time,
298 : Cube<Float> &gain, Cube<Bool> &gain_flag) {
299 :
300 : // radius of the central region
301 0 : Double radius = getRadius();
302 :
303 0 : LOG << "radius = " << radius << POSTLOG;
304 :
305 : // select data within radius
306 0 : auto const data_shape = data.shape();
307 0 : assert(time.nelements() == uInt(data_shape[2]));
308 0 : assert(direction.shape()[1] == data_shape[2]);
309 0 : assert(direction.shape()[0] == 2);
310 0 : findDataWithinRadius(radius, time, data, flag, direction, gain_time, gain,
311 : gain_flag);
312 0 : doCalibrate(gain_time, gain, gain_flag);
313 0 : }
314 :
315 0 : void SDDoubleCircleGainCalImpl::doCalibrate(Vector<Double> &gain_time,
316 : Cube<Float> &gain, Cube<Bool> &gain_flag) {
317 :
318 0 : size_t const num_gain = gain_time.nelements();
319 0 : LOG << "num_gain = " << num_gain << POSTLOG;
320 :
321 0 : size_t const num_pol = ::toUnsigned(gain.shape()[0]);
322 0 : size_t const num_chan = ::toUnsigned(gain.shape()[1]);
323 0 : assert(gain.shape()[2] == num_gain);
324 :
325 : //LOG << "indices within radius: " << within_radius << POSTLOG;
326 :
327 0 : if (num_gain < 100) {
328 0 : LOG << LogIO::WARN << "Probably not enough points for gain calibration: "
329 0 : << num_gain << endl << "Skipping..." << POSTLOG;
330 0 : gain_time.resize();
331 0 : gain.resize();
332 0 : gain_flag.resize();
333 0 : return;
334 : }
335 :
336 : // for each spectral data
337 0 : Vector<Double> work_data(num_gain);
338 0 : Vector<Bool> work_flag(num_gain);
339 0 : Vector<Double> smoothed_data(num_gain);
340 0 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
341 0 : for (size_t ichan = 0; ichan < num_chan; ++ichan) {
342 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
343 0 : work_data[idata] = gain(ipol, ichan, idata);
344 0 : work_flag[idata] = gain_flag(ipol, ichan, idata);
345 : }
346 :
347 : // LOG << "work_data[" << ipol << "," << ichan << "]=" << work_data
348 : // << POSTLOG;
349 :
350 : // smoothing if necessary
351 0 : if (do_smooth_) {
352 0 : Int smooth_size = getEffectiveSmoothingSize();
353 0 : if (smooth_size < 2 || static_cast<size_t>(smooth_size) >= num_gain) {
354 0 : LOG << LogIO::WARN
355 : << "data is not smoothed since smoothing size is invalid: "
356 : << smooth_size << " (number of data " << num_gain << ")"
357 0 : << POSTLOG;
358 : }
359 0 : smooth(smooth_size, work_data, work_flag, smoothed_data);
360 : } else {
361 0 : LOG << "no smoothing" << POSTLOG;
362 0 : smoothed_data = work_data;
363 : }
364 :
365 : // LOG << LogIO::DEBUGGING << "smoothed_data[" << ipol << "," << ichan
366 : // << "]=" << smoothed_data << POSTLOG;
367 :
368 : // derive gain factor: mean(smoothed_data) / smoothed_data
369 : //work_data = mean(smoothed_data) / smoothed_data;
370 0 : auto mean_value = ::average(0, num_gain, smoothed_data, work_flag);
371 :
372 : //work_data = mean_value / smoothed_data;
373 0 : if (mean_value < 0.0) {
374 0 : LOG << LogIO::WARN
375 : << "Negative reference value for gain calibration is found. "
376 0 : << "No valid calibration solution will be provided" << POSTLOG;
377 0 : work_data = 0.0;
378 0 : work_flag = True;
379 : } else {
380 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
381 0 : if (work_data[idata] != 0.0) {
382 0 : work_data[idata] = mean_value / smoothed_data[idata];
383 : }
384 : else {
385 0 : work_data[idata] = 0.0;
386 0 : work_flag[idata] = True;
387 : }
388 : }
389 :
390 0 : LOG << LogIO::DEBUGGING << "mean value = " << mean_value
391 0 : << " (simple mean " << mean(smoothed_data) << ")" << POSTLOG;
392 : }
393 :
394 : // LOG << LogIO::DEBUGGING << "gfactor[" << ipol << "," << ichan << "]=" << work_data
395 : // << POSTLOG;
396 :
397 : // conversion for G type calibration
398 : //work_data = 1.0 / sqrt(work_data);
399 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
400 0 : if (work_data[idata] > 0.0) {
401 0 : work_data[idata] = 1.0 / sqrt(work_data[idata]);
402 : } else {
403 0 : work_data[idata] = 0.0;
404 0 : work_flag[idata] = True;
405 : }
406 : }
407 :
408 : // LOG << LogIO::DEBUGGING << "fparam[" << ipol << "," << ichan << "]=" << work_data
409 : // << POSTLOG;
410 :
411 0 : for (size_t idata = 0; idata < num_gain; ++idata) {
412 0 : gain(ipol, ichan, idata) = work_data[idata];
413 0 : gain_flag(ipol, ichan, idata) = work_flag[idata];
414 : }
415 : }
416 : }
417 0 : }
418 :
419 0 : Double SDDoubleCircleGainCalImpl::getRadius() {
420 : // radius of the central region
421 0 : Double radius = central_region_;
422 0 : if (radius <= 0.0) {
423 : // use default value: primary beam size
424 0 : radius = 0.5 * getPrimaryBeamSize();
425 : }
426 0 : if (radius <= 0.0) {
427 0 : ostringstream ss;
428 0 : ss << "Size of central region is not properly set: " << radius;
429 0 : LOG << LogIO::SEVERE << ss.str() << POSTLOG;
430 0 : throw AipsError(ss.str());
431 0 : }
432 :
433 0 : return radius;
434 : }
435 :
436 0 : Int SDDoubleCircleGainCalImpl::getEffectiveSmoothingSize() {
437 0 : LOG << "do smoothing with size " << smooth_size_ << POSTLOG;
438 :
439 0 : Int smooth_size = smooth_size_;
440 0 : if (smooth_size_ < 0) {
441 0 : smooth_size = getDefaultSmoothingSize();
442 0 : LOG << "default smoothing size will be used: " << smooth_size
443 0 : << POSTLOG;
444 : }
445 :
446 0 : return smooth_size;
447 : }
448 :
449 0 : void SDDoubleCircleGainCalImpl::findDataWithinRadius(Double const radius,
450 : Vector<Double> const &time, Cube<Float> const &data,
451 : Matrix<Double> const &direction, Vector<Double> &gain_time,
452 : Cube<Float> &gain) {
453 0 : Cube<Bool> flag(data.shape(), False);
454 0 : Cube<Bool> gain_flag;
455 0 : findDataWithinRadius(radius, time, data, flag, direction, gain_time, gain,
456 : gain_flag);
457 0 : }
458 :
459 0 : void SDDoubleCircleGainCalImpl::findDataWithinRadius(Double const radius,
460 : Vector<Double> const &time, Cube<Float> const &data, Cube<Bool> const &flag,
461 : Matrix<Double> const &direction, Vector<Double> &gain_time,
462 : Cube<Float> &gain, Cube<Bool> &gain_flag) {
463 0 : size_t num_data = ::toUnsigned(direction.shape()[1]);
464 : // find data within radius
465 0 : Vector<size_t> data_index(num_data);
466 0 : size_t num_gain = 0;
467 0 : for (size_t i = 0; i < num_data; ++i) {
468 0 : Double x = direction(0, i);
469 0 : Double y = direction(1, i);
470 0 : Double r2 = x * x + y * y;
471 0 : if (r2 <= radius * radius) {
472 0 : data_index[num_gain] = i;
473 0 : num_gain++;
474 : }
475 : }
476 :
477 : // store data for calibration
478 0 : gain_time.resize(num_gain);
479 0 : IPosition gain_shape(data.shape());
480 0 : gain_shape[2] = num_gain;
481 0 : gain.resize(gain_shape);
482 0 : gain_flag.resize(gain_shape);
483 0 : for (size_t i = 0; i < num_gain; ++i) {
484 0 : size_t j = data_index[i];
485 0 : gain_time[i] = time[j];
486 0 : gain.xyPlane(i) = data.xyPlane(j);
487 0 : gain_flag.xyPlane(i) = flag.xyPlane(j);
488 : }
489 0 : }
490 :
491 0 : bool SDDoubleCircleGainCalImpl::findTimeRange(Vector<Double> const &time,
492 : Vector<Double> const &interval, Matrix<Double> const &direction,
493 : TimeRangeList &timerange) {
494 : // radius of the central region
495 0 : Double radius = getRadius();
496 :
497 : // find data within radius
498 0 : size_t num_all_data = time.nelements();
499 0 : Vector<size_t> data_index(num_all_data);
500 0 : size_t num_data = 0;
501 0 : for (size_t i = 0; i < num_all_data; ++i) {
502 0 : Double x = direction(0, i);
503 0 : Double y = direction(1, i);
504 0 : Double r2 = x * x + y * y;
505 0 : if (r2 <= radius * radius) {
506 0 : data_index[num_data] = i;
507 0 : num_data++;
508 : }
509 : }
510 :
511 : // list of time and interval
512 0 : Vector<Double> time_sel(num_data);
513 0 : Vector<Double> interval_sel(num_data);
514 0 : for (size_t i = 0; i < num_data; ++i) {
515 0 : time_sel[i] = time[data_index[i]];
516 0 : interval_sel[i] = interval[data_index[i]];
517 : }
518 :
519 : // sort
520 0 : Sort s(time_sel.data(), num_data);
521 0 : Vector<uInt> index_vector;
522 0 : uInt num_sorted = s.sort(index_vector, num_data);
523 0 : if (num_sorted != index_vector.nelements()) {
524 0 : return false;
525 : }
526 :
527 : // conversion from the list of timestamp to TimeRangeList
528 0 : timerange.clear();
529 0 : Double time_from = time_sel[index_vector[0]] - 0.5 * interval_sel[index_vector[0]];
530 0 : Double time_to = time_sel[index_vector[0]] + 0.5 * interval_sel[index_vector[0]];
531 0 : auto iter = index_vector.begin();
532 0 : for (++iter; iter != index_vector.end(); ++iter) {
533 0 : Double current_from = time_sel[*iter] - 0.5 * interval_sel[*iter];
534 0 : Double current_to = time_sel[*iter] + 0.5 * interval_sel[*iter];
535 0 : if (abs(current_from - time_to) < time_to * C::dbl_epsilon) {
536 0 : time_to = current_to;
537 : } else {
538 0 : TimeRange range(time_from, time_to);
539 0 : timerange.push_back(range);
540 0 : time_from = current_from;
541 0 : time_to = current_to;
542 : }
543 : }
544 0 : TimeRange range(time_from, time_to);
545 0 : timerange.push_back(range);
546 :
547 0 : return true;
548 0 : }
549 :
550 : //void SDDoubleCircleGainCalImpl::apply(Vector<Double> const &gain_time,
551 : // Cube<Float> const &gain, Vector<Double> const &time, Cube<Float> &data) {
552 : // // TODO implement
553 : // // not necessary to implement? reuse G type application?
554 : //}
555 : }// namespace casa END
|
