adding riemann caching for fair comparison between romhoom

bin/test_romhom.py

@@ -2,28 +2,27 @@ from pykat.utilities.maps import read_map, tiltmap, surfacemap
 from pykat.utilities.knm import *
 from pykat.utilities.optics.gaussian_beams import beam_param
 import time
+import numpy as np
 
-couplings = np.array([[[0,0,0,0], [0,0,1,0], [0,0,0,1]],
-                      [[1,0,0,0], [1,0,1,0], [1,0,0,1]],
-                      [[0,1,0,0], [0,1,1,0], [0,1,0,1]]])
+couplings = makeCouplingMatrix(5)
 
 q1 = beam_param(w0=3e-2, z=0)
 q2 = beam_param(w0=3e-2, z=0)
 
-size = np.array([200,200])
-stepsize = 0.2/(size-1)
+size = np.array([1200, 1200])
+stepsize = 0.3/(size-1)
 
 # This map has points evenly spaced about the axes
-m = tiltmap("tilt", size, stepsize, (1e-6, 0))
+m = tiltmap("tilt", size, stepsize, (1e-6, 1e-6))
 # Shifting the central point changes the results slightly
 # but this is about a 1e-7 change, probably due to extra
 # clipping
-m.center += 20
+m.center += 0
 
 print "Generating weights..."
 t0 = time.time()
 w, EI = m.generateROMWeights(isModeMatched=True)
-#w.writeToFile("testWeights.rom")
+w.writeToFile("testWeights.rom")
 print "Completed in ", time.time()-t0
 
 print "computing Riemann knm..."
@@ -31,15 +30,14 @@ t0 = time.time()
 K1 = knmHG(couplings, q1, q2, surface_map=m)
 tr = time.time()-t0
 print "Completed in ", tr
-print np.abs(K1)
+#print np.abs(K1)
 
 print "computing ROMHOM knm..."
 t0 = time.time()
 K2 = knmHG(couplings, q1, q2, surface_map=m, method="romhom")
 tt = time.time()-t0
 print "Completed in ", tt
-print np.abs(K2)
-
+#print np.abs(K2)
 print "Speed up", tr/tt
 
 

pykat/utilities/knm.py

@@ -8,9 +8,33 @@ import numpy as np
 import pykat
 import collections
 import math
+
 from romhom import u_star_u
+from progressbar import ProgressBar, ETA, Percentage, Bar
+
+def makeCouplingMatrix(max_order):
+    max_order = int(max_order)
+    c = []
+    for n in range(0, max_order+1):
+        for m in range(0, max_order+1):
+            if n+m <= max_order:
+                c.append([n,m])
+
+    M = []
+
+    for i in c:
+        row = []
+    
+        for j in c:
+            e = list(i)
+            e.extend(j)
+            row.append(e)
         
-def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None):
+        M.append(row)
+    
+    return np.array(M)    
+
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None):
 
     if Axy == None:
         Axy == np.ones((len(x), len(y)))
@@ -24,18 +48,61 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
     if len(mode_in) != 2 or len(mode_out) != 2:
         raise BasePyKatException("Both mode in and out should be a container with modes [n m]")        
 
-    Hg_in  = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
-    Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
-    
     dx = abs(x[1] - x[0])
     dy = abs(y[1] - y[0])    
         
+    if cache == None:
+        Hg_in  = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+        Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
+    
         U1 = Hg_in.Unm(x,y)
         U2 = Hg_out.Unm(x,y).conjugate()
         
         return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+    else:
+        
+        strx = "u1[%i,%i]" % (mode_in[0], mode_out[0])
+        stry = "u2[%i,%i]" % (mode_in[1], mode_out[1])
+        
+        return dx * dy * np.einsum('ij,ij', Axy, np.outer(cache[strx], cache[stry]))
+
+    
 
 
+
+def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
+    if q1y == None:
+        q1y = q1
+        
+    if q2y == None:
+        q2y = q2
+        
+    it = np.nditer(couplings, flags=['refs_ok','f_index'])
+    
+    cache = {}
+    
+    while not it.finished:
+        try:
+            mode_in = [int(it.next()), int(it.next())]
+            mode_out = [int(it.next()), int(it.next())]
+            
+            strx = "u1[%i,%i]" % (mode_in[0], mode_out[0])
+            stry = "u2[%i,%i]" % (mode_in[1], mode_out[1])
+            
+            Hg_in  = HG_beam(qx=q1, qy=q1y, n=mode_in[0], m=mode_in[1])
+            Hg_out = HG_beam(qx=q2, qy=q2y, n=mode_out[0], m=mode_out[1])
+    
+            if strx not in cache:
+                cache[strx] = Hg_in.Un(x) * Hg_out.Un(x).conjugate()   
+            
+            if stry not in cache:
+                cache[stry] = Hg_in.Um(y) * Hg_out.Um(y).conjugate()
+            
+        except StopIteration:
+            break
+    
+    return cache
+    
 def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
 
     if q1y == None:
@@ -65,14 +132,10 @@ def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
             stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
             
             if strx not in cache:
-                cache[strx] = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, mode_in[0], mode_out[0], weights.nodes)   
-            
-            if q1 == q1y and q2 == q2y:
-                cache[stry] = cache[strx]
-            elif stry not in cache:
-                cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], weights.nodes)
-            
+                cache[strx] = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, mode_in[0], mode_out[0], weights.EI["xm"].nodes)   
             
