IndexMap #65

JeffFessler · 2019-08-31T18:39:03Z

This is a generalization of the OffsetMap idea that will address #64 when complete.
I have some local tests that pass but it needs more testing, but I thought I'd get comments first.
Eventually this will support block diagonal linear maps and perhaps will simplify hvcat.

catch up to master

coveralls · 2019-08-31T18:51:15Z

Coverage decreased (-1.6%) to 93.105% when pulling d1dc812 on JeffFessler:master into 98cc86a on Jutho:master.

codecov · 2019-08-31T18:52:11Z

Codecov Report

Merging #65 into master will decrease coverage by 8.11%.
The diff coverage is 0%.

@@            Coverage Diff            @@
##           master     #65      +/-   ##
=========================================
- Coverage   98.62%   90.5%   -8.12%     
=========================================
  Files           8       9       +1     
  Lines         435     474      +39     
=========================================
  Hits          429     429              
- Misses          6      45      +39

Impacted Files	Coverage Δ
src/LinearMaps.jl	`100% <ø> (ø)`	⬆️
src/indexmap.jl	`0% <0%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 03960a9...eb1d1e3. Read the comment docs.

codecov · 2019-08-31T18:52:11Z

Codecov Report

Merging #65 into master will increase coverage by 0.12%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master      #65      +/-   ##
==========================================
+ Coverage   97.67%   97.79%   +0.12%     
==========================================
  Files           8        9       +1     
  Lines         473      545      +72     
==========================================
+ Hits          462      533      +71     
- Misses         11       12       +1

Impacted Files	Coverage Δ
src/LinearMaps.jl	`98.61% <ø> (-1.39%)`	⬇️
src/linearcombination.jl	`96.92% <100%> (ø)`	⬆️
src/indexmap.jl	`100% <100%> (ø)`
src/composition.jl	`98.38% <0%> (ø)`	⬆️
src/blockmap.jl	`98.58% <0%> (+0.03%)`	⬆️
src/uniformscalingmap.jl	`96.36% <0%> (+0.06%)`	⬆️
src/functionmap.jl	`97.22% <0%> (+0.07%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 98cc86a...d1dc812. Read the comment docs.

dkarrasch

Looks good, in principle. What I like is the implicit permutation of indices, that doesn't require writing a permutation factor in a CompositeMap. On the downside, I think we overlooked one major drawback of this approach in #64. The surrounding zero blocks don't come without cost. You have to set all entries to zero first, before you do any meaningful computation on the (potentially much) smaller block. Also, in the linear combination of IndexMaps, you would add entire vectors, as opposed to adding only the relevant part. This could be probably worked around partially by writing a specific 5-arg mul! method. That one, however, is not employed in the LinearCombination A_mul_B!, and there are non-trivial obstacles to switching to it. For those reasons, without significant rework of those other parts of the code, I don't think this LinearCombination-of-IndexMaps-approach can currently replace the BlockMap (and similarly an anticipated BlockDiagonalMap) implementation, unfortunately.

dkarrasch · 2019-09-02T11:48:13Z

src/indexmap.jl

+"""
+struct IndexMap{T, As <: LinearMap, Rs <: AbstractVector{Int},
+		Cs <: AbstractVector{Int}} <: LinearMap{T}
+	map::As


Any indention throughout should be 4 spaces, not aligned with stuff from the previous line.

hmm, then how does one make it clear to the reader that this is a continuation line rather than a new line of code?

But in this case, this is a new line of code, isn't it? You have the (split) line for the struct, and then the new line for its field map.

Oh, perhaps I misunderstood the previous comment.
Perhaps the issue is that I'm using a 4-space tab which looks ok in my editor but different in github where it defaults to an 8-space tab.

src/indexmap.jl

JeffFessler · 2019-09-03T18:16:15Z

I had/have the same concern about the zeros. I have some ideas but it will take time.

JeffFessler · 2019-09-03T18:51:00Z

My first attempt was to try to have only the first of a set of blocks zero out the output vector but that did not work for multiple reasons, including the fact that the addition order seems not to be preserved by the current implementation of LinearCombinationMap:
https://github.com/Jutho/LinearMaps.jl/blob/master/src/linearcombination.jl#L33
Now I really am going to put this on the back-burner.

dkarrasch · 2019-09-03T19:35:39Z

including the fact that the addition order seems not to be preserved by the current implementation of LinearCombinationMap:
https://github.com/Jutho/LinearMaps.jl/blob/master/src/linearcombination.jl#L33

This is just code saving, and calls the explicit method above, since order of addition is usually not expected to matter. If you need the final order of addition (whatever it comes out to be), you could define a special A_mul_B! or mul! function for A::LinearCombination{T,As} where As <: Tuple{Vararg{IndexMap}}, similarly to what I have done in the 5-arg mul-branch. Then, yes, you could first zero out the vector, and subsequently re-use the LinearCombination code. When you get back to this, you should take a look at the that branch, because I do something very similar to that already in the linear combination mul: the first map fills the y-vector via A_mul_B!, the other multiply and add via mul!(y, A.maps[i], x, true, true). In the latter, you would obviously only use the indices that you need (as you already do in the A_mul_B!).

JeffFessler · 2019-09-03T19:46:46Z

Thanks for the tips. I might wait for the dust to settle on the 5-arg mul-branch.

dkarrasch · 2019-09-03T19:49:54Z

Actually, you would only need to define a special 5-arg mul! function, and take care of that, in the branch alpha = beta = 1, so in the mul-add case, you only work on the specific indices. Otherwise, the LinearCombination code would automatically only zero out the output vector for the first map. For subsequent maps in the linear combination, it's the mul!(.., .., .., true, true) part operating. That also means, that you don't need the set0 field.

syncup

JeffFessler · 2019-09-30T02:28:16Z

Hi @dkarrasch, I'm back to working on IndexMap and now I have a version that I think handles the 5-arg mul! correctly and I've included an initial test set. It seems to be working fine for individual IndexMaps in v1.2. So now I want to be sure that LinearCombinations are efficient.
Here I need your input. It seems that the 5-arg mul! is now part of the master branch (yay!), but the code I see for LinearCombination multiplication seems not to be using the 5-arg mul! so I'm a bit lost. I'm hoping to use this for blockdiag so that efficiency is key.

dkarrasch · 2019-09-30T06:38:50Z

What do you mean? It's probably easy to overlook, but we now have this

https://github.com/Jutho/LinearMaps.jl/blob/95cc89cffff4b980af41aea32483af739a2d168c/src/linearcombination.jl#L86-L93

and this

https://github.com/Jutho/LinearMaps.jl/blob/95cc89cffff4b980af41aea32483af739a2d168c/src/linearcombination.jl#L102-L104

Note that 5-arg mul! is available only in Julia v1.3.

JeffFessler · 2019-09-30T11:00:51Z

Oops, I did overlook that because it is not executed in v1.2.
Four questions now:
Can you look at block_diag and see if you can see why it is failing @inferred?
Do we need to add IndexMap to FreeMap?
Is this line missing some A*x ingredients:
https://github.com/Jutho/LinearMaps.jl/blob/95cc89cffff4b980af41aea32483af739a2d168c/src/linearcombination.jl#L130
akin to what this line has?:
https://github.com/Jutho/LinearMaps.jl/blob/95cc89cffff4b980af41aea32483af739a2d168c/src/LinearMaps.jl#L49-L50
Any general comments about this PR? Especially about how to handle VERSION?

JeffFessler · 2019-10-08T12:36:00Z

Hi @dkarrasch, I know you have been busy with kron; when you are done with that will you be able to take a look at this PR?

dkarrasch · 2019-10-08T16:05:58Z

Yes, currently it's more about grading, proposal writing, paper reviewing and paper revisions, and less kron, but I'll certainly take a look when I find the time. 😂

dkarrasch

Okay, so here's a bunch of comments (sorry for the mess!). As for the linear combination, you could use the current FreeMap multiplication code:

function A_mul_B!(y::AbstractVector, A::LinearCombination{T,As}, x::AbstractVector) where {T, As<:Tuple{Vararg{IndexMap}}}
    # no size checking, will be done by individual maps
    A_mul_B!(y, A.maps[1], x) # this fills with zeros, and applies first map
    for n in 2:length(A.maps)
        mul!(y, A.maps[n], x, true, true) # these act only on the respective indices and adds the result
    end
    return y
end

And for the 5-arg mul!, a slight modification of the one I wrote in #72:

function mul!(y::AbstractVector, A::LinearCombination{T,As}, x::AbstractVector, α::Number=true, β::Number=false) where {T, As<:Tuple{Vararg{IndexMap}}}
    require_one_based_indexing(y, x)
    m, n = size(A)
    @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!"))
    iszero(β) && (A_mul_B!(y, A, x); rmul!(y, α); return y)
    if isone(β)
        fill!(y, zero(eltype(y)))
    else
        rmul!(y, β)
    end
    mapind = 0
    @views @inbounds for i in 1:length(A.maps)
        mul!(y, A.maps[i], x, α, true)
    end
    return y
end

At the end of the day, it's pretty much the exact same code running for BlockMaps or the BlockDiagonalMaps, except for the additional zero filling (and a few more generation of views, which is negligible).

dkarrasch · 2019-10-16T12:15:39Z

Manifest.toml

@@ -0,0 +1,19 @@
+# This file is machine-generated - editing it directly is not advised


This file is unnecessary and should be deleted.

dkarrasch · 2019-10-16T12:16:53Z

Project.toml

@@ -4,6 +4,7 @@ version = "2.5.2"

 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"


I'm not sure what exactly you need to set the random seed for. Are your tests sensitive to the actual random matrix that is generated?

dkarrasch · 2019-10-16T12:18:22Z

src/indexmap.jl

+
+
+"""
+`check_index(row::AbstractVector{Int}, dims::Dims{2})`


If this is indented by four spaces, then it is typeset as code even without the ` marks.

src/indexmap.jl

dkarrasch · 2019-10-16T15:52:07Z

src/indexmap.jl

+Base.:(==)(A::IndexMap, B::IndexMap) = (eltype(A) == eltype(B)) &&
+    (A.map == B.map) && (A.dims == B.dims) &&
+    (A.rows == B.rows) && (A.cols == B.cols)


I'm not sure, but I believe this would be the fallback comparison for custom structs anyway, like comparing all parameter types and all fields. If this is so, then one could remove this definition.

dkarrasch · 2019-10-16T16:45:59Z

test/block_diag.jl

+using LinearMaps: LinearMap
+using SparseArrays: sparse, blockdiag
+using Random: seed!
+using Test: @test, @testset


As for the tests, I'd recommend to follow the style in the other test files. I don't think there is need to use every single command. Most likely, users will simply load these packages as a whole.

dkarrasch · 2019-10-16T16:46:57Z

test/block_diag.jl

+    x = randn(size(Md,2))
+#   Bd = @inferred block_diag(L1, L2, L3, L2, L1) # todo: type instability
+    Bd = block_diag(L1, L2, L3, L2, L1)
+    @test Bd isa LinearMaps.LinearCombination


I don't think we would want to hardcode this type. Whatever the type of Bd is, it should yield the correct y = A*x, and this is what we test with, e.g., Matrix(Bd) etc.

dkarrasch · 2019-10-16T16:48:10Z

test/indexmap.jl

+    @test Matrix(Bh') == Matrix(Bh)'
+    @test Mh * x ≈ Bh * x
+
+end


All of these tests look like fairly generic tests. I don't think it's worth it to set the random seed, and therefore add another dependency, even if it's a stdlib.

dkarrasch · 2019-10-16T16:49:58Z

test/block_diag.jl

+    # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse:
+    Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
+    x = randn(size(Md,2))
+#   Bd = @inferred block_diag(L1, L2, L3, L2, L1) # todo: type instability


One would have to take a closer look to identify the cause for type instability. This must be somewhere in the construction of IndexMap, since this comes first and seems to be type-unstable, according to your comment in the other test file.

dkarrasch · 2019-10-17T07:29:43Z

src/indexmap.jl

+            fill!(y_in, zero(eltype(y_in)))
+            A_mul_B!(yv, A, xv)
+        elseif isone(β) # β = 1
+            yv .+= A * xv


This should be mul!(yv, A, xv, true, true), otherwise the right hand side will allocate, even if A is matrix.

JeffFessler and others added 5 commits August 24, 2019 17:21

Fix index in BlockMap adjoint/transpose multiply

19161f9

Add hvcat test for "misaligned" blockmap

65738fa

minor code style changes

1f2501b

Merge pull request #1 from Jutho/master

15a5934

catch up to master

Initial version of IndexMap

eb1d1e3

dkarrasch reviewed Sep 2, 2019

View reviewed changes

Add block_diag

a87f231

JeffFessler and others added 5 commits September 26, 2019 08:37

Merge pull request #2 from Jutho/master

98bdfa2

syncup

Flag possible errors in mul! for LinearCombination

6395500

Add indexmap tests

4bdc059

Add Random for tests

f2e6e9a

Update for 5-arg mul!

7e8b447

JeffFessler added 2 commits September 29, 2019 23:06

Rename block_diag

905ab6d

Update tests

9cd7857

JeffFessler added 3 commits September 30, 2019 06:16

Tab to 4-space

5ff80ff

Try to fix type instability hcat_new block_diag

fd8a30d

Update tests

d767d07

Test complex case

d1dc812

JeffFessler mentioned this pull request Oct 6, 2019

add Kronecker products for LinearMaps #61

Merged

dkarrasch reviewed Oct 17, 2019

View reviewed changes

dkarrasch mentioned this pull request Oct 22, 2019

BlockMap enhancements #72

Merged

dkarrasch mentioned this pull request Mar 26, 2020

Only store non-zero blocks in BlockMap #84

Draft

JeffFessler mentioned this pull request Jun 8, 2020

WIP: get rid of A*_mul_B! #74

Closed

dkarrasch mentioned this pull request Jan 27, 2021

RestrictionMap/SamplingMap? #134

Closed

dkarrasch mentioned this pull request Apr 4, 2022

Add EmbeddedMap #174

Merged

dkarrasch closed this in #174 May 19, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

IndexMap #65

IndexMap #65

JeffFessler commented Aug 31, 2019

coveralls commented Aug 31, 2019 •

edited

Loading

codecov bot commented Aug 31, 2019

codecov bot commented Aug 31, 2019 •

edited

Loading

dkarrasch left a comment

dkarrasch Sep 2, 2019

JeffFessler Sep 3, 2019

dkarrasch Sep 3, 2019

JeffFessler Sep 3, 2019

JeffFessler commented Sep 3, 2019

JeffFessler commented Sep 3, 2019

dkarrasch commented Sep 3, 2019

JeffFessler commented Sep 3, 2019

dkarrasch commented Sep 3, 2019

JeffFessler commented Sep 30, 2019

dkarrasch commented Sep 30, 2019

JeffFessler commented Sep 30, 2019

JeffFessler commented Oct 8, 2019

dkarrasch commented Oct 8, 2019

dkarrasch left a comment

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 16, 2019

dkarrasch Oct 17, 2019

		@@ -0,0 +1,19 @@
		# This file is machine-generated - editing it directly is not advised



		"""
		`check_index(row::AbstractVector{Int}, dims::Dims{2})`

IndexMap #65

IndexMap #65

Conversation

JeffFessler commented Aug 31, 2019

coveralls commented Aug 31, 2019 • edited Loading

codecov bot commented Aug 31, 2019

Codecov Report

codecov bot commented Aug 31, 2019 • edited Loading

Codecov Report

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

JeffFessler commented Sep 3, 2019

JeffFessler commented Sep 3, 2019

dkarrasch commented Sep 3, 2019

JeffFessler commented Sep 3, 2019

dkarrasch commented Sep 3, 2019

JeffFessler commented Sep 30, 2019

dkarrasch commented Sep 30, 2019

JeffFessler commented Sep 30, 2019

JeffFessler commented Oct 8, 2019

dkarrasch commented Oct 8, 2019

dkarrasch left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Aug 31, 2019 •

edited

Loading

codecov bot commented Aug 31, 2019 •

edited

Loading