Blog › ICP guides

Julia developer on retainer: scientific computing, DifferentialEquations.jl, Turing.jl, and GPU acceleration on monthly retainer

September 2, 2026 · ~22 min read

A computational biology group had a simulation performance problem. Their Julia ODE model of a biochemical reaction network took 847 seconds per simulation run on a 24-hour time span — the stiffness of the system (eigenvalue ratio across twelve orders of magnitude) was causing Tsit5, the default explicit solver, to reject 99.3% of proposed integration steps and reduce the step size to picoseconds in order to track the fastest kinetic timescale. The simulation had been running on the explicit solver for eight months because no one had diagnosed why it was slow. A Julia developer on monthly retainer identified the stiffness problem in the first session by examining the Jacobian eigenvalue range using DiffEqSensitivity.jl’s Jacobian utility and confirming that the system was moderately to highly stiff with a condition number in the range that makes explicit methods impractical.

The fix was a solver switch plus an analytical Jacobian. Replacing Tsit5 with Rodas4P — a stiffly stable Rosenbrock method that uses an implicit integration formula — reduced the step rejection rate from 99.3% to 4.1%. The analytical Jacobian, provided via the jac= keyword to ODEFunction, eliminated the finite-difference Jacobian approximation that was spending 63% of the previous solver’s time computing partial derivatives numerically at each step. The tolerances were adjusted to match the physical scale of the species concentrations — reltol=1e-8 and abstol=1e-10 — rather than the default reltol=1e-3 that was insufficient for the accuracy requirements of a pharmacokinetic publication. After these changes, simulation time dropped from 847 seconds to 11 seconds on the same hardware.

No algorithm changed. The ODE system was identical. The reaction kinetics were unchanged. The 77× speedup came from switching to a solver appropriate for the system’s stiffness class and eliminating redundant Jacobian computation. The Julia developer’s retainer invoice logged 16 hours across two sessions: stiffness diagnosis and Jacobian eigenvalue analysis in session one, Rodas4P implementation with analytical Jacobian authoring and tolerance calibration against a known subsystem solution in session two.

The type-stability work that retainers fund

Type stability is the single most important performance concept in Julia, and it is the category of retainer work most frequently invisible in a pull request diff. A type-unstable function — one whose return type cannot be determined at compile time from the argument types alone — forces Julia’s compiler to emit a dynamic dispatch on every call, inhibits LLVM from generating SIMD instructions, and introduces garbage collection pressure when intermediate results cannot be allocated on the stack. The consequence is that a Julia function written by someone accustomed to Python’s dynamic semantics often runs at Python speed, not C speed, even though both the data is statically typed at the call site and the function body contains no inherently dynamic operations.

The diagnostic is @code_warntype f(args...), which prints the compiler’s type inference result for the function body. Any red ::Any annotation indicates a site where the compiler lost type information and will use dynamic dispatch. Common sources include: functions that return different types in different branches (x > 0 ? 1.0 : nothing has return type Union{Float64, Nothing}); global variables that can be reassigned and therefore have unknown type at compile time; container types parameterized on abstract types (Vector{Number} instead of Vector{Float64}) that require dynamic dispatch on every element access; and closures that capture global variables rather than local variables.

A retainer engagement covering type stability typically involves profiling a slow function with BenchmarkTools.jl’s @benchmark f($args) (dollar-sign interpolation captures variables from global scope to avoid the global-variable measurement artifact), identifying functions whose measured throughput is 10× to 100× below the theoretical FLOP-count limit, running @code_warntype on those functions, tracing each ::Any annotation to its source, and restructuring the code to use concrete types throughout. The restructuring typically involves splitting a function into multiple dispatch methods (one per concrete input type), parameterizing structs with type parameters rather than abstract fields (struct Solver{T<:AbstractFloat}; dt::T; end instead of struct Solver; dt::AbstractFloat; end), and replacing global constants that could theoretically be reassigned with const declarations or function arguments. Each fix is a one-to-five line change. Finding the fix requires understanding the Julia compiler’s type inference algorithm, which is a matter of expert knowledge, not code complexity.

DifferentialEquations.jl: solver selection and callbacks

DifferentialEquations.jl’s solver library has over 300 ODE solvers, each appropriate for a different class of problem. The selection decision is invisible to a client reviewing a pull request: the change from solve(prob, Tsit5()) to solve(prob, Rodas4P(), jac=jac_fn) is a one-line diff that produces an 80× speedup on a stiff system and a 0.3× slowdown on a non-stiff system. A retainer Julia developer contributes the expertise to make this selection correctly without trial-and-error across a week of compute time.

The practical selection rules are: for non-stiff ODEs with moderate accuracy requirements (relative tolerance 1e-3 to 1e-6), Tsit5 is the default choice; for high-precision non-stiff problems (relative tolerance below 1e-8), Vern9 or Vern7 provides adaptive step size control with 9th- or 7th-order local error estimation; for stiff systems (eigenvalue ratio above 1000), Rodas4P or QNDF are starting points, with CVODE_BDF from Sundials preferred for very large stiff systems (millions of equations from semi-discretized PDEs) where the sparse Jacobian structure can be exploited; for stiff systems with known Jacobian sparsity patterns, providing the sparsity structure via jac_prototype enables automatic sparse Jacobian coloring that reduces Jacobian computation cost from O(n) to O(k) where k is the number of colors in the graph coloring; and for systems with both stiff and non-stiff components (partitioned ODEs arising in atmospheric chemistry or power grid modeling), implicit-explicit (IMEX) methods like KenCarp4 apply an explicit method to the non-stiff component and an implicit method to the stiff component within each time step.

