I am interested in using the python sky model to simulate galactic foregrounds at high spatial resolution (~1 arcmin) for the preparation of upcoming sub-mm surveys. This would require to upgrade the maps to a healpix nside of 4096 or even 8192. The PySM paper shows that small scale features can be added to the map by extrapolation of the power spectrum to high ell and by modifying the synfast proposal distribution for generating random a_lm. However, it seems like this part of the code is not public at the moment. Would it be feasible for you to share the code for adding small scale features to the maps?

Hi,
I've created the following instrument:

instrument_noisy = pysm.Instrument({
                     'nside': d['nside'],
                     'frequencies' : np.array([143]),
                     'use_smoothing' : False,
                     'beams' : None, 
                     'add_noise' : True,
                     'noise_seed' : 0,
                     'sens_I': np.array([15]), 
                     'sens_P': np.array([30]), 
                     'use_bandpass' : False, 
                     'channel_names' : np.ones(1),
                     'channels' : None,
                     'output_units' : 'uK_RJ',
                     'output_directory' : "./",
                     'output_prefix' : "noisy",
                     'pixel_indices' : None})

Then, I observe my sky with:

instrument_noisy.observe(pysm.Sky(sky_config))

and I read the maps:

hp.read_map('noisy_nu0143p00GHz_total_nside0128.fits', field=(0,1,2))
array([[-31.023283  ,  11.80905437,  21.81184769, ...,  76.90827179,
         82.10586548,  70.3621521 ],
      [  2.53387022,   1.2787329 ,   2.67105889, ...,  -1.31416452,
        -0.67028081,   2.02661443],
      [  0.53173047,   0.92604077,  -0.27713001, ...,  -0.65416229,
        -2.46461463,   0.93581396]])

hp.read_map('noisy_nu0143p00GHz_noise_nside0128.fits', field=(0,1,2))
array([[ 0.58221883,  0.13207038,  0.3230288 , ..., -0.03905814,
        -0.37743953, -0.55901247],    
       [ 1.17966449, -1.23902667,  1.01518738, ..., -1.45075607,
        -0.27817139,  0.76594335],
       [-0.36754149, -0.32368761, -0.06424909, ...,  0.80313766,
        -0.30539998,  0.36780828]])

Now, I do the same with the flag write_outputs=False and I expected to obtain the same maps. Instead, while the noise maps are the same, the total maps are not since the last one doesn't contain noise:

total, noise = instrument_noisy.observe(pysm.Sky(sky_config), write_outputs=False)
total
array([[[-31.60550231,  11.67698437,  21.4888182 , ...,  76.94732776,
          82.48330691,  70.92116465],
       [  1.35420567,   2.51775958,   1.65587158, ...,   0.13659159,
         -0.39210943,   1.26067109],
       [  0.89927195,   1.24972841,  -0.21288092, ...,  -1.45729993,
         -2.15921475,   0.56800566]]])

So, if you add the noise:

total+noise
array([[[-31.02328346,  11.80905476,  21.811847  , ...,  76.90826962,
          82.10586738,  70.36215217],
       [  2.53387016,   1.27873293,   2.67105895, ...,  -1.31416448,
         -0.67028081,   2.02661443],
       [  0.53173045,   0.92604079,  -0.27713001, ...,  -0.6541623 ,
         -2.46461474,   0.93581394]]])

It might be possible to make the output maps independet on the value of the flag write_outputs or, at least, to update the documentation of the method instrument.observe?

