diff --git a/bin/test_romhom.py b/bin/test_romhom.py
new file mode 100644
index 0000000000000000000000000000000000000000..2482c81ff50e091d59fd5ca260cc5fefb3638a33
--- /dev/null
+++ b/bin/test_romhom.py
@@ -0,0 +1,35 @@
+from pykat.utilities.maps import read_map
+from pykat.utilities.knm import *
+from pykat.utilities.optics.gaussian_beams import beam_param
+
+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)
+
+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.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()
+
+print "computing Riemann knm..."
+K1 = knmHG(couplings, q1, q2, surface_map=m)
+print np.abs(K1)
+
+print "computing ROMHOM knm..."
+K2 = knmHG(couplings, q1, q2, surface_map=m, method="romhom")
+print np.abs(K2)
+
+#w.writeToFile("testWeights.rom")
+
+
+
+
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
index 53ca9988938eae101f4d8e3f1f3638d1c4ee7805..f878050f59e0aa136ad455faa4af44956cb942a3 100644
--- a/pykat/utilities/knm.py
+++ b/pykat/utilities/knm.py
@@ -8,10 +8,7 @@ import numpy as np
import pykat
import collections
import math
-
-def makeCouplingMatrix(max_order):
- pass
-
+from romhom import u_star_u
def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None):
@@ -39,11 +36,124 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None
return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
-def ROM_HG_knm(m, mode_in, mode_out, q1, q2, q1y=None, q2y=None):
- pass
+def __gen_ROM_HG_knm_cache(weights, 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 = {}
+
+ cache["w_ij_Q1Q3"] = weights.w_ij_Q1 + weights.w_ij_Q3
+ cache["w_ij_Q2Q4"] = weights.w_ij_Q2 + weights.w_ij_Q4
+ cache["w_ij_Q1Q2"] = weights.w_ij_Q1 + weights.w_ij_Q2
+ cache["w_ij_Q1Q4"] = weights.w_ij_Q1 + weights.w_ij_Q4
+ cache["w_ij_Q2Q3"] = weights.w_ij_Q2 + weights.w_ij_Q3
+ cache["w_ij_Q3Q4"] = weights.w_ij_Q3 + weights.w_ij_Q4
+ cache["w_ij_Q1Q2Q3Q4"] = weights.w_ij_Q1 + weights.w_ij_Q3 + weights.w_ij_Q2 + weights.w_ij_Q4
+ while not it.finished:
+ try:
+ mode_in = [int(it.next()), int(it.next())]
+ mode_out = [int(it.next()), int(it.next())]
+
+ strx = "x[%i,%i]" % (mode_in[0], mode_out[0])
+ 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)
+
+
+
+ except StopIteration:
+ break
+
+ return cache
+
+def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=None):
+ if q1y == None:
+ q1y = q1
-def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
+ if q2y == None:
+ q2y = q2
+
+ # x modes
+ n = mode_in[0]
+ m = mode_out[0]
+
+ # 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)
+
+ return k_ROQ
+
+
+
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann"):
+ if q1y == None:
+ q1y = q1
+
+ if q2y == None:
+ q2y = q2
+
+ assert q1.wavelength == q2.wavelength and q1y.wavelength == q2y.wavelength and q1y.wavelength == q1.wavelength
couplings = np.array(couplings)
@@ -52,8 +162,7 @@ def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
if int(a) - a != 0:
raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")
- if "phase" in surface_map.type:
- Axy = np.exp(2j * 2 * math.pi / q1.wavelength * surface_map.data * surface_map.scaling)
+ Axy = surface_map.z_xy(q1.wavelength)
x = surface_map.x
y = surface_map.y
@@ -64,16 +173,32 @@ def knmHG(couplings, surface_map, q1, q2, q1y=None, q2y=None, method="riemann"):
i = 0
+ if method == "romhom":
+ if surface_map == None:
+ raise BasePyKatException("Using 'romhom' method requires a surface map to be specified")
+
+ weights = surface_map.ROMWeights
+
+ if weights == None:
+ raise BasePyKatException("The ROM weights need to be generated for this map before use.")
+
+ romcache = None#__gen_ROM_HG_knm_cache(weights, couplings, q1=q1, q2=q2, q1y=q1y, q2y=q2y)
+ else:
+ romcache = None
+ weights = None
+
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)
+ elif method == "romhom":
+ K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=romcache)
else:
- raise BasePyKatException("method value not accepted")
+ raise BasePyKatException("method value '%s' not accepted" % method)
i +=1
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index e9f1ab23d19591cf67cfce8ee18f89af990fd370..e115f4c58e2c53eefb3c19be33faed8747f748f2 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -1,5 +1,6 @@
from pykat.utilities.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
import numpy
+import math
class surfacemap:
def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
@@ -7,7 +8,6 @@ class surfacemap:
self.name = name
self.type = maptype
self.center = center
- self.size = size
self.step_size = step_size
self.scaling = scaling
self.data = data
@@ -19,7 +19,7 @@ class surfacemap:
mapfile.write("% Surface map\n")
mapfile.write("% Name: {0}\n".format(self.name))
mapfile.write("% Type: {0}\n".format(self.type))
- mapfile.write("% Size: {0} {1}\n".format(self.size[0], self.size[1]))
+ mapfile.write("% Size: {0} {1}\n".format(self.data.shape[0], self.data.shape[1]))
mapfile.write("% Optical center (x,y): {0} {1}\n".format(self.center[0], self.center[1]))
mapfile.write("% Step size (x,y): {0} {1}\n".format(self.step_size[0], self.step_size[1]))
mapfile.write("% Scaling: {0}\n".format(float(self.scaling)))
@@ -37,7 +37,11 @@ class surfacemap:
@property
def y(self):
return self.step_size[1] * (numpy.array(range(1, self.data.shape[1]+1))- self.center[1])
-
+
+ @property
+ def size(self):
+ return self.data.shape
+
@property
def offset(self):
return numpy.array(self.step_size)*(self.center - numpy.array(self.size)/2)
@@ -45,13 +49,22 @@ class surfacemap:
@property
def ROMWeights(self):
return self._rom_weights
-
+
+ def z_xy(self, wavelength=1064e-9):
+
+ if "phase" in self.type:
+ k = math.pi * 2 / wavelength
+ return numpy.exp(2j * k * self.scaling * self.data)
+ 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)
self._rom_weights = makeWeights(self, EI)
- return self.ROMWeights
+ return self.ROMWeights, EI
def plot(self, show=True, clabel=None):
diff --git a/pykat/utilities/optics/gaussian_beams.py b/pykat/utilities/optics/gaussian_beams.py
index 02ecfbfcfd28a7fbab6133b5f0230e39ba9d02eb..cd1f554349b352a682df73336048b734f9693f81 100644
--- a/pykat/utilities/optics/gaussian_beams.py
+++ b/pykat/utilities/optics/gaussian_beams.py
@@ -82,7 +82,7 @@ class gauss_param(object):
def w(self):
return abs(self.__q)*math.sqrt(self.__lambda / (self.__nr * math.pi * self.__q.imag))
- def w(self, z=None, wavelength=None, nr=None, w0=None):
+ def beamsize(self, z=None, wavelength=None, nr=None, w0=None):
if z == None:
z = self.z
@@ -144,7 +144,7 @@ class gauss_param(object):
else:
return float("inf")
- def Rc(self, z=None, wavelength=None, nr=None, w0=None):
+ def curvature(self, z=None, wavelength=None, nr=None, w0=None):
if z == None:
z = self.z
else:
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
index 840e7aee915dd4ff155b76f6f72939e76aa9a44c..a4677b2781a74d40d0322a7596cb62176f1a90fe 100644
--- a/pykat/utilities/romhom.py
+++ b/pykat/utilities/romhom.py
@@ -1,14 +1,16 @@
+import math
import maps
import os.path
-import numpy as np
import pykat
import collections
+
from itertools import combinations_with_replacement as combinations
from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
from scipy.linalg import inv
from math import factorial
from hermite import *
-import math
+
+import numpy as np
EmpiricalInterpolant = collections.namedtuple('EmpiricalInterpolant', 'B nodes node_indices limits x')
ReducedBasis = collections.namedtuple('ReducedBasis', 'RB limits x')
@@ -16,11 +18,12 @@ 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, limits):
+ def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, nodes, 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.limits = limits
def writeToFile(self, filename):
@@ -101,7 +104,7 @@ def u(re_q1, w0_1, n1, x):
wz1 = w(w0_1, im_q1, re_q1)
- return A_n1 * hermite(n1, np.sqrt(2.)*x / wz1) * np.exp(-1j*(2*math.pi/(1064e-9))* x**2 /(2.*q_z1))
+ return A_n1 * hermite(n1, np.sqrt(2.)*x / wz1) * np.exp(np.array(-1j*(2*math.pi/(1064e-9))* x**2 /(2.*q_z1)))
def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x):
@@ -240,14 +243,14 @@ def makeEmpiricalInterpolant(RB):
invV = inv(V[0:len(node_indices), 0:len(node_indices)])
B = B_matrix(invV, e_x)
-
- return EmpiricalInterpolant(B=np.array(B), nodes=np.array(x_nodes), node_indices=np.array(node_indices), limits=RB.limits, x=RB.x)
+
+ 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)
def makeWeights(smap, EI):
# get full A_xy
- A_xy = smap.data
+ A_xy = smap.z_xy()
hx = len(smap.x)/2
hy = len(smap.y)/2
@@ -289,83 +292,6 @@ def makeWeights(smap, EI):
B_ij_Q4 = np.outer(B[i], B[j])
w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
-
- 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, limits=EI.limits)
-
-def ROMKnm(W, max_order, q1, q2, q1y=None, q2y=None):
-
- if q1y == None:
- q1y = q1
-
- if q2y == None:
- q2y = q2
-
- # get symmetric and anti-symmetric w_ij's
- w_ij_Q1Q3 = W.w_ij_Q1 + W.w_ij_Q3
- w_ij_Q2Q4 = W.w_ij_Q2 + W.w_ij_Q4
- w_ij_Q1Q2 = W.w_ij_Q1 + W.w_ij_Q2
- w_ij_Q1Q4 = W.w_ij_Q1 + W.w_ij_Q4
- w_ij_Q2Q3 = W.w_ij_Q2 + W.w_ij_Q3
- w_ij_Q3Q4 = W.w_ij_Q3 + W.w_ij_Q4
-
- w_ij_Q1Q2Q3Q4 = W.w_ij_Q1 + W.w_ij_Q3 + W.w_ij_Q2 + W.w_ij_Q4
-
- num_fields = int((max_order + 1) * (max_order + 2) / 2);
-
- K_ROQ = np.array((num_fields, num_fields))
-
- for i in range(len(nm_pairs)):
- for j in range(len(nm_pairs)):
-
- # full quadrature
- n = nm_pairs[i][0]
- m = nm_pairs[i][1]
- npr = nm_pairs[j][0]
- mpr = nm_pairs[j][1]
- u_xy = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, full_x), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, full_x))
- k = dx*dy*np.einsum('ij,ij', u_xy, A_xy)
-
- # ROQ
- if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:
- u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
- elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
- u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
- elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 and nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
- u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
- elif nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1 and nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0:
- u_xy_nodes = np.outer(u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv), u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv))
- k_ROQ = np.einsum('ij,ij', u_xy_nodes, w_ij_Q1Q2Q3Q4)
-
- elif nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 1 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 0:
- u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n,m, nodes_nv)
- u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)
-
- if nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 0 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 1:
- u_xy_nodes_Q1Q4 = np.outer(u_x_nodes, u_y_nodes)
- u_xy_nodes_Q2Q3 = - u_xy_nodes_Q1Q4
- k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q1Q4, w_ij_Q1Q4) + np.einsum('ij,ij', u_xy_nodes_Q2Q3, w_ij_Q2Q3)
- else:
- u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
- u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
- k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3)
-
- elif nm_pairs[j][0] % 2 == 0 and nm_pairs[j][1] % 2 == 1 or nm_pairs[j][0] % 2 == 1 and nm_pairs[j][1] % 2 == 0:
- u_x_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, n, m, nodes_nv)
- u_y_nodes = u_star_u(re_q1, re_q2, w0_1, w0_2, npr, mpr, nodes_nv)
-
- if nm_pairs[i][0] % 2 == 0 and nm_pairs[i][1] % 2 == 0 or nm_pairs[i][0] % 2 == 1 and nm_pairs[i][1] % 2 == 1:
- u_xy_nodes_Q3Q4 = np.outer(u_x_nodes, u_y_nodes)
- u_xy_nodes_Q1Q2 = - u_xy_nodes_Q3Q4
- k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q3Q4, w_ij_Q3Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q2, w_ij_Q1Q2)
- else:
- u_xy_nodes_Q2Q4 = np.outer(u_x_nodes, u_y_nodes)
- u_xy_nodes_Q1Q3 = - u_xy_nodes_Q2Q4
- k_ROQ = np.einsum('ij,ij', u_xy_nodes_Q2Q4, w_ij_Q2Q4) + np.einsum('ij,ij', u_xy_nodes_Q1Q3, w_ij_Q1Q3)
-
+ 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