lofar-astron / factor

Facet calibration for LOFAR
http://www.astron.nl/citt/facet-doc
GNU General Public License v2.0
19 stars 12 forks source link

polarization signal effecting DDE selfcal #49

Open rvweeren opened 8 years ago

rvweeren commented 8 years ago

Discovered an "interesting" issue about polarized signal in the data (I thought about it before but never looked into it, until now)

For very bright sources I think we need to also calibrate the slow 0:1 and 1:0 gains.

See attached example of a low-res facet images (made manually, identical brightness scaling). The red circle is 3C295 (70 Jy), this source is quite nicely removed in the DDE selfcal before. However, it only removes the source from the XX and YY correlations (so Stokes I and Stokes Q, as it only calibrates Gain 1:1 and 0:0) and a significant amount of flux remains in the U,V images (and XY and YX correlations). In XY and XY 3C295 flux still dominates over all other sources in the FoV. (and even over XX and YY for the source in the green facet)

The image was made for the facet with the green circle (a ~1 Jy source). The Q, V signal from 3C295 still dominates over the Stokes I flux in the green circle. Since BBS minimizes the -full- visibility correlation set (during solve) our solutions should be affected by the strong "spurious" signal in XY and YX

Interested in hearing some opinions...

Solving for the slow 0:1 and 1:0 Gains for extremely bright sources in the DDE selfcal, typically not more than 1-3 sources per field I think, should solve this (and we assume Q, U are zero). As far as I can see only a few modification are needed in the code. (e.g., the parmdbs need to be created including these entries). So if a source is above ~5 Jy or so we can trigger the 0:1 and 1:0 Gains (optionally the user can change this threshold). This will likely need my new spline outlier rejection smoothing code (as the average gains for 0:1 and 1:0 are not 1.0).

This will likely also benefit the overall polarization effort as well (as at least the brightest sources are removed, Vibor also did that for his work, otherwise he could not study the fainter polarized diffuse galactic emission)

stokesiquv
rvweeren commented 8 years ago

To put some more context, I added the facet layout in this image

stokesiquv2
darafferty commented 8 years ago

I've now added an option to the reimage branch to solve for all 4 correlations during the slow gain solve. It's called "solve_all_correlations_flux_Jy" (under [calibration]) and sets the flux level above which 0:1 and 1:0 will be solve for (default is 5 Jy). The "spline_smooth2D" option must be enabled as well.

rvweeren commented 8 years ago

It does not (yet) like the 0:1 and 1:0. Things proceed normally until it tries the amplitude solve with all polarizations.

I think I spot the mistake. Below is the parset. The "parset" key is incorrect as the file does not exists. As you can see it pastes the full directory path twice! (/usr/local/lib/python2.7/dist-packages/FACTOR-1.0-py2.7.egg/factor/parsets/ needs to be stripped from the name of the parset)

Step.solve.Solve.CellSize.Freq=1 Step.solve.Solve.CellSize.Time=1 Step.solve.Solve.UVRange=[80.0] parset=/usr/local/lib/python2.7/dist-packages/FACTOR-1.0-py2.7.egg/factor/parsets//usr/local/lib/python2.7/dist-packages/FACTOR-1.0-py2.7.egg/factor/parsets/facet_dirdep_amponly_solve_allcorr.parset

When running that BBS simply crashes (as it cannot find the parset)

mhardcastle commented 8 years ago

Fixed the problem in facet_ops.py, I believe.

rvweeren commented 8 years ago

Ah great, I was looking where the problem occurred but didn't manage to find it.

mhardcastle commented 8 years ago

.... but for my dataset it fails immediately after the calibration with error messages about not being able to parse NaN values.

rvweeren commented 8 years ago

For my reduction your fix works and it is proceeding without crashing so far. (successfully solved, smoothed, applied and cleaned)

mhardcastle commented 8 years ago

OK, I'll pursue this further tomorrow.

rvweeren commented 8 years ago

Example for the Gain:0:1 solutions for the RS stations. They have pretty good S/N (well this is 3C295 at 70 Jy). I noticed that the gains are higher for the RS stations compared to the CS and are highest for the most distant stations. There is a very clear dependence with distant to the core. This must be ionosphere related….The phases for 0:1 are rather constant. Trying to understand why it looks like this… (dFaraday Rotation maybe, although I would have expected that to be much smaller…hmm…more thinking required). Anyway the four correlation solve part does seems to do its job in factor.

01sols
rvweeren commented 8 years ago

Found out that the "merge_selfcal_parmdbs" gain 0:1 and gain 1:0 values are wrong (very weird values that do not make much sense). Probably a bug in the code that creates merge_selfcal_parmdbs. So I retract my statement that things work as expected (maybe this is related to the problem Martin has).

I first thought it was the spline smoothing, but I tested that and the smoothing worked and the results look good. Also merge_amp_parmdbs2 and merge_amp_parmdbs1 are fine. So the problems seems to occur when it tries to merge merge_amp_parmdbs2/1 and merge_phase_parmdbs into merge_selfcal_parmdbs.

darafferty commented 8 years ago

Strange, I don't see any reason why the merging shouldn't work. The merge script is very generic, with nothing specific to the 0:0 or 1:1 gains.

