diff --git a/py/helper/__init__.py b/py/helper/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/py/helper/basis/__init__.py b/py/helper/basis/__init__.py new file mode 100644 index 0000000..d2df164 --- /dev/null +++ b/py/helper/basis/__init__.py @@ -0,0 +1,27 @@ +#!/usr/bin/python3 + +import os + +capitalizationList = [ + "BSpline", + "ClenshawCurtis", + "NotAKnot", + "SGpp", +] + +capitalizationListLower = [word.lower() for word in capitalizationList] +dirName = os.path.dirname(__file__) + +for f in os.listdir(dirName): + if (os.path.isfile(os.path.join(dirName, f)) and + (f.endswith(".py")) and (f != "__init__.py")): + moduleName = f[:-3] + words = moduleName.split("_") + words = [word.title() for word in words] + words = [(word if word.lower() not in capitalizationListLower else + capitalizationList[capitalizationListLower.index(word.lower())]) + for word in words] + className = "".join(words) + result = __import__(moduleName, globals=globals(), locals=locals(), + fromlist=[className], level=1) + vars()[className] = getattr(result, className) diff --git a/py/helper/basis/cardinal_bspline.py b/py/helper/basis/cardinal_bspline.py new file mode 100644 index 0000000..0beda74 --- /dev/null +++ b/py/helper/basis/cardinal_bspline.py @@ -0,0 +1,7 @@ +#!/usr/bin/python3 + +from .non_uniform_bspline import NonUniformBSpline + +class CardinalBSpline(NonUniformBSpline): + def __init__(self, p, nu=0): + super().__init__(p, list(range(0, p+2)), nu=nu) diff --git a/py/helper/basis/centralized_cardinal_bspline.py b/py/helper/basis/centralized_cardinal_bspline.py new file mode 100644 index 0000000..d6f10fd --- /dev/null +++ b/py/helper/basis/centralized_cardinal_bspline.py @@ -0,0 +1,9 @@ +#!/usr/bin/python3 + +import numpy as np + +from .non_uniform_bspline import NonUniformBSpline + +class CentralizedCardinalBSpline(NonUniformBSpline): + def __init__(self, p, nu=0): + super().__init__(p, np.linspace(-(p+1)/2, (p+1)/2, p+2), nu=nu) diff --git a/py/helper/figure.py b/py/helper/figure.py new file mode 100644 index 0000000..eb78699 --- /dev/null +++ b/py/helper/figure.py @@ -0,0 +1,221 @@ +import hashlib +import lzma +import os +import pathlib +import pickle +import platform +import subprocess +import sys + +import cycler +import matplotlib as mpl +import matplotlib.pyplot as plt +plt.switch_backend("pgf") +from mpl_toolkits.mplot3d import Axes3D + + + +class Figure(mpl.figure.Figure): + COLORS = { + "anthrazit" : ( 62/255, 68/255, 76/255), + "mittelblau" : ( 0/255, 81/255, 158/255), + "hellblau" : ( 0/255, 190/255, 255/255), + } + + _TEX_PREAMBLE_COMMON = r""" +\usepackage[ngerman,american]{babel} +\usepackage{mathtools} +\usepackage[utf8]{luainputenc} +\usepackage[T1]{fontenc} + +\usepackage{xcolor} + +\definecolor{anthrazit}{RGB}{62,68,76} +\definecolor{mittelblau}{RGB}{0,81,158} +\definecolor{hellblau}{RGB}{0,190,255} + +\definecolor{C0}{rgb}{0.000,0.447,0.741} +\definecolor{C1}{rgb}{0.850,0.325,0.098} +\definecolor{C2}{rgb}{0.929,0.694,0.125} +\definecolor{C3}{rgb}{0.494,0.184,0.556} +\definecolor{C4}{rgb}{0.466,0.674,0.188} +\definecolor{C5}{rgb}{0.301,0.745,0.933} +\definecolor{C6}{rgb}{0.635,0.078,0.184} +\definecolor{C7}{rgb}{0.887,0.465,0.758} +\definecolor{C8}{rgb}{0.496,0.496,0.496} + +% prevent Matplotlib from loading fontspec as this resets the default text font +% (the default math font will be alright with the no-math option to fontspec, +% but the text font will still be CM...) +\makeatletter +\@namedef{ver@fontspec.sty}{} +\makeatother + +%\PassOptionsToPackage{no-math}{fontspec} +""" + + _TEX_PREAMBLE_SPECIAL = { + "beamer" : r""" +\usepackage[scaled]{helvet} +\renewcommand*{\familydefault}{\sfdefault} +\usepackage{sfmath} +""", + "paper" : "", + "thesis" : r""" +\usepackage[bitstream-charter]{mathdesign} +\renewcommand*{\vec}[1]{{\boldsymbol{#1}}} +\renewcommand*{\*}[1]{\vec{#1}} +""", + } + + _LINE_COLORS = [ + (0.000, 0.447, 0.741), + (0.850, 0.325, 0.098), + (0.929, 0.694, 0.125), + (0.494, 0.184, 0.556), + (0.466, 0.674, 0.188), + (0.301, 0.745, 0.933), + (0.635, 0.078, 0.184), + (0.887, 0.465, 0.758), + (0.496, 0.496, 0.496), + ] + + graphicsCounter = 0 + + def __init__(self, *args, fontSize=11, preamble="", mode=None, **kwargs): + super(Figure, self).__init__(*args, **kwargs) + + if mode is None: + if "talk" in Figure._getBuildDir(): + self.mode = "beamer" + elif "thesis" in Figure._getBuildDir(): + self.mode = "thesis" + else: + self.mode = "paper" + else: + self.mode = mode + + fontFamily = ("sans-serif" if self.mode == "beamer" else "serif") + preamble = (Figure._TEX_PREAMBLE_COMMON + + Figure._TEX_PREAMBLE_SPECIAL[self.mode] + + preamble) + + mpl.rcParams.update({ + "axes.prop_cycle" : cycler.cycler(color=Figure._LINE_COLORS), + "font.family" : fontFamily, + "font.size" : fontSize, + "lines.linewidth" : 1, + "pgf.texsystem" : "lualatex", + "pgf.rcfonts" : False, + "pgf.preamble" : preamble.splitlines(), + "text.usetex" : True, + }) + + self._saveDisabled = (platform.node() == "neon") + + @staticmethod + def create(*args, scale=1, **kwargs): + if "figsize" in kwargs: kwargs["figsize"] = [scale * x for x in kwargs["figsize"]] + return plt.figure(*args, FigureClass=Figure, **kwargs) + + @staticmethod + def _getBuildDir(): + return os.path.realpath(os.environ["BUILD_DIR"]) + + @staticmethod + def _getGraphicsBasename(): + return os.path.splitext(os.path.split(os.path.realpath(sys.argv[0]))[1])[0] + + @staticmethod + def _computeHash(path): + try: + with open(path, "rb") as f: return hashlib.md5(f.read()).digest() + except: + return None + + @staticmethod + def load(datPath=None): + Figure.graphicsCounter += 1 + graphicsNumber = Figure.graphicsCounter + + if path is None: + buildDir = Figure._getBuildDir() + graphicsBasename = Figure._getGraphicsBasename() + basename = os.path.join(buildDir, "{}_{}".format(graphicsBasename, graphicsNumber)) + datPath = "{}.pickle.dat".format(basename) + + print("Loading {}...".format(datPath)) + with lzma.open(datPath, "rb") as f: return pickle.load(f) + + def disableSave(): + self._saveDisabled = True + + def enableSave(): + self._saveDisabled = False + + def save(self, graphicsNumber=None, appendGraphicsNumber=True, + hideSpines=True, tightLayout=True, crop=True, close=True, + transparent=True): + plt.figure(self.number) + + if graphicsNumber is None: + Figure.graphicsCounter += 1 + graphicsNumber = Figure.graphicsCounter + else: + Figure.graphicsCounter = graphicsNumber + + if self._saveDisabled: + if close: plt.close(self) + return + + if tightLayout is not False: + if tightLayout is True: tightLayout = {} + self.tight_layout(**tightLayout) + + if hideSpines: + for ax in self.axes: + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + + if graphicsNumber is None: graphicsNumber = self.number + + buildDir = Figure._getBuildDir() + graphicsBasename = Figure._getGraphicsBasename() + + basename = graphicsBasename + if appendGraphicsNumber: basename += "_{}".format(graphicsNumber) + basename = os.path.join(buildDir, basename) + + if transparent: + savefigFcn = ( + lambda path: plt.savefig(path, facecolor="none", transparent=True)) + else: + savefigFcn = ( + lambda path: plt.savefig(path, facecolor=self.get_facecolor())) + + pgfPath = "{}.pgf".format(basename) + print("Saving {}...".format(os.path.split(pgfPath)[1])) + savefigFcn(pgfPath) + + pgfXzPath = "{}.pgf.xz".format(basename) + oldHash = Figure._computeHash(pgfXzPath) + with open(pgfPath, "rb") as f: pgf = f.read() + with lzma.open(pgfXzPath, "wb") as f: pickle.dump(pgf, f) + os.remove(pgfPath) + newHash = Figure._computeHash(pgfXzPath) + + pdfPath = "{}.pdf".format(basename) + + if (oldHash == newHash) and os.path.isfile(pdfPath): + print("No changes since last run.") + pathlib.Path(pdfPath).touch() + else: + print("Compiling to {}...".format(os.path.split(pdfPath)[1])) + savefigFcn(pdfPath) + if crop: subprocess.run(["pdfcrop", pdfPath, pdfPath], check=True) + + datPath = "{}.pickle.xz".format(basename) + print("Saving {}...".format(os.path.split(datPath)[1])) + with lzma.open(datPath, "wb") as f: pickle.dump(self, f) + + if close: plt.close(self) diff --git a/py/helper/grid.py b/py/helper/grid.py new file mode 100644 index 0000000..678ab82 --- /dev/null +++ b/py/helper/grid.py @@ -0,0 +1,244 @@ +#!/usr/bin/python3 + +import numpy as np +import scipy.misc + + + +def getNodalIndices(l): + i1D = [getNodalIndices1D(l1D) for l1D in l] + meshGrid = np.meshgrid(*i1D, indexing="ij") + I = np.column_stack([grid.flatten() for grid in meshGrid]) + return I + +def getNodalIndices1D(l): + i = list(range(2**l + 1)) + return i + +def getHierarchicalIndices(l): + i1D = [getHierarchicalIndices1D(l1D) for l1D in l] + meshGrid = np.meshgrid(*i1D, indexing="ij") + I = np.column_stack([grid.flatten() for grid in meshGrid]) + return I + +def getHierarchicalIndices1D(l): + i = (list(range(1, 2**l, 2)) if l > 0 else [0, 1]) + return i + +def getCoordinates(L, I, distribution="uniform"): + L, I = np.array(L), np.array(I) + if L.size == 1: L = L * np.ones_like(I) + + if distribution == "uniform": + X = I / 2**L + elif distribution == "clenshawCurtis": + HInv = 2**L + K1 = (I < 0) + K2 = (I > HInv) + K = np.logical_not(np.logical_or(K1, K2)) + X = np.zeros(L.shape) + X[K] = (1 - np.cos(np.pi * I[K] / HInv[K])) / 2 + X[K1] = I[K1] * ((1 - np.cos(np.pi / HInv[K1])) / 2) + X[K2] = 1 + (I[K2] - HInv[K2]) * ((1 - np.cos(np.pi / HInv[K2])) / 2) + else: + raise NotImplementedError("Unknown grid point distribution.") + + return X + +def generateMeshGrid(nns): + assert len(nns) == 2 + xx0 = np.linspace(0, 1, nns[0]) + xx1 = np.linspace(0, 1, nns[1]) + XX0, XX1 = np.meshgrid(xx0, xx1) + XX = flattenMeshGrid((XX0, XX1)) + return XX0, XX1, XX + +def flattenMeshGrid(XXs): + XX = np.stack([XXt.flatten() for XXt in XXs], axis=1) + return XX + + + +class RegularSparse(object): + def __init__(self, n, d): + self.n = n + self.d = d + + def getSize(self): + if self.n <= 0: + return (1 if self.d == 0 else 0) + elif self.d == 0: + return 1 + elif self.d < 0: + return 0 + else: + return sum([2**q * scipy.misc.comb(self.d - 1 + q, self.d - 1, exact=True) + for q in range(self.n - self.d + 1)]) + + def generate(self): + n = self.n + d = self.d + N = self.getSize() + L = np.zeros((N, d)) + I = np.zeros((N, d)) + l = np.ones((d,)) + l1 = d + i = np.ones((d,)) + t = 0 + m = 0 + + def generateRecursively(): + nonlocal n, d, L, I, l, l1, i, t, m + + if t == d: + L[m,:] = l + I[m,:] = i + m += 1 + elif l1 <= n: + l[t] += 1 + l1 += 1 + i[t] = 2 * i[t] - 1 + generateRecursively() + i[t] += 2 + generateRecursively() + l[t] -= 1 + l1 -= 1 + i[t] = (i[t] - 1) / 2 + + t += 1 + generateRecursively() + t -= 1 + + generateRecursively() + X = I * 2**(-L) + + return X, L, I + + + +class RegularSparseBoundary(object): + def __init__(self, n, d, b): + self.n = n + self.d = d + self.b = b + + def getSize(self): + if self.b == 0: + return sum([2**q * scipy.misc.comb(self.d, q, exact=True) * + RegularSparse(self.n, self.d - q).getSize() + for q in range(self.d + 1)]) + elif self.b >= 1: + return (RegularSparse(self.n, self.d).getSize() + + sum([2**q * scipy.misc.comb(self.d, q, exact=True) * + RegularSparse(self.n - q - self.b + 1, self.d - q).getSize() + for q in range(1, self.d + 1)])) + else: + raise ValueError("Invalid value for b.") + + def generate(self): + n, d, b = self.n, self.d, self.b + + if self.b == 0: + oldL = [[l] for l in range(n+1)] + newL = list(oldL) + + for t in range(2, d+1): + newL = [] + + for l in oldL: + lNorm = sum(l) + lStar = n - lNorm + newL.extend([(*l, lt) for lt in range(0, lStar+1)]) + + oldL = list(newL) + + L = np.array(newL) + + elif self.b >= 1: + oldL = [[l] for l in range(n-d+2)] + newL = list(oldL) + + for t in range(2, d+1): + newL = [] + + for l in oldL: + lNorm = sum(l) + Nl = sum([lt == 0 for lt in l]) + + if (lNorm + Nl <= n - d + t - b) or (Nl == t-1): + newL.append((*l, 0)) + + if Nl == 0: lStar = n - d + t - lNorm + else: lStar = n - d + t - b + 1 - lNorm - Nl + newL.extend([(*l, lt) for lt in range(1, lStar+1)]) + + oldL = list(newL) + + L = np.array(newL) + + else: + raise ValueError("Invalid value for b.") + + return DimensionallyAdaptiveSparse(L).generate() + + + +class DimensionallyAdaptiveSparse(object): + def __init__(self, L): + self.L = L + + def generate(self): + I = [getHierarchicalIndices(self.L[k,:]) for k in range(self.L.shape[0])] + + if len(I) > 0: + L = np.vstack([np.tile(self.L[k,:], [I[k].shape[0], 1]) + for k in range(self.L.shape[0])]) + I = np.vstack(I) + else: + L = [] + + X = getCoordinates(L, I) + return X, L, I + + + +class SGppGrid(object): + def __init__(self, grid, *args): + if type(grid) is str: + import pysgpp + self.label = "{}({})".format(grid, ", ".join([str(arg) for arg in args])) + gridTypes = { + "bSpline" : pysgpp.Grid.createBsplineBoundaryGrid, + "bSplineNoBoundary" : pysgpp.Grid.createBsplineGrid, + "modifiedBSpline" : pysgpp.Grid.createModBsplineGrid, + "modifiedNotAKnotBSpline" : pysgpp.Grid.createModNotAKnotBsplineGrid, + "notAKnotBSpline" : pysgpp.Grid.createNotAKnotBsplineBoundaryGrid, + } + grid = gridTypes[grid](*args) + + self.grid = grid + + def generateRegular(self, n): + self.grid.getStorage().clear() + self.grid.getGenerator().regular(n) + return self.getPoints() + + def getPoints(self): + gridStorage = self.grid.getStorage() + N = gridStorage.getSize() + d = gridStorage.getDimension() + L = np.zeros((N, d), dtype=np.uint64) + I = np.zeros((N, d), dtype=np.uint64) + + for k in range(N): + gp = gridStorage.getPoint(k) + + for t in range(d): + L[k,t] = gp.getLevel(t) + I[k,t] = gp.getIndex(t) + + X = getCoordinates(L, I) + return X, L, I + + def __str__(self): + return self.label