Authors: S. Farrens
Language: Python 2.7
Download: GitHub
Description: A Python code designed for PSF deconvolution using a low-rank approximation and sparsity. The code can handle a fixed PSF for the entire field or a stack of PSFs for each galaxy position.
Notes: This code was used to produce the results presented in the paper Space variant deconvolution of galaxy survey images. Sample Euclid-like PSF data can be downloaded from here [63Mb]


The following sections provide some details for how to run sf_deconvolve.

The directory lib contains all of the primary functions and classes used for optimisation and analysis. functions contains some additional generic functions and tools.


In order to run the code in this repository the following packages must be installed:

  • Python 2.7 [Tested with v 2.7.11]
  • Numpy [Tested with v 1.11.3]
  • Scipy [Tested with v 0.18.1]
  • Astropy [Tested with v 1.1.2]
  • Matplotlib [Tested with v 1.5.3]
  • Termcolor [Tested with v 1.1.0]

The current implementation of wavelet transformations additionally requires the C++ script from the Sparse2D library in the iSap package [Tested with v 3.1]. These C++ scripts will be need to be compiled in order to run (see iSap Documentation for details).

The low-rank approximation method can be run purely in Python.


The primary code is which is designed to take an observed (i.e. with PSF effects and noise) stack of galaxy images and a known PSF, and attempt to reconstruct the original images. The input format are Numpy binary files (.npy) or FITS image files (.fits).

The code can be run as follows:


 Where INPUT_IMAGES.npy denotes the Numpy binary file containing the stack of observed galaxy images, PSF.npy denotes the PSF corresponding to each galaxy image and OUTPUT_NAME specifies the output path and file name.

Alternatively the code arguments can be stored in a configuration file (with any name) and the code can be run by providing the file name preceded by a @.

$ @config.ini


The following example can be run on the sample data provided in the example directory.

This example takes a sample of 100 galaxy images (with PSF effects and added noise) and the corresponding PSFs, and recovers the original images using low-rank approximation via Condat-Vu optimisation.

$ -i example_image_stack.npy -p example_psf.npy -o example_output --mode lowr

The example can also be run using the configuration file provided.

The result will be two Numpy binary files called example_output_primal.npy and example_output_dual.npy corresponding to the primal and dual variables in the splitting algorithm. The reconstructed images will be in the example_output_primal.npy file.

Code Options

Required Arguments

-i INPUT, --input INPUT: Input data file name. File should be a Numpy binary containing a stack of noisy galaxy images with PSF effects (i.e. a 3D array).

-p PSF, --psf PSF: PSF file name. File should be a Numpy binary containing either: (a) a single PSF (i.e. a 2D array for fixed format) or (b) a stack of PSFs corresponding to each of the galaxy images (i.e. a 3D array for obj_var format).

Optional Arguments

-h, --help: Show the help message and exit.

-v, --version: Show the program's version number and exit.

-q, --quiet: Suppress verbose for each iteration.

-o, --output: Output file name. If not specified output files will placed in input file path.

--output_format Output file format [npy or fits].


-k, --current_res: Current deconvolution results file name (i.e. the file containing the primal results from a previous run).

--noise_est: Initial estimate of the noise standard deviation in the observed galaxy images. If not specified this quantity is automatically calculated using the median absolute deviation of the input image(s).


-m, --mode {all,sparse,lowr,grad}: Option to specify the optimisation mode [all, sparse, lowr or grad]. all performs optimisation using both low-rank approximation and sparsity, sparse using only sparsity, lowr uses only low-rank and grad uses only gradient descent. (default: lowr)

--opt_type {condat,fwbw,gfwbw}: Option to specify the optimisation method to be implemented [condat, fwbw or gfwbw]. condat implements the Condat-Vu proximal splitting method, fwbw implements Forward-Backward splitting with FISTA speed-up and gfwbw implements the generalised Forward-Backward splitting method. (default: condat)

--n_iter: Number of iterations. (default: 150)

--cost_window: Window to measure cost function (i.e. interval of iterations for which cost should be calculated). (default: 1)

--convergence: Convergence tolerance. (default: 0.0001)

--no_pos: Option to turn off positivity constraint.

--no_grad: Option to turn off gradient calculation.

Low-Rank Aproximation

--lowr_thresh_factor: Low rank threshold factor. (default: 1)

--lowr_type: Type of low-rank regularisation [standard or ngole]. (default: standard)

--lowr_thresh_type: Low rank threshold type [soft or hard]. (default: hard)


--wavelet_type: Type of Wavelet to be used (see iSap Documentation). (default: 1)

--wave_thresh_factor: Wavelet threshold factor. (default: [3.0, 3.0, 4.0])

--n_reweights: Number of reweightings. (default: 1)

Condat Algorithm

--relax: Relaxation parameter (rho_n in Condat-Vu method). (default: 0.8)

--condat_sigma: Condat proximal dual parameter. (default: 0.5)

--condat_tau: Condat proximal primal parameter. (default: 0.5)


-c, --clean_data: Clean data file name.

-r, --random_seed: Random seed. Use this option if the input data is a randomly selected subset (with known seed) of the full sample of clean data.

--kernel: Standard deviation of pixels for Gaussian kernel. This option will multiply the deconvolution results by a Gaussian kernel.

--metric: Metric to average errors [median or mean]. (default: median)


If you get the following error:

ERROR: svd() got an unexpected keyword argument 'lapack_driver'

Update your Numpy and Scipy installations

$ pip install --upgrade numpy
$ pip install --upgrade scipy



Paper accepted : New inpainting method to handle colored-noise data to test the weak equivalence principle

