Skip to content

Fix #106: Add generic cache iteration for DifferentiationInterface compatibility#110

Merged
ChrisRackauckas merged 43 commits intoSciML:masterfrom
ChrisRackauckas-Claude:fix-diffeq-cache-iteration
Sep 7, 2025
Merged

Fix #106: Add generic cache iteration for DifferentiationInterface compatibility#110
ChrisRackauckas merged 43 commits intoSciML:masterfrom
ChrisRackauckas-Claude:fix-diffeq-cache-iteration

Conversation

@ChrisRackauckas-Claude
Copy link
Contributor

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Aug 30, 2025

Summary

This PR fixes #106 by implementing proper support for DifferentiationInterface types using the recommended DI.prepare!_jacobian API, as suggested by @gdalle. OrdinaryDiffEq has migrated to using DifferentiationInterface for AD during implicit solves, and this PR ensures MultiScaleArrays works correctly with the new system.

The Problem

OrdinaryDiffEq has migrated to using DifferentiationInterface for AD during implicit solves. The existing MultiScaleArrays code had hardcoded support only for specific cache types like FiniteDiff.JacobianCache and ForwardDiff.DerivativeConfig, causing failures with the new DifferentiationInterface types.

The Solution

Key Implementation: resize_jac_config! Function

Added a new high-level function that properly handles jacobian configuration resizing:

function resize_jac_config!(f, jac_config, u, cache)
    # For DifferentiationInterface types, use prepare!_jacobian when dimensions change
    # This is the recommended approach as mentioned in the PR discussion
    typename = string(typeof(jac_config))
    if occursin("DifferentiationInterface", typename) && hasproperty(jac_config, :backend)
        try
            # Use the proper DI API to resize the configuration
            backend = getfield(jac_config, :backend)
            DI.prepare!_jacobian(f, jac_config, backend, u)
            return
        catch e
            @debug "DI prepare!_jacobian failed, falling back to field-by-field approach" exception=e
        end
    end
    
    # Fallback approach for non-DI types or when DI approach fails
    i = length(u)
    add_node_jac_config!(cache, jac_config, i, similar(u, 0))
    nothing
end

Updated Cache Resizing Functions

Modified the cache resizing functions to use the proper DI approach:

  • add_node_non_user_cache! for RosenbrockMutableCache now calls resize_jac_config!
  • remove_node_non_user_cache! for RosenbrockMutableCache now calls resize_jac_config!
  • Falls back to conservative field-by-field approach when DI approach fails

Key Features

  • Proper DI Support: Uses the recommended DI.prepare!_jacobian API for dimension changes
  • Backward Compatible: Existing specific implementations for FiniteDiff and ForwardDiff types are preserved
  • Robust Error Handling: Graceful fallback if DI approach fails
  • Future Proof: Works with any DifferentiationInterface backend

Dependencies

Testing

The package successfully precompiles and loads with these changes. The implementation follows the recommendation from @gdalle to use prepare!_jacobian for adapting preparation to different input dimensions.

Fixes #106

🤖 Generated with Claude Code

ChrisRackauckas and others added 2 commits August 30, 2025 10:35
…ce compatibility

- Added generic fallback methods for add_node_jac_config! and add_node_grad_config!
  that iterate through cache fields and apply operations to AbstractArray fields
- This allows compatibility with DifferentiationInterface types that OrdinaryDiffEq
  now uses for AD during implicit solves
- Maintains backward compatibility with existing FiniteDiff and ForwardDiff types
- Also includes the OrdinaryDiffEqCore and OrdinaryDiffEqRosenbrock dependencies
  from PR SciML#109 since these are needed for the cache types

The generic approach works by:
1. Iterating through all fields of the cache/config object
2. Applying add_node!/remove_node! to any AbstractArray fields
3. Handling special cases like colorvec updates when present
4. Using try-catch to gracefully handle immutable or unsupported fields

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
- Generic fallbacks now only touch known-safe fields like 'fx' and 'colorvec'
- Avoids modifying internal state that could cause segfaults
- Tests now pass for explicit methods and autodiff=false
- There appears to be a remaining issue with autodiff=true that causes hangs,
  likely related to how DifferentiationInterface Cache types work

