Polarisation imaging

[twshimwell]

Create Q and U cubes with 480 channels about 3 degrees across. This is ~7000 pixels with 1.5arcsec cellsize (current high imaging cellsize) or ~2400 pixels with 4.5arcsec cellsize.

Testing so far has shown that ~1000 pixels is fine with the memory constraints but its unclear what the maximum is.

It would be also nice to produce a stokes V image. My thoughts would be that this can be done with the 4.5arcsec cellsize at 20arcsec resolution.

[twshimwell]

A full bandwidth low resolution Stokes V image takes about 1hr. There is some structure in V around some sources (attached -- V is left I is right), is this to be expected do you think?

We have 2 pulsars in the HETDEX region so ill try those fields and see if we detect either of those.

[sposullivan]

Yes the structure in V gives an indication of the instrumental polarisation. Marco Iacobelli mentioned yesterday that he was interested in imaging Stokes V across the HETDEX area to get a handle on the instrumental polarisation variation across the beam. I think he was planning to use wsclean on the prefactor data, so maybe you can chat to him (he's not here the rest of the week I believe).

[sposullivan]

Cameron Van Eck doesn't mention imaging Stokes V in his paper but it might be worth checking with him in case they are just sitting on a disk somewhere.

[twshimwell]

I can redo them all. Its just an hour per field processing if I just do this step (so no DD). The main time will be downloading and unpacking (especially as SARA is down at the moment).

[twshimwell]

OK and for the QU cubes I chatted a little with @sposullivan and given there is no deconvolution in QU the resolved RM structures might be hard to accurately measure due to sidelobes. So, perhaps for DR2 processing we stick with low resolution imaging of QU even though we will probably loose some sources due to the low resolution.

For the full field of view we need ~2500 pixels but trying this on 256GB ram caused a memory error. 2000 pixels also crashed. 1000 works.

I guess what we will have to do is do 4x1250 pixel images with different phase centers.

I can put this into the pipeline if people are happy with this idea

[cyriltasse]

sounds excelent to me...

[twshimwell]

OK ive changed pipeline.py and created do_polcubes.py so that now we should get polarisation images (4 cubes of 480 channels and 1250x1250 pixels distributed to cover ~3x3 deg). An example polarisation cube is on https://www.dropbox.com/s/rpqg7ujli9wdyls/image_full_low_QU_Cube0.cube.dirty.corr.fits?dl=0 if anyone wants to check this looks ok.

[mhardcastle]

Made this an option -- not everyone using the pipeline will want it.

[mhardcastle]

Tim, wouldn't it be better to make the cubes over the full field but with fewer channels and then stitch them together in frequency? Then if we get QU deconvolution at some point in the future we can deconvolve them properly.

[twshimwell]

Yeh good point. Maybe that would be better indeed. I can look at changing it on Monday.

[twshimwell]

OK ive pushed a change that now does this 1/4 of the bandwidth at a time.

I guess we also now need something to glue the cubes together.

[twshimwell]

Imaging the 2nd part of the bandwidth I get an error in DDFacet.

Traceback (most recent call last): File "/home/shimwell/.local/bin/DDF.py", line 6, in exec(compile(open(file).read(), file, 'exec')) File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/DDF.py", line 387, in main(OP, messages) File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/DDF.py", line 228, in main Imager.Init() File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Imager/ClassDeconvMachine.py", line 169, in Init GD=self.GD) File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Data/ClassVisServer.py", line 91, in init self.Init() File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Data/ClassVisServer.py", line 137, in Init first_ms=self.ListMS[0] if self.ListMS else None) File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Data/ClassMS.py", line 120, in init self.ReadMSInfo(first_ms=first_ms,DoPrint=DoPrint) File "/net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Data/ClassMS.py", line 1040, in ReadMSInfo self.dFreq = ta_spectral.getcol("CHAN_WIDTH")[self._spwid,self.ChanSlice].flatten()[0] IndexError: index 0 is out of bounds for axis 0 with size 0

• 13:25:19 - ClearSHM | Clear shared memory
• 13:25:19 - ClearSHM | Clear Semaphores
• 13:25:19 - ClearSHM | Clear shared dictionaries FAILED to run DDF.py --Output-Name=image_full_low_QU_Cube1 --Data-MS=big-mslist.txt --Deconv-PeakFactor 0.001000 --Data-ColName DATA_DI_CORRECTED --Parallel-NCPU=32 --Beam-CenterNorm=1 --Deconv-CycleFactor=0 --Deconv-MaxMinorIter=1000000 --Deconv-MaxMajorIter=0 --Deconv-Mode SSD --Beam-Model=LOFAR --Beam-LOFARBeamMode=A --Weight-Robust -0.250000 --Image-NPix=2500 --CF-wmax 50000 --CF-Nw 100 --Output-Also onNeds --Image-Cell 4.500000 --Facets-NFacets=11 --SSDClean-NEnlargeData 0 --Freq-NDegridBand 1 --Beam-NBand 1 --Facets-DiamMax 1.5 --Facets-DiamMin 0.1 --Deconv-RMSFactor=3.000000 --SSDClean-ConvFFTSwitch 10000 --Data-Sort 1 --Cache-Dir=. --Log-Memory 1 --GAClean-RMSFactorInitHMP 1.000000 --GAClean-MaxMinorIterInitHMP 10000.000000 --GAClean-AllowNegativeInitHMP True --DDESolutions-SolsDir=SOLSDIR --Cache-Weight=reset --Output-Mode=Clean --Output-RestoringBeam 20.000000 --Weight-ColName="IMAGING_WEIGHT" --Output-Cubes=dD --RIME-PolMode=QU --Output-Mode=Dirty --Freq-NBand=120 --Selection-ChanStart=120 --Selection-ChanEnd=240 --RIME-DecorrMode=FT --SSDClean-SSDSolvePars [S,Alpha] --SSDClean-BICFactor 0 --Mask-Auto=1 --Mask-SigTh=5.00 --DDESolutions-GlobalNorm=None --DDESolutions-DDModeGrid=AP --DDESolutions-DDModeDeGrid=AP --DDESolutions-DDSols=DDS3_full_smoothed --Selection-UVRangeKm=[0.100000,25.750000] --GAClean-MinSizeInit=10 --Beam-Smooth=1: return value is 1

Any ideas?

[cyriltasse]

These keywords are for the individual MSs...

[twshimwell]

Ah yes so it breaks when ChanStart falls outside the first msfile (so from 20 onwards). Is this a DDFacet bug or feature?

[twshimwell]

I've pushed changes that add in StokesV imaging as a option. This also uses the low resolution imaging parameters.

[mhardcastle]

Why don't you do it N MS at a time, then you don't need any selection at all?

I can write something to cubify it later.

[cyriltasse]

That's also what I do to do the shifted restored cubes...

[mhardcastle]

So you could copy Cyril's code, even better (-:

[twshimwell]

Ok, now one MS at a time. Each low resolution QU cube of 20 chans takes about 10mins.

[mhardcastle]

Or 5 mins on a proper computer (-: .

I'll write something to put these together.

[twshimwell]

Wonder how long a high res one would take on a proper computer.

[mhardcastle]

Guess we should see how much use the low-res ones are and then think again.

Cube-making code added to polcube and pushed.

[twshimwell]

Made a stokes-v map with and without the DD corrections to show the difference we are getting. Left with DD (so in the pipeline) and right without DD. Certainly cleans it up a lot and the noise is about 0.08mJy/beam. I'm not sure if we would expect DD to also do something to real stokes V sources but maybe that would require some simulations as im not sure we have any real stokes V sources to use for testing.

[mhardcastle]

Do you mean left with, right without?

[twshimwell]

Ah yes I did indeed mean that.

[cyriltasse]

You've scared me...

[twshimwell]

I've put a cube on ftp://ftp.strw.leidenuniv.nl/pub/shimwell/image_full_low_QU.cube.dirty.corr.fits if you want to check it looks ok @sposullivan

A few changes are still needed to make this work for observations observed multiple times but then this issue is done assuming the cube looks ok.

[sposullivan]

Thanks @twshimwell, downloading now. I will check this tomorrow, and report back. So neither of the two pulsars have Stokes V?

[twshimwell]

No neither yet have obvious excess Stokes V emission but we have only looked in the DI Stokes-V images (DD are not available yet) and we still need to write something that does the Stokes-V source finding and calculates the fractional stokes-V to look at the fractional polarisation. So its not very thorough yet.

[sposullivan]

Nice results from a known polarised source in the field (from image_full_low_QU.cube.dirty.corr.fits). Same RM as from WSClean imaging of prefactor data: I had to make a sub-cube around this source because I couldn't load entire cube into memory for RM synthesis. I will try again soon on computer with more memory (256 GB). Otherwise will split into smaller cubes for RM synthesis.

[sposullivan]

@twshimwell I found 9 (real) polarised sources in total in that cube, so that would imply a source density ~6 times higher than in Van Eck+18. The rms noise in Q and U is ~84 uJy/beam across the centre of the field and increases to ~100 uJy/beam in some of the noisier facets. The green circle is 2.5 deg diameter. As you can see, leakage from the brighter sources can still be problematic... It takes ~2 hours to do the RM synthesis and RMclean (if I use a mask, only including pixels with Stokes I > 1 mJy/beam)

[cyriltasse]

Would be ice to cross check if these sources are also seen in adjacent pointings... just to be sure it's not some weird effect. If indeed much higher it's great... 2 hours is really not very much for the pipeline, if you want to run RMclean within.

[sposullivan]

Going by the wiki, there seems to be a large overlap with P206+50/340810 so we could cross-check with that. It would be great to include RM synthesis/clean if feasible. Have you played around with pyrmsynth (https://github.com/mrbell/pyrmsynth) @cyriltasse? It has some issue reading the QU cube so I write out image planes, then remove the noisy/blank ones and run pyrmsynth using a mask (where the masked values in the fits file need to be zeros).

[cyriltasse]

No i did not have the chance to do a run yet. Not sure I will if it only takes a few hours it's very feasible this way...

[mhardcastle]

If we want to put pyrmsynth in the pipeline then adapting it to be able to read in our cubes should not be too hard.

Is the reason the facets show up in @sposullivan 's images that we don't use a smooth beam for the QU imaging?

[twshimwell]

Yeh I was just checking the smooth-beam thing too. Very odd noise structure.

These were run with beam-smooth=1 so im a bit surprised to see the facets like that.

[twshimwell]

Similarly the Stokes V images dont seem to have the smooth beam either. The command for the stokes V is. Does smooth-beam work for imaging polarisation?

DDF.py --Output-Name=image_full_low_stokesV --Data-MS=big-mslist.txt --Deconv-PeakFactor 0.001000 --Data-ColName DATA_DI_CORRECTED --Parallel-NCPU=32 --Beam-CenterNorm=1 --Deconv-CycleFactor=0 --Deconv-MaxMinorIter=1000000 --Deconv-MaxMajorIter=0 --Deconv-Mode SSD --Beam-Model=LOFAR --Beam-LOFARBeamMode=A --Weight-Robust -0.250000 --Image-NPix=6000 --CF-wmax 50000 --CF-Nw 100 --Output-Also onNeds --Image-Cell 4.500000 --Facets-NFacets=11 --SSDClean-NEnlargeData 0 --Freq-NDegridBand 1 --Beam-NBand 1 --Facets-DiamMax 1.5 --Facets-DiamMin 0.1 --Deconv-RMSFactor=3.000000 --SSDClean-ConvFFTSwitch 10000 --Data-Sort 1 --Cache-Dir=. --Log-Memory 1 --GAClean-RMSFactorInitHMP 1.000000 --GAClean-MaxMinorIterInitHMP 10000.000000 --GAClean-AllowNegativeInitHMP True --DDESolutions-SolsDir=SOLSDIR --Cache-Weight=reset --Output-Mode=Clean --Output-RestoringBeam 20.000000 --Weight-ColName="IMAGING_WEIGHT" --Freq-NBand=2 --RIME-PolMode=IV --Output-Mode=Dirty --RIME-DecorrMode=FT --SSDClean-SSDSolvePars [S,Alpha] --SSDClean-BICFactor 0 --Mask-Auto=1 --Mask-SigTh=5.00 --DDESolutions-GlobalNorm=None --DDESolutions-DDModeGrid=AP --DDESolutions-DDModeDeGrid=AP --DDESolutions-DDSols=DDS3_full_smoothed --Selection-UVRangeKm=[0.100000,25.750000] --GAClean-MinSizeInit=10 --Beam-Smooth=1

[twshimwell]

Hmm weird images when attempting to make very low resolution cubes... not seen this type of behavior before. The bottom right bit of the image is all nans...

[cyriltasse]

I've implemented a mask to remove the interesting science.... Can you redo in I-stokes?

[twshimwell]

Remaking in stokes I (by simply changing RIME-PolMode=QU to =I) shows the same science removing mask.

During these runs I see this in the log

================ Making Dirty Image and/or PSF ==================== /net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Imager/ClassFacetMachine.py:719: RuntimeWarning: invalid value encountered in divide sw /= np.max(sw) /net/para10/data1/shimwell/software/killmsddf/new-install/DDFacet/DDFacet/Imager/ClassFacetMachine.py:721: RuntimeWarning: invalid value encountered in less

[cyriltasse]

Can you look at you .tessel.reg?

[twshimwell]

Looks normal to me

Ive put an example log on https://www.dropbox.com/s/pkkmpdozaw25t3m/DDF-image_full_vlow_QU_Cube12.log?dl=0

[cyriltasse]

F21_S32 look devil

[cyriltasse]

Can you pull again on the DDFacet active branch (ddf_pipeline_fixes_may18) and try again? Perhaps remove your *.ddfcache before running

[cyriltasse]

Hold on - this might make things bad

[cyriltasse]

can you just try to increase --Facets-DiamMin = 0.2 for example

[twshimwell]

Changing that parameter didnt fix it. However, just doubling the image size did so I think its all fine.

[twshimwell]

With these polarisation images the main thing left is to get the smooth beam. @cyriltasse did you have any ideas why the beam might not be smooth for these images?

[twshimwell]

Ive pushed changes to parset.py, do_polcubes.py and pipeline.py so that we now have very low resolution (uvmax 1.6km) polarisation imaging in the pipeline.

[twshimwell]

As Martin mentioned in #73 the largest products we are archiving are by far the QU cubes... Do we really need to store these? Can we not do rmsynthesis in the pipeline and then store smaller products? Or does that take a lot of time to run?