new 3D plot and running the tests in paralell

Author: micha7a

.gitignore

@@ -2,3 +2,4 @@
 __pycache__
 .ipynb*
 surface-images/*
+*.npy

pattern.ipynb

@@ -289,7 +289,7 @@
  "metadata": {
   "anaconda-cloud": {},
   "kernelspec": {
-   "display_name": "Python [default]",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -303,7 +303,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.5.2"
+   "version": "3.6.0"
   }
  },
  "nbformat": 4,

results.ipynb

@@ -92,7 +92,7 @@ class ConstrainedALS(AlternatingLS):
 
     def __init__(self, samples, model_size, model_type, start_pos,
                  show_plots=False, hold_edges=True, stopping_error=1.0e-6, beta=0.01, interval_length=1, max_iter=10000,
-                 fl=1.0, angle=0):
+                 fl=1.0, angle=0, verbose=True):
         super(ConstrainedALS, self).__init__(samples, model_size, model_type, show_plots,
                                              hold_edges, stopping_error, beta, interval_length)
         assert len(samples) == len(start_pos)
@@ -104,6 +104,7 @@ def __init__(self, samples, model_size, model_type, start_pos,
             self.tr_param[2] = fl
             self.tr_param[0] = angle
         self.max_iterations = max_iter
+        self.verb = verbose
 
     def solve(self):
         sign = -1
@@ -120,15 +121,17 @@ def solve(self):
             self.error = np.linalg.norm(np.dot(x, self.parameter_estimate) - self.samples) / self.number_samples
 
             if self.error < self.stopping_error:
-                print("error small enough after fitting parameters")
+                if self.verb:
+                    print("error small enough after fitting parameters")
                 break
 
             g = self.model_type.compute_ls_gradient(self.start_positions, self.parameter_estimate, self.samples,
                                                     self.tr_param)
             # print "gradient: %f" % g
             # print "alpha: %f" % self.tr_param
             if np.max(np.abs(g*self.beta)) < np.finfo(float).eps:
-                print("converged to local minimum after", k, "steps")
+                if self.verb:
+                    print("converged to local minimum after", k, "steps")
                 break
 
             # normalize the gradient so it does not explode for many samples
@@ -137,11 +140,13 @@ def solve(self):
             error = np.linalg.norm(np.dot(x, self.parameter_estimate) - self.samples) / self.number_samples
 
             if error < self.stopping_error:
-                print("error small enough after fitting positions")
+                if self.verb:
+                    print("error small enough after fitting positions")
                 break
 
             if self.error < error:
-                print("error:", self.error, "beta:", self.beta)
+                if self.verb:
+                    print("error:", self.error, "beta:", self.beta)
                 if self.beta > 10*np.finfo(float).eps:
                     self.beta *= 0.5
 
@@ -161,7 +166,8 @@ def solve(self):
                 pylab.pause(0.01)
 
             if k == self.max_iterations - 1:
-                print('force stop after', self.max_iterations, 'steps')
+                if self.verb:
+                    print('force stop after', self.max_iterations, 'steps')
 
 
 class InvertedLS(OrdinaryLS):

solvers.py

@@ -1,60 +1,86 @@
 from samplers import *
 from solvers import *
 from plots import *
+from multiprocessing import Pool
 
 # set parameters
 tests = 50  # number of tests
-n = 4       # number of parameters of the polynomial (degree + 1)
-f = 1.0     # distance between the origin and the image plane
-b = 1.0     # intersection between camera axis and the surface
-slopes = np.linspace(-np.pi/6, np.pi/6, 13)
+pools = 5
+n = 6  # number of parameters of the polynomial (degree + 1)
+ovs = 1  # oversampling
+f = 1.0  # distance between the origin and the image plane
+b = 1.0  # intersection between camera axis and the surface
+slopes = np.linspace(-np.pi / 9, np.pi / 9, 13)
+# slopes = slopes[1:-1]
+noise_scale = 0
+noise_ampl = 0.0
+if noise_scale is not 0:
+    noise_ampl = 10.0 ** (-noise_scale)
+
 
 # containers for results:
-errors = []
-nsr = []
-
-for test in range(tests):
-    print("test: ", test)
-
-    tmp_err = []
-    start_param = nr.randn(n, 1)
-
-    print(start_param)
-
-    for slope in slopes:
-        polynomial = SecondSurfacePolynomial(start_param)
-        sampler = SurfaceSampler(polynomial, 2*n, [slope, b, f], interval_length=2, sigma=0.0, beg=-1)
-        noise = 1e-1*nr.randn(2*n)
-        nsr.append(np.linalg.norm(noise) / np.linalg.norm(sampler.sample_values))
-        sample_values = sampler.sample_values + noise
-
-        solver = ConstrainedALS(
-            sample_values,
-            polynomial.model_size,
-            SecondSurfacePolynomial,
-            start_pos=sampler.sample_positions,
-            stopping_error=1e-10,
-            beta=0.1,
-            show_plots=False,
-            max_iter=10000,
-            fl=f)
-
-        true_error = 0
-        try:
-            solver.solve()
-            true_error = abs(slope - solver.tr_param[0])
-        except AssertionError as as_err:
-            print("assertion error:", as_err.args[0])
-            true_error = np.NAN
-
-        tmp_err.append(true_error)
-    errors.append(tmp_err)
-
-errors = np.array(errors)
-nsr = np.array(nsr)
-
-print("mean:", np.nanmean(errors)*180/np.pi)
-print("median:", np.nanmedian(errors))
-print("std:", np.nanstd(errors))
-print("NANS:", str(np.count_nonzero(np.isnan(errors))/len(errors.flatten())*100)+"%")
+
+def test_block(beginning):
+    errors = []
+    nsr = []
+    params = []
+    np.random.seed(beginning)
+    for test in range(int(tests / pools)):
+        print("test: ", int(tests / pools) * beginning + test)
+
+        tmp_err = []
+        start_param = nr.randn(n, 1)
+        start_param[0] = 1
+
+        # print(start_param)
+        params.append(start_param)
+
+        for slope in slopes:
+            polynomial = SecondSurfacePolynomial(start_param)
+            sampler = SurfaceSampler(polynomial, 2 * ovs * n, [slope, b, f], interval_length=2, sigma=0.0, beg=0)
+            noise = noise_ampl * nr.randn(2 * ovs * n)
+            nsr.append(np.linalg.norm(noise) / np.linalg.norm(sampler.sample_values))
+            sample_values = sampler.sample_values + noise
+
+            solver = ConstrainedALS(
+                sample_values,
+                polynomial.model_size,
+                SecondSurfacePolynomial,
+                start_pos=sampler.sample_positions,
+                stopping_error=1e-10,
+                beta=0.1,
+                show_plots=False,
+                max_iter=10000,
+                fl=f,
+                verbose=False)
+
+            true_error = 0
+            try:
+                solver.solve()
+                true_error = abs(slope - solver.tr_param[0])
+            except AssertionError as as_err:
+                print("assertion error:", as_err.args[0])
+                true_error = np.NAN
+
+            tmp_err.append(true_error)
+        errors.append(np.array(tmp_err))
+    return (errors, nsr, params)
+
+
+p = Pool(pools)
+sth = p.map(test_block, range(pools))
+sth = np.array(sth)
+errors = np.concatenate(sth[:, 0])
+nsr = np.concatenate(sth[:, 1]).reshape(tests, len(slopes))
+params = np.concatenate(sth[:, 2]).reshape(tests, n)
+
+print("mean:", np.degrees(np.nanmean(errors)))
+print("median:", np.degrees(np.nanmedian(errors)))
+print("std:", np.degrees(np.nanstd(errors)))
+print("NANS:", str(np.count_nonzero(np.isnan(errors)) / len(errors.flatten()) * 100) + "%")
 print("noise to signal", np.nanmean(nsr))
+
+version = str(n)+"_"+str(ovs)+"_"+str(noise_scale)
+np.save("errors_"+version, errors)
+np.save("nsr_"+version, nsr)
+np.save("params_"+version, params)