diff --git a/bin/test_romhom.py b/bin/test_romhom.py
index 2482c81ff50e091d59fd5ca260cc5fefb3638a33..0b19cab2525613a5d7b80f5abd3c453e14c173e9 100644
--- a/bin/test_romhom.py
+++ b/bin/test_romhom.py
@@ -1,33 +1,44 @@
from pykat.utilities.maps import read_map
from pykat.utilities.knm import *
from pykat.utilities.optics.gaussian_beams import beam_param
+import time
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)
+q1 = beam_param(w0=4e-2, z=0)
+q2 = beam_param(w0=4e-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.data = np.zeros((100, 100))
+m.x = np.linspace(-0.15, 0.15, m.size[0])
+m.y = np.linspace(-0.15, 0.15, 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()
+t0 = time.time()
+w, EI = m.generateROMWeights(isModeMatched=True)
+print "Completed in ", time.time()-t0
print "computing Riemann knm..."
+t0 = time.time()
K1 = knmHG(couplings, q1, q2, surface_map=m)
+tr = time.time()-t0
+print "Completed in ", tr
print np.abs(K1)
print "computing ROMHOM knm..."
+t0 = time.time()
K2 = knmHG(couplings, q1, q2, surface_map=m, method="romhom")
+tt = time.time()-t0
+print "Completed in ", tt
print np.abs(K2)
+print "Speed up", tr/tt
+
#w.writeToFile("testWeights.rom")
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
index f878050f59e0aa136ad455faa4af44956cb942a3..2303d6babf6c75be8fa3da8c3eab225ba2e44e61 100644
--- a/pykat/utilities/knm.py
+++ b/pykat/utilities/knm.py
@@ -93,55 +93,49 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
# 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)
+ # if cache == None:
+ 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)
+ u_ym_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes-0.1e-2)
+ u_yp_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["yp"].nodes-0.1e-2)
+
+ # 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]
+
+ 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)
+
+ # 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)
+ 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
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index e115f4c58e2c53eefb3c19be33faed8747f748f2..d19b032788b842971bf1d91e4979c988d13ca6cd 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -2,6 +2,22 @@ from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, m
import numpy
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:
def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
@@ -58,10 +74,25 @@ class surfacemap:
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)
+
+ def generateROMWeights(self, isModeMatched=True):
+ xp = self.x[self.x>0]
+ yp = self.y[self.y>0]
+ 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)
return self.ROMWeights, EI
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
index a4677b2781a74d40d0322a7596cb62176f1a90fe..30f9e4bd9056b0176fc7cfd78d77064a6c5fd97e 100644
--- a/pykat/utilities/romhom.py
+++ b/pykat/utilities/romhom.py
@@ -18,12 +18,13 @@ 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, nodes, limits):
+ def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, EI, 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.EI = EI
self.limits = limits
def writeToFile(self, filename):
@@ -111,7 +112,7 @@ 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 makeReducedBasis(x, offset=np.array([0,0]), isModeMatched=True, tolerance = 1e-12, sigma = 1):
+def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
if isModeMatched:
greedypts = 'matched20.txt'
@@ -246,52 +247,107 @@ def makeEmpiricalInterpolant(RB):
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)
+
+from progressbar import ProgressBar, ETA, Percentage, Bar
+
+def makeWeights(smap, EI, verbose=True):
-def makeWeights(smap, EI):
# get full A_xy
A_xy = smap.z_xy()
- hx = len(smap.x)/2
- hy = len(smap.y)/2
- fx = hx*2
- fy = hy*2
+ xm = smap.x[smap.x < 0]
+ xp = smap.x[smap.x > 0]
+
+ ym = smap.y[smap.y < 0]
+ yp = smap.y[smap.y > 0]
+
+ Q1xy = np.ix_(smap.x < 0, smap.y > 0)
+ Q2xy = np.ix_(smap.x > 0, smap.y > 0)
+ Q3xy = np.ix_(smap.x > 0, smap.y < 0)
+ Q4xy = np.ix_(smap.x < 0, smap.y < 0)
# get A_xy in the four quadrants of the x-y plane
- A_xy_Q1 = A_xy[0:hx][0:hy, hy:fy]
- A_xy_Q2 = A_xy[hx:fx][0:hy, hy:fy]
- A_xy_Q3 = A_xy[hx:fx][0:hy, 0:hy]
- A_xy_Q4 = A_xy[0:hx][0:hy, 0:hy]
-
+ A_xy_Q1 = A_xy[Q1xy]
+ A_xy_Q2 = A_xy[Q2xy]
+ A_xy_Q3 = A_xy[Q3xy]
+ A_xy_Q4 = A_xy[Q4xy]
+
+ # Not sure if there is a neater way to select the x=0,y=0 cross elements
+ A_xy_0 = np.hstack([A_xy[:,smap.x==0].flatten(), A_xy[smap.y==0,:].flatten()])
+
full_x = smap.x
full_y = smap.y
dx = full_x[1] - full_x[0]
dy = full_y[1] - full_y[0]
-
- B = EI.B
- nodes_nv = EI.nodes
- nodes_pv = - EI.nodes
+ 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)
+ p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
+
+ n = 0
+
# make integration weights
- w_ij_Q1 = np.zeros((len(B), len(B)), dtype = complex)
- w_ij_Q2 = np.zeros((len(B), len(B)), dtype = complex)
- w_ij_Q3 = np.zeros((len(B), len(B)), dtype = complex)
- w_ij_Q4 = np.zeros((len(B), len(B)), dtype = complex)
-
- for i in range(len(B)):
- for j in range(len(B)):
- B_ij_Q1 = np.outer(B[i], B[j][::-1])
+ Bx = EI["xm"].B
+ By = EI["yp"].B
+ w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q1 = np.outer(Bx[i], By[j])
w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)
- B_ij_Q2 = np.outer(B[i][::-1], B[j][::-1])
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["xp"].B
+ By = EI["yp"].B
+ w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q2 = np.outer(Bx[i], By[j])
w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
- B_ij_Q3 = np.outer(B[i][::-1], B[j])
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["xp"].B
+ By = EI["ym"].B
+ w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q3 = np.outer(Bx[i], By[j])
w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
- B_ij_Q4 = np.outer(B[i], B[j])
+ if verbose:
+ p.update(n)
+ n+=1
+
+ Bx = EI["xm"].B
+ By = EI["ym"].B
+ w_ij_Q4 = np.zeros((len(Bx),len(By)), dtype = complex)
+
+ for i in range(len(Bx)):
+ for j in range(len(By)):
+ B_ij_Q4 = np.outer(Bx[i], By[j])
w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
+
+ if verbose:
+ p.update(n)
+ n+=1
+ if verbose:
+ p.finish()
+
+ 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, EI=EI, limits=EI["xm"].limits)
+
+
- 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
diff --git a/requirements.txt b/requirements.txt
index fd2476d4abd3185e6ada7f0ef6a85c64620faf27..1ec8c6d2b728c9c4d3e65d3c1c94a0caa6c7ffc6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,2 +1,3 @@
numpy >= 1.6.2
-flask >= 0.10.1
\ No newline at end of file
+flask >= 0.10.1
+progressbar >= 2.2