Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correctly initializes a SparseConnectivityTracer.GradientTracer cache #117

Closed
franckgaga opened this issue Mar 13, 2025 · 7 comments
Closed
Labels
bug Something isn't working

Comments

@franckgaga
Copy link

franckgaga commented Mar 13, 2025

Describe the bug 🐞

Correctly initializes a SparseConnectivityTracer.GradientTracer cache at the value specified at construction.

Expected behavior

When differentiating sparse Jacobians with DifferentiationInterface.jl and SparseConnectivityTracer, the floats and duals in a DiffCache are correctly intializes with the specified value at construction, but not the SparseConnectivityTracer.GradientTracer (they start undefined).

Minimal Reproducible Example 👇

using DifferentiationInterface
using SparseConnectivityTracer
using SparseMatrixColorings
using PreallocationTools
x = rand(10)
y = Vector{eltype(x)}(undef, 10)
cache = DiffCache(zeros(eltype(x), length(y))) # should be 0s for all floats, duals, etc.
function f!(y, x::Vector{T}) where T
    y .= get_tmp(cache, T)
    println(y)
    return y
end
sparse_backend = AutoSparse(
    AutoForwardDiff();
    sparsity_detector=TracerSparsityDetector(),
    coloring_algorithm=GreedyColoringAlgorithm(),
)
sparse_prep = prepare_jacobian(f!, y, sparse_backend, x)

Error & Stacktrace ⚠️