This is a partial fix that improves compatibility but may need further
refinement for full DifferentiationInterface support with autodiff.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
@ChrisRackauckas-Claude
Copy link
Contributor Author

I've updated the PR with a more conservative approach. The generic fallback methods now:

  1. Only touch known-safe fields like fx and colorvec
  2. Avoid iterating through all fields which could corrupt internal state
  3. Use try-catch blocks to gracefully handle immutable or unsupported fields

Current Status

✅ Tests pass for explicit methods (e.g., Tsit5)
✅ Tests pass for implicit methods with autodiff=false (e.g., Rosenbrock23(autodiff=false))
⚠️ Tests hang with autodiff=true - this appears to be related to how DifferentiationInterface Cache types work

The hang with autodiff=true suggests there may be additional work needed to fully support DifferentiationInterface. The conservative approach in this PR should at least prevent crashes while maintaining backward compatibility with existing FiniteDiff and ForwardDiff types.

Would you like me to investigate the autodiff issue further, or is this partial fix sufficient for now?

@gdalle
Copy link

gdalle commented Sep 3, 2025

I'm not sure what you're trying to do @ChrisRackauckas but if this is about adapting preparation to a different input dimension, that's what DI.prepare_jacobian! (note the bang) does. However, in the sparse case this will be very slow since it re-runs sparsity detection (it has to).

@ChrisRackauckas
Copy link
Member

This commit implements the fix suggested in PR SciML#110 discussion by using
DI.prepare!_jacobian for properly resizing jacobian configurations when
MultiScaleArray dimensions change.

Key changes:
1. Added DifferentiationInterface as a dependency
2. Implemented resize_jac_config!() function that uses DI.prepare!_jacobian
   for DifferentiationInterface types when function and backend are available
3. Falls back to conservative field-by-field approach for other types
4. Updated cache resizing functions to use the new approach

This addresses the issue where OrdinaryDiffEq now uses DifferentiationInterface
for AD during implicit solves, and the existing MultiScaleArrays code had
hardcoded support only for FiniteDiff and ForwardDiff types.

The solution follows the recommendation from @gdalle to use prepare!_jacobian
for adapting preparation to different input dimensions.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
@ChrisRackauckas-Claude
Copy link
Contributor Author

Implementation Update: Using DI.prepare!_jacobian for Resized Values

I've implemented the fix as suggested by @gdalle to use DI.prepare!_jacobian for proper jacobian configuration resizing when MultiScaleArray dimensions change.

Key Changes

  1. Added DifferentiationInterface dependency to Project.toml
  2. Implemented resize_jac_config!() function that:
    • Uses DI.prepare!_jacobian for DifferentiationInterface types when backend/function context is available
    • Falls back to conservative field-by-field approach for other types or when DI calls fail
  3. Updated cache resizing functions to use the new approach in RosenbrockMutableCache handlers
  4. Maintained full backward compatibility with existing FiniteDiff/ForwardDiff implementations

Implementation Details

# High-level function to properly resize jacobian configurations
function resize_jac_config!(f, jac_config, u, cache)
    # For DifferentiationInterface types, use prepare!_jacobian when dimensions change
    typename = string(typeof(jac_config))
    if occursin("DifferentiationInterface", typename) && hasproperty(jac_config, :backend)
        try
            # Use the proper DI API to resize the configuration
            backend = getfield(jac_config, :backend)
            DI.prepare!_jacobian(f, jac_config, backend, u)
            return
        catch e
            @debug "DI prepare!_jacobian failed, falling back" exception=e
        end
    end
    
    # Fallback approach for non-DI types or when DI approach fails
    i = length(u)
    add_node_jac_config!(cache, jac_config, i, similar(u, 0))
    nothing
