xzackli / bolt.jl

differentiable boltzmann code

Bolt

⚡ Bolt is a pure-Julia integrator for the Boltzmann equations in cosmology. It can accurately compute the gradient of the CMB power spectrum, with respect to cosmological parameters, using forward-mode automatic differentiation.

Contributors: Jamie Sullivan, Zack Li, Marius Millea

Install

Bolt requires Julia 1.5+. To install, from the package prompt, run:

pkg> add https://github.com/xzackli/Bolt.jl

Gallery

A CMB temperature power spectrum and gradient from Bolt.jl.

A linear matter power spectrum and gradient from Bolt.jl.

bolt.jl's Issues

Clean up type tracking in background

There is no reason for us to enforce that things like the quadrature points/weights should be parameterized as the same type "T" that gets used for the dual. This happens in background.jl, quad_pts, quad_wts and in the ie branch where the splines are used for the iteration. Not sure if it will happen in the updated recfast.jl.
I can't imagine that doing this extra tracking is good for code performance when taking derivatives, though I am not sure how much of a difference it will make.
In any event, this is a quick fix, just putting this here so I don't forget

speed of sound term

From @jmsull,

When I look at the MB expression (eqn 67) for the baryon velocity, there is a speed of sound term that is proportional to delta_b that is missing in Callin/Dodelson (and Bolt). Shouldn't we have that?

I should get to the bottom of this...

Lensing

Marco has graciously provided an initial halofit implementation, we should also make some progress toward CMB lensing

AD backend comparison

Have long needed to do this in a quantitative way, so making an issue for it

why is RECFAST slow?

The RECFAST implementation currently in Bolt should be on the order of milliseconds. What's happening?

Photon-baryon decoupling at very early times?

This may just be basic thing I don't understand, but when I was making the movies for my COSMO talk, I noticed that baryons and photons actually decouple initially, then recouple. Is that right physically? I would have though at those early times theyre kept in extremely good equilibrium. Could that indicate an error in the code or the solver? Or maybe I set up the initial conditions wrong? (I changed \Phi initial per-mode, but everything else is consistent with that)

photons_neutrinos_baryons.mp4

My code is at https://github.com/xzackli/Bolt.jl/compare/main...marius311:cosmo21tweaks?expand=1 (not meant to be merged ever)

RECFAST

I believe this is a one-time painful ordeal that must be done -- I have to port RECFAST and test it. It's one-time because so-called "fudged RECFAST" parameters agree extremely well with codes like cosmorec or hyrec (@cmbant), and therefore this RECFAST implementation can supply adjoints for other recombination codes that I don't have to re-implement, just wrap. I think the Fortran is pretty readable, and the paper is too....

Method Error

I am new to Julia and am trying to use this code, but I keep getting the following error when I try to run RECFAST in all of the example code files. Please advise.

`MethodError: no method matching Bolt.RECFAST(::Background{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, Interpolations.ScaledInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, 1, Interpolations.BSplineInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, 1, OffsetArrays.OffsetVector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float6...
Closest candidates are:
Bolt.RECFAST(::AB, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::Any, ::Any, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T) where {T, AB<:(AbstractBackground{T, IT} where IT<:(Interpolations.AbstractInterpolation{T, 1}))} at ~/.julia/packages/Parameters/MK0O4/src/Parameters.jl:526
Bolt.RECFAST(::AB; kws...) where {T, AB<:(AbstractBackground{T, IT} where IT<:(Interpolations.AbstractInterpolation{T, 1}))} at ~/.julia/packages/Bolt/fGczu/src/ionization/recfast.jl:124

Stacktrace:
[1] Bolt.RECFAST(; bg::Background{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, Interpolations.ScaledInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, 1, Interpolations.BSplineInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, 1, OffsetArrays.OffsetVector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, C::Float64, k_B::Float64, h_P::Float64, m_e::Float64, m_H::Float64, not4::Float64, sigma::Float64, a::Float64, G::Float64, Lambda::Float64, Lambda_He::Float64, L_H_ion::Float64, L_H_alpha::Float64, L_He1_ion::Float64, L_He2_ion::Float64, L_He_2s::Float64, L_He_2p::Float64, A2P_s::Float64, A2P_t::Float64, L_He_2Pt::Float64, L_He_2St::Float64, L_He2St_ion::Float64, sigma_He_2Ps::Float64, sigma_He_2Pt::Float64, AGauss1::Float64, AGauss2::Float64, zGauss1::Float64, zGauss2::Float64, wGauss1::Float64, wGauss2::Float64, a_PPB::Float64, b_PPB::Float64, c_PPB::Float64, d_PPB::Float64, a_VF::Float64, b_VF::Float64, ...
@ Bolt ~/.julia/packages/Parameters/MK0O4/src/Parameters.jl:545
[2] pkc(Ω_c::ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1}, k_grid::Vector{Float64})
@ Main ./In[7]:4
[3] fc(Ω_c::ForwardDiff.Dual{ForwardDiff.Tag{typeof(fc), Float64}, Float64, 1})
@ Main ./In[8]:7
[4] derivative(f::typeof(fc), x::Float64)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/derivative.jl:14
[5] top-level scope
@ In[12]:1`

