add some debugging

src/DualTVDD.jl

@@ -47,11 +47,14 @@ function rundd()
     f = zeros(2,2)
     f[1,:] .= 1
     #g = [0. 2; 1 0.]
-    A = diagm(vcat(fill(2, length(f)÷2), fill(1, length(f)÷2)))
-    A = rand(length(f), length(f))
+    A = diagm(vcat(fill(1, length(f)÷2), fill(1/1000, length(f)÷2)))
+    #A = rand(length(f), length(f))
     display(A)
+    println(cond(A))
+    display(eigen(A))
     #A = diagm(fill(1/2, length(f)))
     B = inv(A'*A + β*I)
+    println(norm(sqrt(B)))
 
     #println(norm(sqrt(B)))
 
@@ -61,7 +64,7 @@ function rundd()
     α = .25
 
     md = DualTVDD.DualTVDDModel(f, A, α, 0., 0.)
-    alg = DualTVDD.DualTVDDAlgorithm(M=(1,1), overlap=(1,1), σ=0.25)
+    alg = DualTVDD.DualTVDDAlgorithm(M=(1,1), overlap=(1,1), σ=1)
     ctx = DualTVDD.init(md, alg)
 
     md2 = DualTVDD.OpROFModel(g, B, α)
@@ -69,10 +72,10 @@ function rundd()
     ctx2 = DualTVDD.init(md2, alg2)
 
 
-    for i in 1:1000
+    for i in 1:1
         step!(ctx)
     end
-    for i in 1:10000
+    for i in 1:100000
         step!(ctx2)
     end
 
@@ -99,6 +102,9 @@ function rundd()
     display(recover_u!(ctx))
     display(recover_u!(ctx2))
 
+    println(energy(ctx))
+    println(energy(ctx2))
+
     ctx, ctx2
 end
 

src/dualtvdd.jl

@@ -109,11 +109,11 @@ function step!(ctx::DualTVDDContext)
         subIB = I - ctx.B[vec(li[sax...]), vec(li[sax...])]./λ
         subB = ctx.B[vec(li[sax...]), vec(li[sax...])]./λ
 
-        for j in 1:50
+        for j in 1:1000000
             subΛp = map(kΛ, ctx.subctx[i].p)
             vec(ctx.subg[i]) .= subIB * vec(subΛp) .+ subB * vec(gloc)
 
-            for k in 1:10
+            for k in 1:10000
                 step!(ctx.subctx[i])
             end
         end