fix to jacobian HCAT

Author: nantonel

src/calculus/Jacobian.jl

@@ -32,23 +32,26 @@ Jacobian(S::T2,x::AbstractArray) where {T,L,T2<:Scale{T,L}} = Scale(S.coeff,Jaco
 function Jacobian(H::DCAT{N,L,P1,P2},x) where {N,L,P1,P2} 
 	A = ()
 	c = 0
-	for i = 1:N
-		A = length(H.idxD[i]) == 1 ?
-		(A...,jacobian(H.A[i],x[c+1])) :
-		(A...,jacobian(H.A[i],x[c+1:c+length(H.idxD[i])])) 
-		c += length(H.idxD[i])
+    for (k, idx) in enumerate(H.idxD)
+		if length(idx) == 1 
+            A = (A...,jacobian(H.A[k],x[idx])) 
+        else
+            xx = ([x[i] for i in idx]...)
+            A = (A...,jacobian(H.A[k],xx)) 
+        end
 	end
 	DCAT(A,H.idxD,H.idxC)
 end
 #Jacobian of HCAT
 function Jacobian(H::HCAT{M,N,L,P,C},x::D) where {M,N,L,P,C,D} 
 	A = ()
-	c = 0
-	for i = 1:N
-		A = length(H.idxs[i]) == 1 ?
-		(A...,jacobian(H.A[i],x[c+1])) :
-		(A...,jacobian(H.A[i],x[c+1:c+length(H.idxs[i])])) 
-		c += length(H.idxs[i])
+    for (k, idx) in enumerate(H.idxs)
+		if length(idx) == 1 
+            A = (A...,jacobian(H.A[k],x[idx])) 
+        else
+            xx = ([x[i] for i in idx]...)
+            A = (A...,jacobian(H.A[k],xx)) 
+        end
 	end
 	HCAT(A,H.idxs,H.buf,M)
 end

src/syntax.jl

@@ -35,7 +35,7 @@ function getindex(A::AbstractOperator,idx...)
 	elseif length(idx) == 1  && ndoms(A,2) == length(idx[1])
 		return permute(A,idx[1])
 	else
-		error("cannot split operator of type $(typeof(H.A[i]))")
+		error("cannot split operator of type $(typeof(A))")
 	end
 end
 

test/test_linear_operators_calculus.jl

@@ -190,9 +190,9 @@ y1 = remove_displacement(opD)*(x1,x2,x3)
 y2 = (A1*x1, A2*x2, A3*x3)
 @test all(vecnorm.(y1 .- y2) .<= 1e-12)
 
-#######################
-## test HCAT    #######
-#######################
+########################
+### test HCAT    #######
+########################
 
 m, n1, n2 = 4, 7, 5
 A1 = randn(m, n1)
@@ -206,6 +206,12 @@ y1 = test_op(opH, (x1, x2), randn(m), verb)
 y2 = A1*x1 + A2*x2
 @test vecnorm(y1-y2) <= 1e-12
 
+#permuatation 
+p = [2;1]
+opHp = opH[p]
+y1 = test_op(opHp, (x2, x1), randn(m), verb)
+@test vecnorm(y1-y2) <= 1e-12
+
 # test HCAT longer
 
 m, n1, n2, n3 = 4, 7, 5, 6

test/test_nonlinear_operators_calculus.jl

@@ -11,6 +11,17 @@ y, grad = test_NLop(op,x,r,verb)
 Y = 30*(A*x)
 @test vecnorm(Y-y) <1e-8
 
+m = 3
+x = randn(m)
+r = randn(m)
+A = Pow(Float64,(m,),2)
+op = -A
+
+y, grad = test_NLop(op,x,r,verb)
+
+Y = -A*x
+@test vecnorm(Y-y) <1e-8
+
 #testing DCAT
 n,m = 4,3
 x = (randn(n),randn(m))
@@ -38,6 +49,25 @@ y, grad = test_NLop(op,x,r,verb)
 Y = A*x[1]+B*x[2]
 @test vecnorm(Y-y) <1e-8
 
+m,n = 3,5
+x = (randn(m),randn(n))
+r = randn(m)
+A= Sin(Float64,(m,))
+M = randn(m,n)
+B = MatrixOp(M)
+op = HCAT(A,B)
+
+y, grad = test_NLop(op,x,r,verb)
+
+Y = A*x[1]+M*x[2]
+@test vecnorm(Y-y) <1e-8
+
+p = [2,1]
+opP = permute(op,p)
+J = Jacobian(opP,x[p])'
+println(size(J,1))
+y, grad = test_NLop(opP,x[p],r,verb)
+
 #testing VCAT
 n,m = 4,3
 x = randn(m)