updating knm computations and adding romhom tests

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
+
+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)
+
+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.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()
+
+print "computing Riemann knm..."
+K1 = knmHG(couplings, q1, q2, surface_map=m)
+print np.abs(K1)
+
+print "computing ROMHOM knm..."
+K2 = knmHG(couplings, q1, q2, surface_map=m, method="romhom")
+print np.abs(K2)
+
+#w.writeToFile("testWeights.rom")
+
+
+
+

pykat/utilities/knm.py

@@ -8,10 +8,7 @@ import numpy as np
 import pykat
 import collections
 import math
-
-def makeCouplingMatrix(max_order):
-    pass
-
+from romhom import u_star_u
 
 def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None):
 
@@ -39,11 +36,124 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
     return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
 
 
-def ROM_HG_knm(m, mode_in, mode_out, q1, q2, q1y=None, q2y=None):
-    pass
+def __gen_ROM_HG_knm_cache(weights, 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 = {}
+    
+    cache["w_ij_Q1Q3"] = weights.w_ij_Q1 + weights.w_ij_Q3
+    cache["w_ij_Q2Q4"] = weights.w_ij_Q2 + weights.w_ij_Q4
+    cache["w_ij_Q1Q2"] = weights.w_ij_Q1 + weights.w_ij_Q2
+    cache["w_ij_Q1Q4"] = weights.w_ij_Q1 + weights.w_ij_Q4
+    cache["w_ij_Q2Q3"] = weights.w_ij_Q2 + weights.w_ij_Q3
+    cache["w_ij_Q3Q4"] = weights.w_ij_Q3 + weights.w_ij_Q4
+    cache["w_ij_Q1Q2Q3Q4"] = weights.w_ij_Q1 + weights.w_ij_Q3 + weights.w_ij_Q2 + weights.w_ij_Q4
+    
+    while not it.finished:
+        try:
+            mode_in = [int(it.next()), int(it.next())]
+            mode_out = [int(it.next()), int(it.next())]
+            
+            strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+            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)
+            
+            
+            
+        except StopIteration:
+            break
+    
+    return cache
+    
+def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=None):
+    if q1y == None:
+        q1y = q1
+
+    if q2y == None:
+        q2y = q2
+    
+    # x modes
+    n = mode_in[0]
+    m = mode_out[0]
+
+    # y modes
+    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)
+    
+    return k_ROQ
+    
+
+
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann"):
+    if q1y == None:
+        q1y = q1
         
+    if q2y == None:
+        q2y = q2
         
-def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
+    assert q1.wavelength == q2.wavelength and q1y.wavelength == q2y.wavelength and q1y.wavelength == q1.wavelength
     
     couplings = np.array(couplings)
     
@@ -52,8 +162,7 @@ def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
     if int(a) - a != 0:
         raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")
         
-    if "phase" in surface_map.type:
-        Axy = np.exp(2j * 2 * math.pi / q1.wavelength * surface_map.data * surface_map.scaling)
+    Axy = surface_map.z_xy(q1.wavelength)
     
     x = surface_map.x
     y = surface_map.y
@@ -64,6 +173,20 @@ def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
     
     i = 0
     
+    if method == "romhom":
+        if surface_map == None:
+            raise BasePyKatException("Using 'romhom' method requires a surface map to be specified")
+            
+        weights = surface_map.ROMWeights
+        
+        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)
+    else:
+        romcache = None
+        weights = None
+        
     while not it.finished:
         try:
             
@@ -72,8 +195,10 @@ def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
             
             if method == "riemann":
                 K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy)
+            elif method == "romhom":
+                K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=romcache)
             else:
-                raise BasePyKatException("method value not accepted")
+                raise BasePyKatException("method value '%s' not accepted" % method)
         
             i +=1
                  

pykat/utilities/maps.py

 from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
 import numpy
+import math
 
 class surfacemap:
     def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
@@ -7,7 +8,6 @@ class surfacemap:
         self.name = name
         self.type = maptype
         self.center = center
-        self.size = size
         self.step_size = step_size
         self.scaling = scaling
         self.data = data
@@ -19,7 +19,7 @@ class surfacemap:
             mapfile.write("% Surface map\n")
             mapfile.write("% Name: {0}\n".format(self.name))
             mapfile.write("% Type: {0}\n".format(self.type))
-            mapfile.write("% Size: {0} {1}\n".format(self.size[0], self.size[1]))
+            mapfile.write("% Size: {0} {1}\n".format(self.data.shape[0], self.data.shape[1]))
             mapfile.write("% Optical center (x,y): {0} {1}\n".format(self.center[0], self.center[1]))
             mapfile.write("% Step size (x,y): {0} {1}\n".format(self.step_size[0], self.step_size[1]))
             mapfile.write("% Scaling: {0}\n".format(float(self.scaling)))
@@ -38,6 +38,10 @@ class surfacemap:
     def y(self):
         return self.step_size[1] * (numpy.array(range(1, self.data.shape[1]+1))- self.center[1])
 
+    @property
+    def size(self):
+        return self.data.shape
+            
     @property
     def offset(self):
         return numpy.array(self.step_size)*(self.center - numpy.array(self.size)/2)
@@ -46,12 +50,21 @@ class surfacemap:
     def ROMWeights(self):
         return self._rom_weights
     
+    def z_xy(self, wavelength=1064e-9):
+        
+        if "phase" in self.type:
+            k = math.pi * 2 / wavelength
+            return numpy.exp(2j * k * self.scaling * self.data)
+        else:
+            raise BasePyKatException("Map type needs handling")
+        
+    
     def generateROMWeights(self):
         b = makeReducedBasis(self.x[0:(len(self.x)/2)], offset=self.offset)
         EI = makeEmpiricalInterpolant(b)
         self._rom_weights = makeWeights(self, EI)
         
