Memory overwrite

I am working towards creating an MWE of a UODE problem based on systems of equations. In the course of whittling down my problem, I came across the following:

     println("++> trajectories: ", trajectories)    # Train the model from scratch
    println("--> trajectories: ", trajectories)     # 1

These are two print statements, one after the other, printing the same variable, and I get different OUTPUTS. If I were programming in C++ I would think I overwrote memory. If my code was multithreaded, I could understand this as well. However, the code is single-threaded, and the problem is repeatable. My code is still too long to post, but if anybody has any preliminary ideas, I look forward to hearing them. Thanks.

I take it back. The code uses EnsembleSolve, which is multi-threaded.

This seems very strange. We would really need to see more code to figure out what’s going on (and to check if this is reproducible).

I suggest that you try to systematically reduce your code to a minimal working example (MWE), which should be as minimal as possible but should still display the suspicious behavior. (Note that very often the process of creating the MWE is already very instructive and reveals what the issue is.)

Is it possible that trajectories is of a type with a strange method that mutates the state of trajectories?

Here is my first attempt at an MWE (52 lines), with the same error that started me down this path. I will list the code and the error. The error generates using standard Julia 1.8.5 with and without the DaemonMode.jl module. I have removed the neural network simplified my solvers to one-liners. I solve the matrix equation dx/dt = 0 where x is a 3x3 static matrix with zero initial conditions. One thing I have not done is remove all unneeded libraries from the toml file. But I doubt that is the issue.
Any insight is appreciated.

function dudt_univ_opt(u, p, t)
    du = @SMatrix [0. 0. 0.; 0. 0. 0.; 0. 0. 0.]  # Static array

function single_run()
    # Define the simple shear deformation protocol (all components except v21(t))
    v11(t) = 0

    # Initial conditions and time span
    tspan = (0.0f0, 1.0f0)
    σ0 = SA[0.  0.  0.; 0.  0.  0.; 0.  0.  0.]  # Generates set_index! error

    p_system = Float32[0., 0.]  # change 2023-02-26_15:40
    θi = [0.]

    # Loss function closure (first parameter: concatenate all parameters)
    loss_fn(θ) = loss_univ([θ; p_system], tspan, σ0, 1)
        adtype = Optimization.AutoZygote()
        optf = Optimization.OptimizationFunction((x,p)->loss_fn(x),  adtype)
        optprob = Optimization.OptimizationProblem(optf, θi)
        parameter_res = Optimization.solve(optprob, Optimisers.AMSGrad(), sensealg=ReverseDiffVJP(true), allow_f_increases=false, maxiters=20)

function ensemble_solve(θ, ensemble, tspans, σ0, trajectories)
    θ = [1.]
    prob = ODEProblem(dudt_univ_opt, σ0, tspans[1], θ)

    function prob_func(prob, i, repeat)
            dudt_remade(u, p, t) = dudt_univ_opt(u, p ,t)
        remake(prob, f=dudt_remade, tspan=tspans[i])

    ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
    sim = solve(ensemble_prob, Tsit5(), ensemble, trajectories=trajectories, saveat=0.2)

function loss_univ(θ, tspans, σ0, trajectories)
    loss = 0
    results = ensemble_solve(θ, EnsembleThreads(), tspans, σ0, trajectories)
    xxx = results[1][2,1,:]  # 2nd row, 1st column: σ12 = σ21
    σ12_pred = results[1][2,1,:]  # 2nd row, 1st column: σ12 = σ21
    loss += sum(abs2, σ12_pred )  # Remove this line and error changes
    return loss


and the error:

