Commits
107 107 | |
108 108 | CoordinateSystem cs=templateimage.coordinates(); |
109 109 | if(initialized_p && nx_p==templateimage.shape()(0) && ny_p==templateimage.shape()(1)){ |
110 110 | |
111 111 | Double freq=cs.toWorld(IPosition(4,0,0,0,0))[3]; |
112 112 | if(freq==refFreq_p) |
113 113 | return imWgtColName_p; |
114 114 | } |
115 115 | |
116 116 | visWgt_p=vi.getImagingWeightGenerator(); |
117 - | VisImagingWeight vWghtNat("natural"); |
117 + | VisImagingWeight vWghtNat("natural"); |
118 118 | vi.useImagingWeight(vWghtNat); |
119 119 | vi::VisBuffer2 *vb=vi.getVisBuffer(); |
120 120 | std::vector<pair<Int, Int> > fieldsToUse; |
121 121 | std::set<Int> msInUse; |
122 122 | rownr_t nrows=0; |
123 123 | String ephemtab(""); |
124 124 | if(inRec.isDefined("ephemeristable")){ |
125 125 | inRec.get("ephemeristable", ephemtab); |
126 126 | } |
127 127 | Double freqbeg=cs.toWorld(IPosition(4,0,0,0,0))[3]; |
158 158 | } |
159 159 | //cerr << "CREATING " << endl; |
160 160 | vbrowms2wgtrow_p.clear(); |
161 161 | if(fieldsToUse.size()==0) |
162 162 | fieldsToUse.push_back(make_pair(Int(-1),Int(-1))); |
163 163 | //cerr << "FIELDs to use " << Vector<pair<Int,Int> >(fieldsToUse) << endl; |
164 164 | |
165 165 | ///Lets process the ms independently as swingpad can become very large for MSs seperated by large epochs |
166 166 | for (auto msiter=msInUse.begin(); msiter != msInUse.end(); ++msiter){ |
167 167 | uInt swingpad=estimateSwingChanPad(vi, *msiter, cs, templateimage.shape()[3], |
168 - | 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 - | TempImage<Complex> newTemplate(shp, cs); |
192 - | Matrix<Double> sumWgts; |
193 - | |
194 - | initializeFTMachine(0, newTemplate, inRec); |
195 - | Matrix<Float> dummy; |
196 - | //cerr << "new template shape " << newTemplate.shape() << endl; |
197 - | ft_p[0]->initializeToSky(newTemplate, dummy, *vb); |
198 - | Vector<Double> convFunc(2+superUniformBox_p, 1.0); |
199 - | //cerr << "superuniform box " << superUniformBox_p << endl; |
200 - | ft_p[0]->modifyConvFunc(convFunc, superUniformBox_p, 1); |
201 - | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
202 - | for (vi.origin(); vi.more(); vi.next()) { |
203 - | //cerr << "key and index "<< key << " " << index << " " << multiFieldMap_p[key] << endl; |
204 - | if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){ |
205 - | ft_p[0]->put(*vb, -1, true, FTMachine::PSF); |
206 - | } |
207 - | |
208 - | } |
209 - | } |
210 - | Array<Float> griddedWeight; |
211 - | ft_p[0]->getGrid(griddedWeight); |
212 - | //cerr << " griddedWeight Shape " << griddedWeight.shape() << endl; |
213 - | //grids_p[index]->put(griddedWeight.reform(newTemplate.shape())); |
214 - | sumWgts=ft_p[0]->getSumWeights(); |
215 - | //cerr << "sumweight " << sumWgts[index] << endl; |
216 - | //clear the ftmachine |
217 - | ft_p[0]->finalizeToSky(); |
218 - | |
219 - | ////Lets reset vi before returning |
220 - | vi.originChunks(); |
221 - | vi.origin(); |
222 - | |
223 - | Int nchan=newTemplate.shape()(3); |
224 - | //cerr << "rmode " << rmode_p << endl; |
225 - | |
226 - | for (uInt chan=0; chan < uInt(nchan); ++ chan){ |
227 - | IPosition start(4,0,0,0,chan); |
228 - | IPosition end(4,nx_p-1,ny_p-1,0,chan); |
229 - | Matrix<Float> gwt(griddedWeight(start, end).reform(IPosition(2, nx_p, ny_p))); |
230 - | if ((rmode_p=="norm" || rmode_p=="bwtaper") && (sumWgts(0,chan)> 0.0)) { //See CAS-13021 for bwtaper algorithm details |
231 - | //os << "Normal robustness, robust = " << robust << LogIO::POST; |
232 - | Double sumlocwt = 0.; |
233 - | for(Int vgrid=0;vgrid<ny_p;vgrid++) { |
234 - | for(Int ugrid=0;ugrid<nx_p;ugrid++) { |
235 - | if(gwt(ugrid, vgrid)>0.0) sumlocwt+=square(gwt(ugrid,vgrid)); |
236 - | } |
237 - | } |
238 - | f2_p[0][chan] = square(5.0*pow(10.0,Double(-robust_p))) / (sumlocwt / sumWgts(0,chan)); |
239 - | d2_p[0][chan] = 1.0; |
240 - | |
241 - | } |
242 - | else if (rmode_p=="abs") { |
243 - | //os << "Absolute robustness, robust = " << robust << ", noise = " |
244 - | // << noise.get("Jy").getValue() << "Jy" << LogIO::POST; |
245 - | f2_p[0][chan] = square(robust_p); |
246 - | d2_p[0][chan] = 2.0 * square(noise_p.get("Jy").getValue()); |
247 - | |
248 - | } |
249 - | else { |
250 - | f2_p[0][chan] = 1.0; |
251 - | d2_p[0][chan] = 0.0; |
252 - | } |
253 - | |
254 - | }//chan |
255 - | |
256 - | |
257 - | //std::ofstream myfile; |
258 - | //myfile.open (wgtTab_p->tableName()+".txt"); |
259 - | |
260 - | |
261 - | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
262 - | for (vi.origin(); vi.more(); vi.next()) { |
263 - | if((vb->msId()) == (*msiter)){ |
264 - | if((vb->fieldId()[0] == fieldsToUse[k].second && vb->msId()== fieldsToUse[k].first) || fieldsToUse[k].first==-1){ |
265 - | Matrix<Float> imweight; |
266 - | /*/////////// |
267 - | std::vector<Int> cmap=(ft_p[0]->channelMap(*vb)).tovector(); |
268 - | int n=0; |
269 - | std::for_each(cmap.begin(), cmap.end(), [&n, &myfile](int &val){if(val > -1)myfile << " " << n; n++; }); |
270 - | myfile << "----msid " << vb->msId() << endl; |
271 - | /////////////////////////////*/ |
272 - | getWeightUniform(griddedWeight, imweight, *vb); |
273 - | rownr_t nRows=vb->nRows(); |
274 - | //Int nChans=vb->nChannels(); |
275 - | Vector<uInt> msId(nRows, uInt(vb->msId())); |
276 - | RowNumbers rowidsnr=vb->rowIds(); |
277 - | Vector<uInt>rowids(rowidsnr.nelements()); |
278 - | convertArray(rowids, rowidsnr); |
279 - | wgtTab_p->addRow(nRows, False); |
280 - | rownr_t endrow=wgtTab_p->nrow()-1; |
281 - | rownr_t beginrow=endrow-nRows+1; |
282 - | //Slicer sl(IPosition(2,beginrow,0), IPosition(2,endrow,nChans-1), Slicer::endIsLast); |
283 - | ArrayColumn<Float> col(*wgtTab_p, "IMAGING_WEIGHT"); |
284 - | //cerr << "sl length " << sl.length() << " array shape " << fakeweight.shape() << " col length " << col.nrow() << endl; |
285 - | //col.putColumnRange(sl, fakeweight); |
286 - | //cerr << "nrows " << nRows << " imweight.shape " << imweight.shape() << endl; |
287 - | for (rownr_t row=0; row < nRows; ++row) |
288 - | col.put(beginrow+row, imweight.column(row)); |
289 - | |
290 - | Slicer sl2(IPosition(1, endrow-nRows+1), IPosition(1, endrow), Slicer::endIsLast); |
291 - | ScalarColumn<uInt> col2(*wgtTab_p,"MSID"); |
292 - | col2.putColumnRange(sl2, msId); |
293 - | ScalarColumn<uInt> col3(*wgtTab_p, "ROWID"); |
294 - | |
295 - | col3.putColumnRange(sl2, rowids); |
296 - | } |
297 - | } |
298 - | } |
299 - | } |
168 + | 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; |
194 + | |
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(); |
219 + | |
220 + | |
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; |
250 + | } |
251 + | |
252 + | }//chan |
253 + | |
254 + | |
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; |
269 + | /////////////////////////////*/ |
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 + | } |
297 + | } |
300 298 | //myfile.close(); |
301 299 | |
302 - | |
303 - | } |
300 + | |
301 + | } |
304 302 | } |
305 - | |
303 + | ////Lets reset vi before returning |
304 + | vi.originChunks(); |
305 + | vi.origin(); |
306 + | |
306 307 | initialized_p=True; |
307 308 | wgtTab_p->unlock(); |
308 309 | return imWgtColName_p; |
309 310 | |
310 311 | } |
311 312 | |
312 313 | Int BriggsCubeWeightor::estimateSwingChanPad(vi::VisibilityIterator2& vi, const Int msid, const CoordinateSystem& cs, const Int imNChan, const String& ephemtab){ |
313 314 | |
314 - | vi::VisBuffer2 *vb=vi.getVisBuffer(); |
315 + | |
315 316 | vi.originChunks(); |
316 317 | vi.origin(); |
318 + | vi::VisBuffer2 *vb=vi.getVisBuffer(); |
317 319 | Double freqbeg=cs.toWorld(IPosition(4,0,0,0,0))[3]; |
318 320 | Double freqend=cs.toWorld(IPosition(4,0,0,0,imNChan-1))[3]; |
319 321 | Double freqincr=fabs(cs.increment()[3]); |
320 322 | SpectralCoordinate spCoord=cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)); |
321 323 | MFrequency::Types freqframe=spCoord.frequencySystem(True); |
322 324 | uInt swingpad=16; |
323 325 | Double swingFreq=0.0; |
324 326 | Double minFreq=1e99; |
325 327 | Double maxFreq=0.0; |
326 328 | std::vector<Double> localminfreq, localmaxfreq, firstchanfreq; |
327 329 | Int msID=-1; |
328 330 | Int fieldID=-1; |
329 331 | Int spwID=-1; |
330 332 | ////TESTOO |
331 333 | //int my_cpu_id; |
332 334 | //MPI_Comm_rank(MPI_COMM_WORLD, &my_cpu_id); |
333 335 | /////////////////// |
334 336 | /////// |
335 337 | for (vi.originChunks();vi.moreChunks();vi.nextChunk()) { |
336 338 | for (vi.origin(); vi.more(); vi.next()) { |
337 339 | //process for required msid |
338 - | if(msid==vb->msId()){ |
339 - | |
340 + | //if(msid==vb->msId()){ |
341 + | { |
340 342 | if((msID != vb->msId()) || (fieldID != vb->fieldId()(0)) || (spwID!=vb->spectralWindows()(0))){ |
341 343 | msID=vb->msId(); |
342 344 | fieldID = vb->fieldId()(0); |
343 345 | spwID = vb->spectralWindows()(0); |
344 346 | Double localBeg, localEnd; |
345 347 | Double localNchan=imNChan >1 ? Double(imNChan-1) : 1.0; |
346 348 | Double localStep=abs(freqend-freqbeg)/localNchan; |
347 349 | if(freqbeg < freqend){ |
348 350 | localBeg=freqbeg; |
349 351 | localEnd=freqend; |
403 405 | swingFreq=abs(*itmin -minFreq); |
404 406 | if(swingFreq < abs(*itmax -maxFreq)) |
405 407 | swingFreq=abs(*itmax -maxFreq); |
406 408 | if(firstchanshift < abs(*itf-minfirstchan)) |
407 409 | firstchanshift=abs(*itf-minfirstchan); |
408 410 | itf++; |
409 411 | itmax++; |
410 412 | } |
411 413 | Int extrapad=max(min(4, Int(imNChan/10)),1); |
412 414 | swingpad=2*(Int(std::ceil((swingFreq+firstchanshift)/freqincr))+extrapad); |
413 - | //cerr <<" swingfreq " << (swingFreq/freqincr) << " firstchanshift " << (firstchanshift/freqincr) << " SWINGPAD " << swingpad << endl; |
415 + | cerr <<" swingfreq " << (swingFreq/freqincr) << " firstchanshift " << (firstchanshift/freqincr) << " SWINGPAD " << swingpad << endl; |
414 416 | //////////////// |
415 417 | return swingpad; |
416 418 | |
417 419 | } |
418 420 | |
419 421 | String BriggsCubeWeightor::makeScratchImagingWeightTable(CountedPtr<Table>& weightTable, const String& filetag){ |
420 422 | |
421 423 | //String wgtname=File::newUniqueName(".", "IMAGING_WEIGHT").absoluteName(); |
422 424 | String wgtname=Path("IMAGING_WEIGHT_"+filetag).absoluteName(); |
423 425 | //cerr << "NAME " << wgtname << endl; |