more romhom fixes

bin/test_romhom.py

 from pykat.utilities.maps import read_map
 from pykat.utilities.knm import *
 from pykat.utilities.optics.gaussian_beams import beam_param
+import time
 
 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]]])
 
-q1 = beam_param(w0=2e-2, z=0)
-q2 = beam_param(w0=2e-2, z=0)
+q1 = beam_param(w0=4e-2, z=0)
+q2 = beam_param(w0=4e-2, z=0)
 
 m = read_map('etm08_virtual.txt')
 
-m.data = np.zeros((200, 200))
-m.x = np.linspace(-0.151, 0.149, m.size[0])
-m.y = np.linspace(-0.151, 0.149, m.size[1])
+m.data = np.zeros((100, 100))
+m.x = np.linspace(-0.15, 0.15, m.size[0])
+m.y = np.linspace(-0.15, 0.15, m.size[1])
 m.center = np.array(m.size)/2
 m.step_size = (m.x[1]-m.x[0], m.y[1]-m.y[0])
 
 print "generating weights..."
-w = m.generateROMWeights()
+t0 = time.time()
+w, EI = m.generateROMWeights(isModeMatched=True)
+print "Completed in ", time.time()-t0
 
 print "computing Riemann knm..."
+t0 = time.time()
 K1 = knmHG(couplings, q1, q2, surface_map=m)
+tr = time.time()-t0
+print "Completed in ", tr
 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 "Speed up", tr/tt
+
 #w.writeToFile("testWeights.rom")
 
 

pykat/utilities/knm.py

@@ -94,53 +94,47 @@ 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]
     
-    if cache == None:
-        u_x_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.nodes)   
-        u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.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
-        w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4
-        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_Q3 + weights.w_ij_Q2 + 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])
-        
-        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
-    m_mod_2 = m % 2
-    npr_mod_2 = npr % 2
-    mpr_mod_2 = mpr % 2
-    
-    if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
-        k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4) 	
-
-    elif n_mod_2 != m_mod_2:
-        if npr_mod_2 == mpr_mod_2:
-            k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q2Q3)
-        else:
-            k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3) 
-            
-    elif npr_mod_2 != mpr_mod_2:
-        if n_mod_2 == m_mod_2:
-            k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes,  w_ij_Q1Q2)
-        else:
-            k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
+    # if cache == None:
+    u_xm_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.EI["xm"].nodes-0.1e-2)
+    u_xp_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.EI["xp"].nodes-0.1e-2)
+    u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes-0.1e-2)
+    u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes-0.1e-2)
+    
+    # else:
+    #     strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+    #     stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
+    #
+    #     u_x_nodes = cache[strx]
+    #     u_y_nodes = cache[stry]
+    
+    u_xy_Q1 = np.outer(u_xm_nodes, u_yp_nodes)
+    u_xy_Q2 = np.outer(u_xp_nodes, u_yp_nodes)
+    u_xy_Q3 = np.outer(u_xp_nodes, u_ym_nodes)
+    u_xy_Q4 = np.outer(u_xm_nodes, u_ym_nodes)
+    
+    # n_mod_2 = n % 2
+    # m_mod_2 = m % 2
+    # npr_mod_2 = npr % 2
+    # mpr_mod_2 = mpr % 2
+    #
+    # if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2:
+    #     k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
+    #
+    # elif n_mod_2 != m_mod_2:
+    #     if npr_mod_2 == mpr_mod_2:
+    #         k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q2Q3)
+    #     else:
+    #         k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
+    #
+    # elif npr_mod_2 != mpr_mod_2:
+    #     if n_mod_2 == m_mod_2:
+    #         k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes,  w_ij_Q1Q2)
+    #     else:
+    #         k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
+    k_ROQ  = np.einsum('ij,ij', u_xy_Q1, weights.w_ij_Q1)
+    k_ROQ += np.einsum('ij,ij', u_xy_Q2, weights.w_ij_Q2)
+    k_ROQ += np.einsum('ij,ij', u_xy_Q3, weights.w_ij_Q3)
+    k_ROQ += np.einsum('ij,ij', u_xy_Q4, weights.w_ij_Q4)
      
     return k_ROQ
     

