Skip to content

GitLab

Commit 54f24c7b authored by Daniel Brown's avatar Daniel Brown
Browse files

Removing utlility naming and romhom updates

parent a70516b3
Hide whitespace changes
Inline Side-by-side
bin/test_square_knm.py
View file @ 54f24c7b
......@@ -4,17 +4,17 @@ from pykat.optics.knm import square_aperture_HG_knm, riemann_HG_knm
import math
R = 0.15
q = pykat.beam_param(w0=6e-3, z=2000)
N = 1201
q = pykat.beam_param(w0=0.05, z=0)
N = 1200
mode_i = (0,0)
mode_o = (0,0)
mode_i = (1,1)
mode_o = (1,1)
k = square_aperture_HG_knm(mode_i, mode_o, q, R)
x = y = np.linspace(-R, R, N)
Axy = np.ones((N,N))
k_ = riemann_HG_knm(x, y, mode_i, mode_o, q, q, Axy=Axy, newtonCotesOrder=2)
k_ = riemann_HG_knm(x, y, mode_i, mode_o, q, q, Axy=Axy, newtonCotesOrder=1)
print "%15.15f + %15.15fi" % (k.real, k.imag)
print "%15.15f + %15.15fi" % (k_.real, k_.imag), abs(k-k_)/1e-6, "ppm"
\ No newline at end of file
pykat/__init__.py
View file @ 54f24c7b
......@@ -12,7 +12,7 @@ import components
import detectors
import commands
from pykat.utilities.optics.gaussian_beams import beam_param
from pykat.optics.gaussian_beams import beam_param
......
pykat/commands.py
View file @ 54f24c7b
......@@ -12,7 +12,7 @@ from structs import *
from pykat.param import Param, putter
import pykat.exceptions as pkex
from collections import namedtuple
from pykat.utilities.optics.gaussian_beams import beam_param
from pykat.optics.gaussian_beams import beam_param
class Command(object):
__metaclass__ = abc.ABCMeta
......
pykat/maths/__init__.py
View file @ 54f24c7b
import numpy as np
from scipy.integrate import newton_cotes
def newton_weights(x, newtonCotesOrder):
"""
Constructs the weights for a composite Newton-Cotes rule for integration.
These weights should be multipled by the step-size, equally spaced steps
are assumed. The argument x should be a numpy array of the x samples points.
If the order Newton-Cotes order specified does not produce a rule that fits
into len(x), the order is decremented and that is used to fill gaps that do no fit.
"""
# ensure we have a 1xN array
x = np.array(x).squeeze()
if newtonCotesOrder == 0:
return np.ones(x.shape)
W = newton_cotes(newtonCotesOrder, 1)[0]
wx = np.zeros(x.shape, dtype=np.float64)
N = len(W)
i = 0
while i < len(x) - 1:
try:
wx[i:(i+len(W))] += W
i += len(W)-1
except ValueError as ex:
if len(wx[i:(i+len(W))]) < len(W):
newtonCotesOrder -= 1
W = newton_cotes(newtonCotesOrder, 1)[0]
return wx
\ No newline at end of file
pykat/maths/hermite.py
View file @ 54f24c7b
import numpy as np
def hermite(n, x):
if n == 0:
return 1.0
elif n == 1:
return (2.0 * x)
elif n == 2:
return (4.0 * x * x - 2.0)
elif n == 3:
return (8.0 * x * x * x - 12.0 * x)
elif n == 4:
x_sq = x * x
return (16.0 * x_sq*x_sq - 48.0 * x_sq + 12.0)
elif n == 5:
x_sq = x * x
x_cb = x * x_sq
return (32.0 * x_cb * x_sq - 160.0 * x_cb + 120.0 * x)
elif n == 6:
x_sq = x * x
x_four = x_sq * x_sq
return (64.0 * x_four*x_sq - 480.0 * x_four + 720.0 * x_sq - 120.0)
elif n == 7:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
return (128.0 * x_cb*x_four - 1344.0 * x_cb*x_sq + 3360.0 * x_cb - 1680.0 * x)
elif n == 8:
x_sq = x*x
x_four = x_sq*x_sq
x_six = x_four*x_sq
return (256.0 * x_four*x_four - 3584.0 * x_six + 13440.0 * x_four - 13440.0 * x_sq + 1680.0)
elif n == 9:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
return (512.0 * x_cb*x_cb*x_cb - 9216.0 * x_four*x_cb + 48384.0 * x_cb*x_sq - 80640.0 * x_cb + 30240.0 * x)
elif n == 10:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
x_five = x_four * x
return (1024.0 * x_five*x_five - 23040.0 * x_four*x_four + 161280.0 * x_cb*x_cb - 403200.0 * x_four + 302400.0 * x_sq - 30240.0)
else :
return (2 * x * hermite(n - 1, x) - 2 * (n - 1) * hermite(n - 2, x))
if n == 0:
return 1.0
elif n == 1:
return (2.0 * x)
elif n == 2:
return (4.0 * x * x - 2.0)
elif n == 3:
return (8.0 * x * x * x - 12.0 * x)
elif n == 4:
x_sq = x * x
return (16.0 * x_sq*x_sq - 48.0 * x_sq + 12.0)
elif n == 5:
x_sq = x * x
x_cb = x * x_sq
return (32.0 * x_cb * x_sq - 160.0 * x_cb + 120.0 * x)
elif n == 6:
x_sq = x * x
x_four = x_sq * x_sq
return (64.0 * x_four*x_sq - 480.0 * x_four + 720.0 * x_sq - 120.0)
elif n == 7:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
return (128.0 * x_cb*x_four - 1344.0 * x_cb*x_sq + 3360.0 * x_cb - 1680.0 * x)
elif n == 8:
x_sq = x*x
x_four = x_sq*x_sq
x_six = x_four*x_sq
return (256.0 * x_four*x_four - 3584.0 * x_six + 13440.0 * x_four - 13440.0 * x_sq + 1680.0)
elif n == 9:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
return (512.0 * x_cb*x_cb*x_cb - 9216.0 * x_four*x_cb + 48384.0 * x_cb*x_sq - 80640.0 * x_cb + 30240.0 * x)
elif n == 10:
x_sq = x*x
x_cb = x_sq*x
x_four = x_cb*x
x_five = x_four * x
return (1024.0 * x_five*x_five - 23040.0 * x_four*x_four + 161280.0 * x_cb*x_cb - 403200.0 * x_four + 302400.0 * x_sq - 30240.0)
elif n == 11:
return -665280 * x + 2217600 * x**3 - 1774080 * x**5 + 506880 * x**7 - 56320 * x**9 + 2048 * x**11
elif n == 12:
return 4096 * x**12-135168 * x**10+1520640 * x**8-7096320 * x**6+13305600 * x**4-7983360 * x**2+665280
elif n == 13:
return 8192 * x**13-319488 * x**11+4392960 * x**9-26357760 * x**7+69189120 * x**5-69189120 * x**3+17297280 * x
elif n == 14:
return 16384 * x**14-745472 * x**12+12300288 * x**10-92252160 * x**8+322882560 * x**6-484323840 * x**4+242161920 * x**2-17297280
elif n == 15:
return 32768 * x**15-1720320 * x**13+33546240 * x**11-307507200 * x**9+1383782400 * x**7-2905943040 * x**5+2421619200 * x**3-518918400 * x
elif n == 16:
return 65536 * x**16-3932160 * x**14+89456640 * x**12-984023040 * x**10+5535129600 * x**8-15498362880 * x**6+19372953600 * x**4-8302694400 * x**2+518918400
elif n == 17:
return 131072 * x**17-8912896 * x**15+233963520 * x**13-3041525760 * x**11+20910489600 * x**9-75277762560 * x**7+131736084480 * x**5-94097203200 * x**3+17643225600 * x
elif n == 18:
return 262144 * x**18-20054016 * x**16+601620480 * x**14-9124577280 * x**12+75277762560 * x**10-338749931520 * x**8+790416506880 * x**6-846874828800 * x**4+317578060800 * x**2-17643225600
elif n == 19:
return 524288 * x**19-44826624 * x**17+1524105216 * x**15-26671841280 * x**13+260050452480 * x**11-1430277488640 * x**9+4290832465920 * x**7-6436248698880 * x**5+4022655436800 * x**3-670442572800 * x
elif n == 20:
return 1048576 * x**20-99614720 * x**18+3810263040 * x**16-76205260800 * x**14+866834841600 * x**12-5721109954560 * x**10+21454162329600 * x**8-42908324659200 * x**6+40226554368000 * x**4-13408851456000 * x**2+670442572800
else :
return (2 * x * hermite(n - 1, x) - 2 * (n - 1) * hermite(n - 2, x))
pykat/node_network.py
View file @ 54f24c7b
......@@ -14,7 +14,7 @@ import pykat.exceptions as pkex
from pykat.components import Component, NodeGaussSetter
from pykat.detectors import BaseDetector as Detector
from pykat import beam_param
from pykat.optics.gaussian_beams import beam_param
from copy import deepcopy
class NodeNetwork(object):
......
pykat/optics/knm.py
View file @ 54f24c7b
......@@ -10,6 +10,7 @@ from math import factorial
from pykat.maths.hermite import hermite
from scipy.misc import comb
from scipy.integrate import newton_cotes
from pykat.maths import newton_weights
import time
import pykat.optics.maps
......@@ -303,8 +304,8 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
if isinstance(weights, ROMWeights):
if cache == None:
u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes)
u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes)
u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EIx.nodes)
u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EIy.nodes)
w_ij_Q1Q3 = weights.w_ij_Q1 + weights.w_ij_Q3
w_ij_Q2Q4 = weights.w_ij_Q2 + weights.w_ij_Q4
......@@ -354,8 +355,8 @@ def ROM_HG_knm(weights, mode_in, mode_out, q1, q2, q1y=None, q2y=None, cache=Non
else:
if cache == None:
u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EI["x"].nodes)
u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EI["y"].nodes)
u_x_nodes = u_star_u(q1.z, q2.z, q1.w0, q2.w0, n, m, weights.EIx.nodes)
u_y_nodes = u_star_u(q1y.z, q2y.z, q1y.w0, q2y.w0, npr, mpr, weights.EIy.nodes)
w_ij = weights.w_ij
else:
......@@ -510,9 +511,7 @@ def square_aperture_HG_knm(mode_in, mode_out, q, R):
kx = hg1.constant_x * hg2.constant_x.conjugate()
ky = hg1.constant_y * hg2.constant_y.conjugate()
print hg1.constant_x, hg2.constant_x, hg1.constant_y, hg2.constant_y
f = q.w / math.sqrt(2)
R = R / (q.w / math.sqrt(2))
......
pykat/optics/maps.py
View file @ 54f24c7b
......@@ -13,11 +13,13 @@ from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
from pykat.optics.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights
from scipy.interpolate import interp2d
from pykat.optics.romhom import makeReducedBasis, makeEmpiricalInterpolant, makeWeights, makeWeightsNew
from scipy.interpolate import interp2d, interp1d
from pykat.maths.zernike import *
import numpy as np
import math
from pykat.maths.zernike import *
import pickle
class surfacemap(object):
def __init__(self, name, maptype, size, center, step_size, scaling, data=None):
......@@ -202,7 +204,56 @@ class surfacemap(object):
self._rom_weights = makeWeights(self, EI, verbose=verbose, useSymmetry=useSymmetry)
return self.ROMWeights, EI
def generateROMWeightsNew(self, EIxFilename, EIyFilename=None, verbose=False, interpolate=False):
if interpolate == True:
# Use EI nodes to interpolate if we
with open(EIxFilename, 'rb') as f:
EIx = pickle.load(f)
if EIyFilename is None:
EIy = EIx
else:
with open(EIyFilename, 'rb') as f:
EIy = pickle.load(f)
x = EIx.x
x.sort()
nx = np.unique(np.hstack((x, -x[::-1])))
y = EIy.x
y.sort()
ny = np.unique(np.hstack((y, -y[::-1])))
self.interpolate(nx, ny)
self._rom_weights = makeWeightsNew(self, EIxFilename, EIyFilename, verbose=verbose)
return self.ROMWeights
def interpolate(self, nx, ny, **kwargs):
"""
Interpolates the map for some new x and y values.
Uses scipy.interpolate.interp2d and any keywords arguments are
passed on to it, thus settings like interpolation type and
fill values can be set.
The range of nx and ny must contain the value zero so that the
center point of the map can be set.
"""
D = interp2d(self.x, self.y, self.data, **kwargs)
data = D(nx-self.offset[0], ny-self.offset[0])
Dx = interp1d(nx, np.arange(1,len(nx)+1))
Dy = interp1d(ny, np.arange(1,len(ny)+1))
self.center = (Dx(0), Dy(0))
self.step_size = (nx[1]-nx[0], ny[1]-ny[0])
self.data = data
def plot(self, show=True, clabel=None, xlim=None, ylim=None):
import pylab
......
pykat/optics/romhom.py
View file @ 54f24c7b
......@@ -2,30 +2,40 @@ import math
import os.path
import pykat
import collections
import numpy as np
import multiprocessing
import h5py
import time
import datetime
import pickle
import itertools
from copy import copy
from pykat.external.progressbar import ProgressBar, ETA, Percentage, Bar
from itertools import combinations_with_replacement as combinations
from pykat.optics.gaussian_beams import beam_param, HG_beam
from scipy.linalg import inv
from math import factorial
from pykat.maths.hermite import *
from pykat.maths import newton_weights
from scipy.integrate import newton_cotes
from multiprocessing import Process, Queue, Array, Value, Event
import numpy as np
EmpiricalInterpolant = collections.namedtuple('EmpiricalInterpolant', 'B nodes node_indices limits x')
ReducedBasis = collections.namedtuple('ReducedBasis', 'RB limits x')
ROMLimits = collections.namedtuple('ROMLimits', 'zmin zmax w0min w0max max_order')
ROMLimits = collections.namedtuple('ROMLimits', 'zmin zmax w0min w0max R mapSamples newtonCotesOrder max_order')
class ROMWeights:
def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, EI, limits):
def __init__(self, w_ij_Q1, w_ij_Q2, w_ij_Q3, w_ij_Q4, EIx, EIy):
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.EI = EI
self.limits = limits
self.EIx = EIx
self.EIy = EIy
def writeToFile(self, filename):
"""
......@@ -35,20 +45,20 @@ class ROMWeights:
"""
f = open(filename + ".rom", 'w+')
f.write("zmin=%16.16e\n" % self.limits.zmin)
f.write("zmax=%16.16e\n" % self.limits.zmax)
f.write("w0min=%16.16e\n" % self.limits.w0min)
f.write("w0max=%16.16e\n" % self.limits.w0max)
f.write("maxorder=%i\n" % self.limits.max_order)
f.write("zmin=%16.16e\n" % self.EIx.limits.zmin)
f.write("zmax=%16.16e\n" % self.EIx.limits.zmax)
f.write("w0min=%16.16e\n" % self.EIx.limits.w0min)
f.write("w0max=%16.16e\n" % self.EIx.limits.w0max)
f.write("maxorder=%i\n" % self.EIx.limits.max_order)
f.write("xnodes=%i\n" % len(self.EI["x"].nodes))
f.write("xnodes=%i\n" % len(self.EIx.nodes))
for v in self.EI["x"].nodes.flatten():
for v in self.EIx.nodes.flatten():
f.write("%s\n" % repr(float(v)))
f.write("ynodes=%i\n" % len(self.EI["y"].nodes))
f.write("ynodes=%i\n" % len(self.EIy.nodes))
for v in self.EI["y"].nodes.flatten():
for v in self.EIy.nodes.flatten():
f.write("%s\n" % repr(float(v)))
f.write(repr(self.w_ij_Q1.shape) + "\n")
......@@ -164,7 +174,13 @@ def u_star_u(re_q1, re_q2, w0_1, w0_2, n1, n2, x, x2=None):
return u(re_q1, w0_1, n1, x) * u(re_q2, w0_2, n2, x2).conjugate()
def u_star_u_mm(z, w0, n1, n2, x):
return u(z, w0, n1, x) * u(z, w0, n2, x).conjugate()
################################################################################################
################################################################################################
################################################################################################
def makeReducedBasis(x, isModeMatched=True, tolerance = 1e-12, sigma = 1, greedyfile=None):
if greedyfile != None:
......@@ -316,11 +332,11 @@ def makeEmpiricalInterpolant(RB, sort=False):
return EmpiricalInterpolant(B=B, nodes=nodes, node_indices=node_indices, limits=RB.limits, x=RB.x)
def makeWeights(smap, EI, verbose=True, useSymmetry=True):
def makeWeights(smap, EI, verbose=True, useSymmetry=True, newtonCotesOrder=1):
if useSymmetry:
# get full A_xy
A_xy = smap.z_xy().transpose()
A_xy = smap.z_xy()#.transpose()
xm = smap.x[smap.x < 0]
xp = smap.x[smap.x > 0]
......@@ -358,10 +374,15 @@ def makeWeights(smap, EI, verbose=True, useSymmetry=True):
Bx = EI["x"].B
By = EI["y"].B[:,::-1]
w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex)
from pykat.optics.knm import newton_weights
wx = newton_weights(xm, newtonCotesOrder)
wy = newton_weights(ym, newtonCotesOrder)
for i in range(len(Bx)):
for j in range(len(By)):
B_ij_Q1 = np.outer(Bx[i], By[j])
B_ij_Q1 = np.outer(wx*Bx[i], wy*By[j])
w_ij_Q1[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q1, A_xy_Q1)
if verbose:
......@@ -374,7 +395,7 @@ def makeWeights(smap, EI, verbose=True, useSymmetry=True):
for i in range(len(Bx)):
for j in range(len(By)):
B_ij_Q2 = np.outer(Bx[i], By[j])
B_ij_Q2 = np.outer(wx*Bx[i], wy*By[j])
w_ij_Q2[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q2, A_xy_Q2)
if verbose:
......@@ -387,7 +408,7 @@ def makeWeights(smap, EI, verbose=True, useSymmetry=True):
for i in range(len(Bx)):
for j in range(len(By)):
B_ij_Q3 = np.outer(Bx[i], By[j])
B_ij_Q3 = np.outer(wx*Bx[i], wy*By[j])
w_ij_Q3[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q3, A_xy_Q3)
if verbose:
......@@ -400,7 +421,7 @@ def makeWeights(smap, EI, verbose=True, useSymmetry=True):
for i in range(len(Bx)):
for j in range(len(By)):
B_ij_Q4 = np.outer(Bx[i], By[j])
B_ij_Q4 = np.outer(wx*Bx[i], wy*By[j])
w_ij_Q4[i][j] = dx*dy*np.einsum('ij,ij', B_ij_Q4, A_xy_Q4)
if verbose:
......@@ -447,5 +468,450 @@ def makeWeights(smap, EI, verbose=True, useSymmetry=True):
return ROMWeightsFull(w_ij=w_ij, EI=EI, limits=EI["x"].limits)
###################################################################################################
###################################################################################################
###################################################################################################
# !!! New ROM code below that doesn't need supercomputer
###################################################################################################
###################################################################################################
###################################################################################################
def _compute_TS(queue, oqueue, x, w):
while True:
msg = queue.get()
if msg is None:
break
else:
tmp = u_star_u_mm(msg[0], msg[1], msg[2], msg[3], x)
# includes normalisation with quadrature rule weights
norm = np.sqrt(1/(abs(np.vdot(w*tmp,tmp))))
oqueue.put((msg, tmp*norm))
def _write_TS(queue, filename, tssize):
hfile = h5py.File("%s.h5" % filename, 'a')
i = 0
try:
while True:
msg = queue.get()
if msg is None:
break
else:
# Dump each TS into a group
key = 'TS/%i' % msg[0][4]
hfile[key+"/data"] = msg[1]
hfile[key+"/z"] = msg[0][0]
hfile[key+"/w0"] = msg[0][1]
hfile[key+"/n1"] = msg[0][2]
hfile[key+"/n2"] = msg[0][3]
i += 1
if i % 1000 == 0:
print i/float(tssize)
hfile.flush()
finally:
hfile.close()
def CreateTrainingSetHDF5(filename, maxOrder, z, w0, R, halfMapSamples, newtonCotesOrder=1, NProcesses=1):
"""
newtonCotesOrder: Order of integration to use
0 - Midpoint
1 - Trapezoid
2 - Simpsons
Or higher orders
"""
iq = Queue()
oq = Queue()
Ns = halfMapSamples
h = R / float(Ns-1) # step size
# Close newton-cotes quadrature goes right upto the boundary
# unlike previous midpoint rule.
x = np.linspace(-R, 0, Ns, dtype=np.float64)
w = newton_weights(x, newtonCotesOrder)
nModes = 0
for n in xrange(0, maxOrder+1):
for m in xrange(0, maxOrder+1):
if n+m <= maxOrder and n <= m:
nModes += 1
tssize = len(w0) * len(z) * nModes
hfile = h5py.File("%s.h5" % filename, 'w')
hfile['TSSize'] = tssize
hfile['x'] = x
hfile['zRange'] = (min(z), max(z))
hfile['w0Range'] = (min(w0), max(w0))
hfile['R'] = R
hfile['halfMapSamples'] = halfMapSamples
hfile['maxOrder'] = maxOrder
hfile['newtonCotesOrder'] = newtonCotesOrder
hfile['weights'] = w
hfile.close() # make sure it's closed before
print "Starting processes..."
# Have a bunch of processes doing the computation and one doing the writing
iprocesses = [Process(target=_compute_TS, name="irom%i" % i, args=(iq, oq, x, w)) for i in range(NProcesses)]
oprocess = Process(target=_write_TS, name="orom%i" % i, args=(oq, filename, tssize))
oprocess.start()