How to output selective elements from Matrix?

Amro · August 1, 2022, 6:20pm

I have the below configuration. I am trying to output selective elements from matrix L.sl based on values of L.Seg. More specifically, I want to write syntax such as L.sl[L.Seg] to give me the same results as S. Any simple way to do that?

using SparseArrays, LinearAlgebra

Base.@kwdef mutable struct Ls
   tol::Vector{Float64} = [0,0,0]
   Seg::Vector{Int8} = [1,2,3]
   sl = []
end


L = Ls[];
push!(L, Ls());
L[1].sl = [1 2 3; 4 5 6; 7 8 9];
push!(L, Ls());
L[2].sl = [1 2 3; 4 5 6; 7 8 9];

@eval L = Ls((reduce(vcat, getfield.(L, pn)) for pn in fieldnames(Ls))...);

S = zeros(length(L.tol));
for j in eachindex(L.tol)
   S[j] = L.sl[j,L.Seg[j]];
end
#########################################

julia> L.sl
6×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9
 1  2  3
 4  5  6
 7  8  9

julia> S
6-element Vector{Float64}:
 1.0
 5.0
 9.0
 1.0
 5.0
 9.0

julia> L.sl[L.Seg]
6-element Vector{Int64}:
 1
 4
 7
 1
 4
 7

Amro · August 2, 2022, 1:15pm

Any comment here please?

mikmoore · August 2, 2022, 2:15pm

I’m not sure I completely understand what you’re after. But it sounds like you might be looking for CartesianIndex? For example,

inds = CartesianIndex.(eachindex(L.Seg), L.Seg)
S = L.sl[inds]

You could also accomplish this with broadcasted getindex, although it’s a bit less typical and can be trickier to understand

S = getindex.(Ref(L.sl), eachindex(L.Seg), L.Seg)

Amro · August 2, 2022, 3:13pm

Thank you very much. Exactly this is what I need.
I have another issue please. Below is an MWE, in which I am trying to update the sparse matrix S. However, non of the two methods are working. Could you please suggest a solution?

using SparseArrays, LinearAlgebra

function dump(R)
  R.nzval .= [7,8,9];
end
S = sparse([1 0 0 0 0; 0 0 0 0 0; 0 0 1 0 0; 0 0 0 2 0; 0 0 0 0 5])
# Method-1
dump(S[3:end,3:end])
5×5 SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  2  ⋅
 ⋅  ⋅  ⋅  ⋅  5
# Method-2
dump( @view( S[3:end,3:end]) )
ERROR: type SubArray has no field nzval

mikmoore · August 2, 2022, 3:54pm

You should think very hard before mucking with the internal fields of objects, like you’re doing with R.nzval .= [7,8,9]. It’s generally a bad idea that can break things or at least be difficult to reason about and maintain.

In the dump(S[3:end,3:end]) version, S[3:end,3:end] copies the indexed part of the array, so changing it’s nzval doesn’t do anything to the original S. You need to return R if you want this modified array to be accessible outside of dump (although note that the first two rows/columns of S were never a part of R and will not be returned).

The @views version puts a wrapper object around the original array. The wrapper does not have a nzvals field to access, hence your error.

Amro · August 2, 2022, 4:08pm

Thank you very much. It is clear now.
For performance purposes, what is the best way to mutate S, as below?

function dump(R)
  R.nzval .= [7,8,9];
return R
end
S = sparse([1 0 0 0 0; 0 0 0 0 0; 0 0 1 0 0; 0 0 0 2 0; 0 0 0 0 5])
# Method-1
S[3:end,3:end] .= dump(S[3:end,3:end])

mikmoore · August 2, 2022, 5:04pm

Nothing wrong with loops in Julia:

using SparseArrays
k = 0
newvals = [7,8,9]
ii,jj,vv = findnz(S)
for (i,j) in zip(ii,jj)
  if i >= 3 && j >= 3
    S[i,j] = newvals[k+=1]
  end
end
S

If you’re just trying to mutate the diagonal of the matrix:

newvals = [7,8,9]
for (i,v) in zip(3:5,newvals)
  S[i,i] = v
end

Amro · August 2, 2022, 5:24pm

Thank you for your note.
Actually, my main issue is how to integrate this mutation inside the function dump (which exist in my main code)?

mikmoore · August 2, 2022, 5:55pm

using SparseArrays
function dump(R,startat)
	k = 0
	newvals = 7:9
	ii,jj,vv = findnz(R)
	for (i,j) in zip(ii,jj)
		if i >= startat && j >= startat
			R[i,j] = newvals[k+=1]
		end
	end
	return R
end

# S will be mutated since dump mutates the first argument
dump(S,3)

Since dump mutates its argument, we’d typically name it dump! instead.

Amro · August 2, 2022, 6:45pm

Perfect. There are two things:
1- Since dump(S,3) is not assigned to any parameters, so what will return R do?
2- I deleted the return R and the results are not effected.

mikmoore · August 2, 2022, 7:03pm

The return R allows you to write X = dump(S,3) and now get X === S. Often times this can be useful, but it’s certainly not necessary. Without this, dump returns the value of it’s last line. Since the last “line” here is a loop body, it will return nothing if no explicit return is written.

Amro · August 2, 2022, 7:29pm

Thank you very much for your details!