Line data Source code
1 : //# UVContSubTVI.h: This file contains the implementation of the UVContSubTVI 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/UVContSubTVI.h>
24 : #include <casacore/ms/MeasurementSets/MSFieldColumns.h>
25 :
26 : // OpenMP
27 : #ifdef _OPENMP
28 : #include <omp.h>
29 : #endif
30 :
31 : using namespace casacore;
32 :
33 : namespace casa { //# NAMESPACE CASA - BEGIN
34 :
35 : namespace vi { //# NAMESPACE VI - BEGIN
36 :
37 0 : FitSpec:: FitSpec(Vector<bool> mask, unsigned int order){
38 0 : lineFreeChannelMask = mask;
39 0 : fitOrder = order;
40 0 : };
41 :
42 : //////////////////////////////////////////////////////////////////////////
43 : // UVContSubTVI class
44 : //////////////////////////////////////////////////////////////////////////
45 :
46 : // -----------------------------------------------------------------------
47 : //
48 : // -----------------------------------------------------------------------
49 54 : UVContSubTVI::UVContSubTVI( ViImplementation2 * inputVii,
50 54 : const Record &configuration):
51 54 : FreqAxisTVI (inputVii)
52 : {
53 54 : fitOrder_p = 0;
54 54 : want_cont_p = False;
55 54 : withDenoisingLib_p = true;
56 54 : nThreads_p = 1;
57 54 : niter_p = 1;
58 :
59 54 : inputFrequencyMap_p.clear();
60 :
61 : // Parse and check configuration parameters
62 : // Note: if a constructor finishes by throwing an exception, the memory
63 : // associated with the object itself is cleaned up — there is no memory leak.
64 54 : const auto fitspec = parseConfiguration(configuration);
65 :
66 44 : initialize(fitspec);
67 104 : }
68 :
69 88 : UVContSubTVI::~UVContSubTVI()
70 : {
71 44 : inputFrequencyMap_p.clear();
72 88 : }
73 :
74 :
75 : /**
76 : * Check if this is a valid field index
77 : *
78 : * @param str string from a (multi)field spec string
79 : *
80 : * @return whether this is a valid field index
81 : */
82 74 : bool isFieldIndex(const std::string &str) {
83 74 : return !str.empty() && std::all_of(str.begin(), str.end(), ::isdigit);
84 : }
85 :
86 : /**
87 : * For a field string (example: '3', '0,1,4') get the integer indices
88 : * specified in the string.
89 : * Tokenizes the string by ',' separator and trims spaces.
90 : *
91 : * @param a string with field indices as used in uvcontsub/fitspec parameter
92 : *
93 : * @return vector of field indices
94 : */
95 38 : std::vector<int> stringToFieldIdxs(const std::string &str)
96 : {
97 38 : stringstream stream(str);
98 38 : vector<int> result;
99 77 : while(stream.good()) {
100 41 : string token;
101 41 : getline(stream, token, ',' );
102 41 : token.erase(std::remove(token.begin(), token.end(), ' '), token.end());
103 41 : if (!isFieldIndex(token)) {
104 2 : throw AipsError("Invalid field index: " + token + " (in " + str + ").");
105 :
106 : }
107 : try {
108 39 : auto idx = std::stoi(token);
109 39 : result.push_back(idx);
110 0 : } catch (const std::exception &exc) {
111 0 : throw AipsError("Error trying to parse field ID: " + token + ", " + exc.what());
112 0 : }
113 41 : }
114 72 : return result;
115 40 : }
116 :
117 54 : InFitSpecMap UVContSubTVI::parseFitSpec(const Record &configuration) const
118 : {
119 54 : InFitSpecMap fitspec;
120 :
121 54 : const auto exists = configuration.fieldNumber ("fitspec");
122 54 : if (exists < 0) {
123 0 : return fitspec;
124 : }
125 :
126 : try {
127 : // fitspec is a simple string (spw:chan MSSelection syntax)
128 54 : String specStr;
129 75 : configuration.get(exists, specStr);
130 :
131 99 : std::vector <InFitSpec> specs = { InFitSpec(specStr, fitOrder_p) };
132 : // -1 (all fields): global fit spw:chan string and global fitOrder
133 33 : fitspec.insert({-1, specs});
134 75 : } catch(AipsError &exc) {
135 : // alternatively, fitspec must be a record with field/spw specs
136 21 : fitspec = parseDictFitSpec(configuration);
137 21 : }
138 44 : return fitspec;
139 10 : }
140 :
141 : /**
142 : * For an spw string given in a fitspec per-field subdictionary
143 : * (example: '0', '1,2') get the integer indices specified in the
144 : * string. Tokenizes the string by ',' separator and trims
145 : * spaces. Similar to stringToFieldIDxs but with spw-specific checks.
146 : *
147 : * @param a string with SPW indices as used in uvcontsub/fitspec parameter
148 : * @param allowSemicolon whether to accept the ':' spw:chan separator
149 : *
150 : * @return vector of SPW indices
151 : */
152 75 : std::vector<unsigned int> stringToSPWIdxs(const std::string &spwStr,
153 : bool allowSemicolon=false) {
154 75 : if (0 == spwStr.length()) {
155 0 : throw AipsError("Unexpected empty SPW IDs string");
156 : }
157 75 : if (!allowSemicolon && spwStr.find(':') != std::string::npos) {
158 0 : throw AipsError("Unexpected spw:chan separator character ':' found in SPW: " +
159 0 : spwStr + ". Not allowed in fitspec dictionaries");
160 : }
161 :
162 75 : vector<unsigned int> result;
163 75 : stringstream stream(spwStr);
164 150 : while(stream.good()) {
165 76 : string token;
166 76 : getline(stream, token, ',' );
167 76 : token.erase(std::remove(token.begin(), token.end(), ' '), token.end());
168 76 : if (!allowSemicolon && !isFieldIndex(token)) {
169 1 : throw AipsError("Invalid SPW index: " + token + " (in " + spwStr + ")");
170 : }
171 : try {
172 75 : auto idx = std::stoi(token);
173 75 : result.push_back(idx);
174 0 : } catch (const std::exception &exc) {
175 0 : throw AipsError("Error trying to parse SPW: " + token + ", " + exc.what());
176 0 : }
177 76 : }
178 148 : return result;
179 76 : }
180 :
181 : /**
182 : * Get max valid FIELD IDs for this iterator. MSv2 uses the FIELD
183 : * table row index as FIELD ID.
184 : */
185 21 : rownr_t UVContSubTVI::getMaxMSFieldID() const
186 : {
187 21 : const auto &fieldsTable = getVii()->fieldSubtablecols();
188 21 : return fieldsTable.nrow() - 1;
189 : }
190 :
191 28 : void UVContSubTVI::insertToFieldSpecMap(const std::vector<int> &fieldIdxs,
192 : const InFitSpec &spec,
193 : InFitSpecMap &fitspec) const
194 : {
195 57 : for (const auto fid: fieldIdxs) {
196 29 : const auto fieldFound = fitspec.find(fid);
197 29 : if (fieldFound == fitspec.end()) {
198 69 : std::vector<InFitSpec> firstSpw = { spec };
199 23 : fitspec[fid] = firstSpw;
200 23 : } else {
201 6 : fitspec[fid].emplace_back(spec);
202 : }
203 :
204 : }
205 28 : }
206 :
207 : /**
208 : * Takes one per-field subdict, parses it and populates a map field IDs->fit specs.
209 : *
210 : * @param fieldRec subdict from the configuration dict/record passed from the task/tool
211 : * interface, with uvcontsub task parameters
212 : * @param fieldIdx IDs of the fields this spec applies to
213 : * @param fitspec spec map to populate for this(ese) field(s)
214 : */
215 24 : void UVContSubTVI::parseFieldSubDict(const Record &fieldRec,
216 : const std::vector<int> &fieldIdxs,
217 : InFitSpecMap &fitspec) const
218 : {
219 24 : std::set<unsigned int> spwsSeen;
220 52 : for (unsigned int spw=0; spw < fieldRec.nfields(); ++spw) {
221 32 : const auto spwRec = fieldRec.subRecord(RecordFieldId(spw));
222 64 : std::string spwStr = fieldRec.name(RecordFieldId(spw));
223 32 : const auto spwIdxs = stringToSPWIdxs(spwStr, false);
224 61 : for (const auto sid : spwIdxs) {
225 31 : if (spwsSeen.insert(sid).second == false) {
226 3 : throw AipsError("Error in fitspec. SPW " + to_string(sid) + " is given "
227 2 : "multiple times. Found for SPW subdict '" + spwStr +
228 2 : "' and other subdicts(s).");
229 : }
230 : }
231 : try {
232 31 : const std::string chanStr = spwRec.asString(RecordFieldId("chan"));
233 29 : if (!chanStr.empty()) {
234 26 : spwStr = "";
235 52 : for (const auto sid : spwIdxs) {
236 52 : const auto spw_chan = to_string(sid) + ":" + chanStr;
237 26 : if (spwStr.empty()) {
238 26 : spwStr = spw_chan;
239 : } else {
240 0 : spwStr += "," + spw_chan;
241 : }
242 26 : }
243 : }
244 30 : } catch (const AipsError &exc) {
245 1 : throw AipsError("Error trying to get chan value in subdict for spw "
246 2 : + spwStr + ", " + exc.what());
247 1 : }
248 :
249 29 : int polOrder = fitOrder_p;
250 : try {
251 30 : polOrder = spwRec.asInt(RecordFieldId("fitorder"));
252 1 : } catch (const AipsError &exc) {
253 1 : throw AipsError("Error trying to get fitorder value in subdict for spw "
254 2 : + spwStr + ", " + exc.what());
255 1 : }
256 28 : if (polOrder < 0) {
257 0 : throw AipsError("Fit order cannot be negative. Value given: " +
258 0 : to_string(polOrder));
259 : }
260 :
261 28 : const auto spec = std::make_pair(spwStr, polOrder);
262 28 : insertToFieldSpecMap(fieldIdxs, spec, fitspec);
263 39 : }
264 24 : }
265 :
266 : /**
267 : * Takes the input dict of fieldID: SPW: {chan/fitorder}, parses it
268 : * and populates a map that is returned.
269 : *
270 : * @param configuration dictionary/record passed from the task/tool interface, with
271 : * uvcontsub task parameters
272 : *
273 : * @return InFitSpecMap map populated from an input dict
274 : */
275 21 : InFitSpecMap UVContSubTVI::parseDictFitSpec(const Record &configuration) const
276 : {
277 21 : const Record rawFieldFitspec = configuration.asRecord("fitspec");
278 21 : if (rawFieldFitspec.empty()) {
279 0 : throw AipsError("The dictionary 'fitspec' of fit specifications is empty!");
280 : }
281 :
282 21 : const auto maxMSField = getMaxMSFieldID();
283 21 : InFitSpecMap fitspec;
284 21 : std::set<unsigned int> fieldsSeen;
285 : // iterate through fields given in fitspec
286 49 : for (unsigned int rid=0; rid < rawFieldFitspec.nfields(); ++rid) {
287 76 : const std::string fieldsStr = rawFieldFitspec.name(RecordFieldId(rid));
288 38 : const auto fieldIdxs = stringToFieldIdxs(fieldsStr);
289 72 : for (const auto fid: fieldIdxs) {
290 : // check that the field is in the MS
291 39 : if (fid < 0 || static_cast<unsigned int>(fid) > maxMSField) {
292 4 : throw AipsError("Wrong field ID given: " +
293 6 : std::to_string(fid) +
294 : ". This MeasurementSet has field IDs between"
295 8 : " 0 and " + std::to_string(maxMSField));
296 : }
297 37 : if (fieldsSeen.insert(fid).second == false) {
298 3 : throw AipsError("Error in fitspec. Field " + to_string(fid) + " is given "
299 2 : "multiple times. Found for field subdict '" + fieldsStr +
300 2 : "' and other subdict(s).");
301 : }
302 : }
303 :
304 33 : std::string noneStr;
305 33 : bool subRecordIsString = true;
306 : try {
307 57 : noneStr = rawFieldFitspec.asString(RecordFieldId(rid));
308 24 : } catch (const AipsError &exc) {
309 24 : subRecordIsString = false;
310 24 : }
311 33 : if (subRecordIsString){
312 9 : if (noneStr != "NONE") {
313 1 : throw AipsError("Wrong value found in fitspec, field subdict for field: '"
314 2 : + fieldsStr + "', value: '" +
315 2 : noneStr + "'. Only 'NONE' is accepted (= do not fit for " +
316 2 : "any SPW).");
317 : }
318 : } else {
319 24 : const Record fieldRec = rawFieldFitspec.subRecord(RecordFieldId(rid));
320 24 : parseFieldSubDict(fieldRec, fieldIdxs, fitspec);
321 24 : }
322 51 : }
323 :
324 22 : return fitspec;
325 41 : }
326 :
327 118 : std::string fieldTextFromId(int idx)
328 : {
329 118 : std::string name;
330 118 : if (idx >= 0) {
331 52 : name = std::to_string(idx);
332 : } else {
333 66 : name = "all fields";
334 : }
335 118 : return name;
336 0 : }
337 :
338 44 : void UVContSubTVI::printInputFitSpec(const InFitSpecMap &fitspec) const
339 : {
340 44 : if (fitspec.size() > 0) {
341 88 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
342 : << "Fitspec (order and spw:line-free channel) specification is: "
343 88 : << LogIO::POST;
344 :
345 97 : for (const auto &elem: fitspec) {
346 53 : const auto &specs = elem.second;
347 112 : for (const auto &oneSpw : specs) {
348 118 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
349 59 : << " field: " << fieldTextFromId(elem.first)
350 118 : << ". " << oneSpw.second << ", '" << oneSpw.first << "'"
351 118 : << LogIO::POST;
352 : }
353 : }
354 : } else {
355 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
356 0 : << "Line-free channel selection: not given" << LogIO::POST;
357 : }
358 44 : }
359 :
360 : // -----------------------------------------------------------------------
361 : //
362 : // -----------------------------------------------------------------------
363 54 : InFitSpecMap UVContSubTVI::parseConfiguration(const Record &configuration)
364 : {
365 54 : auto exists = configuration.fieldNumber ("fitorder");
366 54 : if (exists >= 0)
367 : {
368 54 : configuration.get (exists, fitOrder_p);
369 :
370 54 : if (fitOrder_p > 0)
371 : {
372 18 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
373 18 : << "Global/default fit order is " << fitOrder_p << LogIO::POST;
374 : }
375 : }
376 :
377 54 : exists = configuration.fieldNumber ("want_cont");
378 54 : if (exists >= 0)
379 : {
380 54 : configuration.get (exists, want_cont_p);
381 :
382 54 : if (want_cont_p)
383 : {
384 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
385 : << "Producing continuum estimate instead of continuum subtracted data "
386 0 : << LogIO::POST;
387 : }
388 : }
389 :
390 54 : const auto fitspec = parseFitSpec(configuration);
391 44 : printInputFitSpec(fitspec);
392 :
393 44 : exists = configuration.fieldNumber ("writemodel");
394 44 : if (exists >= 0)
395 : {
396 44 : configuration.get(exists, precalcModel_p);
397 44 : if (precalcModel_p)
398 : {
399 6 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
400 : << "Producing continuum estimate in the MODEL_DATA column"
401 6 : << LogIO::POST;
402 : }
403 : }
404 :
405 44 : exists = configuration.fieldNumber ("denoising_lib");
406 44 : if (exists >= 0)
407 : {
408 44 : configuration.get (exists, withDenoisingLib_p);
409 :
410 44 : if (withDenoisingLib_p)
411 : {
412 86 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
413 : << "Using GSL based multiparameter regression with linear "
414 86 : << "least-squares fitting" << LogIO::POST;
415 : }
416 : }
417 :
418 44 : exists = configuration.fieldNumber ("nthreads");
419 44 : if (exists >= 0)
420 : {
421 0 : configuration.get (exists, nThreads_p);
422 :
423 0 : if (nThreads_p > 1)
424 : {
425 : #ifdef _OPENMP
426 0 : if (omp_get_max_threads() < (int)nThreads_p)
427 : {
428 0 : logger_p << LogIO::WARN << LogOrigin("UVContSubTVI", __FUNCTION__)
429 : << "Requested " << nThreads_p
430 : << " OMP threads but maximum possible is " << omp_get_max_threads()<< endl
431 : << "Setting number of OMP threads to " << omp_get_max_threads() << endl
432 : << "Check OMP_NUM_THREADS environmental variable and number of cores in your system"
433 0 : << LogIO::POST;
434 0 : nThreads_p = omp_get_max_threads();
435 : }
436 : else
437 : {
438 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
439 0 : << "Numer of OMP threads set to " << nThreads_p << LogIO::POST;
440 : }
441 : #else
442 : logger_p << LogIO::WARN << LogOrigin("UVContSubTVI", __FUNCTION__)
443 : << "Requested " << nThreads_p << " threads but OMP is not available in your system"
444 : << LogIO::POST;
445 : nThreads_p = 1;
446 : #endif
447 : }
448 : }
449 :
450 44 : exists = configuration.fieldNumber ("niter");
451 44 : if (exists >= 0)
452 : {
453 44 : configuration.get (exists, niter_p);
454 :
455 44 : if (niter_p > 1)
456 : {
457 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
458 0 : << "Number of iterations for re-weighted linear fit: " << niter_p << LogIO::POST;
459 : }
460 : }
461 :
462 : // If the user gives SPWs in fitspec that are not in this list, give a warning
463 44 : exists = configuration.fieldNumber ("allowed_spws");
464 44 : if (exists >= 0)
465 : {
466 44 : String allowed;
467 44 : configuration.get (exists, allowed);
468 44 : allowedSpws_p = allowed;
469 44 : }
470 :
471 44 : return fitspec;
472 0 : }
473 :
474 : /**
475 : * Check consistency of spw related parameters (spw and fitspec). Prints a warning
476 : * if there is an spw selection and fitspec uses SPWs not included in the selection.
477 : * This condition was found to be confusing for users trying to set parameters as
478 : * in the old uvcontsub with combine='spw'.
479 : *
480 : * @param msVii MS object attached to the Vi (remember - not the right
481 : * thing to use in a TVI)
482 : * @param spwChan MSSelection with the spw/chan (spwExpr) selection set
483 : */
484 37 : void UVContSubTVI::spwInputChecks(const MeasurementSet *msVii,
485 : MSSelection &spwChan) const {
486 : // note: spwChan should be const (if getSpwList() were const)
487 37 : if (!allowedSpws_p.empty()) {
488 5 : const auto usedSpws = spwChan.getSpwList(msVii).tovector();
489 5 : MSSelection checker;
490 5 : checker.setSpwExpr(allowedSpws_p);
491 5 : const auto allowed = checker.getSpwList(msVii).tovector();
492 10 : for (const auto spw : usedSpws ) {
493 5 : if (std::find(allowed.begin(), allowed.end(), spw) == allowed.end()) {
494 8 : logger_p << LogIO::WARN << LogOrigin("UVContSubTVI", __FUNCTION__)
495 : << "SPW " << spw << " is used in fitspec but is not included in "
496 4 : << "the SPW selection ('" << allowedSpws_p << "'), please "
497 8 : << "double-check the spw and fitspec parameters." << LogIO::POST;
498 : }
499 : }
500 5 : }
501 37 : }
502 :
503 : /**
504 : * Produces a map spw -> vector of selected channels, to be used as an
505 : * intermediate step to produce a channel mask for the fitting
506 : * routines.
507 : *
508 : * @param spwChanStr An SPW:chan selection string (MSSelection syntax)
509 : *
510 : * @return Map from spw (int ID) -> vector of channel indices
511 : */
512 37 : unordered_map<int, vector<int>> UVContSubTVI::makeLineFreeChannelSelMap(std::string
513 : spwChanStr) const
514 : {
515 37 : MSSelection mssel;
516 37 : mssel.setSpwExpr(spwChanStr);
517 : // This will access the MS directly: far from ideal (CAS-11638) but no easy solution
518 37 : const auto msVii = &(inputVii_p->ms());
519 :
520 37 : spwInputChecks(msVii, mssel);
521 :
522 : // Create line-free channel map based on MSSelection channel ranges
523 37 : const auto spwchan = mssel.getChanList(msVii);
524 37 : const auto nSelections = spwchan.shape()[0];
525 37 : unordered_map<int,vector<int>> lineFreeChannelSelMap;
526 109 : for (uInt selection_i=0; selection_i<nSelections; ++selection_i)
527 : {
528 72 : auto spw = spwchan(selection_i, 0);
529 72 : if (lineFreeChannelSelMap.find(spw) == lineFreeChannelSelMap.end())
530 : {
531 43 : lineFreeChannelSelMap[spw].clear(); // Accessing the vector creates it
532 : }
533 :
534 72 : const auto channelStart = spwchan(selection_i, 1);
535 72 : const auto channelStop = spwchan(selection_i, 2);
536 72 : const auto channelStep = spwchan(selection_i, 3);
537 28128 : for (auto inpChan=channelStart; inpChan<=channelStop; inpChan += channelStep)
538 : {
539 28056 : lineFreeChannelSelMap[spw].push_back(inpChan);
540 : }
541 : }
542 :
543 74 : return lineFreeChannelSelMap;
544 37 : }
545 :
546 : // -----------------------------------------------------------------------
547 : //
548 : // -----------------------------------------------------------------------
549 37 : void UVContSubTVI::populatePerFieldSpec(int fieldID,
550 : const std::vector<InFitSpec> &fitSpecs)
551 : {
552 37 : std::unordered_map<int, int> perSpwFitOrder;
553 37 : std::string fullSpwStr = "";
554 80 : for (const auto &spec : fitSpecs) {
555 43 : if (fullSpwStr == "NONE") {
556 0 : break;
557 43 : } else if (fullSpwStr.length() == 0) {
558 37 : fullSpwStr = spec.first;
559 : } else {
560 : // MSSelection syntax spw:chan, SPWs separated by commas
561 6 : fullSpwStr += "," + spec.first;
562 : }
563 43 : auto spws = stringToSPWIdxs(spec.first, true);
564 86 : for (const auto sid : spws) {
565 43 : perSpwFitOrder[sid] = spec.second;
566 : }
567 43 : }
568 :
569 37 : auto lineFreeChannelSelMap = makeLineFreeChannelSelMap(fullSpwStr);
570 :
571 : // Create line-free channel mask, spw->(channel_mask, fit_order)
572 37 : unordered_map<int, FitSpec> lineFreeChannelMaskMap; // rename: fitSpecMap
573 410 : for (auto const spwInp: spwInpChanIdxMap_p)
574 : {
575 373 : const auto spw = spwInp.first;
576 373 : if (lineFreeChannelMaskMap.find(spw) == lineFreeChannelMaskMap.end())
577 : {
578 373 : if (lineFreeChannelSelMap.size() > 0 && 0 == lineFreeChannelSelMap[spw].size()) {
579 : // Some SPWs have been selected, but this SPW doesn't have any selection
580 : // => a 0-elements mask says there is nothing to mask or fit here, or
581 : // otherwise we'd make the effort to prepare and call the fit routine for an
582 : // "all masked/True" channel mask which will not even start minimizing
583 330 : lineFreeChannelMaskMap[spw].lineFreeChannelMask = Vector<bool>();
584 : } else {
585 43 : auto &mask = lineFreeChannelMaskMap[spw].lineFreeChannelMask;
586 43 : mask = Vector<bool>(spwInp.second.size(), true);
587 28099 : for (uInt selChanIdx=0; selChanIdx<lineFreeChannelSelMap[spw].size();
588 : ++selChanIdx)
589 : {
590 28056 : const auto selChan = lineFreeChannelSelMap[spw][selChanIdx];
591 28056 : mask(selChan) = False;
592 : }
593 : }
594 :
595 373 : unsigned int fitOrder = fitOrder_p;
596 373 : const auto find = perSpwFitOrder.find(spw);
597 373 : if (find != perSpwFitOrder.end()) {
598 43 : fitOrder = find->second;
599 : }
600 373 : lineFreeChannelMaskMap[spw].fitOrder = fitOrder;
601 : }
602 373 : }
603 :
604 : // Poppulate per field map: field -> spw -> fit spec
605 : // emplace struct (linefreechannelmaskmap, + fitorder)
606 37 : perFieldSpecMap_p.emplace(fieldID, lineFreeChannelMaskMap);
607 37 : }
608 :
609 : /**
610 : * Takes the map from fitspec and produces the per field map to be
611 : * used in the iteration, with channel masks and fitorder per field
612 : * per SPW.
613 : */
614 44 : void UVContSubTVI::fitSpecToPerFieldMap(const InFitSpecMap &fitspec)
615 : {
616 : // Process line-free channel specifications
617 97 : for (const auto item: fitspec) {
618 53 : unordered_map<int, FitSpec> fieldSpecMap;
619 : // Parse line-free channel selection using MSSelection syntax
620 53 : const auto fieldID = item.first;
621 53 : const auto &specs = item.second;
622 53 : bool noneFound = false;
623 112 : for (const auto spwSpec : specs) {
624 59 : const auto spwStr = spwSpec.first;
625 59 : const auto order = spwSpec.second;
626 59 : const auto fieldName = fieldTextFromId(fieldID);
627 118 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
628 : << "Parsing fit spw:chan string and order for field: " << fieldName
629 118 : << ", spw/chan: '" << spwStr << "', order: " << order << LogIO::POST;
630 59 : if (spwStr == "NONE") {
631 0 : noneFound = true;
632 0 : if (specs.size() != 1) {
633 0 : throw AipsError("For field '" + fieldName + "' A \"NONE\" fit specification "
634 0 : "has been given but additional fit specs have been given");
635 : }
636 : }
637 59 : }
638 :
639 53 : if (noneFound) {
640 : // Not inserting any entries in perFieldSpecMap_p[fieldID]
641 : // implies no transformation for that field (inputVis -> outputVis)
642 0 : continue;
643 : }
644 :
645 53 : if (1 == specs.size() && specs[0].first.empty()) {
646 : // -1 is the "all-fields-included" field pseudo-index
647 : // empty fieldSpecMap string -> all SPW, channels, leave all SPW masks unset
648 16 : if (-1 == fieldID) {
649 16 : perFieldSpecMap_p.emplace(fieldID, fieldSpecMap);
650 : }
651 16 : continue;
652 : }
653 :
654 37 : populatePerFieldSpec(fieldID, specs);
655 69 : }
656 44 : }
657 :
658 : // -----------------------------------------------------------------------
659 : //
660 : // -----------------------------------------------------------------------
661 44 : void UVContSubTVI::initialize(const InFitSpecMap &fitspec)
662 : {
663 : // Populate nchan input-output maps
664 490 : for (auto spwInp: spwInpChanIdxMap_p)
665 : {
666 446 : auto spw = spwInp.first;
667 446 : spwOutChanNumMap_p[spw] = spwInp.second.size();
668 446 : }
669 :
670 44 : fitSpecToPerFieldMap(fitspec);
671 44 : }
672 :
673 : // -----------------------------------------------------------------------
674 : //
675 : // -----------------------------------------------------------------------
676 0 : void UVContSubTVI::floatData (Cube<Float> & vis) const
677 : {
678 : // Get input VisBuffer
679 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
680 :
681 : // Transform data
682 0 : transformDataCube(vb->visCubeFloat(),vb->weightSpectrum(),vis);
683 :
684 0 : return;
685 : }
686 :
687 : /**
688 : * To save the original_data - continumm_subtracted_data as model data
689 : * for visibilityModel.
690 : *
691 : * @param origVis visibilities before continuum subtraction
692 : * @param contSubvis visibilities after being cont-subtracted
693 : */
694 2875 : void UVContSubTVI::savePrecalcModel(const Cube<Complex> & origVis,
695 : const Cube<Complex> & contSubVis) const
696 : {
697 : // save it for visibilityModel
698 2875 : if (precalcModel_p) {
699 : // Or otherwise the next assignment would fail a conform check
700 53 : if (precalcModelVis_p.shape() != contSubVis.shape()) {
701 3 : precalcModelVis_p.resize(contSubVis.shape());
702 : }
703 :
704 53 : precalcModelVis_p = origVis - contSubVis;
705 : }
706 2875 : }
707 :
708 : // -----------------------------------------------------------------------
709 : //
710 : // -----------------------------------------------------------------------
711 2855 : void UVContSubTVI::visibilityObserved (Cube<Complex> & vis) const
712 : {
713 : // Get input VisBuffer
714 2855 : VisBuffer2 *vb = getVii()->getVisBuffer();
715 :
716 : // Get weightSpectrum from sigmaSpectrum
717 2855 : Cube<Float> weightSpFromSigmaSp;
718 2855 : weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),False);
719 2855 : weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
720 2855 : arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
721 :
722 : // Transform data
723 2855 : transformDataCube(vb->visCube(),weightSpFromSigmaSp,vis);
724 :
725 2855 : savePrecalcModel(vb->visCube(), vis);
726 2855 : }
727 :
728 : // -----------------------------------------------------------------------
729 : //
730 : // -----------------------------------------------------------------------
731 20 : void UVContSubTVI::visibilityCorrected (Cube<Complex> & vis) const
732 : {
733 : // Get input VisBuffer
734 20 : VisBuffer2 *vb = getVii()->getVisBuffer();
735 :
736 : // Transform data
737 20 : transformDataCube(vb->visCubeCorrected(),vb->weightSpectrum(),vis);
738 :
739 20 : savePrecalcModel(vb->visCubeCorrected(), vis);
740 20 : }
741 :
742 : // -----------------------------------------------------------------------
743 : //
744 : // -----------------------------------------------------------------------
745 53 : void UVContSubTVI::visibilityModel (Cube<Complex> & vis) const
746 : {
747 : // Get input VisBuffer
748 53 : VisBuffer2 *vb = getVii()->getVisBuffer();
749 :
750 : // from visiblityObserved we have calculated the polynomial subtraction
751 53 : if (precalcModel_p && precalcModelVis_p.shape() == vb->visCubeModel().shape()) {
752 53 : vis = precalcModelVis_p;
753 53 : return;
754 : }
755 :
756 : // Transform data
757 0 : transformDataCube(vb->visCubeModel(),vb->weightSpectrum(),vis);
758 :
759 0 : return;
760 : }
761 :
762 : // -----------------------------------------------------------------------
763 : //
764 : // -----------------------------------------------------------------------
765 44 : void UVContSubTVI::result(Record &res) const
766 : {
767 44 : auto acc = result_p.getAccumulatedResult();
768 44 : res.defineRecord("uvcontsub", acc);
769 44 : }
770 :
771 : // -----------------------------------------------------------------------
772 : //
773 : // -----------------------------------------------------------------------
774 2875 : template<class T> void UVContSubTVI::transformDataCube( const Cube<T> &inputVis,
775 : const Cube<Float> &inputWeight,
776 : Cube<T> &outputVis) const
777 : {
778 : // Get input VisBuffer
779 2875 : VisBuffer2 *vb = getVii()->getVisBuffer();
780 :
781 2875 : auto fieldID = vb->fieldId()[0];
782 : // First check the "all fields" (no per-individual field list in fitspec)
783 2875 : auto fieldMapIt = perFieldSpecMap_p.find(-1);
784 2875 : if (fieldMapIt == perFieldSpecMap_p.end()) {
785 : // otherwise, check for individual fields (coming from a fitspec list)
786 1205 : fieldMapIt = perFieldSpecMap_p.find(fieldID);
787 1205 : if (fieldMapIt == perFieldSpecMap_p.end()) {
788 : // This was a "NONE" = no-subtraction for this field ==> no-op
789 120 : outputVis = inputVis;
790 1201 : return;
791 : }
792 : }
793 :
794 2755 : Int spwId = vb->spectralWindows()(0);
795 :
796 : // Get input line-free channel mask and fitorder
797 2755 : Vector<bool> *lineFreeChannelMask = nullptr;
798 2755 : auto fitOrder = fitOrder_p;
799 2755 : auto spwMap = fieldMapIt->second;
800 2755 : auto maskLookup = spwMap.find(spwId);
801 2755 : if (maskLookup != spwMap.end())
802 : {
803 1834 : if (maskLookup->second.lineFreeChannelMask.nelements() == 0) {
804 : // This was a non-selected SPW in a non-empty SPW selection string ==> no-op
805 1081 : outputVis = inputVis;
806 1081 : return;
807 : } else {
808 753 : lineFreeChannelMask = &(maskLookup->second.lineFreeChannelMask);
809 753 : fitOrder = maskLookup->second.fitOrder;
810 : }
811 : }
812 :
813 : // Get polynomial model for this SPW (depends on number of channels and gridding)
814 1674 : const auto freqIter = inputFrequencyMap_p.find(spwId);
815 1674 : const Vector<Double> &inputFrequencies = vb->getFrequencies(0);
816 : // Could perhaps have a per field-spw pair map to avoid re-allocations - But can be big
817 : // (n_fields X n_spws X n_chans) for many fields, many SPWs MSs
818 : inputFrequencyMap_p[spwId].
819 1674 : reset(new denoising::GslPolynomialModel<double>(inputFrequencies, fitOrder));
820 :
821 : // Reshape output data before passing it to the DataCubeHolder
822 1674 : outputVis.resize(getVisBuffer()->getShape(),False);
823 :
824 : // Get input flag Cube
825 1674 : const Cube<bool> &flagCube = vb->flagCube();
826 :
827 : // Transform data
828 1674 : if (nThreads_p > 1)
829 : {
830 : #ifdef _OPENMP
831 :
832 0 : uInt nCorrs = vb->getShape()(0);
833 0 : if (nCorrs < nThreads_p)
834 : {
835 0 : omp_set_num_threads(nCorrs);
836 : }
837 : else
838 : {
839 0 : omp_set_num_threads(nThreads_p);
840 : }
841 :
842 0 : #pragma omp parallel for
843 : for (uInt corrIdx=0; corrIdx < nCorrs; corrIdx++)
844 : {
845 : transformDataCore(inputFrequencyMap_p[spwId].get(),lineFreeChannelMask,
846 : inputVis,flagCube,inputWeight,outputVis,corrIdx);
847 : }
848 :
849 0 : omp_set_num_threads(nThreads_p);
850 : #endif
851 : }
852 : else
853 : {
854 1674 : auto scanID = vb->scan()[0];
855 1674 : uInt nCorrs = vb->getShape()(0);
856 5096 : for (uInt corrIdx=0; corrIdx < nCorrs; corrIdx++)
857 : {
858 : Complex chisq =
859 3422 : transformDataCore(inputFrequencyMap_p[spwId].get(), lineFreeChannelMask,
860 : inputVis, flagCube, inputWeight, outputVis, corrIdx);
861 3422 : result_p.addOneFit(fieldID, scanID, spwId, (int)corrIdx, chisq);
862 : }
863 : }
864 2755 : }
865 :
866 : // -----------------------------------------------------------------------
867 : //
868 : // -----------------------------------------------------------------------
869 3422 : template<class T> Complex UVContSubTVI::transformDataCore(denoising::GslPolynomialModel<Double>* model,
870 : Vector<bool> *lineFreeChannelMask,
871 : const Cube<T> &inputVis,
872 : const Cube<bool> &inputFlags,
873 : const Cube<Float> &inputWeight,
874 : Cube<T> &outputVis,
875 : Int parallelCorrAxis) const
876 : {
877 : // Gather input data
878 3422 : DataCubeMap inputData;
879 3422 : DataCubeHolder<T> inputVisCubeHolder(inputVis);
880 3422 : DataCubeHolder<bool> inputFlagCubeHolder(inputFlags);
881 3422 : DataCubeHolder<Float> inputWeightsCubeHolder(inputWeight);
882 3422 : inputData.add(MS::DATA,inputVisCubeHolder);
883 3422 : inputData.add(MS::FLAG,inputFlagCubeHolder);
884 3422 : inputData.add(MS::WEIGHT_SPECTRUM,inputWeightsCubeHolder);
885 :
886 : // Gather output data
887 3422 : DataCubeMap outputData;
888 3422 : DataCubeHolder<T> outputVisCubeHolder(outputVis);
889 3422 : outputData.add(MS::DATA,outputVisCubeHolder);
890 :
891 3422 : Complex chisq;
892 3422 : if (want_cont_p)
893 : {
894 0 : if (withDenoisingLib_p)
895 : {
896 0 : UVContEstimationDenoisingKernel<T> kernel(model,niter_p,lineFreeChannelMask);
897 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
898 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
899 0 : chisq = kernel.getChiSquared();
900 0 : }
901 : else
902 : {
903 0 : UVContEstimationKernel<T> kernel(model,lineFreeChannelMask);
904 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
905 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
906 0 : chisq = kernel.getChiSquared();
907 0 : }
908 : }
909 : else
910 : {
911 3422 : if (withDenoisingLib_p)
912 : {
913 3354 : UVContSubtractionDenoisingKernel<T> kernel(model,niter_p,lineFreeChannelMask);
914 3354 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
915 3354 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
916 3354 : chisq = kernel.getChiSquared();
917 3354 : }
918 : else
919 : {
920 68 : UVContSubtractionKernel<T> kernel(model,lineFreeChannelMask);
921 68 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
922 68 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
923 68 : chisq = kernel.getChiSquared();
924 68 : }
925 : }
926 :
927 3422 : return chisq;
928 3422 : }
929 :
930 : //////////////////////////////////////////////////////////////////////////
931 : // UVContSubTVIFactory class
932 : //////////////////////////////////////////////////////////////////////////
933 :
934 : // -----------------------------------------------------------------------
935 : //
936 : // -----------------------------------------------------------------------
937 0 : UVContSubTVIFactory::UVContSubTVIFactory ( Record &configuration,
938 0 : ViImplementation2 *inputVii)
939 : {
940 0 : inputVii_p = inputVii;
941 0 : configuration_p = configuration;
942 0 : }
943 :
944 : // -----------------------------------------------------------------------
945 : //
946 : // -----------------------------------------------------------------------
947 0 : vi::ViImplementation2 * UVContSubTVIFactory::createVi(VisibilityIterator2 *) const
948 : {
949 0 : return new UVContSubTVI(inputVii_p,configuration_p);
950 : }
951 :
952 : // -----------------------------------------------------------------------
953 : //
954 : // -----------------------------------------------------------------------
955 0 : vi::ViImplementation2 * UVContSubTVIFactory::createVi() const
956 : {
957 0 : return new UVContSubTVI(inputVii_p,configuration_p);
958 : }
959 :
960 : //////////////////////////////////////////////////////////////////////////
961 : // UVContSubTVILayerFactory class
962 : //////////////////////////////////////////////////////////////////////////
963 :
964 54 : UVContSubTVILayerFactory::UVContSubTVILayerFactory(Record &configuration) :
965 54 : ViiLayerFactory()
966 : {
967 54 : configuration_p = configuration;
968 54 : }
969 :
970 : ViImplementation2*
971 54 : UVContSubTVILayerFactory::createInstance(ViImplementation2* vii0) const
972 : {
973 : // Make the UVContSubTVI, using supplied ViImplementation2, and return it
974 54 : ViImplementation2 *vii = new UVContSubTVI(vii0, configuration_p);
975 44 : return vii;
976 : }
977 :
978 : //////////////////////////////////////////////////////////////////////////
979 : // UVContSubTransformEngine class
980 : //////////////////////////////////////////////////////////////////////////
981 :
982 : // -----------------------------------------------------------------------
983 : //
984 : // -----------------------------------------------------------------------
985 3422 : template<class T> UVContSubTransformEngine<T>::UVContSubTransformEngine( UVContSubKernel<T> *kernel,
986 : DataCubeMap *inputData,
987 : DataCubeMap *outputData):
988 3422 : FreqAxisTransformEngine2<T>(inputData,outputData)
989 : {
990 3422 : uvContSubKernel_p = kernel;
991 3422 : }
992 :
993 : // -----------------------------------------------------------------------
994 : //
995 : // -----------------------------------------------------------------------
996 69280 : template<class T> void UVContSubTransformEngine<T>::transform( )
997 : {
998 69280 : uvContSubKernel_p->setDebug(debug_p);
999 69280 : uvContSubKernel_p->kernel(inputData_p,outputData_p);
1000 69280 : }
1001 :
1002 : //////////////////////////////////////////////////////////////////////////
1003 : // UVContSubKernel class
1004 : //////////////////////////////////////////////////////////////////////////
1005 :
1006 : // -----------------------------------------------------------------------
1007 : //
1008 : // -----------------------------------------------------------------------
1009 3422 : template<class T> UVContSubKernel<T>::UVContSubKernel( denoising::GslPolynomialModel<Double> *model,
1010 3422 : Vector<bool> *lineFreeChannelMask)
1011 : {
1012 3422 : model_p = model;
1013 3422 : fitOrder_p = model_p->ncomponents()-1;
1014 3422 : freqPows_p.reference(model_p->getModelMatrix());
1015 3422 : frequencies_p.reference(model_p->getLinearComponentFloat());
1016 :
1017 3422 : lineFreeChannelMask_p = lineFreeChannelMask ? lineFreeChannelMask : nullptr;
1018 3422 : debug_p = False;
1019 3422 : }
1020 :
1021 : // -----------------------------------------------------------------------
1022 : //
1023 : // -----------------------------------------------------------------------
1024 69280 : template<class T> void UVContSubKernel<T>::kernel(DataCubeMap *inputData,
1025 : DataCubeMap *outputData)
1026 : {
1027 : // Get input/output data
1028 69280 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1029 69280 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1030 :
1031 : // Apply line-free channel mask
1032 69280 : Vector<bool> &inputFlags = inputData->getVector<bool>(MS::FLAG);
1033 69280 : if (lineFreeChannelMask_p != nullptr)
1034 40186 : inputFlags |= *lineFreeChannelMask_p;
1035 :
1036 : // Calculate number of valid data points and adapt fit
1037 69280 : size_t validPoints = nfalse(inputFlags);
1038 69280 : if (validPoints > 0)
1039 : {
1040 68782 : bool restoreDefaultPoly = False;
1041 68782 : uInt tmpFitOrder = fitOrder_p;
1042 :
1043 : // Reduce fit order to match number of valid points
1044 68782 : if (validPoints <= fitOrder_p)
1045 : {
1046 1800 : changeFitOrder(validPoints-1);
1047 1800 : restoreDefaultPoly = true;
1048 : }
1049 :
1050 : // Get weights
1051 68782 : Vector<Float> &inputWeight = inputData->getVector<Float>(MS::WEIGHT_SPECTRUM);
1052 :
1053 : // Convert flags to mask
1054 68782 : Vector<bool> mask = !inputFlags;
1055 :
1056 : // Calculate and subtract continuum
1057 68782 : chisq_p = kernelCore(inputVector,mask,inputWeight,outputVector);
1058 :
1059 : // Go back to default fit order to match number of valid points
1060 68782 : if (restoreDefaultPoly)
1061 : {
1062 1800 : changeFitOrder(tmpFitOrder);
1063 : }
1064 68782 : }
1065 : else
1066 : {
1067 498 : defaultKernel(inputVector,outputVector);
1068 : }
1069 69280 : }
1070 :
1071 :
1072 : //////////////////////////////////////////////////////////////////////////
1073 : // UVContSubtractionKernel class
1074 : //////////////////////////////////////////////////////////////////////////
1075 :
1076 : // -----------------------------------------------------------------------
1077 : //
1078 : // -----------------------------------------------------------------------
1079 68 : template<class T> UVContSubtractionKernel<T>::UVContSubtractionKernel( denoising::GslPolynomialModel<Double>* model,
1080 : Vector<bool> *lineFreeChannelMask):
1081 68 : UVContSubKernel<T>(model,lineFreeChannelMask)
1082 : {
1083 68 : changeFitOrder(fitOrder_p);
1084 68 : }
1085 :
1086 : // -----------------------------------------------------------------------
1087 : //
1088 : // -----------------------------------------------------------------------
1089 68 : template<class T> void UVContSubtractionKernel<T>::changeFitOrder(size_t order)
1090 : {
1091 68 : fitOrder_p = order;
1092 68 : Polynomial<AutoDiff<Float> > poly(order);
1093 68 : fitter_p.setFunction(poly); // poly It is cloned
1094 68 : fitter_p.asWeight(true);
1095 :
1096 136 : return;
1097 68 : }
1098 :
1099 :
1100 : // -----------------------------------------------------------------------
1101 : //
1102 : // -----------------------------------------------------------------------
1103 0 : template<class T> void UVContSubtractionKernel<T>::defaultKernel( Vector<Complex> &inputVector,
1104 : Vector<Complex> &outputVector)
1105 : {
1106 0 : outputVector = inputVector;
1107 0 : return;
1108 : }
1109 :
1110 : // -----------------------------------------------------------------------
1111 : //
1112 : // -----------------------------------------------------------------------
1113 0 : template<class T> void UVContSubtractionKernel<T>::defaultKernel( Vector<Float> &inputVector,
1114 : Vector<Float> &outputVector)
1115 : {
1116 0 : outputVector = inputVector;
1117 0 : return;
1118 : }
1119 :
1120 : // -----------------------------------------------------------------------
1121 : //
1122 : // -----------------------------------------------------------------------
1123 680 : template<class T> Complex UVContSubtractionKernel<T>::kernelCore( Vector<Complex> &inputVector,
1124 : Vector<bool> &inputFlags,
1125 : Vector<Float> &inputWeights,
1126 : Vector<Complex> &outputVector)
1127 : {
1128 : // Fit for imaginary and real components separately
1129 680 : Vector<Float> realCoeff;
1130 680 : Vector<Float> imagCoeff;
1131 680 : realCoeff = fitter_p.fit(frequencies_p, real(inputVector), inputWeights, &inputFlags);
1132 680 : double realChisq = fitter_p.chiSquare();
1133 680 : imagCoeff = fitter_p.fit(frequencies_p, imag(inputVector), inputWeights, &inputFlags);
1134 680 : double imagChisq = fitter_p.chiSquare();
1135 :
1136 : // Fill output data
1137 680 : outputVector = inputVector;
1138 680 : outputVector -= Complex(realCoeff(0),imagCoeff(0));
1139 680 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1140 : {
1141 0 : Complex coeff(realCoeff(order_idx),imagCoeff(order_idx));
1142 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1143 : {
1144 0 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff;
1145 : }
1146 : }
1147 :
1148 680 : if (debug_p)
1149 : {
1150 0 : LogIO logger;
1151 0 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1152 0 : logger << "realCoeff =" << realCoeff << LogIO::POST;
1153 0 : logger << "imagCoeff =" << imagCoeff << LogIO::POST;
1154 0 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1155 0 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1156 0 : logger << "inputVector =" << inputVector << LogIO::POST;
1157 0 : logger << "outputVector =" << outputVector << LogIO::POST;
1158 0 : }
1159 :
1160 1360 : return Complex(realChisq, imagChisq);
1161 680 : }
1162 :
1163 : // -----------------------------------------------------------------------
1164 : //
1165 : // -----------------------------------------------------------------------
1166 0 : template<class T> Complex UVContSubtractionKernel<T>::kernelCore( Vector<Float> &inputVector,
1167 : Vector<bool> &inputFlags,
1168 : Vector<Float> &inputWeights,
1169 : Vector<Float> &outputVector)
1170 : {
1171 : // Fit model
1172 0 : Vector<Float> coeff;
1173 0 : coeff = fitter_p.fit(frequencies_p, inputVector, inputWeights, &inputFlags);
1174 0 : const double chisq = fitter_p.chiSquare();
1175 :
1176 : // Fill output data
1177 0 : outputVector = inputVector;
1178 0 : outputVector -= coeff(0);
1179 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1180 : {
1181 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1182 : {
1183 0 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1184 : }
1185 : }
1186 :
1187 0 : if (debug_p)
1188 : {
1189 0 : LogIO logger;
1190 0 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1191 0 : logger << "coeff =" << coeff << LogIO::POST;
1192 0 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1193 0 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1194 0 : logger << "inputVector =" << inputVector << LogIO::POST;
1195 0 : logger << "outputVector =" << outputVector << LogIO::POST;
1196 0 : }
1197 :
1198 0 : return Complex(chisq, chisq);
1199 0 : }
1200 :
1201 :
1202 : //////////////////////////////////////////////////////////////////////////
1203 : // UVContEstimationKernel class
1204 : //////////////////////////////////////////////////////////////////////////
1205 :
1206 : // -----------------------------------------------------------------------
1207 : //
1208 : // -----------------------------------------------------------------------
1209 0 : template<class T> UVContEstimationKernel<T>::UVContEstimationKernel( denoising::GslPolynomialModel<Double>* model,
1210 : Vector<bool> *lineFreeChannelMask):
1211 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1212 : {
1213 0 : changeFitOrder(fitOrder_p);
1214 0 : }
1215 :
1216 : // -----------------------------------------------------------------------
1217 : //
1218 : // -----------------------------------------------------------------------
1219 0 : template<class T> void UVContEstimationKernel<T>::changeFitOrder(size_t order)
1220 : {
1221 0 : fitOrder_p = order;
1222 0 : Polynomial<AutoDiff<Float> > poly(order);
1223 0 : fitter_p.setFunction(poly); // poly It is cloned
1224 0 : fitter_p.asWeight(true);
1225 :
1226 0 : return;
1227 0 : }
1228 :
1229 : // -----------------------------------------------------------------------
1230 : //
1231 : // -----------------------------------------------------------------------
1232 0 : template<class T> void UVContEstimationKernel<T>::defaultKernel(Vector<Complex> &,
1233 : Vector<Complex> &outputVector)
1234 : {
1235 0 : outputVector = 0;
1236 0 : return;
1237 : }
1238 :
1239 : // -----------------------------------------------------------------------
1240 : //
1241 : // -----------------------------------------------------------------------
1242 0 : template<class T> void UVContEstimationKernel<T>::defaultKernel(Vector<Float> &,
1243 : Vector<Float> &outputVector)
1244 : {
1245 0 : outputVector = 0;
1246 0 : return;
1247 : }
1248 :
1249 : // -----------------------------------------------------------------------
1250 : //
1251 : // -----------------------------------------------------------------------
1252 0 : template<class T> Complex UVContEstimationKernel<T>::kernelCore(Vector<Complex> &inputVector,
1253 : Vector<bool> &inputFlags,
1254 : Vector<Float> &inputWeights,
1255 : Vector<Complex> &outputVector
1256 : )
1257 : {
1258 : // Fit for imaginary and real components separately
1259 0 : Vector<Float> realCoeff;
1260 0 : Vector<Float> imagCoeff;
1261 0 : realCoeff = fitter_p.fit(frequencies_p, real(inputVector), inputWeights, &inputFlags);
1262 0 : double realChisq = fitter_p.chiSquare();
1263 :
1264 0 : imagCoeff = fitter_p.fit(frequencies_p, imag(inputVector), inputWeights, &inputFlags);
1265 0 : double imagChisq = fitter_p.chiSquare();
1266 :
1267 : // Fill output data
1268 0 : outputVector = Complex(realCoeff(0),imagCoeff(0));
1269 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1270 : {
1271 0 : Complex coeff(realCoeff(order_idx),imagCoeff(order_idx));
1272 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1273 : {
1274 0 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff;
1275 : }
1276 : }
1277 :
1278 0 : return Complex(realChisq, imagChisq);
1279 0 : }
1280 :
1281 : // -----------------------------------------------------------------------
1282 : //
1283 : // -----------------------------------------------------------------------
1284 0 : template<class T> Complex UVContEstimationKernel<T>::kernelCore( Vector<Float> &inputVector,
1285 : Vector<bool> &inputFlags,
1286 : Vector<Float> &inputWeights,
1287 : Vector<Float> &outputVector)
1288 : {
1289 : // Fit model
1290 0 : Vector<Float> coeff;
1291 0 : coeff = fitter_p.fit(frequencies_p, inputVector, inputWeights, &inputFlags);
1292 0 : double chisq = fitter_p.chiSquare();
1293 :
1294 : // Fill output data
1295 0 : outputVector = coeff(0);
1296 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1297 : {
1298 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1299 : {
1300 0 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1301 : }
1302 : }
1303 :
1304 0 : return Complex(chisq, chisq);
1305 0 : }
1306 :
1307 : //////////////////////////////////////////////////////////////////////////
1308 : // UVContSubtractionDenoisingKernel class
1309 : //////////////////////////////////////////////////////////////////////////
1310 :
1311 : // -----------------------------------------------------------------------
1312 : //
1313 : // -----------------------------------------------------------------------
1314 3354 : template<class T> UVContSubtractionDenoisingKernel<T>::UVContSubtractionDenoisingKernel(denoising::GslPolynomialModel<Double>* model,
1315 : size_t nIter,
1316 : Vector<bool> *lineFreeChannelMask):
1317 3354 : UVContSubKernel<T>(model,lineFreeChannelMask)
1318 : {
1319 3354 : fitter_p.resetModel(*model);
1320 3354 : fitter_p.setNIter(nIter);
1321 3354 : }
1322 :
1323 : // -----------------------------------------------------------------------
1324 : //
1325 : // -----------------------------------------------------------------------
1326 3600 : template<class T> void UVContSubtractionDenoisingKernel<T>::changeFitOrder(size_t order)
1327 : {
1328 3600 : fitOrder_p = order;
1329 3600 : fitter_p.resetNComponents(order+1);
1330 3600 : return;
1331 : }
1332 :
1333 :
1334 : // -----------------------------------------------------------------------
1335 : //
1336 : // -----------------------------------------------------------------------
1337 498 : template<class T> void UVContSubtractionDenoisingKernel<T>::defaultKernel( Vector<T> &inputVector,
1338 : Vector<T> &outputVector)
1339 : {
1340 498 : outputVector = inputVector;
1341 498 : return;
1342 : }
1343 :
1344 : // -----------------------------------------------------------------------
1345 : //
1346 : // -----------------------------------------------------------------------
1347 68102 : template<class T> Complex UVContSubtractionDenoisingKernel<T>::kernelCore( Vector<T> &inputVector,
1348 : Vector<bool> &inputFlags,
1349 : Vector<Float> &inputWeights,
1350 : Vector<T> &outputVector)
1351 68102 : { fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1352 68102 : vector<T> coeff;
1353 68102 : Complex chisq;
1354 68102 : tie(coeff, chisq) = fitter_p.calcFitCoeff(inputVector);
1355 :
1356 68102 : Vector<T> model(outputVector.size());
1357 68102 : fitter_p.calcFitModel(model);
1358 :
1359 68102 : outputVector = inputVector;
1360 68102 : outputVector -= model;
1361 :
1362 : /*
1363 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1364 : Vector<T> coeff = fitter_p.calcFitCoeff(inputVector);
1365 :
1366 : // Fill output data
1367 : outputVector = inputVector;
1368 : outputVector -= coeff(0);
1369 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1370 : {
1371 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1372 : {
1373 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1374 : }
1375 : }
1376 : */
1377 :
1378 : /*
1379 : if (debug_p)
1380 : {
1381 : LogIO logger;
1382 : logger << "freqPows_p =" << freqPows_p << LogIO::POST;
1383 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1384 : logger << "coeff =" << coeff << LogIO::POST;
1385 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1386 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1387 : logger << "inputVector =" << inputVector << LogIO::POST;
1388 : logger << "outputVector =" << outputVector << LogIO::POST;
1389 : }
1390 : */
1391 :
1392 68102 : return chisq;
1393 68102 : }
1394 :
1395 : //////////////////////////////////////////////////////////////////////////
1396 : // UVContEstimationDenoisingKernel class
1397 : //////////////////////////////////////////////////////////////////////////
1398 :
1399 : // -----------------------------------------------------------------------
1400 : //
1401 : // -----------------------------------------------------------------------
1402 0 : template<class T> UVContEstimationDenoisingKernel<T>::UVContEstimationDenoisingKernel( denoising::GslPolynomialModel<Double>* model,
1403 : size_t nIter,
1404 : Vector<bool> *lineFreeChannelMask):
1405 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1406 : {
1407 0 : fitter_p.resetModel(*model);
1408 0 : fitter_p.setNIter(nIter);
1409 0 : }
1410 :
1411 : // -----------------------------------------------------------------------
1412 : //
1413 : // -----------------------------------------------------------------------
1414 0 : template<class T> void UVContEstimationDenoisingKernel<T>::changeFitOrder(size_t order)
1415 : {
1416 0 : fitOrder_p = order;
1417 0 : fitter_p.resetNComponents(order+1);
1418 0 : return;
1419 : }
1420 :
1421 : // -----------------------------------------------------------------------
1422 : //
1423 : // -----------------------------------------------------------------------
1424 0 : template<class T> void UVContEstimationDenoisingKernel<T>::defaultKernel( Vector<T> &,
1425 : Vector<T> &outputVector)
1426 : {
1427 0 : outputVector = 0;
1428 0 : return;
1429 : }
1430 :
1431 : // -----------------------------------------------------------------------
1432 : //
1433 : // -----------------------------------------------------------------------
1434 0 : template<class T> Complex UVContEstimationDenoisingKernel<T>::kernelCore( Vector<T> &inputVector,
1435 : Vector<bool> &inputFlags,
1436 : Vector<Float> &inputWeights,
1437 : Vector<T> &outputVector)
1438 : {
1439 0 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1440 0 : std::vector<T> coeff;
1441 0 : Complex chisq;
1442 0 : tie(coeff, chisq) = fitter_p.calcFitCoeff(inputVector);
1443 0 : fitter_p.calcFitModel(outputVector);
1444 :
1445 : /*
1446 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1447 : Vector<T> coeff = fitter_p.calcFitCoeff(inputVector);
1448 :
1449 : // Fill output data
1450 : outputVector = coeff(0);
1451 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1452 : {
1453 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1454 : {
1455 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1456 : }
1457 : }
1458 : */
1459 :
1460 0 : return chisq;
1461 0 : }
1462 :
1463 : } //# NAMESPACE VI - END
1464 :
1465 : } //# NAMESPACE CASA - END
1466 :
1467 :
|
