added a plot showing progression of a sample with the iteratons of CG

demos/krylov_sampling.py

@@ -32,8 +32,10 @@ R_p = R * FFT * A
 j = R_p.adjoint(N.inverse(d))
 D_inv = R_p.adjoint * N.inverse * R_p + S.inverse
 
+history = []
 
 def sample(D_inv, S, j,  N_samps, N_iter):
+    global history
     space = D_inv.domain
     x = ift.Field.zeros(space)
     r = j.copy()
@@ -43,14 +45,16 @@ def sample(D_inv, S, j,  N_samps, N_iter):
     for i in range(N_samps):
         y += [S.draw_sample()]
     for k in range(1, 1 + N_iter):
+        history += [y[0].copy()]
         gamma = r.vdot(r) / d
         if gamma == 0.:
             break
         x += gamma * p
-        print(p.vdot(D_inv(j)))
+        #print(p.vdot(D_inv(j)))
         for i in range(N_samps):
             y[i] -= p.vdot(D_inv(y[i])) * p / d
             y[i] += np.random.randn() / np.sqrt(d) * p
+        print("variance iteration "+str(k)+":", np.sqrt(p.vdot(p)/d))
         #r_new = j - D_inv(x)
         r_new = r - gamma * D_inv(p)
         beta = r_new.vdot(r_new) / (r.vdot(r))
@@ -106,3 +110,10 @@ plt.plot((s_x - m_x).val, color='k', label='signal - mean')
 plt.legend()
 plt.savefig('Krylov_uncertainty.png')
 plt.close()
+
+for i in range(min(6, len(history))):
+    plt.plot(sky(history[i]).val, label="step " + str(i+1))
+plt.plot(s_x.val-m_x.val, 'k-', label="residual")
+plt.legend()
+plt.savefig('iterations.png')
+plt.close()