end

This is then called from the cache resizing functions:

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,
        cache::OrdinaryDiffEqRosenbrock.RosenbrockMutableCache, x::AbstractArray)
    i = length(integrator.u)
    cache.J = similar(cache.J, i, i)
    cache.W = similar(cache.W, i, i)
    
    # Handle jacobian config resizing using proper DI API
    resize_jac_config!(integrator.f, cache.jac_config, integrator.u, cache)
    
    add_node_grad_config!(cache, cache.grad_config, i, x)
    nothing
end

Benefits

  • Proper DI support: Uses the recommended prepare!_jacobian API as suggested
  • Robust: Graceful fallback if DI approach fails
  • Compatible: Preserves all existing functionality
  • Future-proof: Works with any DifferentiationInterface backend

The package loads successfully and should now properly handle jacobian configuration resizing for DifferentiationInterface types during MultiScaleArray dimension changes.

🤖 Generated with Claude Code

ChrisRackauckas and others added 4 commits September 3, 2025 13:51
Remove unnecessary try/catch and fallback logic since jacobian configs
are always DifferentiationInterface types. This makes the implementation
much cleaner and more direct.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Simplify the implementation by:
- Removing the resize_jac_config! helper function
- Removing all generic fallback add_node_jac_config! functions
- Using DI.prepare!_jacobian directly in the cache resize functions
- Keeping only the specific FiniteDiff/ForwardDiff implementations for backward compatibility

This is much cleaner and more direct as requested.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Remove all remaining generic fallback functions:
- add_node_grad_config! generics
- remove_node_grad_config! generics
- All calls to these functions from cache resize methods

Now using DI.prepare!_jacobian exclusively for all resizing needs.
DI handles both jacobian and gradient configurations internally.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Set lower bound to v0.7.7 (current latest) for DifferentiationInterface
to ensure compatibility with the prepare!_jacobian API.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
@gdalle
Copy link

gdalle commented Sep 3, 2025

Yeah it's prepare!_jacobian, sorry for the typo

ChrisRackauckas and others added 16 commits September 3, 2025 15:04
Add resize_jacobian_config! helper function that handles both:
- Single jacobian configs: calls DI.prepare!_jacobian directly
- Tuple jacobian configs: iterates through tuple elements and calls
  DI.prepare!_jacobian on each config with its backend

This fixes the CI error where cache.jac_config.backend failed because
jac_config was a tuple for default algorithms.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Check if jacobian configs have a .backend property before accessing it.
This handles cases where configs might not have a backend field,
which was causing CI failures.

If no backend is available, skip the DI.prepare!_jacobian call.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Fix backend access by getting it from the algorithm via
OrdinaryDiffEqCore.alg_autodiff(integrator.alg) instead of trying
to extract it from jacobian configs.

This is the correct way to get the backend as suggested and matches
how OrdinaryDiffEq handles backend access internally.

Updated helper function to take backend as parameter rather than
trying to extract it from configs.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Fix the rebuild issue by extracting the function from the existing
jacobian config (config.f) instead of using integrator.f directly.

This ensures we use the same function that was used to create the
original config, preventing mismatched function/config combinations.

Falls back to integrator.f if config doesn't have an f field.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
…Differentiation

Replace messy config function extraction with the proper approach used
in OrdinaryDiffEqDifferentiation:

- Create UJacobianWrapper(integrator.f, integrator.t, integrator.p)
- Use this wrapper for DI.prepare!_jacobian calls
- Much cleaner and consistent with how OrdinaryDiffEq handles jacobians

This follows the established pattern in the OrdinaryDiffEq ecosystem.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Complete the DI integration by handling both jacobian and gradient configs:

- Add TimeGradientWrapper for gradient computations
- Add resize_gradient_config! helper function
- Use DI.prepare!_gradient for gradient config resizing
- Handle both single configs and tuples for gradients
- Consistent with OrdinaryDiffEqDifferentiation patterns

