Line data Source code
1 : //# VisImagingWeight.cc: imaging weight calculation for a give buffer
2 : //# Copyright (C) 2011
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be adressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $Id$
28 :
29 :
30 : #include <casacore/casa/Utilities/CountedPtr.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/Vector.h>
33 : #include <casacore/casa/OS/Timer.h>
34 : #include <casacore/casa/Containers/Record.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 : #include <casacore/tables/Tables/ScaRecordColDesc.h>
37 : #include <casacore/tables/Tables/TableUtil.h>
38 : #include <components/ComponentModels/ComponentList.h>
39 : #include <casacore/ms/MSSel/MSSelection.h>
40 : #include <casacore/ms/MSSel/MSSelectionTools.h>
41 : #include <casacore/ms/MeasurementSets/MSSource.h>
42 : #include <casacore/ms/MSSel/MSSourceIndex.h>
43 : #include <casacore/ms/MeasurementSets/MSSourceColumns.h>
44 :
45 : #include <msvis/MSVis/VisBuffer.h>
46 : #include <msvis/MSVis/VisBuffer2.h>
47 : #include <msvis/MSVis/VisBuffer2Adapter.h>
48 : #include <msvis/MSVis/VisBufferUtil.h>
49 : #include <synthesis/TransformMachines/VisModelData.h>
50 : #include <synthesis/TransformMachines/FTMachine.h>
51 : #include <synthesis/TransformMachines/SimpleComponentFTMachine.h>
52 : #include <synthesis/TransformMachines/GridFT.h>
53 : #include <synthesis/TransformMachines/rGridFT.h>
54 : #include <synthesis/TransformMachines/MosaicFT.h>
55 : #include <synthesis/TransformMachines/MosaicFTNew.h>
56 : #include <synthesis/TransformMachines/WProjectFT.h>
57 : #include <synthesis/TransformMachines/MultiTermFT.h>
58 : #include <synthesis/TransformMachines/MultiTermFTNew.h>
59 : #include <synthesis/TransformMachines/SetJyGridFT.h>
60 :
61 : namespace {
62 :
63 1 : casa::VisModelDataI * createVisModelData (){
64 1 : return new casa::VisModelData ();
65 : }
66 : // FIXME initializeVisModelDataFactory defined but not used triggers compiler warning
67 : bool initializeVisModelDataFactory = casa::VisModelDataI::setFactory (createVisModelData);
68 :
69 : }
70 :
71 :
72 : using namespace casacore;
73 : namespace casa { //# NAMESPACE CASA - BEGIN
74 :
75 :
76 1 : VisModelData::VisModelData(): clholder_p(0), ftholder_p(0), flatholder_p(0), version_p(-1){
77 :
78 1 : cft_p=new SimpleComponentFTMachine();
79 1 : }
80 :
81 0 : VisModelData::~VisModelData(){
82 :
83 :
84 0 : }
85 :
86 : VisModelDataI *
87 0 : VisModelData::clone ()
88 : {
89 0 : return new VisModelData (* this);
90 : }
91 :
92 0 : void VisModelData::listModel(const MeasurementSet& thems){
93 :
94 : //Table newTab(thems);
95 :
96 0 : MSColumns msc(thems);
97 0 : Vector<String> fldnames=msc.field().name().getColumn();
98 0 : Vector<Int> fields=msc.fieldId().getColumn();
99 0 : const Sort::Order order=Sort::Ascending;
100 0 : const Int option=Sort::HeapSort | Sort::NoDuplicates;
101 0 : Int nfields=GenSort<Int>::sort (fields, order, option);
102 :
103 0 : LogIO logio;
104 0 : if (nfields>0) {
105 :
106 : logio << "MS Header field records:"
107 0 : << LogIO::POST;
108 :
109 0 : Int nlis(0);
110 0 : for (Int j=0; j< 2; ++j){
111 0 : const Table* thetab=&thems;
112 0 : if (j==1)
113 0 : thetab=&(thems.source());
114 0 : for (Int k=0; k< nfields; ++k){
115 0 : if(thetab->keywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
116 : {
117 0 : String elkey=thetab->keywordSet().asString("definedmodel_field_"+String::toString(fields[k]));
118 0 : if(thetab->keywordSet().isDefined(elkey))
119 0 : logio << " " << fldnames[fields[k]] << " (id = " << fields[k] << ")" << LogIO::POST;
120 0 : ++nlis;
121 : }
122 : }
123 : }
124 0 : if (nlis==0)
125 0 : logio << " None." << LogIO::POST;
126 : }
127 :
128 0 : }
129 :
130 0 : void VisModelData::deleteDiskImage(MeasurementSet& theMS, const String& theKey){
131 0 : TableRecord theRec;
132 0 : if(!VisModelData::getModelRecord(theKey, theRec, theMS))
133 0 : return;
134 0 : if(theRec.isDefined("numft")){
135 0 : Int numft=0;
136 0 : numft=theRec.asInt("numft");
137 0 : for (Int k=0; k < numft; ++k){
138 0 : String ftname=String("ft_")+String::toString(k);
139 0 : const Record& ftrec=theRec.asRecord(ftname);
140 0 : if(ftrec.asRecord("container").isDefined("diskimage")){
141 0 : String diskim=ftrec.asRecord("container").asString("diskimage");
142 0 : if(TableUtil::canDeleteTable(diskim))
143 0 : TableUtil::deleteTable(diskim);
144 : }
145 :
146 : }
147 : }
148 :
149 :
150 : }
151 0 : void VisModelData::removeRecordByKey(MeasurementSet& theMS, const String& theKey){
152 :
153 :
154 0 : deleteDiskImage(theMS, theKey);
155 0 : if(Table::isReadable(theMS.sourceTableName()) &&theMS.source().nrow() > 0 ){
156 0 : if(theMS.source().keywordSet().isDefined(theKey)){
157 0 : Int rowid=theMS.source().keywordSet().asInt(theKey);
158 0 : TableRecord elrec;
159 : //Replace the model with an empty record
160 0 : MSSourceColumns srcCol(theMS.source());
161 0 : srcCol.sourceModel().put(rowid, elrec);
162 0 : theMS.source().rwKeywordSet().removeField(theKey);
163 : }
164 : }
165 : //Remove from Main table if it is there
166 0 : if(theMS.rwKeywordSet().isDefined(theKey))
167 0 : theMS.rwKeywordSet().removeField(theKey);
168 0 : }
169 :
170 0 : void VisModelData::clearModel(const MeasurementSet& thems){
171 : // Table newTab(thems);
172 0 : MeasurementSet& newTab=const_cast<MeasurementSet& >(thems);
173 0 : if(!newTab.isWritable())
174 0 : return;
175 0 : auto part_block = newTab.getPartNames(true);
176 0 : Vector<String> theParts(part_block.begin( ),part_block.end( ));
177 0 : if(theParts.nelements() > 1){
178 0 : for (uInt k=0; k < theParts.nelements(); ++k){
179 0 : MeasurementSet subms(theParts[k], newTab.lockOptions(), Table::Update);
180 0 : clearModel(subms);
181 : }
182 0 : return;
183 : }
184 0 : Bool alreadyLocked=newTab.hasLock(True);
185 0 : if(!alreadyLocked)
186 0 : newTab.lock(True);
187 0 : if(Table::isReadable(newTab.sourceTableName())){
188 0 : newTab.source().lock(True);
189 : }
190 0 : LogIO logio;
191 : logio << "Clearing all model records in MS header."
192 0 : << LogIO::POST;
193 :
194 0 : MSColumns msc(thems);
195 0 : Vector<Int> fields=msc.fieldId().getColumn();
196 0 : Vector<String> fldnames=msc.field().name().getColumn();
197 0 : const Sort::Order order=Sort::Ascending;
198 0 : const Int option=Sort::HeapSort | Sort::NoDuplicates;
199 0 : Int nfields=GenSort<Int>::sort (fields, order, option);
200 0 : for (Int k=0; k< nfields; ++k){
201 0 : String elkey;
202 : Int srow;
203 : //if(newTab.rwKeywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
204 0 : if(isModelDefined(fields[k], thems, elkey, srow))
205 : {
206 0 : logio << " " << fldnames[fields[k]] << " (id = " << fields[k] << ") deleted." << LogIO::POST;
207 : //Remove from Source table
208 0 : removeRecordByKey(newTab, elkey);
209 0 : if(srow > -1){
210 0 : if(thems.source().keywordSet().isDefined("definedmodel_field_"+String::toString(fields[k]))){
211 0 : newTab.source().rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
212 : }
213 : }
214 : //Remove from Main table if it is there
215 0 : if(newTab.rwKeywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
216 0 : newTab.rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
217 : }
218 : else{
219 : ///remove keys that point to nowehere
220 : //cerr << "The key " << elkey << endl;
221 0 : if(!newTab.source().isNull() && newTab.source().rwKeywordSet().isDefined(elkey))
222 0 : newTab.source().rwKeywordSet().removeField(elkey);
223 0 : if(newTab.rwKeywordSet().isDefined(elkey))
224 0 : newTab.rwKeywordSet().removeField(elkey);
225 0 : if( !newTab.source().isNull() && newTab.source().rwKeywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
226 0 : newTab.source().rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
227 0 : if(newTab.rwKeywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
228 0 : newTab.rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
229 : }
230 : }
231 0 : if(!alreadyLocked)
232 0 : newTab.unlock();
233 0 : if(Table::isReadable(newTab.sourceTableName())){
234 0 : newTab.source().unlock();
235 : }
236 : ////Cleaning out orphaned image disk models
237 0 : String srctable=thems.source().isNull() ? "" : thems.source().tableName();
238 0 : Vector<String> possibleFT(2);
239 0 : possibleFT(0)=srctable+"/FT_MODEL*";
240 0 : possibleFT(1)=theParts[0]+"/FT_MODEL*";
241 :
242 0 : Vector<String> ftmods=Directory::shellExpand(possibleFT);
243 0 : for (uInt kk=0; kk < ftmods.nelements(); ++kk){
244 0 : if(TableUtil::canDeleteTable(ftmods[kk])){
245 0 : TableUtil::deleteTable(ftmods[kk]);
246 : }
247 : }
248 :
249 : }
250 0 : void VisModelData::clearModel(const MeasurementSet& thems, const String field, const String specwindows){
251 0 : MeasurementSet& newTab=const_cast<MeasurementSet& >(thems);
252 0 : auto part_block = newTab.getPartNames(true);
253 0 : Vector<String> theParts(part_block.begin( ),part_block.end( ));
254 0 : if(theParts.nelements() > 1){
255 0 : for (uInt k =0; k < theParts.nelements(); ++k){
256 0 : MeasurementSet subms(theParts[k], newTab.lockOptions(), Table::Update);
257 0 : clearModel(subms, field, specwindows);
258 : }
259 0 : return;
260 : }
261 0 : if(!newTab.isWritable())
262 0 : return;
263 :
264 0 : Bool alreadyLocked=newTab.hasLock(True);
265 0 : if(!alreadyLocked)
266 0 : newTab.lock(True);
267 0 : if(Table::isReadable(newTab.sourceTableName())){
268 0 : newTab.source().lock(True);
269 : }
270 :
271 0 : MSColumns msc(thems);
272 0 : Vector<String> fldnames=msc.field().name().getColumn();
273 0 : Int nfields=0;
274 0 : Vector<Int> fields(0);
275 : {
276 : // Parse field specification
277 0 : MSSelection mss;
278 0 : mss.setFieldExpr(field);
279 0 : TableExprNode ten=mss.toTableExprNode(&thems);
280 0 : fields=mss.getFieldList();
281 0 : nfields=fields.nelements();
282 : }
283 :
284 :
285 0 : if (nfields==0)
286 : // Call the method that deletes them all
287 0 : VisModelData::clearModel(thems);
288 : else {
289 :
290 : //Now we have the two cases the whole field or specific spws
291 : // only delete the specified ones
292 0 : Vector<Int> spws(0);
293 0 : Int nspws=0;
294 : {
295 : // Parse field specification
296 0 : MSSelection mss;
297 0 : mss.setFieldExpr(field);
298 0 : mss.setSpwExpr(specwindows);
299 0 : TableExprNode ten=mss.toTableExprNode(&thems);
300 0 : spws=mss.getSpwList();
301 0 : nspws=spws.nelements();
302 : }
303 :
304 :
305 0 : LogIO logio;
306 : logio << "Clearing model records in MS header for selected fields."
307 0 : << LogIO::POST;
308 :
309 0 : for (Int k=0; k< nfields; ++k){
310 0 : String elkey;
311 : Int srow;
312 0 : if(isModelDefined(fields[k], newTab, elkey, srow))
313 :
314 : {
315 :
316 0 : if(nspws==0){
317 0 : logio << " " << fldnames[fields[k]] << " (id = " << fields[k] << ") deleted." << LogIO::POST;
318 0 : removeRecordByKey(newTab, elkey);
319 0 : if(srow > -1 ){
320 0 : if(newTab.source().keywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
321 0 : newTab.source().rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
322 : }
323 0 : if(newTab.rwKeywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
324 0 : newTab.rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
325 : }
326 : else{
327 : //if(newTab.rwKeywordSet().isDefined(elkey)){
328 :
329 0 : TableRecord conteneur;
330 : //cerr << "Just before getModelRecord " << elkey << " fields " << fields << endl;
331 0 : getModelRecord(elkey, conteneur, newTab);
332 : //cerr << "keys " ;
333 : //for (uInt k=0; k< conteneur.nfields() ; ++k)
334 : // cerr << conteneur. name(k) ;
335 : //cerr << endl;
336 : //=newTab.rwKeywordSet().asRecord(elkey);
337 0 : if(removeSpw(conteneur, spws, fields)){
338 0 : putRecordByKey(newTab, elkey, conteneur, srow);
339 : }
340 : else{
341 0 : removeRecordByKey(newTab, elkey);
342 0 : if(newTab.source().keywordSet().isDefined("definedmodel_field_"+String::toString(fields[k])))
343 0 : newTab.source().rwKeywordSet().removeField("definedmodel_field_"+String::toString(fields[k]));
344 : }
345 : /*
346 : Vector<Int> defspws=conteneur.asArrayInt("spws");
347 : Vector<Int> newdefspw(defspws.nelements(), -1);
348 : Int counter=0;
349 : for(uInt k=0; k < defspws.nelements(); ++k){
350 : for (Int j=0; j < nspws; ++j){
351 : if(defspws[k] != spws[j]){
352 : newdefspw[counter]=defspws[k];
353 : ++counter;
354 : }
355 : }
356 : }
357 : if(counter==0){
358 : //Now we have to remove this ftm or cft
359 : newTab.rwKeywordSet().removeField(elkey);
360 : }
361 : else{
362 : conteneur.define("spws", newdefspw);
363 : updatespw(conteneur,
364 : }
365 : */
366 : }
367 :
368 :
369 :
370 : }
371 : else
372 0 : logio << " " << fldnames[fields[k]] << " (id = " << fields[k] << ") not found." << LogIO::POST;
373 : }
374 :
375 : }
376 0 : if(!alreadyLocked)
377 0 : newTab.unlock();
378 0 : if(Table::isReadable(newTab.sourceTableName())){
379 0 : newTab.source().unlock();
380 : }
381 :
382 : }
383 :
384 0 : Bool VisModelData::removeSpw(TableRecord& therec, const Vector<Int>& spws, const Vector<Int>& fields){
385 0 : Int numft=0;
386 0 : Vector<String> numtype(2);
387 0 : Vector<String> keyval(2);
388 0 : numtype[0]=String("numft");
389 0 : numtype[1]=String("numcl");
390 0 : keyval[0]="ft_";
391 0 : keyval[1]="cl_";
392 :
393 0 : for (Int j=0; j< 2; ++j){
394 0 : if(therec.isDefined(numtype[j])){
395 0 : numft=therec.asInt(numtype[j]);
396 0 : Vector<String> toberemoved(numft);
397 0 : Int numrem=0;
398 0 : for (Int k=0; k < numft; ++k){
399 0 : RecordInterface& ftrec=therec.asrwRecord(keyval[j]+String::toString(k));
400 0 : if(!ftrec.isDefined("version")){
401 0 : Bool hasField=true;
402 0 : if(fields.nelements() >0){
403 0 : hasField=false;
404 0 : Vector<Int> fieldsinrec;
405 0 : ftrec.get("fields", fieldsinrec);
406 0 : if(anyEQ(fieldsinrec, fields))
407 0 : hasField=true;
408 : }
409 :
410 0 : if(hasField && !removeSpwFromMachineRec(ftrec, spws)){
411 0 : toberemoved[numrem]=String(keyval[j]+String::toString(k));
412 0 : ++numrem;
413 : }
414 : }
415 : else{
416 : //Version 2
417 0 : const Matrix<Int>& ftcombind=ftrec.asArrayInt("indexcombination");
418 0 : Matrix<Int> newComb;
419 0 : for(uInt row=0; row < ftcombind.nrow(); ++row){
420 0 : Vector<Int> rowcomb=ftcombind.row(row);
421 0 : Bool rowTobeRemoved=False;
422 0 : for (uInt ff=0; ff < fields.nelements() ; ++ff){
423 0 : for (uInt ss=0; ss < spws.nelements(); ++ss){
424 0 : if(rowcomb[0]==fields[ff] && rowcomb[1] == spws[ss] )
425 0 : rowTobeRemoved=True;
426 :
427 :
428 : }
429 : }
430 0 : if(!rowTobeRemoved){
431 0 : newComb.resize(newComb.nrow()+1,4, True);
432 0 : newComb.row(newComb.nrow()-1)=rowcomb;
433 : }
434 : }
435 0 : if(newComb.nrow()==0){
436 0 : toberemoved[numrem]=String(keyval[j]+String::toString(k));
437 0 : ++numrem;
438 : }
439 : else{
440 0 : ftrec.define("indexcombination", newComb);
441 : }
442 :
443 : }
444 : }
445 0 : if(numrem > 0){
446 0 : for (Int k=0; k < numrem; ++k)
447 0 : removeFTFromRec(therec, toberemoved[k], k==(numrem-1));
448 : }
449 : }
450 : }
451 0 : numft=0; Int numcl=0;
452 0 : if(therec.isDefined("numft")) numft=therec.asInt("numft");
453 0 : if(therec.isDefined("numcl")) numcl=therec.asInt("numcl");
454 0 : if (numft==0 && numcl==0)
455 0 : return false;
456 0 : return true;
457 : }
458 :
459 0 : Bool VisModelData::removeFTFromRec(TableRecord& therec, const String& keyval, const Bool relabel){
460 :
461 :
462 0 : String *splitkey=new String[2];
463 0 : Int nsep=split(keyval, splitkey, 2, String("_"));
464 0 : if (nsep <1 || !therec.isDefined(keyval))
465 0 : return false;
466 0 : String eltype=splitkey[0];
467 : //Int modInd=String::toInt(splitkey[1]);
468 0 : Int numcomp= (eltype==String("ft")) ? therec.asInt("numft"): therec.asInt("numcl");
469 0 : therec.removeField(keyval);
470 :
471 0 : numcomp=numcomp-1;
472 0 : if(eltype=="ft")
473 0 : therec.define("numft", numcomp);
474 : else
475 0 : therec.define("numcl", numcomp);
476 0 : if(relabel){
477 :
478 :
479 0 : eltype=eltype+String("_");
480 0 : Int id=0;
481 0 : for(uInt k=0; k < therec.nfields(); ++k){
482 0 : if(therec.name(k).contains(eltype)){
483 0 : therec.renameField(eltype+String::toString(id), k);
484 0 : ++id;
485 : }
486 : }
487 : }
488 :
489 0 : delete [] splitkey;
490 0 : return true;
491 : }
492 :
493 0 : Bool VisModelData::removeSpwFromMachineRec(RecordInterface& ftclrec, const Vector<Int>& spws){
494 0 : Vector<Int> defspws=ftclrec.asArrayInt("spws");
495 0 : Vector<Int> newdefspw(defspws.nelements(), -1);
496 0 : Int counter=0;
497 0 : for(uInt k=0; k < defspws.nelements(); ++k){
498 0 : for (uInt j=0; j < spws.nelements(); ++j){
499 0 : if(defspws[k] == spws[j]){
500 0 : defspws[k]=-1;
501 0 : ++counter;
502 : }
503 : }
504 : }
505 0 : if(defspws.nelements() == uInt(counter)){
506 : //Now we have to remove this ft or cl model
507 :
508 0 : return false;
509 :
510 : }
511 0 : newdefspw.resize(defspws.nelements()-counter);
512 0 : counter=0;
513 0 : for (uInt k=0; k < defspws.nelements(); ++k){
514 0 : if(defspws[k] != -1){
515 0 : newdefspw[counter]=defspws[k];
516 0 : ++counter;
517 : }
518 : }
519 :
520 0 : ftclrec.define("spws", newdefspw);
521 0 : return true;
522 : }
523 :
524 0 : Bool VisModelData::addToRec(TableRecord& therec, const Vector<Int>& spws){
525 :
526 0 : Int numft=0;
527 0 : Int numcl=0;
528 0 : Vector<Bool> hasSpw(spws.nelements(), false);
529 0 : if(therec.isDefined("numft")){
530 0 : numft=therec.asInt("numft");
531 0 : Vector<Int> ft_toremove(numft, 0);
532 0 : for(Int k=0; k < numft; ++k){
533 0 : const Record& ftrec=therec.asRecord("ft_"+String::toString(k));
534 0 : const Vector<Int>& ftspws=ftrec.asArrayInt("spws");
535 0 : for (uInt i=0; i<spws.nelements(); ++i){
536 0 : for (uInt j=0; j<ftspws.nelements(); ++j){
537 0 : if(spws[i]==ftspws[j]){
538 0 : hasSpw[i]=true;
539 0 : ft_toremove[k]=1;
540 : }
541 : }
542 : }
543 : }
544 0 : if(sum(ft_toremove) >0){
545 0 : for(Int k=0; k < numft; ++k){
546 0 : if(ft_toremove[k]==1)
547 0 : therec.removeField("ft_"+String::toString(k));
548 : }
549 0 : numft=numft-sum(ft_toremove);
550 0 : therec.define("numft", numft);
551 0 : Int id=0;
552 0 : for(uInt k=0; k < therec.nfields(); ++k){
553 0 : if(therec.name(k).contains("ft_")){
554 0 : therec.renameField("ft_"+String::toString(id), k);
555 0 : ++id;
556 : }
557 : }
558 : }
559 : }
560 0 : if(therec.isDefined("numcl")){
561 0 : numcl=therec.asInt("numcl");
562 0 : Vector<Int> cl_toremove(numcl, 0);
563 0 : for(Int k=0; k < numcl; ++k){
564 0 : const Record& clrec=therec.asRecord("cl_"+String::toString(k));
565 0 : const Vector<Int>& clspws=clrec.asArrayInt("spws");
566 0 : for (uInt i=0; i<spws.nelements(); ++i){
567 0 : for (uInt j=0; j<clspws.nelements(); ++j){
568 0 : if(spws[i]==clspws[j]){
569 0 : hasSpw[i]=true;
570 0 : cl_toremove[k]=1;
571 : }
572 : }
573 : }
574 : }
575 0 : if(sum(cl_toremove) >0){
576 0 : for(Int k=0; k < numcl; ++k){
577 0 : if(cl_toremove[k]==1)
578 0 : therec.removeField("cl_"+String::toString(k));
579 : }
580 0 : numcl=numcl-sum(cl_toremove);
581 0 : therec.define("numcl", numcl);
582 0 : Int id=0;
583 0 : for(uInt k=0; k < therec.nfields(); ++k){
584 0 : if(therec.name(k).contains("cl_")){
585 0 : therec.renameField("cl_"+String::toString(id), k);
586 0 : ++id;
587 : }
588 : }
589 : }
590 : }
591 0 : return (!allTrue(hasSpw) || ((numft+numcl)>0));
592 : }
593 :
594 0 : Bool VisModelData::addToRec(TableRecord& therec, const Matrix<Int>& combind){
595 0 : if(therec.nfields() >0 && ((!therec.isDefined("version")) || (therec.asInt("version") !=2)))
596 0 : throw(AipsError("cannot combine with older version virtual model column;\n please clear it before proceeding"));
597 0 : Int numft=0;
598 0 : Int numcl=0;
599 0 : Vector<Bool> hasComb(combind.nrow(), false);
600 0 : if(therec.isDefined("numft")){
601 0 : numft=therec.asInt("numft");
602 0 : Vector<Int> ft_toremove(numft, 0);
603 0 : for(Int k=0; k < numft; ++k){
604 0 : const Record& ftrec=therec.asRecord("ft_"+String::toString(k));
605 0 : const Matrix<Int>& ftcombind=ftrec.asArrayInt("indexcombination");
606 0 : for (uInt i=0; i<combind.nrow(); ++i){
607 0 : for (uInt j=0; j<ftcombind.nrow(); ++j){
608 0 : if(allEQ(combind.row(i), ftcombind.row(j))){
609 : //We could upgrade this to just delete the intent from the ft
610 0 : hasComb[i]=true;
611 0 : ft_toremove[k]=1;
612 : }
613 : }
614 : }
615 : }
616 0 : if(sum(ft_toremove) >0){
617 0 : for(Int k=0; k < numft; ++k){
618 0 : if(ft_toremove[k]==1)
619 0 : therec.removeField("ft_"+String::toString(k));
620 : }
621 0 : numft=numft-sum(ft_toremove);
622 0 : therec.define("numft", numft);
623 0 : Int id=0;
624 0 : for(uInt k=0; k < therec.nfields(); ++k){
625 0 : if(therec.name(k).contains("ft_")){
626 0 : therec.renameField("ft_"+String::toString(id), k);
627 0 : ++id;
628 : }
629 : }
630 : }
631 : }
632 0 : if(therec.isDefined("numcl")){
633 0 : numcl=therec.asInt("numcl");
634 0 : Vector<Int> cl_toremove(numcl, 0);
635 0 : for(Int k=0; k < numcl; ++k){
636 0 : const Record& clrec=therec.asRecord("cl_"+String::toString(k));
637 0 : const Matrix<Int>& clcombind=clrec.asArrayInt("indexcombination");
638 0 : for (uInt i=0; i<combind.nrow(); ++i){
639 0 : for (uInt j=0; j<clcombind.nrow(); ++j){
640 0 : if(allEQ(combind.row(i),clcombind.row(j))){
641 0 : hasComb[i]=true;
642 0 : cl_toremove[k]=1;
643 : }
644 : }
645 : }
646 : }
647 0 : if(sum(cl_toremove) >0){
648 0 : for(Int k=0; k < numcl; ++k){
649 0 : if(cl_toremove[k]==1)
650 0 : therec.removeField("cl_"+String::toString(k));
651 : }
652 0 : numcl=numcl-sum(cl_toremove);
653 0 : therec.define("numcl", numcl);
654 0 : Int id=0;
655 0 : for(uInt k=0; k < therec.nfields(); ++k){
656 0 : if(therec.name(k).contains("cl_")){
657 0 : therec.renameField("cl_"+String::toString(id), k);
658 0 : ++id;
659 : }
660 : }
661 : }
662 : }
663 0 : return (!allTrue(hasComb) || ((numft+numcl)>0));
664 : }
665 :
666 :
667 200 : Bool VisModelData::isModelDefined(const Int fieldId, const MeasurementSet& thems, String& thekey, Int& sourceRow){
668 200 : sourceRow=-1;
669 600 : String modelkey=String("definedmodel_field_")+String::toString(fieldId);
670 200 : thekey="";
671 200 : if(Table::isReadable(thems.sourceTableName()) && thems.source().nrow() > 0 && (thems.source().keywordSet().isDefined(modelkey))){
672 : {
673 0 : thekey=thems.source().keywordSet().asString(modelkey);
674 0 : if(thems.source().keywordSet().isDefined(thekey))
675 0 : sourceRow=thems.source().keywordSet().asInt(thekey);
676 : }
677 : }else{
678 200 : if(thems.keywordSet().isDefined(modelkey)){
679 0 : thekey=thems.keywordSet().asString(modelkey);
680 : }
681 : }
682 200 : if(thekey != "" )
683 0 : return VisModelData::isModelDefined(thekey, thems);
684 200 : return false;
685 :
686 : }
687 :
688 0 : Bool VisModelData::isModelDefined(const String& elkey, const MeasurementSet& thems){
689 : //Let's try the Source table
690 0 : if(Table::isReadable(thems.sourceTableName()) &&thems.source().nrow() > 0 ){
691 0 : if(thems.source().keywordSet().isDefined(elkey))
692 0 : return true;
693 : }
694 : //Let's try the Main table
695 0 : if(thems.keywordSet().isDefined(elkey))
696 0 : return true;
697 0 : return false;
698 : }
699 :
700 0 : Bool VisModelData::getModelRecord(const String& theKey, TableRecord& theRec, const MeasurementSet& theMs){
701 : //Let's try the Source table
702 0 : if(Table::isReadable(theMs.sourceTableName()) &&theMs.source().nrow() > 0 && (theMs.source().keywordSet().isDefined(theKey))){
703 : {
704 : //Get the row for the model
705 0 : Int row=theMs.source().keywordSet().asInt(theKey);
706 : //MSSourceColumns srcCol(theMs.source());
707 :
708 0 : ScalarColumn<TableRecord> scol(theMs.source(), "SOURCE_MODEL");
709 0 : scol.get(row, theRec);
710 : }
711 0 : return true;
712 : }
713 : //Let's try the Main table
714 0 : if(theMs.keywordSet().isDefined(theKey)){
715 0 : theRec=theMs.keywordSet().asRecord(theKey);
716 0 : return true;
717 : }
718 0 : return false;
719 :
720 :
721 : }
722 :
723 0 : void VisModelData::putRecordByKey(MeasurementSet& theMS, const String& theKey, const TableRecord& theRec, const Int sourceRowNum){
724 : //Prefer the Source table first
725 0 : if( (sourceRowNum> -1) && Table::isReadable(theMS.sourceTableName()) && Int(theMS.source().nrow()) > sourceRowNum ){
726 0 : MSSource& mss=theMS.source();
727 0 : if(!mss.isColumn(MSSource::SOURCE_MODEL) ){
728 0 : mss.addColumn(ScalarRecordColumnDesc("SOURCE_MODEL"), true);
729 : }
730 0 : if(mss.rwKeywordSet().isDefined(theKey))
731 0 : mss.rwKeywordSet().removeField(theKey);
732 0 : mss.rwKeywordSet().define(theKey, sourceRowNum);
733 0 : MSSourceColumns srcCol(mss);
734 0 : srcCol.sourceModel().put(sourceRowNum, theRec);
735 0 : return;
736 : }
737 : //Oh well no source table so will add it to the main table
738 0 : theMS.rwKeywordSet().defineRecord(theKey, theRec);
739 :
740 :
741 : }
742 :
743 0 : Bool VisModelData::putModelRecord(const Vector<Int>& fieldIds, TableRecord& theRec, MeasurementSet& theMS){
744 0 : auto part_block = theMS.getPartNames(true);
745 0 : Vector<String> theParts(part_block.begin( ),part_block.end( ));
746 0 : if(theParts.nelements() > 1){
747 0 : Bool retval=true;
748 0 : for (uInt k =0; k < theParts.nelements(); ++k){
749 0 : MeasurementSet subms(theParts[k], theMS.lockOptions(), Table::Update);
750 0 : retval= retval && putModelRecord(fieldIds, theRec, subms);
751 : }
752 0 : return retval;
753 : }
754 :
755 0 : String elkey="model";
756 0 : for (uInt k=0; k < fieldIds.nelements(); ++k){
757 0 : elkey=elkey+"_"+String::toString(fieldIds[k]);
758 : }
759 0 : if(theMS.rwKeywordSet().isDefined(elkey))
760 0 : theMS.rwKeywordSet().removeField(elkey);
761 0 : Int row=-1;
762 : //Prefer the Source table first
763 0 : MSFieldColumns fCol(theMS.field());
764 0 : if(Table::isReadable(theMS.sourceTableName()) && (theMS.source().nrow() > 0) && (!fCol.sourceId().isNull()) && (fCol.sourceId().get(fieldIds[0]) != -1) ){
765 : //
766 0 : row=0;
767 0 : MSSource& mss=theMS.source();
768 :
769 0 : Int sid=fCol.sourceId().get(fieldIds[0]);
770 0 : auto rows=MSSourceIndex(mss).getRowNumbersOfSourceID(sid);
771 0 : if(rows.nelements() > 0)
772 0 : row=rows[0];
773 : else{
774 0 : LogIO logio;
775 0 : logio << "Invalid Source_id "+String::toString(sid)+" found in FIELD table\n"
776 : <<"Model is being written at Source ID 0 position which will be erased\n"
777 : << "Fix the FIELD table before proceeding "
778 0 : << LogIO::WARN << LogIO::POST;
779 :
780 : }
781 :
782 :
783 :
784 0 : putRecordByKey(theMS, elkey, theRec, row);
785 0 : for (uInt k=0; k < fieldIds.nelements(); ++k){
786 0 : mss.rwKeywordSet().define("definedmodel_field_"+String::toString(fieldIds[k]), elkey);
787 : }
788 0 : return true;
789 :
790 : }
791 : //Oh well no source table so will add it to the main table
792 0 : putRecordByKey(theMS, elkey, theRec, -1);
793 0 : for (uInt k=0; k < fieldIds.nelements(); ++k){
794 0 : theMS.rwKeywordSet().define("definedmodel_field_"+String::toString(fieldIds[k]), elkey);
795 : }
796 :
797 0 : return true;
798 : }
799 :
800 :
801 0 : void VisModelData::putModel(const MeasurementSet& thems, const RecordInterface& rec, const Vector<Int>& validfieldids, const Vector<Int>& spws, const Vector<Int>& starts, const Vector<Int>& nchan, const Vector<Int>& incr, Bool iscomponentlist, Bool incremental){
802 :
803 0 : LogIO logio;
804 :
805 : try{
806 : //A field can have multiple FTmachines and ComponentList associated with it
807 : //For example having many flanking images for the model
808 : //For componentlist it may have multiple componentlist ...for different spw
809 : //Timer tim;
810 : //tim.mark();
811 0 : Int counter=0;
812 0 : Record modrec;
813 0 : modrec.define("fields", validfieldids);
814 0 : modrec.define("spws", spws);
815 0 : modrec.define("start", starts);
816 0 : modrec.define("nchan", nchan);
817 0 : modrec.define("incr", incr);
818 0 : modrec.defineRecord("container", rec);
819 0 : String elkey="model";
820 0 : for (uInt k=0; k < validfieldids.nelements(); ++k){
821 0 : elkey=elkey+"_"+String::toString(validfieldids[k]);
822 : }
823 0 : TableRecord outRec;
824 0 : Bool addtorec=false;
825 0 : MeasurementSet& newTab=const_cast<MeasurementSet& >(thems);
826 : //cerr << elkey << " incr " << incremental << endl;
827 0 : if(isModelDefined(elkey, newTab)){
828 0 : getModelRecord(elkey, outRec, thems);
829 : //if incremental no need to check & remove what is in the record
830 0 : if(!incremental)
831 0 : addtorec=addToRec(outRec, spws);
832 : //cerr << "addToRec " << addtorec << endl;
833 : }
834 : else{
835 : ///////even if it is not defined some other field model might be sitting on that
836 : //////model key
837 0 : Int hasSourceRecord=firstSourceRowRecord(validfieldids[0], thems, outRec);
838 0 : if(hasSourceRecord > -1 && outRec.nfields() > 0)
839 0 : addtorec=true;
840 : //cerr << "has Source " << hasSourceRecord << " addToRec " << addtorec << endl;
841 : ////
842 : }
843 0 : incremental=incremental || addtorec;
844 0 : if(iscomponentlist){
845 0 : modrec.define("type", "componentlist");
846 0 : if(outRec.isDefined("numcl"))
847 0 : counter=incremental ? outRec.asInt("numcl") : 0;
848 :
849 : }
850 : else{
851 0 : modrec.define("type", "ftmachine");
852 0 : if(outRec.isDefined("numft"))
853 0 : counter=incremental ? outRec.asInt("numft") : 0;
854 : }
855 0 : iscomponentlist ? outRec.define("numcl", counter+1) : outRec.define("numft", counter+1);
856 :
857 : //for (uInt k=0; k < validfieldids.nelements(); ++k){
858 : // newTab.rwKeywordSet().define("definedmodel_field_"+String::toString(validfieldids[k]), elkey);
859 :
860 : // }
861 0 : iscomponentlist ? outRec.defineRecord("cl_"+String::toString(counter), modrec):
862 0 : outRec.defineRecord("ft_"+String::toString(counter), modrec);
863 : //////////////////;
864 : //for (uInt k=0; k < newTab.rwKeywordSet().nfields() ; ++k){
865 : // cerr << "keys " << k << " is " << newTab.rwKeywordSet().name(k) << " type " << newTab.rwKeywordSet().dataType(k) << endl;
866 : //}
867 : ////////////////////////
868 : //if image for a given key is on disk and not incrementing ...lets remove it
869 0 : if(!incremental)
870 0 : deleteDiskImage(newTab, elkey);
871 0 : putModelRecord(validfieldids, outRec, newTab);
872 : //if(newTab.rwKeywordSet().isDefined(elkey))
873 : // newTab.rwKeywordSet().removeField(elkey);
874 : //newTab.rwKeywordSet().defineRecord(elkey, outRec);
875 :
876 : // tim.show("Time taken to save record ");
877 :
878 : /*
879 : String subName=newTab.tableName()+"/"+elkey;
880 : if(Table::isReadable(subName) && !Table::canDeleteTable(subName, true))
881 : throw(AipsError("Cannot save model into MS"));
882 : else if ((Table::isReadable(subName) && Table::canDeleteTable(subName, true)))
883 : Table::deleteTable(subName, true);
884 : TableDesc td1 ("td1", TableDesc::New);
885 : //td1.addColumn (ArrayColumnDesc<Int> ("SPECTRAL_WINDOW_ID"));
886 : //td1.addColumn(ScalarColumnDesc<TableRecord>("MODEL"));
887 : SetupNewTable newtab1 (subName, td1, Table::New);
888 : Table tab1 (newtab1);
889 : newTab.rwKeywordSet().defineTable(elkey, tab1);
890 : tab1.rwKeywordSet().defineRecord(elkey, outRec);
891 : */
892 : //ArrayColumn<Int> spwCol(tab1, "SPECTRAL_WINDOW_ID");
893 : //ScalarColumn<TableRecord> modCol(tab1, "MODEL");
894 : //newTab.addRow(1,false);
895 : //spwCol.put(0,spws);
896 : //modCol.put(0,outRec);
897 : //newTab.flush();
898 : // MSSource& mss=newTab.source();
899 : // cerr << "has model_source " << mss.isColumn(MSSource::SOURCE_MODEL) << endl;
900 : // if(!mss.isColumn(MSSource::SOURCE_MODEL) ){
901 : //mss.addColumn(ScalarRecordColumnDesc("SOURCE_MODEL"), true);
902 : // }
903 : // MSSourceColumns srcCol(mss);
904 : // srcCol.sourceModel().put(0, outRec);
905 : }
906 0 : catch(...){
907 0 : logio << "Could not save virtual model data column due to an artificial virtual model size limit. \nYou may need to use the scratch column if you need model visibilities" << LogIO::WARN << LogIO::POST ;
908 :
909 : }
910 :
911 0 : }
912 :
913 0 : void VisModelData::putModel(const MeasurementSet& thems,const RecordInterface& rec, const Matrix<Int>& indexComb, const Matrix<Int>& chanSel,Bool iscomponentlist, Bool incremental){
914 0 : LogIO logio;
915 :
916 : try{
917 : //A field can have multiple FTmachines and ComponentList associated with it
918 : //For example having many flanking images for the model
919 : //For componentlist it may have multiple componentlist ...for different spw
920 : //Timer tim;
921 : //tim.mark();
922 0 : Int counter=0;
923 0 : Record modrec;
924 0 : modrec.define("version", Int(2));
925 0 : Vector<Int> validfieldids;
926 0 : validfieldids=indexComb.column(0);
927 0 : Int nfields=GenSort<Int>::sort (validfieldids, Sort::Ascending,
928 0 : Sort::QuickSort|Sort::NoDuplicates);
929 0 : validfieldids.resize(nfields, True);
930 0 : modrec.define("indexcombination", indexComb);
931 0 : modrec.define("channelselection", chanSel);
932 0 : modrec.defineRecord("container", rec);
933 0 : String elkey="model";
934 0 : for (uInt k=0; k < validfieldids.nelements(); ++k){
935 0 : elkey=elkey+"_"+String::toString(validfieldids[k]);
936 : }
937 0 : TableRecord outRec;
938 0 : Bool addtorec=false;
939 0 : MeasurementSet& newTab=const_cast<MeasurementSet& >(thems);
940 0 : newTab.lock(True);
941 0 : if(Table::isReadable(newTab.sourceTableName())){
942 0 : newTab.source().lock(True);
943 : }
944 : ////TESTOO
945 : //Int CPUID;
946 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
947 : //cerr <<"VISMOD-ftm1 ver2 " << CPUID << " ELKEY " << elkey << "has lock "<< newTab.hasLock(True) << endl;
948 : //
949 0 : if(isModelDefined(elkey, newTab)){
950 0 : getModelRecord(elkey, outRec, thems);
951 : //if incremental no need to check & remove what is in the record
952 0 : if(!incremental)
953 0 : addtorec=addToRec(outRec, indexComb);
954 : //cerr << "addToRec " << addtorec << endl;
955 : }
956 : else{
957 : ///////even if it is not defined some other field model might be sitting on that
958 : //////model key
959 0 : Int hasSourceRecord=firstSourceRowRecord(validfieldids[0], thems, outRec);
960 0 : if(hasSourceRecord > -1 && outRec.nfields() > 0)
961 0 : addtorec=true;
962 : //cerr << "has Source " << hasSourceRecord << " addToRec " << addtorec << endl;
963 : ////
964 : }
965 0 : if(outRec.nfields() >0 && ((!outRec.isDefined("version")) || (outRec.asInt("version") !=2)))
966 0 : throw(AipsError("cannot combine with older version virtual model column;\n please clear it before proceeding"));
967 0 : outRec.define("version", Int(2));
968 0 : incremental=incremental || addtorec;
969 0 : if(iscomponentlist){
970 0 : modrec.define("type", "componentlist");
971 0 : if(outRec.isDefined("numcl"))
972 0 : counter=incremental ? outRec.asInt("numcl") : 0;
973 :
974 : }
975 : else{
976 0 : modrec.define("type", "ftmachine");
977 0 : if(outRec.isDefined("numft"))
978 0 : counter=incremental ? outRec.asInt("numft") : 0;
979 : }
980 0 : iscomponentlist ? outRec.define("numcl", counter+1) : outRec.define("numft", counter+1);
981 :
982 : //for (uInt k=0; k < validfieldids.nelements(); ++k){
983 : // newTab.rwKeywordSet().define("definedmodel_field_"+String::toString(validfieldids[k]), elkey);
984 :
985 : // }
986 0 : iscomponentlist ? outRec.defineRecord("cl_"+String::toString(counter), modrec):
987 0 : outRec.defineRecord("ft_"+String::toString(counter), modrec);
988 : //////////////////;
989 : //for (uInt k=0; k < newTab.rwKeywordSet().nfields() ; ++k){
990 : // cerr << "keys " << k << " is " << newTab.rwKeywordSet().name(k) << " type " << newTab.rwKeywordSet().dataType(k) << endl;
991 : //}
992 : ////////////////////////
993 : //if image for a given key is on disk and not incrementing ...lets remove it
994 0 : if(!incremental)
995 0 : deleteDiskImage(newTab, elkey);
996 0 : putModelRecord(validfieldids, outRec, newTab);
997 0 : newTab.unlock();
998 0 : if(Table::isReadable(newTab.sourceTableName())){
999 0 : newTab.source().unlock();
1000 : }
1001 :
1002 : }
1003 0 : catch(...){
1004 0 : logio << "Could not save virtual model data for some reason \nYou may need clear the model and redo or use the scratch column if you need model visibilities" << LogIO::WARN << LogIO::POST ;
1005 0 : const_cast<MeasurementSet& >(thems).unlock();
1006 : }
1007 :
1008 0 : }
1009 0 : void VisModelData::addModel(const RecordInterface& rec, const Vector<Int>& /*msids*/, const vi::VisBuffer2& vb)
1010 : {
1011 :
1012 0 : Vector<Int> dum;
1013 0 : addModel(rec, dum, vi::VisBuffer2Adapter(&vb));
1014 :
1015 :
1016 0 : }
1017 0 : void VisModelData::modifyDiskImagePath(Record& rec, const VisBuffer& vb){
1018 0 : Record& ftmac= rec.rwSubRecord("container");
1019 0 : if(ftmac.isDefined("diskimage")){
1020 0 : String theDiskImage=ftmac.asString("diskimage");
1021 0 : if(File(theDiskImage).exists())
1022 0 : return;
1023 0 : String subPathname[30];
1024 0 : String sep = "/";
1025 0 : uInt nw = split(theDiskImage, subPathname, 20, sep);
1026 0 : String theposs=(subPathname[nw-1]);
1027 0 : String msname=vb.msName(false);
1028 0 : for (uInt i=nw-2 ; i>0; --i){
1029 0 : if(File(msname+"/"+theposs).exists()){
1030 0 : theDiskImage=msname+"/"+theposs;
1031 0 : ftmac.define("diskmage", theDiskImage);
1032 0 : return;
1033 : }
1034 0 : theposs=subPathname[i]+"/"+theposs;
1035 : }
1036 : }
1037 : }
1038 :
1039 0 : void VisModelData::addModel(const RecordInterface& rec, const Vector<Int>& /*msids*/, const VisBuffer& vb){
1040 :
1041 :
1042 :
1043 0 : Int indexft=-1;
1044 0 : if(rec.isDefined("numft")){
1045 0 : Int numft=rec.asInt("numft");
1046 0 : if(numft >0){
1047 0 : for(Int ftk=0; ftk < numft; ++ftk){
1048 0 : Record ftrec(rec.asRecord("ft_"+String::toString(ftk)));
1049 0 : modifyDiskImagePath(ftrec, vb);
1050 0 : Vector<Int>fields;
1051 0 : Vector<Int> spws;
1052 0 : ftrec.get("fields", fields);
1053 0 : ftrec.get("spws", spws);
1054 0 : if(anyEQ(spws, vb.spectralWindow())){
1055 0 : indexft=ftholder_p.nelements();
1056 0 : ftholder_p.resize(indexft+1, false, true);
1057 0 : ftholder_p[indexft].resize(1);
1058 0 : ftholder_p[indexft][0]=NEW_FT(ftrec.asRecord("container"));
1059 0 : if(!( ftholder_p[indexft][0]))
1060 0 : throw(AipsError("Unsupported FTMachine found in virtual MODEL_DATA column"));
1061 0 : ftholder_p[indexft][0]->initMaps(vb);
1062 :
1063 0 : for( uInt fi=0; fi < fields.nelements(); ++fi){
1064 0 : for(uInt spi=0; spi < spws.nelements(); ++spi){
1065 0 : Int indx=-1;
1066 0 : Int ftindx=-1;
1067 0 : if(hasModel(vb.msId(), fields[fi], spws[spi]) && (ftindex_p(spws[spi], fields[fi], vb.msId()) >= 0 )){
1068 :
1069 0 : indx=ftindex_p(spws[spi], fields[fi], vb.msId());
1070 0 : ftindx=ftholder_p[indx].nelements();
1071 0 : Bool alreadyAdded=false;
1072 0 : for (Int kk=1; kk < ftindx; ++kk){
1073 0 : alreadyAdded= alreadyAdded || (ftholder_p[indexft][0]==ftholder_p[indx][kk]);
1074 : }
1075 0 : if(!alreadyAdded){
1076 0 : ftholder_p[indx].resize(ftindx+1, true);
1077 0 : ftholder_p[indx][ftindx]=ftholder_p[indexft][0];
1078 : }
1079 : }
1080 : else{
1081 0 : ftindex_p(spws[spi], fields[fi], vb.msId())=indexft;
1082 : }
1083 : }
1084 : }
1085 : }
1086 : else{
1087 0 : if(hasModel(vb.msId(), vb.fieldId(), vb.spectralWindow()) < 0)
1088 0 : ftindex_p(vb.spectralWindow(), vb.fieldId(), vb.msId())=-2;
1089 : }
1090 :
1091 :
1092 : }
1093 : }
1094 : }
1095 0 : Int indexcl=-1;
1096 0 : if(rec.isDefined("numcl")){
1097 0 : Int numcl=rec.asInt("numcl");
1098 0 : if(numcl >0){
1099 0 : for(Int clk=0; clk < numcl; ++clk){
1100 0 : Vector<Int>fields;
1101 0 : Vector<Int> spws;
1102 0 : Record clrec(rec.asRecord("cl_"+String::toString(clk)));
1103 0 : clrec.get("fields", fields);
1104 0 : clrec.get("spws", spws);
1105 0 : if(anyEQ(spws, vb.spectralWindow())){
1106 0 : indexcl=clholder_p.nelements();
1107 0 : clholder_p.resize(indexcl+1, false, true);
1108 0 : clholder_p[indexcl].resize(1);
1109 0 : clholder_p[indexcl][0]=new ComponentList();
1110 0 : String err;
1111 0 : if(!((clholder_p[indexcl][0])->fromRecord(err, clrec.asRecord("container"))))
1112 0 : throw(AipsError("Component model failed to load for field "+String::toString(fields)));
1113 0 : for( uInt fi=0; fi < fields.nelements(); ++fi){
1114 0 : for(uInt spi=0; spi < spws.nelements(); ++spi){
1115 0 : Int indx=-1;
1116 0 : Int clindx=-1;
1117 0 : if(hasModel(vb.msId(), fields[fi], spws[spi]) && (clindex_p(spws[spi], fields[fi], vb.msId()) >= 0 )){
1118 0 : indx=clindex_p(spws[spi], fields[fi], vb.msId());
1119 0 : clindx=clholder_p[indx].nelements();
1120 0 : Bool alreadyAdded=false;
1121 0 : for (Int kk=1; kk < clindx; ++kk){
1122 0 : alreadyAdded= alreadyAdded || (clholder_p[indexcl][0]==clholder_p[indx][kk]);
1123 : }
1124 0 : if(!alreadyAdded){
1125 0 : clholder_p[indx].resize(clindx+1, true);
1126 0 : clholder_p[indx][clindx]=clholder_p[indexcl][0];
1127 : }
1128 : }
1129 : else{
1130 0 : clindex_p(spws[spi], fields[fi], vb.msId())=indexcl;
1131 : }
1132 : }
1133 : }
1134 : }
1135 : else{
1136 0 : if(hasModel(vb.msId(), vb.fieldId(), vb.spectralWindow()) < 0)
1137 0 : clindex_p(vb.spectralWindow(), vb.fieldId(), vb.msId())=-2;
1138 : }
1139 :
1140 : }
1141 : }
1142 : }
1143 :
1144 :
1145 0 : }
1146 :
1147 0 : FTMachine* VisModelData::NEW_FT(const Record& ftrec){
1148 0 : String name=ftrec.asString("name");
1149 : //cerr << "NEW_FT " << name << endl;
1150 0 : if(name=="GridFT")
1151 0 : return new GridFT(ftrec);
1152 0 : if(name=="rGridFT")
1153 0 : return new rGridFT(ftrec);
1154 0 : if(name=="WProjectFT")
1155 0 : return new WProjectFT(ftrec);
1156 0 : if(name=="MultiTermFT")
1157 0 : return new MultiTermFT(ftrec);
1158 0 : if(name=="MosaicFT")
1159 0 : return new MosaicFT(ftrec);
1160 0 : if(name=="MosaicFTNew")
1161 0 : return new MosaicFTNew(ftrec);
1162 0 : if(name=="SetJyGridFT")
1163 0 : return new SetJyGridFT(ftrec);
1164 0 : if(name=="MultiTermFTNew")
1165 0 : return new MultiTermFTNew(ftrec);
1166 : //When the following have constructors from Record they should be uncommented
1167 : // if(name=="AWProjectFT")
1168 : // return new AWProjectFT(ftrec);
1169 : //if(name=="AWProjectWBFT")
1170 : // return new AWProjectWBFT(ftrec);
1171 : //if(name=="MultiTermAWProjectWBFT")
1172 : // return new MultiTermAWProjectWBFT(ftrec);
1173 0 : return NULL;
1174 : }
1175 :
1176 0 : Int VisModelData::hasModel(Int msid, Int field, Int spw){
1177 :
1178 0 : IPosition oldcubeShape=ftindex_p.shape();
1179 0 : if(oldcubeShape(0) <(spw+1) || oldcubeShape(1) < (field+1) || oldcubeShape(2) < (msid+1)){
1180 0 : Cube<Int> newind(max((spw+1), oldcubeShape(0)), max((field+1),oldcubeShape(1)) , max((msid+1), oldcubeShape(2)));
1181 0 : newind.set(-1);
1182 0 : newind(IPosition(3, 0,0,0), (oldcubeShape-1))=ftindex_p;
1183 0 : ftindex_p.assign(newind);
1184 0 : newind.set(-1);
1185 0 : newind(IPosition(3, 0,0,0), (oldcubeShape-1))=clindex_p;
1186 0 : clindex_p.assign(newind);
1187 : }
1188 :
1189 0 : if( (clindex_p(spw, field, msid) + ftindex_p(spw, field, msid)) < -2)
1190 0 : return -2;
1191 0 : else if( (clindex_p(spw, field, msid) ==-1) && (ftindex_p(spw, field, msid) ==-1))
1192 0 : return -1;
1193 0 : return 1;
1194 :
1195 :
1196 : }
1197 :
1198 0 : void VisModelData::initializeToVis(){
1199 :
1200 :
1201 0 : }
1202 :
1203 0 : Bool VisModelData::getModelVis(vi::VisBuffer2& vb){
1204 0 : vi::VisBuffer2Adapter vba(&vb);
1205 0 : return getModelVis(vba);
1206 :
1207 : }
1208 0 : void VisModelData::getMatchingMachines(Vector<CountedPtr<FTMachine> >& ft, Vector<CountedPtr<ComponentList> >& cl, const VisBuffer& vb){
1209 0 : if(isVersion2()){
1210 0 : Vector<Vector<Int> > combindx;
1211 0 : getUniqueIndicesComb(vb, combindx);
1212 0 : if(combindx.nelements() != 1)
1213 0 : throw(AipsError("Cannot deal with multiple intent per visbuffer "));
1214 0 : std::vector<Int> indexInBuf=combindx[0].tovector();
1215 0 : indexInBuf.push_back(vb.msId());
1216 0 : auto itFTMap = ftindex2_p.find(indexInBuf);
1217 0 : if(itFTMap != ftindex2_p.end() ) {
1218 0 : if((itFTMap->second) < 0){
1219 0 : ft.resize(0);
1220 : }
1221 : else{
1222 0 : ft=ftholder_p[itFTMap->second];
1223 : }
1224 : }
1225 0 : auto itCLMap = clindex2_p.find(indexInBuf);
1226 0 : if(itCLMap != clindex2_p.end() ) {
1227 0 : if((itCLMap->second)< 0){
1228 0 : cl.resize(0);
1229 : }
1230 : else{
1231 0 : cl=clholder_p[itCLMap->second];
1232 : }
1233 : }
1234 : ///Now let's deal with a key that has not been visited before
1235 0 : if((itCLMap == clindex2_p.end()) && (itFTMap == ftindex2_p.end() )){
1236 0 : cerr << "no matching holder " << Vector<Int>(indexInBuf) << " num of cl " << clholder_p.nelements() << endl;
1237 0 : updateHolders(vb, indexInBuf);
1238 0 : getMatchingMachines(ft, cl, vb);
1239 :
1240 : }
1241 : }
1242 : else{
1243 0 : cl=getCL(vb.msId(), vb.fieldId(), vb.spectralWindow());
1244 0 : ft=getFT(vb.msId(), vb.fieldId(), vb.spectralWindow());
1245 :
1246 : }
1247 0 : return;
1248 :
1249 : }
1250 0 : void VisModelData::updateHolders(const VisBuffer& vb, const std::vector<Int>& indexInBuf){
1251 0 : Int fieldId=vb.fieldId();
1252 0 : const MeasurementSet& thems= (vb.getVisibilityIterator())->ms();
1253 0 : Int snum=-1;
1254 0 : Bool hasmodkey=False;
1255 0 : String modelkey;
1256 0 : hasmodkey=isModelDefined(fieldId, thems, modelkey, snum);
1257 0 : if(!hasmodkey){
1258 0 : clindex2_p[indexInBuf]=-1;
1259 0 : ftindex2_p[indexInBuf]=-1;
1260 0 : return;
1261 : }
1262 : //if we have already filled for this field
1263 0 : for (auto it=ftindex2_p.begin(); it != ftindex2_p.end(); ++it){
1264 : // cerr << Vector<Int>(it->first) << " val " << it->second << endl;
1265 0 : if((it->first)[0]==fieldId){
1266 0 : clindex2_p[indexInBuf]=-2;
1267 0 : ftindex2_p[indexInBuf]=-2;
1268 0 : return;
1269 : }
1270 :
1271 : }
1272 : //We do have this key
1273 0 : TableRecord therec;
1274 0 : getModelRecord(modelkey, therec, thems);
1275 0 : if(!therec.isDefined("version") || (therec.asInt("version")!=2)){
1276 : ///Set something versioning so as it can go to do version 1 way
1277 0 : version_p=1;
1278 0 : clindex2_p[indexInBuf]=-1;
1279 0 : ftindex2_p[indexInBuf]=-1;
1280 0 : return;
1281 : }
1282 : else{
1283 0 : version_p=2;
1284 0 : Int indexft=-1;
1285 0 : if(therec.isDefined("numft")){
1286 0 : Int numft=therec.asInt("numft");
1287 0 : if(numft >0){
1288 0 : for(Int ftk=0; ftk < numft; ++ftk){
1289 0 : Record ftrec(therec.asRecord("ft_"+String::toString(ftk)));
1290 0 : modifyDiskImagePath(ftrec, vb);
1291 0 : const Matrix<Int>& ftcombind=ftrec.asArrayInt("indexcombination");
1292 :
1293 0 : indexft=ftholder_p.nelements();
1294 0 : ftholder_p.resize(indexft+1, false, true);
1295 0 : ftholder_p[indexft].resize(1);
1296 0 : ftholder_p[indexft][0]=NEW_FT(ftrec.asRecord("container"));
1297 0 : if(!( ftholder_p[indexft][0]))
1298 0 : throw(AipsError("Unsupported FTMachine found in virtual MODEL_DATA column"));
1299 0 : ftholder_p[indexft][0]->initMaps(vb);
1300 0 : for(uInt row=0; row < ftcombind.nrow(); ++row){
1301 0 : std::vector<int> key=ftcombind.row(row).tovector();
1302 0 : key.push_back(int(vb.msId()));
1303 0 : if(ftindex2_p.count(key) >0){
1304 0 : Int numftforkey=ftholder_p[ftindex2_p[key]].nelements();
1305 0 : Int indx=ftindex2_p[key];
1306 0 : Bool alreadyAdded=false;
1307 0 : for (Int kk=1; kk < numftforkey; ++kk){
1308 0 : alreadyAdded= alreadyAdded || (ftholder_p[indexft][0]==ftholder_p[indx][kk]);
1309 : }
1310 0 : if(!alreadyAdded){
1311 0 : ftholder_p[indx].resize(numftforkey+1, true);
1312 0 : ftholder_p[indx][numftforkey]=ftholder_p[indexft][0];
1313 : }
1314 : }
1315 : else{
1316 0 : ftindex2_p[key]=indexft;
1317 : }
1318 : }
1319 :
1320 : }
1321 : }
1322 : }
1323 0 : if(ftindex2_p.count(indexInBuf)==0)
1324 0 : ftindex2_p[indexInBuf]=-2;
1325 :
1326 0 : Int indexcl=-1;
1327 0 : if(therec.isDefined("numcl")){
1328 0 : Int numcl=therec.asInt("numcl");
1329 0 : if(numcl >0){
1330 0 : for(Int clk=0; clk < numcl; ++clk){
1331 :
1332 0 : Record clrec(therec.asRecord("cl_"+String::toString(clk)));
1333 0 : const Matrix<Int>& clcombind=clrec.asArrayInt("indexcombination");
1334 :
1335 :
1336 0 : indexcl=clholder_p.nelements();
1337 0 : clholder_p.resize(indexcl+1, false, true);
1338 0 : clholder_p[indexcl].resize(1);
1339 0 : clholder_p[indexcl][0]=new ComponentList();
1340 0 : String err;
1341 0 : if(!((clholder_p[indexcl][0])->fromRecord(err, clrec.asRecord("container"))))
1342 0 : throw(AipsError("Component model failed to load for field "+String::toString(clcombind.column(0))));
1343 0 : for(uInt row=0; row < clcombind.nrow(); ++row){
1344 0 : std::vector<int> key=clcombind.row(row).tovector();
1345 0 : key.push_back(int(vb.msId()));
1346 0 : if(clindex2_p.count(key) >0){
1347 0 : Int numclforkey=clholder_p[clindex2_p[key]].nelements();
1348 0 : Int indx=clindex2_p[key];
1349 0 : Bool alreadyAdded=false;
1350 0 : for (Int kk=1; kk < numclforkey; ++kk){
1351 0 : alreadyAdded= alreadyAdded || (clholder_p[indexcl][0]==clholder_p[indx][kk]);
1352 : }
1353 0 : if(!alreadyAdded){
1354 0 : clholder_p[indx].resize(numclforkey+1, true);
1355 0 : clholder_p[indx][numclforkey]=clholder_p[indexcl][0];
1356 : }
1357 : }
1358 : else{
1359 0 : clindex2_p[key]=indexcl;
1360 : }
1361 : }
1362 :
1363 : }
1364 : }
1365 : }
1366 0 : if(clindex2_p.count(indexInBuf)==0)
1367 0 : clindex2_p[indexInBuf]=-2;
1368 : }
1369 : }
1370 :
1371 0 : void VisModelData::init(const VisBuffer& vb){
1372 0 : if(version_p < 1){
1373 0 : updateHolders(vb, std::vector<Int>({vb.fieldId(), vb.spectralWindow(), vb.scan()(0), vb.stateId()(0), vb.msId()}));
1374 :
1375 : }
1376 :
1377 0 : }
1378 :
1379 0 : void VisModelData::getUniqueIndicesComb(const VisBuffer& vb, Vector< Vector<Int> >& retval){
1380 0 : const Vector<Int>& state = vb.stateId();
1381 0 : const Vector<Int>& scan=vb.scan();
1382 0 : const Vector<Double>& t=vb.time();
1383 0 : Int fldid=vb.fieldId();
1384 0 : Int spwid=vb.spectralWindow();
1385 0 : Vector<uInt> uniqIndx;
1386 0 : uInt nTimes=GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates);
1387 0 : Vector<Int> comb(4);
1388 0 : comb(0)=fldid;
1389 0 : comb(1)=spwid;
1390 0 : for (uInt k=0; k < nTimes; ++k){
1391 0 : comb(2)=scan(uniqIndx[k]);
1392 0 : comb(3)=state(uniqIndx[k]);
1393 0 : Bool hasComb=False;
1394 0 : if(retval.nelements() >0){
1395 0 : for (uInt j=0; j < retval.nelements(); ++j){
1396 0 : if(allEQ(retval[j], comb))
1397 0 : hasComb=True;
1398 : }
1399 :
1400 : }
1401 0 : if(!hasComb){
1402 0 : retval.resize(retval.nelements()+1, True);
1403 0 : retval[retval.nelements()-1]=comb;
1404 : }
1405 : }
1406 :
1407 0 : }
1408 0 : Bool VisModelData::getModelVis(VisBuffer& vb){
1409 :
1410 : //Vector<CountedPtr<ComponentList> >cl=getCL(vb.msId(), vb.fieldId(), vb.spectralWindow());
1411 : //Vector<CountedPtr<FTMachine> > ft=getFT(vb.msId(), vb.fieldId(), vb.spectralWindow());
1412 0 : Vector<CountedPtr<ComponentList> >cl(0);
1413 0 : Vector<CountedPtr<FTMachine> > ft(0);
1414 0 : getMatchingMachines(ft, cl, vb);
1415 : //Fill the buffer with 0.0; also prevents reading from disk if MODEL_DATA exists
1416 : ///Oh boy this is really dangerous...
1417 : //nCorr etc are public..who know who changed these values before reaching here.
1418 0 : Cube<Complex> mod(vb.nCorr(), vb.nChannel(), vb.nRow(), Complex(0.0));
1419 0 : vb.setModelVisCube(mod);
1420 0 : Bool incremental=false;
1421 0 : if( cl.nelements()>0){
1422 : //cerr << "In cft " << cl.nelements() << endl;
1423 0 : for (uInt k=0; k < cl.nelements(); ++k)
1424 0 : if(!cl[k].null()){
1425 0 : cft_p->get(vb, *(cl[k]), -1);
1426 : //cerr << "max " << max(vb.modelVisCube()) << endl;
1427 0 : incremental=true;
1428 : }
1429 : }
1430 0 : if(ft.nelements()>0){
1431 0 : Cube<Complex> tmpModel;
1432 0 : if(incremental || ft.nelements() >1)
1433 0 : tmpModel.assign(vb.modelVisCube());
1434 0 : Bool allnull=true;
1435 0 : for (uInt k=0; k < ft.nelements(); ++k){
1436 0 : if(!ft[k].null()){
1437 0 : if(k >0) vb.setModelVisCube(Cube<Complex> (vb.nCorr(), vb.nChannel(), vb.nRow(), Complex(0.0)));
1438 :
1439 0 : ft[k]->get(vb, -1);
1440 0 : if(ft.nelements()>1 || incremental){
1441 0 : tmpModel+=vb.modelVisCube();
1442 : }
1443 0 : allnull=false;
1444 : }
1445 : }
1446 : //cerr << "min max after ft " << min(vb.modelVisCube()) << max(vb.modelVisCube()) << endl;
1447 0 : if(!allnull){
1448 0 : if(ft.nelements()>1 || incremental)
1449 0 : vb.modelVisCube()=tmpModel;
1450 0 : incremental=true;
1451 : }
1452 : }
1453 0 : if(!incremental){
1454 : //No model was set so....
1455 : ///Set the Model to 1.0 for parallel hand and 0.0 for x-hand
1456 :
1457 0 : vb.modelVisCube().set(Complex(1.0));
1458 0 : Vector<Int> corrType=vb.corrType();
1459 0 : uInt nCorr = corrType.nelements();
1460 0 : for (uInt i=0; i<nCorr; i++) {
1461 0 : if(corrType[i]==Stokes::RL || corrType[i]==Stokes::LR ||
1462 0 : corrType[i]==Stokes::XY || corrType[i]==Stokes::YX){
1463 0 : vb.modelVisCube().yzPlane(i)=0.0;
1464 : }
1465 : }
1466 : }
1467 :
1468 0 : return true;
1469 :
1470 : }
1471 :
1472 :
1473 0 : Vector<CountedPtr<ComponentList> > VisModelData::getCL(const Int msId, const Int fieldId, const Int spwId){
1474 0 : if(!hasModel(msId, fieldId, spwId))
1475 0 : return Vector<CountedPtr<ComponentList> >(0);
1476 0 : Int indx=clindex_p(spwId, fieldId, msId);
1477 : //cerr << "indx " << indx << " " << clholder_p[indx].nelements() << " spw " << spwId << endl;
1478 0 : if(indx <0)
1479 0 : return Vector<CountedPtr<ComponentList> >(0);
1480 0 : return clholder_p[indx];
1481 :
1482 :
1483 : }
1484 :
1485 0 : Vector<CountedPtr<FTMachine> >VisModelData::getFT(const Int msId, const Int fieldId, Int spwId){
1486 :
1487 0 : if(!hasModel(msId, fieldId, spwId))
1488 0 : return Vector<CountedPtr<FTMachine> >(0);
1489 0 : Int indx=ftindex_p(spwId, fieldId, msId);
1490 : //cerr << "indx " << indx << endl;
1491 0 : if(indx <0)
1492 0 : return Vector<CountedPtr<FTMachine> >(0);
1493 0 : return ftholder_p[indx];
1494 : }
1495 :
1496 0 : Int VisModelData::firstSourceRowRecord(const Int field, const MeasurementSet& theMS, TableRecord& rec){
1497 0 : Int row=-1;
1498 :
1499 : //Prefer the Source table first
1500 0 : MSFieldColumns fCol(theMS.field());
1501 0 : if(Table::isReadable(theMS.sourceTableName()) && (theMS.source().nrow() > 0) && (!fCol.sourceId().isNull()) && (fCol.sourceId().get(field) != -1) ){
1502 :
1503 0 : const MSSource& mss=theMS.source();
1504 :
1505 0 : Int sid=fCol.sourceId().get(field);
1506 0 : auto rows=MSSourceIndex(mss).getRowNumbersOfSourceID(sid);
1507 0 : if(rows.nelements() > 0)
1508 0 : row=rows[0];
1509 0 : const TableRecord& keywords=mss.keywordSet();
1510 0 : Bool rowIsUsed=false;
1511 0 : for (uInt n=0; n< keywords.nfields(); ++n){
1512 0 : if(keywords.dataType(n) == TpInt){
1513 0 : if(row==keywords.asInt(n)){
1514 0 : rowIsUsed=true;
1515 : //cerr << "has record pointed to" << endl;
1516 : }
1517 : }
1518 : }
1519 : //nobody is using that row ..any record there, if any, must be dangling
1520 0 : if(!rowIsUsed)
1521 0 : return -1;
1522 0 : if(row >-1 && mss.isColumn(MSSource::SOURCE_MODEL)){
1523 0 : ScalarColumn<TableRecord> scol(theMS.source(), "SOURCE_MODEL");
1524 0 : scol.get(row, rec);
1525 : }
1526 : }
1527 0 : return row;
1528 : }
1529 :
1530 :
1531 : }//# NAMESPACE CASA - END
1532 :
|