-        return self.ROMWeights
+        return self.ROMWeights, EI
         
     def plot(self, show=True, clabel=None):
         

pykat/utilities/optics/gaussian_beams.py

@@ -82,7 +82,7 @@ class gauss_param(object):
     def w(self):
         return abs(self.__q)*math.sqrt(self.__lambda / (self.__nr * math.pi * self.__q.imag))
     
-    def w(self, z=None, wavelength=None, nr=None, w0=None):
+    def beamsize(self, z=None, wavelength=None, nr=None, w0=None):
 
         if z == None:
             z = self.z
@@ -144,7 +144,7 @@ class gauss_param(object):
         else:
             return float("inf")
     
-    def Rc(self, z=None, wavelength=None, nr=None, w0=None):
+    def curvature(self, z=None, wavelength=None, nr=None, w0=None):
         if z == None:
             z = self.z
         else:

pykat/utilities/romhom.py

+import math
 import maps
 import os.path
-import numpy as np
 import pykat
 import collections
+
 from itertools import combinations_with_replacement as combinations
 from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
 from scipy.linalg import inv
 from math import factorial
 from hermite import *
-import math
+
+import numpy as np
 
 EmpiricalInterpolant = collections.namedtuple('EmpiricalInterpolant', 'B nodes node_indices limits x')
 ReducedBasis = collections.namedtuple('ReducedBasis', 'RB limits x')
@@ -16,11 +18,12 @@ 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, limits):
+    def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, nodes, 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.limits = limits
         
     def writeToFile(self, filename):
@@ -101,7 +104,7 @@ def u(re_q1, w0_1, n1, x):
 
     wz1 = w(w0_1, im_q1, re_q1)
 
-    return A_n1 * hermite(n1, np.sqrt(2.)*x / wz1) * np.exp(-1j*(2*math.pi/(1064e-9))* x**2 /(2.*q_z1))
+    return A_n1 * hermite(n1, np.sqrt(2.)*x / wz1) * np.exp(np.array(-1j*(2*math.pi/(1064e-9))* x**2 /(2.*q_z1)))
 
 
 def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x):
@@ -241,13 +244,13 @@ def makeEmpiricalInterpolant(RB):
     invV = inv(V[0:len(node_indices), 0:len(node_indices)])
     B = B_matrix(invV, e_x)
     
-    return EmpiricalInterpolant(B=np.array(B), nodes=np.array(x_nodes), node_indices=np.array(node_indices),  limits=RB.limits, x=RB.x)
+    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):
     
     # get full A_xy
-    A_xy = smap.data
+    A_xy = smap.z_xy()
     
     hx = len(smap.x)/2
     hy = len(smap.y)/2
@@ -290,82 +293,5 @@ def makeWeights(smap, EI):
             B_ij_Q4 = np.outer(B[i], B[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, limits=EI.limits)
-    
-    
-def ROMKnm(W, max_order, q1, q2, q1y=None, q2y=None):
-    
-    if q1y == None:
-        q1y = q1
-        
-    if q2y == None:
-        q2y = q2
-        
-    # get symmetric and anti-symmetric w_ij's 
-    w_ij_Q1Q3 = W.w_ij_Q1 + W.w_ij_Q3
-    w_ij_Q2Q4 = W.w_ij_Q2 + W.w_ij_Q4
-    w_ij_Q1Q2 = W.w_ij_Q1 + W.w_ij_Q2
-    w_ij_Q1Q4 = W.w_ij_Q1 + W.w_ij_Q4
-    w_ij_Q2Q3 = W.w_ij_Q2 + W.w_ij_Q3
-    w_ij_Q3Q4 = W.w_ij_Q3 + W.w_ij_Q4
-
-    w_ij_Q1Q2Q3Q4 = W.w_ij_Q1 + W.w_ij_Q3 + W.w_ij_Q2 + W.w_ij_Q4
-
-    num_fields = int((max_order + 1) * (max_order + 2) / 2);
-
-    K_ROQ = np.array((num_fields, num_fields))
-
-    for i in range(len(nm_pairs)):
-        for j in range(len(nm_pairs)):
-
-            # full quadrature
-            n = nm_pairs[i][0]
-            m = nm_pairs[i][1]
-            npr = nm_pairs[j][0]
-            mpr = nm_pairs[j][1]
-            u_xy =  np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, full_x), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, full_x))
-            k = dx*dy*np.einsum('ij,ij', u_xy, A_xy)	
-            
-            # ROQ
-            if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:
-                u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
-                k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4) 	
-	
-            elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
-                u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
-                k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
-            elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
-                u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
-                k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
-            elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:	
-                u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
-                k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-                
-            elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 1 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 0:
-                u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n,m, nodes_nv)	
-                u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)   
-
-                if nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
-                    u_xy_nodes_Q1Q4 = np.outer(u_x_nodes, u_y_nodes)
-                    u_xy_nodes_Q2Q3 = - u_xy_nodes_Q1Q4
-                    k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q1Q4, w_ij_Q1Q4) + np.einsum('ij,ij', u_xy_nodes_Q2Q3, w_ij_Q2Q3)
-                else:
-                    u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
-                    u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
-                    k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3) 
-
-            elif nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 1 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 0:
-                u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv)   
-                u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)  
-                
-                if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1:
-                    u_xy_nodes_Q3Q4 = np.outer(u_x_nodes, u_y_nodes)
-                    u_xy_nodes_Q1Q2 = - u_xy_nodes_Q3Q4
-                    k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q3Q4, w_ij_Q3Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q2,  w_ij_Q1Q2)
-                else:
-                    u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
-                    u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
-                    k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3)
+    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)
     
\ No newline at end of file