Commit 3a418d7c by Daniel Brown

### more romhom and map fixes

parent 16c0700e
 from pykat.utilities.maps import read_map from pykat.utilities.maps import read_map, tiltmap, surfacemap from pykat.utilities.knm import * from pykat.utilities.knm import * from pykat.utilities.optics.gaussian_beams import beam_param from pykat.utilities.optics.gaussian_beams import beam_param import time import time ... @@ -7,20 +7,23 @@ couplings = np.array([[[0,0,0,0], [0,0,1,0], [0,0,0,1]], ... @@ -7,20 +7,23 @@ 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]], [[1,0,0,0], [1,0,1,0], [1,0,0,1]], [[0,1,0,0], [0,1,1,0], [0,1,0,1]]]) [[0,1,0,0], [0,1,1,0], [0,1,0,1]]]) q1 = beam_param(w0=4e-2, z=0) q1 = beam_param(w0=3e-2, z=0) q2 = beam_param(w0=4e-2, z=0) q2 = beam_param(w0=3e-2, z=0) m = read_map('etm08_virtual.txt') size = np.array([200,200]) stepsize = 0.2/(size-1) m.data = np.zeros((100, 100)) # This map has points evenly spaced about the axes m.x = np.linspace(-0.15, 0.15, m.size[0]) m = tiltmap("tilt", size, stepsize, (1e-6, 0)) m.y = np.linspace(-0.15, 0.15, m.size[1]) # Shifting the central point changes the results slightly m.center = np.array(m.size)/2 # but this is about a 1e-7 change, probably due to extra m.step_size = (m.x[1]-m.x[0], m.y[1]-m.y[0]) # clipping m.center += 20 print "generating weights..." print "Generating weights..." t0 = time.time() t0 = time.time() w, EI = m.generateROMWeights(isModeMatched=True) w, EI = m.generateROMWeights(isModeMatched=True) #w.writeToFile("testWeights.rom") print "Completed in ", time.time()-t0 print "Completed in ", time.time()-t0 print "computing Riemann knm..." print "computing Riemann knm..." ... @@ -39,7 +42,7 @@ print np.abs(K2) ... @@ -39,7 +42,7 @@ print np.abs(K2) print "Speed up", tr/tt print "Speed up", tr/tt #w.writeToFile("testWeights.rom") ... ...
 ... @@ -94,47 +94,64 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non ... @@ -94,47 +94,64 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non npr = mode_in[1] npr = mode_in[1] mpr = mode_out[1] mpr = mode_out[1] # if cache == None: foundSymmetry = np.all(weights.EI["xp"].nodes == -weights.EI["xm"].nodes) and np.all(weights.EI["yp"].nodes == -weights.EI["ym"].nodes) 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) if foundSymmetry: u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes-0.1e-2) if cache == None: u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes-0.1e-2) 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(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["ym"].nodes) # else: w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3 # strx = "x[%i,%i]" % (mode_in[0], mode_out[0]) w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4 # stry = "y[%i,%i]" % (mode_in[1], mode_out[1]) w_ij_Q1Q2 = weights.w_ij_Q1 + weights.w_ij_Q2 # w_ij_Q1Q4 = weights.w_ij_Q1 + weights.w_ij_Q4 # u_x_nodes = cache[strx] w_ij_Q2Q3 = weights.w_ij_Q2 + weights.w_ij_Q3 # u_y_nodes = cache[stry] w_ij_Q3Q4 = weights.w_ij_Q3 + weights.w_ij_Q4 w_ij_Q1Q2Q3Q4 = weights.w_ij_Q1 + weights.w_ij_Q2 + weights.w_ij_Q3 + weights.w_ij_Q4 u_xy_Q1 = np.outer(u_xm_nodes, u_yp_nodes) else: u_xy_Q2 = np.outer(u_xp_nodes, u_yp_nodes) strx = "x[%i,%i]" % (mode_in[0], mode_out[0]) u_xy_Q3 = np.outer(u_xp_nodes, u_ym_nodes) stry = "y[%i,%i]" % (mode_in[1], mode_out[1]) u_xy_Q4 = np.outer(u_xm_nodes, u_ym_nodes) u_x_nodes = cache[strx] # n_mod_2 = n % 2 u_y_nodes = cache[stry] # m_mod_2 = m % 2 # npr_mod_2 = npr % 2 u_xy_nodes = np.outer(u_x_nodes, u_y_nodes) # mpr_mod_2 = mpr % 2 # n_mod_2 = n % 2 # if n_mod_2 == m_mod_2 and npr_mod_2 == mpr_mod_2: m_mod_2 = m % 2 # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4) npr_mod_2 = npr % 2 # mpr_mod_2 = mpr % 2 # elif n_mod_2 != m_mod_2: # if npr_mod_2 == mpr_mod_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_Q1Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q2Q3) k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4) # else: # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3) elif n_mod_2 != m_mod_2: # if npr_mod_2 == mpr_mod_2: # elif 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) # if n_mod_2 == m_mod_2: else: # k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q2) k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q2Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q3) # 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: k_ROQ = np.einsum('ij,ij', u_xy_Q1, weights.w_ij_Q1) if n_mod_2 == m_mod_2: k_ROQ += np.einsum('ij,ij', u_xy_Q2, weights.w_ij_Q2) k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q3Q4) + np.einsum('ij,ij', -u_xy_nodes, w_ij_Q1Q2) k_ROQ += np.einsum('ij,ij', u_xy_Q3, weights.w_ij_Q3) else: k_ROQ += np.einsum('ij,ij', u_xy_Q4, weights.w_ij_Q4) 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 return k_ROQ ... ...
 from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights import numpy import numpy as np import math 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: class surfacemap(object): def __init__(self, name, maptype, size, center, step_size, scaling, data=None): def __init__(self, name, maptype, size, center, step_size, scaling, data=None): self.name = name self.name = name ... @@ -26,7 +11,12 @@ class surfacemap: ... @@ -26,7 +11,12 @@ class surfacemap: self.center = center self.center = center self.step_size = step_size self.step_size = step_size self.scaling = scaling self.scaling = scaling self.data = data if data == None: self.data = np.zeros(size) else: self.data = data self._rom_weights = None self._rom_weights = None def write_map(self, filename): def write_map(self, filename): ... @@ -48,11 +38,11 @@ class surfacemap: ... @@ -48,11 +38,11 @@ class surfacemap: @property @property def x(self): def x(self): return self.step_size[0] * (numpy.array(range(1, self.data.shape[0]+1)) - self.center[0]) return self.step_size[0] * (np.array(range(1, self.data.shape[0]+1)) - self.center[0]) @property @property def y(self): def y(self): return self.step_size[1] * (numpy.array(range(1, self.data.shape[1]+1))- self.center[1]) return self.step_size[1] * (np.array(range(1, self.data.shape[1]+1))- self.center[1]) @property @property def size(self): def size(self): ... @@ -60,7 +50,7 @@ class surfacemap: ... @@ -60,7 +50,7 @@ class surfacemap: @property @property def offset(self): def offset(self): return numpy.array(self.step_size)*(self.center - numpy.array(self.size)/2) return np.array(self.step_size)*(self.center - np.array(self.size)/2) @property @property def ROMWeights(self): def ROMWeights(self): ... @@ -70,7 +60,7 @@ class surfacemap: ... @@ -70,7 +60,7 @@ class surfacemap: if "phase" in self.type: if "phase" in self.type: k = math.pi * 2 / wavelength k = math.pi * 2 / wavelength return numpy.exp(2j * k * self.scaling * self.data) return np.exp(2j * k * self.scaling * self.data) else: else: raise BasePyKatException("Map type needs handling") raise BasePyKatException("Map type needs handling") ... @@ -87,7 +77,7 @@ class surfacemap: ... @@ -87,7 +77,7 @@ class surfacemap: if len(xp) > 0: EI["xp"] = makeEmpiricalInterpolant(makeReducedBasis(xp, 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(ym) > 0: EI["ym"] = makeEmpiricalInterpolant(makeReducedBasis(ym, isModeMatched=isModeMatched)) if len(yp) > 0: EI["yp"] = makeEmpiricalInterpolant(makeReducedBasis(yp, 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(x0) > 0: EI["x0"] = makeEmpiricalInterpolant(makeReducedBasis(x0, isModeMatched=isModeMatched)) #if len(y0) > 0: EI["y0"] = makeEmpiricalInterpolant(makeReducedBasis(y0, isModeMatched=isModeMatched)) #if len(y0) > 0: EI["y0"] = makeEmpiricalInterpolant(makeReducedBasis(y0, isModeMatched=isModeMatched)) ... @@ -122,8 +112,27 @@ class surfacemap: ... @@ -122,8 +112,27 @@ class surfacemap: return fig return fig class tiltmap(surfacemap): def __init__(self, name, size, step_size, tilt): surfacemap.__init__(self, name, "phase", size, (np.array(size)+1)/2.0, step_size, 1) self.tilt = tilt @property def tilt(self): return self.__tilt @tilt.setter def tilt(self, value): self.__tilt = value xx, yy = np.meshgrid(self.x, self.y) self.data = xx * self.tilt[0] + yy * self.tilt[1] def read_map(filename): def read_map(filename): with open(filename, 'r') as f: with open(filename, 'r') as f: ... @@ -137,7 +146,7 @@ def read_map(filename): ... @@ -137,7 +146,7 @@ def read_map(filename): data = numpy.loadtxt(filename, dtype=numpy.float64,ndmin=2,comments='%') data = np.loadtxt(filename, dtype=np.float64,ndmin=2,comments='%') return surfacemap(name,maptype,size,center,step,scaling,data) return surfacemap(name,maptype,size,center,step,scaling,data) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!