Julia Implementation Research Articles

Codebase release r0.3 for SmoQyDQMC.jl

We introduce the SmoQyDQMC.jl package, a Julia implementation of the determinant quantum Monte Carlo algorithm. SmoQyDQMC.jl supports generalized tight-binding Hamiltonians with on-site Hubbard and generalized electron-phonon (ee-ph) interactions, including non-linear ee-ph coupling and anharmonic lattice potentials. Our implementation uses hybrid Monte Carlo methods with exact forces for sampling the phonon fields, enabling efficient simulation of low-energy phonon branches, including acoustic phonons. The SmoQyDQMC.jl package also uses a flexible scripting interface, allowing users to adapt it to different workflows and interface with other software packages in the Julia ecosystem. The code for this package can be downloaded from our GitHub repository at https://github.com/SmoQySuite/SmoQyDQMC.jl or installed using the Julia package manager. The online documentation, including examples, can be obtained from our document page at https://smoqysuite.github.io/SmoQyDQMC.jl/stable/.

SmoQyDQMC.jl: A flexible implementation of determinant quantum Monte Carlo for Hubbard and electron-phonon interactions

We introduce the SmoQyDQMC.jl package, a Julia implementation of the determinant quantum Monte Carlo algorithm. SmoQyDQMC.jl supports generalized tight-binding Hamiltonians with on-site Hubbard and generalized electron-phonon (ee-ph) interactions, including non-linear ee-ph coupling and anharmonic lattice potentials. Our implementation uses hybrid Monte Carlo methods with exact forces for sampling the phonon fields, enabling efficient simulation of low-energy phonon branches, including acoustic phonons. The SmoQyDQMC.jl package also uses a flexible scripting interface, allowing users to adapt it to different workflows and interface with other software packages in the Julia ecosystem. The code for this package can be downloaded from our GitHub repository at https://github.com/SmoQySuite/SmoQyDQMC.jl or installed using the Julia package manager. The online documentation, including examples, can be obtained from our document page at https://smoqysuite.github.io/SmoQyDQMC.jl/stable/.

Singularity swap quadrature for nearly singular line integrals on closed curves in two dimensions

This paper presents a quadrature method for evaluating layer potentials in two dimensions close to periodic boundaries, discretized using the trapezoidal rule. It is an extension of the method of singularity swap quadrature, which recently was introduced for boundaries discretized using composite Gauss–Legendre quadrature. The original method builds on swapping the target singularity for its preimage in the complexified space of the curve parametrization, where the source panel is flat. This allows the integral to be efficiently evaluated using an interpolatory quadrature with a monomial basis. In this extension, we use the target preimage to swap the singularity to a point close to the unit circle. This allows us to evaluate the integral using an interpolatory quadrature with complex exponential basis functions. This is well-conditioned, and can be efficiently evaluated using the fast Fourier transform. The resulting method has exponential convergence, and can be used to accurately evaluate layer potentials close to the source geometry. We report experimental results on a simple test geometry, and provide a baseline Julia implementation that can be used for further experimentation.

Codebase release 1.0 for Nevanlinna.jl

We introduce a Julia implementation of the recently proposed Nevanlinna analytic continuation method. The method is based on Nevanlinna interpolants and inherently preserves the causality of a response function due to its construction. For theoretical calculations without statistical noise, this continuation method is a powerful tool to extract real-frequency information from numerical input data on the Matsubara axis. This method has been applied to first-principles calculations of correlated materials. This paper presents its efficient and full-featured open-source implementation of the method including the Hamburger moment problem and smoothing.

Nevanlinna.jl: A Julia implementation of Nevanlinna analytic continuation

We introduce a Julia implementation of the recently proposed Nevanlinna analytic continuation method. The method is based on Nevanlinna interpolants and inherently preserves the causality of a response function due to its construction. For theoretical calculations without statistical noise, this continuation method is a powerful tool to extract real-frequency information from numerical input data on the Matsubara axis. This method has been applied to first-principles calculations of correlated materials. This paper presents its efficient and full-featured open-source implementation of the method including the Hamburger moment problem and smoothing.