DifferentialEquations.jl’s callback interface is the second category of solver expertise that retainer engagements typically cover. A ContinuousCallback detects a zero crossing of a condition function during integration and applies a custom effect: stopping at an impact event in a bouncing ball model, resetting a state variable when a membrane potential reaches a threshold in a neuron model, or switching a differential equation’s right-hand side when a control input changes sign. A DiscreteCallback fires at prescribed time points: administering a drug dose at regular intervals in a pharmacokinetic model, recording a state snapshot for a data assimilation checkpoint, or updating a parameter value from an external control signal. Writing these callbacks correctly requires understanding the affect! function signature (affect!(integrator) with access to integrator.u for state, integrator.p for parameters, and integrator.t for current time) and the save_positions argument that controls whether solution points are saved immediately before and after the callback fires. A retainer engagement covering callbacks typically involves translating a physical event description into a correctly functioning callback with appropriate condition and affect functions, verifying that the event detection does not cause step rejection spikes (which it does when the condition function is poorly scaled relative to the ODE’s right-hand side), and testing the callback against a known analytical solution for the simple case before applying it to the full system.

Turing.jl probabilistic programming and posterior inference

Turing.jl enables Bayesian statistical modeling in Julia using a probabilistic programming language embedded via Julia macros. A Turing model is a Julia function annotated with @model where ~ denotes a stochastic relationship: @model function growth_model(y, t) mu = Alpha * exp.(Beta * t); y .~ Normal.(mu, Sigma) end. The ~ operator with a scalar distribution draws from a prior when sampling parameters and evaluates the log density when conditioning on observed data. The .~ broadcasting form applies element-wise for vector observations.

The NUTS sampler (No-U-Turn Sampler) is the default choice for continuous-parameter models because it automatically tunes step size and trajectory length via dual averaging and the U-turn criterion, eliminating the manual hyperparameter tuning that basic HMC requires. Parallelized sampling with sample(model, NUTS(), MCMCThreads(), 2000, 4) runs four chains in parallel across Julia threads, enabling convergence diagnosis via the Gelman-Rubin R-hat statistic (computed with MCMCChains.gelmandiag(chains); values below 1.01 indicate convergence) and effective sample size estimates that account for within-chain autocorrelation. A retainer engagement covering Turing.jl typically involves: designing the model structure to match the generative process of the data rather than the statistical convenience of conjugate priors; identifying parameters whose posterior geometry (funnel shapes from hierarchical models, strong correlations from collinear covariates) causes NUTS divergences, and reparameterizing the model using the non-centered parameterization to improve the posterior geometry; and writing the generated_quantities block for posterior predictive checks by running the model forward with sampled parameter values and comparing the predictive distribution to the observed data distribution.

GPU acceleration with CUDA.jl and kernel design

CUDA.jl exposes NVIDIA GPU computation in Julia through a CuArray type that mirrors Julia’s Array API: CUDA.zeros(Float32, n, m) allocates a GPU array, and most built-in Julia array operations (map, reduce, broadcasting with dot operators, linear algebra via cuBLAS) work on CuArray without code changes. The first GPU acceleration step in a retainer engagement is typically identifying which operations in the client’s compute pipeline are dominated by linear algebra or element-wise transforms (both of which map directly to GPU acceleration) versus which operations have irregular memory access or branching that makes GPU execution inefficient.

Custom GPU kernels with the @cuda macro are appropriate when the required operation is not covered by cuBLAS, cuDNN, or CUDA.jl’s broadcast implementation. A kernel function decorated with function my_kernel(A, B, C) and called with @cuda threads=256 blocks=ceil(Int, n/256) my_kernel(A, B, C) executes on the GPU with each thread processing one element using i = (blockIdx().x - 1) * blockDim().x + threadIdx().x for thread-to-element mapping. The most common performance pitfall is non-coalesced global memory access: a kernel that accesses a 2D array column-by-column in row-major memory layout reads non-contiguous memory addresses across a warp of 32 threads, causing 32 separate cache line fetches instead of 1. Identifying this pattern requires the Nsight Compute profiler (accessible from Julia via NVTX.jl range annotations and CUDA.@profile) and fixing it requires either transposing the memory layout or using shared memory tiling to stage data in fast on-chip memory before the irregular access.

How HourTab tracks Julia developer retainer hours

Julia developer retainers present a specific hour-visibility problem: a nine-hour type stability session produces a five-line diff (three type parameter additions, one dispatch split, one const annotation) and a 12× speedup. The diff is visible; the 847 lines of @code_warntype output that were read to find the five-line fix are not. A time log entry that says “performance optimization, 9h” is technically accurate and completely uninformative to a client who does not read LLVM IR.

HourTab gives Julia consultants a public retainer-hours URL that shows the client the burn-down without a status email. The work log entry is where the diagnostic process is documented: the function name, the @code_warntype annotation observed, the source of the type instability, the dispatch redesign applied, and the @benchmark numbers before and after. That entry takes five minutes to write after the session and converts a nine-hour black box into a nine-hour analysis with documented findings. The next client check-in is a two-sentence summary (“type instability in compute_rates was causing 12× slowdown; fixed with dispatch split, now 0.38s per run”) rather than a twenty-minute tutorial on Julia’s type inference system.

The retainer model fits Julia scientific computing because the performance work is continuous and unpredictable. A Julia upgrade from 1.9 to 1.10 can change the inference behavior for edge-case type patterns, introducing new type instabilities or resolving existing ones. A new package version can change the performance characteristics of a solver or kernel. A new dataset with different stiffness properties can require a solver switch. A project contract closes when the analysis is done. A retainer stays open for the next Julia version, the next stiff ODE, and the next @code_warntype annotation that explains why a simulation that ran fast last quarter is slow this quarter.

FAQ: Julia developer retainers