I am fitting a statistical model and I refactored the code to reduce the allocations it generated. I morally halved them, with respect to the first version, but there are still some unclear sections which I see from running my code with the trackallocation=user
option, where I cant understand why they allocate.
The sections involved are the following. In the fitting functon:
 # lC = @MVector zeros(2) # prellocated before this loop
0 for k in 1:nclus_red
 # filter the covariates of the units of label k
0 aux_idxs = findall_faster(jj > Si_red[jj] == k, 1:n_red)
0 if sPPM
0 copy!(s1o, sp1_red[aux_idxs])
0 copy!(s2o, sp2_red[aux_idxs])
0 copy!(s1n,s1o); push!(s1n, sp1[j])
0 copy!(s2n,s2o); push!(s2n, sp2[j])
1012512 spatial_cohesion!(spatial_cohesion_idx, s1o, s2o, sp_params_real, true, M_dp, S,1,false,lC)
1012512 spatial_cohesion!(spatial_cohesion_idx, s1n, s2n, sp_params_real, true, M_dp, S,2,false,lC)
 end
In the utils file:
 function spatial_cohesion!(idx::Real, s1::AbstractVector{Float64}, s2::AbstractVector{Float64}, sp_params::Vector, lg::Bool, M::Real, S=@MMatrix(zeros(2, 2)), case::Int=1, add::Bool=false, lC=@MVector(zeros(2)))
0 if idx==1.0 cohesion1!(s1,s2,sp_params[1],lg,M,case,add,lC); return; end
0 if idx==2.0 cohesion2!(s1,s2,sp_params[1],lg,M,case,add,lC); return; end
4525664 if idx==3.0 cohesion3!(s1,s2,sp_params[1],sp_params[2],sp_params[3],sp_params[4],lg,M,S,case,add,lC); return; end
0 if idx==4.0 cohesion4!(s1,s2,sp_params[1],sp_params[2],sp_params[3],sp_params[4],lg,M,S,case,add,lC); return; end
0 if idx==5.0 cohesion5!(s1,s2,sp_params[1],lg,M,case,add,lC); return; end
0 if idx==6.0 cohesion6!(s1,s2,sp_params[1],lg,M,case,add,lC); return; end
 end
 function cohesion3!(s1::AbstractVector{Float64}, s2::AbstractVector{Float64}, mu_0::AbstractVector{Float64}, k0::Real, v0::Real, Psi::AbstractMatrix{Float64}, lg::Bool, M::Real=1.0, S=@MMatrix(zeros(2, 2)), case::Int=1, add::Bool=false, lC=@MVector(zeros(2)))::Float64
0 sdim = length(s1)
 # Compute sample means
0 sbar1 = mean(s1)
0 sbar2 = mean(s2)
 # Compute deviations from the sample mean
0 S .= 0.
0 @inbounds for i in 1:sdim
0 s_sbar1 = s1[i]  sbar1
0 s_sbar2 = s2[i]  sbar2

0 S[1, 1] += s_sbar1 * s_sbar1
0 S[2, 2] += s_sbar2 * s_sbar2
0 S[2, 1] += s_sbar1 * s_sbar2
0 end
0 S[1, 2] = S[2, 1] # to avoid repeating computations
 # Updated parameters for cohesion 3
0 kn = k0 + sdim
0 vn = v0 + sdim

 sbar = SVector((sbar1, sbar2))
0 auxvec1 = sbar . mu_0
0 auxmat1 = auxvec1 * auxvec1'

0 auxconst1 = k0 * sdim
0 auxconst2 = k0 + sdim
0 Psi_n = Psi .+ S .+ auxconst1 / (auxconst2) .* auxmat1

0 out = sdim * logpi + G2a(0.5 * vn, true)  G2a(0.5 * v0, true) + 0.5 * v0 * logdet(Psi)  0.5 * vn * logdet(Psi_n) + log(k0)  log(kn)
0 if add
0 lC[case] += lg ? out : exp(out)
 else
0 lC[case] = lg ? out : exp(out)
 end
 end
 function G2a(a::Real, lg::Bool)
 out = logpi + lgamma(a) + lgamma(a  0.5)
 return lg ? out : exp(out)
 end
I dont provide a minimal working example since itâs a bit complex to simulate and since I believe/hope itâs just a theoretical/syntactical detail that I missed (e.g. maybe the overhead of the intermediate function call and not directly cohesion3!
). Because working with inplace mutating functions I expected to have no memory issues.