Computing Characteristic Polynomials of Hyperplane Arrangements with Symmetries

We introduce a new algorithm computing the characteristic polynomials of hyperplane arrangements which exploits their underlying symmetry groups. Our algorithm counts the chambers of an arrangement as a byproduct of computing its characteristic polynomial. We showcase our julia implementation, based on OSCAR, on examples coming from hyperplane arrangements with applications to physics and computer science.

Three-Phase Power Flow Tool for Electric Distribution Grids: A Julia Implementation for Electrical Engineering Students

Most of the concepts involved in Electrical Engineering have been developed from the 20th century to the present day. One of the most popular studies in this field is the power flow problem for electrical systems composed of multiple nonlinear loads (i.e., loads with constant power consumption) [1].

Branched Latent Neural Maps

We introduce Branched Latent Neural Maps (BLNMs) to learn finite dimensional input–output maps encoding complex physical processes. A BLNM is defined by a simple and compact feedforward partially-connected neural network that structurally disentangles inputs with different intrinsic roles, such as the time variable from model parameters of a differential equation, while transferring them into a generic field of interest. BLNMs leverage latent outputs to enhance the learned dynamics and break the curse of dimensionality by showing excellent in-distribution generalization properties with small training datasets and short training times on a single processor. Indeed, their in-distribution generalization error remains comparable regardless of the adopted discretization during the testing phase. Moreover, the partial connections, in place of a fully-connected structure, significantly reduce the number of tunable parameters. We show the capabilities of BLNMs in a challenging test case involving biophysically detailed electrophysiology simulations in a biventricular cardiac model of a pediatric patient with hypoplastic left heart syndrome. The model includes a 1D Purkinje network for fast conduction and a 3D heart-torso geometry. Specifically, we trained BLNMs on 150 in silico generated 12-lead electrocardiograms (ECGs) while spanning 7 model parameters, covering cell-scale, organ-level and electrical dyssynchrony. Although the 12-lead ECGs manifest very fast dynamics with sharp gradients, after automatic hyperparameter tuning the optimal BLNM, trained in less than 3 h on a single CPU, retains just 7 hidden layers and 19 neurons per layer. The resulting mean square error is on the order of 10−4 on an independent test dataset comprised of 50 additional electrophysiology simulations. In the online phase, the BLNM allows for 5000x faster real-time simulations of cardiac electrophysiology on a single core standard computer and can be employed to solve inverse problems via global optimization in a few seconds of computational time. This paper provides a novel computational tool to build reliable and efficient reduced-order models for digital twinning in engineering applications. The Julia implementation is publicly available under MIT License at https://github.com/StanfordCBCL/BLNM.jl.

BondGraphs.jl: Composable energy-based modelling in systems biology.

BondGraphs.jl is a Julia implementation of bond graphs. Bond graphs provide a modelling framework that describes energy flow through a physical system and by construction enforce thermodynamic constraints. The framework is widely used in engineering and has recently been shown to be a powerful approach for modelling biology. Models are mutable, hierarchical, multi-scale, multi-physics, and BondGraphs.jl is compatible with the Julia modelling ecosystem. BondGraphs.jl is freely available under the MIT license. Source code and documentation can be found at https://github.com/jedforrest/BondGraphs.jl. Supplementary data are available at Bioinformatics online.

A framework for synthetic power system dynamics.

We present a modular framework for generating synthetic power grids that consider the heterogeneity of real power grid dynamics but remain simple and tractable. This enables the generation of large sets of synthetic grids for a wide range of applications. For the first time, our synthetic model also includes the major drivers of fluctuations on short-time scales and a set of validators that ensure the resulting system dynamics are plausible. The synthetic grids generated are robust and show good synchronization under all evaluated scenarios, as should be expected for realistic power grids. A software package that includes an efficient Julia implementation of the framework is released as a companion to the paper.

