Commits

Kumar Golap authored acbe7be224a
remove printouts
No tags

casa5/code/synthesis/TransformMachines2/BriggsCubeWeightor.cc

Modified
135 135
136 136 pair<Int, Int> ms_field=make_pair(vb->msId(), vb->fieldId()[0]);
137 137 if (std::find(fieldsToUse.begin(), fieldsToUse.end(), ms_field) == fieldsToUse.end()) {
138 138 fieldsToUse.push_back(ms_field);
139 139 }
140 140
141 141 }
142 142 }
143 143 }
144 144 wgtTab_p=nullptr;
145 + Int tmpInt;
146 + inRec.get("freqinterpmethod", tmpInt);
147 + freqInterpMethod_p=static_cast<InterpolateArray1D<Double, Complex >::InterpolationMethod>(tmpInt);
145 148 ostringstream oss;
146 - oss << std::setprecision(12) << nrows << "_" << freqbeg << "_" << freqend << "_"<< rmode_p << "_" << robust_p <<"_interp_"<< freqInterpMethod_p;
149 + oss << std::setprecision(12) << nrows << "_" << freqbeg << "_" << freqend << "_"<< rmode_p << "_" << robust_p << "_interp_"<< tmpInt ;
147 150
148 151 //cerr << "STRING " << oss.str() << endl;;
149 152 imWgtColName_p=makeScratchImagingWeightTable(wgtTab_p, oss.str());
150 153 if(wgtTab_p->nrow()==nrows){
151 154 //cerr << "REUSING " << endl;
152 155 return imWgtColName_p;
153 156 }
154 157 else{
155 158 wgtTab_p->lock(True,100);
156 159 wgtTab_p->removeRow(wgtTab_p->rowNumbers());
157 160
158 161 }
159 162 //cerr << "CREATING " << endl;
160 163 vbrowms2wgtrow_p.clear();
161 164 if(fieldsToUse.size()==0)
162 165 fieldsToUse.push_back(make_pair(Int(-1),Int(-1)));
163 166 //cerr << "FIELDs to use " << Vector<pair<Int,Int> >(fieldsToUse) << endl;
164 167
168 + Bool inOneGo=True;
169 + uInt allSwingPad=estimateSwingChanPad(vi, -1, cs, templateimage.shape()[3],
170 + ephemtab);
171 + if(msInUse.size() > 1){
172 +
173 + if((allSwingPad > 4) && (allSwingPad > uInt(templateimage.shape()[3]/10)))
174 + inOneGo=False;
175 + //cerr << "allSwingPad " << allSwingPad << " inOneGo " << inOneGo << endl;
176 +
177 + }
178 +
179 +
165 180 ///Lets process the ms independently as swingpad can become very large for MSs seperated by large epochs
166 - for (auto msiter=msInUse.begin(); msiter != msInUse.end(); ++msiter){
167 - uInt swingpad=estimateSwingChanPad(vi, *msiter, cs, templateimage.shape()[3],
181 + if(inOneGo){
182 + fillImgWeightCol(vi, inRec, -1, fieldsToUse, allSwingPad, templateimage.shape(), cs);
183 + }
184 + else{
185 +
186 + for (auto msiter=msInUse.begin(); msiter != msInUse.end(); ++msiter){
187 + uInt swingpad=estimateSwingChanPad(vi, *msiter, cs, templateimage.shape()[3],
168 188 ephemtab);
169 - for (uInt k=0; k < fieldsToUse.size(); ++k){
170 - vi.originChunks();
171 - vi.origin();
172 - //cerr << "####nchannels " << vb->nChannels() << " swingpad " << swingpad << endl;
173 - IPosition shp=templateimage.shape();
174 - nx_p=shp[0];
175 - ny_p=shp[1];
176 - CoordinateSystem cs=templateimage.coordinates();
177 - refFreq_p=cs.toWorld(IPosition(4,0,0,0,0))[3];
178 - Vector<String> units = cs.worldAxisUnits();
179 - units[0]="rad"; units[1]="rad";
180 - cs.setWorldAxisUnits(units);
181 - Vector<Double> incr=cs.increment();
182 - uscale_p=(nx_p*incr[0]);
183 - vscale_p=(ny_p*incr[1]);
184 - uorigin_p=nx_p/2;
185 - vorigin_p=ny_p/2;
186 - //IPosition shp=templateimage.shape();
187 - shp[3]=shp[3]+swingpad; //add extra channels at begining and end;
188 - Vector<Double> refpix=cs.referencePixel();
189 - refpix[3]+=swingpad/2;
190 - cs.setReferencePixel(refpix);
191 - cerr << "REFPIX " << refpix << " new SHAPE " << shp << " swigpad " << swingpad << " msid " <<fieldsToUse[k].first << " fieldid " << fieldsToUse[k].second << endl;
192 - TempImage<Complex> newTemplate(shp, cs);
193 - Matrix<Double> sumWgts;
189 + fillImgWeightCol(vi, inRec, *msiter, fieldsToUse, swingpad, templateimage.shape(), cs);
190 +
191 + }
192 + }
193 + initialized_p=True;
194 + wgtTab_p->unlock();
195 + return imWgtColName_p;
196 +
197 +}
198 +
199 + void BriggsCubeWeightor::fillImgWeightCol(vi::VisibilityIterator2& vi, const Record& inRec, const Int msid, std::vector<pair<Int, Int> >& fieldsToUse, const uInt swingpad, const IPosition& origShp, CoordinateSystem cs){
200 + vi::VisBuffer2 *vb=vi.getVisBuffer();
201 +
202 + for (uInt k=0; k < fieldsToUse.size(); ++k){
203 + vi.originChunks();
204 + vi.origin();
205 + //cerr << "####nchannels " << vb->nChannels() << " swingpad " << swingpad << endl;
206 + IPosition shp=origShp;
207 + nx_p=shp[0];
208 + ny_p=shp[1];
209 + //CoordinateSystem cs=templateimage.coordinates();
210 + refFreq_p=cs.toWorld(IPosition(4,0,0,0,0))[3];
211 + Vector<String> units = cs.worldAxisUnits();
212 + units[0]="rad"; units[1]="rad";
213 + cs.setWorldAxisUnits(units);
214 + Vector<Double> incr=cs.increment();
215 + uscale_p=(nx_p*incr[0]);
216 + vscale_p=(ny_p*incr[1]);
217 + uorigin_p=nx_p/2;
218 + vorigin_p=ny_p/2;
219 + //IPosition shp=templateimage.shape();
220 + shp[3]=shp[3]+swingpad; //add extra channels at begining and end;
221 + Vector<Double> refpix=cs.referencePixel();
222 + refpix[3]+=swingpad/2;
223 + cs.setReferencePixel(refpix);
224 + //cerr << "REFPIX " << refpix << " new SHAPE " << shp << " swigpad " << swingpad << " msid " <<fieldsToUse[k].first << " fieldid " << fieldsToUse[k].second << endl;
225 + TempImage<Complex> newTemplate(shp, cs);
226 + Matrix<Double> sumWgts;
194 227
195 - initializeFTMachine(0, newTemplate, inRec);
196 - Matrix<Float> dummy;
197 - //cerr << "new template shape " << newTemplate.shape() << endl;
198 - ft_p[0]->initializeToSky(newTemplate, dummy, *vb);
199 - Vector<Double> convFunc(2+superUniformBox_p, 1.0);
200 - //cerr << "superuniform box " << superUniformBox_p << endl;
201 - ft_p[0]->modifyConvFunc(convFunc, superUniformBox_p, 1);
202 - for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
203 - for (vi.origin(); vi.more(); vi.next()) {
204 - //cerr << "key and index "<< key << " " << index << " " << multiFieldMap_p[key] << endl;
205 - if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){
206 - ft_p[0]->put(*vb, -1, true, FTMachine::PSF);
207 - }
208 -
209 - }
210 - }
211 - Array<Float> griddedWeight;
212 - ft_p[0]->getGrid(griddedWeight);
213 - //cerr << " griddedWeight Shape " << griddedWeight.shape() << endl;
214 - //grids_p[index]->put(griddedWeight.reform(newTemplate.shape()));
215 - sumWgts=ft_p[0]->getSumWeights();
216 - //cerr << "sumweight " << sumWgts[index] << endl;
217 - //clear the ftmachine
218 - ft_p[0]->finalizeToSky();
228 + initializeFTMachine(0, newTemplate, inRec);
229 + Matrix<Float> dummy;
230 + //cerr << "new template shape " << newTemplate.shape() << endl;
231 + ft_p[0]->initializeToSky(newTemplate, dummy, *vb);
232 + Vector<Double> convFunc(2+superUniformBox_p, 1.0);
233 + //cerr << "superuniform box " << superUniformBox_p << endl;
234 + ft_p[0]->modifyConvFunc(convFunc, superUniformBox_p, 1);
235 + for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
236 + for (vi.origin(); vi.more(); vi.next()) {
237 + //cerr << "key and index "<< key << " " << index << " " << multiFieldMap_p[key] << endl;
238 + if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){
239 + ft_p[0]->put(*vb, -1, true, FTMachine::PSF);
240 + }
241 +
242 + }
243 + }
244 + Array<Float> griddedWeight;
245 + ft_p[0]->getGrid(griddedWeight);
246 + //cerr << " griddedWeight Shape " << griddedWeight.shape() << endl;
247 + //grids_p[index]->put(griddedWeight.reform(newTemplate.shape()));
248 + sumWgts=ft_p[0]->getSumWeights();
249 + //cerr << "sumweight " << sumWgts[index] << endl;
250 + //clear the ftmachine
251 + ft_p[0]->finalizeToSky();
219 252
220 253
221 - Int nchan=newTemplate.shape()(3);
222 - //cerr << "rmode " << rmode_p << endl;
223 -
224 - for (uInt chan=0; chan < uInt(nchan); ++ chan){
225 - IPosition start(4,0,0,0,chan);
226 - IPosition end(4,nx_p-1,ny_p-1,0,chan);
227 - Matrix<Float> gwt(griddedWeight(start, end).reform(IPosition(2, nx_p, ny_p)));
228 - if ((rmode_p=="norm" || rmode_p=="bwtaper") && (sumWgts(0,chan)> 0.0)) { //See CAS-13021 for bwtaper algorithm details
229 - //os << "Normal robustness, robust = " << robust << LogIO::POST;
230 - Double sumlocwt = 0.;
231 - for(Int vgrid=0;vgrid<ny_p;vgrid++) {
232 - for(Int ugrid=0;ugrid<nx_p;ugrid++) {
233 - if(gwt(ugrid, vgrid)>0.0) sumlocwt+=square(gwt(ugrid,vgrid));
234 - }
235 - }
236 - f2_p[0][chan] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumWgts(0,chan));
237 - d2_p[0][chan] = 1.0;
238 -
239 - }
240 - else if (rmode_p=="abs") {
241 - //os << "Absolute robustness, robust = " << robust << ", noise = "
242 - // << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
243 - f2_p[0][chan] = square(robust_p);
244 - d2_p[0][chan] = 2.0 * square(noise_p.get("Jy").getValue());
245 -
246 - }
247 - else {
248 - f2_p[0][chan] = 1.0;
249 - d2_p[0][chan] = 0.0;
254 + Int nchan=newTemplate.shape()(3);
255 + //cerr << "rmode " << rmode_p << endl;
256 +
257 + for (uInt chan=0; chan < uInt(nchan); ++ chan){
258 + IPosition start(4,0,0,0,chan);
259 + IPosition end(4,nx_p-1,ny_p-1,0,chan);
260 + Matrix<Float> gwt(griddedWeight(start, end).reform(IPosition(2, nx_p, ny_p)));
261 + if ((rmode_p=="norm" || rmode_p=="bwtaper") && (sumWgts(0,chan)> 0.0)) { //See CAS-13021 for bwtaper algorithm details
262 + //os << "Normal robustness, robust = " << robust << LogIO::POST;
263 + Double sumlocwt = 0.;
264 + for(Int vgrid=0;vgrid<ny_p;vgrid++) {
265 + for(Int ugrid=0;ugrid<nx_p;ugrid++) {
266 + if(gwt(ugrid, vgrid)>0.0) sumlocwt+=square(gwt(ugrid,vgrid));
267 + }
268 + }
269 + f2_p[0][chan] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumWgts(0,chan));
270 + d2_p[0][chan] = 1.0;
271 +
272 + }
273 + else if (rmode_p=="abs") {
274 + //os << "Absolute robustness, robust = " << robust << ", noise = "
275 + // << noise.get("Jy").getValue() << "Jy" << LogIO::POST;
276 + f2_p[0][chan] = square(robust_p);
277 + d2_p[0][chan] = 2.0 * square(noise_p.get("Jy").getValue());
278 +
279 + }
280 + else {
281 + f2_p[0][chan] = 1.0;
282 + d2_p[0][chan] = 0.0;
250 283 }
251 -
252 - }//chan
284 +
285 + }//chan
253 286
254 287
255 - //std::ofstream myfile;
256 - //myfile.open (wgtTab_p->tableName()+".txt");
257 - vi.originChunks();
258 - vi.origin();
259 - for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
260 - for (vi.origin(); vi.more(); vi.next()) {
261 - if((vb->msId()) == (*msiter)){
262 - if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){
263 - Matrix<Float> imweight;
264 - /*///////////
265 - std::vector<Int> cmap=(ft_p[0]->channelMap(*vb)).tovector();
266 - int n=0;
267 - std::for_each(cmap.begin(), cmap.end(), [&n, &myfile](int &val){if(val > -1)myfile << " " << n; n++; });
268 - myfile << "----msid " << vb->msId() << endl;
288 + //std::ofstream myfile;
289 + //myfile.open (wgtTab_p->tableName()+".txt");
290 + vi.originChunks();
291 + vi.origin();
292 + for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
293 + for (vi.origin(); vi.more(); vi.next()) {
294 + if((msid < 0) || ((vb->msId()) == (msid))){
295 + if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){
296 + Matrix<Float> imweight;
297 + /*///////////
298 + std::vector<Int> cmap=(ft_p[0]->channelMap(*vb)).tovector();
299 + int n=0;
300 + std::for_each(cmap.begin(), cmap.end(), [&n, &myfile](int &val){if(val > -1)myfile << " " << n; n++; });
301 + myfile << "----msid " << vb->msId() << endl;
269 302 /////////////////////////////*/
270 - getWeightUniform(griddedWeight, imweight, *vb);
271 - rownr_t nRows=vb->nRows();
272 - //Int nChans=vb->nChannels();
273 - Vector<uInt> msId(nRows, uInt(vb->msId()));
274 - RowNumbers rowidsnr=vb->rowIds();
275 - Vector<uInt>rowids(rowidsnr.nelements());
276 - convertArray(rowids, rowidsnr);
277 - wgtTab_p->addRow(nRows, False);
278 - rownr_t endrow=wgtTab_p->nrow()-1;
279 - rownr_t beginrow=endrow-nRows+1;
280 - //Slicer sl(IPosition(2,beginrow,0), IPosition(2,endrow,nChans-1), Slicer::endIsLast);
281 - ArrayColumn<Float> col(*wgtTab_p, "IMAGING_WEIGHT");
282 - //cerr << "sl length " << sl.length() << " array shape " << fakeweight.shape() << " col length " << col.nrow() << endl;
283 - //col.putColumnRange(sl, fakeweight);
284 - //cerr << "nrows " << nRows << " imweight.shape " << imweight.shape() << endl;
285 - for (rownr_t row=0; row < nRows; ++row)
286 - col.put(beginrow+row, imweight.column(row));
287 -
288 - Slicer sl2(IPosition(1, endrow-nRows+1), IPosition(1, endrow), Slicer::endIsLast);
289 - ScalarColumn<uInt> col2(*wgtTab_p,"MSID");
290 - col2.putColumnRange(sl2, msId);
291 - ScalarColumn<uInt> col3(*wgtTab_p, "ROWID");
292 -
293 - col3.putColumnRange(sl2, rowids);
294 - }
295 - }
296 - }
303 + getWeightUniform(griddedWeight, imweight, *vb);
304 + rownr_t nRows=vb->nRows();
305 + //Int nChans=vb->nChannels();
306 + Vector<uInt> msId(nRows, uInt(vb->msId()));
307 + RowNumbers rowidsnr=vb->rowIds();
308 + Vector<uInt>rowids(rowidsnr.nelements());
309 + convertArray(rowids, rowidsnr);
310 + wgtTab_p->addRow(nRows, False);
311 + rownr_t endrow=wgtTab_p->nrow()-1;
312 + rownr_t beginrow=endrow-nRows+1;
313 + //Slicer sl(IPosition(2,beginrow,0), IPosition(2,endrow,nChans-1), Slicer::endIsLast);
314 + ArrayColumn<Float> col(*wgtTab_p, "IMAGING_WEIGHT");
315 + //cerr << "sl length " << sl.length() << " array shape " << fakeweight.shape() << " col length " << col.nrow() << endl;
316 + //col.putColumnRange(sl, fakeweight);
317 + //cerr << "nrows " << nRows << " imweight.shape " << imweight.shape() << endl;
318 + for (rownr_t row=0; row < nRows; ++row)
319 + col.put(beginrow+row, imweight.column(row));
320 +
321 + Slicer sl2(IPosition(1, endrow-nRows+1), IPosition(1, endrow), Slicer::endIsLast);
322 + ScalarColumn<uInt> col2(*wgtTab_p,"MSID");
323 + col2.putColumnRange(sl2, msId);
324 + ScalarColumn<uInt> col3(*wgtTab_p, "ROWID");
325 +
326 + col3.putColumnRange(sl2, rowids);
297 327 }
328 + }
329 + }
330 + }
298 331 //myfile.close();
299 332
300 333
301 - }
302 - }
334 + }
335 +
303 336 ////Lets reset vi before returning
304 - vi.originChunks();
305 - vi.origin();
306 -
307 - initialized_p=True;
308 - wgtTab_p->unlock();
309 - return imWgtColName_p;
310 -
311 -}
337 + vi.originChunks();
338 + vi.origin();
339 +
340 +
341 +
342 +
343 +
344 + }
345 +
312 346
313 347 Int BriggsCubeWeightor::estimateSwingChanPad(vi::VisibilityIterator2& vi, const Int msid, const CoordinateSystem& cs, const Int imNChan, const String& ephemtab){
314 348
315 349
316 350 vi.originChunks();
317 351 vi.origin();
318 352 vi::VisBuffer2 *vb=vi.getVisBuffer();
319 353 Double freqbeg=cs.toWorld(IPosition(4,0,0,0,0))[3];
320 354 Double freqend=cs.toWorld(IPosition(4,0,0,0,imNChan-1))[3];
321 355 Double freqincr=fabs(cs.increment()[3]);
322 356 SpectralCoordinate spCoord=cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL));
323 357 MFrequency::Types freqframe=spCoord.frequencySystem(True);
324 358 uInt swingpad=16;
325 359 Double swingFreq=0.0;
326 360 Double minFreq=1e99;
327 361 Double maxFreq=0.0;
328 362 std::vector<Double> localminfreq, localmaxfreq, firstchanfreq;
329 - Int msID=-1;
363 + //Int msID=-1;
330 364 Int fieldID=-1;
331 365 Int spwID=-1;
332 366 ////TESTOO
333 367 //int my_cpu_id;
334 368 //MPI_Comm_rank(MPI_COMM_WORLD, &my_cpu_id);
335 369 ///////////////////
336 370 ///////
337 371 for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
338 372 for (vi.origin(); vi.more(); vi.next()) {
339 373 //process for required msid
340 - //if(msid==vb->msId()){
341 - {
342 - if((msID != vb->msId()) || (fieldID != vb->fieldId()(0)) || (spwID!=vb->spectralWindows()(0))){
343 - msID=vb->msId();
374 + if((msid < 0) || (msid==vb->msId()))
375 + {
376 + //if((msID != vb->msId()) || (fieldID != vb->fieldId()(0)) || (spwID!=vb->spectralWindows()(0))){
377 + {
378 + //msID=vb->msId();
344 379 fieldID = vb->fieldId()(0);
345 380 spwID = vb->spectralWindows()(0);
346 381 Double localBeg, localEnd;
347 382 Double localNchan=imNChan >1 ? Double(imNChan-1) : 1.0;
348 383 Double localStep=abs(freqend-freqbeg)/localNchan;
349 384 if(freqbeg < freqend){
350 385 localBeg=freqbeg;
351 386 localEnd=freqend;
352 387 }
353 388 else{
405 440 swingFreq=abs(*itmin -minFreq);
406 441 if(swingFreq < abs(*itmax -maxFreq))
407 442 swingFreq=abs(*itmax -maxFreq);
408 443 if(firstchanshift < abs(*itf-minfirstchan))
409 444 firstchanshift=abs(*itf-minfirstchan);
410 445 itf++;
411 446 itmax++;
412 447 }
413 448 Int extrapad=max(min(4, Int(imNChan/10)),1);
414 449 swingpad=2*(Int(std::ceil((swingFreq+firstchanshift)/freqincr))+extrapad);
415 - cerr <<" swingfreq " << (swingFreq/freqincr) << " firstchanshift " << (firstchanshift/freqincr) << " SWINGPAD " << swingpad << endl;
450 + //cerr <<" swingfreq " << (swingFreq/freqincr) << " firstchanshift " << (firstchanshift/freqincr) << " SWINGPAD " << swingpad << endl;
416 451 ////////////////
417 452 return swingpad;
418 453
419 454 }
420 455
421 456 String BriggsCubeWeightor::makeScratchImagingWeightTable(CountedPtr<Table>& weightTable, const String& filetag){
422 457
423 458 //String wgtname=File::newUniqueName(".", "IMAGING_WEIGHT").absoluteName();
424 459 String wgtname=Path("IMAGING_WEIGHT_"+filetag).absoluteName();
425 460 //cerr << "NAME " << wgtname << endl;

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut