galsci / pysm Goto Github PK

Based on d10

d11 will have templates just at large scale saved at low nside, C_ell of small scales and map of modulation at low nside loaded from network.
synalm + almxfl + alm2map to generate templates on the fly.
Using a specific seed should reproduce d10 but on the fly.
If seed is not specified, the small scales are stochastic.

Apparently importlib_resources doesn't get installed in python 3.6 but seems to be required.
Is it something to be fixed in the install process or is it a problem on my side?

I recently updated my Mac OS to Catalina, and the first thing post-upgrade was getting my Python packages back up and running. Unfortunately, I was unable to get pysm3 to install via conda following the docs, though the pip install worked without issue. Perhaps that should have been a warning sign.

When I try to run the simple code:

import pysm3
import astropy.units as u

nside = 64
sky = pysm3.Sky(nside = nside, preset_strings = ['d0'])
test_map = sky.get_emission(250.*u.GHz)
smooth_map = pysm3.apply_smoothing_and_coord_transform(test_map, fwhm=5./60.*u.deg)

I get an ominous OMP error:

OMP: Error #15: Initializing libomp.dylib, but found libomp.dylib already initialized.
OMP: Hint This means that multiple copies of the OpenMP runtime have been linked into the program. That is dangerous, since it can degrade performance or cause incorrect results. The best thing to do is to ensure that only a single OpenMP runtime is linked into the process, e.g. by avoiding static linking of the OpenMP runtime in any library. As an unsafe, unsupported, undocumented workaround you can set the environment variable KMP_DUPLICATE_LIB_OK=TRUE to allow the program to continue to execute, but that may cause crashes or silently produce incorrect results. For more information, please see http://openmp.llvm.org/
Abort trap: 6

I will try the suggested workaround on my machine, but it's not an ideal fix. Any insights or suggestions would be greatly appreciated.

My most common case of use is the production of a multi-frequency observation with delta band passes. I have the frequencies of my experiment in an array of shape (n_freq) and I'd like to get an array of shape (n_freq, n_stokes, n_pixels) containing the observation.

From a quick look at the doc and the code, it seems to me that Sky.get_emission can only return (n_stokes, n_pixels): arrays of frequencies are processed as band-pass samples and integrated over.