I cannot find details about the license of the code, is it open source? (https://github.com/blog/1530-choosing-an-open-source-license)

Hi,
I am confused about the units of the noise levels in Instrument.
The doc says that they are expected to be in uK_RJarcmin
https://github.com/bthorne93/PySM_public/blob/b3e23bb8f88182bf440b731a2bab7100a255af0d/pysm/pysm.py#L147

However, from here
https://github.com/bthorne93/PySM_public/blob/b3e23bb8f88182bf440b731a2bab7100a255af0d/pysm/pysm.py#L475
and
https://github.com/bthorne93/PySM_public/blob/b3e23bb8f88182bf440b731a2bab7100a255af0d/pysm/pysm.py#L532
I understand that the noise levels are expected in uK_CMBarcmin.

If you confirm that sens_I and sens_P are actually expected in CMB units I can open a pull request that changes the doc of Instrument
Thanks

@zonca These are the changes I think will need to be made to allow the code to consider only partial maps. Most likely an incomplete list, but I hope it helps! Thanks so much!

• Initialization:
  • Any use of the code is initialized using the models function in nominal.py. I imagine just adding a kwarg for pixel_indices to that would be the thing to do.
• Inputs:
  • Input templates are read-in in nominal.py. Modifying the models function in this file will do a lot of the work.
• Processing: Several models generate maps of random values on the fly (e.g. CMB, HD17). These will need to be checked to make sure they are either generated for the required number of pixels when possible, or for the whole sky and then reduced for e.g. the CMB.
  • components.py:
    • CMB class, taylens function.
    • Dust class, hensley_draine_2017 function through the draw_uval function, which draws a map of random values. Unfortunately, due to a work around we had to use in order to speed it up this model also relies on the initialise_hd_dust_model_bandpass function in pysm.py, which initializes the model when bandpass-integration is requested.
  • instrument.py:
    • noiser method of instrument class.
• Outputs:
  • Maps are only saved to file when the Instrument.observe is used on an instance of the Sky object. So Instrument.writer will have to be able to pad the chosen pixels with zeros when doing only partial skies.

Need to write unit tests for bandpass integrated noise unit conversion.

Hi!
I'm starting to learn how to use PySM. I just created a sky object:

sky_config = {
'dust' : models('d1', nside),
'synchrotron' : models('s1', nside),
 'ame' : models('a1', nside),
'freefree' : models('f1', nside),
'cmb' : models('c1', nside)
}

and then:

sky = pysm.Sky(skyconfig)
T, Q, U = sky.signal(nu=43)

it's probably a trivial question, but it's not clear to me which are the units of T, Q and U.
It's u_K_RJ? or u_K_CMB? or something else?
I read the documentation but I couldn't find this information explicitly reported anywhere.

Thank you in advance!

Silvia

As there are no current Pull Requests open I think it would be a good time to tag a release on Github, maybe 2.1.0 (as I think the first version of this new codebase was tagged 2.0).

It would be useful for our work on TOAST so we can install a specific release.

See https://help.github.com/articles/creating-releases/

Hi, I am wondering if there is a method to calculate the E and B mode signal from Q and U map in the pysm?

K_RJ bandpass unit conversion seems to be causing a ~10% depression, even when simulating delta bandpasses with very narrow frequency ranges.

Hello,
I have a problem using pysm.Instrument.

I use this standard sky configuration:

nside = 128
sky_config = {
'synchrotron': models('s1', nside),
'dust': models('d1', nside),
'freefree': models('f1', nside),
'cmb': models('c1', nside),
'ame': models('a1', nside)} 

then I create a sky with:

sky = pysm.Sky(sky_config)

Now, I create the two following different instruments:

instrument_mono_noiseless = pysm.Instrument({
    'nside': nside,
    'frequencies' : np.array([143]), 
    'use_smoothing' : True,
    'beams' : np.array([7.1]),  
    'add_noise' : False, 
    'noise_seed' : None,  
    'sens_I': None, 
    'sens_P': None, 
    'use_bandpass' : False,  
    'channel_names' : ['143_GHz'],
    'channels' : None,
    'output_units' : 'uK_RJ',
    'output_directory' : "./",
    'output_prefix' : "noiseless_mono",
    'pixel_indices' : None})

instrument_mono = pysm.Instrument({
    'nside': nside,
    'frequencies' : np.array([143]), 
    'use_smoothing' : True,
    'beams' : np.array([7.1]),  
    'add_noise' : True, 
    'noise_seed' : 0,  
    'sens_I': np.ones(1), 
    'sens_P': np.ones(1), 
    'use_bandpass' : False,  
    'channel_names' : ['143_GHz'],
    'channels' : None,
    'output_units' : 'uK_RJ',
    'output_directory' : "./",
    'output_prefix' : "mono",
    'pixel_indices' : None})

I observe the sky through these two instruments with instrument_*.observe(sky) and then I read the corresponding maps that have been produced.

In the first case the created map does not have the polarization components but only the intensity one. In fact, you will get an error message if you try:

hp.read_map('noiseless_mono_nu0143p00GHz_total_nside0128.fits', field=(0, 1, 2))

while the right command would be:

hp.read_map('noiseless_mono_nu0143p00GHz_total_nside0128.fits', field=(0))

but this means that Q and U maps have not been generated.

In the second case, the three maps (I, Q, U) are generated but in a wrong way since the Q and U maps are in the same order of magnitude of the I map:

hp.read_map('mono_nu0143p00GHz_total_nside0128.fits', field=(0, 1, 2))

array([[-36.27698898,  14.9692564 ,  25.86756897, ...,  76.63456726, 
         81.74451447,  67.1904068 ],
       [-36.27647781,  14.91915131,  25.87987328, ...,  76.58881378, 
         81.76040649,  67.25320435],
       [-36.32805252,  14.94966221,  25.84389114, ...,  76.66394043, 
         81.7594986 ,  67.2399292 ]])

Is there something I'm doing wrongly somewhere?

Hello,
I think there is something misleading in the documentation of pysm.Instrument. Here in fact, it is written that sens_I and sens_P should be expressed in uK_RJ arcmin while in pysm.Instrument.noiser they are expected in uK_CMB arcmin. Is that right? This should be fix in the documentation, otherwise it could lead to underestimate the noise levels in the maps.

Thank you!

Hi ,
i 've just updated my current version of PySM and found that 2 Nosetests fail cause file not found , i.e.
test_data/test_nu0023p00GHz_noise_nside1024.fits

I'm trying to use your code to simulate galactic foreground around 1.4 GHz, and I find that most your model components are high frequency, except the 408 MHz all sky survey. So is it OK to use your code for 1.4 GHz?

hi!
the current 2.1.0 release tgz file at https://github.com/bthorne93/PySM_public/releases is 380MB.
This is quite big and is causing some random unit tests of TOAST to fail because the download from Github takes too long.

I wonder if it would be a good idea to store the maps in a separate Github repository and then use astropy.utils.data functionality to download it only when it is necessary (then cache it).

I recently implemented something similar in healpy to download pixel weights when needed: https://github.com/healpy/healpy/pull/442/files

Cell 13 in the example notebook (https://github.com/bthorne93/PySM_public/blob/master/example.ipynb) seems to suggest that one could make all spectral parameters of a given component scalars instead of maps. This however seems to work only for spectral indices, not for e.g. dust temperature (the code will crash if you replace 'spectral_index' with 'temp' in that same cell)

I also think some aspects of the example notebooks are outdated (e.g. bare 'print' statements a la python2.7)

Calculation of the unit conversion factor does not normalise the bandpass, however calculation of the integrated map does. Therefore there will be a normalization mismatch.

@rainwoodman recommended to run tests with runtests, which has the neat feature of supporting MPI tests.

@rainwoodman:

• does runtests allow to run mixed MPI and non-MPI tests?
• is the test collection different than py.test? I tried to vendor run-tests.py in the root folder of the package, rename the tested package and run it but it doesn't find any test, py.test instead finds them.