Note: this is a alpha software.
Cosmology.jl is a code written in Julia to compute various cosmological quantities like angular diameter distances, (WIP:) matter power spectra, or (TODO:) the CMB anisotropy spectra.
Requirements:
Installation (from the package prompt):
pkg> add https://github.com/marius311/Cosmology.jlExample:
julia> using Cosmology
julia> p = new_params(H0=67, ωb=0.0225, ωc=0.12, mν=0.06eV, Nν_massive=1, Nν_massless=2.046, reltol=1e-4);
julia> rs(p,zdrag(p))/Mpc # sound horizon at baryon drag redshift
146.94167579179634
# etc..Things currently handled correctly:
The field of Cosmology already has some pretty mature and widely used codes which perform these calculations, CAMB and CLASS. So why Cosmology.jl?
Cosmology.jl can be faster. Julia is JIT compiled to native code, so it can be just as fast as C or Fortran, or faster. As a simple example where the codes perform essentially identical mathematical computations, the current calculation of the angular diameter distance is 7x faster than CAMB. This is encouraging and I see no reason more complex calculations can't have similar performance.
The code is far more readable. These calculations shouldn't be a black box only a few experts understand. Julia is built to write clear scientific code (for example, the aforemention angular diameter distance calculation is roughly 10x fewer lines of code than either CAMB or CLASS's, and still faster). A few specific things that really help:
Unicode variable names. This may seem silly but it makes a big differnce in reading code when you're familiar with the equations but not with the code-base. And Julia has many more Unicode characters available than Python, so this works really well. E.g., a code snippet:
function add_derived!()
ργ₀ = (π²/15)*Tγ₀^4
h² = (H0/100)^2
ωk = Ωk*h²
if mν == 0
Nν_massless += Nν_massive
Nν_massive = 0
end
ων = ρν(0)/ρx_over_ωx
ωΛ = h² - ωk - ωb - ωc - ων - ργ₀/ρx_over_ωx
endDAE solver. Lots of these calculations involve solving ordinary differential equations (ODEs). Traditional ODE solvers require you to manipulate your differential equations into the form y′ = F(y, x), which often puts them in an unfamiliar form far removed from the physically motivated way in which they are usually written down. We use instead a "differential algebraic solver" (DAE) which lets you write any general system of differential equations F(y′,y,x)=0. This (coupled with Unicode variable names) allows for writing very clear code, for example from the matter power spectrum calculation:
@eqs begin
#Photons
δγ′ == -4/3*θγ + 4ϕ′
θγ′ == k²/4*δγ + k²*ψ + a*nₑ*σT*(θb-θγ)
#CDM
δc′ == -θc + 3ϕ′
θc′ == -(a′/a)*θc + k²*ψ
#Baryons
δb′ == -θb + 3ϕ′
θb′ == -(a′/a)*θb + cₛ²*k²*δb + a*nₑ*σT/R*(θγ-θb) + k²*ψ
#Einstein equations
k²*ϕ + 3*(a′/a)*(ϕ′+(a′/a)*ψ) == -4π*a^2*(ρc(z)*δc + ρb(z)*δb + (ργ(z)+ρν(z))*δγ)
endPlanck units Absolutely every quantity in Cosmology.jl is in Planck units, so there's never any confusion, and equations are as simple as possible. (Other codes don't do this in an attempt to keep native 32bit support, since 32bit floats don't have a big enough exponent for expressing cosmological quantities in Planck units; we drop native 32bit support, i.e. you can run on 32bit, but its slower)
Opens up the possibility GPU acceleration and automatic differentiation. Julia has great support for automatic differentiation, which works through almost any arbitrary Julia code. Something like ForwardDiff.jl should work almost out-of-the-box and provide exact gradients in O(1) times the cost of running the calculation itself. Not only can this be used to compute Fisher matrices or for analysis techniques that require likelihood gradients (HMC, VI, etc...), but exact gradients can also be used to speed up the internal ODE solvers too. Additionally, Julia has great GPU integration, and using something like CUDA.jl or KernelAbstractions.jl one could offload much of the computation to GPU almost transparently, and without needing to write anything beside Julia code.