ERROR: LoadError: setindex!(::StaticArraysCore.SMatrix{3, 3, Float64, 9}, value, ::Int) is not defined.
 [1] error at ./error.jl:35
 [2] setindex! at /Users/erlebach/.julia/packages/StaticArrays/pTgFe/src/indexing.jl:3
 [3] macro expansion at /Users/erlebach/.julia/packages/StaticArrays/pTgFe/src/indexing.jl:66
 [4] _setindex!_scalar at /Users/erlebach/.julia/packages/StaticArrays/pTgFe/src/indexing.jl:46
 [5] setindex! at /Users/erlebach/.julia/packages/StaticArrays/pTgFe/src/indexing.jl:42
 [6] setindex! at /Users/erlebach/.julia/packages/RecursiveArrayTools/VzH8Y/src/vector_of_array.jl:336
 [7] macro expansion at ./multidimensional.jl:946
 [8] macro expansion at ./cartesian.jl:64
 [9] macro expansion at ./multidimensional.jl:941
 [10] _unsafe_setindex! at ./multidimensional.jl:953
 [11] _setindex! at ./multidimensional.jl:930
 [12] setindex! at ./abstractarray.jl:1344
 [13] AbstractVectorOfArray_getindex_adjoint at /Users/erlebach/.julia/packages/RecursiveArrayTools/VzH8Y/src/zygote.jl:121
 [14] #108#back at /Users/erlebach/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [15] Pullback at /Users/erlebach/src/2022/rude/giesekus/GE_rude.jl/rude_optimized/mwe/rude_functions.jl:48
 [16] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [17] Pullback at /Users/erlebach/src/2022/rude/giesekus/GE_rude.jl/rude_optimized/mwe/rude_functions.jl:24
 [18] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [19] Pullback at /Users/erlebach/src/2022/rude/giesekus/GE_rude.jl/rude_optimized/mwe/rude_functions.jl:26
 [20] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [21] #208 at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/lib/lib.jl:206
 [22] #2066#back at /Users/erlebach/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [23] Pullback at /Users/erlebach/.julia/packages/SciMLBase/gTrkJ/src/scimlfunctions.jl:3904
 [24] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [25] #208 at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/lib/lib.jl:206
 [26] #2066#back at /Users/erlebach/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [27] Pullback at /Users/erlebach/.julia/packages/Optimization/XjqVZ/src/function/zygote.jl:30
 [28] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [29] #208 at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/lib/lib.jl:206
 [30] #2066#back at /Users/erlebach/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [31] Pullback at /Users/erlebach/.julia/packages/Optimization/XjqVZ/src/function/zygote.jl:34
 [32] Pullback at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface2.jl:0
 [33] #60 at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface.jl:45
 [34] gradient at /Users/erlebach/.julia/packages/Zygote/g2w9o/src/compiler/interface.jl:97
 [35] #157 at /Users/erlebach/.julia/packages/Optimization/XjqVZ/src/function/zygote.jl:32
 [36] macro expansion at /Users/erlebach/.julia/packages/OptimizationOptimisers/KGKWE/src/OptimizationOptimisers.jl:36
 [37] macro expansion at /Users/erlebach/.julia/packages/Optimization/XjqVZ/src/utils.jl:37
 [38] #__solve#1 at /Users/erlebach/.julia/packages/OptimizationOptimisers/KGKWE/src/OptimizationOptimisers.jl:35
 [39] #solve#552 at /Users/erlebach/.julia/packages/SciMLBase/gTrkJ/src/solve.jl:85
 [40] single_run at /Users/erlebach/src/2022/rude/giesekus/GE_rude.jl/rude_optimized/mwe/rude_functions.jl:28
 [41] top-level scope at /Users/erlebach/src/2022/rude/giesekus/GE_rude.jl/rude_optimized/mwe/rude_functions.jl:53
 [42] eval at ./boot.jl:368

And finally, the project.toml file:

ArgMacros = "dbc42088-9de8-42a0-8ec8-2cd114e1ea3e"
BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DaemonMode = "d749ddd5-2b29-4920-8305-6ff5a704e36e"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
DiffEqFlux = "aae7a2af-3d4f-5e19-a356-7da93b79d9d0"
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
NPZ = "15e1cf62-19b3-5cfa-8e77-841668bca605"
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Pickle = "fbb45041-c46e-462f-888f-7c521cafbc2c"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

I’m not sure I understand how this relates to your original issue? Here you seem to show some code which produces an error related to trying to mutate a static matrix, and without any printing?

1 Like

I discovered the to problem as I started to reduce my code towards the MWE. Perhaps I should have started a new thread.