adding bayerhelms knm calculator, fixing bits with maps and romhom

pykat/__init__.py

@@ -5,3 +5,4 @@ import components
 import detectors
 import commands
 
+from pykat.utilities.optics.gaussian_beams import beam_param

pykat/utilities/knm.py

@@ -2,12 +2,14 @@ from itertools import combinations_with_replacement as combinations
 from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
 from pykat.exceptions import BasePyKatException
 
+import time
 import maps
 import os.path
 import numpy as np
 import pykat
 import collections
 import math
+import cmath
 
 from romhom import u_star_u
 from progressbar import ProgressBar, ETA, Percentage, Bar
@@ -34,7 +36,7 @@ def makeCouplingMatrix(max_order):
     
     return np.array(M)
 
-def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None):
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0)):
 
     if Axy == None:
         Axy == np.ones((len(x), len(y)))
@@ -55,7 +57,7 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=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)
+        U1 = Hg_in.Unm(x+delta[0], y+delta[1])
         U2 = Hg_out.Unm(x,y).conjugate()
         
         return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
@@ -69,8 +71,7 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
 
 
     
-
-def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
+def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0)):
     if q1y == None:
         q1y = q1
         
@@ -89,20 +90,24 @@ def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
             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])
+            #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()   
+                cache[strx] = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, mode_in[0], mode_out[0], x, x+delta[0])    
+                #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()
+                cache[stry] = u_star_u(q1y.z,   q2y.z,  q1y.w0,  q2y.w0, mode_in[1], mode_out[1], y, y+delta[1])    
+                #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:
@@ -142,6 +147,8 @@ def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
     
     return cache
 
+
+
 def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=None):
     if q1y == None:
         q1y = q1
@@ -157,10 +164,6 @@ 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)
-    
-    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(q1y.z,   q2y.z,  q1y.w0,  q2y.w0, npr, mpr,   weights.EI["ym"].nodes)
@@ -211,28 +214,104 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
         else:
             k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3)
     
-    else:
-        # If there is no symmetry about the axes then each quadrant is different
-        u_xm_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.EI["xm"].nodes)
-        u_xp_nodes = u_star_u(q1.z,   q2.z,  q1.w0,  q2.w0, n,   m,   weights.EI["xp"].nodes)
-        u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes)
-        u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes)
-        
-        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)
-        
-        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
 
+__fac_cache = []
+
+def fac(n):
+    global __fac_cache
+    return __fac_cache[n]
+    #return math.factorial(int(n))
+
+def m_1_pow(n):
+    if n % 2 == 0:
+        return 1
+    else:
+        return -1
+
+
+def __Ss(u, _u, F, _F, d=0):
+    r = 0
+    
+    for s in range(0, min(u,_u)+1):
+        r += m_1_pow(s) * _F ** (u-s) * _F ** (_u-s) / (fac(2*s+d)*fac(u-s)*fac(_u-s))
+        
+    return r
 
 
-def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=True):
+def __S(m, _m, X, _X, F, _F, d=0):
+    if m % 2 == 1:
+        lim1 = (m-1)/2
+    else:
+        lim1 = m/2 
+
+    if _m % 2 == 1:
+        lim2 = (_m-1)/2
+    else:
+        lim2 = _m/2
+    
+    r = 0
+    
+    for u in range(0, lim1+1):
+        for _u in range(0, lim2+1):
+            r += m_1_pow(u) * _X**(m-2*u) * X**(_m-2*_u) / ( fac(m-2*u)*fac(_m-2*_u) )   * __Ss(u, _u, F, _F, d=d)
+    
+    return r
+           
+
+def __bayerhelms_kn(n, _n, q1, q2, gamma=0.0):
+    
+    K0 = (q1.zr - q2.zr)/q2.zr
+    K2 = (q1.z - q2.z)/q2.zr
+    K = (K0 + 1j*K2)/2.0
+    
+    Ktilde = abs(K / (1+K))
+
+    if gamma != 0:
+        a  = q2.zr * math.sin(gamma) / (cmath.sqrt(1+K.conjugate()) * q2.w0)
+
+        _X = - a * (q2.z/q2.zr - 1j)
+        X  = - a * (q2.z/q2.zr + 1j*(1+2*K.conjugate()))
+        Ex = cmath.exp(-_X*X / 2.0)
+    else:
+        _X = 0.0
+        X  = 0.0
+        Ex = 1.0
+    
+    _F  = K / (2.0 * (1.0+K0))
+    F = K.conjugate() / 2.0 
+
+    Sg = __S(n, _n, X, _X, F, _F)
+
+    if n > 0 and _n > 0:
+        Su = __S(n-1, _n-1, X, _X, F, _F, 1)
+    else:
+        Su = 0
+    
+    b = m_1_pow(_n) * cmath.sqrt(fac(n) * fac(_n) * (1.0 + K.real)**(n+0.5) * (1.0 + K.conjugate()) ** (-(n+_n+1)))
+    
+    return b * Ex * (Sg - Su)
+
+
+def bayerhelms_HG_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, gamma=(0,0)):
+    if q1y == None:
+        q1y = q1
+
+    if q2y == None:
+        q2y = q2
+
+    # x modes
+    n = mode_in[0]
+    _n = mode_out[0]
+
+    # y modes
+    m = mode_in[1]
+    _m = mode_out[1]
+
+    return __bayerhelms_kn(n,_n, q1, q2, 2*gamma[0]) * __bayerhelms_kn(m, _m, q1y, q2y, 2*gamma[1])
+
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0)):
     if q1y == None:
         q1y = q1
         
