get_map() method of orphics.maps.MapGen class is broken

[EEmGuzman]

Hi Mat,

The get_map() method in the orphics.maps.MapGen class is broken when using Pixell Version 0.10.3.

Running the command:
unlensed,kappa,lensed,beamed,noise_map,observed = flsims.get_sim(return_intermediate=True)

Results in an error of:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-1017c01d7b16> in <module>
     17 # generate maps
     18 print("Pixell Version {}".format(pixell.__version__))
---> 19 unlensed,kappa,lensed,beamed,noise_map,observed = flsims.get_sim(return_intermediate=True)

~/.local/lib/python3.8/site-packages/orphics/lensing.py in get_sim(self, seed_cmb, seed_kappa, seed_noise, lens_order, return_intermediate, skip_lensing, cfrac)
    269         else:
    270             if not(self._fixed):
--> 271                 kappa = self.get_kappa(seed_kappa)
    272                 self.kappa = kappa
    273                 self.alpha = alpha_from_kappa(kappa,posmap=self.posmap)

~/.local/lib/python3.8/site-packages/orphics/lensing.py in get_kappa(self, seed)
    261         return self.mgen.get_map(seed=seed)
    262     def get_kappa(self,seed=None):
--> 263         return self.kgen.get_map(seed=seed)
    264     def get_sim(self,seed_cmb=None,seed_kappa=None,seed_noise=None,lens_order=5,return_intermediate=False,skip_lensing=False,cfrac=None):
    265         unlensed = self.get_unlensed(seed_cmb)

~/.local/lib/python3.8/site-packages/orphics/maps.py in get_map(self, seed, scalar, iau, real, harm)
    452                 print("covsqrt shape {}".format(np.shape(self.covsqrt)))
    453                 print("rand shape {}".format(np.shape(rand)))
--> 454                 data = enmap.map_mul(self.covsqrt, rand)
    455                 print("map mul shape {}".format(np.shape(data)))
    456                 kmap = enmap.ndmap(data, self.wcs)

~/.local/lib/python3.8/site-packages/pixell/enmap.py in map_mul(mat, vec)
   1099         if mat.ndim == 2 and vec.ndim == 2: return mat*vec
   1100         # Otherwise we do a matrix product along the last axes
-> 1101         ovec = samewcs(np.einsum("...abyx,...byx->...ayx", mat, vec), mat, vec)
   1102         return ovec
   1103 

<__array_function__ internals> in einsum(*args, **kwargs)

~/bin/anaconda3/envs/standard_dpipeline/lib/python3.8/site-packages/numpy/core/einsumfunc.py in einsum(out, optimize, *operands, **kwargs)
   1348         if specified_out:
   1349             kwargs['out'] = out
-> 1350         return c_einsum(*operands, **kwargs)
   1351 
   1352     # Check the kwargs to avoid a more cryptic error later, without having to

ValueError: einstein sum subscripts string contains too many subscripts for operand 1
---------------------------------------------------------------------------

I printed the dimensions of what was being multiplied when it failed:

covsqrt shape (1, 1, 128, 128)
rand shape (128, 128)

Expanding the dimension of the rand variable in the get_map() method when rand.ndim==2 then reducing the dimension again after the multiplication fixed the issue on my end, so I assume the problem is isolated to the case of rand.ndim==2.

[debuadakiitk]

Facing same issue. Can I request you to part of the code that you modified to resolve this issue? or it is related to numpy version. Tried with tensorflow routine. Could not resolve

[EEmGuzman]

@debuadakiitk I am still using orphics but have not kept up with the latest versions, so hopefully my fix still works.

To fix the issue I changed the get_map() function of the MapGen class in maps.py to:

def get_map(self,seed=None,scalar=False,iau=False,real=False,harm=False):
                if seed is not None: np.random.seed(seed)
                rand = enmap.fft(enmap.rand_gauss(self.shape, self.wcs)) if real else enmap.rand_gauss_harm(self.shape, self.wcs)
                # Changes begin here
                if rand.ndim == 2:
                    rand = np.expand_dims(rand, axis=0)
                    data = (enmap.map_mul(self.covsqrt, rand))[0]
                else:
                    data = enmap.map_mul(self.covsqrt, rand)
                # Changes end here
                kmap = enmap.ndmap(data, self.wcs)
                if harm: 
                    return kmap
                else:
                    if scalar:
                            return enmap.ifft(kmap).real
                    else:
                            return enmap.harm2map(kmap,iau=iau)

If I remember correctly, just modifying the axes of the array fixes the issue with the numpy.einsum() in pixell. I can't quite remember if I made any other changes anywhere else in the code, but I think that was all I did. Ill check later today if just making these few changes fixes the current version of orphics as well.

[debuadakiitk]

Thanks a lot. If you find some other modifications, kindly let me know. It would be helpful for me. thanks, Debabrata

…

On Tue, Nov 8, 2022 at 10:55 PM EEmGuzman ***@***.***> wrote: @debuadakiitk <https://github.com/debuadakiitk> I am still using orphics but have not kept up with the latest versions, so hopefully my fix still works. To fix the issue I changed the get_map() function of the MapGen class in maps.py to: def get_map(self,seed=None,scalar=False,iau=False,real=False,harm=False): if seed is not None: np.random.seed(seed) rand = enmap.fft(enmap.rand_gauss(self.shape, self.wcs)) if real else enmap.rand_gauss_harm(self.shape, self.wcs) # Changes begin here if rand.ndim == 2: rand = np.expand_dims(rand, axis=0) data = (enmap.map_mul(self.covsqrt, rand))[0] else: data = enmap.map_mul(self.covsqrt, rand) # Changes end here kmap = enmap.ndmap(data, self.wcs) if harm: return kmap else: if scalar: return enmap.ifft(kmap).real else: return enmap.harm2map(kmap,iau=iau) If I remember correctly, just modifying the axes of the array fixes the issue with the numpy.einsum() in pixell. I can't quite remember if I made any other changes anywhere else in the code, but I think that was all I did. Ill check later today if just making these few changes fixes the current version of orphics as well. — Reply to this email directly, view it on GitHub <#25 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A4CNW34VRVAMXB6XKLLGGSTWHKEHBANCNFSM4W6LDRFA> . You are receiving this because you were mentioned.Message ID: ***@***.***>