Is there a way of getting all the frequencies (as in pysm2's instrument.observe(sky)) or I'm expected to iterate over the frequencies and stack the result in my application code?
Thanks!

Recently, an issue has been raised on the transformation of the intensity component within the Log Pol tens formalism . I address it here

below the maps how they look like

Finally I asserted that :

• E[i'] = E[i] +b: True up to a relative tolerance of 1e-7
• Var[i' ] =Var[i] : True up to a relative tolerance of 1e-4
• Cl^i'i' = Cl^ii : True up to a rel. tol of 1e-5

i would conclude that this issue shouldn't be worrisome for time being. what do you think ?

I am running an SO-specific TOAST simulation with this preset:

SO_d0,SO_s0,SO_a0,SO_f0,SO_x1_cib,SO_x1_tsz,SO_x1_ksz,SO_x1_cmb_lensed_solardipole

or in the case of LF1 band:

SO_d0,SO_s0,SO_a0,SO_f0,SO_x1_tsz,SO_x1_ksz,SO_x1_cmb_lensed_solardipole

The simulation resolution is nside=512.

In the LF1 case, simulating roughly 12 days of data at 37Hz for 65 detectors took about 9 minutes using 32 MPI tasks distributed over 32 nodes. Each MPI task was given 8 threads. Each node was running a total of 8 independent PySM calculations simultaneously.

In the MFS1 case, simulating same 12 days of data at same 37Hz for 43 detectors, using 16 processes on 16 nodes and 8 threads each did not finish in 2 hours. The test has only been run once.

Either the overall size of the jobs (32 processes vs. 256 processes, 518 threads vs. 5558 threads) causes a scaling issue, or the CIB component is very expensive to evaluate. I'll soon have numbers for an LF2 job that is nearly identical with the LF1 case but with CIB.

Hi,
@yaojian95 has been working a bit on poltens dust maps produced by @giuspugl in #97, in order to estimate non-Gaussianity. He has computed power spectra for different fsky. In particular, power spectra on full sky look a bit weird, with bumps around \ell~100 (see below). Have you already discussed this point? I'm not sure whether this is critical or not, as spectra look ok once we move away from the Galactic plane, but we just want to make sure that you are aware of this features.

When executing the following code:

import numpy as np
import pysm
freqs = np.linspace(10, 400, 100) * u.GHz
sky = pysm.Sky(nside=nside, preset_strings=['d1'])
emission_d1 = sky.get_emission(freqs)
print(emission_d1.shape)

I get an array of shape (3, npix), rather than (nfreqs, 3, npix). This is not the expected behavior according to the docstring. It looks like the code automatically defines a 'bandpass' with equal weight at each of these frequencies, and averages over them.

Using master branch, with a conda-forge python3.6 stack and running the toast unit tests I get:

Traceback (most recent call last):
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/tests/ops_sim_pysm.py", line 192, in test_op_pysm_nosmooth
    units="uK_RJ",
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/timing.py", line 43, in df
    return f(*args, **kwargs)
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/todmap/sim_det_pysm.py", line 141, in __init__
    map_dist=map_dist,
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/todmap/pysm.py", line 76, in __init__
    if init_sky
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/timing.py", line 43, in df
    return f(*args, **kwargs)
  File "/home/kisner/software/toast/lib/python3.6/site-packages/toast/todmap/pysm.py", line 124, in init_sky
    output_unit=self._units,
  File "/home/kisner/software/cmbenv/cmbenv_aux/lib/python3.6/site-packages/pysm-3.0.0-py3.6.egg/pysm/sky.py", line 101, in __init__
    component_config, nside=nside, map_dist=map_dist
  File "/home/kisner/software/cmbenv/cmbenv_aux/lib/python3.6/site-packages/pysm-3.0.0-py3.6.egg/pysm/sky.py", line 37, in create_components_from_config
    map_dist=map_dist
  File "/home/kisner/software/cmbenv/cmbenv_aux/lib/python3.6/site-packages/pysm-3.0.0-py3.6.egg/pysm/models/spdust.py", line 53, in __init__
    self.I_ref = self.read_map(map_I, unit=unit_I)
  File "/home/kisner/software/cmbenv/cmbenv_aux/lib/python3.6/site-packages/pysm-3.0.0-py3.6.egg/pysm/models/template.py", line 65, in read_map
    dataurl=self.dataurl,
  File "/home/kisner/software/cmbenv/cmbenv_aux/lib/python3.6/site-packages/pysm-3.0.0-py3.6.egg/pysm/models/template.py", line 200, in read_map
    filename = data.get_pkg_data_filename(path)
  File "/home/kisner/software/cmbenv/cmbenv_python/lib/python3.6/site-packages/astropy/utils/data.py", line 588, in get_pkg_data_filename
    timeout=remote_timeout)
  File "/home/kisner/software/cmbenv/cmbenv_python/lib/python3.6/site-packages/astropy/utils/data.py", line 1012, in download_file
    with shelve.open(urlmapfn) as url2hash:
  File "/home/kisner/software/cmbenv/cmbenv_python/lib/python3.6/shelve.py", line 243, in open
    return DbfilenameShelf(filename, flag, protocol, writeback)
  File "/home/kisner/software/cmbenv/cmbenv_python/lib/python3.6/shelve.py", line 227, in __init__
    Shelf.__init__(self, dbm.open(filename, flag), protocol, writeback)
  File "/home/kisner/software/cmbenv/cmbenv_python/lib/python3.6/dbm/__init__.py", line 91, in open
    "available".format(result))
dbm.error: db type is dbm.gnu, but the module is not available

Is the dbm format python version dependent? Installing the "gdbm" package did not help.

In the LiteBIRD Simulation Framework, we generate a set of CMB maps and pass them to PySM3, instead of using the maps produced by PySM3 itself (see the code here).

Since we generate the CMB maps in our code, they are saved outside the directory where PySM3 keeps its templates. However, pysm3.CMBMap expects relative paths rather than absolute. To let pysm3.CMBMap find our maps, our code overrides PYSM_LOCAL_DATA every time we need to include the CMB in the simulated sky. This is not optimal, as in our CI pipeline we need PYSM_LOCAL_DATA to point to the folder where we store a cache of the templates: obviously, this folder is always overwritten when we include any CMB signal in the sky!

I and @NicolettaK have thought about possible way to solve this:

1. Avoid including the CMB signal when running the tests. We have decided to follow this route, but we fear that the problem might appear again in the future, if somebody wants to run our pipeline on a server where PySM3 templates are supposed to be kept locally.

2. Save the value of PYSM_LOCAL_DATA before overriding it, and then restore the old value once the CMB sky has been generated. However, this is tricky to implement, if one wants to be robust when running several MPI processes.

3. Prepend/append the directory where we have saved our CMB maps to PYSM_LOCAL_DATA instead of overwriting the variable completely. This would work if this variable worked similarly to the UNIX PATH variable, i.e., a colon-separated list of directories like /storage/PySM3/templates:/litebird_CMB_maps. Unfortunately, it seems that this behavior is not supported, as pysm3/utils/data.py contains the following line:

    self.data_folders.append(os.environ["PYSM_LOCAL_DATA"])

   instead of

    self.data_folders += os.environ["PYSM_LOCAL_DATA"].split(os.pathsep)

Does any of this sound reasonable? Or are there other simpler and viable solutions we are missing?

Running PySM3 standalone using modifications of this script:

https://github.com/healpy/pysm/blob/master/mpi_examples/pysm3_mpi_so_pysm_models.py

I have all the 8 components requested for the TOD sims, which load a total of 30 maps. Plus bandpass integration of 5 points each for 10 channels and distributed smoothing with libsharp.

Based on the slides emailed on Oct. 19, @delabrou recommends we, "[m]ultiply dust template maps by 0.92 to get a template at 353 GHz rather than integrated in the Planck 353 GHz band." Note that a similar color correction factor appears in Planck 2018 XI, Table 2 as 1/1.098 = 0.911. I suggest we use this value so that we can cite the Planck table. @giuspugl

when using the <<= operator, need to add a with statement with the equivalency I want to use:

http://docs.astropy.org/en/stable/units/equivalencies.html#equivalency-context

Get the script in TOAST that checks memory usage across the MPI communicator: https://github.com/hpc4cmb/toast/blob/master/pipelines/toast_ground_sim.py#L54

@bthorne93 the pysm.units and from pysm.component_models.extragalactic.cib import InterpolatedCIB modules used in the tests are missing, can you please add them?

In #11 implemented MPI support and pixel distribution based on rings as expected by libsharp.
Now should do the actual distributed smoothing with libsharp.

Import the code from https://github.com/hpc4cmb/toast/blob/master/src/python/map/smooth.py

I updated to the latest vesion 3.2.2, but I cannot run a model that includes dust. For example, if I run this example

import pysm3
import pysm3.units as u
import healpy as hp
import numpy as np
import warnings
warnings.filterwarnings("ignore")
sky = pysm3.Sky(nside=128, preset_strings=["d1", "s1"])
map_100GHz = sky.get_emission(100 * u.GHz)

I get this error

Traceback (most recent call last):
  File "Test.py", line 11, in <module>
    map_100GHz = sky.get_emission(100 * u.GHz)
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/sky.py", line 151, in get_emission
    output = self.components[0].get_emission(freq, weights=weights)
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/astropy/units/decorators.py", line 234, in wrapper
    return_ = wrapped_function(*func_args, **func_kwargs)
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/models/dust.py", line 99, in get_emission
    outputs = get_emission_numba(
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/numba/core/dispatcher.py", line 415, in _compile_for_args
    error_rewrite(e, 'typing')
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/numba/core/dispatcher.py", line 358, in error_rewrite
    reraise(type(e), e, None)
  File "/home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/numba/core/utils.py", line 80, in reraise
    raise value.with_traceback(tb)
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython mode backend)
Failed in nopython mode pipeline (step: nopython frontend)
No implementation of function Function(<built-in function getitem>) found for signature:
 
 >>> getitem(array(float64, 1d, C), UniTuple(uint64 x 2))
 
There are 22 candidate implementations:
  - Of which 20 did not match due to:
  Overload of function 'getitem': File: <numerous>: Line N/A.
    With argument(s): '(array(float64, 1d, C), UniTuple(uint64 x 2))':
   No match.
  - Of which 2 did not match due to:
  Overload in function 'GetItemBuffer.generic': File: numba/core/typing/arraydecl.py: Line 162.
    With argument(s): '(array(float64, 1d, C), UniTuple(uint64 x 2))':
   Rejected as the implementation raised a specific error:
     TypeError: cannot index array(float64, 1d, C) with 2 indices: UniTuple(uint64 x 2)
  raised from /home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/numba/core/typing/arraydecl.py:84

During: typing of intrinsic-call at /home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/models/dust.py (137)

File "../../../Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/models/dust.py", line 137:
def get_emission_numba(
    <source elided>
        temp[I] *= (freq / freq_ref_I) ** (mbb_index - 2.0)
        temp[Q:] *= (freq / freq_ref_P) ** (mbb_index - 2.0)
        ^

During: lowering "id=8[LoopNest(index_variable = parfor_index.121, range = (0, $190binary_subscr.52_size0.50, 1)), LoopNest(index_variable = parfor_index.122, range = (0, mbb_index_size0.19, 1))]{306: <ir.Block at /home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/models/dust.py (137)>}Var($parfor_index_tuple_var.123, dust.py:137)" at /home/chervias/Software/anaconda3/envs/healpy/lib/python3.8/site-packages/pysm3/models/dust.py (137)

I believe it is a good idea, as it is now, to have the Model subclasses take care of the bandpass integration, see #9

However, I think it is convenient to have a Instrument class with a similar interface to PySM 2, it is a good way of grouping the instrument properties together that are otherwise just independent arrays like fwhms and bandpasses.

We can still provide the same observe method that calls get_integrated_emission of Sky and passes bandpasses.

Smoothing instead can just be implemented inside Instrument, I don't think there is any advantage in implementing it inside Model.

Finally Instrument will be able to generate white noise maps and also optionally write the outputs to disk. So it should process one channel at a time.

In order to run the Simons Observatory simulations, I would like to have all the models used by https://github.com/simonsobs/so_pysm_models/blob/master/so_pysm_models/so_models.py

In particular, modified_black_body and power_law we should already have. Are power_law and SynchrotronPowerLaw the same thing?

I need also spdust and spdust_pol for AME

@NicolettaK has an updated version of the notebook to generate synchrotron beta maps (apologies Nico if you've already posted this elsewhere!). A standalone implementation (i.e., with various parameters hard coded) begins in Cell 31. @giuspugl does this have everything you need to make an implementation that parallels yours for dust beta?

Originally posted by @brandonshensley in #73 (comment)

The dust code has properties and setters, check if we can simplify the code or they are really the best way to implement that code.

Current sky models at high resolution are very expensive to run, especially memory-wise.

A complete sky model with CMB, foregrounds and extragalactic can have 8-9 components, each needs to load a number of high resolution maps, generally 2 or 3, but sometimes 5 or 6.

Running PySM on the fly for time domain simulations is very expensive. Some pre-computation is most probably needed.
The simplest model would be to only simulate the impact of shifting the center of the bandpass, keeping bandwidth and shape the same for all detectors.

My proposal for this is:

• Choose a representative band shape, simplest is top-hat, but it would be better if we have something more realistic (i.e. lower towards the edges)
• Precompute with PySM complete sky maps densely around the center frequency. E.g. we are simulating 40GHz, the tolerance of center frequency is +-1 GHz, we create 1 PySM map which includes all components for 39GHz, 39.2, 39.4...40.8 41GHz, so ~11 maps.
• the time domain simulation tool (TOAST) doesn't run PySM, it loads the 11 maps, then for each detector draws a center frequency 39-40GHz and interpolates linearly the 2 closest maps

Next level:

• Also modulate bandwidth, based on center frequency - the PySM maps have a 30% bandwidth, so the 39GHz channel has smaller bandwidth
• Also modulate bandwidth, independently - we precompute PySM for 5 points in frequency and 4 different bandwidths, produce 20 maps

Trying to run LAT LF simulation at Nside=4096. There are 222 detectors. I can only allocate two MPI processes on each KNL node not to exhaust the memory during PySM call. I split the communicator into 256 groups, each with one node and two processes. There are two rank communicators with 256 MPI processes, each trying to simulate 111 detectors. The job did not complete in two hours. The medium and high frequency tubes have 10X as many detectors.

I'm having problems at installing pysm3 at nersc.
I'm using the pip install pysm3 command with the --prefix option and getting the following error:

ERROR: Command errored out with exit status 1: command: /usr/common/software/tensorflow/gpu-tensorflow/1.15.0-py37/bin/python -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-aftj8io3/pysm3/setup.py'"'"'; __file__='"'"'/tmp/pip-install-aftj8io3/pysm3/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' egg_info --egg-base /tmp/pip-install-aftj8io3/pysm3/pip-egg-info cwd: /tmp/pip-install-aftj8io3/pysm3/ Complete output (5 lines): Traceback (most recent call last): File "<string>", line 1, in <module> File "/tmp/pip-install-aftj8io3/pysm3/setup.py", line 51, in <module> from astropy_helpers.setup_helpers import ( ModuleNotFoundError: No module named 'astropy_helpers.setup_helpers' ---------------------------------------- ERROR: Command errored out with exit status 1: python setup.py egg_info Check the logs for full command output.

Having variable spectral parameters are of course one of the challenges for many foreground cleaning methodologies. personally, I am in favor of having another model with spatially variable beta and T_d , e.g. the ones from GNILC maps (see Planck 2016 ).

Originally posted by @giuspugl in #83 (comment)

I can adapt that code to the new PySM and then integrate it directly inside PySM

https://github.com/zonca/pysm_interpolating_component

Originally posted by @zonca in #4 (comment)

I would like to optimize how bandpass integration is handled.

However, even before that, I would like to optimize how emissions are calculated when we get an input array of frequency.

For example now in the synchrotron power law, we have a loop on the frequencies:

https://github.com/healpy/pysm/blob/584898cad071e591eeb74b633d946502341aa6d9/pysm/component_models/galactic/synchrotron.py#L59

First I'd like to replace those loops with vectorized operations with numpy. This should be a significant performance improvement.

Second, I would like each component model to have a get_integrated_emission that gets an array of frequency and an array of weights. This function should perform bandpass integration without having in memory all of the maps at each frequency. For example if you have a tophat with 10 points, in PySM 2 and also in the current implementation of PySM 3, you get 10 IQU maps for all the frequencies and just after perform bandpass integration. I'd like to skip a step and just accumulate a single IQU output map.

Problem

PySM needs to store the templates at very different resolutions, we want people to be able to do maps at low nside without the need of downloading and processing very high resolution maps.

Current implementation

Currently PySM 3 implements the same method of PySM 2.
Input templates are stored at a fixed N_side (512 for PySM 2 models), higher resolution models implemented in so_pysm_models have both 512 and 4096. Any other N_side is computed via ud_grade from the higher resolution map.
ud_grade is quick but heavily smoothes the spectrum, in fact both SO and S4 run simulations at 512 and 4096, so effectively never suffering this.

See for example impact of ud_grade on a power law signal with a slope of -2:

"Provide all the maps"

The simplest option is to pre-compute maps using alm2map after clipping the Alms to 3 * N_side - 1:

It is a lot of maps but the N_side 8192 is dominating the storage space anyway, the sum of all other maps is less than half of it.
We should provide maps at all N_sides from 8192 down to 128 (Lower resolution for SO)

In this case, we download, cache and read in memory the maps at the requested N_side, then do bandpass integration, finally do beam smoothing. So if we are simulating 1 channel at a time we do 1 map2alm and 1 alm2map. Just as a reference, at 2048 on Jupyter@NERSC map2alm on a IQU map takes 80s, alm2map takes 40s, so about 120s.

"Provide the Alms"

Alternatively we could provide the Alms instead of maps, that would be at variable resolution, for example from 8192 (quoting N_side, I assume ell_max is 3*N_side-1) down to 512.
Then we first apply the beam smoothing to all templates with almxfl, then transform to maps with alm2map at the target N_side. So the Sky object now would also have a specific beam instead of having only a specific N_side, but I guess most people run PySM for 1 component and 1 channel at a time.
Most components need to multiply maps in pixel space, so we need to transform the templates to maps before running the model.

Then we do bandpass integration on those templates, the output map is already smoothed.
For a model like GNILC dust, using N_side 2048 on Jupyter@NERSC, it would be 40s to apply alm2map to the template, 10s and 10s for spectral index and dust temperature.

File sizes

Size for each component in single precision, multiply by 3 for IQU, multiply by 2 for double precision.

Got in touch with Ginés Martínez-Solaeche and Jacques Delabrouille about implementing the model: "A 3-D model of polarised dust emission in the Milky Way" published in https://arxiv.org/abs/1706.04162

Would like first to complete and merge #37 that we could use as a reference.

It would be nice to know:

• is there already a Python implementation of the model?
• would a frequency-based interpolation of existing maps be a suitable implementation?
• does the model have any specific requirements? e.g. needs to load huge input datasets

Hi,

I need to run pysm3 in parallel, but I'm encountering a problem I can't solve.

I have this simplified version of my code:

 1 import healpy as hp
  2 import numpy as np
  3 import pysm3
  4 import pysm3.units as u
  5 import os
  6
  7
  8 from mpi4py import MPI
  9 comm = MPI.COMM_WORLD
 10 rank = comm.Get_rank()
 11 size = comm.Get_size()
 12
 13 write_dir =f'./test/{rank}/'
 14 if not os.path.exists(write_dir):
 15     os.makedirs(write_dir)
 16 test_map = np.arange(hp.nside2npix(128))
 17 hp.write_map(f'{write_dir}/test_map_{rank}.fits', test_map, overwrite=True, dtype=np.float32)
 18 sky = pysm3.Sky(nside=128, preset_strings=["s1"])
 19 hp.write_map(f'{write_dir}/test_map_{rank}_after_sky.fits', test_map, overwrite=True, dtype=np.float32)
 20 sky_extrap = sky.get_emission(145.*u.GHz)
 21 hp.write_map(f'{write_dir}/test_map_{rank}_after_getem.fits', test_map, overwrite=True, dtype=np.float32)
 22 hp.write_map(f'{write_dir}/sky_extrep_{rank}.fits', sky_extrap, overwrite=True, dtype=np.float32)

I'm trying to run it on a interactive job at NERSC: salloc -N 2 -C knl -q interactive -t 04:00:00
I then set export OMP_NUM_THREADS=2
and run the code with: mpirun -np 100 python test_parallel.py

what happens is the following:
the code correctly writes in 100 different folders the 100 test maps at line 17 and 19. But it writes the maps at line 21 an 22 only for a subset of processes (between 4 and 6 depending on the run).

Not that this happens for "s0","d0" ,"d1" (I haven't tried the others) but not for "c1"!

Any idea why this could happen?

Thanks a lot!

Hi ! I'm running PySM at NERSC on CORI and I have an error that appeared with version 3.2.0 and was not there with 3.1.0. This is when converting units. The simple script below triggers the error on Cori with 3.2.0 and not with 3.1.0:

`import numpy as np
import pysm
import pysm.units as u
from pysm import utils

print('PySM version:',pysm.version)
nus = np.linspace(130., 170., 12)
print('nus: ',nus)
filter_uK_CMB = np.ones(len(nus), dtype=np.double)
print('filter Init: ', filter_uK_CMB)
filter_uK_CMB_normalized = utils.normalize_weights(nus * u.GHz, filter_uK_CMB)
print('filter Normalized: ', filter_uK_CMB_normalized)
print ('All done')
`

With 3.1.0:
PySM version: 3.1.0 nus: [130. 133.63636364 137.27272727 140.90909091 144.54545455 148.18181818 151.81818182 155.45454545 159.09090909 162.72727273 166.36363636 170. ] filter Init: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] filter Normalized: [0.01866534 0.01972416 0.02081219 0.02192942 0.02307587 0.02425152 0.02545639 0.02669046 0.02795374 0.02924623 0.03056793 0.03191884] All done

while with 3.2.0:
PySM version: 3.2.0 nus: [130. 133.63636364 137.27272727 140.90909091 144.54545455 148.18181818 151.81818182 155.45454545 159.09090909 162.72727273 166.36363636 170. ] filter Init: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] Traceback (most recent call last): File "simpletest.py", line 11, in <module> filter_uK_CMB_normalized = utils.normalize_weights(nus * u.GHz, filter_uK_CMB) File "/global/homes/h/hamilton/.conda/envs/qubic_env/lib/python3.8/site-packages/pysm/utils/__init__.py", line 38, in normalize_weights (u.Jy / u.sr), equivalencies=u.cmb_equivalencies(freqs * u.GHz) File "/global/homes/h/hamilton/.conda/envs/qubic_env/lib/python3.8/site-packages/astropy/units/decorators.py", line 229, in wrapper _validate_arg_value(param.name, wrapped_function.__name__, File "/global/homes/h/hamilton/.conda/envs/qubic_env/lib/python3.8/site-packages/astropy/units/decorators.py", line 79, in _validate_arg_value raise UnitsError("Argument '{}' to function '{}' must be in units" astropy.units.core.UnitsError: Argument 'spec' to function 'cmb_equivalencies' must be in units convertible to 'GHz'.

Am I doing something wrong ?

Thanks in advance

The implementation of the GNILC-based dust templates is complete, see details and the executed notebook at:

#97 (comment)

Now it would be useful to have some independent validation, load the maps, take spectra, compare known regions, check everything is reasonable (see also #100)

See comparison between GNILC and PySM 2 here:
https://nbviewer.jupyter.org/gist/zonca/c6aa3e83f0151666e33a28cc84c0ba20

Suggestion by Hans Kristian (referring to the PySM 2 templates) was:

1. threshold beta_dust between 1.4 and 1.7, and smooth to ~2 degrees
   FWHM, to suppress the worst CO, point source and CIB effects
2. Smooth T_dust also to ~2 degrees, to reduce the worst Planck noise
   sensitivity

However no idea how to proceed. @brandonshensley @NicolettaK?

I think I must be doing something wrong with the synchrotron maps since I am unable to reproduce the Haslam 408MHz map; see attached notebook.

CheckHaslam_pysm3.ipynb.gz

As suggested by @zonca, I am starting an issue where we can discuss how the new models could be used as a baseline simulation suite, with the goal being to have models that span optimistic/baseline/pessimistic (or similar spirit). My initial suggestion is:

1. d9s4f1a1
2. d10s4f1a1
3. d12s3f1a2

where I have assumed s4 corresponds to the new synchrotron model with small scales in both amplitude and spectral index. Another point of discussion is whether some or all of these should include CO emission and polarization. Thoughts welcome.

• We need to provide guidelines on how to properly use the code, for example if people want maps at N_side 512, they should run the models at N_side 1024 and transform to target N_side. As we are smoothing at the end, we can easily provide a target_nside parameter, so people can smooth a 1024 map and get a 512 output without transforming twice.

Originally posted by @zonca in #90 (comment)

All PySM 2 models except HD2017 dust are implemented.

Looking for beta-testers, please look at the docs:
http://pysm3.readthedocs.org/

Try to install and run some models, if you find bugs, please open a dedicated issue.
If you have general feedback, please reply here. In particular I am looking for high level feedback:

• is there any feature you need that is missing?
• can the interface be simplified?
• do you have suggestions for better function/classes names?
• can the docs be improved?

Port local folder lookup from so_pysm_models to PySM:

https://github.com/simonsobs/so_pysm_models/blob/master/so_pysm_models/utils/__init__.py#L47

but also support an environment variable.

Add also docs on how to clone the data.

Related to hpc4cmb/toast#330

• mpi4py from conda-forge pulls mpich, it gives "cannot create executable" error
• testing now openmpi from conda-forge + mpi4py with pip

For the validation, I'd like to make some different realisations of foregrounds assuming Gaussian variance. Are there any way to implement this in PySM?
I can realise using healpy for each output frequency map. However, foreground maps should be correlated among frequencies. I think I have to create different realisations before the colour corrections.

PySM relies on libsharp to smooth maps over MPI.

Libsharp is not maintained anymore: https://gitlab.mpcdf.mpg.de/mtr/libsharp

Is anyone using PySM over MPI? Should we simplify the code and remove MPI support? or find an alternative for libsharp (@mreineck suggestions?) ?

@NicolettaK @giuspugl @bthorne93

It would be great if we can include this emission into pysm, similarly to what we've done in https://github.com/simonsobs/so_pysm_models/blob/master/so_pysm_models/co_lines.py

@zonca do you think that the interface changed that much ?

Longer term, it would be nice to import the hensley_draine_2017 dust model, which is possibly the most complex model available in PySM. This could even require changing the architecture of PySM 3.

As of astropy v3.1 , thermodynamic units are supported. Switch from our version defined in pysm.units.

MSS-001 production run

• 700 nodes
• 10 nodes/group
• 70 groups
• 8 MPI processes/node
• ~ 7 observation per group
• 2580 channels
• 1/20th year
• TOD channels per node 258
• PySM channels per group 2580/70 = 37, per node 4

Opening a new issue to discuss the synchrotron curvature parameter, which may affect the design of the Simulation Suite (#103). We have gotten some feedback that synchrotron spectral complexity is a high priority concern.

The s3 model has curvature as a single global parameter. I propose considering a new synchrotron model that has the amplitude and beta improvements of s5 but that also includes curvature. This could be implemented as a global mean value in agreement with s3, but with fluctuations following the same power spectrum as synchrotron beta in s5. That leaves only the amplitude of fluctuations to be decided, and unfortunately current data are not informative here. We could iterate starting with something modest like a 10% standard deviation around the mean at ~1 degree scales and then refine if that produces something either too extreme or completely unnoticeable.

Thoughts very much welcome @NicolettaK, @seclark, @zonca, @giuspugl, @bthorne93, others.

I've published to PyPI PySM 3.4.0 beta 1, you can install it with:

pip install --pre pysm3

(--pre grabs pre-releases)

Alternatively you can test in your browser using Google Colab (or use the notebook code locally):

https://colab.research.google.com/drive/1LI7zGC2Exd8R0HiYMl6b3Bnw-0bRp0Ws?usp=sharing

Please report any bugs here or open dedicated issue.
Feedback on the docs is particualrly appreciated, see https://pysm3.readthedocs.io/en/latest/models.html#dust, you can also click on "Edit this file on Github" to create a pull request directly.

Hi,
I'm working with bandpass integration. In order to check that results are correct I'm doing the following:
I integrate a CMB map in two different bands and then check that the integrated maps in thermodynamic units are still equal. This is what I expect for CMB in thermodynamic units also after band integration.

The code I'm using is the following:

 def check_unit_conversion(nuc, deltanu):
    bandpass_frequencies = np.linspace(nuc-deltanu/2., nuc+deltanu/2., 51) * u.GHz
    weights = np.ones(len(bandpass_frequencies))
    sky = pysm3.Sky(nside=128, preset_strings=["c1"])
    CMB_rj_int = sky.get_emission(bandpass_frequencies, weights)
    CMB_thermo_int = CMB_rj_int*pysm3.bandpass_unit_conversion(bandpass_frequencies, weights, u.uK_CMB)
    return CMB_thermo_int

CMB_thermo_int_100 = check_unit_conversion(100, 50)
CMB_thermo_int_400 = check_unit_conversion(400, 50)

if then I plot CMB_thermo_int_100[0]-CMB_thermo_int_400[0] this is what I get:

while what I expect is a map of zeros.

I've done some calculations and I think the problem is in bandpass_unit_conversion, I can share them with you, but first we should agree that this is a good test to do and that I'm actually implementing it in the right way.

I propose to continue the PySM 2 numbering, the first model should be d9.

Templates

The new GNILC dust IQU templates are being implemented in #82
What should we use for spectral index and dust temperature?

Resolution

The model will be natively at 8192, so we want to also have a version of it at 512, so it keeps being usable on laptops.
My proposal is to generate all templates both at 8192 and 512 from the same alms. Then if the user requests a NSIDE=512 or lower, the 512 map is loaded, otherwise the 8192 map.

Model implementation

Use ModifiedBlackBody as other standard dust components:

@giuspugl @seclark @bthorne93 @brandonshensley any feedback?

I am simulating pure CMB and because I want to update my cosmology easily, I use the pysm.CMBMap(nside, filename) with a set of I,Q,U maps I have simulated with synfast from my model (Note that I don't care about having the non-gaussian features from taylens here, so doing what I do is sufficient).

I do it like that:

cmbmap = pysm.CMBMap(nside, map_IQU=filename)
sky.add_component(cmbmap)

and then when I want to get this map from PySM, I do:

themaps_iqu = sky.get_emission([nus_edge[i], nus_edge[i + 1]] * u.GHz)

where nus_edges defin my bandpass, I then convert to uK_CMB with:

mysky[i, :, :] = np.array(themaps_iqu.to(u.uK_CMB, equivalencies=u.cmb_equivalencies(nus_in[i] * u.GHz))).T

where nus_in[i] is the center of my bandpass

Doing so gives pixels with a few percent difference from one band to another and they don;'t match the initial CMB maps. So I probably do something wrong here...

I could get consistent results by using single frequency instead of integration but I'm afraid this is not so good when I will be including something else but CMB...