fixing more python 3 bugs with romhom adding in merged map objecting, for merging surface maps

bin/test_2d_plot.py

 from pykat import finesse
-from pykat.utilities.optics.gaussian_beams import gauss_param
 from pykat.detectors import *
 from pykat.components import *
 from pykat.commands import *
@@ -26,7 +25,7 @@ x2axis b1 y lin -6 6 100
 kat = finesse.kat()
 
 kat.parseCommands(code)
-kat.s1.n1.q = gauss_param(w0=1e-3, z=0)
+kat.s1.n1.q = pykat.beam_param(w0=1e-3, z=0)
 
 out = kat.run(printout=0,printerr=0)
 

pykat/finesse.py

@@ -249,7 +249,7 @@ class katRun(object):
       
 class katRun2D(object):
     def __init__(self):
-        self.runtime
+        self.runtime = None
         self.startDateTime = datetime.datetime.now()
         self.x = None
         self.y = None
@@ -1692,10 +1692,13 @@ class kat(object):
         if self.noxaxis == True:
             out.append("noxaxis\n")
             
-        if self.yaxis != None:
+        if self.yaxis is not None:
             out.append("yaxis {0}\n".format(self.yaxis))
          
-        if self.printmatrix != None and self.printmatrix == True:
+        if self.retrace is not None:
+            out.append("retrace {0}\n".format(str(self.retrace)))
+            
+        if self.printmatrix is not None and self.printmatrix == True:
             out.append("printmatrix\n")
             
         if self.lambda0 != 1064e-9:

pykat/maths/hermite.py

@@ -10,36 +10,36 @@ def hermite(n, x):
     elif n == 3:		
         return (8.0 * x * x * x - 12.0 * x) 
     elif n == 4:
-    	x_sq = x * x	
+        x_sq = x * x	
         return (16.0 * x_sq*x_sq - 48.0 * x_sq + 12.0) 
     elif n == 5:	
-    	x_sq = x * x 
-    	x_cb = x * x_sq
+        x_sq = x * x 
+        x_cb = x * x_sq
         return (32.0 * x_cb * x_sq - 160.0 * x_cb + 120.0 * x)
     elif n == 6:
-    	x_sq = x * x	
-    	x_four = x_sq * x_sq
+        x_sq = x * x	
+        x_four = x_sq * x_sq
         return (64.0 * x_four*x_sq - 480.0 * x_four + 720.0 * x_sq - 120.0) 
     elif n == 7:
         x_sq = x*x
-    	x_cb = x_sq*x
-    	x_four = x_cb*x
+        x_cb = x_sq*x
+        x_four = x_cb*x
         return (128.0 * x_cb*x_four - 1344.0 * x_cb*x_sq + 3360.0 * x_cb - 1680.0 * x) 
     elif n == 8:
-    	x_sq = x*x
-    	x_four = x_sq*x_sq
-    	x_six = x_four*x_sq	
+        x_sq = x*x
+        x_four = x_sq*x_sq
+        x_six = x_four*x_sq	
         return (256.0 * x_four*x_four - 3584.0 * x_six + 13440.0 * x_four - 13440.0 * x_sq + 1680.0) 
     elif n == 9:
-    	x_sq = x*x
-    	x_cb = x_sq*x
-    	x_four = x_cb*x	
+        x_sq = x*x
+        x_cb = x_sq*x
+        x_four = x_cb*x	
         return (512.0 * x_cb*x_cb*x_cb - 9216.0 * x_four*x_cb + 48384.0 * x_cb*x_sq - 80640.0 * x_cb + 30240.0 * x) 
     elif n == 10:
-    	x_sq = x*x
-    	x_cb = x_sq*x
-    	x_four = x_cb*x
-    	x_five = x_four * x
+        x_sq = x*x
+        x_cb = x_sq*x
+        x_four = x_cb*x
+        x_five = x_four * x
         return (1024.0 * x_five*x_five - 23040.0 * x_four*x_four + 161280.0 * x_cb*x_cb - 403200.0 * x_four + 302400.0 * x_sq - 30240.0) 
     elif n == 11:
         return -665280 * x + 2217600 * x**3 - 1774080 * x**5 + 506880 * x**7 - 56320 * x**9 + 2048 * x**11

pykat/optics/maps.py

@@ -303,7 +303,182 @@ class surfacemap(object):
             pylab.show()
         
         return fig
