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 0 : UVContSubTVI::UVContSubTVI( ViImplementation2 * inputVii,
50 0 : const Record &configuration):
51 0 : FreqAxisTVI (inputVii)
52 : {
53 0 : fitOrder_p = 0;
54 0 : want_cont_p = False;
55 0 : withDenoisingLib_p = true;
56 0 : nThreads_p = 1;
57 0 : niter_p = 1;
58 :
59 0 : 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 0 : const auto fitspec = parseConfiguration(configuration);
65 :
66 0 : initialize(fitspec);
67 0 : }
68 :
69 0 : UVContSubTVI::~UVContSubTVI()
70 : {
71 0 : inputFrequencyMap_p.clear();
72 0 : }
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 0 : bool isFieldIndex(const std::string &str) {
83 0 : 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 0 : std::vector<int> stringToFieldIdxs(const std::string &str)
96 : {
97 0 : stringstream stream(str);
98 0 : vector<int> result;
99 0 : while(stream.good()) {
100 0 : string token;
101 0 : getline(stream, token, ',' );
102 0 : token.erase(std::remove(token.begin(), token.end(), ' '), token.end());
103 0 : if (!isFieldIndex(token)) {
104 0 : throw AipsError("Invalid field index: " + token + " (in " + str + ").");
105 :
106 : }
107 : try {
108 0 : auto idx = std::stoi(token);
109 0 : 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 0 : }
114 0 : return result;
115 0 : }
116 :
117 0 : InFitSpecMap UVContSubTVI::parseFitSpec(const Record &configuration) const
118 : {
119 0 : InFitSpecMap fitspec;
120 :
121 0 : const auto exists = configuration.fieldNumber ("fitspec");
122 0 : if (exists < 0) {
123 0 : return fitspec;
124 : }
125 :
126 : try {
127 : // fitspec is a simple string (spw:chan MSSelection syntax)
128 0 : String specStr;
129 0 : configuration.get(exists, specStr);
130 :
131 0 : std::vector <InFitSpec> specs = { InFitSpec(specStr, fitOrder_p) };
132 : // -1 (all fields): global fit spw:chan string and global fitOrder
133 0 : fitspec.insert({-1, specs});
134 0 : } catch(AipsError &exc) {
135 : // alternatively, fitspec must be a record with field/spw specs
136 0 : fitspec = parseDictFitSpec(configuration);
137 0 : }
138 0 : return fitspec;
139 0 : }
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 0 : std::vector<unsigned int> stringToSPWIdxs(const std::string &spwStr,
153 : bool allowSemicolon=false) {
154 0 : if (0 == spwStr.length()) {
155 0 : throw AipsError("Unexpected empty SPW IDs string");
156 : }
157 0 : 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 0 : vector<unsigned int> result;
163 0 : stringstream stream(spwStr);
164 0 : while(stream.good()) {
165 0 : string token;
166 0 : getline(stream, token, ',' );
167 0 : token.erase(std::remove(token.begin(), token.end(), ' '), token.end());
168 0 : if (!allowSemicolon && !isFieldIndex(token)) {
169 0 : throw AipsError("Invalid SPW index: " + token + " (in " + spwStr + ")");
170 : }
171 : try {
172 0 : auto idx = std::stoi(token);
173 0 : 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 0 : }
178 0 : return result;
179 0 : }
180 :
181 : /**
182 : * Get max valid FIELD IDs for this iterator. MSv2 uses the FIELD
183 : * table row index as FIELD ID.
184 : */
185 0 : rownr_t UVContSubTVI::getMaxMSFieldID() const
186 : {
187 0 : const auto &fieldsTable = getVii()->fieldSubtablecols();
188 0 : return fieldsTable.nrow() - 1;
189 : }
190 :
191 0 : void UVContSubTVI::insertToFieldSpecMap(const std::vector<int> &fieldIdxs,
192 : const InFitSpec &spec,
193 : InFitSpecMap &fitspec) const
194 : {
195 0 : for (const auto fid: fieldIdxs) {
196 0 : const auto fieldFound = fitspec.find(fid);
197 0 : if (fieldFound == fitspec.end()) {
198 0 : std::vector<InFitSpec> firstSpw = { spec };
199 0 : fitspec[fid] = firstSpw;
200 0 : } else {
201 0 : fitspec[fid].emplace_back(spec);
202 : }
203 :
204 : }
205 0 : }
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 0 : void UVContSubTVI::parseFieldSubDict(const Record &fieldRec,
216 : const std::vector<int> &fieldIdxs,
217 : InFitSpecMap &fitspec) const
218 : {
219 0 : std::set<unsigned int> spwsSeen;
220 0 : for (unsigned int spw=0; spw < fieldRec.nfields(); ++spw) {
221 0 : const auto spwRec = fieldRec.subRecord(RecordFieldId(spw));
222 0 : std::string spwStr = fieldRec.name(RecordFieldId(spw));
223 0 : const auto spwIdxs = stringToSPWIdxs(spwStr, false);
224 0 : for (const auto sid : spwIdxs) {
225 0 : if (spwsSeen.insert(sid).second == false) {
226 0 : throw AipsError("Error in fitspec. SPW " + to_string(sid) + " is given "
227 0 : "multiple times. Found for SPW subdict '" + spwStr +
228 0 : "' and other subdicts(s).");
229 : }
230 : }
231 : try {
232 0 : const std::string chanStr = spwRec.asString(RecordFieldId("chan"));
233 0 : if (!chanStr.empty()) {
234 0 : spwStr = "";
235 0 : for (const auto sid : spwIdxs) {
236 0 : const auto spw_chan = to_string(sid) + ":" + chanStr;
237 0 : if (spwStr.empty()) {
238 0 : spwStr = spw_chan;
239 : } else {
240 0 : spwStr += "," + spw_chan;
241 : }
242 0 : }
243 : }
244 0 : } catch (const AipsError &exc) {
245 0 : throw AipsError("Error trying to get chan value in subdict for spw "
246 0 : + spwStr + ", " + exc.what());
247 0 : }
248 :
249 0 : int polOrder = fitOrder_p;
250 : try {
251 0 : polOrder = spwRec.asInt(RecordFieldId("fitorder"));
252 0 : } catch (const AipsError &exc) {
253 0 : throw AipsError("Error trying to get fitorder value in subdict for spw "
254 0 : + spwStr + ", " + exc.what());
255 0 : }
256 0 : if (polOrder < 0) {
257 0 : throw AipsError("Fit order cannot be negative. Value given: " +
258 0 : to_string(polOrder));
259 : }
260 :
261 0 : const auto spec = std::make_pair(spwStr, polOrder);
262 0 : insertToFieldSpecMap(fieldIdxs, spec, fitspec);
263 0 : }
264 0 : }
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 0 : InFitSpecMap UVContSubTVI::parseDictFitSpec(const Record &configuration) const
276 : {
277 0 : const Record rawFieldFitspec = configuration.asRecord("fitspec");
278 0 : if (rawFieldFitspec.empty()) {
279 0 : throw AipsError("The dictionary 'fitspec' of fit specifications is empty!");
280 : }
281 :
282 0 : const auto maxMSField = getMaxMSFieldID();
283 0 : InFitSpecMap fitspec;
284 0 : std::set<unsigned int> fieldsSeen;
285 : // iterate through fields given in fitspec
286 0 : for (unsigned int rid=0; rid < rawFieldFitspec.nfields(); ++rid) {
287 0 : const std::string fieldsStr = rawFieldFitspec.name(RecordFieldId(rid));
288 0 : const auto fieldIdxs = stringToFieldIdxs(fieldsStr);
289 0 : for (const auto fid: fieldIdxs) {
290 : // check that the field is in the MS
291 0 : if (fid < 0 || static_cast<unsigned int>(fid) > maxMSField) {
292 0 : throw AipsError("Wrong field ID given: " +
293 0 : std::to_string(fid) +
294 : ". This MeasurementSet has field IDs between"
295 0 : " 0 and " + std::to_string(maxMSField));
296 : }
297 0 : if (fieldsSeen.insert(fid).second == false) {
298 0 : throw AipsError("Error in fitspec. Field " + to_string(fid) + " is given "
299 0 : "multiple times. Found for field subdict '" + fieldsStr +
300 0 : "' and other subdict(s).");
301 : }
302 : }
303 :
304 0 : std::string noneStr;
305 0 : bool subRecordIsString = true;
306 : try {
307 0 : noneStr = rawFieldFitspec.asString(RecordFieldId(rid));
308 0 : } catch (const AipsError &exc) {
309 0 : subRecordIsString = false;
310 0 : }
311 0 : if (subRecordIsString){
312 0 : if (noneStr != "NONE") {
313 0 : throw AipsError("Wrong value found in fitspec, field subdict for field: '"
314 0 : + fieldsStr + "', value: '" +
315 0 : noneStr + "'. Only 'NONE' is accepted (= do not fit for " +
316 0 : "any SPW).");
317 : }
318 : } else {
319 0 : const Record fieldRec = rawFieldFitspec.subRecord(RecordFieldId(rid));
320 0 : parseFieldSubDict(fieldRec, fieldIdxs, fitspec);
321 0 : }
322 0 : }
323 :
324 0 : return fitspec;
325 0 : }
326 :
327 0 : std::string fieldTextFromId(int idx)
328 : {
329 0 : std::string name;
330 0 : if (idx >= 0) {
331 0 : name = std::to_string(idx);
332 : } else {
333 0 : name = "all fields";
334 : }
335 0 : return name;
336 0 : }
337 :
338 0 : void UVContSubTVI::printInputFitSpec(const InFitSpecMap &fitspec) const
339 : {
340 0 : if (fitspec.size() > 0) {
341 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
342 : << "Fitspec (order and spw:line-free channel) specification is: "
343 0 : << LogIO::POST;
344 :
345 0 : for (const auto &elem: fitspec) {
346 0 : const auto &specs = elem.second;
347 0 : for (const auto &oneSpw : specs) {
348 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
349 0 : << " field: " << fieldTextFromId(elem.first)
350 0 : << ". " << oneSpw.second << ", '" << oneSpw.first << "'"
351 0 : << 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 0 : }
359 :
360 : // -----------------------------------------------------------------------
361 : //
362 : // -----------------------------------------------------------------------
363 0 : InFitSpecMap UVContSubTVI::parseConfiguration(const Record &configuration)
364 : {
365 0 : auto exists = configuration.fieldNumber ("fitorder");
366 0 : if (exists >= 0)
367 : {
368 0 : configuration.get (exists, fitOrder_p);
369 :
370 0 : if (fitOrder_p > 0)
371 : {
372 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
373 0 : << "Global/default fit order is " << fitOrder_p << LogIO::POST;
374 : }
375 : }
376 :
377 0 : exists = configuration.fieldNumber ("want_cont");
378 0 : if (exists >= 0)
379 : {
380 0 : configuration.get (exists, want_cont_p);
381 :
382 0 : 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 0 : const auto fitspec = parseFitSpec(configuration);
391 0 : printInputFitSpec(fitspec);
392 :
393 0 : exists = configuration.fieldNumber ("writemodel");
394 0 : if (exists >= 0)
395 : {
396 0 : configuration.get(exists, precalcModel_p);
397 0 : if (precalcModel_p)
398 : {
399 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
400 : << "Producing continuum estimate in the MODEL_DATA column"
401 0 : << LogIO::POST;
402 : }
403 : }
404 :
405 0 : exists = configuration.fieldNumber ("denoising_lib");
406 0 : if (exists >= 0)
407 : {
408 0 : configuration.get (exists, withDenoisingLib_p);
409 :
410 0 : if (withDenoisingLib_p)
411 : {
412 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
413 : << "Using GSL based multiparameter regression with linear "
414 0 : << "least-squares fitting" << LogIO::POST;
415 : }
416 : }
417 :
418 0 : exists = configuration.fieldNumber ("nthreads");
419 0 : 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 0 : exists = configuration.fieldNumber ("niter");
451 0 : if (exists >= 0)
452 : {
453 0 : configuration.get (exists, niter_p);
454 :
455 0 : 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 0 : exists = configuration.fieldNumber ("allowed_spws");
464 0 : if (exists >= 0)
465 : {
466 0 : String allowed;
467 0 : configuration.get (exists, allowed);
468 0 : allowedSpws_p = allowed;
469 0 : }
470 :
471 0 : 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 0 : void UVContSubTVI::spwInputChecks(const MeasurementSet *msVii,
485 : MSSelection &spwChan) const {
486 : // note: spwChan should be const (if getSpwList() were const)
487 0 : if (!allowedSpws_p.empty()) {
488 0 : const auto usedSpws = spwChan.getSpwList(msVii).tovector();
489 0 : MSSelection checker;
490 0 : checker.setSpwExpr(allowedSpws_p);
491 0 : const auto allowed = checker.getSpwList(msVii).tovector();
492 0 : for (const auto spw : usedSpws ) {
493 0 : if (std::find(allowed.begin(), allowed.end(), spw) == allowed.end()) {
494 0 : logger_p << LogIO::WARN << LogOrigin("UVContSubTVI", __FUNCTION__)
495 : << "SPW " << spw << " is used in fitspec but is not included in "
496 0 : << "the SPW selection ('" << allowedSpws_p << "'), please "
497 0 : << "double-check the spw and fitspec parameters." << LogIO::POST;
498 : }
499 : }
500 0 : }
501 0 : }
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 0 : unordered_map<int, vector<int>> UVContSubTVI::makeLineFreeChannelSelMap(std::string
513 : spwChanStr) const
514 : {
515 0 : MSSelection mssel;
516 0 : mssel.setSpwExpr(spwChanStr);
517 : // This will access the MS directly: far from ideal (CAS-11638) but no easy solution
518 0 : const auto msVii = &(inputVii_p->ms());
519 :
520 0 : spwInputChecks(msVii, mssel);
521 :
522 : // Create line-free channel map based on MSSelection channel ranges
523 0 : const auto spwchan = mssel.getChanList(msVii);
524 0 : const auto nSelections = spwchan.shape()[0];
525 0 : unordered_map<int,vector<int>> lineFreeChannelSelMap;
526 0 : for (uInt selection_i=0; selection_i<nSelections; ++selection_i)
527 : {
528 0 : auto spw = spwchan(selection_i, 0);
529 0 : if (lineFreeChannelSelMap.find(spw) == lineFreeChannelSelMap.end())
530 : {
531 0 : lineFreeChannelSelMap[spw].clear(); // Accessing the vector creates it
532 : }
533 :
534 0 : const auto channelStart = spwchan(selection_i, 1);
535 0 : const auto channelStop = spwchan(selection_i, 2);
536 0 : const auto channelStep = spwchan(selection_i, 3);
537 0 : for (auto inpChan=channelStart; inpChan<=channelStop; inpChan += channelStep)
538 : {
539 0 : lineFreeChannelSelMap[spw].push_back(inpChan);
540 : }
541 : }
542 :
543 0 : return lineFreeChannelSelMap;
544 0 : }
545 :
546 : // -----------------------------------------------------------------------
547 : //
548 : // -----------------------------------------------------------------------
549 0 : void UVContSubTVI::populatePerFieldSpec(int fieldID,
550 : const std::vector<InFitSpec> &fitSpecs)
551 : {
552 0 : std::unordered_map<int, int> perSpwFitOrder;
553 0 : std::string fullSpwStr = "";
554 0 : for (const auto &spec : fitSpecs) {
555 0 : if (fullSpwStr == "NONE") {
556 0 : break;
557 0 : } else if (fullSpwStr.length() == 0) {
558 0 : fullSpwStr = spec.first;
559 : } else {
560 : // MSSelection syntax spw:chan, SPWs separated by commas
561 0 : fullSpwStr += "," + spec.first;
562 : }
563 0 : auto spws = stringToSPWIdxs(spec.first, true);
564 0 : for (const auto sid : spws) {
565 0 : perSpwFitOrder[sid] = spec.second;
566 : }
567 0 : }
568 :
569 0 : auto lineFreeChannelSelMap = makeLineFreeChannelSelMap(fullSpwStr);
570 :
571 : // Create line-free channel mask, spw->(channel_mask, fit_order)
572 0 : unordered_map<int, FitSpec> lineFreeChannelMaskMap; // rename: fitSpecMap
573 0 : for (auto const spwInp: spwInpChanIdxMap_p)
574 : {
575 0 : const auto spw = spwInp.first;
576 0 : if (lineFreeChannelMaskMap.find(spw) == lineFreeChannelMaskMap.end())
577 : {
578 0 : 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 0 : lineFreeChannelMaskMap[spw].lineFreeChannelMask = Vector<bool>();
584 : } else {
585 0 : auto &mask = lineFreeChannelMaskMap[spw].lineFreeChannelMask;
586 0 : mask = Vector<bool>(spwInp.second.size(), true);
587 0 : for (uInt selChanIdx=0; selChanIdx<lineFreeChannelSelMap[spw].size();
588 : ++selChanIdx)
589 : {
590 0 : const auto selChan = lineFreeChannelSelMap[spw][selChanIdx];
591 0 : mask(selChan) = False;
592 : }
593 : }
594 :
595 0 : unsigned int fitOrder = fitOrder_p;
596 0 : const auto find = perSpwFitOrder.find(spw);
597 0 : if (find != perSpwFitOrder.end()) {
598 0 : fitOrder = find->second;
599 : }
600 0 : lineFreeChannelMaskMap[spw].fitOrder = fitOrder;
601 : }
602 0 : }
603 :
604 : // Poppulate per field map: field -> spw -> fit spec
605 : // emplace struct (linefreechannelmaskmap, + fitorder)
606 0 : perFieldSpecMap_p.emplace(fieldID, lineFreeChannelMaskMap);
607 0 : }
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 0 : void UVContSubTVI::fitSpecToPerFieldMap(const InFitSpecMap &fitspec)
615 : {
616 : // Process line-free channel specifications
617 0 : for (const auto item: fitspec) {
618 0 : unordered_map<int, FitSpec> fieldSpecMap;
619 : // Parse line-free channel selection using MSSelection syntax
620 0 : const auto fieldID = item.first;
621 0 : const auto &specs = item.second;
622 0 : bool noneFound = false;
623 0 : for (const auto spwSpec : specs) {
624 0 : const auto spwStr = spwSpec.first;
625 0 : const auto order = spwSpec.second;
626 0 : const auto fieldName = fieldTextFromId(fieldID);
627 0 : logger_p << LogIO::NORMAL << LogOrigin("UVContSubTVI", __FUNCTION__)
628 : << "Parsing fit spw:chan string and order for field: " << fieldName
629 0 : << ", spw/chan: '" << spwStr << "', order: " << order << LogIO::POST;
630 0 : 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 0 : }
638 :
639 0 : 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 0 : 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 0 : if (-1 == fieldID) {
649 0 : perFieldSpecMap_p.emplace(fieldID, fieldSpecMap);
650 : }
651 0 : continue;
652 : }
653 :
654 0 : populatePerFieldSpec(fieldID, specs);
655 0 : }
656 0 : }
657 :
658 : // -----------------------------------------------------------------------
659 : //
660 : // -----------------------------------------------------------------------
661 0 : void UVContSubTVI::initialize(const InFitSpecMap &fitspec)
662 : {
663 : // Populate nchan input-output maps
664 0 : for (auto spwInp: spwInpChanIdxMap_p)
665 : {
666 0 : auto spw = spwInp.first;
667 0 : spwOutChanNumMap_p[spw] = spwInp.second.size();
668 0 : }
669 :
670 0 : fitSpecToPerFieldMap(fitspec);
671 0 : }
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 0 : void UVContSubTVI::savePrecalcModel(const Cube<Complex> & origVis,
695 : const Cube<Complex> & contSubVis) const
696 : {
697 : // save it for visibilityModel
698 0 : if (precalcModel_p) {
699 : // Or otherwise the next assignment would fail a conform check
700 0 : if (precalcModelVis_p.shape() != contSubVis.shape()) {
701 0 : precalcModelVis_p.resize(contSubVis.shape());
702 : }
703 :
704 0 : precalcModelVis_p = origVis - contSubVis;
705 : }
706 0 : }
707 :
708 : // -----------------------------------------------------------------------
709 : //
710 : // -----------------------------------------------------------------------
711 0 : void UVContSubTVI::visibilityObserved (Cube<Complex> & vis) const
712 : {
713 : // Get input VisBuffer
714 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
715 :
716 : // Get weightSpectrum from sigmaSpectrum
717 0 : Cube<Float> weightSpFromSigmaSp;
718 0 : weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),False);
719 0 : weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
720 0 : arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
721 :
722 : // Transform data
723 0 : transformDataCube(vb->visCube(),weightSpFromSigmaSp,vis);
724 :
725 0 : savePrecalcModel(vb->visCube(), vis);
726 0 : }
727 :
728 : // -----------------------------------------------------------------------
729 : //
730 : // -----------------------------------------------------------------------
731 0 : void UVContSubTVI::visibilityCorrected (Cube<Complex> & vis) const
732 : {
733 : // Get input VisBuffer
734 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
735 :
736 : // Transform data
737 0 : transformDataCube(vb->visCubeCorrected(),vb->weightSpectrum(),vis);
738 :
739 0 : savePrecalcModel(vb->visCubeCorrected(), vis);
740 0 : }
741 :
742 : // -----------------------------------------------------------------------
743 : //
744 : // -----------------------------------------------------------------------
745 0 : void UVContSubTVI::visibilityModel (Cube<Complex> & vis) const
746 : {
747 : // Get input VisBuffer
748 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
749 :
750 : // from visiblityObserved we have calculated the polynomial subtraction
751 0 : if (precalcModel_p && precalcModelVis_p.shape() == vb->visCubeModel().shape()) {
752 0 : vis = precalcModelVis_p;
753 0 : return;
754 : }
755 :
756 : // Transform data
757 0 : transformDataCube(vb->visCubeModel(),vb->weightSpectrum(),vis);
758 :
759 0 : return;
760 : }
761 :
762 : // -----------------------------------------------------------------------
763 : //
764 : // -----------------------------------------------------------------------
765 0 : void UVContSubTVI::result(Record &res) const
766 : {
767 0 : auto acc = result_p.getAccumulatedResult();
768 0 : res.defineRecord("uvcontsub", acc);
769 0 : }
770 :
771 : // -----------------------------------------------------------------------
772 : //
773 : // -----------------------------------------------------------------------
774 0 : 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 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
780 :
781 0 : auto fieldID = vb->fieldId()[0];
782 : // First check the "all fields" (no per-individual field list in fitspec)
783 0 : auto fieldMapIt = perFieldSpecMap_p.find(-1);
784 0 : if (fieldMapIt == perFieldSpecMap_p.end()) {
785 : // otherwise, check for individual fields (coming from a fitspec list)
786 0 : fieldMapIt = perFieldSpecMap_p.find(fieldID);
787 0 : if (fieldMapIt == perFieldSpecMap_p.end()) {
788 : // This was a "NONE" = no-subtraction for this field ==> no-op
789 0 : outputVis = inputVis;
790 0 : return;
791 : }
792 : }
793 :
794 0 : Int spwId = vb->spectralWindows()(0);
795 :
796 : // Get input line-free channel mask and fitorder
797 0 : Vector<bool> *lineFreeChannelMask = nullptr;
798 0 : auto fitOrder = fitOrder_p;
799 0 : auto spwMap = fieldMapIt->second;
800 0 : auto maskLookup = spwMap.find(spwId);
801 0 : if (maskLookup != spwMap.end())
802 : {
803 0 : if (maskLookup->second.lineFreeChannelMask.nelements() == 0) {
804 : // This was a non-selected SPW in a non-empty SPW selection string ==> no-op
805 0 : outputVis = inputVis;
806 0 : return;
807 : } else {
808 0 : lineFreeChannelMask = &(maskLookup->second.lineFreeChannelMask);
809 0 : fitOrder = maskLookup->second.fitOrder;
810 : }
811 : }
812 :
813 : // Get polynomial model for this SPW (depends on number of channels and gridding)
814 0 : const auto freqIter = inputFrequencyMap_p.find(spwId);
815 0 : 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 0 : reset(new denoising::GslPolynomialModel<double>(inputFrequencies, fitOrder));
820 :
821 : // Reshape output data before passing it to the DataCubeHolder
822 0 : outputVis.resize(getVisBuffer()->getShape(),False);
823 :
824 : // Get input flag Cube
825 0 : const Cube<bool> &flagCube = vb->flagCube();
826 :
827 : // Transform data
828 0 : 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 0 : auto scanID = vb->scan()[0];
855 0 : uInt nCorrs = vb->getShape()(0);
856 0 : for (uInt corrIdx=0; corrIdx < nCorrs; corrIdx++)
857 : {
858 : Complex chisq =
859 0 : transformDataCore(inputFrequencyMap_p[spwId].get(), lineFreeChannelMask,
860 : inputVis, flagCube, inputWeight, outputVis, corrIdx);
861 0 : result_p.addOneFit(fieldID, scanID, spwId, (int)corrIdx, chisq);
862 : }
863 : }
864 0 : }
865 :
866 : // -----------------------------------------------------------------------
867 : //
868 : // -----------------------------------------------------------------------
869 0 : 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 0 : DataCubeMap inputData;
879 0 : DataCubeHolder<T> inputVisCubeHolder(inputVis);
880 0 : DataCubeHolder<bool> inputFlagCubeHolder(inputFlags);
881 0 : DataCubeHolder<Float> inputWeightsCubeHolder(inputWeight);
882 0 : inputData.add(MS::DATA,inputVisCubeHolder);
883 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
884 0 : inputData.add(MS::WEIGHT_SPECTRUM,inputWeightsCubeHolder);
885 :
886 : // Gather output data
887 0 : DataCubeMap outputData;
888 0 : DataCubeHolder<T> outputVisCubeHolder(outputVis);
889 0 : outputData.add(MS::DATA,outputVisCubeHolder);
890 :
891 0 : Complex chisq;
892 0 : 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 0 : if (withDenoisingLib_p)
912 : {
913 0 : UVContSubtractionDenoisingKernel<T> kernel(model,niter_p,lineFreeChannelMask);
914 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
915 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
916 0 : chisq = kernel.getChiSquared();
917 0 : }
918 : else
919 : {
920 0 : UVContSubtractionKernel<T> kernel(model,lineFreeChannelMask);
921 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
922 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
923 0 : chisq = kernel.getChiSquared();
924 0 : }
925 : }
926 :
927 0 : return chisq;
928 0 : }
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 0 : UVContSubTVILayerFactory::UVContSubTVILayerFactory(Record &configuration) :
965 0 : ViiLayerFactory()
966 : {
967 0 : configuration_p = configuration;
968 0 : }
969 :
970 : ViImplementation2*
971 0 : UVContSubTVILayerFactory::createInstance(ViImplementation2* vii0) const
972 : {
973 : // Make the UVContSubTVI, using supplied ViImplementation2, and return it
974 0 : ViImplementation2 *vii = new UVContSubTVI(vii0, configuration_p);
975 0 : return vii;
976 : }
977 :
978 : //////////////////////////////////////////////////////////////////////////
979 : // UVContSubTransformEngine class
980 : //////////////////////////////////////////////////////////////////////////
981 :
982 : // -----------------------------------------------------------------------
983 : //
984 : // -----------------------------------------------------------------------
985 0 : template<class T> UVContSubTransformEngine<T>::UVContSubTransformEngine( UVContSubKernel<T> *kernel,
986 : DataCubeMap *inputData,
987 : DataCubeMap *outputData):
988 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
989 : {
990 0 : uvContSubKernel_p = kernel;
991 0 : }
992 :
993 : // -----------------------------------------------------------------------
994 : //
995 : // -----------------------------------------------------------------------
996 0 : template<class T> void UVContSubTransformEngine<T>::transform( )
997 : {
998 0 : uvContSubKernel_p->setDebug(debug_p);
999 0 : uvContSubKernel_p->kernel(inputData_p,outputData_p);
1000 0 : }
1001 :
1002 : //////////////////////////////////////////////////////////////////////////
1003 : // UVContSubKernel class
1004 : //////////////////////////////////////////////////////////////////////////
1005 :
1006 : // -----------------------------------------------------------------------
1007 : //
1008 : // -----------------------------------------------------------------------
1009 0 : template<class T> UVContSubKernel<T>::UVContSubKernel( denoising::GslPolynomialModel<Double> *model,
1010 0 : Vector<bool> *lineFreeChannelMask)
1011 : {
1012 0 : model_p = model;
1013 0 : fitOrder_p = model_p->ncomponents()-1;
1014 0 : freqPows_p.reference(model_p->getModelMatrix());
1015 0 : frequencies_p.reference(model_p->getLinearComponentFloat());
1016 :
1017 0 : lineFreeChannelMask_p = lineFreeChannelMask ? lineFreeChannelMask : nullptr;
1018 0 : debug_p = False;
1019 0 : }
1020 :
1021 : // -----------------------------------------------------------------------
1022 : //
1023 : // -----------------------------------------------------------------------
1024 0 : template<class T> void UVContSubKernel<T>::kernel(DataCubeMap *inputData,
1025 : DataCubeMap *outputData)
1026 : {
1027 : // Get input/output data
1028 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1029 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1030 :
1031 : // Apply line-free channel mask
1032 0 : Vector<bool> &inputFlags = inputData->getVector<bool>(MS::FLAG);
1033 0 : if (lineFreeChannelMask_p != nullptr)
1034 0 : inputFlags |= *lineFreeChannelMask_p;
1035 :
1036 : // Calculate number of valid data points and adapt fit
1037 0 : size_t validPoints = nfalse(inputFlags);
1038 0 : if (validPoints > 0)
1039 : {
1040 0 : bool restoreDefaultPoly = False;
1041 0 : uInt tmpFitOrder = fitOrder_p;
1042 :
1043 : // Reduce fit order to match number of valid points
1044 0 : if (validPoints <= fitOrder_p)
1045 : {
1046 0 : changeFitOrder(validPoints-1);
1047 0 : restoreDefaultPoly = true;
1048 : }
1049 :
1050 : // Get weights
1051 0 : Vector<Float> &inputWeight = inputData->getVector<Float>(MS::WEIGHT_SPECTRUM);
1052 :
1053 : // Convert flags to mask
1054 0 : Vector<bool> mask = !inputFlags;
1055 :
1056 : // Calculate and subtract continuum
1057 0 : chisq_p = kernelCore(inputVector,mask,inputWeight,outputVector);
1058 :
1059 : // Go back to default fit order to match number of valid points
1060 0 : if (restoreDefaultPoly)
1061 : {
1062 0 : changeFitOrder(tmpFitOrder);
1063 : }
1064 0 : }
1065 : else
1066 : {
1067 0 : defaultKernel(inputVector,outputVector);
1068 : }
1069 0 : }
1070 :
1071 :
1072 : //////////////////////////////////////////////////////////////////////////
1073 : // UVContSubtractionKernel class
1074 : //////////////////////////////////////////////////////////////////////////
1075 :
1076 : // -----------------------------------------------------------------------
1077 : //
1078 : // -----------------------------------------------------------------------
1079 0 : template<class T> UVContSubtractionKernel<T>::UVContSubtractionKernel( denoising::GslPolynomialModel<Double>* model,
1080 : Vector<bool> *lineFreeChannelMask):
1081 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1082 : {
1083 0 : changeFitOrder(fitOrder_p);
1084 0 : }
1085 :
1086 : // -----------------------------------------------------------------------
1087 : //
1088 : // -----------------------------------------------------------------------
1089 0 : template<class T> void UVContSubtractionKernel<T>::changeFitOrder(size_t order)
1090 : {
1091 0 : fitOrder_p = order;
1092 0 : Polynomial<AutoDiff<Float> > poly(order);
1093 0 : fitter_p.setFunction(poly); // poly It is cloned
1094 0 : fitter_p.asWeight(true);
1095 :
1096 0 : return;
1097 0 : }
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 0 : 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 0 : Vector<Float> realCoeff;
1130 0 : Vector<Float> imagCoeff;
1131 0 : realCoeff = fitter_p.fit(frequencies_p, real(inputVector), inputWeights, &inputFlags);
1132 0 : double realChisq = fitter_p.chiSquare();
1133 0 : imagCoeff = fitter_p.fit(frequencies_p, imag(inputVector), inputWeights, &inputFlags);
1134 0 : double imagChisq = fitter_p.chiSquare();
1135 :
1136 : // Fill output data
1137 0 : outputVector = inputVector;
1138 0 : outputVector -= Complex(realCoeff(0),imagCoeff(0));
1139 0 : 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 0 : 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 0 : return Complex(realChisq, imagChisq);
1161 0 : }
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 0 : template<class T> UVContSubtractionDenoisingKernel<T>::UVContSubtractionDenoisingKernel(denoising::GslPolynomialModel<Double>* model,
1315 : size_t nIter,
1316 : Vector<bool> *lineFreeChannelMask):
1317 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1318 : {
1319 0 : fitter_p.resetModel(*model);
1320 0 : fitter_p.setNIter(nIter);
1321 0 : }
1322 :
1323 : // -----------------------------------------------------------------------
1324 : //
1325 : // -----------------------------------------------------------------------
1326 0 : template<class T> void UVContSubtractionDenoisingKernel<T>::changeFitOrder(size_t order)
1327 : {
1328 0 : fitOrder_p = order;
1329 0 : fitter_p.resetNComponents(order+1);
1330 0 : return;
1331 : }
1332 :
1333 :
1334 : // -----------------------------------------------------------------------
1335 : //
1336 : // -----------------------------------------------------------------------
1337 0 : template<class T> void UVContSubtractionDenoisingKernel<T>::defaultKernel( Vector<T> &inputVector,
1338 : Vector<T> &outputVector)
1339 : {
1340 0 : outputVector = inputVector;
1341 0 : return;
1342 : }
1343 :
1344 : // -----------------------------------------------------------------------
1345 : //
1346 : // -----------------------------------------------------------------------
1347 0 : template<class T> Complex UVContSubtractionDenoisingKernel<T>::kernelCore( Vector<T> &inputVector,
1348 : Vector<bool> &inputFlags,
1349 : Vector<Float> &inputWeights,
1350 : Vector<T> &outputVector)
1351 0 : { fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1352 0 : vector<T> coeff;
1353 0 : Complex chisq;
1354 0 : tie(coeff, chisq) = fitter_p.calcFitCoeff(inputVector);
1355 :
1356 0 : Vector<T> model(outputVector.size());
1357 0 : fitter_p.calcFitModel(model);
1358 :
1359 0 : outputVector = inputVector;
1360 0 : 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 0 : return chisq;
1393 0 : }
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 :
|