@@ -248,6 +327,18 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
     if int(a) - a != 0:
         raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")
     
+    maxtem = 0
+    c = couplings.flatten()
+    
+    for i in range(0, c.size/2):
+        maxtem = max(sum(c[i*2:(i*2+2)]), maxtem)
+    
+    global __fac_cache
+    
+    for n in range(0, maxtem+1):
+        __fac_cache.append(math.factorial(n))
+    
+    if surface_map != None:  
         Axy = surface_map.z_xy(q1.wavelength)
     
         x = surface_map.x
@@ -259,6 +350,9 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
     
     i = 0
     
+    if profile:
+        t0 = time.time()
+        
     if method == "romhom":
         if surface_map == None:
             raise BasePyKatException("Using 'romhom' method requires a surface map to be specified")
@@ -268,32 +362,44 @@ 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.")
 
-        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)
+        if surface_map == None:
+            raise BasePyKatException("Using 'riemann' method requires a surface map to be specified")
+            
+        cache = __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=delta)
     else:
         cache = None
         weights = None
     
+    if profile:
+        cache_time = time.time() - t0
+        Ktime = np.zeros((couplings.size/4,), dtype=np.float64)
+    
     if verbose:
         p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])
     
     while not it.finished:
         try:
+            if profile:
+                t0 = time.time()
+                
             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, cache=cache)
+                K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache, delta=delta)
             elif method == "romhom":
                 K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache)
+            elif method == "bayerhelms":
+                K[i] = bayerhelms_HG_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, gamma=gamma)
             else:
                 raise BasePyKatException("method value '%s' not accepted" % method)
             
+            if profile:
+                Ktime[i] = time.time() - t0
+            
             i +=1
             
             if verbose:
@@ -303,19 +409,41 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
         except StopIteration:
             break
 
+    if profile:
+        return K.reshape(couplings.shape[:-1]), Ktime.reshape(couplings.shape[:-1]), cache_time
+    else:
         return K.reshape(couplings.shape[:-1])
 
 
+def plot_knm_matrix(couplings, knm):
+    import pylab as plt
     
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
+    cax = ax.imshow(knm, interpolation='nearest')
+    fig.colorbar(cax)
     
+    numrows, numcols = knm.shape
     
+    c = couplings[:, 0, :2]
+    c_ = []
     
+    for d in c:
+        c_.append("%i,%i"%(d[0], d[1]))
     
+    ax.set_xticklabels(c_)
+    ax.set_yticklabels(c_)
     
+    def format_coord(x, y):
+        col = int(x+0.5)
+        row = int(y+0.5)
         
