Added code to verify result: read u1.mat containing matlab result matlab to compare results

TestSuite/Phase_JWST_demo.py

@@ -5,6 +5,6 @@ import JWST_data_processor
 import JWST_in
 from phase import Phase
 
-JWST = Phase(JWST_in.new_config);
-JWST.solve();
+JWST = Phase(JWST_in.new_config)
+JWST.solve()
 JWST.show()

proxtoolbox/Problems/Phase/phase.py

@@ -6,6 +6,8 @@ from proxtoolbox import ProxOperators
 from proxtoolbox.ProxOperators.proxoperators import ProxOperator
 from proxtoolbox.Problems.Phase.JWST_graphics import JWST_graphics
 from numpy.linalg import norm
+import numpy as np
+import h5py
 from numpy import square, sqrt, nonzero
 
 
@@ -211,7 +213,33 @@ class Phase(Problem):
         #print(self.u1[71,46]);
         #print(self.gap);
         #print(self.change);
-    
+        f = h5py.File('u1.mat')
+        u1 = f['u1'].value.view(np.complex)
+        u1 =np.array(u1)
+        u1 = u1.T
+        print("Compare u1:")
+        print("Nonzero indices matlab:")
+        print(nonzero(u1))
+        print("Nonzero indices python:")
+        print(nonzero(self.output['u1']))
+        print("Nonzero indices equal:")
+        print(np.array_equal(nonzero(u1),nonzero(self.output['u1'])))
+        print("Nonzero values matlab:")
+        print(u1[nonzero(u1)])
+        print("Nonzero values python:")
+        print(self.output['u1'][nonzero(self.output['u1'])])
+        print("Difference at nonzero values:")
+        nonz = nonzero(u1)
+        diff = u1 - self.output['u1']
+        print(diff[nonz])
+        print("Maximum of absolute value of difference:")
+        print(np.amax(abs(diff)));
+        print("Frobenius of difference:")
+        print(norm(diff))
+        print("Frobenius norm of matlab u1:")
+        print(norm(u1))
+        print("Frobenius norm of python u1:")
+        print(norm(self.output['u1']))
     
     def _postsolve(self):
         """
@@ -224,7 +252,8 @@ class Phase(Problem):
         """
         Generates graphical output from the solution
         """
-
+        
+        print("Calculation time:")
         print(self.elapsed_time)
         JWST_graphics(self.config,self.output)