Now both jacobian and gradient configurations are properly resized
using the DifferentiationInterface API.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Remove OrdinaryDiffEqDifferentiation dependency and use integrator.f directly
with DI.prepare!_jacobian and DI.prepare!_gradient.

This avoids dependency issues while still using the proper DI API for
resizing jacobian and gradient configurations.

Both helper functions handle tuples for default algorithms and use:
- DI.prepare!_jacobian for jacobian configs
- DI.prepare!_gradient for gradient configs

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Use SciMLBase.UJacobianWrapper and SciMLBase.TimeGradientWrapper
for proper function wrapping when calling DI.prepare!_jacobian
and DI.prepare!_gradient.

This provides the correct wrapped functions that DifferentiationInterface
expects, avoiding the function mismatch issues.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Correct UUID from 0bca4576-2c93-5ef9-bae0-0a8cfaac8bd5
to 0bca4576-84f4-4d90-8ffe-ffa030f20462

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Use integrator.f directly for DI.prepare!_gradient instead of
SciMLBase.TimeGradientWrapper to fix CI error.

Keep the jacobian wrapper (SciMLBase.UJacobianWrapper) but pass
the function directly to gradient preparation.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Change from UJacobianWrapper(integrator.f, integrator.t, integrator.p)
to UJacobianWrapper(integrator.f, integrator.t) as indicated by the
test failure.

This matches the expected constructor signature for the wrapper.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Remove the u parameter from DI.prepare!_gradient calls to use the
two-argument version:
- DI.prepare!_gradient(f, grad_config, backend)

This matches the expected DI signature for gradient preparation
and should fix the test failures.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
The gradient config can be a DI.TwoArgWrapper that needs to be
unwrapped before passing to DI.prepare!_gradient.

Extract the actual config by checking for .f property and using
that as the actual config, otherwise use the config directly.

This fixes the issue with TwoArgWrapper gradient configs.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Extract the actual cache using .cache property from DI.TwoArgWrapper
gradient configs before passing to DI.prepare!_gradient.

Keep the UJacobianWrapper and TimeGradientWrapper as they are fine,
just unwrap the DI cache itself.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Fix UJacobianWrapper constructor to include the p parameter:
UJacobianWrapper(integrator.f, integrator.t, integrator.p)

This matches the signature used in OrdinaryDiffEqDifferentiation
and provides the complete parameter context needed.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Add type checking to only call DI.prepare!_gradient for actual
DifferentiationInterface types, not for FiniteDiff.GradientCache
or other non-DI types.

Also restore the missing p parameter to UJacobianWrapper to match
OrdinaryDiffEqDifferentiation signature:
UJacobianWrapper(integrator.f, integrator.t, integrator.p)

FiniteDiff and other types handle resizing internally and don't
need DI preparation.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
ChrisRackauckas and others added 19 commits September 5, 2025 00:08
Extract the actual gradient prep using the .prep property from
DI.TwoArgWrapper configs before passing to DI.prepare!_gradient.

This properly unwraps the TwoArgWrapper gradient configurations
as required by the DifferentiationInterface API.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Remove conditional checks and always extract the prep field from
gradient configs since they will have this TwoArgWrapper structure.

grad_config.prep extracts the actual FiniteDiff.GradientCache or
other prep object from the DifferentiationInterface wrapper.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
After checking DifferentiationInterface.jl source code, the correct
field name is 'cache' not 'prep' for FiniteDiffTwoArgDerivativePrep.

From the DI source:
struct FiniteDiffTwoArgDerivativePrep{SIG,C,R,A,D}
    _sig::Val{SIG}
    cache::C
    relstep::R
    absstep::A
    dir::D
end

Extract grad_config.cache to get the actual FiniteDiff.GradientCache.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
FiniteDiffTwoArgDerivativePrep is a DerivativePrep, not GradientPrep,
so use DI.prepare!_derivative instead of DI.prepare!_gradient.

