Automatic differentiation

This page shows how to use Automatic Differentiation in combination with an EPG simulation.

The AD package tested is ForwardDiff.jl, maybe it works with others with some minor modification to the following code.

Load package

using EPGsim, ForwardDiff, CairoMakie

Building signal function

function MESE_EPG(T2,T1,TE,ETL,delta)
  T = eltype(complex(T2))
  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
  echo_vec = Vector{Complex{eltype(T2)}}()

  E = epgRotation(E,pi/2*delta, pi/2)
  # loop over refocusing-pulses
  for i = 1:ETL
    E = epgDephasing(E,1)
    E = epgRelaxation(E,TE,T1,T2)
    E  = epgRotation(E,pi*delta,0.0)
    E  = epgDephasing(E,1)
    push!(echo_vec,E.Fp[1])
  end

  return abs.(echo_vec)
end;
MESE_EPG (generic function with 1 method)
Specific types with AD

ForwardDiff use a specific type : Dual <: Real. The target function must be written generically enough to accept numbers of type T<:Real as input (or arrays of these numbers).

We also need to create an EPGStates that is of that type. We need to force it to be complex :

T = eltype(complex(T2))
E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])

Define parameters for simulation and run it

T2 = 60.0
T1 = 1000.0
TE = 7
ETL = 50
deltaB1 = 1

TE_vec = range(TE,TE*ETL,ETL)

amp = MESE_EPG(T2,T1,TE,ETL,deltaB1)
lines(TE_vec,abs.(amp),axis =(;title = "MESE Signal", xlabel="TE [ms]"))

As expected, we get a standard T2 decaying exponential curve :

\[S(TE) = exp(-TE/T_2)\]

We can analytically derive the equation according to $T_2$ :

\[\frac{\partial S}{\partial T_2} = \frac{TE}{T_2^2} exp(-TE/T_2)\]

which give the following curves:

df = TE_vec .* exp.(-TE_vec./T2)./(T2^2)

lines(TE_vec,abs.(df),axis =(;title = "dS/dT2", xlabel="TE [ms]"))

Find the derivative with Automatic Differentiation

Because we want to obtain the derivate at multiple time points (TE), we will use ForwardDiff.jacobian :

j = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])

Let's compare it to the analytical equation :

f=Figure()
ax = Axis(f[1,1],title ="Analytic vs Automatic Differentiation")

lines!(ax,TE_vec,abs.(df),label = "Analytic Differentiation",linewidth=3)
lines!(ax,TE_vec,abs.(vec(j)),label = "Automatic Differentiation",linestyle=:dash,linewidth=3)
axislegend(ax)
f

Of course, in that case we don't really need the AD possibility. But if we reduce the B1+ value the equation becomes complicated enough and might lead to error during derivation if we don't use AD.

deltaB1 = 0.8

amp = MESE_EPG(T2,T1,TE,ETL,deltaB1)
j = ForwardDiff.jacobian(x -> MESE_EPG(x,T1,TE,ETL,deltaB1),[T2])

f = Figure()
ax = Axis(f[1,1], title = "MESE signal with B1 = $(deltaB1)",xlabel="TE [ms]")
lines!(ax,TE_vec,abs.(amp))
ax = Axis(f[1,2], title = "AD of MESE signal with B1 = $(deltaB1)",xlabel="TE [ms]")
lines!(ax,TE_vec,df)
f

Differentiation along multiple variables

If we want to obtain the derivation along T1 and T2 we need to change the EPG_MESE function. The function should take as input a vector containing T2 and T1 (here noted T2/T1) :

function MESE_EPG2(T2T1,TE,ETL,delta)
  T2,T1 = T2T1
  T = complex(eltype(T2))
  E = EPGStates([T(0.0)],[T(0.0)],[T(1.0)])
  echo_vec = Vector{Complex{eltype(T2)}}()

  E = epgRotation(E,pi/2*delta, pi/2)
  ##loop over refocusing-pulses
  for i = 1:ETL
    E = epgDephasing(E,1)
    E = epgRelaxation(E,TE,T1,T2)
    E  = epgRotation(E,pi*delta,0.0)
    E  = epgDephasing(E,1)
    push!(echo_vec,E.Fp[1])
  end

  return abs.(echo_vec)
end

j2 = ForwardDiff.jacobian(x -> MESE_EPG2(x,TE,ETL,deltaB1),[T2,T1])

Here we can see that the second column corresponding to T1 is equal to 0 which is expected for a MESE sequence and the derivative along T2 gives the same results :

j2[:,1] ≈ vec(j)

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils
io = IOBuffer();
versioninfo(io);
split(String(take!(io)), '\n')
16-element Vector{SubString{String}}:
 "Julia Version 1.9.2"
 "Commit e4ee485e909 (2023-07-05 09:39 UTC)"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 8 × Intel(R) Xeon(R) Platinum 8375C CPU @ 2.90GHz"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-14.0.6 (ORCJIT, icelake-server)"
 "  Threads: 1 on 8 virtual cores"
 "Environment:"
 "  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495"
 "  JULIA_PKG_UNPACK_REGISTRY = true"
 "  JULIA_VERSION = 1.9.2"
 "  JULIA_PATH = /usr/local/julia"
 "  JULIA_PKG_SERVER = "
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/.julia/packages/EPGsim/gBRuj/docs/Project.toml`
  [13f3f980] CairoMakie v0.10.8
  [e30172f5] Documenter v0.27.25
  [b6a82cc1] EPGsim v0.1.1 `~/.julia/packages/EPGsim/gBRuj`
  [f6369f11] ForwardDiff v0.10.36
  [98b081ad] Literate v2.14.1
  [b77e0a4c] InteractiveUtils