rvweeren commented 8 years ago

Hmm, very weird. I tested the merge script standalone and it indeed works. However, I still have this wrong merge_selfcal_parmdbs, see image (merge_selfcal_parmdbs on top, merge_amp_parmdbs2 at the bottom). Any other code modifies merge_selfcal_parmdbs ?

parmdbs
darafferty commented 8 years ago

Does the smooth_amp2 parmdb look OK? The merge_selfcal_parmdbs script merges smooth_amp2 with merge_phase_parmdbs.

rvweeren commented 8 years ago

Yes. (& see above plot which shows Gain 0:1 for smooth_amp2).

And when I merge smooth_amp2 with merge_phase_parmdbs the output looks good.

rvweeren commented 8 years ago

Maybe I should start from scratch and see if I hit the problem again if there is no other factor code that touches/modifies merge_phase_parmdbs ?

darafferty commented 8 years ago

No, the parmdb should not be modified after merging.

rvweeren commented 8 years ago

Ah I see where it goes wrong. The normalization is not correct when we have four correlations smooth_amps_spline.py.

Clip extremely low amplitude solutions to prevent very high

                # amplitudes in the corrected data
                low_ind = numpy.where(amp < 0.2)
                amp[low_ind] = 0.2
rvweeren commented 8 years ago

I think this should be modified to

                if pol == '1:1' or pol == '0:0':
                  low_ind = numpy.where(amp < 0.2)
                  amp[low_ind] = 0.2
darafferty commented 8 years ago

Oh, yeah - that would do it. Perhaps this check isn't needed anyway? (I think I added it during some of my testing when I had noisy solutions, but I think they were noisy due to poor choice of solution interval, too much flagged data, etc., all of which are handled better now.)

rvweeren commented 8 years ago

Maybe we should indeed remove this. The reason is that if such values still occur either the user is doing something wrong or there is really bad data present. In that case it is maybe better not to "hide" the fact that there are amps with values < 0.2. Maybe we should throw a warning to the user if it encounters such values.

rvweeren commented 8 years ago

I am working on updating plot_selfcal_solutions.py to handle 4 polarizations. Amplitudes are done, still working on the phases.

test_amp_channel0

rvweeren commented 8 years ago

And the phases.

test_phase_channel3

BTW Is there a way for me to push this updated script to the reimage branch ? "git push origin reimage" gives me a permission denied error...

mhardcastle commented 8 years ago

If you have selected the reimage branch with 'git checkout -b reimage' then later 'git push' should just work. (It did for me the other day.)

rvweeren commented 8 years ago

Hmmm, does someone need to give me read/write access to the repository ? (git checkout -b reimage did not fix the problem)

remote: Permission to lofar-astron/factor.git denied to user. fatal: unable to access 'https://github.com/lofar-astron/factor.git/': The requested URL returned error: 403

AHorneffer commented 8 years ago

Yes, you need write permission to push directly to the repository. Looking at the existing list of collaborators I guess it is O.K. to add you too, so that's what I just did.

Or you do a fork of the repository to your own account, push the new code there, and then do a pull request on GitHub.

darafferty commented 8 years ago

I've added the arguments to the pipeline steps to activate the 4-pol plotting when needed.

rvweeren commented 8 years ago

I think it is working properly. I'll do one more test to double check.

rvweeren commented 8 years ago

I am still trying to confirm the Q,U,V images are more "empty" now. Should have results by tomorrow.

rvweeren commented 8 years ago

BTW I did some tests and I do not see improvements for the all four correlation solve (results change but do not clean up Q, U, V, not fully understand this yet). For now I suggest to put solve_all_correlations_flux_Jy to a very high value by default (999 or so). More testing and thinking is needed. I do not see obvious problems with factor itself. The code seems to do what it is asked to do.

This might be a more fundamental issue of how/what one solves. Maybe we should solve for leakage terms on the slow-fains timescales, and not for gain:0:1 and gain:1:0. It could be that gain:0:1 and gain:1:0 are inherently unstable in combination with a correct step. Also with the WSRT, JVLA, even at low-freqs one never solves for gain:0:1 and gain:1:0 and only leakage. (for example WSRT has done 150 MHz polarization, and again one solves for leakage in that case). The must have been a reason why one never solves for gain:1:0 and gain:0:1 in casa/aips/miriad…. (someone has the anwer?)

I can try testing solving for leakage manually (taking some of the factor output and running casa), but it will take quite some time to delve into this.

mhardcastle commented 8 years ago

On Tue, May 03, 2016 at 09:50:21AM -0700, Reinout van Weeren wrote:

Also with the WSRT, JVLA, even at low-freqs one never solves for gain:0:1 and gain:1:0 and only leakage. (for example WSRT has done 150 MHz polarization, and again one solves for leakage in that case). The must have been a reason why one never solves for gain:1:0 and gain:0:1 in casa/aips/miriada*|. (someone has the anwer?)

I think it's what you said: it's more stable to fit for a gain in 0:0 and 1:1 and a small leakage term than to fit for 0:1 and 1:0 directly. Especially if the leakage term can be parametrized as a constant per antenna.

darafferty commented 8 years ago

I've set the default to 1000.0 Jy.