A New ODE-Based Julia Implementation of the Anaerobic Digestion Model No. 1 Greatly Outperforms Existing DAE-Based Java and Python Implementations

The Anaerobic Digestion Model 1 is the quasi-industry standard for modelling anaerobic digestion, and it has seen several new implementations in recent years. It is assumed that these implementations would give the same results; however, a thorough comparison of these implementations has never been reported. This paper considers four different implementations of ADM1: one in Julia, one in Java, and two in Python. The Julia code is a de novo implementation of the ODE formulation of ADM1 that is reported here for the first time. The existing Java and Python codes implement the more common DAE formulation. Therefore, this paper also examines how DAE implementations compare to ODE implementations in terms of computational speed as well as solutions returned. As expected, the ODE and DAE forms both return comparable solutions. However, contrary to popular belief, the Julia ODE implementation is faster than the DAE implementations, namely by one to three orders of magnitude of compute time, depending on the simulation scenario and the reference implementation used for comparison.

Training End-to-End Unrolled Iterative Neural Networks for SPECT Image Reconstruction.

Training end-to-end unrolled iterative neural networks for SPECT image reconstruction requires a memory-efficient forward-backward projector for efficient backpropagation. This paper describes an open-source, high performance Julia implementation of a SPECT forward-backward projector that supports memory-efficient backpropagation with an exact adjoint. Our Julia projector uses only ~5% of the memory of an existing Matlab-based projector. We compare unrolling a CNN-regularized expectation-maximization (EM) algorithm with end-to-end training using our Julia projector with other training methods such as gradient truncation (ignoring gradients involving the projector) and sequential training, using XCAT phantoms and virtual patient (VP) phantoms generated from SIMIND Monte Carlo (MC) simulations. Simulation results with two different radionuclides (90Y and 177Lu) show that: 1) For 177Lu XCAT phantoms and 90Y VP phantoms, training unrolled EM algorithm in end-to-end fashion with our Julia projector yields the best reconstruction quality compared to other training methods and OSEM, both qualitatively and quantitatively. For VP phantoms with 177Lu radionuclide, the reconstructed images using end-to-end training are in higher quality than using sequential training and OSEM, but are comparable with using gradient truncation. We also find there exists a trade-off between computational cost and reconstruction accuracy for different training methods. End-to-end training has the highest accuracy because the correct gradient is used in backpropagation; sequential training yields worse reconstruction accuracy, but is significantly faster and uses much less memory.

Reliable genotyping of recombinant genomes using a robust hidden Markov model.

Meiotic recombination is an essential mechanism during sexual reproduction and includes the exchange of chromosome segments between homologous chromosomes. New allelic combinations are transmitted to the new generation, introducing novel genetic variation in the offspring genomes. With the improvement of high-throughput whole-genome sequencing technologies, large numbers of recombinant individuals can now be sequenced with low sequencing depth at low costs, necessitating computational methods for reconstructing their haplotypes. The main challenge is the uncertainty in haplotype calling that arises from the low information content of a single genomic position. Straightforward sliding window-based approaches are difficult to tune and fail to place recombination breakpoints precisely. Hidden Markov model (HMM)-based approaches, on the other hand, tend to over-segment the genome. Here, we present RTIGER, an HMM-based model that exploits in a mathematically precise way the fact that true chromosome segments typically have a certain minimum length. We further separate the task of identifying the correct haplotype sequence from the accurate placement of haplotype borders, thereby maximizing the accuracy of border positions. By comparing segmentations based on simulated data with known underlying haplotypes, we highlight the reasons for RTIGER outperforming traditional segmentation approaches. We then analyze the meiotic recombination pattern of segregants of 2 Arabidopsis (Arabidopsis thaliana) accessions and a previously described hyper-recombining mutant. RTIGER is available as an R package with an efficient Julia implementation of the core algorithm.

Unsupervised discovery of ancestry-informative markers and genetic admixture proportions in biobank-scale datasets

