adding adaptive knm solver, it's pretty slow though. Adding aperture map

.gitignore

 *.pyc
+.DS_store
+*.dat
+*.txt
+*.log
+*.kat
+*.pdf
+*.out
+*.bak
\ No newline at end of file

bin/test_knm.py

+from pykat import *
+from pykat.utilities.knm import knmHG, makeCouplingMatrix, plot_knm_matrix
+from pykat.utilities.maps import aperturemap
+import numpy as np
+import time
+
+q1 = beam_param(w0=5e-2, z=0)
+q2 = beam_param(w0=5e-2, z=0)
+
+aperture = 0.1
+s = 1000
+size = np.array([s, s])
+stepsize = 0.3/(size-1)
+smap = aperturemap("tilt", size, stepsize, aperture)
+
+couplings = makeCouplingMatrix(1)
+
+params = {"usepolar":True, "aperture":aperture, "epsabs": 1e-3, "epsrel": 1e-3}
+
+t0 = time.time()
+kbh = knmHG(couplings, q1, q2, method="adaptive", verbose=True, params=params)
+print time.time() - t0
+
+t0 = time.time()
+kr = knmHG(couplings, q1, q2, surface_map=smap, method="riemann", verbose=True)
+print time.time() - t0
+
+smap.generateROMWeights(isModeMatched=True, verbose=True)
+t0 = time.time()
+krm = knmHG(couplings, q1, q2, surface_map=smap, method="romhom")
+print time.time() - t0
+
+print kbh
+print kr
+print krm
+
+plot_knm_matrix(couplings, np.log10(np.abs(kbh - kr)))
+
+print params
\ No newline at end of file

pykat/utilities/knm.py

 from itertools import combinations_with_replacement as combinations
 from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
 from pykat.exceptions import BasePyKatException
+from romhom import u_star_u
+from progressbar import ProgressBar, ETA, Percentage, Bar
+from scipy.interpolate import interp2d
+from scipy.integrate import dblquad
 
 import time
 import maps
@@ -11,10 +15,7 @@ import collections
 import math
 import cmath
 
-from romhom import u_star_u
-from progressbar import ProgressBar, ETA, Percentage, Bar
-
-def makeCouplingMatrix(max_order):
+def makeCouplingMatrix(max_order, Neven=True, Nodd=True, Meven=True, Modd=True):
     max_order = int(max_order)
     c = []
     for n in range(0, max_order+1):
@@ -30,13 +31,119 @@ def makeCouplingMatrix(max_order):
         for j in c:
             e = list(i)
             e.extend(j)
+            
+            if not Neven and (e[0]-e[2]) % 2 == 0: continue
+            if not Nodd and (e[0]-e[2]) % 2 == 1: continue
+            if not Meven and (e[1]-e[3]) % 2 == 0: continue
+            if not Modd and (e[1]-e[3]) % 2 == 1: continue
+            
             row.append(e)
         
-        M.append(row)
+        
+        M.append(np.array(row).squeeze())
     
     return np.array(M)
 
-def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0)):
+def adaptive_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, smap=None, delta=(0,0), params={}):
+    
+    if q1y == None:
+        q1y = q1
+        
+    if q2y == None:
+        q2y = q2
+    
+    if "epsabs" not in params: params["epsabs"] = 1e-6
+    if "epsrel" not in params: params["epsrel"] = 1e-6
+    if "usepolar" not in params: params["usepolar"] = False
+        
+    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])
+    
+    Nfuncs = []
+    Nfuncs.append(0)
+    
+    if smap != None:
+        
+        if not params["usepolar"]:
+            xlims = (min(smap.x), max(smap.x))
+            ylims = (min(smap.y), max(smap.y))
+            
+            def Rfunc(y,x):
+                Nfuncs[-1] += len(x)
+                return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real
+                
+            def Ifunc(y,x):
+                Nfuncs[-1] += len(x)
+                return (Hg_in.Unm(x+delta[0], y+delta[1]) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag
+            
+        else:
+            xlims = (0, 2*math.pi)
+            ylims = (0, params["aperture"])
+            
+            def Rfunc(r, phi):
+                Nfuncs[-1] += len(x)
+                x = r*np.cos(phi)
+                y = r*np.sin(phi)
+                return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).real
+                
+            def Ifunc(r, phi):
+                Nfuncs[-1] += len(x)
+                x = r*np.cos(phi)
+                y = r*np.sin(phi)
+                return (r * Hg_in.Unm(x, y) * smap.z_xy(x=x,y=y) * Hg_out.Unm(x, y).conjugate()).imag
+            
+    else:
+        if not params["usepolar"]:
+            _x = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1.w
+            _y = 4 * math.sqrt(1+max(mode_in[0],mode_in[1])) * q1y.w
+        
+            xlims = (-_x, _x)
+            ylims = (-_y, _y)
+        
+            def Rfunc(y, x):
+                Nfuncs[-1] += len(r)
+                return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).real
+                
+            def Ifunc(y,x):
+                Nfuncs[-1] += len(r)
+                return (Hg_in.Unm(x+delta[0], y+delta[1]) * Hg_out.Unm(x, y).conjugate()).imag
+        else:
+            xlims = (0, 2*math.pi)
+            ylims = (0, params["aperture"])
+            
+            def Rfunc(r, phi):
+                
+                if hasattr(r, "__len__"):
+                    Nfuncs[-1] += len(r)
+                else:
+                    Nfuncs[-1] += 1
+                    
+                x = r*np.cos(phi)
+                y = r*np.sin(phi)
+                return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).real
+                
+            def Ifunc(r, phi):
+                if hasattr(r, "__len__"):
+                    Nfuncs[-1] += len(r)
+                else:
+                    Nfuncs[-1] += 1
+                    
+                x = r*np.cos(phi)
+                y = r*np.sin(phi)
+                return (r * Hg_in.Unm(x, y) * Hg_out.Unm(x, y).conjugate()).imag
+    
+    R, errR = dblquad(Rfunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"])
+    I, errI = dblquad(Ifunc, xlims[0], xlims[1], lambda y: ylims[0], lambda y: ylims[1], epsabs=params["epsabs"], epsrel=params["epsrel"])
+    
+    params["Nfuncs"] = Nfuncs[0]
+    params["errors"] = (errR, errI)
+    
+    return R + 1j * I
+    
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0), params={}):
 
     if Axy == None:
         Axy == np.ones((len(x), len(y)))