The context

The MICROSCOPE space mission, launched on April 25, 2016, aims to test the weak equivalence principle (WEP) with a 1015 precision. To reach this performance requires an accurate and robust data analysis method, especially since the possible WEP violation signal will be dominated by a strongly colored noise. An important complication is brought by the fact that some values will be missing –therefore, the measured time series will not be strictly regularly sampled. Those missing values induce a spectral leakage that significantly increases the noise in Fourier space, where the WEP violation signal is looked for, thereby complicating scientific returns.

FIG. 1 (from Pires et al, 2016): The black curve shows the MICROSCOPE PSD es- timate for a 120 orbits simulation. An example of a possible EPV signal of 3 × 10−15 in the inertial mode is shown by the peak at 1.8 × 10−4 Hz. The grey curve shows the spectral leakage affecting the PSD estimate when gaps are present in the data.















The results

Recently, we developed an inpainting algorithm to correct the MICROSCOPE data for missing values (in red, fig 4). This code has been integrated in the official MICROSCOPE data processing and analysis pipeline because it enables us to significantly measure an equivalence principle violation (EPV) signal in a model-independent way, in the inertial satellite configuration. In this work, we present several improvements to the method that may allow us now to reach the MICROSCOPE requirements for both inertial and spin satellite configurations (green curve, fig. 4).



FIG. 4 (from Pires et al, 2016): MICROSCOPE differential acceleration PSD estimates averaged over 100 simulations in the inertial mode (upper panel) and in the spin mode (lower panel). The black lines show the PSD estimated when all the data is available, the red lines show the PSD estimated from data filled with the inpainting method developed in Paper I and the green lines show the PSD estimated from data filled with the new inpainting method (ICON) presented in this paper.

The code ICON

The code corresponding to the paper is available for download here.

Although, the inpainting method presented in this paper has been optimized to the MICROSCOPE data, it remains sufficiently general to be used in the general context of missing data in time series dominated by an unknown colored-noise.


 Dealing with missing data in the MICROSCOPE space mission: An adaptation of inpainting to handle colored-noise data, S. Pires, J. Bergé, Q. Baghi, P. Touboul, G. Métris, accepted in Physical Review D, December 2016

Dealing with missing data: An inpainting application to the MICROSCOPE space mission, J. Bergé, S. Pires, Q. Baghi, P. Touboul, G. Metris, Physical Review D, 92, 11, December 2015 


CFIS proposal accepted

On the day of the Brexit outcome, so disastrous for Europe and the UK, there is at least good news for the cosmological community: CFIS, the Canada-France Imaging Survey, has been accepted! This survey consists of two parts. The WIQD​(Wide ­ Image Quality ­ Deep) part will cover 5,000 deg2 of the Northern sky, observed in the r-band with the CFHT (Canada-France Hawai'i telescope). The u-band will cover 10,000 deg2 to a lower depth, and is part of LUAU (Legacy for the U­-band All­-sky Universe). 271 nights are granted, observations will start in 8 months from now.CFIS Logo

CFIS will allow us to study properties of dark-matter structures, including filaments between galaxy clusters and groups, stripping of dark-matter halos of satellite galaxies in clusters, and the shapes of dark-matter halos. In addition, the laws of gravity on large scales will be tested, and modifications to Einstein's theory of general relativity will be looked for. CFIS will observe a very large number of distant, high-redshift galaxies, and will use techniques of galaxy clustering and weak gravitational lensing to achieve its goals.

In addition, CFIS will create synergie with other ongoing and planned surveys: CFIS will provide ground-based optical data for Euclid photometric-redshifts. It will produce a very useful imaging data set for target selection for spectroscopic surveys such as DESI, WEAVE, and MSE. It will further provide optical data of galaxy clusters that will enhance the science outcome of the X-ray mission eRosita.

PIs: Jean-Charles Cuillandre (CEA Saclay/France) & Alan McConnachie (Victoria/Canada).
CosmoStat participants: Martin Kilbinger, Jean-Luc Starck. Sandrine Pires.
Irfu participants: Monique Arnaud, Hervé Aussel, Olivier Boulade, Pierre-Alain Duc, David Elbaz, Christophe Magneville, Yannick Mellier, Marguerite Pierre, Anand Raichoor, Jim Rich, Vanina Ruhlmann-Kleider, Marc Sauvage, Christophe Yèche.


École d'été 2016 en cosmologie

C'est l'annonce d'une école d'été en cosmologie qui va avoir lieu cette août, avec participation des membres de CosmoStat.

Dates: 22 au 27 août 2016
Lieu: Narbonne, France
Site web:

Date limite d'inscription: 3 juillet 2016
Comité d'organisation: Alain Blanchard, pour la communité Euclid-France
L'objectif prioritaire est d'offrir à la communauté une initiation et un approfondissement
aux problématiques des grands relevés cosmologiques appropriés aux scientifiques qui
s'impliquent fortement dans le projet Euclid comme dans les autres grands projets dédiés
à l'étude de l'énergie noire.
Programme et intervenants:
* Modèle cosmologique standard, énergie noire (Martin Kunz)
* Clustering des galaxies (Francis Bernardeau, Olivier LeFèvre, Stéphanie Escoffier)
* Lentillage gravitationnel (Martin Kilbinger)
* Amas de galaxies (Alain Blanchard)
* Paramètres de nuisance (Sandrine Codis)

* Outils: Python (Yannick Copin)
* Outils: CLASS (Julien Lesgourgues)

Les cours sont complémentés par des TD.

Slides for weak-lensing (no figure to reduce file size):
Day 1