+
+class mergedmap:
+    """
+    A merged map combines multiple surfaces map to form one. Such a map can be used
+    for computations of coupling coefficients but it cannot be written to a file to 
+    be used with Finesse. For this you must output each map separately.
+    
+    """
+    
+    def __init__(self, name, size, center, step_size, scaling):
+        
+        self.name = name
+        self.center = center
+        self.step_size = step_size
+        self.scaling = scaling
+        self.__interp = None
+        self._rom_weights = None
+        self.__maps = []
+      
+    def addMap(self, m):
+        self.__maps.append(m)
+    
+    @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.size[0]+1)) - self.center[0])
+        
+    @property
+    def y(self):
+        return self.step_size[1] * (np.array(range(1, self.size[1]+1))- self.center[1])
+
+    @property
+    def size(self):
+        return self.__maps[0].data.shape
+            
+    @property
+    def offset(self):
+        return np.array(self.step_size)*(np.array(self.center) - 1/2. - np.array(self.size)/2.0)
+    
+    @property
+    def ROMWeights(self):
+        return self._rom_weights
+    
+    def z_xy(self, wavelength=1064e-9, direction="reflection", nr1=1.0, nr2=1.0):
+        
+        z_xy = np.ones(self.size, dtype=np.complex128)
+        
+        for m in self.__maps:
+            z_xy *= m.z_xy(wavelength=wavelength, direction=direction, nr1=nr1, nr2=nr2)
+            
+        return z_xy
+        
+    def generateROMWeights(self, EIxFilename, EIyFilename=None, verbose=False, interpolate=False):
+        if interpolate == True:
+            # Use EI nodes to interpolate if we
+            with open(EIxFilename, 'rb') as f:
+                EIx = pickle.load(f)
+
+            if EIyFilename is None:
+                EIy = EIx
+            else:
+                with open(EIyFilename, 'rb') as f:
+                    EIy = pickle.load(f)
+
+            x = EIx.x
+            x.sort()
+            nx = np.unique(np.hstack((x, -x[::-1])))
+        
+            y = EIy.x
+            y.sort()
+            ny = np.unique(np.hstack((y, -y[::-1])))
+            
+            self.interpolate(nx, ny)
+        
+        self._rom_weights = makeWeightsNew(self, EIxFilename, EIyFilename, verbose=verbose)
+        return self.ROMWeights
+
+    def interpolate(self, nx, ny, **kwargs):
+        """
+        Interpolates all the maps that are used to fc
+        
+        Uses scipy.interpolate.interp2d and any keywords arguments are
+        passed on to it, thus settings like interpolation type and
+        fill values can be set.
+        
+        The range of nx and ny must contain the value zero so that the
+        center point of the map can be set.
+        """
+
+        for m in self.__maps:
+            m.interpolate(nx, ny)
+
+    def plot(self, mode="absorption", show=True, clabel=None, xlim=None, ylim=None, wavelength=1064e-9):
+        
+        import pylab
+        
+        if xlim is not None:
+            _x = np.logical_and(self.x<=max(xlim)/100.0, self.x>=min(xlim)/100.0)
+            xmin = np.min(np.where(_x == True))
+            xmax = np.max(np.where(_x == True))
+        else:
+            xmin = 0
+            xmax = len(self.x)-1
+            xlim = [self.x.min()*100, self.x.max()*100]
+    
+        if ylim is not None:
+            _y = np.logical_and(self.y<=max(ylim)/100.0, self.y>=min(ylim)/100.0)
+            ymin = np.min(np.where(_y == True))
+            ymax = np.max(np.where(_y == True))
+        else:
+            ymin = 0
+            ymax = len(self.y)-1
+            ylim = [self.y.min()*100, self.y.max()*100]
+
+        if mode == "absorption":
+            # plots how much of field is absorbed
+            data = 1-np.abs(self.z_xy())
+        elif mode == "meter":
+            # plot the phase in terms of meters of displacement
+            k = 2*np.pi/wavelength
+            data = np.angle(self.z_xy()) / (2*k)
+            
+        zmin = data[xmin:xmax,ymin:ymax].min()
+        zmax = data[xmin:xmax,ymin:ymax].max()
+
+        # 100 factor for scaling to cm
+        xrange = 100*self.x
+        yrange = 100*self.y
+
+        fig = pylab.figure()
+        axes = pylab.pcolormesh(xrange, yrange, data, vmin=zmin, vmax=zmax)
+        pylab.xlabel('x [cm]')
+        pylab.ylabel('y [cm]')
+
+        if xlim is not None: pylab.xlim(xlim)
+        if ylim is not None: pylab.ylim(ylim)
+
+        pylab.title('Merged map {0}, mode {1}'.format(self.name, mode))
+
+        cbar = fig.colorbar(axes)
+        cbar.set_clim(zmin, zmax)
         