pykat/utilities/maps.py

@@ -2,6 +2,22 @@ from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, m
 import numpy
 import math
 
+class tiltmap(surfacemap):
+    
+    def __init__(self, name, size, step_size, tilt, center=(0,0) ):
+        surfacemap.__init__(name, "phase", size, center, step_sizes)
+        self.tilt = tilt
+        
+    @property
+    def tilt(self):
+        return self.__tilt
+    
+    @tilt.setter
+    def tilt(self, value):
+        
+        self.__tilt = value
+        
+        
 class surfacemap:
     def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
         
@@ -59,9 +75,24 @@ class surfacemap:
             raise BasePyKatException("Map type needs handling")
         
 
-    def generateROMWeights(self):
-        b = makeReducedBasis(self.x[0:(len(self.x)/2)], offset=self.offset)
-        EI = makeEmpiricalInterpolant(b)
+    def generateROMWeights(self, isModeMatched=True):
+        xp = self.x[self.x>0]
+        yp = self.y[self.y>0]
+        xm = self.x[self.x<0]
+        ym = self.y[self.y<0]
+        
+        EI = {}
+        
+        if len(xm) > 0: EI["xm"] = makeEmpiricalInterpolant(makeReducedBasis(xm, isModeMatched=isModeMatched))
+        if len(xp) > 0: EI["xp"] = makeEmpiricalInterpolant(makeReducedBasis(xp, isModeMatched=isModeMatched))
+        if len(ym) > 0: EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched))
+        if len(yp) > 0: EI["yp"] = makeEmpiricalInterpolant(makeReducedBasis(yp, isModeMatched=isModeMatched))
+
+        #if len(x0) > 0: EI["x0"] = makeEmpiricalInterpolant(makeReducedBasis(x0, isModeMatched=isModeMatched))
+        #if len(y0) > 0: EI["y0"] = makeEmpiricalInterpolant(makeReducedBasis(y0, isModeMatched=isModeMatched))
+        
+        EI["limits"] = EI["xm"].limits
+        
         self._rom_weights = makeWeights(self, EI)
         
         return self.ROMWeights, EI

pykat/utilities/romhom.py

@@ -18,12 +18,13 @@ ROMLimits = collections.namedtuple('ROMLimits', 'zmin zmax w0min w0max max_order
 
 class ROMWeights:
     
-    def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, nodes, limits):
+    def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, EI, limits):
         self.w_ij_Q1 = w_ij_Q1
         self.w_ij_Q2 = w_ij_Q2
         self.w_ij_Q3 = w_ij_Q3
         self.w_ij_Q4 = w_ij_Q4
-        self.nodes = nodes
+        
+        self.EI = EI
         self.limits = limits
         
     def writeToFile(self, filename):
@@ -111,7 +112,7 @@ def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x):
     return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x).conjugate()
 
     
-def makeReducedBasis(x, offset=np.array([0,0]), isModeMatched=True, tolerance = 1e-12, sigma = 1):
+def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
     
     if isModeMatched:
         greedypts = 'matched20.txt'
@@ -247,21 +248,33 @@ 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)
     
 
-def makeWeights(smap, EI):
+from progressbar import ProgressBar, ETA, Percentage, Bar
+
+def makeWeights(smap, EI, verbose=True):
+    
     
     # get full A_xy
     A_xy = smap.z_xy()
     
-    hx = len(smap.x)/2
-    hy = len(smap.y)/2
-    fx = hx*2
-    fy = hy*2
+    xm = smap.x[smap.x < 0]
+    xp = smap.x[smap.x > 0]
+    
+    ym = smap.y[smap.y < 0]
+    yp = smap.y[smap.y > 0]
+    
+    Q1xy = np.ix_(smap.x < 0, smap.y > 0)
+    Q2xy = np.ix_(smap.x > 0, smap.y > 0)
+    Q3xy = np.ix_(smap.x > 0, smap.y < 0)
+    Q4xy = np.ix_(smap.x < 0, smap.y < 0)
     
     # get A_xy in the four quadrants of the x-y plane
