diff --git a/examples/diapason/pod.py b/examples/diapason/pod.py
index 659f67c..89eb30d 100644
--- a/examples/diapason/pod.py
+++ b/examples/diapason/pod.py
@@ -1,150 +1,149 @@
import numpy as np
from diapason_engine import DiapasonEngine, DiapasonEngineDamped
from rrompy.reduction_methods.distributed import RationalInterpolant as Pade
from rrompy.reduction_methods.distributed import RBDistributed as RB
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
-
verb = 100
sol = "single"
sol = "sweep"
algo = "Pade"
#algo = "RB"
polyBasis = "LEGENDRE"
polyBasis = "CHEBYSHEV"
#polyBasis = "MONOMIAL"
dampingEta = 0. * 1e4 / 2. / np.pi
ktar = 1.e4 # [Hz]
k0s = [2.5e2, 1.0e4]
#k0s = np.array([2.5e3, 1.5e4])
#k0s = np.array([5.0e4, 1.0e5])
k0s = [2.0e5, 3.0e5]
k0 = np.mean(np.power(k0s, 2.)) ** .5
theta = 20. * np.pi / 180.
phi = 10. * np.pi / 180.
c = 3.e2
rho = 8e3 * (2. * np.pi) ** 2.
E = 1.93e11
nu = .3
T = 1e6
###
if np.isclose(dampingEta, 0.):
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
solver = DiapasonEngine(kappa = k0, c = c, rho = rho, E = E, nu = nu,
T = T, theta = theta, phi = phi, meshNo = 1,
degree_threshold = 8, verbosity = 0)
else:
rescaling = lambda x: x
rescalingInv = lambda x: x
solver = DiapasonEngineDamped(kappa = k0, c = c, rho = rho, E = E, nu = nu,
T = T, theta = theta, phi = phi,
dampingEta = dampingEta, meshNo = 1,
degree_threshold = 8, verbosity = 0)
params = {'N':39, 'M':39, 'S':40, 'POD':True, 'polybasis':polyBasis,
'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)}#,
# 'robustTol':1e-16}
if algo == "Pade":
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
else:
params.pop("N")
params.pop("M")
params.pop("polybasis")
# params.pop("robustTol")
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if sol == "single":
approx.outParaviewTimeDomainSamples(
filename = "out/outSamples{}".format(dampingEta),
forceNewFile = False, folders = True)
nameBase = "{}_{}".format(ktar, dampingEta)
approx.outParaviewTimeDomainApprox(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTApp{}".format(nameBase),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainHF(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTHF{}".format(nameBase),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainErr(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTErr{}".format(nameBase),
forceNewFile = False, folder = True)
approx.outParaviewTimeDomainRes(ktar, omega = 2. * np.pi * ktar,
filename = "out/outTRes{}".format(nameBase),
forceNewFile = False, folder = True)
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
poles = approx.getPoles()
print('Poles:', poles)
if sol == "sweep":
k0s = np.linspace(k0s[0], k0s[1], 100)
kl, kr = min(k0s), max(k0s)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
kl, kr = np.real(kl), np.real(kr)
from matplotlib import pyplot as plt
normApp = np.zeros(len(k0s))
norm = np.zeros_like(normApp)
err = np.zeros_like(normApp)
res = np.zeros_like(normApp)
# errApp = np.zeros_like(normApp)
fNorm = approx.normRHS(k0s[0])
for j in range(len(k0s)):
normApp[j] = approx.normApprox(k0s[j])
norm[j] = approx.normHF(k0s[j])
err[j] = approx.normErr(k0s[j]) / norm[j]
res[j] = approx.normRes(k0s[j]) / fNorm
# errApp[j] = res[j] / np.min(np.abs(k0s[j] - poles))
# errApp *= np.mean(err) / np.mean(errApp)
plt.figure()
plt.semilogy(k0s, norm)
plt.semilogy(k0s, normApp, '--')
plt.semilogy(np.real(approx.mus),
1.05*np.max(norm)*np.ones_like(approx.mus, dtype = float),
'rx')
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, res)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(k0s, err)
# plt.semilogy(k0s, errApp)
plt.xlim([kl, kr])
plt.grid()
plt.show()
plt.close()
polesApp = approx.getPoles()
mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr)
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()
diff --git a/examples/from_papers/pod_membrane_centered.py b/examples/from_papers/pod_membrane_centered.py
index 965c695..18f3748 100644
--- a/examples/from_papers/pod_membrane_centered.py
+++ b/examples/from_papers/pod_membrane_centered.py
@@ -1,96 +1,72 @@
import fenics as fen
import numpy as np
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
from rrompy.reduction_methods.centered import RationalPade as TP
-from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
verb = 0
-test = "poles"
-test = "error"
-
k0 = 10
ktars = np.linspace(78**.5, 122**.5, 50)
def boundaryNeumann(x, on_boundary):
return on_boundary and x[1] > .25 and x[0] > 0.995 and x[0] < 1.005
meshname = '../data/mesh/crack_coarse.xml'
#meshname = '../data/mesh/crack_fine.xml'
mesh = fen.Mesh(meshname)
x, y = fen.SpatialCoordinate(mesh)[:]
x0, y0 = .5, .5
Rr, Ri = .1, .1
forcingTerm = fen.exp(- ((x - x0)**2 + (y - y0)**2) / 2 / Rr**2)
solver = HPE(verbosity = verb)
solver.omega = np.real(k0)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.forcingTerm = forcingTerm
solver.NeumannBoundary = boundaryNeumann
solver.DirichletBoundary = 'rest'
-if test == "poles":
- appPoles = {}
- Emax = 13
- params = {'N':6, 'M':0, 'E':6, 'sampleType':'Arnoldi',
- 'POD':True}
+appPoles = {}
+Emax = 13
+params = {'N':6, 'M':0, 'E':6, 'sampleType':'Arnoldi',
+ 'POD':True}
+
+approxPade = TP(solver, mu0 = k0, approxParameters = params,
+ verbosity = verb)
+for E in range(6, Emax + 1):
+ approxPade.E = E
+ appPoles[E] = np.sort(approxPade.getPoles())
- approxPade = TP(solver, mu0 = k0, approxParameters = params,
- verbosity = verb)
- for E in range(6, Emax + 1):
- approxPade.E = E
- appPoles[E] = np.sort(approxPade.getPoles())
+a = fen.dot(fen.grad(solver.u), fen.grad(solver.v)) * fen.dx
+A = fen.assemble(a)
+fen.DirichletBC(solver.V, fen.Constant(0.),
+ solver.DirichletBoundary).apply(A)
+AMat = fen.as_backend_type(A).mat()
+Ar, Ac, Av = AMat.getValuesCSR()
+import scipy.sparse as scsp
+A = scsp.csr_matrix((Av, Ac, Ar), shape = AMat.size)
- a = fen.dot(fen.grad(solver.u), fen.grad(solver.v)) * fen.dx
- A = fen.assemble(a)
- fen.DirichletBC(solver.V, fen.Constant(0.),
- solver.DirichletBoundary).apply(A)
- AMat = fen.as_backend_type(A).mat()
- Ar, Ac, Av = AMat.getValuesCSR()
- import scipy.sparse as scsp
- A = scsp.csr_matrix((Av, Ac, Ar), shape = AMat.size)
+m = fen.dot(solver.u, solver.v) * fen.dx
+M = fen.assemble(m)
+fen.DirichletBC(solver.V, fen.Constant(0.),
+ solver.DirichletBoundary).apply(M)
+MMat = fen.as_backend_type(M).mat()
+Mr, Mc, Mv = MMat.getValuesCSR()
+import scipy.sparse as scsp
+M = scsp.csr_matrix((Mv, Mc, Mr), shape = MMat.size)
- m = fen.dot(solver.u, solver.v) * fen.dx
- M = fen.assemble(m)
- fen.DirichletBC(solver.V, fen.Constant(0.),
- solver.DirichletBoundary).apply(M)
- MMat = fen.as_backend_type(M).mat()
- Mr, Mc, Mv = MMat.getValuesCSR()
- import scipy.sparse as scsp
- M = scsp.csr_matrix((Mv, Mc, Mr), shape = MMat.size)
+poles = scsp.linalg.eigs(A, k = 7, M = M, sigma = 100.,
+ return_eigenvectors = False)
+II = np.argsort(np.abs(poles - k0))
+poles = poles[II]
+print('Exact', end = ': ')
+[print('{},{}'.format(np.real(x), np.imag(x)), end = ',') for x in poles]
+print()
- poles = scsp.linalg.eigs(A, k = 7, M = M, sigma = 100.,
- return_eigenvectors = False)
- II = np.argsort(np.abs(poles - k0))
- poles = poles[II]
- print('Exact', end = ': ')
- [print('{},{}'.format(np.real(x), np.imag(x)), end = ',') for x in poles]
+for E in range(6, Emax + 1):
+ print(E, end = ': ')
+ [print('{},{}'.format(np.real(x), np.imag(x)), end = ',')\
+ for x in np.sort(appPoles[E])]
print()
-
- for E in range(6, Emax + 1):
- print(E, end = ': ')
- [print('{},{}'.format(np.real(x), np.imag(x)), end = ',')\
- for x in np.sort(appPoles[E])]
- print()
-
-elif test == "error":
- M0 = 5
- Emax = 8
- params = {'E':Emax, 'sampleType':'Arnoldi', 'POD':True}
- paramsSetsPade = [None] * (Emax - M0 + 1)
- for M in range(M0, Emax + 1):
- paramsSetsPade[M - M0] = {'N':M0 + 1, 'M':M, 'E':max(M, M0 + 1)}
- approxPade = TP(solver, mu0 = k0, approxParameters = params,
- verbosity = verb)
-
- sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
- sweeper.ROMEngine = approxPade
- sweeper.params = paramsSetsPade
- filenamePade = sweeper.sweep('membrane_error.dat',
- outputs = 'ALL')
- sweeper.plot(filenamePade, ['muRe'], ['normHF', 'normApprox'], ['M'],
- onePlot = True)
- sweeper.plot(filenamePade, ['muRe'], ['normErr'], ['M'])
diff --git a/rrompy/utilities/parameter_sweeper/__init__.py b/rrompy/utilities/parameter_sweeper/__init__.py
deleted file mode 100644
index 03ca6b3..0000000
--- a/rrompy/utilities/parameter_sweeper/__init__.py
+++ /dev/null
@@ -1,25 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from .parameter_sweeper import ParameterSweeper
-
-__all__ = [
- 'ParameterSweeper'
- ]
-
-
diff --git a/rrompy/utilities/parameter_sweeper/parameter_sweeper.py b/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
deleted file mode 100644
index d82cf6c..0000000
--- a/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
+++ /dev/null
@@ -1,503 +0,0 @@
-# Copyright (C) 2018 by the RROMPy authors
-#
-# This file is part of RROMPy.
-#
-# RROMPy is free software: you can redistribute it and/or modify
-# it under the terms of the GNU Lesser General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# RROMPy is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public License
-# along with RROMPy. If not, see .
-#
-
-from copy import copy
-import itertools
-import csv
-import numpy as np
-from matplotlib import pyplot as plt
-from rrompy.utilities.base.types import Np1D, DictAny, List, ROMEng
-from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
-from rrompy.utilities.exception_manager import RROMPyWarning
-from rrompy.utilities.exception_manager import RROMPyException
-
-__all__ = ['ParameterSweeper']
-
-def C2R2csv(x):
- x = np.ravel(x)
- y = np.concatenate((np.real(x), np.imag(x)))
- z = np.ravel(np.reshape(y, [2, np.size(x)]).T)
- return np.array2string(z, separator = '_', suppress_small = False,
- max_line_width = np.inf, sign = '+',
- formatter = {'all' : lambda x : "{:.15E}".format(x)}
- )[1 : -1]
-
-class ParameterSweeper:
- """
- ROM approximant parameter sweeper.
-
- Args:
- ROMEngine(optional): Generic approximant class. Defaults to None.
- mutars(optional): Array of parameter values to sweep. Defaults to empty
- array.
- params(optional): List of parameter settings (each as a dict) to
- explore. Defaults to single empty set.
- mostExpensive(optional): String containing label of most expensive
- step, to be executed fewer times. Allowed options are 'HF' and
- 'Approx'. Defaults to 'HF'.
-
- Attributes:
- ROMEngine: Generic approximant class.
- mutars: Array of parameter values to sweep.
- params: List of parameter settings (each as a dict) to explore.
- mostExpensive: String containing label of most expensive step, to be
- executed fewer times.
- """
-
- allowedOutputsStandard = ["normHF", "normApprox", "normRes", "normResRel",
- "normErr", "normErrRel"]
- allowedOutputs = allowedOutputsStandard + ["HFFunc", "ApproxFunc",
- "ErrFunc", "ErrFuncRel"]
- allowedOutputsFull = allowedOutputs + ["poles"]
-
- def __init__(self, ROMEngine : ROMEng = None, mutars : Np1D = np.array([]),
- params : List[DictAny] = [{}], mostExpensive : str = "HF"):
- self.ROMEngine = ROMEngine
- self.mutars = mutars
- self.params = params
- self.mostExpensive = mostExpensive
-
- def name(self) -> str:
- return self.__class__.__name__
-
- def __str__(self) -> str:
- return self.name()
-
- def __repr__(self) -> str:
- return self.__str__() + " at " + hex(id(self))
-
- @property
- def mostExpensive(self):
- """Value of mostExpensive."""
- return self._mostExpensive
- @mostExpensive.setter
- def mostExpensive(self, mostExpensive:str):
- mostExpensive = mostExpensive.upper()
- if mostExpensive not in ["HF", "APPROX"]:
- RROMPyWarning(("Value of mostExpensive not recognized. Overriding "
- "to 'APPROX'."))
- mostExpensive = "APPROX"
- self._mostExpensive = mostExpensive
-
- def checkValues(self) -> bool:
- """Check if sweep can be performed."""
- if self.ROMEngine is None:
- raise RROMPyException("ROMEngine is missing. Aborting.")
- if len(self.mutars) == 0:
- raise RROMPyException("Empty target parameter vector. Aborting.")
- if len(self.params) == 0:
- raise RROMPyException("Empty method parameters vector. Aborting.")
-
- def sweep(self, filename : str = "out.dat", outputs : List[str] = [],
- verbose : int = 10, timestamp : bool = True):
- self.checkValues()
- try:
- if outputs.upper() == "ALL":
- outputs = self.allowedOutputsFull
- except:
- if len(outputs) == 0:
- outputs = self.allowedOutputsStandard
- outputs = purgeList(outputs, self.allowedOutputsFull,
- listname = self.name() + ".outputs",
- baselevel = 1)
- poles = ("poles" in outputs)
- if len(outputs) == 0:
- raise RROMPyException("Empty outputs. Aborting.")
- outParList = self.ROMEngine.parameterList
- Nparams = len(self.params)
- if poles: polesCheckList = []
- allowedParams = self.ROMEngine.parameterList
-
- dotPos = filename.rfind('.')
- if dotPos in [-1, len(filename) - 1]:
- filename = getNewFilename(filename[:dotPos])
- else:
- filename = getNewFilename(filename[:dotPos], filename[dotPos + 1:])
-
- append_write = "w"
- initial_row = (outParList + ["muRe", "muIm"]
- + [x for x in self.allowedOutputs if x in outputs]
- + ["type"] + ["poles"] * poles)
- with open(filename, append_write, buffering = 1) as fout:
- writer = csv.writer(fout, delimiter=",")
- writer.writerow(initial_row)
-
- if self.mostExpensive == "HF":
- outerSet = self.mutars
- innerSet = self.params
- elif self.mostExpensive == "APPROX":
- outerSet = self.params
- innerSet = self.mutars
-
- for outerIdx, outerPar in enumerate(outerSet):
- if self.mostExpensive == "HF":
- i, mutar = outerIdx, outerPar
- elif self.mostExpensive == "APPROX":
- j, par = outerIdx, outerPar
- self.ROMEngine.approxParameters = {k: par[k] for k in\
- par.keys() & allowedParams}
- self.ROMEngine.setupApprox()
-
- for innerIdx, innerPar in enumerate(innerSet):
- if self.mostExpensive == "APPROX":
- i, mutar = innerIdx, innerPar
- elif self.mostExpensive == "HF":
- j, par = innerIdx, innerPar
- self.ROMEngine.approxParameters = {k: par[k] for k in\
- par.keys() & allowedParams}
- self.ROMEngine.setupApprox()
-
- if verbose >= 5:
- verbtype = "MAIN"
- if innerIdx == 0:
- verbtype = "INIT"
- verbosityDepth(verbtype, ("Set {}/{}\tmu_{} = "
- "{:.10f}").format(j + 1,
- Nparams, i,
- mutar),
- timestamp = timestamp)
-
- outData = []
- if "normHF" in outputs:
- valNorm = self.ROMEngine.normHF(mutar)
- outData = outData + [valNorm]
- if "normApprox" in outputs:
- val = self.ROMEngine.normApprox(mutar)
- outData = outData + [val]
- if "normRes" in outputs:
- valNRes = self.ROMEngine.normRes(mutar)
- outData = outData + [valNRes]
- if "normResRel" in outputs:
- if "normRes" not in outputs:
- valNRes = self.ROMEngine.normRes(mutar)
- val = self.ROMEngine.normRHS(mutar)
- outData = outData + [valNRes / val]
- if "normErr" in outputs:
- valNErr = self.ROMEngine.normErr(mutar)
- outData = outData + [valNErr]
- if "normErrRel" in outputs:
- if "normHF" not in outputs:
- valNorm = self.ROMEngine.normHF(mutar)
- if "normErr" not in outputs:
- valNErr = self.ROMEngine.normErr(mutar)
- outData = outData + [valNErr / valNorm]
- if "HFFunc" in outputs:
- valFunc = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getHF(mutar))
- outData = outData + [valFunc]
- if "ApproxFunc" in outputs:
- valFApp = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getApprox(mutar))
- outData = outData + [valFApp]
- if "ErrFunc" in outputs:
- if "HFFunc" not in outputs:
- valFunc = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getHF(mutar))
- if "ApproxFunc" not in outputs:
- valFApp = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getApprox(mutar))
- valFErr = np.abs(valFApp - valFunc)
- outData = outData + [valFErr]
- if "ErrFuncRel" in outputs:
- if not ("HFFunc" in outputs or "ErrFunc" in outputs):
- valFunc = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getHF(mutar))
- if not ("AppFunc" in outputs or "ErrFunc" in outputs):
- valFApp = self.ROMEngine.HFEngine.functional(
- self.ROMEngine.getApprox(mutar))
- val = np.nan
- if not np.isclose(valFunc, 0.):
- val = valFApp / valFunc
- outData = outData + [val]
- writeData = []
- for parn in outParList:
- writeData = (writeData
- + [self.ROMEngine.approxParameters[parn]])
- writeData = (writeData + [mutar.real, mutar.imag]
- + outData + [self.ROMEngine.name()])
- if poles:
- if j not in polesCheckList:
- polesCheckList += [j]
- writeData = writeData + [C2R2csv(
- self.ROMEngine.getPoles())]
- else:
- writeData = writeData + [""]
- writer.writerow(str(x) for x in writeData)
- if verbose >= 5:
- if self.mostExpensive == "APPROX":
- out = "Set {}/{}\tdone.".format(j + 1, Nparams)
- elif self.mostExpensive == "HF":
- out = "Point mu_{} = {:.10f}\tdone.".format(i, mutar)
- verbosityDepth("DEL", out, timestamp = timestamp)
- self.filename = filename
- return self.filename
-
- def read(self, filename:str, restrictions : DictAny = {},
- outputs : List[str] = []) -> DictAny:
- """
- Execute a query on a custom format CSV.
-
- Args:
- filename: CSV filename.
- restrictions(optional): Parameter configurations to output.
- Defaults to empty dictionary, i.e. output all.
- outputs(optional): Values to output. Defaults to empty list, i.e.
- no output.
-
- Returns:
- Dictionary of desired results, with a key for each entry of
- outputs, and a numpy 1D array as corresponding value.
- """
- with open(filename, 'r') as f:
- reader = csv.reader(f, delimiter=',')
- header = next(reader)
- restrIndices, outputIndices, outputData = {}, {}, {}
- for key in restrictions.keys():
- try:
- restrIndices[key] = header.index(key)
- if not isinstance(restrictions[key], list):
- restrictions[key] = [restrictions[key]]
- restrictions[key] = copy(restrictions[key])
- except:
- RROMPyWarning("Ignoring key {} from restrictions.".format(
- key))
- for key in outputs:
- try:
- outputIndices[key] = header.index(key)
- outputData[key] = np.array([])
- except:
- RROMPyWarning("Ignoring key {} from outputs.".format(key))
-
- for row in reader:
- restrTrue = True
- for key in restrictions.keys():
- if row[restrIndices[key]] == restrictions[key]:
- continue
- try:
- if np.any(np.isclose(float(row[restrIndices[key]]),
- [float(x) for x in restrictions[key]])):
- continue
- except: pass
- restrTrue = False
- if restrTrue:
- for key in outputIndices.keys():
- try:
- val = row[outputIndices[key]]
- val = float(val)
- finally:
- outputData[key] = np.append(outputData[key], val)
- return outputData
-
- def plot(self, filename:str, xs:List[str], ys:List[str], zs:List[str],
- onePlot : bool = False, save : str = None,
- saveFormat : str = "eps", saveDPI : int = 100, **figspecs):
- """
- Perform plots from data in filename.
-
- Args:
- filename: CSV filename.
- xs: Values to put on x axes.
- ys: Values to put on y axes.
- zs: Meta-values for constraints.
- onePlot(optional): Whether to create a single figure per x.
- Defaults to False.
- save(optional): Where to save plot(s). Defaults to None, i.e. no
- saving.
- saveFormat(optional): Format for saved plot(s). Defaults to "eps".
- saveDPI(optional): DPI for saved plot(s). Defaults to 100.
- figspecs(optional key args): Optional arguments for matplotlib
- figure creation.
- """
- if save is not None:
- save = save.strip()
- zsVals = self.read(filename, outputs = zs)
- zs = list(zsVals.keys())
- zss = None
- for key in zs:
- vals = np.unique(zsVals[key])
- if zss is None:
- zss = copy(vals)
- else:
- zss = list(itertools.product(zss, vals))
- lzs = len(zs)
- for z in zss:
- if lzs <= 1:
- constr = {zs[0] : z}
- else:
- constr = {zs[j] : z[j] for j in range(len(zs))}
- data = self.read(filename, restrictions = constr, outputs = xs+ys)
- if onePlot:
- for x in xs:
- xVals = data[x]
- p = plt.figure(**figspecs)
- logScale = False
- for y in ys:
- yVals = data[y]
- label = '{} vs {} for {}'.format(y, x, constr)
- if np.min(yVals) <= - np.finfo(float).eps:
- plt.plot(xVals, yVals, label = label)
- else:
- plt.plot(xVals, yVals, label = label)
- if np.log10(np.max(yVals) / np.min(yVals)) > 1.:
- logScale = True
- if logScale:
- ax = p.get_axes()[0]
- ax.set_yscale('log')
- plt.legend()
- plt.grid()
- if save is not None:
- prefix = "{}_{}_vs_{}_{}".format(save, ys, x, constr)
- plt.savefig(getNewFilename(prefix, saveFormat),
- format = saveFormat, dpi = saveDPI)
- plt.show()
- plt.close()
- else:
- for x, y in itertools.product(xs, ys):
- xVals, yVals = data[x], data[y]
- label = '{} vs {} for {}'.format(y, x, constr)
- p = plt.figure(**figspecs)
- if np.min(yVals) <= - np.finfo(float).eps:
- plt.plot(xVals, yVals, label = label)
- else:
- plt.plot(xVals, yVals, label = label)
- if np.log10(np.max(yVals) / np.min(yVals)) > 1.:
- ax = p.get_axes()[0]
- ax.set_yscale('log')
- plt.legend()
- plt.grid()
- if save is not None:
- prefix = "{}_{}_vs_{}_{}".format(save, y, x, constr)
- plt.savefig(getNewFilename(prefix, saveFormat),
- format = saveFormat, dpi = saveDPI)
- plt.show()
- plt.close()
-
- def plotCompare(self, filenames:List[str], xs:List[str], ys:List[str],
- zs:List[str], onePlot : bool = False, save : str = None,
- ylims : dict = None, saveFormat : str = "eps",
- saveDPI : int = 100, labels : List[str] = None,
- **figspecs):
- """
- Perform plots from data in filename1 and filename2.
-
- Args:
- filenames: CSV filenames.
- xs: Values to put on x axes.
- ys: Values to put on y axes.
- zs: Meta-values for constraints.
- onePlot(optional): Whether to create a single figure per x.
- Defaults to False.
- save(optional): Where to save plot(s). Defaults to None, i.e. no
- saving.
- clip(optional): Custom y axis limits. If None, automatic values are
- kept. Defaults to None.
- saveFormat(optional): Format for saved plot(s). Defaults to "eps".
- saveDPI(optional): DPI for saved plot(s). Defaults to 100.
- labels: Label for each dataset.
- figspecs(optional key args): Optional arguments for matplotlib
- figure creation.
- """
- nfiles = len(filenames)
- if save is not None:
- save = save.strip()
- if labels is None:
- labels = ["{}".format(j + 1) for j in range(nfiles)]
- zsVals = self.read(filenames[0], outputs = zs)
- zs = list(zsVals.keys())
- zss = None
- for key in zs:
- vals = np.unique(zsVals[key])
- if zss is None:
- zss = copy(vals)
- else:
- zss = list(itertools.product(zss, vals))
- lzs = len(zs)
- for z in zss:
- if lzs <= 1:
- constr = {zs[0] : z}
- else:
- constr = {zs[j] : z[j] for j in range(len(zs))}
- data = [None] * nfiles
- for j in range(nfiles):
- data[j] = self.read(filenames[j], restrictions = constr,
- outputs = xs + ys)
- if onePlot:
- for x in xs:
- xVals = [None] * nfiles
- for j in range(nfiles):
- try:
- xVals[j] = data[j][x]
- except:
- pass
- p = plt.figure(**figspecs)
- logScale = False
- for y in ys:
- for j in range(nfiles):
- try:
- yVals = data[j][y]
- except:
- pass
- l = '{} vs {} for {}, {}'.format(y, x, constr,
- labels[j])
- if np.min(yVals) <= - np.finfo(float).eps:
- plt.plot(xVals[j], yVals, label = l)
- else:
- plt.plot(xVals[j], yVals, label = l)
- if np.log10(np.max(yVals)/np.min(yVals)) > 1.:
- logScale = True
- if logScale:
- ax = p.get_axes()[0]
- ax.set_yscale('log')
- if ylims is not None:
- plt.ylim(**ylims)
- plt.legend()
- plt.grid()
- if save is not None:
- prefix = "{}_{}_vs_{}_{}".format(save, ys, x, constr)
- plt.savefig(getNewFilename(prefix, saveFormat),
- format = saveFormat, dpi = saveDPI)
- plt.show()
- plt.close()
- else:
- for x, y in itertools.product(xs, ys):
- p = plt.figure(**figspecs)
- logScale = False
- for j in range(nfiles):
- xVals, yVals = data[j][x], data[j][y]
- l = '{} vs {} for {}, {}'.format(y, x, constr,
- labels[j])
- if np.min(yVals) <= - np.finfo(float).eps:
- plt.plot(xVals, yVals, label = l)
- else:
- plt.plot(xVals, yVals, label = l)
- if np.log10(np.max(yVals)/np.min(yVals)) > 1.:
- logScale = True
- if logScale:
- ax = p.get_axes()[0]
- ax.set_yscale('log')
- if ylims is not None:
- plt.ylim(**ylims)
- plt.legend()
- plt.grid()
- if save is not None:
- prefix = "{}_{}_vs_{}_{}".format(save, y, x, constr)
- plt.savefig(getNewFilename(prefix, saveFormat),
- format = saveFormat, dpi = saveDPI)
- plt.show()
- plt.close()
-