Admixture estimation plays a crucial role in ancestry inference and genome-wide association studies (GWASs). Computer programs such as ADMIXTURE and STRUCTURE are commonly employed to estimate the admixture proportions of sample individuals. However, these programs can be overwhelmed by the computational burdens imposed by the 105 to 106 samples and millions of markers commonly found in modern biobanks. An attractive strategy is to run these programs on a set of ancestry-informative SNP markers (AIMs) that exhibit substantially different frequencies across populations. Unfortunately, existing methods for identifying AIMs require knowing ancestry labels for a subset of the sample. This supervised learning approach creates a chicken and the egg scenario. In this paper, wepresent an unsupervised, scalable framework that seamlessly carries out AIM selection and likelihood-based estimation of admixture proportions.Our simulated and real data examples show that this approach is scalable to modern biobank datasets. OpenADMIXTURE, our Julia implementation of the method, is open source and available for free.

Libdlr: Efficient imaginary time calculations using the discrete Lehmann representation

We introduce libdlr, a library implementing the recently introduced discrete Lehmann representation (DLR) of imaginary time Green's functions. The DLR basis consists of a collection of exponentials chosen by the interpolative decomposition to ensure stable and efficient recovery of Green's functions from imaginary time or Matsubara frequency samples. The library provides subroutines to build the DLR basis and grids, and to carry out various standard operations. The simplicity of the DLR makes it straightforward to incorporate into existing codes as a replacement for less efficient representations of imaginary time Green's functions, and libdlr is intended to facilitate this process. libdlr is written in Fortran, provides a C header interface, and contains a Python module pydlr. We also introduce a stand-alone Julia implementation, Lehmann.jl. Program summaryProgram Title: libdlrCPC Library link to program files:https://doi.org/10.17632/56z594pzsj.1Developer's repository link:https://github.com/jasonkaye/libdlrLicensing provisions: Apache-2.0Programming language: Fortran, C, Python, JuliaNature of problem: Discretization and compression of functions (Green's functions and self-energies) with an imaginary time variable.Solution method: Explicit basis functions and discretization points obtained by low rank compression of the analytical continuation kernel.

Accelerating advances in landscape connectivity modelling with the ConScape library

AbstractIncreasingly precise spatial data (e.g. high‐resolution imagery from remote sensing) allow for improved representations of the landscape network for assessing the combined effects of habitat loss and connectivity declines on biodiversity. However, evaluating large landscape networks presents a major computational challenge both in terms of working memory and computation time. We present theConScape(i.e. “connected landscapes”) software library implemented in the high‐performance open‐source Julia language to compute metrics for connected habitat and movement flow on high‐resolution landscapes.The combination of Julia's ‘just‐in‐time’ compiler, efficient algorithms and ‘landmarks’ to reduce the computational load allowsConScapeto compute landscape ecological metrics—originally developed in metapopulation ecology (such as ‘metapopulation capacity’ and ‘probability of connectivity’)—for large landscapes. An additional major innovation inConScapeis the adoption of the randomized shortest paths framework to represent connectivity along the continuum from optimal to random movements, instead of only those extremes.We demonstrateConScape's potential for using large datasets in sustainable land planning by modelling landscape connectivity based on remote‐sensing data paired with GPS tracking of wild reindeer in Norway. To guide users, we discuss other applications, and provide a series of worked examples to showcase allConScape's functionalities in Supplementary Material.Built by a team of ecologists, network scientists and software developers,ConScapeis able to efficiently compute landscape metrics for high‐resolution landscape representations to leverage the availability of large data for sustainable land use and biodiversity conservation. As a Julia implementation,ConScapecombines computational efficiency with a transparent code base, which facilitates continued innovation through contributions from the rapidly growing community of landscape and connectivity modellers using Julia.

Equations for GL Invariant Families of Polynomials

