playing around, definetly working!

paper/1d_separation.py

@@ -15,7 +15,6 @@ if __name__ == '__main__':
     h_space = s_space.get_default_codomain()
     FFT = ift.FFTOperator(h_space)
     p_spec = lambda k: (1./(1+k)**4)
-    mod_p_spec = lambda k: p_spec(k*1024)*1024**4
     binbounds = ift.PowerSpace.useful_binbounds(h_space,logarithmic=False)#, nbin=100)
     p_space = ift.PowerSpace(h_space,binbounds=binbounds)
     k_lengths=p_space.k_lengths
@@ -39,26 +38,30 @@ if __name__ == '__main__':
     sextracted8 = fits.open('check8.fits')
     sextracted8 = sextracted8[0].data
     sextracted8 = sextracted8.clip(0.001)
+    dcnn_diffuse = np.load('tf_diffuse.npy')
+    dcnn_points = np.load('tf_points.npy')
+    dcnn_diffuse = dcnn_diffuse.clip(0.001)
+    dcnn_points = dcnn_points.clip(0.001)
 
-    energy1 = sb.build_starblade(data,1.0, newton_iterations=200, cg_iterations=3,
+    energy1 = sb.build_starblade(data,1.0, newton_iterations=200, cg_iterations=10,
                                  q=q)#, manual_power_spectrum= p_spec)
-    energy2 = sb.build_starblade(data,1.5,  newton_iterations=200, cg_iterations=3,
+    energy2 = sb.build_starblade(data,1.5,  newton_iterations=200, cg_iterations=10,
                                  q=q)#, manual_power_spectrum= p_spec)
-    energy3 = sb.build_starblade(data,3., newton_iterations=200, cg_iterations=3,
+    energy3 = sb.build_starblade(data,3., newton_iterations=200, cg_iterations=10,
                                  q=q)#, manual_power_spectrum= p_spec)
     # ift.extra.check_value_gradient_consistency(energy1, tol=1e-3)
-    for i in range(10):
-        energy1 = sb.starblade_iteration(energy1, samples=5)
-        energy2 = sb.starblade_iteration(energy2, samples=5)
-        energy3 = sb.starblade_iteration(energy3, samples=5)
-        print "error energy1:", np.sqrt(((ift.exp(energy1.s).val - ift.exp(s).val) ** 2).mean())
-        print "error energy2:", np.sqrt(((ift.exp(energy2.s).val - ift.exp(s).val) ** 2).mean())
-        print "error energy3:", np.sqrt(((ift.exp(energy3.s).val - ift.exp(s).val) ** 2).mean())
+    for i in range(5):
+        energy1 = sb.starblade_iteration(energy1, samples=3)
+        energy2 = sb.starblade_iteration(energy2, samples=3)
+        energy3 = sb.starblade_iteration(energy3, samples=3)
+        print "error energy1:", np.sqrt(((1 - ift.exp(s).val/ift.exp(energy1.s).val) ** 2).mean())
+        print "error energy2:", np.sqrt(((1 - ift.exp(s).val/ift.exp(energy2.s).val) ** 2).mean())
+        print "error energy3:", np.sqrt(((1 - ift.exp(s).val/ift.exp(energy3.s).val) ** 2).mean())
 
     samples = []
     n=30
     for i in range(n):
-        samples.append(energy1.curvature.inverse.draw_sample())
+        samples.append(energy2.curvature.inverse.draw_sample())
     m = 0
     v = 0
     s_s=0
@@ -66,7 +69,7 @@ if __name__ == '__main__':
     pos_tanh = ift.library.nonlinearities.PositiveTanh()
     p=0
     for sample in samples:
-        a_s = pos_tanh(sample)
+        a_s = pos_tanh(energy2.position+sample)
         m +=  a_s
         v += a_s**2
         s_s += ift.log(d*(1-a_s))
@@ -94,7 +97,7 @@ if __name__ == '__main__':
     plt.savefig('2d_data.pdf')
     # plt.figure(figsize=(8, 8))
     # plt.plot(data, 'k-')
-    f, (ax0, ax1,ax2,ax3,ax4) = plt.subplots(5, sharex=True, sharey=True,figsize=(8, 10))
+    f, (ax0, ax1,ax2,ax3,ax4,ax5) = plt.subplots(6, sharex=True, sharey=True,figsize=(8, 12))
     plt.suptitle('diffuse components', size=size)
 
     for i in range(energy1.s.val.shape[0]):
@@ -119,12 +122,17 @@ if __name__ == '__main__':
     for i in range(energy1.s.val.shape[0]):
         ax4.plot((sextracted8)[i].clip(0.0001), 'k-', alpha=(0.15 / (energy1.s.val.shape[0]) * i))
     ax4.yaxis.set_label_position("right")
-    ax4.set_ylabel(r'sextractor'+'\n' + r'\textsf{BACK\_SIZE}'+r'$=8$', size=size)
+    ax5.set_ylabel(r'sextractor'+'\n' + r'\textsf{BACK\_SIZE}'+r'$=8$', size=size)
+    for i in range(energy1.s.val.shape[0]):
+        ax5.plot((dcnn_diffuse)[i].clip(0.0001), 'k-', alpha=(0.15 / (energy1.s.val.shape[0]) * i))
+    ax5.yaxis.set_label_position("right")
+
+    ax5.set_ylabel(r'DCAE', size=size)
 
     plt.savefig('1d_diffuse.pdf')
 
     plt.figure()
-    f, (ax0, ax1,ax2,ax3,ax4) = plt.subplots(5, sharex=True, sharey=True,figsize=(8, 10))
+    f, (ax0, ax1,ax2,ax3,ax4,ax5) = plt.subplots(6, sharex=True, sharey=True,figsize=(8, 12))
 
     plt.suptitle('point-like components', size=size)
     for i in range(energy1.s.val.shape[0]):
@@ -149,6 +157,10 @@ if __name__ == '__main__':
         ax4.plot((data-sextracted8)[i].clip(0.0001), 'k-', alpha=(0.15 / (energy1.s.val.shape[0]) * i))
     ax4.yaxis.set_label_position("right")
     ax4.set_ylabel(r'sextractor'+'\n' + r'\textsf{BACK\_SIZE}'+r'$=8$', size=size)
+    for i in range(energy1.s.val.shape[0]):
+        ax5.plot((dcnn_points)[i].clip(0.0001), 'k-', alpha=(0.15 / (energy1.s.val.shape[0]) * i))
+    ax5.set_ylabel(r'DCAE', size=size)
+    ax5.yaxis.set_label_position("right")
 
     ax0.set_yscale("log")
 
@@ -157,7 +169,7 @@ if __name__ == '__main__':
     # plt.ylim(1e-0)
     plt.savefig('1d_points.pdf')
     plt.figure()
-    f, (ax0, ax1,ax2) = plt.subplots(3, sharex=True, sharey=True)
+    f, (ax0, ax1,ax2) = plt.subplots(3, sharex=True, sharey=True,figsize=(8, 6))
     plt.suptitle('data and true components', size=size)
     for i in range(energy1.s.val.shape[0]):
         ax0.plot(data[i], 'k-', alpha=(0.15 / (energy1.s.val.shape[0]) * i))
@@ -207,10 +219,23 @@ if __name__ == '__main__':
     plt.ylabel('power',size=15)
     plt.xlabel('harmonic mode',size=15)
     plt.savefig('1d_power.pdf')
+    plt.figure()
+    plt.scatter(s.val, np.log(dcnn_diffuse), alpha = 0.2, c='k', label='DDCAE')
+    plt.scatter(s.val, energy2.s.val, alpha = 0.2, c='r', label = 'starblade')
+    plt.title('diffuse truth vs DDCAE and starblade')
+    plt.legend()
+    plt.savefig('scatter.pdf')
     plt.close('all')
-    print "error energy1:", np.sqrt(((ift.exp(energy1.s).val-ift.exp(s).val)**2).mean())
-    print "error energy2:", np.sqrt(((ift.exp(energy2.s).val-ift.exp(s).val)**2).mean())
-    print "error energy3:", np.sqrt(((ift.exp(energy3.s).val-ift.exp(s).val)**2).mean())
-    print "error back_size8:", np.sqrt(((sextracted8-ift.exp(s).val)**2).mean())
-    print "error back_size64:", np.sqrt(((sextracted64-ift.exp(s).val)**2).mean())
 
+    print "flux error energy1:", np.sqrt((((1 - ift.exp(s).val / ift.exp(energy1.s).val) ** 2)*d.val/d.sum()).mean())
+    print "flux error energy2:", np.sqrt((((1 - ift.exp(s).val / ift.exp(energy2.s).val) ** 2)*d.val/d.sum()).mean())
+    print "flux error energy3:", np.sqrt((((1 - ift.exp(s).val / ift.exp(energy3.s).val) ** 2)*d.val/d.sum()).mean())
+    print "flux error back_size8:", (np.sqrt(((1-sextracted8/ift.exp(s).val)**2)*d.val/d.sum()).mean())
+    print "flux error back_size64:", (np.sqrt(((1-sextracted64/ift.exp(s).val)**2)*d.val/d.sum()).mean())
+    print "flux error dcnn:", (np.sqrt(((1-dcnn_diffuse/ift.exp(s).val)**2)*d.val/d.sum()).mean())
+    print "class error energy1:", np.sqrt(((1 -  ift.exp(energy1.s).val/ift.exp(s).val) ** 2).mean())
+    print "class error energy2:", np.sqrt(((1 - ift.exp(energy2.s).val/ift.exp(s).val) ** 2).mean())
+    print "class error energy3:", np.sqrt(((1 - ift.exp(energy3.s).val/ift.exp(s).val) ** 2).mean())
+    print "class error back_size8:", np.sqrt(((1-sextracted8/ift.exp(s).val)**2).mean())
+    print "class error back_size64:", np.sqrt((((1-sextracted64/ift.exp(s).val)**2)).mean())
+    print "class error dcnn:", np.sqrt((((1-dcnn_diffuse/ift.exp(s).val)**2)).mean())

starblade/sugar.py

@@ -52,8 +52,8 @@ def build_starblade(data, alpha=1.5, q=1e-40, cg_iterations=100, newton_iteratio
     if manual_power_spectrum is None:
         initial_spectrum = ift.power_analyze(FFT.adjoint(ift.log(data)),
                                              binbounds=p_space.binbounds)
-        initial_spectrum /= 100*(p_space.k_lengths+1.)**4
-        # initial_spectrum = ift.Field(p_space,val=1e-3)
+        # initial_spectrum /= 100*(p_space.k_lengths+1.)**4
+        initial_spectrum = ift.Field(p_space,val=1e-8)
         update_power = True