-    A_xy_Q1 = A_xy[0:hx][0:hy, hy:fy]
-    A_xy_Q2 = A_xy[hx:fx][0:hy, hy:fy]
-    A_xy_Q3 =  A_xy[hx:fx][0:hy, 0:hy]
-    A_xy_Q4 = A_xy[0:hx][0:hy, 0:hy]
+    A_xy_Q1 = A_xy[Q1xy]
+    A_xy_Q2 = A_xy[Q2xy]
+    A_xy_Q3 = A_xy[Q3xy]
+    A_xy_Q4 = A_xy[Q4xy]
+    
+    # Not sure if there is a neater way to select the x=0,y=0 cross elements
+    A_xy_0  = np.hstack([A_xy[:,smap.x==0].flatten(), A_xy[smap.y==0,:].flatten()])
     
     full_x = smap.x
     full_y = smap.y
@@ -269,29 +282,72 @@ def makeWeights(smap, EI):
     dx = full_x[1] - full_x[0]
     dy = full_y[1] - full_y[0]
 
-    B = EI.B
-    nodes_nv = EI.nodes 
-    nodes_pv = - EI.nodes
+    if verbose:
+        count  = len(EI["xm"].B) * len(EI["yp"].B)
+        count += len(EI["xp"].B) * len(EI["yp"].B)
+        count += len(EI["xp"].B) * len(EI["ym"].B)
+        count += len(EI["xm"].B) * len(EI["ym"].B)
+        p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
+    
+    n = 0
     
     # make integration weights
-    w_ij_Q1 = np.zeros((len(B), len(B)), dtype = complex)
-    w_ij_Q2 = np.zeros((len(B), len(B)), dtype = complex)
-    w_ij_Q3 = np.zeros((len(B), len(B)), dtype = complex)
-    w_ij_Q4 = np.zeros((len(B), len(B)), dtype = complex)
-
-    for i in range(len(B)):
-        for j in range(len(B)):
-            B_ij_Q1  = np.outer(B[i], B[j][::-1])		
+    Bx = EI["xm"].B
+    By = EI["yp"].B
+    w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
+    
+    for i in range(len(Bx)):
+        for j in range(len(By)):
+            B_ij_Q1 = np.outer(Bx[i], By[j])
             w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)	
 
-            B_ij_Q2 = np.outer(B[i][::-1], B[j][::-1])
+            if verbose:
+                p.update(n)
+                n+=1
+
+    Bx = EI["xp"].B
+    By = EI["yp"].B
+    w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex)
+    
+    for i in range(len(Bx)):
+        for j in range(len(By)):
+            B_ij_Q2 = np.outer(Bx[i], By[j])
             w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
 
-            B_ij_Q3 = np.outer(B[i][::-1], B[j])
+            if verbose:
+                p.update(n)
+                n+=1
+
+    Bx = EI["xp"].B
+    By = EI["ym"].B
+    w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)
+    
+    for i in range(len(Bx)):
+        for j in range(len(By)):
+            B_ij_Q3 = np.outer(Bx[i], By[j])
             w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
 
-            B_ij_Q4 = np.outer(B[i], B[j])
+            if verbose:
+                p.update(n)
+                n+=1
+
+    Bx = EI["xm"].B
+    By = EI["ym"].B
+    w_ij_Q4 = np.zeros((len(Bx),len(By)), dtype = complex)
+    
+    for i in range(len(Bx)):
+        for j in range(len(By)):
+            B_ij_Q4 = np.outer(Bx[i], By[j])
             w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
 
-    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, nodes=EI.nodes, limits=EI.limits)
+            if verbose:
+                p.update(n)
+                n+=1
+    if verbose:
+        p.finish()
+        
+    return ROMWeights(w_ij_Q1=w_ij_Q1, w_ij_Q2=w_ij_Q2, w_ij_Q3=w_ij_Q3, w_ij_Q4=w_ij_Q4, EI=EI, limits=EI["xm"].limits)
+    
+    
+    
     
\ No newline at end of file

requirements.txt

 numpy >= 1.6.2
 flask >= 0.10.1
+progressbar >= 2.2