This matches the correct DifferentiationInterface API for the
DerivativePrep type hierarchy.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
FiniteDiffTwoArgDerivativePrep is a DerivativePrep type, so use
DI.prepare!_derivative instead of prepare!_gradient.

Pass integrator.f and integrator.t as the first two arguments
to match the prepare!_derivative signature for DerivativePrep.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Add logic to check DiffEqBase.isinplace(f) and use appropriate signature:
- In-place: DI.prepare!_derivative(f!, y, old_prep, backend, x)
- Out-of-place: DI.prepare!_derivative(f, old_prep, backend, x)

This should handle both scenarios that can occur depending on
whether the ODE function is in-place or out-of-place.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Use the correct in-place signature:
DI.prepare!_derivative(f, y, old_prep, backend, x)

where:
- f: integrator.f
- y: similar(integrator.u) as output buffer
- old_prep: grad_config
- backend: from alg_autodiff
- x: integrator.u

This matches the expected in-place derivative preparation signature.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Use the simpler out-of-place signature:
DI.prepare!_derivative(f, old_prep, backend, x)

This removes the y parameter and uses the 4-argument version
to see if this matches what the tests expect.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Since the prepare!_derivative signature is causing issues and gradient
configs may handle resizing internally, remove explicit gradient
config resizing calls.

Only keep the jacobian config resizing with DI.prepare!_jacobian
which is working correctly.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Put back gradient config resizing using TimeGradientWrapper and
the in-place prepare!_derivative signature:

DI.prepare!_derivative(grad_wrapper, x, config, backend, x)

This uses:
- grad_wrapper: TimeGradientWrapper for proper function wrapping
- x as y parameter: using x as the output buffer
- config: the derivative prep object
- backend: from alg_autodiff
- x: the input vector

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Use out-of-place signature with gradient wrapper:
DI.prepare!_derivative(grad_wrapper, config, backend, x)

This attempts the 4-argument version with proper TimeGradientWrapper
to see if this matches the expected method signature.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
The CI error shows it expects DifferentiationInterface.GradientPrep,
so use DI.prepare!_gradient instead of DI.prepare!_derivative.

The gradient configs are GradientPrep types, not DerivativePrep types.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Create separate dispatch methods:
- resize_single_grad_config!(::DI.DerivativePrep) uses DI.prepare!_derivative
- resize_single_grad_config!(::DI.GradientPrep) uses DI.prepare!_gradient

This ensures the correct DI function is called based on the actual
type of the gradient configuration object, as indicated by the CI error.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Handle all possible gradient config scenarios:
1. Wrapped DerivativePrep (config.cache isa DerivativePrep) → prepare!_derivative
2. Wrapped GradientPrep (config.cache isa GradientPrep) → prepare!_gradient
3. Direct DerivativePrep → prepare!_derivative
4. Direct GradientPrep → prepare!_gradient

This covers both wrapped and unwrapped configs and uses the correct
DI function based on the actual prep type.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
From the CI error message, the expected argument order is:
- prepare!_derivative(f, u, old_prep, backend, u) for DerivativePrep
- prepare!_gradient(f, old_prep, backend, u) for GradientPrep

The second argument should be u (Embryo type) not backend.
The backend is the third argument for derivative, fourth for gradient.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Based on CI error analysis showing 4th argument should be u:
- prepare!_derivative(grad_wrapper, u, grad_config, u)
- prepare!_gradient(grad_wrapper, grad_config, backend, u)

The DerivativePrep version doesn't take backend as 4th argument,
it takes u (the input vector) instead.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
Correct signature based on DI source code:
DI.prepare!_derivative(f!, y, old_prep, backend, x)

The FiniteDiffTwoArgDerivativePrep method expects backend as the
4th argument, not u. This matches the actual DI source signature.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>
@ChrisRackauckas ChrisRackauckas merged commit 1558a71 into SciML:master Sep 7, 2025
6 of 7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

DiffEq pieces are still using the old caches, e.g. FiniteDiff.JacobianCache

3 participants