+            if stry not in cache:
+                cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], weights.EI["ym"].nodes)
             
         except StopIteration:
             break
@@ -94,12 +157,14 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
     npr = mode_in[1]
     mpr = mode_out[1]
     
-    foundSymmetry = np.all(weights.EI["ym"].nodes == -weights.EI["xm"].nodes) and np.all(weights.EI["xp"].nodes == -weights.EI["xm"].nodes) and np.all(weights.EI["yp"].nodes == -weights.EI["ym"].nodes)
+    foundSymmetry = np.all(weights.EI["ym"].nodes == weights.EI["xm"].nodes) and np.all(weights.EI["xp"].nodes == -weights.EI["xm"].nodes) and np.all(weights.EI["yp"].nodes == -weights.EI["ym"].nodes)
     
     if foundSymmetry:
+        
         if cache == None:
             u_x_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,     m,   weights.EI["xm"].nodes)
-            u_y_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.EI["ym"].nodes)
+            u_y_nodes = u_star_u(q1y.z,   q2y.z,  q1y.w0,  q2y.w0, npr, mpr,   weights.EI["ym"].nodes)
+            
             w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3
             w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4
             w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2
@@ -107,6 +172,7 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
             w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3
             w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4
             w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4
+            
         else:
             strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
             stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
@@ -114,6 +180,15 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
             u_x_nodes = cache[strx]
             u_y_nodes = cache[stry]
             
+            w_ij_Q1Q3 = cache["w_ij_Q1Q3"]
+            w_ij_Q2Q4 = cache["w_ij_Q2Q4"]
+            w_ij_Q1Q2 = cache["w_ij_Q1Q2"]
+            w_ij_Q1Q4 = cache["w_ij_Q1Q4"]
+            w_ij_Q2Q3 = cache["w_ij_Q2Q3"]
+            w_ij_Q3Q4 = cache["w_ij_Q3Q4"]
+            w_ij_Q1Q2Q3Q4 = cache["w_ij_Q1Q2Q3Q4"]
+            
+
         u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
     
         n_mod_2 = n % 2
@@ -157,7 +232,7 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
     
 
 
-def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann"):
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=True):
     if q1y == None:
         q1y = q1
         
@@ -193,27 +268,54 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
         if weights == None:
             raise BasePyKatException("The ROM weights need to be generated for this map before use.")
 
-        romcache = None#__gen_ROM_HG_knm_cache(weights, couplings, q1=q1, q2=q2, q1y=q1y, q2y=q2y)
+        if verbose: print "Computing caches"
+        
+        cache = __gen_ROM_HG_knm_cache(weights, couplings, q1=q1, q2=q2, q1y=q1y, q2y=q2y)
+        
+    elif method == "riemann":
+        if verbose: print "Computing caches"
+        cache = __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None)
     else:
-        romcache = None
+        cache = None
         weights = None
     
+    if verbose:
+        p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])
+    
     while not it.finished:
         try:
-            
             mode_in = [int(it.next()), int(it.next())]
             mode_out = [int(it.next()), int(it.next())]
             
             if method == "riemann":
-                K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy)
+                K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache)
             elif method == "romhom":
-                K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=romcache)
+                K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache)
             else:
                 raise BasePyKatException("method value '%s' not accepted" % method)
                 
             i +=1
             
+            if verbose:
+                p.update(i*4)
+                
+                 
         except StopIteration:
             break
 
     return K.reshape(couplings.shape[:-1])
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

pykat/utilities/romhom.py

@@ -4,6 +4,7 @@ import os.path
 import pykat
 import collections
 
+from progressbar import ProgressBar, ETA, Percentage, Bar
 from itertools import combinations_with_replacement as combinations
 from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
 from scipy.linalg import inv
@@ -150,21 +151,22 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
 
     for i in range(len(params)):
 
-        q1 = beam_param(w0=params[i][0], z=params[i][2])
+        q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
 
         if isModeMatched:
             q2 = q1
             n = int(params[i][2])
             m = int(params[i][3])
-            w0_1 = params[i][0]
-            w0_2 = w0_1
-            re_q1 = params[i][1]
-            re_q2 = re_q1
         else:
-            q2 = beam_param(w0=params[i][1], z=params[i][3])            
+            q2 = beam_param(w0=float(params[i][1]), z=float(params[i][3]))            
             n = int(params[i][4])
             m = int(params[i][5])
             
+        w0_1 = q1.w0
+        w0_2 = q2.w0
+        re_q1 = q1.z
+        re_q2 = q2.z
+            
         TS[IDx] = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, x)
         
         # normalize
@@ -248,11 +250,8 @@ def makeEmpiricalInterpolant(RB):
     return EmpiricalInterpolant(B=np.array(B), nodes=np.array(x_nodes, dtype=np.float64), node_indices=np.array(node_indices, dtype=np.int32), limits=RB.limits, x=RB.x)
     
 
-from progressbar import ProgressBar, ETA, Percentage, Bar
-
 def makeWeights(smap, EI, verbose=True):
     
-    
     # get full A_xy
     A_xy = smap.z_xy()