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 : }
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 Vector<Double> &inputFrequencies = vb->getFrequencies(0);
815 : // Could perhaps have a per field-spw pair map to avoid re-allocations - But can be big
816 : // (n_fields X n_spws X n_chans) for many fields, many SPWs MSs
817 : inputFrequencyMap_p[spwId].
818 0 : reset(new denoising::GslPolynomialModel<double>(inputFrequencies, fitOrder));
819 :
820 : // Reshape output data before passing it to the DataCubeHolder
821 0 : outputVis.resize(getVisBuffer()->getShape(),False);
822 :
823 : // Get input flag Cube
824 0 : const Cube<bool> &flagCube = vb->flagCube();
825 :
826 : // Transform data
827 0 : if (nThreads_p > 1)
828 : {
829 : #ifdef _OPENMP
830 :
831 0 : uInt nCorrs = vb->getShape()(0);
832 0 : if (nCorrs < nThreads_p)
833 : {
834 0 : omp_set_num_threads(nCorrs);
835 : }
836 : else
837 : {
838 0 : omp_set_num_threads(nThreads_p);
839 : }
840 :
841 0 : #pragma omp parallel for
842 : for (uInt corrIdx=0; corrIdx < nCorrs; corrIdx++)
843 : {
844 : transformDataCore(inputFrequencyMap_p[spwId].get(),lineFreeChannelMask,
845 : inputVis,flagCube,inputWeight,outputVis,corrIdx);
846 : }
847 :
848 0 : omp_set_num_threads(nThreads_p);
849 : #endif
850 : }
851 : else
852 : {
853 0 : auto scanID = vb->scan()[0];
854 0 : uInt nCorrs = vb->getShape()(0);
855 0 : for (uInt corrIdx=0; corrIdx < nCorrs; corrIdx++)
856 : {
857 : Complex chisq =
858 0 : transformDataCore(inputFrequencyMap_p[spwId].get(), lineFreeChannelMask,
859 : inputVis, flagCube, inputWeight, outputVis, corrIdx);
860 0 : result_p.addOneFit(fieldID, scanID, spwId, (int)corrIdx, chisq);
861 : }
862 : }
863 0 : }
864 :
865 : // -----------------------------------------------------------------------
866 : //
867 : // -----------------------------------------------------------------------
868 0 : template<class T> Complex UVContSubTVI::transformDataCore(denoising::GslPolynomialModel<Double>* model,
869 : Vector<bool> *lineFreeChannelMask,
870 : const Cube<T> &inputVis,
871 : const Cube<bool> &inputFlags,
872 : const Cube<Float> &inputWeight,
873 : Cube<T> &outputVis,
874 : Int parallelCorrAxis) const
875 : {
876 : // Gather input data
877 0 : DataCubeMap inputData;
878 0 : DataCubeHolder<T> inputVisCubeHolder(inputVis);
879 0 : DataCubeHolder<bool> inputFlagCubeHolder(inputFlags);
880 0 : DataCubeHolder<Float> inputWeightsCubeHolder(inputWeight);
881 0 : inputData.add(MS::DATA,inputVisCubeHolder);
882 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
883 0 : inputData.add(MS::WEIGHT_SPECTRUM,inputWeightsCubeHolder);
884 :
885 : // Gather output data
886 0 : DataCubeMap outputData;
887 0 : DataCubeHolder<T> outputVisCubeHolder(outputVis);
888 0 : outputData.add(MS::DATA,outputVisCubeHolder);
889 :
890 0 : Complex chisq;
891 0 : if (want_cont_p)
892 : {
893 0 : if (withDenoisingLib_p)
894 : {
895 0 : UVContEstimationDenoisingKernel<T> kernel(model,niter_p,lineFreeChannelMask);
896 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
897 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
898 0 : chisq = kernel.getChiSquared();
899 0 : }
900 : else
901 : {
902 0 : UVContEstimationKernel<T> kernel(model,lineFreeChannelMask);
903 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
904 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
905 0 : chisq = kernel.getChiSquared();
906 0 : }
907 : }
908 : else
909 : {
910 0 : if (withDenoisingLib_p)
911 : {
912 0 : UVContSubtractionDenoisingKernel<T> kernel(model,niter_p,lineFreeChannelMask);
913 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
914 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
915 0 : chisq = kernel.getChiSquared();
916 0 : }
917 : else
918 : {
919 0 : UVContSubtractionKernel<T> kernel(model,lineFreeChannelMask);
920 0 : UVContSubTransformEngine<T> transformer(&kernel,&inputData,&outputData);
921 0 : transformFreqAxis2(inputVis.shape(),transformer,parallelCorrAxis);
922 0 : chisq = kernel.getChiSquared();
923 0 : }
924 : }
925 :
926 0 : return chisq;
927 0 : }
928 :
929 : //////////////////////////////////////////////////////////////////////////
930 : // UVContSubTVIFactory class
931 : //////////////////////////////////////////////////////////////////////////
932 :
933 : // -----------------------------------------------------------------------
934 : //
935 : // -----------------------------------------------------------------------
936 0 : UVContSubTVIFactory::UVContSubTVIFactory ( Record &configuration,
937 0 : ViImplementation2 *inputVii)
938 : {
939 0 : inputVii_p = inputVii;
940 0 : configuration_p = configuration;
941 0 : }
942 :
943 : // -----------------------------------------------------------------------
944 : //
945 : // -----------------------------------------------------------------------
946 0 : vi::ViImplementation2 * UVContSubTVIFactory::createVi(VisibilityIterator2 *) const
947 : {
948 0 : return new UVContSubTVI(inputVii_p,configuration_p);
949 : }
950 :
951 : // -----------------------------------------------------------------------
952 : //
953 : // -----------------------------------------------------------------------
954 0 : vi::ViImplementation2 * UVContSubTVIFactory::createVi() const
955 : {
956 0 : return new UVContSubTVI(inputVii_p,configuration_p);
957 : }
958 :
959 : //////////////////////////////////////////////////////////////////////////
960 : // UVContSubTVILayerFactory class
961 : //////////////////////////////////////////////////////////////////////////
962 :
963 0 : UVContSubTVILayerFactory::UVContSubTVILayerFactory(Record &configuration) :
964 0 : ViiLayerFactory()
965 : {
966 0 : configuration_p = configuration;
967 0 : }
968 :
969 : ViImplementation2*
970 0 : UVContSubTVILayerFactory::createInstance(ViImplementation2* vii0) const
971 : {
972 : // Make the UVContSubTVI, using supplied ViImplementation2, and return it
973 0 : ViImplementation2 *vii = new UVContSubTVI(vii0, configuration_p);
974 0 : return vii;
975 : }
976 :
977 : //////////////////////////////////////////////////////////////////////////
978 : // UVContSubTransformEngine class
979 : //////////////////////////////////////////////////////////////////////////
980 :
981 : // -----------------------------------------------------------------------
982 : //
983 : // -----------------------------------------------------------------------
984 0 : template<class T> UVContSubTransformEngine<T>::UVContSubTransformEngine( UVContSubKernel<T> *kernel,
985 : DataCubeMap *inputData,
986 : DataCubeMap *outputData):
987 0 : FreqAxisTransformEngine2<T>(inputData,outputData)
988 : {
989 0 : uvContSubKernel_p = kernel;
990 0 : }
991 :
992 : // -----------------------------------------------------------------------
993 : //
994 : // -----------------------------------------------------------------------
995 0 : template<class T> void UVContSubTransformEngine<T>::transform( )
996 : {
997 0 : uvContSubKernel_p->setDebug(debug_p);
998 0 : uvContSubKernel_p->kernel(inputData_p,outputData_p);
999 0 : }
1000 :
1001 : //////////////////////////////////////////////////////////////////////////
1002 : // UVContSubKernel class
1003 : //////////////////////////////////////////////////////////////////////////
1004 :
1005 : // -----------------------------------------------------------------------
1006 : //
1007 : // -----------------------------------------------------------------------
1008 0 : template<class T> UVContSubKernel<T>::UVContSubKernel( denoising::GslPolynomialModel<Double> *model,
1009 0 : Vector<bool> *lineFreeChannelMask)
1010 : {
1011 0 : model_p = model;
1012 0 : fitOrder_p = model_p->ncomponents()-1;
1013 0 : freqPows_p.reference(model_p->getModelMatrix());
1014 0 : frequencies_p.reference(model_p->getLinearComponentFloat());
1015 :
1016 0 : lineFreeChannelMask_p = lineFreeChannelMask ? lineFreeChannelMask : nullptr;
1017 0 : debug_p = False;
1018 0 : }
1019 :
1020 : // -----------------------------------------------------------------------
1021 : //
1022 : // -----------------------------------------------------------------------
1023 0 : template<class T> void UVContSubKernel<T>::kernel(DataCubeMap *inputData,
1024 : DataCubeMap *outputData)
1025 : {
1026 : // Get input/output data
1027 0 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1028 0 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1029 :
1030 : // Apply line-free channel mask
1031 0 : Vector<bool> &inputFlags = inputData->getVector<bool>(MS::FLAG);
1032 0 : if (lineFreeChannelMask_p != nullptr)
1033 0 : inputFlags |= *lineFreeChannelMask_p;
1034 :
1035 : // Calculate number of valid data points and adapt fit
1036 0 : size_t validPoints = nfalse(inputFlags);
1037 0 : if (validPoints > 0)
1038 : {
1039 0 : bool restoreDefaultPoly = False;
1040 0 : uInt tmpFitOrder = fitOrder_p;
1041 :
1042 : // Reduce fit order to match number of valid points
1043 0 : if (validPoints <= fitOrder_p)
1044 : {
1045 0 : changeFitOrder(validPoints-1);
1046 0 : restoreDefaultPoly = true;
1047 : }
1048 :
1049 : // Get weights
1050 0 : Vector<Float> &inputWeight = inputData->getVector<Float>(MS::WEIGHT_SPECTRUM);
1051 :
1052 : // Convert flags to mask
1053 0 : Vector<bool> mask = !inputFlags;
1054 :
1055 : // Calculate and subtract continuum
1056 0 : chisq_p = kernelCore(inputVector,mask,inputWeight,outputVector);
1057 :
1058 : // Go back to default fit order to match number of valid points
1059 0 : if (restoreDefaultPoly)
1060 : {
1061 0 : changeFitOrder(tmpFitOrder);
1062 : }
1063 0 : }
1064 : else
1065 : {
1066 0 : defaultKernel(inputVector,outputVector);
1067 : }
1068 0 : }
1069 :
1070 :
1071 : //////////////////////////////////////////////////////////////////////////
1072 : // UVContSubtractionKernel class
1073 : //////////////////////////////////////////////////////////////////////////
1074 :
1075 : // -----------------------------------------------------------------------
1076 : //
1077 : // -----------------------------------------------------------------------
1078 0 : template<class T> UVContSubtractionKernel<T>::UVContSubtractionKernel( denoising::GslPolynomialModel<Double>* model,
1079 : Vector<bool> *lineFreeChannelMask):
1080 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1081 : {
1082 0 : changeFitOrder(fitOrder_p);
1083 0 : }
1084 :
1085 : // -----------------------------------------------------------------------
1086 : //
1087 : // -----------------------------------------------------------------------
1088 0 : template<class T> void UVContSubtractionKernel<T>::changeFitOrder(size_t order)
1089 : {
1090 0 : fitOrder_p = order;
1091 0 : Polynomial<AutoDiff<Float> > poly(order);
1092 0 : fitter_p.setFunction(poly); // poly It is cloned
1093 0 : fitter_p.asWeight(true);
1094 :
1095 0 : return;
1096 0 : }
1097 :
1098 :
1099 : // -----------------------------------------------------------------------
1100 : //
1101 : // -----------------------------------------------------------------------
1102 0 : template<class T> void UVContSubtractionKernel<T>::defaultKernel( Vector<Complex> &inputVector,
1103 : Vector<Complex> &outputVector)
1104 : {
1105 0 : outputVector = inputVector;
1106 0 : return;
1107 : }
1108 :
1109 : // -----------------------------------------------------------------------
1110 : //
1111 : // -----------------------------------------------------------------------
1112 0 : template<class T> void UVContSubtractionKernel<T>::defaultKernel( Vector<Float> &inputVector,
1113 : Vector<Float> &outputVector)
1114 : {
1115 0 : outputVector = inputVector;
1116 0 : return;
1117 : }
1118 :
1119 : // -----------------------------------------------------------------------
1120 : //
1121 : // -----------------------------------------------------------------------
1122 0 : template<class T> Complex UVContSubtractionKernel<T>::kernelCore( Vector<Complex> &inputVector,
1123 : Vector<bool> &inputFlags,
1124 : Vector<Float> &inputWeights,
1125 : Vector<Complex> &outputVector)
1126 : {
1127 : // Fit for imaginary and real components separately
1128 0 : Vector<Float> realCoeff;
1129 0 : Vector<Float> imagCoeff;
1130 0 : realCoeff = fitter_p.fit(frequencies_p, real(inputVector), inputWeights, &inputFlags);
1131 0 : double realChisq = fitter_p.chiSquare();
1132 0 : imagCoeff = fitter_p.fit(frequencies_p, imag(inputVector), inputWeights, &inputFlags);
1133 0 : double imagChisq = fitter_p.chiSquare();
1134 :
1135 : // Fill output data
1136 0 : outputVector = inputVector;
1137 0 : outputVector -= Complex(realCoeff(0),imagCoeff(0));
1138 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1139 : {
1140 0 : Complex coeff(realCoeff(order_idx),imagCoeff(order_idx));
1141 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1142 : {
1143 0 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff;
1144 : }
1145 : }
1146 :
1147 0 : if (debug_p)
1148 : {
1149 0 : LogIO logger;
1150 0 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1151 0 : logger << "realCoeff =" << realCoeff << LogIO::POST;
1152 0 : logger << "imagCoeff =" << imagCoeff << LogIO::POST;
1153 0 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1154 0 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1155 0 : logger << "inputVector =" << inputVector << LogIO::POST;
1156 0 : logger << "outputVector =" << outputVector << LogIO::POST;
1157 0 : }
1158 :
1159 0 : return Complex(realChisq, imagChisq);
1160 0 : }
1161 :
1162 : // -----------------------------------------------------------------------
1163 : //
1164 : // -----------------------------------------------------------------------
1165 0 : template<class T> Complex UVContSubtractionKernel<T>::kernelCore( Vector<Float> &inputVector,
1166 : Vector<bool> &inputFlags,
1167 : Vector<Float> &inputWeights,
1168 : Vector<Float> &outputVector)
1169 : {
1170 : // Fit model
1171 0 : Vector<Float> coeff;
1172 0 : coeff = fitter_p.fit(frequencies_p, inputVector, inputWeights, &inputFlags);
1173 0 : const double chisq = fitter_p.chiSquare();
1174 :
1175 : // Fill output data
1176 0 : outputVector = inputVector;
1177 0 : outputVector -= coeff(0);
1178 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1179 : {
1180 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1181 : {
1182 0 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1183 : }
1184 : }
1185 :
1186 0 : if (debug_p)
1187 : {
1188 0 : LogIO logger;
1189 0 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1190 0 : logger << "coeff =" << coeff << LogIO::POST;
1191 0 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1192 0 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1193 0 : logger << "inputVector =" << inputVector << LogIO::POST;
1194 0 : logger << "outputVector =" << outputVector << LogIO::POST;
1195 0 : }
1196 :
1197 0 : return Complex(chisq, chisq);
1198 0 : }
1199 :
1200 :
1201 : //////////////////////////////////////////////////////////////////////////
1202 : // UVContEstimationKernel class
1203 : //////////////////////////////////////////////////////////////////////////
1204 :
1205 : // -----------------------------------------------------------------------
1206 : //
1207 : // -----------------------------------------------------------------------
1208 0 : template<class T> UVContEstimationKernel<T>::UVContEstimationKernel( denoising::GslPolynomialModel<Double>* model,
1209 : Vector<bool> *lineFreeChannelMask):
1210 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1211 : {
1212 0 : changeFitOrder(fitOrder_p);
1213 0 : }
1214 :
1215 : // -----------------------------------------------------------------------
1216 : //
1217 : // -----------------------------------------------------------------------
1218 0 : template<class T> void UVContEstimationKernel<T>::changeFitOrder(size_t order)
1219 : {
1220 0 : fitOrder_p = order;
1221 0 : Polynomial<AutoDiff<Float> > poly(order);
1222 0 : fitter_p.setFunction(poly); // poly It is cloned
1223 0 : fitter_p.asWeight(true);
1224 :
1225 0 : return;
1226 0 : }
1227 :
1228 : // -----------------------------------------------------------------------
1229 : //
1230 : // -----------------------------------------------------------------------
1231 0 : template<class T> void UVContEstimationKernel<T>::defaultKernel(Vector<Complex> &,
1232 : Vector<Complex> &outputVector)
1233 : {
1234 0 : outputVector = 0;
1235 0 : return;
1236 : }
1237 :
1238 : // -----------------------------------------------------------------------
1239 : //
1240 : // -----------------------------------------------------------------------
1241 0 : template<class T> void UVContEstimationKernel<T>::defaultKernel(Vector<Float> &,
1242 : Vector<Float> &outputVector)
1243 : {
1244 0 : outputVector = 0;
1245 0 : return;
1246 : }
1247 :
1248 : // -----------------------------------------------------------------------
1249 : //
1250 : // -----------------------------------------------------------------------
1251 0 : template<class T> Complex UVContEstimationKernel<T>::kernelCore(Vector<Complex> &inputVector,
1252 : Vector<bool> &inputFlags,
1253 : Vector<Float> &inputWeights,
1254 : Vector<Complex> &outputVector
1255 : )
1256 : {
1257 : // Fit for imaginary and real components separately
1258 0 : Vector<Float> realCoeff;
1259 0 : Vector<Float> imagCoeff;
1260 0 : realCoeff = fitter_p.fit(frequencies_p, real(inputVector), inputWeights, &inputFlags);
1261 0 : double realChisq = fitter_p.chiSquare();
1262 :
1263 0 : imagCoeff = fitter_p.fit(frequencies_p, imag(inputVector), inputWeights, &inputFlags);
1264 0 : double imagChisq = fitter_p.chiSquare();
1265 :
1266 : // Fill output data
1267 0 : outputVector = Complex(realCoeff(0),imagCoeff(0));
1268 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1269 : {
1270 0 : Complex coeff(realCoeff(order_idx),imagCoeff(order_idx));
1271 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1272 : {
1273 0 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff;
1274 : }
1275 : }
1276 :
1277 0 : return Complex(realChisq, imagChisq);
1278 0 : }
1279 :
1280 : // -----------------------------------------------------------------------
1281 : //
1282 : // -----------------------------------------------------------------------
1283 0 : template<class T> Complex UVContEstimationKernel<T>::kernelCore( Vector<Float> &inputVector,
1284 : Vector<bool> &inputFlags,
1285 : Vector<Float> &inputWeights,
1286 : Vector<Float> &outputVector)
1287 : {
1288 : // Fit model
1289 0 : Vector<Float> coeff;
1290 0 : coeff = fitter_p.fit(frequencies_p, inputVector, inputWeights, &inputFlags);
1291 0 : double chisq = fitter_p.chiSquare();
1292 :
1293 : // Fill output data
1294 0 : outputVector = coeff(0);
1295 0 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1296 : {
1297 0 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1298 : {
1299 0 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1300 : }
1301 : }
1302 :
1303 0 : return Complex(chisq, chisq);
1304 0 : }
1305 :
1306 : //////////////////////////////////////////////////////////////////////////
1307 : // UVContSubtractionDenoisingKernel class
1308 : //////////////////////////////////////////////////////////////////////////
1309 :
1310 : // -----------------------------------------------------------------------
1311 : //
1312 : // -----------------------------------------------------------------------
1313 0 : template<class T> UVContSubtractionDenoisingKernel<T>::UVContSubtractionDenoisingKernel(denoising::GslPolynomialModel<Double>* model,
1314 : size_t nIter,
1315 : Vector<bool> *lineFreeChannelMask):
1316 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1317 : {
1318 0 : fitter_p.resetModel(*model);
1319 0 : fitter_p.setNIter(nIter);
1320 0 : }
1321 :
1322 : // -----------------------------------------------------------------------
1323 : //
1324 : // -----------------------------------------------------------------------
1325 0 : template<class T> void UVContSubtractionDenoisingKernel<T>::changeFitOrder(size_t order)
1326 : {
1327 0 : fitOrder_p = order;
1328 0 : fitter_p.resetNComponents(order+1);
1329 0 : return;
1330 : }
1331 :
1332 :
1333 : // -----------------------------------------------------------------------
1334 : //
1335 : // -----------------------------------------------------------------------
1336 0 : template<class T> void UVContSubtractionDenoisingKernel<T>::defaultKernel( Vector<T> &inputVector,
1337 : Vector<T> &outputVector)
1338 : {
1339 0 : outputVector = inputVector;
1340 0 : return;
1341 : }
1342 :
1343 : // -----------------------------------------------------------------------
1344 : //
1345 : // -----------------------------------------------------------------------
1346 0 : template<class T> Complex UVContSubtractionDenoisingKernel<T>::kernelCore( Vector<T> &inputVector,
1347 : Vector<bool> &inputFlags,
1348 : Vector<Float> &inputWeights,
1349 : Vector<T> &outputVector)
1350 0 : { fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1351 0 : vector<T> coeff;
1352 0 : Complex chisq;
1353 0 : tie(coeff, chisq) = fitter_p.calcFitCoeff(inputVector);
1354 :
1355 0 : Vector<T> model(outputVector.size());
1356 0 : fitter_p.calcFitModel(model);
1357 :
1358 0 : outputVector = inputVector;
1359 0 : outputVector -= model;
1360 :
1361 : /*
1362 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1363 : Vector<T> coeff = fitter_p.calcFitCoeff(inputVector);
1364 :
1365 : // Fill output data
1366 : outputVector = inputVector;
1367 : outputVector -= coeff(0);
1368 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1369 : {
1370 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1371 : {
1372 : outputVector(chan_idx) -= (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1373 : }
1374 : }
1375 : */
1376 :
1377 : /*
1378 : if (debug_p)
1379 : {
1380 : LogIO logger;
1381 : logger << "freqPows_p =" << freqPows_p << LogIO::POST;
1382 : logger << "fit order = " << fitOrder_p << LogIO::POST;
1383 : logger << "coeff =" << coeff << LogIO::POST;
1384 : logger << "inputFlags =" << inputFlags << LogIO::POST;
1385 : logger << "inputWeights =" << inputWeights << LogIO::POST;
1386 : logger << "inputVector =" << inputVector << LogIO::POST;
1387 : logger << "outputVector =" << outputVector << LogIO::POST;
1388 : }
1389 : */
1390 :
1391 0 : return chisq;
1392 0 : }
1393 :
1394 : //////////////////////////////////////////////////////////////////////////
1395 : // UVContEstimationDenoisingKernel class
1396 : //////////////////////////////////////////////////////////////////////////
1397 :
1398 : // -----------------------------------------------------------------------
1399 : //
1400 : // -----------------------------------------------------------------------
1401 0 : template<class T> UVContEstimationDenoisingKernel<T>::UVContEstimationDenoisingKernel( denoising::GslPolynomialModel<Double>* model,
1402 : size_t nIter,
1403 : Vector<bool> *lineFreeChannelMask):
1404 0 : UVContSubKernel<T>(model,lineFreeChannelMask)
1405 : {
1406 0 : fitter_p.resetModel(*model);
1407 0 : fitter_p.setNIter(nIter);
1408 0 : }
1409 :
1410 : // -----------------------------------------------------------------------
1411 : //
1412 : // -----------------------------------------------------------------------
1413 0 : template<class T> void UVContEstimationDenoisingKernel<T>::changeFitOrder(size_t order)
1414 : {
1415 0 : fitOrder_p = order;
1416 0 : fitter_p.resetNComponents(order+1);
1417 0 : return;
1418 : }
1419 :
1420 : // -----------------------------------------------------------------------
1421 : //
1422 : // -----------------------------------------------------------------------
1423 0 : template<class T> void UVContEstimationDenoisingKernel<T>::defaultKernel( Vector<T> &,
1424 : Vector<T> &outputVector)
1425 : {
1426 0 : outputVector = 0;
1427 0 : return;
1428 : }
1429 :
1430 : // -----------------------------------------------------------------------
1431 : //
1432 : // -----------------------------------------------------------------------
1433 0 : template<class T> Complex UVContEstimationDenoisingKernel<T>::kernelCore( Vector<T> &inputVector,
1434 : Vector<bool> &inputFlags,
1435 : Vector<Float> &inputWeights,
1436 : Vector<T> &outputVector)
1437 : {
1438 0 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1439 0 : std::vector<T> coeff;
1440 0 : Complex chisq;
1441 0 : tie(coeff, chisq) = fitter_p.calcFitCoeff(inputVector);
1442 0 : fitter_p.calcFitModel(outputVector);
1443 :
1444 : /*
1445 : fitter_p.setWeightsAndFlags(inputWeights,inputFlags);
1446 : Vector<T> coeff = fitter_p.calcFitCoeff(inputVector);
1447 :
1448 : // Fill output data
1449 : outputVector = coeff(0);
1450 : for (uInt order_idx = 1; order_idx <= fitOrder_p; order_idx++)
1451 : {
1452 : for (uInt chan_idx=0; chan_idx < outputVector.size(); chan_idx++)
1453 : {
1454 : outputVector(chan_idx) += (freqPows_p(order_idx,chan_idx))*coeff(order_idx);
1455 : }
1456 : }
1457 : */
1458 :
1459 0 : return chisq;
1460 0 : }
1461 :
1462 : } //# NAMESPACE VI - END
1463 :
1464 : } //# NAMESPACE CASA - END
1465 :
1466 :
|
