Skip to content
- Low-rank plus sparse decomposition for exoplanet detection in direct-imaging ADI sequences. The LLSG algorithm
- Context. Data processing constitutes a critical component of high-contrast exoplanet imaging. Its role is almost as important as the choice of a coronagraph or a wavefront control system, and it is intertwined with the chosen observing strategy. Among the data processing techniques for angular differential imaging (ADI), the most recent is the family of principal component analysis (PCA) based algorithms. It is a widely used statistical tool developed during the first half of the past century. PCA serves, in this case, as a subspace projection technique for constructing a reference point spread function (PSF) that can be subtracted from the science data for boosting the detectability of potential companions present in the data. Unfortunately, when building this reference PSF from the science data itself, PCA comes with certain limitations such as the sensitivity of the lower dimensional orthogonal subspace to non-Gaussian noise. Aims: Inspired by recent advances in machine learning algorithms such as robust PCA, we aim to propose a localized subspace projection technique that surpasses current PCA-based post-processing algorithms in terms of the detectability of companions at near real-time speed, a quality that will be useful for future direct imaging surveys. Methods: We used randomized low-rank approximation methods recently proposed in the machine learning literature, coupled with entry-wise thresholding to decompose an ADI image sequence locally into low-rank, sparse, and Gaussian noise components (LLSG). This local three-term decomposition separates the starlight and the associated speckle noise from the planetary signal, which mostly remains in the sparse term. We tested the performance of our new algorithm on a long ADI sequence obtained on β Pictoris with VLT/NACO. Results: Compared to a standard PCA approach, LLSG decomposition reaches a higher signal-to-noise ratio and has an overall better performance in the receiver operating characteristic space. This three-term decomposition brings a detectability boost compared to the full-frame standard PCA approach, especially in the small inner working angle region where complex speckle noise prevents PCA from discerning true companions from noise.
- FARGO3D: A New GPU-oriented MHD Code
- We present the FARGO3D code, recently publicly released. It is a magnetohydrodynamics code developed with special emphasis on the physics of protoplanetary disks and planet-disk interactions, and parallelized with MPI. The hydrodynamics algorithms are based on finite-difference upwind, dimensionally split methods. The magnetohydrodynamics algorithms consist of the constrained transport method to preserve the divergence-free property of the magnetic field to machine accuracy, coupled to a method of characteristics for the evaluation of electromotive forces and Lorentz forces. Orbital advection is implemented, and an N-body solver is included to simulate planets or stars interacting with the gas. We present our implementation in detail and present a number of widely known tests for comparison purposes. One strength of FARGO3D is that it can run on either graphical processing units (GPUs) or central processing units (CPUs), achieving large speed-up with respect to CPU cores. We describe our implementation choices, which allow a user with no prior knowledge of GPU programming to develop new routines for CPUs, and have them translated automatically for GPUs.
- AstroBlend: An astrophysical visualization package for Blender
- The rapid growth in scale and complexity of both computational and observational astrophysics over the past decade necessitates efficient and intuitive methods for examining and visualizing large datasets. Here, I present AstroBlend, an open-source Python library for use within the three dimensional modeling software, Blender. While Blender has been a popular open-source software among animators and visual effects artists, in recent years it has also become a tool for visualizing astrophysical datasets. AstroBlend combines the three dimensional capabilities of Blender with the analysis tools of the widely used astrophysical toolset, yt, to afford both computational and observational astrophysicists the ability to simultaneously analyze their data and create informative and appealing visualizations. The introduction of this package includes a description of features, work flow, and various example visualizations.
- Automated detection of extended sources in radio maps: progress from the SCORPIO survey
- Automated source extraction and parametrization represents a crucial challenge for the next-generation radio interferometer surveys, such as those performed with the Square Kilometre Array (SKA) and its precursors. In this paper, we present a new algorithm, called CAESAR (Compact And Extended Source Automated Recognition), to detect and parametrize extended sources in radio interferometric maps. It is based on a pre-filtering stage, allowing image denoising, compact source suppression and enhancement of diffuse emission, followed by an adaptive superpixel clustering stage for final source segmentation. A parametrization stage provides source flux information and a wide range of morphology estimators for post-processing analysis. We developed CAESAR in a modular software library, also including different methods for local background estimation and image filtering, along with alternative algorithms for both compact and diffuse source extraction. The method was applied to real radio continuum data collected at the Australian Telescope Compact Array (ATCA) within the SCORPIO project, a pathfinder of the Evolutionary Map of the Universe (EMU) survey at the Australian Square Kilometre Array Pathfinder (ASKAP). The source reconstruction capabilities were studied over different test fields in the presence of compact sources, imaging artefacts and diffuse emission from the Galactic plane and compared with existing algorithms. When compared to a human-driven analysis, the designed algorithm was found capable of detecting known target sources and regions of diffuse emission, outperforming alternative approaches over the considered fields.
- A Novel Parallel Method for Speckle Masking Reconstruction Using the OpenMP
- High resolution reconstruction technology is developed to help enhance the spatial resolution of observational images for ground-based solar telescopes, such as speckle masking. Near real-time reconstruction performance is achieved on a high performance cluster using the Message Passing Interface (MPI). However, much time is spent in reconstructing solar subimages in such a speckle reconstruction. We design and implement a novel parallel method for speckle masking reconstruction of solar subimage on a shared memory machine using the OpenMP. Real tests are performed to verify the correctness of our codes. We present the details of several parallel reconstruction steps. The parallel implementation between various modules shows a great speed increase as compared to single thread serial implementation, and a speedup of about 2.5 is achieved in one subimage reconstruction. The timing result for reconstructing one subimage with 256×256 pixels shows a clear advantage with greater number of threads. This novel parallel method can be valuable in real-time reconstruction of solar images, especially after porting to a high performance cluster.
- SymPix: A Spherical Grid for Efficient Sampling of Rotationally Invariant Operators
- We present SymPix, a special-purpose spherical grid optimized for efficiently sampling rotationally invariant linear operators. This grid is conceptually similar to the Gauss-Legendre (GL) grid, aligning sample points with iso-latitude rings located on Legendre polynomial zeros. Unlike the GL grid, however, the number of grid points per ring varies as a function of latitude, avoiding expensive oversampling near the poles and ensuring nearly equal sky area per grid point. The ratio between the number of grid points in two neighboring rings is required to be a low-order rational number (3, 2, 1, 4/3, 5/4, or 6/5) to maintain a high degree of symmetries. Our main motivation for this grid is to solve linear systems using multi-grid methods, and to construct efficient preconditioners through pixel-space sampling of the linear operator in question. As a benchmark and representative example, we compute a preconditioner for a linear system that involves the operator \widehat{{\boldsymbol{D}}}+{\widehat{{\boldsymbol{B}}}}T{{\boldsymbol{N}}}-1\widehat{{\boldsymbol{B}}}, where \widehat{{\boldsymbol{B}}} and \widehat{{\boldsymbol{D}}} may be described as both local and rotationally invariant operators, and {\boldsymbol{N}} is diagonal in the pixel domain. For a bandwidth limit of {{\ell }}{max} = 3000, we find that our new SymPix implementation yields average speed-ups of 360 and 23 for {\widehat{{\boldsymbol{B}}}}T{{\boldsymbol{N}}}-1\widehat{{\boldsymbol{B}}} and \widehat{{\boldsymbol{D}}}, respectively, compared with the previous state-of-the-art implementation.
- Radio astronomical image formation using constrained least squares and Krylov subspaces
- Aims: Image formation for radio astronomy can be defined as estimating the spatial intensity distribution of celestial sources throughout the sky, given an array of antennas. One of the challenges with image formation is that the problem becomes ill-posed as the number of pixels becomes large. The introduction of constraints that incorporate a priori knowledge is crucial. Methods: In this paper we show that in addition to non-negativity, the magnitude of each pixel in an image is also bounded from above. Indeed, the classical “dirty image” is an upper bound, but a much tighter upper bound can be formed from the data using array processing techniques. This formulates image formation as a least squares optimization problem with inequality constraints. We propose to solve this constrained least squares problem using active set techniques, and the steps needed to implement it are described. It is shown that the least squares part of the problem can be efficiently implemented with Krylov-subspace-based techniques. We also propose a method for correcting for the possible mismatch between source positions and the pixel grid. This correction improves both the detection of sources and their estimated intensities. The performance of these algorithms is evaluated using simulations. Results: Based on parametric modeling of the astronomical data, a new imaging algorithm based on convex optimization, active sets, and Krylov-subspace-based solvers is presented. The relation between the proposed algorithm and sequential source removing techniques is explained, and it gives a better mathematical framework for analyzing existing algorithms. We show that by using the structure of the algorithm, an efficient implementation that allows massive parallelism and storage reduction is feasible. Simulations are used to compare the new algorithm to classical CLEAN. Results illustrate that for a discrete point model, the proposed algorithm is capable of detecting the correct number of sources and producing highly accurate intensity estimates.
- Fifty Years of Pulsar Candidate Selection: From simple filters to a new principled real-time classification approach
- Improving survey specifications are causing an exponential rise in pulsar candidate numbers and data volumes. We study the candidate filters used to mitigate these problems during the past fifty years. We find that some existing methods such as applying constraints on the total number of candidates collected per observation, may have detrimental effects on the success of pulsar searches. Those methods immune to such effects are found to be ill-equipped to deal with the problems associated with increasing data volumes and candidate numbers, motivating the development of new approaches. We therefore present a new method designed for on-line operation. It selects promising candidates using a purpose-built tree-based machine learning classifier, the Gaussian Hellinger Very Fast Decision Tree (GH-VFDT), and a new set of features for describing candidates. The features have been chosen so as to i) maximise the separation between candidates arising from noise and those of probable astrophysical origin, and ii) be as survey-independent as possible. Using these features our new approach can process millions of candidates in seconds (˜1 million every 15 seconds), with high levels of pulsar recall (90%+). This technique is therefore applicable to the large volumes of data expected to be produced by the Square Kilometre Array (SKA). Use of this approach has assisted in the discovery of 20 new pulsars in data obtained during the LOFAR Tied-Array All-Sky Survey (LOTAAS).
- A package for the automated classification of periodic variable stars
- We present a machine learning package for the classification of periodic variable stars. Our package is intended to be general: it can classify any single band optical light curve comprising at least a few tens of observations covering durations from weeks to years with arbitrary time sampling. We use light curves of periodic variable stars taken from OGLE and EROS-2 to train the model. To make our classifier relatively survey-independent, it is trained on 16 features extracted from the light curves (e.g., period, skewness, Fourier amplitude ratio). The model classifies light curves into one of seven superclasses – δ Scuti, RR Lyrae, Cepheid, Type II Cepheid, eclipsing binary, long-period variable, non-variable – as well as subclasses of these, such as ab, c, d, and e types for RR Lyraes. When trained to give only superclasses, our model achieves 0.98 for both recall and precision as measured on an independent validation dataset (on a scale of 0 to 1). When trained to give subclasses, it achieves 0.81 for both recall and precision. The majority of misclassifications of the subclass model is caused by confusion within a superclass rather than between superclasses. To assess classification performance of the subclass model, we applied it to the MACHO, LINEAR, and ASAS periodic variables, which gave recall/precision of 0.92/0.98, 0.89/0.96, and 0.84/0.88, respectively. We also applied the subclass model to Hipparcos periodic variable stars of many other variability types that do not exist in our training set, in order to examine how much those types degrade the classification performance of our target classes. In addition, we investigate how the performance varies with the number of data points and duration of observations. We find that recall and precision do not vary significantly if there are more than 80 data points and the duration is more than a few weeks. The classifier software of the subclass model is available (in Python) from the GitHub repository (http://https://goo.gl/xmFO6Q).
- CHIPS: The Cosmological H i Power Spectrum Estimator
- Detection of the cosmological neutral hydrogen signal from the Epoch of Reionization (EoR) and estimation of its basic physical parameters are principal scientific aims of many current low-frequency radio telescopes. Here we describe the Cosmological H i Power Spectrum Estimator (CHIPS), an algorithm developed and implemented with data from the Murchison Widefield Array, to compute the two-dimensional and spherically-averaged power spectrum of brightness temperature fluctuations. The principal motivations for CHIPS are the application of realistic instrumental and foreground models to form the optimal estimator, thereby maximizing the likelihood of unbiased signal estimation, and allowing a full covariant understanding of the outputs. CHIPS employs an inverse-covariance weighting of the data through the maximum likelihood estimator, thereby allowing use of the full parameter space for signal estimation (“foreground suppression”). We describe the motivation for the algorithm, implementation, application to real and simulated data, and early outputs. Upon application to a set of 3 hr of data, we set a 2σ upper limit on the EoR dimensionless power at k=0.05 {{h}} Mpc-1 of {{{Δ }}}k2\lt 7.6× {10}4 mK2 in the redshift range z = [6.2-6.6], consistent with previous estimates.
- A Computer-generated Visual Morphology Catalog of ~3,000,000 SDSS Galaxies
- We have applied computer analysis to classify the broad morphological types of ∼3 · 106 Sloan Digital Sky Survey (SDSS) galaxies. For each galaxy, the catalog provides the DR8 object ID, the R.A., the decl., and the certainty for the automatic classification as either spiral or elliptical. The certainty of the classification allows us to control the accuracy of a subset of galaxies by sacrificing some of the least certain classifications. The accuracy of the catalog was tested using galaxies that were classified by the manually annotated Galaxy Zoo catalog. The results show that the catalog contains ∼900,000 spiral galaxies and ∼600,000 elliptical galaxies with classification certainty that has a statistical agreement rate of ∼98% with the Galaxy Zoo debiased “superclean” data set. The catalog also shows that objects assigned by the SDSS pipeline with a relatively high redshift (z > 0.4) can have clear visual spiral morphology. The catalog can be downloaded at http://vfacstaff.ltu.edu/lshamir/data/morph_catalog. The image analysis software that was used to create the catalog is also publicly available.
- WISE Photometry for 400 Million SDSS Sources
- We present photometry of images from the Wide-Field Infrared Survey Explorer (WISE) of over 400 million sources detected by the Sloan Digital Sky Survey (SDSS). We use a “forced photometry” technique, using measured SDSS source positions, star-galaxy classification, and galaxy profiles to define the sources whose fluxes are to be measured in the WISE images. We perform photometry with The Tractor image modeling code, working on our “unWISE” coaddds and taking account of the WISE point-spread function and a noise model. The result is a measurement of the flux of each SDSS source in each WISE band. Many sources have little flux in the WISE bands, so often the measurements we report are consistent with zero given our uncertainties. However, for many sources we get 3σ or 4σ measurements; these sources would not be reported by the “official” WISE pipeline and will not appear in the WISE catalog, yet they can be highly informative for some scientific questions. In addition, these small-signal measurements can be used in stacking analyses at the catalog level. The forced photometry approach has the advantage that we measure a consistent set of sources between SDSS and WISE, taking advantage of the resolution and depth of the SDSS images to interpret the WISE images; objects that are resolved in SDSS but blended together in WISE still have accurate measurements in our photometry. Our results, and the code used to produce them, are publicly available at http://unwise.me.
- Asteroid orbits with Gaia using random-walk statistical ranging
- We describe statistical inverse methods for the computation of initial asteroid orbits within the data processing and analysis pipeline of the ESA Gaia space mission. Given small numbers of astrometric observations across short time intervals, we put forward a random-walk ranging method, in which the orbital-element phase space is uniformly sampled, up to a limiting χ2-value, with the help of the Markov-chain Monte Carlo technique (MCMC). The sample orbits obtain weights from the a posteriori probability density value and the MCMC rejection rate. For the first time, we apply the method to Gaia astrometry of asteroids. The results are nominal in that the method provides realistic estimates for the orbital uncertainties and meets the efficiency requirements for the daily, short-term processing of unknown objects.
- Cosmic Reionization on Computers: Numerical and Physical Convergence
- In this paper I show that simulations of reionization performed under the Cosmic Reionization On Computers project do converge in space and mass, albeit rather slowly. A fully converged solution (for a given star formation and feedback model) can be determined at a level of precision of about 20%, but such a solution is useless in practice, since achieving it in production-grade simulations would require a large set of runs at various mass and spatial resolutions, and computational resources for such an undertaking are not yet readily available. In order to make progress in the interim, I introduce a weak convergence correction factor in the star formation recipe, which allows one to approximate the fully converged solution with finite-resolution simulations. The accuracy of weakly converged simulations approaches a comparable, ~20% level of precision for star formation histories of individual galactic halos and other galactic properties that are directly related to star formation rates, such as stellar masses and metallicities. Yet other properties of model galaxies, for example, their H I masses, are recovered in the weakly converged runs only within a factor of 2.
- Odyssey: A Public GPU-based Code for General Relativistic Radiative Transfer in Kerr Spacetime
- General relativistic radiative transfer calculations coupled with the calculation of geodesics in the Kerr spacetime are an essential tool for determining the images, spectra, and light curves from matter in the vicinity of black holes. Such studies are especially important for ongoing and upcoming millimeter/submillimeter very long baseline interferometry observations of the supermassive black holes at the centers of Sgr A* and M87. To this end we introduce Odyssey, a graphics processing unit (GPU) based code for ray tracing and radiative transfer in the Kerr spacetime. On a single GPU, the performance of Odyssey can exceed 1 ns per photon, per Runge-Kutta integration step. Odyssey is publicly available, fast, accurate, and flexible enough to be modified to suit the specific needs of new users. Along with a Graphical User Interface powered by a video-accelerated display architecture, we also present an educational software tool, Odyssey_Edu, for showing in real time how null geodesics around a Kerr black hole vary as a function of black hole spin and angle of incidence onto the black hole.
- A chemical reaction network solver for the astrophysics code NIRVANA
- Context. Chemistry often plays an important role in astrophysical gases. It regulates thermal properties by changing species abundances and via ionization processes. This way, time-dependent cooling mechanisms and other chemistry-related energy sources can have a profound influence on the dynamical evolution of an astrophysical system. Modeling those effects with the underlying chemical kinetics in realistic magneto-gasdynamical simulations provide the basis for a better link to observations. Aims: The present work describes the implementation of a chemical reaction network solver into the magneto-gasdynamical code NIRVANA. For this purpose a multispecies structure is installed, and a new module for evolving the rate equations of chemical kinetics is developed and coupled to the dynamical part of the code. A small chemical network for a hydrogen-helium plasma was constructed including associated thermal processes which is used in test problems. Methods: Evolving a chemical network within time-dependent simulations requires the additional solution of a set of coupled advection-reaction equations for species and gas temperature. Second-order Strang-splitting is used to separate the advection part from the reaction part. The ordinary differential equation (ODE) system representing the reaction part is solved with a fourth-order generalized Runge-Kutta method applicable for stiff systems inherent to astrochemistry. Results: A series of tests was performed in order to check the correctness of numerical and technical implementation. Tests include well-known stiff ODE problems from the mathematical literature in order to confirm accuracy properties of the solver used as well as problems combining gasdynamics and chemistry. Overall, very satisfactory results are achieved. Conclusions: The NIRVANA code is now ready to handle astrochemical processes in time-dependent simulations. An easy-to-use interface allows implementation of complex networks including thermal processes. In combination with NIRVANA’s self-gravity solver, its efficient solver for dissipation terms and its adaptive mesh refinement capability challenging astrophysical problems come into reach with the code.
- A panoptic model for planetesimal formation and pebble delivery
- Context. The journey from dust particle to planetesimal involves physical processes acting on scales ranging from micrometers (the sticking and restructuring of aggregates) to hundreds of astronomical units (the size of the turbulent protoplanetary nebula). Considering these processes simultaneously is essential when studying planetesimal formation. Aims: The goal of this work is to quantify where and when planetesimal formation can occur as the result of porous coagulation of icy grains and to understand how the process is influenced by the properties of the protoplanetary disk. Methods: We develop a novel, global, semi-analytical model for the evolution of the mass-dominating dust particles in a turbulent protoplanetary disk that takes into account the evolution of the dust surface density while preserving the essential characteristics of the porous coagulation process. This panoptic model is used to study the growth from sub-micron to planetesimal sizes in disks around Sun-like stars. Results: For highly porous ices, unaffected by collisional fragmentation and erosion, rapid growth to planetesimal sizes is possible in a zone stretching out to ~10 AU for massive disks. When porous coagulation is limited by erosive collisions, the formation of planetesimals through direct coagulation is not possible, but the creation of a large population of aggregates with Stokes numbers close to unity might trigger the streaming instability (SI). However, we find that reaching conditions necessary for SI is difficult and limited to dust-rich disks, (very) cold disks, or disks with weak turbulence. Conclusions: Behind the snow-line, porosity-driven aggregation of icy grains results in rapid (~104 yr) formation of planetesimals. If erosive collisions prevent this, SI might be triggered for specific disk conditions. The numerical approach introduced in this work is ideally suited for studying planetesimal formation and pebble delivery simultaneously and will help build a coherent picture of the start of the planet formation process.
- Deep-HiTS: Rotation Invariant Convolutional Neural Network for Transient Detection
- We introduce Deep-HiTS, a rotation-invariant convolutional neural network (CNN) model for classifying images of transient candidates into artifacts or real sources for the High cadence Transient Survey (HiTS). CNNs have the advantage of learning the features automatically from the data while achieving high performance. We compare our CNN model against a feature engineering approach using random forests (RFs). We show that our CNN significantly outperforms the RF model, reducing the error by almost half. Furthermore, for a fixed number of approximately 2000 allowed false transient candidates per night, we are able to reduce the misclassified real transients by approximately one-fifth. To the best of our knowledge, this is the first time CNNs have been used to detect astronomical transient events. Our approach will be very useful when processing images from next generation instruments such as the Large Synoptic Survey Telescope. We have made all our code and data available to the community for the sake of allowing further developments and comparisons at https://github.com/guille-c/Deep-HiTS.
- An efficient and flexible Abel-inversion method for noisy data
- We propose an efficient and flexible method for solving the Abel integral equation of the first kind, frequently appearing in many fields of astrophysics, physics, chemistry, and applied sciences. This equation represents an ill-posed problem, thus solving it requires some kind of regularization. Our method is based on solving the equation on a so-called compact set of functions and/or using Tikhonov’s regularization. A priori constraints on the unknown function, defining a compact set, are very loose and can be set using simple physical considerations. Tikhonov’s regularization in itself does not require any explicit a priori constraints on the unknown function and can be used independently of such constraints or in combination with them. Various target degrees of smoothness of the unknown function may be set, as required by the problem at hand. The advantage of the method, apart from its flexibility, is that it gives uniform convergence of the approximate solution to the exact solution, as the errors of input data tend to zero. The method is illustrated on several simulated models with known solutions. An example of astrophysical application of the method is also given.
- Fast Algorithms for Demixing Sparse Signals From Nonlinear Observations
- We study the problem of demixing a pair of sparse signals from noisy, nonlinear observations of their superposition. Mathematically, we consider a nonlinear signal observation model, yi = g(aiTx) + ei, i = 1,…,m, where x = Φw + Ψz denotes the superposition signal, Φ and Ψ are orthonormal bases in ℝn, and w, z ∈ ℝn are sparse coefficient vectors of the constituent signals, and ei represents the noise. Moreover, g represents a nonlinear link function, and ai ∈ ℝn is the ith row of the measurement matrix A ∈ ℝm×n. Problems of this nature arise in several applications ranging from astronomy, computer vision, and machine learning. In this paper, we make some concrete algorithmic progress for the above demixing problem. Specifically, we consider two scenarios: first, the case when the demixing procedure has no knowledge of the link function, and second is the case when the demixing algorithm has perfect knowledge of the link function. In both cases, we provide fast algorithms for recovery of the constituents w and z from the observations. Moreover, we support these algorithms with a rigorous theoretical analysis and derive (nearly) tight upper bounds on the sample complexity of the proposed algorithms for achieving stable recovery of the component signals. We also provide a range of numerical simulations to illustrate the performance of the proposed algorithms on both real and synthetic signals and images.