+        if col>=0 and col<numcols and row>=0 and row<numrows:
+            z = knm[row,col]
+            return 'x=%s, y=%s, z=%1.4f' % (c_[col], c_[row], z)
+        else:
+            return 'x=%1.4f, y=%1.4f'%(x, y)
 
+    ax.format_coord = format_coord
 
-
-
-
-
+    plt.show()

pykat/utilities/maps.py

@@ -65,25 +65,18 @@ class surfacemap(object):
             raise BasePyKatException("Map type needs handling")
         
 
-    def generateROMWeights(self, isModeMatched=True):
-        xp = self.x[self.x>0]
-        yp = self.y[self.y>0]
+    def generateROMWeights(self, isModeMatched=True, verbose=False):
         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)
+        self._rom_weights = makeWeights(self, EI, verbose=verbose)
         
         return self.ROMWeights, EI
         
@@ -129,7 +122,7 @@ class tiltmap(surfacemap):
         
         xx, yy = np.meshgrid(self.x, self.y)
         
-        self.data = xx * self.tilt[0] + yy * self.tilt[1]
+        self.data = xx * self.tilt[1] + yy * self.tilt[0]
         
         
         

pykat/utilities/romhom.py

 import math
-import maps
 import os.path
 import pykat
 import collections
@@ -109,8 +108,11 @@ def u(re_q1, w0_1, n1, x):
     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):
-    return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x).conjugate()
+def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x, x2=None):
+    if x2 == None:
+        x2 = x
+        
+    return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x2).conjugate()
 
     
 def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
@@ -150,14 +152,13 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
     IDx = 0 #TS index 
 
     for i in range(len(params)):
-
-        q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
-
         if isModeMatched:
+            q1 = beam_param(w0=float(params[i][0]), z=float(params[i][1]))
             q2 = q1
             n = int(params[i][2])
             m = int(params[i][3])
         else:
+            q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
             q2 = beam_param(w0=float(params[i][1]), z=float(params[i][3]))            
             n = int(params[i][4])
             m = int(params[i][5])
@@ -180,10 +181,6 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
 
     #### Begin greedy: see Field et al. arXiv:1308.3565v2 #### 
 
-    tolerance = 1e-12 # set maximum RB projection error
-
-    sigma = 1 # projection error at 0th iteration
-
     RB_matrix = [TS[0]] # seed greedy algorithm (arbitrary)
 
     proj_coefficients = np.zeros(TS_size*TS_size, dtype = complex).reshape(TS_size, TS_size)
@@ -196,6 +193,7 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
     	# go through training set and find worst projection error
     	projections = project_onto_basis(dx, RB_matrix, TS, projections, proj_coefficients, iter)
     	residual = TS - projections
+        
     	projection_errors = [np.vdot(dx* residual[i], residual[i]) for i in range(len(residual))]
     	sigma = abs(max(projection_errors))
     	index = np.argmax(projection_errors) 
@@ -208,6 +206,7 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
         
     	iter += 1
         
+        
     return ReducedBasis(RB=np.array(RB_matrix), limits=romlimits, x=x)
     
 def makeEmpiricalInterpolant(RB):
@@ -282,17 +281,14 @@ def makeWeights(smap, EI, verbose=True):
     dy = full_y[1] - full_y[0]
 
     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)
+        count  = 4*len(EI["xm"].B) * len(EI["ym"].B)
         p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
     
     n = 0
     
     # make integration weights
     Bx = EI["xm"].B
-    By = EI["yp"].B
+    By = EI["ym"].B[:,::-1]
     w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
     
     for i in range(len(Bx)):
@@ -304,8 +300,8 @@ def makeWeights(smap, EI, verbose=True):
                 p.update(n)
                 n+=1
 
-    Bx = EI["xp"].B
-    By = EI["yp"].B
+    Bx = EI["xm"].B[:,::-1]
+    By = EI["ym"].B[:,::-1]
     w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex)
     
     for i in range(len(Bx)):
@@ -317,7 +313,7 @@ def makeWeights(smap, EI, verbose=True):
                 p.update(n)
                 n+=1
 
-    Bx = EI["xp"].B
+    Bx = EI["xm"].B[:,::-1]
     By = EI["ym"].B
     w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)