Rodas5 instability with current massive neutrinos

Following up from mass_nu PR:
"""
Regarding the instability you encountered -- let's turn that into an issue, I imagine it's easy to make an minimal demonstration. My current set of questions is

1. Does this happen if you remove the effect of massive neutrinos on the metric? (i.e. they're just free streaming the entire time)
2. Does this happen if you truncate the massive neutrinos to just the monopole / dipole?
3. Does this happen when you change the q grid?

However, we can deal with it later.

Originally posted by @xzackli in #39 (comment)
"""

Background speed

The background is clocking in at ~200ms in btime - that's too long. Most of the time is in the FD background integral, so we should think about how to optimize this better.

Code not working with Julia 1.10

It seems that recent updates to Julia and some of its common packages have broken some part of the code, rendering it unusable. I encountered the error while running basic_usage.jl. It seems that the error is originating from the DiffEqBase package, with error:

---

ERROR: LoadError: MethodError: anyeltypedual(::Type{Union{}}) is ambiguous.

Candidates:
anyeltypedual(::Type{T}) where T<:ForwardDiff.Dual
@ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/forwarddiff.jl:147
anyeltypedual(x::Type{T}) where T<:ForwardDiff.AbstractConfig
@ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/forwarddiff.jl:120
anyeltypedual(::Type{T}) where T<:(Tuple{Vararg{T, N}} where {N, T})
@ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/forwarddiff.jl:166
anyeltypedual(::Type{T}) where T<:Union{Set, AbstractArray}
@ DiffEqBase ~/.julia/packages/DiffEqBase/4084R/src/forwarddiff.jl:151

Possible fix, define
anyeltypedual(::Type{Union{}})

Two of my students are getting the same exact error, so it is not something special about my local installation. I am running julia 1.10.1 on MacOS 11.7.4. Reverting to Julia 1.8 fixes the issue, so it indicates that it has to do with some recent updates.

fix source function from relativistic neutrinos

Looks like #28 broke the source function integration, should fix.

need a fast array type for the component vector

Interacting with the component vector could be easier. The vector is [Θ..., Θᵖ..., 𝒩..., Φ, δ, v, δb, vb]. The relativistic components Θ, Θᵖ, 𝒩 are expanded in multipoles, and therefore are parts of the array indexed by multipole (typically 0 through ~8). @marius311 tried using ComponentArrays in main...marius311:componentarrays but incurred a 2x speed hit, maybe it's something related to memory layout?

A simple array type that just does what is currently done (offsetarray of views) would be a nice addition.

Update C_ell functions to use primordial power parameters

The plot_derivs_cl example was written before A and ns were included as parameters, so should be updated to include them - just need to add a line here that has p.A*(k/kpivot)^p.n.

Scalar C\ell accuracy test

We did the power spectrum - we should do C\ell_TT,_TE,_EE (since there are now functions for those). I will take a first pass at this without digging too deeply into precision settings.

Recfast splines complain if final redshift is not zero

If I run the following

𝕡 = CosmoParams()
xstart,xend = -15.0,-1.0
bg = Background(𝕡; x_grid=xstart:0.1:xend)
𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b)
ih = IonizationHistory(𝕣, 𝕡, bg)

when I get to the last line I get the following error:

"ERROR: BoundsError: attempt to access 141-element scale(interpolate(OffsetArray(::Array{Float64,1}, 0:142), BSpline(Cubic(Line(OnGrid())))), (-15.0:0.1:-1.0,)) with element type Float64 at index [-0.9980221350275645]"

which is coming from around here since zfinal=0.0 appears hardcoded.

Can be persuaded otherwise but it would be nice if we allowed arbitrary final redshift for the perturbations? I will continue using 0 for now since this is not a priority, but this is something to maybe come back to later.

Difference between \rho_P0, \rho_\sigma, and \rho0\scrM

As mentioned in #44, I noticed something weird when plotting perturbations.

When normalizing the massive neutrino perturbations for plotting I normalize by the value of the \rho_\sigma function (from perturbations.jl) acting on a vector of all ones at the particular redshift.

I find that this is the same as the density I get from the background spline (bg.\rho0\scrM) of the background \rho_P0 function.

However when I call the background \rho_P0 function from the plotting script I get a different value than the one returned by \rho_\sigma *a^-4 and the spline (both of which have the same value).
I find this especially weird since \rho_P0 is what is used to generate the spline.