@@ -71,7 +178,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, delta=(0,0)):
+def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0), params={}):
     if q1y == None:
         q1y = q1
         
@@ -311,7 +418,9 @@ def bayerhelms_HG_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, gamma=(0,0)
 
     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)):
+
+
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0), params={}):
     if q1y == None:
         q1y = q1
         
@@ -339,7 +448,7 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
         __fac_cache.append(math.factorial(n))
     
     if surface_map != None:  
-        Axy = surface_map.z_xy(q1.wavelength)
+        Axy = surface_map.z_xy(wavelength=q1.wavelength)
     
         x = surface_map.x
         y = surface_map.y
@@ -388,12 +497,15 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
             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, 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)
+            elif method == "adaptive":
+                K[i] = adaptive_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, smap=surface_map, delta=delta, params=params)
             else:
                 raise BasePyKatException("method value '%s' not accepted" % method)
             

pykat/utilities/maps.py

 from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
+from scipy.interpolate import interp2d
 import numpy as np
 import math
         
@@ -11,6 +12,7 @@ class surfacemap(object):
         self.center = center
         self.step_size = step_size
         self.scaling = scaling
+        self.__interp = None
         
         if data == None:
             self.data = np.zeros(size)
@@ -36,6 +38,42 @@ class surfacemap(object):
                     mapfile.write("%.15g " % self.data[i,j])
                 mapfile.write("\n")
     
+    @property
+    def data(self):
+        return self.__data
+    
+    @data.setter
+    def data(self, value):
+        self.__data = value
+        self.__interp = None
+    
+    @property
+    def center(self):
+        return self.__center
+    
+    @center.setter
+    def center(self, value):
+        self.__center = value
+        self.__interp = None
+    
+    @property
+    def step_size(self):
+        return self.__step_size
+    
+    @step_size.setter
+    def step_size(self, value):
+        self.__step_size = value
+        self.__interp = None
+
+    @property
+    def scaling(self):
+        return self.__scaling
+    
+    @scaling.setter
+    def scaling(self, value):
+        self.__scaling = value
+        self.__interp = None
+
     @property
     def x(self):
         return self.step_size[0] * (np.array(range(1, self.data.shape[0]+1)) - self.center[0])
@@ -56,22 +94,34 @@ class surfacemap(object):
     def ROMWeights(self):
         return self._rom_weights
     
-    def z_xy(self, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0):
+    def z_xy(self, x=None, y=None, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0):
+        
+        if x == None and y == None:
+            data = self.scaling * self.data
+        else:
+            if self.__interp == None:
+                self.__interp = interp2d(self.x, self.y, self.data * self.scaling)
+                
+            data = self.__interp(x, y)
             
         if direction == "reflection":
             if "phase" in self.type:
                 k = math.pi * 2 / wavelength
-                return np.exp(2j * k * self.scaling * self.data)
+                return np.exp(2j * k * data)
+                
             elif "absorption" in self.type:
-                return np.sqrt(1.0 - self.scaling * self.data)
+                return np.sqrt(1.0 - data)
+                
             else:
                 raise BasePyKatException("Map type needs handling")
         elif direction == "transmission":
             if "phase" in self.type:
                 k = math.pi * 2 / wavelength
-                return np.exp((nr1-nr2)*k * self.scaling * self.data)
+                return np.exp((nr1-nr2)*k * data)
+                
             elif "absorption" in self.type:
-                return np.sqrt(1.0 - self.scaling * self.data)
+                return np.sqrt(1.0 - data)
+                
             else:
                 raise BasePyKatException("Map type needs handling")
         else:

pykat/utilities/optics/gaussian_beams.py

@@ -314,12 +314,12 @@ class HG_beam(object):
             
     def _calc_constants(self):
         self.__xpre_const = math.pow(2.0/math.pi, 0.25)
-        self.__xpre_const *= math.sqrt(1.0/(self._qx.w*2**self._n * math.factorial(self._n)))
+        self.__xpre_const *= math.sqrt(1.0/(self._qx.w0 * 2**(self._n) * math.factorial(self._n)))
         self.__xpre_const *= cmath.sqrt(1j*self._qx.imag / self._qx.q)
         self.__xpre_const *= ((1j*self._qx.imag * self._qx.q.conjugate())/(-1j*self._qx.imag * self._qx.q)) ** ( self._n/2.0)
         
         self.__ypre_const = math.pow(2.0/math.pi, 0.25)
-        self.__ypre_const *= math.sqrt(1.0/(self._qy.w*2**self._m * math.factorial(self._m)))
+        self.__ypre_const *= math.sqrt(1.0/(self._qy.w0 * 2**(self._m) * math.factorial(self._m)))
         self.__ypre_const *= cmath.sqrt(1j*self._qy.imag / self._qy.q)
         self.__ypre_const *= ((1j*self._qy.imag * self._qy.q.conjugate())/(-1j*self._qy.imag * self._qy.q)) **(self._m/2.0)