Improving performance of dynamic programming problem of households in a labor market model

General Usage Performance
22

Hi. Thanks a lot. I will print the results from the profiler. I am quite bad at interpreting this, but it is clear that interpolations is costly. There may be other concerns as well, however.

I need to truncate the output in order to fit it within the space limit.

 
   75╎    ╎    55961 @Optim\src\maximize.jl:17; (::Optim.var"#127#128"{var"#38#40"{Int64, Int64, Int64, var"#V_fun#43"{F...
    4╎    ╎     4     @Base\float.jl:320; -(x::Float64)
  420╎    ╎     55882 ...Research\LUS heterogeneous agents\Julia\MWE.jl:211; (::var"#38#40"{Int64, Int64, Int64, var"#V_fun#43"{Float64, Vector{Flo...        
  
  
 1465╎    ╎    ╎    ╎    ╎     1465  @Base\math.jl:899; ^
 
  264╎    ╎    ╎    ╎    ╎ 264   @Base\boot.jl:448; Array
    1╎    ╎    ╎    ╎    1     @Base\set.jl:391; allunique(C::Vector{Float64})
    7╎    ╎    ╎    ╎    7     @Base\set.jl:397; allunique(C::Vector{Float64})
   14╎    ╎    ╎    ╎    287   @Base\set.jl:399; allunique(C::Vector{Float64})
    1╎    ╎    ╎    ╎     1     @Base\Base.jl:0; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   37╎    ╎    ╎    ╎     37    @Base\dict.jl:305; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎     11    @Base\dict.jl:307; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    4╎    ╎    ╎    ╎    ╎ 4     @Base\Base.jl:33; getproperty
    7╎    ╎    ╎    ╎    ╎ 7     @Base\array.jl:197; length
     ╎    ╎    ╎    ╎     2     @Base\dict.jl:309; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    2╎    ╎    ╎    ╎    ╎ 2     @Base\Base.jl:33; getproperty
     ╎    ╎    ╎    ╎     126   @Base\dict.jl:310; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 126   @Base\dict.jl:169; hashindex
 
  254╎    ╎    ╎    ╎    ╎     254   @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎   6     @Base\array.jl:504; zeros
    6╎    ╎    ╎    ╎    ╎    6     @Base\array.jl:406; fill!
     ╎    ╎    ╎    ╎    ╎ 588   @Base\dict.jl:193; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  588╎    ╎    ╎    ╎    ╎  588   @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎ 1000  @Base\dict.jl:194; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
 1000╎    ╎    ╎    ╎    ╎  1000  @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎ 2     @Base\dict.jl:199; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  
     ╎    ╎    ╎    ╎    ╎ 192   @Base\dict.jl:210; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  192╎    ╎    ╎    ╎    ╎  192   @Base\array.jl:839; setindex!
 
   23╎    ╎    ╎    ╎  24    @Interpolations\src\gridded\gridded.jl:47; Interpolations.GriddedInterpolation(#unused#::Type{Float64}, knots:...
    1╎    ╎    ╎    ╎   1     @Interpolations\src\gridded\gridded.jl:18; GriddedInterpolation
     ╎    ╎    ╎   8     ...rpolations\src\extrapolation\extrapolation.jl:48; Extrapolation
     ╎    ╎    ╎    3     @Interpolations\src\Interpolations.jl:147; bounds
     ╎    ╎    ╎     2     @Interpolations\src\gridded\gridded.jl:175; lbounds
     ╎    ╎    ╎    ╎ 2     @Base\broadcast.jl:883; materialize
     ╎    ╎    ╎    ╎  2     @Base\broadcast.jl:1098; copy
     ╎    ╎    ╎    ╎   2     @Base\ntuple.jl:48; ntuple
     ╎    ╎    ╎    ╎    2     @Base\broadcast.jl:1098; #19
     ╎    ╎    ╎    ╎     1     @Base\broadcast.jl:620; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎ 1     @Base\broadcast.jl:645; _getindex
     ╎    ╎    ╎    ╎    ╎  1     @Base\broadcast.jl:598; _broadcast_getindex
    1╎    ╎    ╎    ╎    ╎   1     @Base\tuple.jl:29; getindex
     ╎    ╎    ╎    ╎     1     @Base\broadcast.jl:621; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎ 1     @Base\broadcast.jl:648; _broadcast_getindex_evalf
     ╎    ╎    ╎    ╎    ╎  1     @Base\abstractarray.jl:368; first
    1╎    ╎    ╎    ╎    ╎   1     @Base\array.jl:801; getindex
     ╎    ╎    ╎     1     @Interpolations\src\gridded\gridded.jl:176; ubounds
     ╎    ╎    ╎    ╎ 1     @Base\broadcast.jl:883; materialize
     ╎    ╎    ╎    ╎  1     @Base\broadcast.jl:1098; copy
     ╎    ╎    ╎    ╎   1     @Base\ntuple.jl:48; ntuple
     ╎    ╎    ╎    ╎    1     @Base\broadcast.jl:1098; #19
     ╎    ╎    ╎    ╎     1     @Base\broadcast.jl:621; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎ 1     @Base\broadcast.jl:648; _broadcast_getindex_evalf
     ╎    ╎    ╎    ╎    ╎  1     @Base\abstractarray.jl:437; last
    1╎    ╎    ╎    ╎    ╎   1     @Base\array.jl:801; getindex
     ╎    ╎    ╎    5     ...polations\src\extrapolation\extrapolation.jl:102; inbounds_position
     ╎    ╎    ╎     5     ...polations\src\extrapolation\extrapolation.jl:108; inbounds_index
     ╎    ╎    ╎    ╎ 5     ...olations\src\extrapolation\extrapolation.jl:128; inbounds_index
     ╎    ╎    ╎    ╎  2     ...olations\src\extrapolation\extrapolation.jl:131; inbounds_index
     ╎    ╎    ╎    ╎   2     ...olations\src\extrapolation\extrapolation.jl:148; minp
    2╎    ╎    ╎    ╎    2     @Base\math.jl:726; min
     ╎    ╎    ╎    ╎  3     ...olations\src\extrapolation\extrapolation.jl:150; maxp
    3╎    ╎    ╎    ╎   3     @Base\math.jl:722; max
     ╎    ╎    ╎   20    ...rpolations\src\extrapolation\extrapolation.jl:49; Extrapolation
     ╎    ╎    ╎    17    @Interpolations\src\gridded\indexing.jl:4; GriddedInterpolation
     ╎    ╎    ╎     17    @Interpolations\src\b-splines\indexing.jl:66; weightedindexes
     ╎    ╎    ╎    ╎ 17    @Interpolations\src\b-splines\indexing.jl:70; map3argf
     ╎    ╎    ╎    ╎  15    @Interpolations\src\gridded\indexing.jl:55; weightedindex_parts
     ╎    ╎    ╎    ╎   15    @Interpolations\src\gridded\indexing.jl:51; find_knot_index
     ╎    ╎    ╎    ╎    15    @Base\sort.jl:295; searchsortedfirst
    1╎    ╎    ╎    ╎     1     @Base\sort.jl:0; searchsortedfirst
     ╎    ╎    ╎    ╎     1     @Base\sort.jl:183; searchsortedfirst
     ╎    ╎    ╎    ╎    ╎ 1     @Base\sort.jl:170; midpoint
    1╎    ╎    ╎    ╎    ╎  1     @Base\int.jl:86; -
     ╎    ╎    ╎    ╎     13    @Base\sort.jl:184; searchsortedfirst
     ╎    ╎    ╎    ╎    ╎ 13    @Base\ordering.jl:109; lt
   13╎    ╎    ╎    ╎    ╎  13    @Base\float.jl:381; isless
     ╎    ╎    ╎    ╎  2     @Interpolations\src\gridded\indexing.jl:56; weightedindex_parts
     ╎    ╎    ╎    ╎   1     @Interpolations\src\gridded\indexing.jl:61; weightedindex_parts2
    1╎    ╎    ╎    ╎    1     @Base\math.jl:65; clamp
     ╎    ╎    ╎    ╎   1     @Interpolations\src\gridded\indexing.jl:62; weightedindex_parts2
     ╎    ╎    ╎    ╎    1     @Interpolations\src\gridded\indexing.jl:72; weightedindex
     ╎    ╎    ╎    ╎     1     @Interpolations\src\utils.jl:10; fmap
     ╎    ╎    ╎    ╎    ╎ 1     @Interpolations\src\utils.jl:11; _fmap
     ╎    ╎    ╎    ╎    ╎  1     @Interpolations\src\b-splines\linear.jl:54; value_weights
     ╎    ╎    ╎    ╎    ╎   1     @Base\promotion.jl:323; -
    1╎    ╎    ╎    ╎    ╎    1     @Base\float.jl:329; -
     ╎    ╎    ╎    3     @Interpolations\src\gridded\indexing.jl:5; GriddedInterpolation
     ╎    ╎    ╎     1     @Base\abstractarray.jl:1170; getindex
     ╎    ╎    ╎    ╎ 1     @Interpolations\src\Interpolations.jl:285; _getindex
     ╎    ╎    ╎    ╎  1     @Interpolations\src\Interpolations.jl:319; interp_getindex
     ╎    ╎    ╎    ╎   1     @Interpolations\src\Interpolations.jl:329; macro expansion
    1╎    ╎    ╎    ╎    1     @Base\float.jl:326; +
     ╎    ╎    ╎     2     @Interpolations\src\gridded\gridded.jl:141; coefficients
    2╎    ╎    ╎    ╎ 2     @Base\Base.jl:33; getproperty
    5╎    ╎    ╎ 49626 ...esearch\LUS heterogeneous agents\Julia\MWE.jl:147; bellman(a_prime::Float64, ia::Int64, is::Int64, iom::Int64, V_fun::var...        
    2╎    ╎    ╎  2     @Base\array.jl:0; (::var"#V_fun#43"{Float64, Vector{Float64}})(x::Float64, is::Int64, i...
   26╎    ╎    ╎  26    @Base\float.jl:332; *
    6╎    ╎    ╎  6     @Base\float.jl:326; +
    2╎    ╎    ╎  2     @Base\float.jl:335; /
  373╎    ╎    ╎  373   @Base\math.jl:899; ^
    3╎    ╎    ╎  3     @Base\math.jl:900; ^
    7╎    ╎    ╎  11    @Base\math.jl:722; max
    4╎    ╎    ╎   4     @Base\floatfuncs.jl:15; signbit
     ╎    ╎    ╎  10    @Base\range.jl:674; iterate
   10╎    ╎    ╎   10    @Base\promotion.jl:410; ==
    7╎    ╎    ╎  7     ...esearch\LUS heterogeneous agents\Julia\MWE.jl:0; (::var"#V_fun#43"{Float64, Vector{Float64}})(x::Float64, is::Int64, i...
   39╎    ╎    ╎  49181 ...esearch\LUS heterogeneous agents\Julia\MWE.jl:237; (::var"#V_fun#43"{Float64, Vector{Float64}})(x::Float64, is::Int64, i...        
     ╎    ╎    ╎   14904 @Base\broadcast.jl:883; materialize
     ╎    ╎    ╎    14904 @Base\broadcast.jl:908; copy
     ╎    ╎    ╎     13761 @Base\broadcast.jl:936; copyto!
    3╎    ╎    ╎    ╎ 3     @Base\broadcast.jl:0; copyto!
     ╎    ╎    ╎    ╎ 11    @Base\broadcast.jl:972; copyto!
     ╎    ╎    ╎    ╎  11    @Base\tuple.jl:362; ==
     ╎    ╎    ╎    ╎   11    @Base\tuple.jl:366; _eq
    9╎    ╎    ╎    ╎    9     @Base\range.jl:802; ==
     ╎    ╎    ╎    ╎    2     @Base\range.jl:804; ==
    2╎    ╎    ╎    ╎     2     @Base\promotion.jl:410; ==
     ╎    ╎    ╎    ╎ 8     @Base\broadcast.jl:980; copyto!
     ╎    ╎    ╎    ╎  8     @Base\broadcast.jl:963; preprocess
     ╎    ╎    ╎    ╎   8     @Base\broadcast.jl:966; preprocess_args
     ╎    ╎    ╎    ╎    8     @Base\broadcast.jl:967; preprocess_args
     ╎    ╎    ╎    ╎     8     @Base\broadcast.jl:963; preprocess
     ╎    ╎    ╎    ╎    ╎ 8     @Base\broadcast.jl:966; preprocess_args
     ╎    ╎    ╎    ╎    ╎  8     @Base\broadcast.jl:964; preprocess
     ╎    ╎    ╎    ╎    ╎   8     @Base\broadcast.jl:957; broadcast_unalias
     ╎    ╎    ╎    ╎    ╎    8     @Base\abstractarray.jl:1349; unalias
    2╎    ╎    ╎    ╎    ╎     2     @Base\abstractarray.jl:0; mightalias
    5╎    ╎    ╎    ╎    ╎     6     @Base\abstractarray.jl:1384; mightalias
     ╎    ╎    ╎    ╎    ╎    ╎ 1     @Base\abstractarray.jl:1408; dataids
     ╎    ╎    ╎    ╎    ╎    ╎  1     @Base\abstractarray.jl:1116; pointer
    1╎    ╎    ╎    ╎    ╎    ╎   1     @Base\pointer.jl:65; unsafe_convert
     ╎    ╎    ╎    ╎ 13739 @Base\broadcast.jl:983; copyto!
  118╎    ╎    ╎    ╎  118   @Base\simdloop.jl:0; macro expansion
     ╎    ╎    ╎    ╎  3     @Base\simdloop.jl:72; macro expansion
    3╎    ╎    ╎    ╎   3     @Base\int.jl:83; <
     ╎    ╎    ╎    ╎  13618 @Base\simdloop.jl:77; macro expansion
     ╎    ╎    ╎    ╎   13618 @Base\broadcast.jl:984; macro expansion
 1070╎    ╎    ╎    ╎    1070  @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎    12548 @Base\broadcast.jl:575; getindex
     ╎    ╎    ╎    ╎     12449 @Base\broadcast.jl:620; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎ 12449 @Base\broadcast.jl:644; _getindex
     ╎    ╎    ╎    ╎    ╎  12449 @Base\broadcast.jl:645; _getindex
     ╎    ╎    ╎    ╎    ╎   4     @Base\broadcast.jl:620; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎    4     @Base\broadcast.jl:644; _getindex
     ╎    ╎    ╎    ╎    ╎     4     @Base\broadcast.jl:614; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎    ╎ 4     @Base\subarray.jl:309; getindex
    4╎    ╎    ╎    ╎    ╎    ╎  4     @Base\array.jl:801; getindex
     ╎    ╎    ╎    ╎    ╎   12445 @Base\broadcast.jl:621; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎    12445 @Base\broadcast.jl:648; _broadcast_getindex_evalf
12362╎    ╎    ╎    ╎    ╎     12362 @Base\math.jl:899; ^
   82╎    ╎    ╎    ╎    ╎     83    @Base\math.jl:900; ^
     ╎    ╎    ╎    ╎    ╎    ╎ 1     @Base\float.jl:449; isnan
    1╎    ╎    ╎    ╎    ╎    ╎  1     @Base\float.jl:368; !=
     ╎    ╎    ╎    ╎     99    @Base\broadcast.jl:621; _broadcast_getindex
     ╎    ╎    ╎    ╎    ╎ 99    @Base\broadcast.jl:648; _broadcast_getindex_evalf
   99╎    ╎    ╎    ╎    ╎  99    @Base\float.jl:332; *
     ╎    ╎    ╎     1143  @Base\broadcast.jl:196; similar
     ╎    ╎    ╎    ╎ 1143  @Base\broadcast.jl:197; similar
     ╎    ╎    ╎    ╎  1143  @Base\abstractarray.jl:784; similar
     ╎    ╎    ╎    ╎   1143  @Base\abstractarray.jl:785; similar
     ╎    ╎    ╎    ╎    1143  @Base\boot.jl:465; Array
     ╎    ╎    ╎    ╎     1143  @Base\boot.jl:457; Array
 1143╎    ╎    ╎    ╎    ╎ 1143  @Base\boot.jl:448; Array
     ╎    ╎    ╎   4     @Base\subarray.jl:176; view
     ╎    ╎    ╎    4     @Base\indices.jl:324; to_indices
     ╎    ╎    ╎     4     @Base\abstractarray.jl:89; axes
     ╎    ╎    ╎    ╎ 4     @Base\array.jl:135; size
     ╎    ╎    ╎    ╎  4     @Base\ntuple.jl:50; ntuple
     ╎    ╎    ╎    ╎   4     @Base\array.jl:135; #85
    4╎    ╎    ╎    ╎    4     @Base\array.jl:132; size
    9╎    ╎    ╎   12    @Base\subarray.jl:177; view
    1╎    ╎    ╎    1     @Base\abstractarray.jl:0; checkbounds
    2╎    ╎    ╎    2     @Base\abstractarray.jl:616; checkbounds
     ╎    ╎    ╎   33975 @Interpolations\src\convenience-constructors.jl:9; LinearInterpolation##kw
     ╎    ╎    ╎    33975 @Interpolations\src\convenience-constructors.jl:9; #LinearInterpolation#53
    3╎    ╎    ╎     33975 @Interpolations\src\gridded\gridded.jl:166; interpolate
   31╎    ╎    ╎    ╎ 31    @Interpolations\src\gridded\gridded.jl:148; interpolate(#unused#::Type{Float64}, #unused#::Type{Float64}, knots...
 3363╎    ╎    ╎    ╎ 33941 @Interpolations\src\gridded\gridded.jl:149; interpolate(#unused#::Type{Float64}, #unused#::Type{Float64}, knots...
  311╎    ╎    ╎    ╎  311   @Interpolations\src\gridded\gridded.jl:37; Interpolations.GriddedInterpolation(#unused#::Type{Float64}, knots:...
     ╎    ╎    ╎    ╎  29727 @Interpolations\src\gridded\gridded.jl:40; Interpolations.GriddedInterpolation(#unused#::Type{Float64}, knots:...
  
   55╎    ╎    ╎    ╎    ╎ 55    @Base\array.jl:0; issorted
  104╎    ╎    ╎    ╎    ╎ 104   @Base\sort.jl:0; issorted
     ╎    ╎    ╎    ╎    ╎ 2     @Base\sort.jl:57; issorted
     ╎    ╎    ╎    ╎    ╎  2     @Base\array.jl:777; iterate
     ╎    ╎    ╎    ╎    ╎   2     @Base\array.jl:777; iterate
    2╎    ╎    ╎    ╎    ╎    2     @Base\array.jl:801; getindex
   47╎    ╎    ╎    ╎    ╎ 47    @Base\sort.jl:61; issorted
   21╎    ╎    ╎    ╎    ╎ 650   @Base\sort.jl:63; issorted
     ╎    ╎    ╎    ╎    ╎  629   @Base\ordering.jl:109; lt
  629╎    ╎    ╎    ╎    ╎   629   @Base\float.jl:381; isless
     ╎    ╎    ╎    ╎    ╎ 25    @Base\sort.jl:65; issorted
     ╎    ╎    ╎    ╎    ╎  25    @Base\array.jl:777; iterate
     ╎    ╎    ╎    ╎    ╎   25    @Base\int.jl:448; <
   25╎    ╎    ╎    ╎    ╎    25    @Base\int.jl:441; <
    6╎    ╎    ╎    ╎   28818 @Interpolations\src\gridded\gridded.jl:76; check_gridded
    9╎    ╎    ╎    ╎    9     @Base\int.jl:0; allunique(C::Vector{Float64})
   30╎    ╎    ╎    ╎    30    @Base\set.jl:0; allunique(C::Vector{Float64})
   42╎    ╎    ╎    ╎    42    @Base\set.jl:384; allunique(C::Vector{Float64})
   12╎    ╎    ╎    ╎    2999  @Base\set.jl:385; allunique(C::Vector{Float64})
   29╎    ╎    ╎    ╎     29    @Base\dict.jl:88; Dict{Float64, Nothing}()
  389╎    ╎    ╎    ╎     2958  @Base\dict.jl:90; Dict{Float64, Nothing}()
     ╎    ╎    ╎    ╎    ╎ 538   @Base\array.jl:499; zeros
     ╎    ╎    ╎    ╎    ╎  521   @Base\array.jl:503; zeros
     ╎    ╎    ╎    ╎    ╎   521   @Base\boot.jl:457; Array
  521╎    ╎    ╎    ╎    ╎    521   @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎  17    @Base\array.jl:504; zeros
   14╎    ╎    ╎    ╎    ╎   17    @Base\array.jl:406; fill!
    3╎    ╎    ╎    ╎    ╎    3     @Base\array.jl:197; length
 2031╎    ╎    ╎    ╎    ╎ 2031  @Base\boot.jl:448; Array
 
   62╎    ╎    ╎    ╎    62    @Base\set.jl:397; allunique(C::Vector{Float64})
  134╎    ╎    ╎    ╎    2972  @Base\set.jl:399; allunique(C::Vector{Float64})
   68╎    ╎    ╎    ╎     68    @Base\Base.jl:0; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
  382╎    ╎    ╎    ╎     382   @Base\dict.jl:305; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎     73    @Base\dict.jl:307; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   37╎    ╎    ╎    ╎    ╎ 37    @Base\Base.jl:33; getproperty
   36╎    ╎    ╎    ╎    ╎ 36    @Base\array.jl:197; length
     ╎    ╎    ╎    ╎     16    @Base\dict.jl:309; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   16╎    ╎    ╎    ╎    ╎ 16    @Base\Base.jl:33; getproperty
     ╎    ╎    ╎    ╎     1394  @Base\dict.jl:310; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 1394  @Base\dict.jl:169; hashindex
  130╎    ╎    ╎    ╎    ╎  1314  @Base\hashing.jl:18; hash
   64╎    ╎    ╎    ╎    ╎   64    @Base\float.jl:0; hash(x::Float64, h::UInt64)
  103╎    ╎    ╎    ╎    ╎   103   @Base\float.jl:465; hash(x::Float64, h::UInt64)
  109╎    ╎    ╎    ╎    ╎   150   @Base\float.jl:467; hash(x::Float64, h::UInt64)
   41╎    ╎    ╎    ╎    ╎    41    @Base\float.jl:374; <=
   17╎    ╎    ╎    ╎    ╎   17    @Base\float.jl:468; hash(x::Float64, h::UInt64)
  133╎    ╎    ╎    ╎    ╎   216   @Base\float.jl:469; hash(x::Float64, h::UInt64)
     ╎    ╎    ╎    ╎    ╎    83    @Base\operators.jl:129; isequal
     ╎    ╎    ╎    ╎    ╎     83    @Base\operators.jl:125; signequal
   51╎    ╎    ╎    ╎    ╎    ╎ 51    @Base\floatfuncs.jl:15; signbit
   32╎    ╎    ╎    ╎    ╎    ╎ 32    @Base\promotion.jl:410; ==
    1╎    ╎    ╎    ╎    ╎   12    @Base\float.jl:470; hash(x::Float64, h::UInt64)
     ╎    ╎    ╎    ╎    ╎    11    @Base\hashing.jl:71; hash
     ╎    ╎    ╎    ╎    ╎     11    @Base\hashing.jl:62; hash_uint64
     ╎    ╎    ╎    ╎    ╎    ╎ 2     @Base\hashing.jl:29; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  1     @Base\int.jl:464; <<
    1╎    ╎    ╎    ╎    ╎    ╎   1     @Base\int.jl:457; <<
    1╎    ╎    ╎    ╎    ╎    ╎  1     @Base\int.jl:285; ~
     ╎    ╎    ╎    ╎    ╎    ╎ 3     @Base\hashing.jl:30; hash_64_64
    3╎    ╎    ╎    ╎    ╎    ╎  3     @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎    ╎ 2     @Base\hashing.jl:32; hash_64_64
    2╎    ╎    ╎    ╎    ╎    ╎  2     @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎    ╎ 1     @Base\hashing.jl:33; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  1     @Base\operators.jl:560; +
    1╎    ╎    ╎    ╎    ╎    ╎   1     @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    ╎ 3     @Base\hashing.jl:34; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  3     @Base\int.jl:462; >>
    3╎    ╎    ╎    ╎    ╎    ╎   3     @Base\int.jl:456; >>
   45╎    ╎    ╎    ╎    ╎   622   @Base\float.jl:480; hash(x::Float64, h::UInt64)
     ╎    ╎    ╎    ╎    ╎    515   @Base\hashing.jl:62; hash_uint64
     ╎    ╎    ╎    ╎    ╎     106   @Base\hashing.jl:29; hash_64_64
   21╎    ╎    ╎    ╎    ╎    ╎ 21    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    ╎ 5     @Base\int.jl:464; <<
    5╎    ╎    ╎    ╎    ╎    ╎  5     @Base\int.jl:457; <<
   80╎    ╎    ╎    ╎    ╎    ╎ 80    @Base\int.jl:285; ~
     ╎    ╎    ╎    ╎    ╎     73    @Base\hashing.jl:30; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎ 67    @Base\int.jl:462; >>
   67╎    ╎    ╎    ╎    ╎    ╎  67    @Base\int.jl:456; >>
    6╎    ╎    ╎    ╎    ╎    ╎ 6     @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎     25    @Base\hashing.jl:31; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎ 25    @Base\operators.jl:560; +
   25╎    ╎    ╎    ╎    ╎    ╎  25    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎     78    @Base\hashing.jl:32; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎ 76    @Base\int.jl:462; >>
   76╎    ╎    ╎    ╎    ╎    ╎  76    @Base\int.jl:456; >>
    2╎    ╎    ╎    ╎    ╎    ╎ 2     @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎     72    @Base\hashing.jl:33; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎ 72    @Base\operators.jl:560; +
   72╎    ╎    ╎    ╎    ╎    ╎  72    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎     64    @Base\hashing.jl:34; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎ 34    @Base\int.jl:462; >>
   34╎    ╎    ╎    ╎    ╎    ╎  34    @Base\int.jl:456; >>
   30╎    ╎    ╎    ╎    ╎    ╎ 30    @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎     97    @Base\hashing.jl:35; hash_64_64
   97╎    ╎    ╎    ╎    ╎    ╎ 97    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    30    @Base\int.jl:923; *
   30╎    ╎    ╎    ╎    ╎     30    @Base\int.jl:88; *
   32╎    ╎    ╎    ╎    ╎    32    @Base\int.jl:86; -
   27╎    ╎    ╎    ╎    ╎  27    @Base\int.jl:309; &
   30╎    ╎    ╎    ╎    ╎  30    @Base\int.jl:87; +
   23╎    ╎    ╎    ╎    ╎  23    @Base\int.jl:86; -
     ╎    ╎    ╎    ╎     23    @Base\dict.jl:312; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   23╎    ╎    ╎    ╎    ╎ 23    @Base\Base.jl:33; getproperty
  194╎    ╎    ╎    ╎     270   @Base\dict.jl:315; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 76    @Base\dict.jl:171; isslotempty
   13╎    ╎    ╎    ╎    ╎  13    @Base\Base.jl:0; getindex
   24╎    ╎    ╎    ╎    ╎  24    @Base\Base.jl:0; getproperty
   39╎    ╎    ╎    ╎    ╎  39    @Base\array.jl:801; getindex
     ╎    ╎    ╎    ╎     58    @Base\dict.jl:316; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   58╎    ╎    ╎    ╎    ╎ 58    @Base\int.jl:83; <
  208╎    ╎    ╎    ╎     224   @Base\dict.jl:319; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   16╎    ╎    ╎    ╎    ╎ 16    @Base\int.jl:85; -
   36╎    ╎    ╎    ╎     142   @Base\dict.jl:328; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   39╎    ╎    ╎    ╎    ╎ 39    @Base\Base.jl:0; getindex
   37╎    ╎    ╎    ╎    ╎ 37    @Base\array.jl:801; getindex
   30╎    ╎    ╎    ╎    ╎ 30    @Base\float.jl:0; isequal
     ╎    ╎    ╎    ╎     25    @Base\dict.jl:332; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   15╎    ╎    ╎    ╎    ╎ 15    @Base\int.jl:309; &
   10╎    ╎    ╎    ╎    ╎ 10    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎     37    @Base\dict.jl:333; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
   37╎    ╎    ╎    ╎    ╎ 37    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎     16    @Base\dict.jl:334; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 16    @Base\operators.jl:305; >
   16╎    ╎    ╎    ╎    ╎  16    @Base\int.jl:83; <
     ╎    ╎    ╎    ╎     7     @Base\dict.jl:337; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    7╎    ╎    ╎    ╎    ╎ 7     @Base\int.jl:83; <
     ╎    ╎    ╎    ╎     11    @Base\dict.jl:339; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 6     @Base\int.jl:462; >>
    6╎    ╎    ╎    ╎    ╎  6     @Base\int.jl:455; >>
    3╎    ╎    ╎    ╎    ╎ 5     @Base\promotion.jl:421; max
    2╎    ╎    ╎    ╎    ╎  2     @Base\int.jl:83; <
    6╎    ╎    ╎    ╎     14    @Base\dict.jl:341; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    8╎    ╎    ╎    ╎    ╎ 8     @Base\int.jl:83; <
   18╎    ╎    ╎    ╎     34    @Base\dict.jl:342; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
     ╎    ╎    ╎    ╎    ╎ 16    @Base\dict.jl:172; isslotfilled
   16╎    ╎    ╎    ╎    ╎  16    @Base\promotion.jl:410; ==
     ╎    ╎    ╎    ╎     2     @Base\dict.jl:343; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    2╎    ╎    ╎    ╎    ╎ 2     @Base\Base.jl:34; setproperty!
   25╎    ╎    ╎    ╎     32    @Base\dict.jl:344; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    7╎    ╎    ╎    ╎    ╎ 7     @Base\int.jl:85; -
     ╎    ╎    ╎    ╎     7     @Base\dict.jl:346; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    6╎    ╎    ╎    ╎    ╎ 6     @Base\int.jl:309; &
    1╎    ╎    ╎    ╎    ╎ 1     @Base\int.jl:87; +
     ╎    ╎    ╎    ╎     3     @Base\dict.jl:347; ht_keyindex2!(h::Dict{Float64, Nothing}, key::Float64)
    3╎    ╎    ╎    ╎    ╎ 3     @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    15    @Base\set.jl:400; allunique(C::Vector{Float64})
     ╎    ╎    ╎    ╎     15    @Base\operators.jl:305; >
   15╎    ╎    ╎    ╎    ╎ 15    @Base\int.jl:83; <
     ╎    ╎    ╎    ╎    22577 @Base\set.jl:401; allunique(C::Vector{Float64})
     ╎    ╎    ╎    ╎     145   @Base\dict.jl:356; _setindex!
   13╎    ╎    ╎    ╎    ╎ 13    @Base\Base.jl:33; getproperty
  132╎    ╎    ╎    ╎    ╎ 132   @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎     510   @Base\dict.jl:357; _setindex!
  116╎    ╎    ╎    ╎    ╎ 116   @Base\Base.jl:33; getproperty
  394╎    ╎    ╎    ╎    ╎ 394   @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎     89    @Base\dict.jl:358; _setindex!
   66╎    ╎    ╎    ╎    ╎ 66    @Base\Base.jl:33; getproperty
   23╎    ╎    ╎    ╎    ╎ 23    @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎     93    @Base\dict.jl:359; _setindex!
   12╎    ╎    ╎    ╎    ╎ 12    @Base\Base.jl:33; getproperty
   37╎    ╎    ╎    ╎    ╎ 37    @Base\Base.jl:34; setproperty!
   44╎    ╎    ╎    ╎    ╎ 44    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎     7     @Base\dict.jl:360; _setindex!
    7╎    ╎    ╎    ╎    ╎ 7     @Base\Base.jl:34; setproperty!
     ╎    ╎    ╎    ╎     120   @Base\dict.jl:361; _setindex!
   22╎    ╎    ╎    ╎    ╎ 22    @Base\Base.jl:33; getproperty
   98╎    ╎    ╎    ╎    ╎ 98    @Base\int.jl:83; <
     ╎    ╎    ╎    ╎     140   @Base\dict.jl:367; _setindex!
   52╎    ╎    ╎    ╎    ╎ 52    @Base\int.jl:88; *
     ╎    ╎    ╎    ╎    ╎ 27    @Base\int.jl:462; >>
   27╎    ╎    ╎    ╎    ╎  27    @Base\int.jl:455; >>
     ╎    ╎    ╎    ╎    ╎ 26    @Base\operators.jl:305; >
   26╎    ╎    ╎    ╎    ╎  26    @Base\int.jl:83; <
     ╎    ╎    ╎    ╎    ╎ 35    @Base\operators.jl:352; >=
   35╎    ╎    ╎    ╎    ╎  35    @Base\int.jl:442; <=
    7╎    ╎    ╎    ╎     21395 @Base\dict.jl:369; _setindex!
    1╎    ╎    ╎    ╎    ╎ 1     @Base\array.jl:0; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  273╎    ╎    ╎    ╎    ╎ 273   @Base\dict.jl:0; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
   82╎    ╎    ╎    ╎    ╎ 82    @Base\dict.jl:175; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
     ╎    ╎    ╎    ╎    ╎ 15    @Base\dict.jl:176; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
   15╎    ╎    ╎    ╎    ╎  15    @Base\Base.jl:33; getproperty

 3007╎    ╎    ╎    ╎    ╎     3007  @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎   61    @Base\array.jl:504; zeros
   47╎    ╎    ╎    ╎    ╎    61    @Base\array.jl:406; fill!
   14╎    ╎    ╎    ╎    ╎     14    @Base\array.jl:197; length
     ╎    ╎    ╎    ╎    ╎ 5026  @Base\dict.jl:193; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
 5026╎    ╎    ╎    ╎    ╎  5026  @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎ 8834  @Base\dict.jl:194; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
 8834╎    ╎    ╎    ╎    ╎  8834  @Base\boot.jl:448; Array
     ╎    ╎    ╎    ╎    ╎ 10    @Base\dict.jl:199; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
    6╎    ╎    ╎    ╎    ╎  10    @Base\range.jl:670; iterate
     ╎    ╎    ╎    ╎    ╎   4     @Base\range.jl:519; isempty
     ╎    ╎    ╎    ╎    ╎    4     @Base\operators.jl:305; >
    4╎    ╎    ╎    ╎    ╎     4     @Base\int.jl:83; <
   68╎    ╎    ╎    ╎    ╎ 186   @Base\dict.jl:200; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
   29╎    ╎    ╎    ╎    ╎  29    @Base\array.jl:801; getindex
   89╎    ╎    ╎    ╎    ╎  89    @Base\promotion.jl:410; ==
     ╎    ╎    ╎    ╎    ╎ 117   @Base\dict.jl:201; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  117╎    ╎    ╎    ╎    ╎  117   @Base\array.jl:801; getindex
     ╎    ╎    ╎    ╎    ╎ 1471  @Base\dict.jl:203; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
     ╎    ╎    ╎    ╎    ╎  1471  @Base\dict.jl:169; hashindex
    9╎    ╎    ╎    ╎    ╎   9     @Base\dict.jl:0; -
  115╎    ╎    ╎    ╎    ╎   1354  @Base\hashing.jl:18; hash
   75╎    ╎    ╎    ╎    ╎    75    @Base\float.jl:0; hash(x::Float64, h::UInt64)
   88╎    ╎    ╎    ╎    ╎    88    @Base\float.jl:465; hash(x::Float64, h::UInt64)
  127╎    ╎    ╎    ╎    ╎    159   @Base\float.jl:467; hash(x::Float64, h::UInt64)
   32╎    ╎    ╎    ╎    ╎     32    @Base\float.jl:374; <=
   17╎    ╎    ╎    ╎    ╎    17    @Base\float.jl:468; hash(x::Float64, h::UInt64)
  143╎    ╎    ╎    ╎    ╎    258   @Base\float.jl:469; hash(x::Float64, h::UInt64)
     ╎    ╎    ╎    ╎    ╎     115   @Base\operators.jl:129; isequal
     ╎    ╎    ╎    ╎    ╎    ╎ 115   @Base\operators.jl:125; signequal
   61╎    ╎    ╎    ╎    ╎    ╎  61    @Base\floatfuncs.jl:15; signbit
   54╎    ╎    ╎    ╎    ╎    ╎  54    @Base\promotion.jl:410; ==
     ╎    ╎    ╎    ╎    ╎    20    @Base\float.jl:470; hash(x::Float64, h::UInt64)
 
   58╎    ╎    ╎    ╎    ╎    622   @Base\float.jl:480; hash(x::Float64, h::UInt64)
     ╎    ╎    ╎    ╎    ╎     516   @Base\hashing.jl:62; hash_uint64
     ╎    ╎    ╎    ╎    ╎    ╎ 104   @Base\hashing.jl:29; hash_64_64
   16╎    ╎    ╎    ╎    ╎    ╎  16    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    ╎  12    @Base\int.jl:464; <<
   12╎    ╎    ╎    ╎    ╎    ╎   12    @Base\int.jl:457; <<
   76╎    ╎    ╎    ╎    ╎    ╎  76    @Base\int.jl:285; ~
     ╎    ╎    ╎    ╎    ╎    ╎ 73    @Base\hashing.jl:30; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  63    @Base\int.jl:462; >>
   63╎    ╎    ╎    ╎    ╎    ╎   63    @Base\int.jl:456; >>
   10╎    ╎    ╎    ╎    ╎    ╎  10    @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎    ╎ 21    @Base\hashing.jl:31; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  21    @Base\operators.jl:560; +
   21╎    ╎    ╎    ╎    ╎    ╎   21    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    ╎ 99    @Base\hashing.jl:32; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  94    @Base\int.jl:462; >>
   94╎    ╎    ╎    ╎    ╎    ╎   94    @Base\int.jl:456; >>
    5╎    ╎    ╎    ╎    ╎    ╎  5     @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎    ╎ 62    @Base\hashing.jl:33; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  62    @Base\operators.jl:560; +
   62╎    ╎    ╎    ╎    ╎    ╎   62    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎    ╎ 69    @Base\hashing.jl:34; hash_64_64
     ╎    ╎    ╎    ╎    ╎    ╎  52    @Base\int.jl:462; >>
   52╎    ╎    ╎    ╎    ╎    ╎   52    @Base\int.jl:456; >>
   17╎    ╎    ╎    ╎    ╎    ╎  17    @Base\int.jl:333; xor
     ╎    ╎    ╎    ╎    ╎    ╎ 88    @Base\hashing.jl:35; hash_64_64
   88╎    ╎    ╎    ╎    ╎    ╎  88    @Base\int.jl:87; +
     ╎    ╎    ╎    ╎    ╎     31    @Base\int.jl:923; *
   31╎    ╎    ╎    ╎    ╎    ╎ 31    @Base\int.jl:88; *
   17╎    ╎    ╎    ╎    ╎     17    @Base\int.jl:86; -
   23╎    ╎    ╎    ╎    ╎   23    @Base\int.jl:309; &
   85╎    ╎    ╎    ╎    ╎   85    @Base\int.jl:87; +
  223╎    ╎    ╎    ╎    ╎ 283   @Base\dict.jl:204; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
   39╎    ╎    ╎    ╎    ╎  39    @Base\array.jl:801; getindex
   
   78╎    ╎    ╎    ╎    ╎  78    @Base\int.jl:86; -
     ╎    ╎    ╎    ╎    ╎ 31    @Base\dict.jl:208; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
     ╎    ╎    ╎    ╎    ╎  31    @Base\operators.jl:305; >
   31╎    ╎    ╎    ╎    ╎   31    @Base\int.jl:83; <
     ╎    ╎    ╎    ╎    ╎ 103   @Base\dict.jl:209; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
  103╎    ╎    ╎    ╎    ╎  103   @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎    ╎ 1607  @Base\dict.jl:210; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
 1607╎    ╎    ╎    ╎    ╎  1607  @Base\array.jl:839; setindex!
     ╎    ╎    ╎    ╎    ╎ 15    @Base\dict.jl:214; rehash!(h::Dict{Float64, Nothing}, newsz::Int64)
     ╎    ╎    ╎    ╎    ╎  15    @Base\operators.jl:204; !=

  536╎    ╎    ╎    ╎  540   @Interpolations\src\gridded\gridded.jl:47; Interpolations.GriddedInterpolation(#unused#::Type{Float64}, knots:...
23

Hi, to clarify, I deleted may rows corresponding to a smaller number of samples to get the output to fit.

24

This is very concerning if it is true. the @unpack is supposed to lower to identical code as just accessing the properties, which is just “names” and should have no impact on performance?

Since in Julia 1.7 we can now use the (; a, b) = p instead of the @unpack a, b = p maybe that would show something diferent.

25

Hi. I tried LinearInterpolations.jl rather than Interpolations.jl and benchmarked the function update_value!. It runs about twice as fast, and the number of allocations drops from over the (seemingly incredible) 4GB to 950 MB. Now, that still seems like too much, but there is significant progress here.

So, I wonder if LinearInterpolations.jl is the best lightweight interpolator around?

Regards.

26

Another is GitHub - PumasAI/DataInterpolations.jl: A library of data interpolation and smoothing functions which might be better on allocations.

After that, the usual suspect in these things is the optimization algorithm. Perhaps the fortran code is using something far better for these circumstances that requires radically fewer function evaluations. For example, giving bounds , better itinital conditions, etc.

Another thing that sometimes happens comparing code is that the fortran optimizer might have a much looser convergence threshold and ends up being sped up partially because of that.

The last thing I will suggest is that you consider using a better fixed-point algorithm such as anderson iteration (which is in https://github.com/JuliaNLSolvers/NLsolve.jl#fixed-points for example) this can sometimes speed things up significantly.

But all of that is to say that while Julia is just as fast as fortran in principle, in practice it sometimes takes time to get there - so you need to decide how much investment is worthwhile

1 Like
27

That’s great that LinearInterpolations.jl saves some time. It looks like it allocates a tiny bit for a struct, and it would be best if even that were avoided. I also looked at @jlperla’s suggestion of DataInterpolations, and it seems to allocate but in a cache. It may thus be possible to update the cache with each iteration, and avoid new allocations.

Your profile also says that a lot of time is spent on maximize (55961 backtraces), so I would wonder how many iterations it’s spending in there. It looks like the Fortran toolbox also uses Brent method, and I suspect with an initial guess based on the previous iteration. So I’d suggest use the previous solution as initial guess to Julia maximize. Separately from that, if the maximize allocation is taking time, it would be trivial to remove those.

I still wonder if maximize allocations are costly. How much of the total 950 MB is done at the beginning, and how much in loops? It seems like the hot loops should not allocate at all.

Both the optimization and interpolation packages seem not to be intended for repeated runs on same-size but changing data. So unfortunately it may take mucking within those packages, or doing a custom solution, so make things faster. Interestingly, the Fortran book you cite seems to have a lightweight numerics library based on old Numerical Recipes code. Modern libraries tend to add lots of cool features, but there’s probably nothing like old Fortran to run fast with no allocations. Maybe there could be some lightweight Julia packages for that purpose.

2 Likes