diff --git a/bin/test_romhom.py b/bin/test_romhom.py
index 6275e6fbaf06d884fb808e15de2bf1ea61ea9592..d48170610c5d993b3209946dc167991425b855b6 100644
--- a/bin/test_romhom.py
+++ b/bin/test_romhom.py
@@ -2,28 +2,27 @@ from pykat.utilities.maps import read_map, tiltmap, surfacemap
from pykat.utilities.knm import *
from pykat.utilities.optics.gaussian_beams import beam_param
import time
+import numpy as np
-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]]])
+couplings = makeCouplingMatrix(5)
q1 = beam_param(w0=3e-2, z=0)
q2 = beam_param(w0=3e-2, z=0)
-size = np.array([200,200])
-stepsize = 0.2/(size-1)
+size = np.array([1200, 1200])
+stepsize = 0.3/(size-1)
# This map has points evenly spaced about the axes
-m = tiltmap("tilt", size, stepsize, (1e-6, 0))
+m = tiltmap("tilt", size, stepsize, (1e-6, 1e-6))
# Shifting the central point changes the results slightly
# but this is about a 1e-7 change, probably due to extra
# clipping
-m.center += 20
+m.center += 0
print "Generating weights..."
t0 = time.time()
w, EI = m.generateROMWeights(isModeMatched=True)
-#w.writeToFile("testWeights.rom")
+w.writeToFile("testWeights.rom")
print "Completed in ", time.time()-t0
print "computing Riemann knm..."
@@ -31,15 +30,14 @@ t0 = time.time()
K1 = knmHG(couplings, q1, q2, surface_map=m)
tr = time.time()-t0
print "Completed in ", tr
-print np.abs(K1)
+#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 np.abs(K2)
print "Speed up", tr/tt
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
index 0ebb3eaafc783596771b1d600b37dd711654b11f..7d418c4696bb366a26aa2ae43b0c4392de706969 100644
--- a/pykat/utilities/knm.py
+++ b/pykat/utilities/knm.py
@@ -8,9 +8,33 @@ import numpy as np
import pykat
import collections
import math
+
from romhom import u_star_u
+from progressbar import ProgressBar, ETA, Percentage, Bar
+
+def makeCouplingMatrix(max_order):
+ max_order = int(max_order)
+ c = []
+ for n in range(0, max_order+1):
+ for m in range(0, max_order+1):
+ if n+m <= max_order:
+ c.append([n,m])
+
+ M = []
+
+ for i in c:
+ row = []
+
+ for j in c:
+ e = list(i)
+ e.extend(j)
+ row.append(e)
+
+ M.append(row)
+
+ return np.array(M)
-def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None):
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None):
if Axy == None:
Axy == np.ones((len(x), len(y)))
@@ -23,19 +47,62 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
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])
-
+
dx = abs(x[1] - x[0])
- dy = abs(y[1] - y[0])
+ dy = abs(y[1] - y[0])
+
+ if cache == 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)
- U2 = Hg_out.Unm(x,y).conjugate()
+ U1 = Hg_in.Unm(x,y)
+ U2 = Hg_out.Unm(x,y).conjugate()
+
+ return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+ else:
+
+ strx = "u1[%i,%i]" % (mode_in[0], mode_out[0])
+ stry = "u2[%i,%i]" % (mode_in[1], mode_out[1])
+
+ return dx * dy * np.einsum('ij,ij', Axy, np.outer(cache[strx], cache[stry]))
+
- return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
+
+def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
+ if q1y == None:
+ q1y = q1
+
+ if q2y == None:
+ q2y = q2
+
+ it = np.nditer(couplings, flags=['refs_ok','f_index'])
+
+ cache = {}
+
+ while not it.finished:
+ try:
+ mode_in = [int(it.next()), int(it.next())]
+ mode_out = [int(it.next()), int(it.next())]
+
+ 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])
+
+ if strx not in cache:
+ cache[strx] = 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()
+
+ except StopIteration:
+ break
+
+ return cache
+
def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
if q1y == None:
@@ -65,14 +132,10 @@ def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
stry = "y[%i,%i]" % (mode_in[1], mode_out[1])
if strx not in cache:
- cache[strx] = u_star_u(q1.z, q2.z, q1.w0, q2.w0, mode_in[0], mode_out[0], weights.nodes)
-
- if q1 == q1y and q2 == q2y:
- cache[stry] = cache[strx]
- elif stry not in cache:
- cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], weights.nodes)
-
+ cache[strx] = u_star_u(q1.z, q2.z, q1.w0, q2.w0, mode_in[0], mode_out[0], weights.EI["xm"].nodes)
+ if stry not in cache:
+ cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], weights.EI["ym"].nodes)
except StopIteration:
break
@@ -94,12 +157,14 @@ 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)
+ 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(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["ym"].nodes)
+ 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)
+
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
@@ -107,12 +172,22 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
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_Q2 + weights.w_ij_Q3 + 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)
@@ -157,7 +232,7 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
-def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann"):
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=True):
if q1y == None:
q1y = q1
@@ -183,7 +258,7 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
it = np.nditer(couplings, flags=['refs_ok','f_index'])
i = 0
-
+
if method == "romhom":
if surface_map == None:
raise BasePyKatException("Using 'romhom' method requires a surface map to be specified")
@@ -192,28 +267,55 @@ 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)
- romcache = None#__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)
else:
- romcache = None
+ cache = None
weights = None
-
+
+ if verbose:
+ p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])
+
while not it.finished:
try:
-
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)
+ K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache)
elif method == "romhom":
- K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=romcache)
+ K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache)
else:
raise BasePyKatException("method value '%s' not accepted" % method)
-
+
i +=1
+
+ if verbose:
+ p.update(i*4)
+
except StopIteration:
break
- return K.reshape(couplings.shape[:-1])
\ No newline at end of file
+ return K.reshape(couplings.shape[:-1])
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
index 30f9e4bd9056b0176fc7cfd78d77064a6c5fd97e..e134a7b502e1ba6cdea605b8b01ddbafa75d0861 100644
--- a/pykat/utilities/romhom.py
+++ b/pykat/utilities/romhom.py
@@ -4,6 +4,7 @@ import os.path
import pykat
import collections
+from progressbar import ProgressBar, ETA, Percentage, Bar
from itertools import combinations_with_replacement as combinations
from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
from scipy.linalg import inv
@@ -150,21 +151,22 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
for i in range(len(params)):
- q1 = beam_param(w0=params[i][0], z=params[i][2])
+ q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
if isModeMatched:
q2 = q1
n = int(params[i][2])
m = int(params[i][3])
- w0_1 = params[i][0]
- w0_2 = w0_1
- re_q1 = params[i][1]
- re_q2 = re_q1
else:
- q2 = beam_param(w0=params[i][1], z=params[i][3])
+ q2 = beam_param(w0=float(params[i][1]), z=float(params[i][3]))
n = int(params[i][4])
m = int(params[i][5])
+ w0_1 = q1.w0
+ w0_2 = q2.w0
+ re_q1 = q1.z
+ re_q2 = q2.z
+
TS[IDx] = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, x)
# normalize
@@ -248,11 +250,8 @@ 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):
-
# get full A_xy
A_xy = smap.z_xy()