+        if clabel is not None:
+            cbar.set_label(clabel)
+    
+        if show:
+            pylab.show()
+        
+        return fig
+
 class aperturemap(surfacemap):
     
     def __init__(self, name, size, step_size, R):

pykat/optics/romhom.py

@@ -323,7 +323,7 @@ def makeEmpiricalInterpolant(RB, sort=False):
     nodes = np.array(x_nodes, dtype=np.float64)
     
     if sort:
-        print sort
+        print (sort)
         ix = np.argsort(nodes)
         nodes = nodes[ix]
         node_indices = node_indices[ix]
@@ -513,7 +513,7 @@ def _write_TS(queue, filename, tssize):
                 i += 1
                 
                 if i % 1000 == 0:
-                    print i/float(tssize)
+                    print(i/float(tssize))
                     hfile.flush()
     finally:
         hfile.close()
@@ -542,8 +542,8 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo
 
     nModes = 0
     
-    for n in xrange(0, maxOrder+1):
-        for m in xrange(0, maxOrder+1):
+    for n in range(0, maxOrder+1):
+        for m in range(0, maxOrder+1):
             if n+m <= maxOrder and n <= m:
                 nModes += 1
             
@@ -563,10 +563,10 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo
     
     hfile.close() # make sure it's closed before
     
-    print "Starting processes..."
+    print("Starting processes...")
     # Have a bunch of processes doing the computation and one doing the writing
     iprocesses = [Process(target=_compute_TS, name="irom%i" % i, args=(iq, oq, x, w)) for i in range(NProcesses)]
-    oprocess = Process(target=_write_TS, name="orom%i" % i, args=(oq, filename, tssize))
+    oprocess = Process(target=_write_TS, name="orom", args=(oq, filename, tssize))
     
     oprocess.start()
     
@@ -574,10 +574,10 @@ def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCo
         for P in iprocesses:
             P.start()
         
-        print "Filling queue..."
+        print("Filling queue...")
         curr = 0
-        for n in xrange(0, maxOrder+1):
-            for m in xrange(0, maxOrder+1):
+        for n in range(0, maxOrder+1):
+            for m in range(0, maxOrder+1):
                 if n+m <= maxOrder and n <= m:
                     for _z in z:
                         for _w0 in w0:
@@ -760,10 +760,10 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1
                 next_RB_index = max_idx[np.argmax(max_res)]
         
             if worst_error <= tol:
-                print "Final basis size = %d, Final error = %e, Tolerance=%e" % (k, worst_error, tol) 
+                print( "Final basis size = %d, Final error = %e, Tolerance=%e" % (k, worst_error, tol) )
                 break
 
-            print "worst error = %e at %i on iteration %d" % (worst_error, next_RB_index, k)			
+            print ("worst error = %e at %i on iteration %d" % (worst_error, next_RB_index, k))
         
             epsilon = TS['%i/data' % next_RB_index].value
             res   = epsilon - emp_interp(B, epsilon, EI_indices)
@@ -782,7 +782,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1
             _TS = TS[str(next_RB_index)]
             f.write("%15.15e %15.15e %i %i\n" % (_TS["w0"].value, _TS["z"].value, _TS["n1"].value, _TS["n2"].value))
                 
-    print time.time() - start, "Seconds"
+    print (time.time() - start, "Seconds")
     
     if NProcesses > 1:
         for P in procs: P.terminate()
@@ -791,7 +791,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1
 
     greedyFilenameBase = os.path.splitext(greedyFilename)[0]
     
-    print "Writing to %s" % greedyFilename
+    print ("Writing to %s" % greedyFilename)
                        
     EI = EmpiricalInterpolant(B=np.matrix(B).real,
                               nodes=np.array(x_nodes).squeeze(),
@@ -803,7 +803,7 @@ def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-1
         with open("%s.p" % EIFilename, 'wb') as f:
             pickle.dump(EI, f)
         
-        print "Writing to %s.p" % EIFilename
+        print ("Writing to %s.p" % EIFilename)
                               
     return EI