diff --git a/pykat/__init__.py b/pykat/__init__.py
index 0f0686889ad479e6ff6a0ccfe59e682b02375c2f..7d2920ba0635ec19cfc8add4e073e4b3a8ec1d79 100644
--- a/pykat/__init__.py
+++ b/pykat/__init__.py
@@ -5,3 +5,4 @@ import components
import detectors
import commands
+from pykat.utilities.optics.gaussian_beams import beam_param
diff --git a/pykat/utilities/knm.py b/pykat/utilities/knm.py
index 7d418c4696bb366a26aa2ae43b0c4392de706969..1607bc26d990bff8813f08f01536d50637f48bbd 100644
--- a/pykat/utilities/knm.py
+++ b/pykat/utilities/knm.py
@@ -2,12 +2,14 @@ from itertools import combinations_with_replacement as combinations
from pykat.utilities.optics.gaussian_beams import beam_param, HG_beam
from pykat.exceptions import BasePyKatException
+import time
import maps
import os.path
import numpy as np
import pykat
import collections
import math
+import cmath
from romhom import u_star_u
from progressbar import ProgressBar, ETA, Percentage, Bar
@@ -32,9 +34,9 @@ def makeCouplingMatrix(max_order):
M.append(row)
- return np.array(M)
+ return np.array(M)
-def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None):
+def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=None, cache=None, delta=(0,0)):
if Axy == None:
Axy == np.ones((len(x), len(y)))
@@ -55,7 +57,7 @@ def riemann_HG_knm(x, y, mode_in, mode_out, q1, q2, q1y=None, q2y=None, Axy=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)
+ U1 = Hg_in.Unm(x+delta[0], y+delta[1])
U2 = Hg_out.Unm(x,y).conjugate()
return dx * dy * np.einsum('ij,ij', Axy, U1*U2)
@@ -66,11 +68,10 @@ 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, np.outer(cache[strx], cache[stry]))
-
-
-def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
+
+def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=(0,0)):
if q1y == None:
q1y = q1
@@ -89,20 +90,24 @@ def __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None):
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])
+ #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()
+ cache[strx] = u_star_u(q1.z, q2.z, q1.w0, q2.w0, mode_in[0], mode_out[0], x, x+delta[0])
+ #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()
+ cache[stry] = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, mode_in[1], mode_out[1], y, y+delta[1])
+ #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:
@@ -141,7 +146,9 @@ def __gen_ROM_HG_knm_cache(weights, couplings, q1, q2, q1y=None, q2y=None):
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
@@ -157,82 +164,154 @@ 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)
-
- 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(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["ym"].nodes)
- 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(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
- 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_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])
+ 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_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_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)
- u_xy_nodes = np.outer(u_x_nodes, u_y_nodes)
+ 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)
- n_mod_2 = n % 2
- m_mod_2 = m % 2
- npr_mod_2 = npr % 2
- mpr_mod_2 = mpr % 2
+
+ return k_ROQ
- 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)
+__fac_cache = []
- 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)
+def fac(n):
+ global __fac_cache
+ return __fac_cache[n]
+ #return math.factorial(int(n))
- 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)
-
+def m_1_pow(n):
+ if n % 2 == 0:
+ return 1
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)
+ return -1
+
+
+def __Ss(u, _u, F, _F, d=0):
+ r = 0
+
+ for s in range(0, min(u,_u)+1):
+ r += m_1_pow(s) * _F ** (u-s) * _F ** (_u-s) / (fac(2*s+d)*fac(u-s)*fac(_u-s))
- 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 r
+
+
+def __S(m, _m, X, _X, F, _F, d=0):
+ if m % 2 == 1:
+ lim1 = (m-1)/2
+ else:
+ lim1 = m/2
+
+ if _m % 2 == 1:
+ lim2 = (_m-1)/2
+ else:
+ lim2 = _m/2
+
+ r = 0
+
+ for u in range(0, lim1+1):
+ for _u in range(0, lim2+1):
+ r += m_1_pow(u) * _X**(m-2*u) * X**(_m-2*_u) / ( fac(m-2*u)*fac(_m-2*_u) ) * __Ss(u, _u, F, _F, d=d)
+
+ return r
+
+
+def __bayerhelms_kn(n, _n, q1, q2, gamma=0.0):
+
+ K0 = (q1.zr - q2.zr)/q2.zr
+ K2 = (q1.z - q2.z)/q2.zr
+ K = (K0 + 1j*K2)/2.0
+
+ Ktilde = abs(K / (1+K))
+
+ if gamma != 0:
+ a = q2.zr * math.sin(gamma) / (cmath.sqrt(1+K.conjugate()) * q2.w0)
+
+ _X = - a * (q2.z/q2.zr - 1j)
+ X = - a * (q2.z/q2.zr + 1j*(1+2*K.conjugate()))
+ Ex = cmath.exp(-_X*X / 2.0)
+ else:
+ _X = 0.0
+ X = 0.0
+ Ex = 1.0
+ _F = K / (2.0 * (1.0+K0))
+ F = K.conjugate() / 2.0
+ Sg = __S(n, _n, X, _X, F, _F)
-def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=True):
+ if n > 0 and _n > 0:
+ Su = __S(n-1, _n-1, X, _X, F, _F, 1)
+ else:
+ Su = 0
+
+ b = m_1_pow(_n) * cmath.sqrt(fac(n) * fac(_n) * (1.0 + K.real)**(n+0.5) * (1.0 + K.conjugate()) ** (-(n+_n+1)))
+
+ return b * Ex * (Sg - Su)
+
+
+def bayerhelms_HG_knm(mode_in, mode_out, q1, q2, q1y=None, q2y=None, gamma=(0,0)):
+ if q1y == None:
+ q1y = q1
+
+ if q2y == None:
+ q2y = q2
+
+ # x modes
+ n = mode_in[0]
+ _n = mode_out[0]
+
+ # y modes
+ m = mode_in[1]
+ _m = mode_out[1]
+
+ return __bayerhelms_kn(n,_n, q1, q2, 2*gamma[0]) * __bayerhelms_kn(m, _m, q1y, q2y, 2*gamma[1])
+
+def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riemann", verbose=False, profile=False, gamma=(0,0), delta=(0,0)):
if q1y == None:
q1y = q1
@@ -247,17 +326,32 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
if int(a) - a != 0:
raise BasePyKatException("Iterator should be product of 4, each element of coupling array should be [n,m,n',m']")
-
- Axy = surface_map.z_xy(q1.wavelength)
- x = surface_map.x
- y = surface_map.y
+ maxtem = 0
+ c = couplings.flatten()
+
+ for i in range(0, c.size/2):
+ maxtem = max(sum(c[i*2:(i*2+2)]), maxtem)
+
+ global __fac_cache
+
+ for n in range(0, maxtem+1):
+ __fac_cache.append(math.factorial(n))
+
+ if surface_map != None:
+ Axy = surface_map.z_xy(q1.wavelength)
+
+ x = surface_map.x
+ y = surface_map.y
K = np.zeros((couplings.size/4,), dtype=np.complex128)
it = np.nditer(couplings, flags=['refs_ok','f_index'])
i = 0
+
+ if profile:
+ t0 = time.time()
if method == "romhom":
if surface_map == None:
@@ -268,32 +362,44 @@ 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)
elif method == "riemann":
- if verbose: print "Computing caches"
- cache = __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None)
+ if surface_map == None:
+ raise BasePyKatException("Using 'riemann' method requires a surface map to be specified")
+
+ cache = __gen_riemann_knm_cache(x, y, couplings, q1, q2, q1y=None, q2y=None, delta=delta)
else:
cache = None
weights = None
+ if profile:
+ cache_time = time.time() - t0
+ Ktime = np.zeros((couplings.size/4,), dtype=np.float64)
+
if verbose:
p = ProgressBar(maxval=couplings.size, widgets=["Knm (%s): " % method, Percentage(), Bar(), ETA()])
while not it.finished:
try:
+ if profile:
+ t0 = time.time()
+
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, cache=cache)
+ K[i] = riemann_HG_knm(x, y, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, Axy=Axy, cache=cache, delta=delta)
elif method == "romhom":
K[i] = ROM_HG_knm(weights, mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, cache=cache)
+ elif method == "bayerhelms":
+ K[i] = bayerhelms_HG_knm(mode_in, mode_out, q1=q1, q2=q2, q1y=q1y, q2y=q2y, gamma=gamma)
else:
raise BasePyKatException("method value '%s' not accepted" % method)
-
+
+ if profile:
+ Ktime[i] = time.time() - t0
+
i +=1
if verbose:
@@ -303,19 +409,41 @@ def knmHG(couplings, q1, q2, surface_map=None, q1y=None, q2y=None, method="riema
except StopIteration:
break
- return K.reshape(couplings.shape[:-1])
-
-
-
-
-
-
-
-
-
-
-
+ if profile:
+ return K.reshape(couplings.shape[:-1]), Ktime.reshape(couplings.shape[:-1]), cache_time
+ else:
+ return K.reshape(couplings.shape[:-1])
+def plot_knm_matrix(couplings, knm):
+ import pylab as plt
+
+ fig = plt.figure()
+ ax = fig.add_subplot(111)
+ cax = ax.imshow(knm, interpolation='nearest')
+ fig.colorbar(cax)
+
+ numrows, numcols = knm.shape
+
+ c = couplings[:, 0, :2]
+ c_ = []
+
+ for d in c:
+ c_.append("%i,%i"%(d[0], d[1]))
+
+ ax.set_xticklabels(c_)
+ ax.set_yticklabels(c_)
+
+ def format_coord(x, y):
+ col = int(x+0.5)
+ row = int(y+0.5)
+
+ if col>=0 and col<numcols and row>=0 and row<numrows:
+ z = knm[row,col]
+ return 'x=%s, y=%s, z=%1.4f' % (c_[col], c_[row], z)
+ else:
+ return 'x=%1.4f, y=%1.4f'%(x, y)
+ ax.format_coord = format_coord
+ plt.show()
diff --git a/pykat/utilities/maps.py b/pykat/utilities/maps.py
index eb663bbc72c1d61d41fc25f9ee95917643c8ae48..8debfd706fa7379b980bbdcd0fc72e6fbece4522 100644
--- a/pykat/utilities/maps.py
+++ b/pykat/utilities/maps.py
@@ -65,25 +65,18 @@ class surfacemap(object):
raise BasePyKatException("Map type needs handling")
- def generateROMWeights(self, isModeMatched=True):
- xp = self.x[self.x>0]
- yp = self.y[self.y>0]
+ def generateROMWeights(self, isModeMatched=True, verbose=False):
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)
+ self._rom_weights = makeWeights(self, EI, verbose=verbose)
return self.ROMWeights, EI
@@ -129,7 +122,7 @@ class tiltmap(surfacemap):
xx, yy = np.meshgrid(self.x, self.y)
- self.data = xx * self.tilt[0] + yy * self.tilt[1]
+ self.data = xx * self.tilt[1] + yy * self.tilt[0]
diff --git a/pykat/utilities/romhom.py b/pykat/utilities/romhom.py
index e134a7b502e1ba6cdea605b8b01ddbafa75d0861..ae77c020c611dc66de1b6865c49823e1fa6d3447 100644
--- a/pykat/utilities/romhom.py
+++ b/pykat/utilities/romhom.py
@@ -1,9 +1,8 @@
import math
-import maps
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
@@ -109,8 +108,11 @@ def u(re_q1, w0_1, n1, x):
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):
- return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x).conjugate()
+def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x, x2=None):
+ if x2 == None:
+ x2 = x
+
+ return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x2).conjugate()
def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
@@ -142,7 +144,7 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
params = np.loadtxt(greedyfile, skiprows=5)
TS_size = len(params) # training space of TS_size number of waveforms
-
+
#### allocate memory for training space ####
TS = np.zeros(TS_size*len(x), dtype = complex).reshape(TS_size, len(x)) # store training space in TS_size X len(x) array
@@ -150,14 +152,13 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
IDx = 0 #TS index
for i in range(len(params)):
-
- q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
-
if isModeMatched:
+ q1 = beam_param(w0=float(params[i][0]), z=float(params[i][1]))
q2 = q1
n = int(params[i][2])
m = int(params[i][3])
else:
+ q1 = beam_param(w0=float(params[i][0]), z=float(params[i][2]))
q2 = beam_param(w0=float(params[i][1]), z=float(params[i][3]))
n = int(params[i][4])
m = int(params[i][5])
@@ -180,34 +181,32 @@ def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1):
#### Begin greedy: see Field et al. arXiv:1308.3565v2 ####
- tolerance = 1e-12 # set maximum RB projection error
-
- sigma = 1 # projection error at 0th iteration
-
RB_matrix = [TS[0]] # seed greedy algorithm (arbitrary)
proj_coefficients = np.zeros(TS_size*TS_size, dtype = complex).reshape(TS_size, TS_size)
projections = np.zeros(TS_size*len(x), dtype = complex).reshape(TS_size, len(x))
iter = 0
-
+
while sigma > tolerance:
#for k in range(TS_size-1):
# go through training set and find worst projection error
projections = project_onto_basis(dx, RB_matrix, TS, projections, proj_coefficients, iter)
residual = TS - projections
+
projection_errors = [np.vdot(dx* residual[i], residual[i]) for i in range(len(residual))]
sigma = abs(max(projection_errors))
index = np.argmax(projection_errors)
-
+
#Gram-Schmidt to get the next basis and normalize
next_basis = TS[index] - projections[index]
next_basis /= np.sqrt(abs(np.vdot(dx* next_basis, next_basis)))
RB_matrix.append(next_basis)
-
- iter += 1
-
+
+ iter += 1
+
+
return ReducedBasis(RB=np.array(RB_matrix), limits=romlimits, x=x)
def makeEmpiricalInterpolant(RB):
@@ -282,42 +281,39 @@ def makeWeights(smap, EI, verbose=True):
dy = full_y[1] - full_y[0]
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)
+ count = 4*len(EI["xm"].B) * len(EI["ym"].B)
p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()])
n = 0
# make integration weights
Bx = EI["xm"].B
- By = EI["yp"].B
+ By = EI["ym"].B[:,::-1]
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)
-
+
if verbose:
p.update(n)
n+=1
- Bx = EI["xp"].B
- By = EI["yp"].B
+ Bx = EI["xm"].B[:,::-1]
+ By = EI["ym"].B[:,::-1]
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)
-
+
if verbose:
p.update(n)
n+=1
- Bx = EI["xp"].B
+ Bx = EI["xm"].B[:,::-1]
By = EI["ym"].B
w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex)