Creating a global map using foscat from a single model#
Context#
This notebook demonstrates the use of FOSCAT packages to synthesize new maps from an initial target map. The synthesized maps retain the same Scattering Covariance (SC) statistics as the target map. The target map is here a Large Scale Structure simulation of weak lensing, from the CosmoGrid data set (Kacprzak et al. 2023; Fluri et al. 2022); utilized in this example is identical to the one in Mousset et al. (submitted 2024), and the synthesis method is nearly the same.
The method#
The idea is to calculate the SC statistic \(\Phi(d)\) for the target map \(d\) and generate a new map \(u\) that minimizes the following loss: \(L = \sum_k{\frac{(\Phi(d)_k - \Phi(u)_k)^{2}}{\sigma_k^{2}}}\)
Here, \(k\) denotes the \(k^{th}\) SC coefficients, and \(\sigma_k^{2}\) represents the variance of \(\Phi(d)_k\) across the target map. Specifically, \(\Phi(d)_k\) values are obtained by integrating over the entire domain wavelet maps, while \(\sigma_k^{2}\) values correspond to the variance within this domain.
The map \(u\) is initially set to random Gaussian noise, and the final result is obtained by iteratively minimizing \(L\) using gradient descent. The resulting map \(u\) is then returned.
For further details on the synthesis domain, please refer to:#
Fluri, J., Kacprzak, T., Lucchi, A., et al. 2022, Phys. Rev. D, 105, 083518
Kacprzak, T., Fluri, J., Schneider, A., Refregier, A., & Stadel, J. 2023, Journal of Cosmology and Astroparticle Physics, 2023, 05
Mousset, L., Allys E., Pric, M. A., . Aumo, J., Delouis, J.-M., Montier, L., McEwen, J. D., Sumitted 2024 Astronomy & Astrophysicwen2 07
Installation of required packages#
!pip install foscat==3.1.3
Collecting foscat==3.1.3
Using cached foscat-3.1.3-py3-none-any.whl.metadata (1.0 kB)
Requirement already satisfied: imageio in /srv/conda/envs/notebook/lib/python3.12/site-packages (from foscat==3.1.3) (2.34.2)
Requirement already satisfied: imagecodecs in /srv/conda/envs/notebook/lib/python3.12/site-packages (from foscat==3.1.3) (2024.6.1)
Requirement already satisfied: matplotlib in /srv/conda/envs/notebook/lib/python3.12/site-packages (from foscat==3.1.3) (3.9.2)
Requirement already satisfied: numpy in /srv/conda/envs/notebook/lib/python3.12/site-packages (from foscat==3.1.3) (1.26.4)
Requirement already satisfied: tensorflow in /srv/conda/envs/notebook/lib/python3.12/site-packages (from foscat==3.1.3) (2.16.2)
Collecting healpy (from foscat==3.1.3)
Using cached healpy-1.17.3-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (8.9 kB)
Requirement already satisfied: astropy in /srv/conda/envs/notebook/lib/python3.12/site-packages (from healpy->foscat==3.1.3) (6.1.2)
Requirement already satisfied: pillow>=8.3.2 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from imageio->foscat==3.1.3) (10.4.0)
Requirement already satisfied: contourpy>=1.0.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (4.53.1)
Requirement already satisfied: kiwisolver>=1.3.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (24.1)
Requirement already satisfied: pyparsing>=2.3.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from matplotlib->foscat==3.1.3) (2.8.2)
Requirement already satisfied: absl-py>=1.0.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (2.1.0)
Requirement already satisfied: astunparse>=1.6.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (1.6.3)
Requirement already satisfied: flatbuffers>=23.5.26 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (24.3.25)
Requirement already satisfied: gast!=0.5.0,!=0.5.1,!=0.5.2,>=0.2.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (0.5.5)
Requirement already satisfied: google-pasta>=0.1.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (0.2.0)
Requirement already satisfied: h5py>=3.10.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (3.11.0)
Requirement already satisfied: ml-dtypes~=0.3.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (0.3.2)
Requirement already satisfied: opt-einsum>=2.3.2 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (3.3.0)
Requirement already satisfied: protobuf!=4.21.0,!=4.21.1,!=4.21.2,!=4.21.3,!=4.21.4,!=4.21.5,<5.0.0dev,>=3.20.3 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (4.25.3)
Requirement already satisfied: requests<3,>=2.21.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (2.32.3)
Requirement already satisfied: setuptools in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (72.1.0)
Requirement already satisfied: six>=1.12.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (1.16.0)
Requirement already satisfied: termcolor>=1.1.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (2.4.0)
Requirement already satisfied: typing-extensions>=3.6.6 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (4.12.2)
Requirement already satisfied: wrapt>=1.11.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (1.16.0)
Requirement already satisfied: grpcio<2.0,>=1.24.3 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (1.62.2)
Requirement already satisfied: tensorboard<2.17,>=2.16 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (2.16.2)
Requirement already satisfied: keras>=3.0.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorflow->foscat==3.1.3) (3.5.0)
Requirement already satisfied: wheel<1.0,>=0.23.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from astunparse>=1.6.0->tensorflow->foscat==3.1.3) (0.44.0)
Requirement already satisfied: rich in /srv/conda/envs/notebook/lib/python3.12/site-packages (from keras>=3.0.0->tensorflow->foscat==3.1.3) (13.7.1)
Requirement already satisfied: namex in /srv/conda/envs/notebook/lib/python3.12/site-packages (from keras>=3.0.0->tensorflow->foscat==3.1.3) (0.0.8)
Requirement already satisfied: optree in /srv/conda/envs/notebook/lib/python3.12/site-packages (from keras>=3.0.0->tensorflow->foscat==3.1.3) (0.12.1)
Requirement already satisfied: charset-normalizer<4,>=2 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from requests<3,>=2.21.0->tensorflow->foscat==3.1.3) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from requests<3,>=2.21.0->tensorflow->foscat==3.1.3) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from requests<3,>=2.21.0->tensorflow->foscat==3.1.3) (1.26.19)
Requirement already satisfied: certifi>=2017.4.17 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from requests<3,>=2.21.0->tensorflow->foscat==3.1.3) (2024.7.4)
Requirement already satisfied: markdown>=2.6.8 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorboard<2.17,>=2.16->tensorflow->foscat==3.1.3) (3.6)
Requirement already satisfied: tensorboard-data-server<0.8.0,>=0.7.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorboard<2.17,>=2.16->tensorflow->foscat==3.1.3) (0.7.0)
Requirement already satisfied: werkzeug>=1.0.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from tensorboard<2.17,>=2.16->tensorflow->foscat==3.1.3) (3.0.3)
Requirement already satisfied: pyerfa>=2.0.1.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from astropy->healpy->foscat==3.1.3) (2.0.1.4)
Requirement already satisfied: astropy-iers-data>=0.2024.7.1.0.34.3 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from astropy->healpy->foscat==3.1.3) (0.2024.8.5.0.32.23)
Requirement already satisfied: PyYAML>=3.13 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from astropy->healpy->foscat==3.1.3) (6.0.2)
Requirement already satisfied: MarkupSafe>=2.1.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from werkzeug>=1.0.1->tensorboard<2.17,>=2.16->tensorflow->foscat==3.1.3) (2.1.5)
Requirement already satisfied: markdown-it-py>=2.2.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from rich->keras>=3.0.0->tensorflow->foscat==3.1.3) (3.0.0)
Requirement already satisfied: pygments<3.0.0,>=2.13.0 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from rich->keras>=3.0.0->tensorflow->foscat==3.1.3) (2.18.0)
Requirement already satisfied: mdurl~=0.1 in /srv/conda/envs/notebook/lib/python3.12/site-packages (from markdown-it-py>=2.2.0->rich->keras>=3.0.0->tensorflow->foscat==3.1.3) (0.1.2)
Using cached foscat-3.1.3-py3-none-any.whl (73 kB)
Using cached healpy-1.17.3-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.5 MB)
Installing collected packages: healpy, foscat
Successfully installed foscat-3.1.3 healpy-1.17.3
import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
import foscat.scat_cov as sc
import foscat.Synthesis as synthe
Convert Data in HEALPix#
The above gap-filling method does not take into account that our subject is a sphere. Because we work with full globe data, we need to treat the grid system that conveys the shape of the sphere. To do that, we use one of the Discrete Global Grid Systems (DGGS), HEALPix.
HEALPix stands for Hierarchical Equal Area isoLatitude Pixelation of a sphere. This pixelation produces a subdivision of a spherical surface in which each pixel covers the same surface area as every other pixel.
See https://healpix.sourceforge.io and/or the HEALPix Primer for more information.
The healpy tutorial is also a very good starting point to understand more about HEALPix.
Resolution#
The resolution of the grid is expressed by the parameter Nside, which defines the number of divisions along the side of a base-resolution pixel that is needed to reach a desired high-resolution partition.
Ordering Systems#
HEALPix supports two pixel ordering systems: nested and ring.
Detailed explanations of the two pixel ordering systems can be found at https://healpix.jpl.nasa.gov/html/intronode4.htm.
In our exam
choose the resolution#
we use \(nside=64\) in this context to prevent excessive computation.ple we use nested
nside=64
Get Large Scale Structure data#
The target map is here a LSS simulation of weak lensing, from the CosmoGrid data set (Kacprzak et al. 2023; Fluri et al. 2022);
!wget -O target_map_lss.npy https://github.com/astro-informatics/s2scat/raw/main/notebooks/data/target_map_lss.npy
--2024-09-09 18:40:02-- https://github.com/astro-informatics/s2scat/raw/main/notebooks/data/target_map_lss.npy
Resolving github.com (github.com)... 140.82.121.4
Connecting to github.com (github.com)|140.82.121.4|:443... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://raw.githubusercontent.com/astro-informatics/s2scat/main/notebooks/data/target_map_lss.npy [following]
--2024-09-09 18:40:02-- https://raw.githubusercontent.com/astro-informatics/s2scat/main/notebooks/data/target_map_lss.npy
Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 2606:50c0:8002::154, 2606:50c0:8000::154, 2606:50c0:8001::154, ...
Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|2606:50c0:8002::154|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1046656 (1022K) [application/octet-stream]
Saving to: ‘target_map_lss.npy’
target_map_lss.npy 100%[===================>] 1022K --.-KB/s in 0.04s
2024-09-09 18:40:03 (23.7 MB/s) - ‘target_map_lss.npy’ saved [1046656/1046656]
from scipy.interpolate import RegularGridInterpolator
# convert the input data in a nside=128 healpix map
l_nside=128
im = np.load('target_map_lss.npy')
xsize,ysize=im.shape
# Define the new row and column to be added to prepare the interpolation
new_row = im[0:1,:] # A new row with N elements (the other longitude)
new_column = np.concatenate([im[:,0:1],im[-2:-1,0:1]],0) # A new column with N+1 elements to add previous latitude
# Add the new row to the array
im = np.vstack([im, new_row])
# Add the new column to the array with the new row
im = np.hstack([im, new_column])
# Create a grid of coordinates corresponding to the array indices
x = np.linspace(0, im.shape[0]-1, im.shape[0])
y = np.linspace(0, im.shape[1]-1, im.shape[1])
# Create an interpolator
interpolator = RegularGridInterpolator((x, y), im)
# List of healpix coordinate to interpol
colatitude,longitude = hp.pix2ang(l_nside,np.arange(12*l_nside**2),nest=True)
coords = np.concatenate([colatitude/np.pi*xsize,longitude/(2*np.pi)*ysize],0).reshape(2,colatitude.shape[0]).T
# Perform the interpolation
heal_im = interpolator(coords)
#reduce the final map to the expected resolution
if nside>128:
th,ph=hp.pix2ang(nside,np.arange(12*nside**2),nest=True)
heal_im=hp.get_interp_val(heal_im,th,ph,nest=True)
else:
heal_im=np.mean(heal_im.reshape(12*nside**2,(l_nside//nside)**2),1)
hp.mollview(heal_im,cmap='plasma',nest=True,title='LSS',min=-3,max=3)
# free memory
del coords
del interpolator
del colatitude
del longitude
Initialise FOSCAT scattering covariance description#
scat_op=sc.funct(NORIENT=4, # define the number of wavelet orientation
KERNELSZ=3, #KERNELSZ, # define the kernel size
JmaxDelta=0, # The used Jmax is Jmax-JmaxDelta
all_type='float64')
================================================
START FOSCAT CONFIGURATION
================================================
2024-09-09 18:40:15.486772: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:10575] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered
2024-09-09 18:40:15.486800: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:479] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered
2024-09-09 18:40:15.488952: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1442] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered
2024-09-09 18:40:15.494428: I tensorflow/core/platform/cpu_feature_guard.cc:210] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: SSE4.1 SSE4.2 AVX AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
Num GPUs Available: 1
1 Physical GPUs, 1 Logical GPUs
============================================================
== ==
== ==
== RUN ON GPU Rank 0 : /device:GPU:0 ==
== ==
== ==
============================================================
2024-09-09 18:40:17.723128: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1928] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 43598 MB memory: -> device: 0, name: NVIDIA A40, pci bus id: 0000:a3:00.0, compute capability: 8.6
Init Wave 2
Init Wave 4
Init Wave 8
Init Wave 16
Init Wave 32
Compute the reference statistics#
ref,sref=scat_op.eval(heal_im,calc_var=True)
ref.plot(name='LSS',color='b')
sref.plot(name='Variance',color='r',hold=False)
Define the loss function#
def The_loss(u,scat_operator,args):
ref = args[0]
sref = args[1]
# compute scattering covariance of the current synthetised map called u
learn=scat_operator.eval(u)
# make the difference withe the reference coordinates
loss=scat_operator.reduce_distance(learn,ref,sigma=sref)
return(loss)
Build the synthesis model#
loss1=synthe.Loss(The_loss,scat_op,ref,sref)
sy = synthe.Synthesis([loss1])
Run the synthesis#
imap=np.random.randn(12*nside**2)*np.std(heal_im)
inscat=scat_op.eval(imap)
omap=sy.run(imap,
EVAL_FREQUENCY=100,
NUM_EPOCHS = 1000).numpy()
Total number of loss 1
Run [PROC=0000] on GPU /device:GPU:0
WARNING:tensorflow:AutoGraph could not transform <bound method Socket.send of <zmq.Socket(zmq.PUSH) at 0x7f57603fa890>> and will run it as-is.
Please report this to the TensorFlow team. When filing the bug, set the verbosity to 10 (on Linux, `export AUTOGRAPH_VERBOSITY=10`) and attach the full output.
Cause: module, class, method, function, traceback, frame, or code object was expected, got cython_function_or_method
To silence this warning, decorate the function with @tf.autograph.experimental.do_not_convert
WARNING: AutoGraph could not transform <bound method Socket.send of <zmq.Socket(zmq.PUSH) at 0x7f57603fa890>> and will run it as-is.
Please report this to the TensorFlow team. When filing the bug, set the verbosity to 10 (on Linux, `export AUTOGRAPH_VERBOSITY=10`) and attach the full output.
Cause: module, class, method, function, traceback, frame, or code object was expected, got cython_function_or_method
To silence this warning, decorate the function with @tf.autograph.experimental.do_not_convert
Itt 0 L= 2.15e+03 ( 2.15e+03 ) 24.433s
Itt 100 L= 6.18 ( 6.18 ) 2.853s
Itt 200 L= 2.36 ( 2.36 ) 2.857s
Itt 300 L= 1.5 ( 1.5 ) 2.967s
Itt 400 L= 1.16 ( 1.16 ) 2.771s
Itt 500 L= 0.982 ( 0.982 ) 2.875s
Itt 600 L= 0.855 ( 0.855 ) 2.831s
Itt 700 L= 0.763 ( 0.763 ) 2.791s
Itt 800 L= 0.701 ( 0.701 ) 2.795s
Itt 900 L= 0.649 ( 0.649 ) 2.790s
Itt 1000 L= 0.616 ( 0.616 ) 2.815s
Plot the results#
plt.figure(figsize=(12,10))
hp.mollview(heal_im,cmap='plasma',nest=True,title='Target LSS',min=-3,max=3,
hold=False,sub=(1,2,1),cbar=False)
hp.mollview(omap,cmap='plasma',nest=True,title='Synthesised LSS',min=-3,max=3,
hold=False,sub=(1,2,2),cbar=False)
synth=scat_op.eval(omap)
inscat=scat_op.eval(imap)
inscat.plot(name='Input White noise',color='grey',norm=True)
ref.plot(name='Target LSS',color='b',norm=True,lw=4,hold=False)
synth.plot(name='Synthesised LSS',color='r',hold=False,norm=True)