Not sure if this issue will affect the evolution, but I will look into it more carefully.

Make hierarchy/integrator types less restrictive

In the boltsolve and related functions, there is no need to restrict to "BasicNewtonian()" - we can instead allow any integrator (and it would be good to do this for ie branch). Will take a pass back through and do this for the hierarchy when that branch is eventually merged (or at least caught up).

source function integrations without quadrature

The Bessel integrals required to integrate the source function are pretty painful. I should implement some of the harmonic approaches that work well at small scales, like

• power law decomposition / FFTLog type stuff (Assassi+17)
• quasi-discrete Hankel transform methods may work better for large scales which need ISW

Update parameter defaults

Following CAMB and CLASS, we should update Neff to the default 3.044.

We should really also change n to be not 1.0 (and pick some fiducial cosmology for the other $\Lambda$ CDM parameters too, like Planck).

This is not a priority, but opening an issue anyways

Bug in plin code

Not sure how I missed it before, but there are some errant factors of h in the plin function here. I will remove them.

w0wa

Easy feature to add - add CLP fluid dark energy w0wa to background.jl.

clean up background integration

The background integration should basically be instant compared to the perturbation evolution. We should remove all the quadgk and replace it with fixed grid quadrature.

Cosmological perturbations without the Boltzmann hierarchy

Might be worth reading and potentially deciding this would be something that could make Bolt even cooler (or maybe a spinoff package): https://arxiv.org/abs/2105.02887

Add reionization

As referenced in #41, we need to put reionization modeling in to get the correct sound speed and photon/baryon perturbations at low redshift.

This does not look that complicated based on looking at CAMB notes section III.

Clean up notation a bit

e.g. Omega_r should be Omega_gamma, which may cause a problem for people using Bolt

Compare to CAMB Pk

So far we have been comparing to CLASS (for perturbations b/c of gauge), but should also compare to CAMB on matter power spectrum.

Clarification on Boltzmann Equations in Cosmology

Hello, I came across your software that uses Boltzmann equations in cosmology, but I am not familiar with this term. Could you please provide some clarification on what Boltzmann equations are and how they are applied in cosmology? Additionally, I am curious about your use of the Bohm-Hatz equation to solve vortex formation interactions and how this approach differs from the conventional one. Specifically, how does it address the challenges posed by large-scale objects like galaxies that do not always obey Newton's second law? Thank you for your help and for creating such interesting software!

New RECFAST doesn't like dual-numbers

Dear Jamie and Zack

After the latest update I get the following error when I try to ForwardDiff.jl through RECFAST:

Got exception outside of a @test
  MethodError: no method matching Bolt.RECFAST(::Bolt.Background{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, Interpolations.ScaledInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, 1, Interpolations.BSplineInterpolation{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, 1, OffsetArrays.OffsetVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Int64, ::Int64, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::ForwardDiff.Dual{ForwardDiff.Tag{var"#lin_Bolt#19"{Vector{Float64}}, Float64}, Float64, 1}, ::Float64, ::Float64, ::Float64)

  Closest candidates are:
    Bolt.RECFAST(::AB, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::Any, ::Any, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T) where {T, AB<:(Bolt.AbstractBackground{T, IT} where IT<:(Interpolations.AbstractInterpolation{T, 1}))}

Code to replicate the error:

function get_Bolt_pk0(Ω_c)
    𝕡 = Bolt.CosmoParams(Ω_c = Ωc)
    bg = Background(𝕡)
    𝕣 = Bolt.RECFAST(bg=bg, Yp=𝕡.Y_p, OmegaB=𝕡.Ω_b,OmegaG=𝕡.Ω_r)
    ih = IonizationHistory(𝕣, 𝕡, bg)
    # bg.H₀ = 0.00023349 [1/Mpc] 
    kmin,kmax,nk = 0.5*0.00023349,30000*0.00023349,70 
    ks_Bolt = log10_k(kmin,kmax,nk)
    pL = [plin(k,𝕡,bg,ih) for k in ks_Bolt]
    return ks_Bolt, pL
end
ForwardDiff.derivative(get_Bolt_pk0, 0.25)

(not exactly what I am running but should the trick)

The speed up is amazing though!!!

Thanks a lot,
Jaime

zero polarization in RSA

Right now there is no implementation of photon polarization multipoles during RSA:

Bolt.jl/src/perturbations.jl

Line 231 in de08e4d

Θᵖ′[:] = zeros(ℓᵧ+1)

This is definitely not right - these should at least be zeroed out so that we neglect them rather than preserving whatever the final values were before turning on RSA.

If not zeroing them, we can somehow relate them to the photon RSA expressions? It is not obvious to me if they do this in the CLASS 2 paper, but we should be able to figure something out.

RSA issues

The RSA trigger is set to the prescription given in the CLASS 2 paper - but I found for k=1h/Mpc, the RSA turns on too late, and the photon monopole has already started exploding. This (and a related bug I fixed in ie branch) probably is causing the residual large scale-dependent disagreements in the radiation transfer functions shown in #44, and (maybe?) could be contributing significantly to the 0.1% scale-dependent deviations in $P_L(k)$.

GPUs

We want this code to run on GPUs. We can try the easiest thing first - just passing CUDA arrays and then get more involved as necessary.

Strange early-time evolution in metric perturbations

In trying to match the CLASS perturbations in #44, we saw before that there was a several percent disagreement pretty much all the time (i.e. even on superhorizon scales / early times). I have looked at this more carefully, and \Phi is clearly doing something strange at very early times (plot attached for k=0.03 h/Mpc) that we were not looking at in that draft PR thread.

The orange dashed line is Bolt and shows that at early times \Phi is evolving in the first 100 steps or so before leveling off until horizon entry. This definitely should not be happening - \Phi (and \Psi) should be constant (and \Phi' negligible) until horizon entry (e.g. MB Fig.2).

You can actually see that initially (very first step) we are in agreement with the initial \Phi value of CLASS, but after this evolution, our \Phi is too high and this affects subsequent evolution, resulting in the percent level differences we see (I am surprised it is not worse, to be honest).

I tried setting Neff to some very small value (1e-4) and this seems to make the problem go away (though late-time evolution won't match CLASS because Neff is wrong):

This is also true if I take out both massive neutrinos and massless neutrinos (by commenting out the lines in the background and einstein equations). The issue is not related to the massive neutrinos, as it persists when only massless neutrinos are included.

I likely made some mistake in the past when setting up the massless neutrinos. I will dig into this in the first few steps to see what about massless neutrinos is making \Phi' non-negligible. I have checked that we get the correct \Psi in the initial conditions, so I guess it is somehow an error in the 00 Einstein eqn (or something less obvious with the neutrinos).

Bessel function integration for zero argument

Not sure this will ever actually be an issue in practice, but right now the bessel function interpolator is constructed using a point at x =0 (were we have j_\ell(x)) when \eta = \eta_final.

At least for the polarization source function, there is a problem because j_2(x)/x^2 |x=0 is not zero, but the way the source function integration is done is by calling the bessel interpolator at x=0 and then multiplying it into the source function.

I caught a NaN in the source function for x=0, but the bessel interpolator returning j_2(0)=0 will still be wrong.

This might be solved by the new Bessel integration?

photon boundary condition / RSA

As discussed in #44 by @jmsull, we appear to need something like CLASS's RSA to avoid unphysical oscillations from reflections off the boundary condition.

CLASS paper II has a description in the synchronous gauge, and Doran+05 has a similar approach for the Newtonian gauge. Taking the equations from Doran is the lowest effort option (no need to transform any variables), but I definitely need to read these papers in more detail.

More tests

Turn at least some of the perturbation-level and Pk comparisons with CLASS in #44 into tests.

Enzyme isn't ready for use with Bolt

I played with Enzyme a little bit, and I suspect it's not ready for use with our package. It can't differentiate simple ODEs at present. There's a fair amount of linear algebra in OrdinaryDiffEq, so it's probably tripping up on BLAS.

# this will crash
using OrdinaryDiffEq
using Enzyme

f(u,p,t) = 1.01*u
function test(u0)
    tspan = (0.0,1.0)
    prob = ODEProblem(f,u0,tspan)
    sol = solve(prob,Rodas4(),reltol=1e-8,abstol=1e-8)
    return sol(1.0)
end

autodiff(test, Active(1.0)) # (g, 1.0)

Scalar polarization C_ells

Add C^TE and C^EE to existing C^TT - should be quick - just need to add the (much simpler) E source function and re-use/generalise the current TT los integration function.

[compat] disabled until first tagged release

I'm disabling the compat helper and removing all compat entries. Let's add them back when we tag our first release, as things are pretty fluid right now. I don't want to worry about [compat] until we have a matter power spectrum we can write papers with.

Minor conformal time integration issue

At least on the ie branch, when I run the conformal time ODE solve, I find
┌ Warning: dt(1.8189894035458565e-12) <= dtmin(1.8189894035458565e-12) at t=0.0006687152338204658. Aborting. There is either an error in your model specification or the true solution is unstable. └ @ SciMLBase ~/.julia/packages/SciMLBase/MXZxL/src/integrator_interface.jl:529

But only the first time - on a second (or any subequent) solve, the solver runs fine. This seems like something we want to avoid?

I saw Chris Rackauckas mention something about this in Julia diff-eq slack a while ago, so I am tempted to think this is not our fault - will dig it up here.