We provide an algorithm that takes as an input a given parametric family of homogeneous polynomials, which is invariant under the action of the general linear group, and an integer d. It outputs the ideal of that family intersected with the space of homogeneous polynomials of degree d. Our motivation comes from Question 7 in Ranestad and Sturmfels (Le Math. 75, 411–424, 2020) and Problem 13 in Sturmfels (2014), which ask to find equations for varieties of cubic and quartic symmetroids. The algorithm relies on a database of specific Young tableaux and highest weight polynomials. We provide the database and the implementation of the database construction algorithm. Moreover, we provide a Julia implementation to run the algorithm using the database, so that more varieties of homogeneous polynomials can easily be treated in the future.

Safe-by-default Concurrency for Modern Programming Languages

Modern “safe” programming languages follow a design principle that we call safety by default and performance by choice . By default, these languages enforce important programming abstractions, such as memory and type safety, but they also provide mechanisms that allow expert programmers to explicitly trade some safety guarantees for increased performance. However, these same languages have adopted the inverse design principle in their support for multithreading. By default, multithreaded programs violate important abstractions, such as program order and atomic access to individual memory locations to admit compiler and hardware optimizations that would otherwise need to be restricted. Not only does this approach conflict with the design philosophy of safe languages, but very little is known about the practical performance cost of providing a stronger default semantics. In this article, we propose a safe-by-default and performance-by-choice multithreading semantics for safe languages, which we call volatile -by-default . Under this semantics, programs have sequential consistency (SC) by default, which is the natural “interleaving” semantics of threads. However, the volatile -by-default design also includes annotations that allow expert programmers to avoid the associated overheads in performance-critical code. We describe the design, implementation, optimization, and evaluation of the volatile -by-default semantics for two different safe languages: Java and Julia. First, we present V BD-HotSpot and V BDA-HotSpot, modifications of Oracle’s HotSpot JVM that enforce the volatile -by-default semantics on Intel x86-64 hardware and ARM-v8 hardware. Second, we present S C-Julia, a modification to the just-in-time compiler within the standard Julia implementation that provides best-effort enforcement of the volatile -by-default semantics on x86-64 hardware for the purpose of performance evaluation. We also detail two different implementation techniques: a baseline approach that simply reuses existing mechanisms in the compilers for handling atomic accesses, and a speculative approach that avoids the overhead of enforcing the volatile -by-default semantics until there is the possibility of an SC violation. Our results show that the cost of enforcing SC is significant but arguably still acceptable for some use cases today. Further, we demonstrate that compiler optimizations as well as programmer annotations can reduce the overhead considerably.

COSMO: A Conic Operator Splitting Method for Convex Conic Problems

This paper describes the conic operator splitting method (COSMO) solver, an operator splitting algorithm and associated software package for convex optimisation problems with quadratic objective function and conic constraints. At each step, the algorithm alternates between solving a quasi-definite linear system with a constant coefficient matrix and a projection onto convex sets. The low per-iteration computational cost makes the method particularly efficient for large problems, e.g. semidefinite programs that arise in portfolio optimisation, graph theory, and robust control. Moreover, the solver uses chordal decomposition techniques and a new clique merging algorithm to effectively exploit sparsity in large, structured semidefinite programs. Numerical comparisons with other state-of-the-art solvers for a variety of benchmark problems show the effectiveness of our approach. Our Julia implementation is open source, designed to be extended and customised by the user, and is integrated into the Julia optimisation ecosystem.

Julia subtyping: a rational reconstruction

Programming languages that support multiple dispatch rely on an expressive notion of subtyping to specify method applicability. In these languages, type annotations on method declarations are used to select, out of a potentially large set of methods, the one that is most appropriate for a particular tuple of arguments. Julia is a language for scientific computing built around multiple dispatch and an expressive subtyping relation. This paper provides the first formal definition of Julia's subtype relation and motivates its design. We validate our specification empirically with an implementation of our definition that we compare against the existing Julia implementation on a collection of real-world programs. Our subtype implementation differs on 122 subtype tests out of 6,014,476. The first 120 differences are due to a bug in Julia that was fixed once reported; the remaining 2 are under discussion.