SparseConnectivityTracer.GradientTracer{SparseConnectivityTracer.IndexSetGradientPattern{Int64, BitSet}}[#undef, #undef, #undef, #undef, #undef, #undef, #undef, #undef, #undef]
ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ ./essentials.jl:917 [inlined]
  [2] iterate (repeats 2 times)
    @ ./array.jl:902 [inlined]
  [3] iterate
    @ ./iterators.jl:206 [inlined]
  [4] iterate
    @ ./iterators.jl:205 [inlined]
  [5] jacobian_tracers_to_matrix(xt::Vector{…}, yt::Vector{…})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/Cscac/src/parse_outputs_to_matrix.jl:17
  [6] _jacobian_sparsity(f!::Function, y::Vector{…}, x::Vector{…}, ::Type{…})
    @ SparseConnectivityTracer ~/.julia/packages/SparseConnectivityTracer/Cscac/src/trace_functions.jl:88
  [7] jacobian_sparsity
    @ ~/.julia/packages/SparseConnectivityTracer/Cscac/src/adtypes_interface.jl:60 [inlined]
  [8] _prepare_sparse_jacobian_aux(::DifferentiationInterface.PushforwardFast, ::Vector{…}, ::Tuple{…}, ::AutoSparse{…}, ::Vector{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/qrWdQ/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:61
  [9] prepare_jacobian(::typeof(f!), ::Vector{…}, ::AutoSparse{…}, ::Vector{…})
    @ DifferentiationInterfaceSparseMatrixColoringsExt ~/.julia/packages/DifferentiationInterface/qrWdQ/ext/DifferentiationInterfaceSparseMatrixColoringsExt/jacobian.jl:49
 [10] top-level scope
    @ ~/Dropbox/Programmation/Julia/TestMPC/src/di.jl:91
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

  • Output of using Pkg; Pkg.status()
Project TestMPC v0.1.0
Status `~/Dropbox/Programmation/Julia/TestMPC/Project.toml`
  [6e4b80f9] BenchmarkTools v1.6.0
  [aaaaaaaa] ControlSystemsBase v1.14.5
  [c47d62df] DAQP v0.6.0
  [163ba53b] DiffResults v1.1.0
  [a0c0ee7d] DifferentiationInterface v0.6.43
  [e30172f5] Documenter v1.8.1
  [f6369f11] ForwardDiff v0.10.38
  [87dc4568] HiGHS v1.14.0
  [b6b21f68] Ipopt v1.7.3
  [c3a54625] JET v0.9.18
  [4076af6c] JuMP v1.24.0
  [67920dd8] KNITRO v0.14.4
  [33e6dc65] MKL v0.8.0
  [2621e9c9] MadNLP v0.8.5
  [61f9bdb8] ModelPredictiveControl v1.4.2 `~/.julia/dev/ModelPredictiveControl`
  [ab2f91bb] OSQP v0.8.1
  [ccf2f8ad] PlotThemes v3.3.0
  [91a5bcdd] Plots v1.40.10
  [d236fae5] PreallocationTools v0.4.25
  [21216c6a] Preferences v1.4.3
  [295af30f] Revise v3.7.2
  [9f842d2f] SparseConnectivityTracer v0.6.14
  [0a514795] SparseMatrixColorings v0.4.14
  [f8b46487] TestItemRunner v1.1.0
  [37e2e46d] LinearAlgebra v1.11.0
  • Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Project TestMPC v0.1.0
Status `~/Dropbox/Programmation/Julia/TestMPC/Manifest.toml`
  [47edcb42] ADTypes v1.14.0
  [14f7f29c] AMD v0.5.3
  [a4c015fc] ANSIColoredPrinters v0.0.1
  [1520ce14] AbstractTrees v0.4.5
  [79e6a3ab] Adapt v4.2.0
  [66dad0bd] AliasTables v1.1.3
  [4fba245c] ArrayInterface v7.18.0
  [6e4b80f9] BenchmarkTools v1.6.0
  [d1d4a3ce] BitFlags v0.1.9
  [da1fd8a2] CodeTracking v1.3.6
  [523fee87] CodecBzip2 v0.8.5
  [944b1d66] CodecZlib v0.7.8
  [35d6a980] ColorSchemes v3.29.0
  [3da002f7] ColorTypes v0.12.0
  [c3611d14] ColorVectorSpace v0.11.0
  [5ae59095] Colors v0.13.0
  [bbf7d656] CommonSubexpressions v0.3.1
  [34da2185] Compat v4.16.0
  [f0e56b4a] ConcurrentUtilities v2.5.0
  [187b0558] ConstructionBase v1.5.8
  [d38c429a] Contour v0.6.3
  [aaaaaaaa] ControlSystemsBase v1.14.5
  [c47d62df] DAQP v0.6.0
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.20
  [8bb1440f] DelimitedFiles v1.9.1
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.6.43
  [ffbed154] DocStringExtensions v0.9.3
  [e30172f5] Documenter v1.8.1
  [460bff9d] ExceptionUnwrapping v0.1.11
  [e2ba6199] ExprTools v0.1.10
  [c87230d0] FFMPEG v0.4.2
  [9aa1b823] FastClosures v0.3.2
  [1a297f60] FillArrays v1.13.0
  [53c48c17] FixedPointNumbers v0.8.5
  [1fa38f19] Format v1.3.7
  [f6369f11] ForwardDiff v0.10.38
  [28b8d3ca] GR v0.73.13
  [d7ba0133] Git v1.3.1
  [42e2da0e] Grisu v1.0.2
  [cd3eb016] HTTP v1.10.15
  [87dc4568] HiGHS v1.14.0
  [b5f81e59] IOCapture v0.2.5
  [b6b21f68] Ipopt v1.7.3
  [92d709cd] IrrationalConstants v0.2.4
  [c3a54625] JET v0.9.18
  [1019f520] JLFzf v0.1.9
  [692b3bcd] JLLWrappers v1.7.0
  [682c06a0] JSON v0.21.4
  [0f8b85d8] JSON3 v1.14.1
  [4076af6c] JuMP v1.24.0
  [aa1ae85d] JuliaInterpreter v0.9.42
  [70703baa] JuliaSyntax v0.4.10
  [67920dd8] KNITRO v0.14.4
  [40e66cde] LDLFactorizations v0.10.1
  [b964fa9f] LaTeXStrings v1.4.0
  [23fbe1c1] Latexify v0.16.6
  [0e77f7df] LazilyInitializedFields v1.3.0
  [7a12625a] LinearMaps v3.11.4
  [5c8ed15e] LinearOperators v2.9.0
  [2ab3a3ac] LogExpFunctions v0.3.29
  [e6f89c97] LoggingExtras v1.1.0
  [6f1432cf] LoweredCodeUtils v3.1.0
  [33e6dc65] MKL v0.8.0
  [1914dd2f] MacroTools v0.5.15
  [2621e9c9] MadNLP v0.8.5
  [d0879d2d] MarkdownAST v0.1.2
  [b8f27783] MathOptInterface v1.37.2
  [99c1a7ee] MatrixEquations v2.4.2
  [48965c70] MatrixPencils v1.8.1
  [739be429] MbedTLS v1.1.9
  [442fdcdd] Measures v0.3.2
  [e1d29d7a] Missings v1.2.0
  [61f9bdb8] ModelPredictiveControl v1.4.2 `~/.julia/dev/ModelPredictiveControl`
  [d8a4904e] MutableArithmetics v1.6.4
  [a4795742] NLPModels v0.21.3
  [77ba4419] NaNMath v1.1.2
  [ab2f91bb] OSQP v0.8.1
  [4d8831e6] OpenSSL v1.4.3
  [bac558e1] OrderedCollections v1.8.0
  [69de0a69] Parsers v2.8.1
  [b98c9c47] Pipe v1.3.0
  [ccf2f8ad] PlotThemes v3.3.0
  [995b91a9] PlotUtils v1.4.3
  [91a5bcdd] Plots v1.40.10
  [f27b6e38] Polynomials v4.0.19
  [d236fae5] PreallocationTools v0.4.25
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [33c8b6b6] ProgressLogging v0.1.4
  [43287f4e] PtrArrays v1.3.0
  [3cdcf5f2] RecipesBase v1.3.4
  [01d81517] RecipesPipeline v0.6.12
  [189a3867] Reexport v1.2.2
  [2792f1a3] RegistryInstances v0.1.0
  [05181044] RelocatableFolders v1.0.1
  [ae029012] Requires v1.3.1
  [295af30f] Revise v3.7.2
  [6c6a2e73] Scratch v1.2.1
  [efcf1570] Setfield v1.1.2
  [992d4aef] Showoff v1.0.3
  [777ac1f9] SimpleBufferStream v1.2.0
  [ff4d7338] SolverCore v0.3.8
  [a2af1166] SortingAlgorithms v1.2.1
  [9f842d2f] SparseConnectivityTracer v0.6.14
  [0a514795] SparseMatrixColorings v0.4.14
  [276daf66] SpecialFunctions v2.5.0
  [860ef19b] StableRNGs v1.0.2
  [1e83bf80] StaticArraysCore v1.4.3
  [10745b16] Statistics v1.11.1
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.4
  [856f2bd8] StructTypes v1.11.0
  [62fd8b95] TensorCore v0.1.1
  [f8b46487] TestItemRunner v1.1.0
  [1c621080] TestItems v1.0.0
  [a759f4b9] TimerOutputs v0.5.28
  [3bb67fe8] TranscodingStreams v0.11.3
  [5c2747f8] URIs v1.5.1
  [1cfade01] UnicodeFun v0.4.1
  [1986cc42] Unitful v1.22.0
  [45397f5d] UnitfulLatexify v1.6.4
  [41fe7b60] Unzip v0.2.0
  [ae81ac8f] ASL_jll v0.1.3+0
  [6e34b625] Bzip2_jll v1.0.9+0
  [83423d85] Cairo_jll v1.18.2+1
  [5c51c916] DAQP_jll v0.6.0+0
  [ee1fde0b] Dbus_jll v1.14.10+0
  [2702e6a9] EpollShim_jll v0.0.20230411+1
  [2e619515] Expat_jll v2.6.5+0
⌅ [b22a6f82] FFMPEG_jll v4.4.4+1
  [a3f928ae] Fontconfig_jll v2.15.0+0
  [d7e528f0] FreeType2_jll v2.13.3+1
  [559328eb] FriBidi_jll v1.0.16+0
  [0656b61e] GLFW_jll v3.4.0+2
  [d2c73de3] GR_jll v0.73.13+0
  [78b55507] Gettext_jll v0.21.0+0
  [f8c6e375] Git_jll v2.47.1+0
  [7746bdde] Glib_jll v2.82.4+0
  [3b182d85] Graphite2_jll v1.3.14+1
  [2e76f6c2] HarfBuzz_jll v8.5.0+0
  [8fd58aa0] HiGHS_jll v1.9.0+0
  [e33a78d0] Hwloc_jll v2.12.0+0
  [1d5cc7b8] IntelOpenMP_jll v2025.0.4+0
  [9cc047cb] Ipopt_jll v300.1400.1700+0
  [aacddb02] JpegTurbo_jll v3.1.1+0
  [c1c5ebd0] LAME_jll v3.100.2+0
  [88015f11] LERC_jll v4.0.1+0
  [1d63c593] LLVMOpenMP_jll v18.1.7+0
  [dd4b983a] LZO_jll v2.10.3+0
⌅ [e9f186c6] Libffi_jll v3.2.2+2
  [d4300ac3] Libgcrypt_jll v1.11.0+0
  [7e76a0d4] Libglvnd_jll v1.7.0+0
  [7add5ba3] Libgpg_error_jll v1.51.1+0
  [94ce4f54] Libiconv_jll v1.18.0+0
  [4b2f31a3] Libmount_jll v2.40.3+0
  [89763e89] Libtiff_jll v4.7.1+0
  [38a345b3] Libuuid_jll v2.40.3+0
  [d00139f3] METIS_jll v5.1.3+0
  [856f044c] MKL_jll v2025.0.1+1
  [d7ed1dd3] MUMPS_seq_jll v500.700.301+0
  [9c4f68bf] OSQP_jll v0.600.200+0
  [e7412a2a] Ogg_jll v1.3.5+1
  [656ef2d0] OpenBLAS32_jll v0.3.29+0
  [458c3c95] OpenSSL_jll v3.0.16+0
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [91d4177d] Opus_jll v1.3.3+0
  [36c8627f] Pango_jll v1.56.1+0
⌅ [30392449] Pixman_jll v0.43.4+0
⌅ [c0090381] Qt6Base_jll v6.7.1+1
⌅ [629bc702] Qt6Declarative_jll v6.7.1+2
⌅ [ce943373] Qt6ShaderTools_jll v6.7.1+1
⌃ [e99dba38] Qt6Wayland_jll v6.7.1+1
⌅ [319450e9] SPRAL_jll v2024.5.8+0
  [a44049a8] Vulkan_Loader_jll v1.3.243+0
  [a2964d1f] Wayland_jll v1.21.0+2
  [2381bf8a] Wayland_protocols_jll v1.36.0+0
  [02c8fc9c] XML2_jll v2.13.6+1
  [aed1982a] XSLT_jll v1.1.42+0
  [ffd25f8a] XZ_jll v5.6.4+1
  [f67eecfb] Xorg_libICE_jll v1.1.1+0
  [c834827a] Xorg_libSM_jll v1.2.4+0
  [4f6342f7] Xorg_libX11_jll v1.8.6+3
  [0c0b7dd1] Xorg_libXau_jll v1.0.12+0
  [935fb764] Xorg_libXcursor_jll v1.2.3+0
  [a3789734] Xorg_libXdmcp_jll v1.1.5+0
  [1082639a] Xorg_libXext_jll v1.3.6+3
  [d091e8ba] Xorg_libXfixes_jll v6.0.0+0
  [a51aa0fd] Xorg_libXi_jll v1.8.2+0
  [d1454406] Xorg_libXinerama_jll v1.1.5+0
  [ec84b674] Xorg_libXrandr_jll v1.5.4+0
  [ea2f1a96] Xorg_libXrender_jll v0.9.11+1
  [14d82f49] Xorg_libpthread_stubs_jll v0.1.2+0
  [c7cfdc94] Xorg_libxcb_jll v1.17.0+3
  [cc61e674] Xorg_libxkbfile_jll v1.1.2+1
  [e920d4aa] Xorg_xcb_util_cursor_jll v0.1.4+0
  [12413925] Xorg_xcb_util_image_jll v0.4.0+1
  [2def613f] Xorg_xcb_util_jll v0.4.0+1
  [975044d2] Xorg_xcb_util_keysyms_jll v0.4.0+1
  [0d47668e] Xorg_xcb_util_renderutil_jll v0.3.9+1
  [c22f9ab0] Xorg_xcb_util_wm_jll v0.4.1+1
  [35661453] Xorg_xkbcomp_jll v1.4.6+1
  [33bec58e] Xorg_xkeyboard_config_jll v2.39.0+0
  [c5fb5394] Xorg_xtrans_jll v1.5.1+0
  [3161d3a3] Zstd_jll v1.5.7+1
  [35ca27e7] eudev_jll v3.2.9+0
  [214eeab7] fzf_jll v0.56.3+0
  [1a1c6b14] gperf_jll v3.1.1+1
  [a4ae2306] libaom_jll v3.11.0+0
  [0ac62f75] libass_jll v0.15.2+0
  [1183f4f0] libdecor_jll v0.2.2+0
  [2db6ffa8] libevdev_jll v1.11.0+0
  [f638f0a6] libfdk_aac_jll v2.0.3+0
  [36db933b] libinput_jll v1.18.0+0
  [b53b4c65] libpng_jll v1.6.47+0
  [f27f6e37] libvorbis_jll v1.3.7+2
  [009596ad] mtdev_jll v1.1.6+0
  [1317d2d5] oneTBB_jll v2022.0.0+0
⌅ [1270edf5] x264_jll v2021.5.5+0
⌅ [dfaa095f] x265_jll v3.5.0+0
  [d8fb68d0] xkbcommon_jll v1.4.1+2
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [9abbd945] Profile v1.11.0
  [3fa0cd96] REPL v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [f489334b] StyledStrings v1.11.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.1+4
  [efcefdf7] PCRE2_jll v10.42.0+1
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
  • Output of versioninfo()
Julia Version 1.11.4
Commit 8561cc3d68d (2025-03-10 11:36 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × 12th Gen Intel(R) Core(TM) i7-1270P
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 16 default, 0 interactive, 8 GC (on 16 virtual cores)
Environment:
  LD_LIBRARY_PATH = /home/francis/Documents/Applications/knitro-14.0.0-Linux64/lib:/home/francis/Documents/Applications/knitro-14.0.0-Linux64/lib:
  JULIA_EDITOR = code

Additional context

It's for the integration of DifferentiationInterface.jl in `ModelPredictiveControl.jl.

@franckgaga
Copy link
Author

franckgaga commented Mar 13, 2025

Also what AD backend DiffCache supports ? I did not found this information in the documentation here. Does all AbstractADType will work? (e.g. AutoEnzyme and AutoFiniteDiff)

edit:
just did some simple tests and the following backend works with DiffCache:

  • AutoFiniteDiff
  • AutoForwardDiff
  • AutoReverseDiff

but not them:

  • AutoZygote
  • AutoEnzyme

@thomvet
Copy link
Contributor

thomvet commented Mar 13, 2025

It's meant to be used with ForwardDiff; it essentially provides caches that swap between Dual-numbers (that carry the value and partial values with them) and "normal" numbers. That was at least the initial intent :)

@franckgaga
Copy link
Author

franckgaga commented Mar 13, 2025

@ChrisRackauckas mentioned on discourse that it should works with any differentiation tools that is based on type. My understanding is that it include all AbstractADType except AutoZygote and AutoEnzyme (and also AutoMooncake?)

Knowing that, I'm wondering if there is a way to avoid PreallocationTools in MPC.jl (while keeping the allocations as small as possible, and in a function that will be differentiated. There always the option to allocate new vectors at each call according to the type of x elements, but allocations impact the performance!)

@ChrisRackauckas
Copy link
Member

Well the issue here is a user issue. The function being tested doesn't make sense.

function f!(y, x::Vector{T}) where T
    y .= get_tmp(cache, T)
    println(y)
    return y
end

This violates the whole purpose of a cache, since there isn't anything being cached. Because of this, you get an uninitialized memory error because nothing has been written to the memory that you're using. You might've intended:

function f!(y, x::Vector{T}) where T
    tmp = get_tmp(cache, T)
    tmp .= 0
    y .= tmp
    println(y)
    return y
end

or

function f!(y, x::Vector{T}) where T
    tmp = get_tmp(cache, T)
    tmp .= x
    y .= tmp
    println(y)
    return y
end

either of which are fine.

@franckgaga
Copy link
Author

franckgaga commented Mar 21, 2025

@ChrisRackauckas

Well the issue here is a user issue. The function being tested doesn't make sense.

That's not very a professional phrasing. Maybe you should read again your SciML Community Standards, but let's put this aside.

That was only a MRE. I am well aware that the function being differentiated is useless. My application is for memoization. Here's a more realistic MRE:

using DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings
using PreallocationTools
use_sparse = false # works with `false`, crashes with `true`
backend = if use_sparse
    AutoSparse(
        AutoForwardDiff(); 
        sparsity_detector  = TracerSparsityDetector(), 
        coloring_algorithm = GreedyColoringAlgorithm()
    )
else 
    AutoForwardDiff()
end
i_counter, n = 0, 3
last_x_cache, last_y_cache = DiffCache(zeros(n)), DiffCache(zeros(n))
function f!(y, x::Vector{T}) where T
    last_x, last_y = get_tmp(last_x_cache, T), get_tmp(last_y_cache, T)
    if last_x != x
        # some heavy computations here:
        y .= x.^2 .+ 3x
        global i_counter += 1
        last_x .= x
        last_y .= y
    else
        y .= last_y
    end
    return y
end
x, y, J = ones(n), zeros(n), zeros(n, n)
map(i -> jacobian!(f!, y, J, backend, x), 1:1000)
display(J)        # J = diag{5, 5, 5}    
@show i_counter;  # executed only once, memoization works

If you set use_spase=true it will throw an error similar as above. I understand now that's not really possible with current DiffCache design, since it contains only floats and duals, not sparsity tracers. But I found a workaround with DifferentiationInterface.Cache, so I no longer need this feature.

edit : maybe you should mention in the Doc that we cannot assume that the initial value of DiffCache will be the same as the one used at its construction.

@ChrisRackauckas
Copy link
Member

That's not very a professional phrasing. Maybe you should read again your SciML Community Standards, but let's put this aside.

I'm not entirely sure what was read into that but there was no ill-intent with it. I just meant that the function does not have a sensical definition because the caches are never written into, so the outputs are random with an undefined probability distribution. That would be considered a user error here and the error message is correctly saying that an uninitialized value is accessed. It's the same issue as x = Vector{Any}(undef, 1); x[1] which will error with the same message for the same reason.

I understand now that's not really possible with current DiffCache design, since it contains only floats and duals, not sparsity tracers

No that is incorrect. There is a fallback Any cache that is used for the other cases, which is how it supports sparsity tracing with symbolics.

https://github.com/adrhill/SparseConnectivityTracer.jl/blob/main/src/tracers.jl#L15

With AbstractTracer defined as real it should work exactly the same way Symbolics does.

If you set use_spase=true it will throw an error similar as above. I understand now that's not really possible

This again is a user error not a package bug, but it's a bit more subtle. The issue is that x_last is not defined on the first call, so last_x != x is not defined. It will throw this same error because x_last is undef memory and you would then be doing the same thing as x = Vector{Any}(undef, 1); x[1]. If you did not error, i.e. used @inbounds here, then last_x != x would be random, though most of the time it would be true, and the function would do what you wanted, but you could randomly get that the first call has last_x == x (for example, if you ran the function in a test loop during a benchmark and it pulled retrieved memory after GC) in which case you'd also get y .= last_y and so the output would be random. So the error is guarding against unsafe memory behavior and is correctly triggering here.

To make your function robust to x_last not being defined on the first call, you need to check for assignment, i.e.:

julia> x = Vector{Any}(undef, 1)
1-element Vector{Any}:
 #undef

julia> isassigned(x,1)
false

In other words, here's the form of the function that would correctly be memoizing without the incorrect edge cases:

function f!(y, x::Vector{T}) where T
    last_x, last_y = get_tmp(last_x_cache, T), get_tmp(last_y_cache, T)
   @asset !isempty(last_x) # need to make sure there is at least one element 
   if !isassigned(last_x, 1) || last_x != x # Short circuit || is required to avoid error access
        # some heavy computations here:
        y .= x.^2 .+ 3x
        global i_counter += 1
        last_x .= x
        last_y .= y
    else
        y .= last_y
    end
    return y
end

And that's a version that should just work fine with sparsity detection since it will not check the contents of last_x before writing to it.

Now again, in "most" scenarios your code would look like it was okay so it might not seem like a user error, but you do have to be careful about uninitialized memory since assuming that it will be zero is incorrect, as can be verified by something like this:

julia> for i in 1:10
           x = Vector{UInt8}(undef, 1)
           @show @inbounds x[1]
       end
#= REPL[54]:3 =# @inbounds(x[1]) = 0xf8
#= REPL[54]:3 =# @inbounds(x[1]) = 0xff
#= REPL[54]:3 =# @inbounds(x[1]) = 0x00
#= REPL[54]:3 =# @inbounds(x[1]) = 0x10
#= REPL[54]:3 =# @inbounds(x[1]) = 0x00
#= REPL[54]:3 =# @inbounds(x[1]) = 0x10
#= REPL[54]:3 =# @inbounds(x[1]) = 0xff
#= REPL[54]:3 =# @inbounds(x[1]) = 0xd0
#= REPL[54]:3 =# @inbounds(x[1]) = 0x40
#= REPL[54]:3 =# @inbounds(x[1]) = 0x00

Did that elaboration help?

@franckgaga
Copy link
Author

franckgaga commented Mar 21, 2025

Yes thanks for the elaboration, it is more clear now. In short, the memory is never initialized with DiffCache, and it is true for floats, duals and sparsity detectors.

So the issue is in the documentation IMO. It should be mentioned in https://docs.sciml.ai/PreallocationTools/stable/#DiffCache or in https://docs.sciml.ai/PreallocationTools/stable/preallocationtools/#PreallocationTools.DiffCache that memory is never initialized, and the supplied values at construction are simply ignored. For exemple, I thought that:

last_x_cache, last_y_cache = DiffCache(fill(NaN, n)), DiffCache(zeros(n))

would be a simple fix but now I know that it would not with your clarification.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

3 participants