Fix weighted normal diag op

Author: nHackel

src/DiagOp.jl

@@ -148,8 +148,8 @@ function LinearOperatorCollection.normalOperator(diag::DiagOp, W=opEye(eltype(di
 
     # we promote the weights to be of the same type as T, which will be required
     # when creating the temporary vector in normalOperator in a later stage
-    opInner = normalOperator(diag.ops[1], WeightingOp(T; weights=T.(weights[diag.yIdx[1]:diag.yIdx[2]-1])); copyOpsFn = copyOpsFn, kwargs...)
-    op = DiagNormalOp([copyOpsFn(opInner) for i=1:length(diag.ops)], size(diag,2), diag.xIdx, S(undef, diag.ncol) )
+    op = DiagNormalOp([normalOperator(copyOpsFn(diag.ops[1]), WeightingOp(T; weights=T.(weights[diag.yIdx[i]:diag.yIdx[i+1]-1])); copyOpsFn = copyOpsFn, kwargs...)
+                     for i in 1:length(diag.ops)], size(diag,2), diag.xIdx, S(undef, diag.ncol) )
   else
     op = DiagNormalOp([normalOperator(diag.ops[i], WeightingOp(T; weights=T.(weights[diag.yIdx[i]:diag.yIdx[i+1]-1])); copyOpsFn = copyOpsFn, kwargs...)
                      for i in 1:length(diag.ops)], size(diag,2), diag.xIdx, S(undef, diag.ncol) )

test/testOperators.jl

@@ -363,6 +363,23 @@ function testDiagOp(N=32,K=2;arrayType = Array)
     @test y2 ≈ y3 rtol = 1e-2
   end
 
+  @testset "Weighted Diag Normal" begin
+    w = rand(eltype(op1), size(op1, 1))
+    wop = WeightingOp(w)
+    prod1 = ProdOp(wop, op1)
+    prod2 = ProdOp(wop, op2)
+    prod3 = ProdOp(wop, op3)
+
+    y = Array(adjoint(F) * adjoint(wop) * wop * F* x)
+    y1 = Array(normalOperator(prod1) * x)
+    y2 = Array(normalOperator(prod2) * x)
+    y3 = Array(normalOperator(prod3) * x)
+
+    @test y ≈ y1 rtol = 1e-2
+    @test y1 ≈ y2 rtol = 1e-2
+    @test y2 ≈ y3 rtol = 1e-2
+  end
+
 end
 
 function testRadonOp(N=32;arrayType = Array)