diff --git a/bin/test_square_knm.py b/bin/test_square_knm.py index c1950c9807926489c0e05d185961065033c0baaf..5afe81d29f42c87b4a6becbaa4abf1f6dab30f08 100644 --- a/bin/test_square_knm.py +++ b/bin/test_square_knm.py @@ -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 diff --git a/pykat/__init__.py b/pykat/__init__.py index 3008f8b854cdea48a960f2090521cf3f60f5616e..9b14f735d2213da106e9b1cd22c2b5be53a8cd38 100644 --- a/pykat/__init__.py +++ b/pykat/__init__.py @@ -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 diff --git a/pykat/commands.py b/pykat/commands.py index 8300b1805ec06a2ac5691e1e0d8c0f3bfd86111a..6f0f1c89dfb111a55d78e3efcda25923b305dc4a 100644 --- a/pykat/commands.py +++ b/pykat/commands.py @@ -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 diff --git a/pykat/maths/__init__.py b/pykat/maths/__init__.py index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0d2742adce2ff0cb2ec965ab001b3ffab962b963 100644 --- a/pykat/maths/__init__.py +++ b/pykat/maths/__init__.py @@ -0,0 +1,40 @@ +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 diff --git a/pykat/maths/hermite.py b/pykat/maths/hermite.py index aa4152fa65f3b3351bdaff5ecaeac98870256294..84f729ccd554046db57bb96cf353a680d2e5d9bc 100644 --- a/pykat/maths/hermite.py +++ b/pykat/maths/hermite.py @@ -1,69 +1,65 @@ 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)) diff --git a/pykat/node_network.py b/pykat/node_network.py index d1c2745a683b1d95f069aa78549342508fafb9b6..1a565b8354cf5cdf27209e81a5c4f0aa83459605 100644 --- a/pykat/node_network.py +++ b/pykat/node_network.py @@ -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): diff --git a/pykat/optics/knm.py b/pykat/optics/knm.py index 881efc4fad864c47603083e44b65e534929c8320..4ecd9b6a2a6551739ed5f81f46110b9b20a12d1a 100644 --- a/pykat/optics/knm.py +++ b/pykat/optics/knm.py @@ -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)) diff --git a/pykat/optics/maps.py b/pykat/optics/maps.py index bee94b8f72b42419307b90c876db87ab1f45685f..d96ae12ced82da0c6d989149f4b1b6d6bc4a11cf 100644 --- a/pykat/optics/maps.py +++ b/pykat/optics/maps.py @@ -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 diff --git a/pykat/optics/romhom.py b/pykat/optics/romhom.py index 1938d249f8e088c5de92cbd9e6933b9f259b418e..9f96d2c858d3eeb724ca939a831150a144641197 100644 --- a/pykat/optics/romhom.py +++ b/pykat/optics/romhom.py @@ -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() + try: + for P in iprocesses: + P.start() + + print "Filling queue..." + curr = 0 + for n in xrange(0, maxOrder+1): + for m in xrange(0, maxOrder+1): + if n+m <= maxOrder and n <= m: + for _z in z: + for _w0 in w0: + iq.put((_z, _w0, n, m, curr)) + # Can't use Queue.qsize() on OSX so count normally... + curr += 1 + + for P in iprocesses: + iq.put(None) # put a None for each thread to catch to end + + for P in iprocesses: + P.join() + + # Put none to stop output process + oq.put(None) + + oprocess.join() + + except: + print("Exception occurred") + + for P in iprocesses: + P.terminate() + + oprocess.terminate() + + print("Completed training set creation") + print("Data written to %s.h5" % filename) + + +def _worker_ROM(hdf5Filename, job_queue, result_err, result_idx, event): + # h5py drivers: 'core', 'sec2', 'stdio', 'mpio' + # Need to use something ot her than sec2, the default on OSX, + # as it doesn't play nice with multiple processes reading files + with h5py.File("%s.h5" % hdf5Filename, driver="core", mode='r') as file: + + TS = file["TS"] + + while True: + + msg = job_queue.get() + + if msg is None: + break + else: + TSidx, B, EI_indices = msg + TSidx = np.array(TSidx) + + max_err = [] + + for ll in TSidx: + a = TS['%i/data' % ll].value + + _err = np.max(np.abs(a - emp_interp(B, a, EI_indices))) + + max_err.append(_err) + + result_err.value = np.max(max_err) + result_idx.value = TSidx[np.argmax(max_err)] + + event.set() + +def MakeROMFromHDF5(hdf5Filename, greedyFilename=None, EIFilename=None, tol=1e-10, NProcesses=1, maxRBsize=50): + start = time.time() + + #### Start reading TS file #### + TSdata = h5py.File("%s.h5" % hdf5Filename, 'r') + TS = TSdata['TS'] + + quadratureWeights = TSdata['weights'][...] + x = TSdata['x'][...] + maxRBsize = maxRBsize + TSsize = TSdata['TSSize'][...] + + #### Set up stuff for greedy #### + tol = tol + rb_errors = [] + x_nodes = [] + + # Initial RB to seed with + next_RB_index = 0 + #EI_indices = [next_RB_index] + #x_nodes.append(x[EI_indices]) + + RB0 = TS['%i/data' % next_RB_index][...] + #RB0 /= RB0[EI_indices[0]] + EI_indices = [np.argmax(RB0)] + x_nodes.extend(x[EI_indices]) + RB0 /= RB0[EI_indices[0]] + RB_matrix = [RB0] + + V = np.zeros(((maxRBsize), (maxRBsize)), dtype=complex) + + V[0][0] = RB_matrix[0][EI_indices[0]] + invV = inv(V[0:len(EI_indices), 0:len(EI_indices)]) + B = B_matrix(invV, np.array(RB_matrix)) + + RBs = [] + + if NProcesses > 1: + queue = Queue() + locks = [Event() for l in range(NProcesses)] + + result_err = [Value('d', np.inf) for l in range(NProcesses)] + result_idx = [Value('i', -1) for l in range(NProcesses)] + + Names = range(NProcesses) + procs = [Process(name="process_%i" % l[0], target=_worker_ROM, + args=(hdf5Filename, queue, l[1], l[2], l[3])) for l in itertools.izip(Names, result_err, result_idx, locks)] + + max_res = np.zeros((NProcesses), dtype='d') + max_idx = np.zeros((NProcesses), dtype='i') + + for P in procs: P.start() + + dstr = datetime.datetime.strftime(datetime.datetime.now(), "%d%m%Y_%H%M%S") + + if greedyFilename is None: + greedyFilename = "GreedyPoints_%s" % dstr + + greedyFilename += ".dat" + + limits = ROMLimits(zmin=min(TSdata['zRange'].value), + zmax=max(TSdata['zRange'].value), + w0min=min(TSdata['w0Range'].value), + w0max=max(TSdata['w0Range'].value), + R=TSdata['R'].value, + mapSamples=TSdata['halfMapSamples'].value, + newtonCotesOrder=int(TSdata['newtonCotesOrder'].value), + max_order=int(TSdata['maxOrder'].value)) + + with open(greedyFilename, "w") as f: + f.write("min w0 = %15.15e\n" % limits.zmin) + f.write("max w0 = %15.15e\n" % limits.zmax) + f.write("min z = %15.15e\n" % limits.w0min) + f.write("min z = %15.15e\n" % limits.w0max) + f.write("R = %15.15e\n" % limits.R) + f.write("max order = %i\n" % limits.max_order) + f.write("NC order = %i\n" % limits.newtonCotesOrder) + f.write("half map samples = %i\n" % limits.mapSamples) + + # write initial RB + _TS = TS[str(next_RB_index)] + f.write("%15.15e %15.15e %i %i\n" % (_TS["z"].value, _TS["w0"].value, _TS["n1"].value, _TS["n2"].value)) + + for k in range(1, maxRBsize): + + if NProcesses == 1: + max_res = [] + + TSidx = np.array(range(TSsize)) + + for ll in TSidx: + a = TS['%i/data' % ll].value + max_res.append(np.max(a - emp_interp(B, a, EI_indices))) + + worst_error = max(np.abs(max_res)) + next_RB_index = np.argmax(max_res) + + else: + + TSs = [range(TSsize)[i::NProcesses] for i in range(NProcesses)] + + for l in TSs: queue.put((l, B, EI_indices)) + + end_locks = copy(locks) + + while len(end_locks) > 0: + for e in end_locks: + if e.wait(): + end_locks.remove(e) + + for e in locks: e.clear() + + for i in range(NProcesses): + max_res[i] = result_err[i].value + max_idx[i] = result_idx[i].value + + worst_error = max(np.abs(max_res)) + next_RB_index = max_idx[np.argmax(max_res)] + + if worst_error <= tol: + print "Final basis size = %d, Final error = %e, Tolerance=%e" % (k, worst_error, tol) + break + + print "worst error = %e at %i on iteration %d" % (worst_error, next_RB_index, k) + + epsilon = TS['%i/data' % next_RB_index].value + res = epsilon - emp_interp(B, epsilon, EI_indices) + index = np.argmax(res) + EI_indices.append(index) + x_nodes.append(TSdata['x'][index]) + RB_matrix.append(res/max(res)) + + for l in range(len(EI_indices)): # Part of (5) of Algorithm 2: making V_{ij} + for m in range(len(EI_indices)): # Part of (5) of Algorithm 2: making V_{ij} + V[m][l] = RB_matrix[l][EI_indices[m]] # Part of (5) of Algorithm 2: making V_{ij} + + invV = inv(V[0:len(EI_indices), 0:len(EI_indices)]) + B = B_matrix(invV, np.array(RB_matrix)) + + _TS = TS[str(next_RB_index)] + f.write("%15.15e %15.15e %i %i\n" % (_TS["w0"].value, _TS["z"].value, _TS["n1"].value, _TS["n2"].value)) + + print time.time() - start, "Seconds" + + if NProcesses > 1: + for P in procs: P.terminate() + + TSdata.close() + + greedyFilenameBase = os.path.splitext(greedyFilename)[0] + + print "Writing to %s" % greedyFilename + + EI = EmpiricalInterpolant(B=np.matrix(B).real, + nodes=np.array(x_nodes).squeeze(), + node_indices=np.array(EI_indices).squeeze(), + limits=limits, + x=x.squeeze()) + + if EIFilename is not None: + with open("%s.p" % EIFilename, 'wb') as f: + pickle.dump(EI, f) + + print "Writing to %s.p" % EIFilename + + return EI + + +def makeWeightsNew(smap, EIxFilename, EIyFilename=None, verbose=True, newtonCotesOrder=1): + with open("%s" % EIxFilename, 'rb') as f: + EIx = pickle.load(f) + + if EIyFilename is None: + EIy = EIx + else: + with open("%s" % EIyFilename, 'rb') as f: + EIy = pickle.load(f) + + # get full A_xy + A_xy = smap.z_xy().transpose() + + 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[Q1xy] + A_xy_Q2 = A_xy[Q2xy] + A_xy_Q3 = A_xy[Q3xy] + A_xy_Q4 = A_xy[Q4xy] + + full_x = smap.x + full_y = smap.y + + dx = full_x[1] - full_x[0] + dy = full_y[1] - full_y[0] + + if verbose: + count = 4*len(EIx.B) * len(EIy.B) + p = ProgressBar(maxval=count, widgets=["Computing weights: ", Percentage(), Bar(), ETA()]) + + n = 0 + + # make integration weights + Bx = EIx.B + By = EIy.B[:,::-1] + w_ij_Q1 = np.zeros((len(Bx),len(By)), dtype = complex) + + wx = newton_weights(Bx[0], EIx.limits.newtonCotesOrder) + wy = newton_weights(By[0], EIy.limits.newtonCotesOrder) + W = np.outer(wx, wy) + + A = A_xy_Q1 * W + + 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) + + if verbose: + p.update(n) + n+=1 + + Bx = EIx.B[:,::-1] + By = EIy.B[:,::-1] + w_ij_Q2 = np.zeros((len(Bx),len(By)), dtype = complex) + A = A_xy_Q2 * W + + 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) + + if verbose: + p.update(n) + n+=1 + + Bx = EIx.B[:,::-1] + By = EIy.B + w_ij_Q3 = np.zeros((len(Bx),len(By)), dtype = complex) + A = A_xy_Q3 * W + + 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) + + if verbose: + p.update(n) + n+=1 + + Bx = EIx.B + By = EIy.B + w_ij_Q4 = np.zeros((len(Bx),len(By)), dtype = complex) + A = A_xy_Q4 * W + + 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) + + 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, EIx=EIx, EIy=EIy) \ No newline at end of file