matrix_gradient(): refactor

* rename from get_matrix_derivative(): it's not a getter, and gradient is a better term
* construct sparse matrix directly, which is much more efficient
* parameters arg is not needed

Author: Alexey Stukalov

src/additional_functions/parameters.jl

@@ -86,15 +86,19 @@ function check_constants(M)
     return false
 end
 
-function get_matrix_derivative(M_indices, parameters, n_long)
-    ∇M = [
-        sparsevec(M_indices[i], ones(length(M_indices[i])), n_long) for
-        i in 1:length(parameters)
-    ]
-
-    ∇M = reduce(hcat, ∇M)
-
-    return ∇M
+# construct length(M)×length(parameters) sparse matrix of 1s at the positions,
+# where the corresponding parameter occurs in the M matrix
+function matrix_gradient(M_indices::ArrayParamsMap, M_length::Integer)
+    rowval = reduce(vcat, M_indices)
+    colptr =
+        pushfirst!(accumulate((ptr, M_ind) -> ptr + length(M_ind), M_indices, init = 1), 1)
+    return SparseMatrixCSC(
+        M_length,
+        length(M_indices),
+        colptr,
+        rowval,
+        ones(length(rowval)),
+    )
 end
 
 # fill M with parameters

src/imply/RAM/generic.jl

@@ -158,8 +158,8 @@ function RAM(;
     I_A = similar(A_pre)
 
     if gradient
-        ∇A = get_matrix_derivative(A_indices, parameters, n_nod^2)
-        ∇S = get_matrix_derivative(S_indices, parameters, n_nod^2)
+        ∇A = matrix_gradient(A_indices, n_nod^2)
+        ∇S = matrix_gradient(S_indices, n_nod^2)
     else
         ∇A = nothing
         ∇S = nothing
@@ -168,15 +168,8 @@ function RAM(;
     # μ
     if meanstructure
         has_meanstructure = Val(true)
-
-        if gradient
-            ∇M = get_matrix_derivative(M_indices, parameters, n_nod)
-        else
-            ∇M = nothing
-        end
-
+        ∇M = gradient ? matrix_gradient(M_indices, n_nod) : nothing
         μ = zeros(n_var)
-
     else
         has_meanstructure = Val(false)
         M_indices = nothing