diff --git a/rrompy/hfengines/base/hfengine_base.py b/rrompy/hfengines/base/hfengine_base.py index 84e099a..93976ae 100644 --- a/rrompy/hfengines/base/hfengine_base.py +++ b/rrompy/hfengines/base/hfengine_base.py @@ -1,271 +1,272 @@ # 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 abc import abstractmethod import numpy as np import scipy.sparse as scsp from numbers import Number from copy import copy as softcopy from rrompy.utilities.base.decorators import nonaffine_construct from rrompy.utilities.base.types import (Np1D, Np2D, List, DictAny, paramVal, paramList, sampList) from rrompy.utilities.numerical import solve as tsolve, dot, customPInv from rrompy.utilities.expression import expressionEvaluator from rrompy.utilities.exception_manager import RROMPyAssert +from rrompy.sampling.sample_list import sampleList from rrompy.parameter import (checkParameter, checkParameterList, - parameterMap as pMap) + parameterList, parameterMap as pMap) from rrompy.solver.linear_solver import setupSolver from rrompy.utilities.parallel import (poolRank, masterCore, listScatter, matrixGatherv, isend, recv) __all__ = ['HFEngineBase'] class HFEngineBase: """Generic solver for parametric problems.""" def __init__(self, verbosity : int = 10, timestamp : bool = True): self.verbosity = verbosity self.timestamp = timestamp self.setSolver("SPSOLVE", {"use_umfpack" : False}) self.npar = 0 self._C = None self.outputNormMatrix = 1. 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)) def __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] def __deepcopy__(self, memo): return softcopy(self) @property def npar(self): """Value of npar.""" return self._npar @npar.setter def npar(self, npar): nparOld = self._npar if hasattr(self, "_npar") else -1 if npar != nparOld: self.parameterMap = pMap(1., npar) self._npar = npar @property def spacedim(self): return 1 def checkParameter(self, mu:paramVal) -> paramVal: muP = checkParameter(mu, self.npar) if self.npar == 0: muP.reset((1, 0), muP.dtype) return muP def checkParameterList(self, mu:paramList, check_if_single : bool = False) -> paramList: muL = checkParameterList(mu, self.npar, check_if_single) return muL def mapParameterList(self, mu:paramList, direct : str = "F", idx : List[int] = None) -> paramList: if idx is None: idx = np.arange(self.npar) muMapped = checkParameterList(mu, len(idx)) for j, d in enumerate(idx): muMapped.data[:, j] = expressionEvaluator( self.parameterMap[direct][d], muMapped(j)).flatten() return muMapped def buildEnergyNormForm(self): """ Build sparse matrix (in CSR format) representative of scalar product. """ self.energyNormMatrix = 1. def buildEnergyNormDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product without duality. """ self.energyNormDualMatrix = 1. def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False, dual : bool = False, is_state : bool = True) -> Np2D: """Scalar product.""" if is_state or self.isCEye: if dual: if not hasattr(self, "energyNormDualMatrix"): self.buildEnergyNormDualForm() energyMat = self.energyNormDualMatrix else: if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() energyMat = self.energyNormMatrix else: energyMat = self.outputNormMatrix - if not isinstance(u, (np.ndarray,)): u = u.data - if not isinstance(v, (np.ndarray,)): v = v.data + if isinstance(u, (parameterList, sampleList)): u = u.data + if isinstance(v, (parameterList, sampleList)): v = v.data if onlyDiag: return np.sum(dot(energyMat, u) * v.conj(), axis = 0) return dot(dot(energyMat, u).T, v.conj()).T def norm(self, u:Np2D, dual : bool = False, is_state : bool = True) -> Np1D: return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual, is_state = is_state)) ** .5 def baselineA(self): """Return 0 of shape consistent with operator of linear system.""" if (hasattr(self, "As") and hasattr(self.As, "__len__") and self.As[0] is not None): d = self.As[0].shape[0] else: d = self.spacedim return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)), shape = (d, d), dtype = np.complex) def baselineb(self): """Return 0 of shape consistent with RHS of linear system.""" return np.zeros(self.spacedim, dtype = np.complex) @nonaffine_construct @abstractmethod def A(self, mu : paramVal = [], der : List[int] = 0) -> Np2D: """ Assemble terms of operator of linear system and return it (or its derivative) at a given parameter. """ return @nonaffine_construct @abstractmethod def b(self, mu : paramVal = [], der : List[int] = 0) -> Np1D: """ Assemble terms of RHS of linear system and return it (or its derivative) at a given parameter. """ return @property def C(self): """Value of C.""" if self._C is None: self._C = 1. return self._C @property def isCEye(self): return isinstance(self.C, Number) def applyC(self, u:sampList): """Apply LHS of linear system.""" return dot(self.C, u) def applyCpInv(self, u:sampList): """Apply pseudoinverse of LHS of linear system.""" return dot(customPInv(self.C), u) def setSolver(self, solverType:str, solverArgs : DictAny = {}): """Choose solver type and parameters.""" self._solver, self._solverArgs = setupSolver(solverType, solverArgs) def solve(self, mu : paramList = [], RHS : sampList = None, return_state : bool = False) -> sampList: """ Find solution of linear system. Args: mu: parameter value. RHS: RHS of linear system. If None, defaults to that of parametric system. Defaults to None. return_state: whether to return state before multiplication by c. Defaults to False. """ from rrompy.sampling import sampleList, emptySampleList if mu == []: mu = self.mu0 mu = self.checkParameterList(mu) if len(mu) == 0: return emptySampleList() mu, idx, sizes = listScatter(mu, return_sizes = True) mu = self.checkParameterList(mu) req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(mu) == 0: uL, uT = recv(source = 0, tag = poolRank()) sol = np.empty((uL, 0), dtype = uT) else: if RHS is None: RHS = sampleList([self.b(m) for m in mu]) else: RHS = sampleList(RHS) if len(RHS) > 1: RHS = sampleList([RHS[i] for i in idx]) mult = 0 if len(RHS) == 1 else 1 RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size") for j, mj in enumerate(mu): u = tsolve(self.A(mj), RHS[mult * j], self._solver, self._solverArgs) if j == 0: sol = np.empty((len(u), len(mu)), dtype = u.dtype) if masterCore(): for dest in emptyCores: req += [isend((len(u), u.dtype), dest = dest, tag = dest)] sol[:, j] = u if not return_state: sol = self.applyC(sol) for r in req: r.wait() return sampleList(matrixGatherv(sol, sizes)) def residual(self, mu : paramList = [], u : sampList = None, post_c : bool = True) -> sampList: """ Find residual of linear system for given approximate solution. Args: mu: parameter value. u: numpy complex array with function dofs. If None, set to 0. post_c: whether to post-process using c. Defaults to True. """ from rrompy.sampling import sampleList, emptySampleList if mu == []: mu = self.mu0 mu = self.checkParameterList(mu) if len(mu) == 0: return emptySampleList() mu, idx, sizes = listScatter(mu, return_sizes = True) mu = self.checkParameterList(mu) req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(mu) == 0: uL, uT = recv(source = 0, tag = poolRank()) res = np.empty((uL, 0), dtype = uT) else: v = sampleList(np.zeros((self.spacedim, len(mu)))) if u is not None: u = sampleList(u) v = v + sampleList([u[i] for i in idx]) for j, (mj, vj) in enumerate(zip(mu, v)): r = self.b(mj) - dot(self.A(mj), vj) if j == 0: res = np.empty((len(r), len(mu)), dtype = r.dtype) if masterCore(): for dest in emptyCores: req += [isend((len(r), r.dtype), dest = dest, tag = dest)] res[:, j] = r if post_c: res = self.applyC(res) for r in req: r.wait() return sampleList(matrixGatherv(res, sizes)) diff --git a/rrompy/hfengines/base/marginal_proxy_engine.py b/rrompy/hfengines/base/marginal_proxy_engine.py index 2e83dd9..4b4d90c 100644 --- a/rrompy/hfengines/base/marginal_proxy_engine.py +++ b/rrompy/hfengines/base/marginal_proxy_engine.py @@ -1,160 +1,158 @@ # 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 . # import inspect import numpy as np from copy import copy as softcopy from rrompy.utilities.base.types import Np1D, paramVal, paramList, HFEng from rrompy.utilities.base import freepar as fp from rrompy.utilities.base.decorators import (affine_construct, nonaffine_construct) from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter, checkParameterList __all__ = ['MarginalProxyEngine'] def MarginalProxyEngine(HFEngine:HFEng, marginalized:Np1D): Aaff = hasattr(HFEngine.A, "is_affine") and HFEngine.A.is_affine baff = hasattr(HFEngine.b, "is_affine") and HFEngine.b.is_affine if Aaff: if baff: return MarginalProxyEngineAffineAb(HFEngine, marginalized) return MarginalProxyEngineAffineA(HFEngine, marginalized) if baff: return MarginalProxyEngineAffineb(HFEngine, marginalized) return MarginalProxyEngineNonAffine(HFEngine, marginalized) class MarginalProxyEngineNonAffine: """ Marginalized should prescribe fixed value for the marginalized parameters and leave freepar/None elsewhere. """ _allowedMuDependencies = ["A", "b", "checkParameter", "checkParameterList", "_assembleObject", "solve", "residual"] def __init__(self, HFEngine:HFEng, marginalized:Np1D): self.baseHF = HFEngine self.marg = marginalized for name in HFEngine.__dir_base__(): att = getattr(HFEngine, name) if inspect.ismethod(att): attargs = inspect.getfullargspec(att).args if "mu" not in attargs: setattr(self.__class__, name, getattr(HFEngine, name)) else: if name not in self._allowedMuDependencies: raise RROMPyException(("Function {} depends on mu " "and was not accounted for. " "Must override.").format(name)) @property def affinePoly(self): return self.nparFixed == 0 and self.baseHF.affinePoly @property def freeLocations(self): return [x for x in range(self.baseHF.npar) if self.marg[x] == fp] @property def fixedLocations(self): return [x for x in range(self.baseHF.npar) if self.marg[x] != fp] @property def _freeLocationsInsert(self): return np.cumsum([m == fp for m in self.marg])[self.fixedLocations] @property def muFixed(self): muF = checkParameter([m for m in self.marg if m != fp]) if self.baseHF.npar - self.nparFree > 0: muF = muF[0] return muF @property def nparFree(self): """Value of nparFree.""" return len(self.freeLocations) @property def nparFixed(self): """Value of nparFixed.""" return len(self.fixedLocations) def name(self) -> str: return "{}-proxy for {}".format(self.freeLocations, self.baseHF.name()) def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] def __deepcopy__(self, memo): return softcopy(self) def completeMu(self, mu:paramVal): - mu = checkParameter(mu, self.nparFree) - return np.insert(mu.data, self._freeLocationsInsert, self.muFixed, - axis = 1) + mu = checkParameter(mu, self.nparFree, return_data = True) + return np.insert(mu, self._freeLocationsInsert, self.muFixed, axis = 1) def completeMuList(self, mu:paramList): - mu = checkParameterList(mu, self.nparFree) - return np.insert(mu.data, self._freeLocationsInsert, self.muFixed, - axis = 1) + mu = checkParameterList(mu, self.nparFree, return_data = True) + return np.insert(mu, self._freeLocationsInsert, self.muFixed, axis = 1) @nonaffine_construct def A(self, mu : paramVal = [], *args, **kwargs): return self.baseHF.A(self.completeMu(mu), *args, **kwargs) @nonaffine_construct def b(self, mu : paramVal = [], *args, **kwargs): return self.baseHF.b(self.completeMu(mu), *args, **kwargs) def checkParameter(self, mu : paramVal = [], *args, **kwargs): return self.baseHF.checkParameter(self.completeMu(mu), *args, **kwargs) def checkParameterList(self, mu : paramList = [], *args, **kwargs): return self.baseHF.checkParameterList(self.completeMuList(mu), *args, **kwargs) def _assembleObject(self, mu : paramVal = [], *args, **kwargs): return self.baseHF._assembleObject(self.completeMu(mu), *args, **kwargs) def solve(self, mu : paramList = [], *args, **kwargs): return self.baseHF.solve(self.completeMuList(mu), *args, **kwargs) def residual(self, mu : paramList = [], *args, **kwargs): return self.baseHF.residual(self.completeMuList(mu), *args, **kwargs) class MarginalProxyEngineAffineA(MarginalProxyEngineNonAffine): @affine_construct def A(self, mu : paramVal = [], *args, **kwargs): return self.baseHF.A(self.completeMu(mu), *args, **kwargs) class MarginalProxyEngineAffineb(MarginalProxyEngineNonAffine): @affine_construct def b(self, mu : paramVal = [], *args, **kwargs): return self.baseHF.b(self.completeMu(mu), *args, **kwargs) class MarginalProxyEngineAffineAb(MarginalProxyEngineAffineA, MarginalProxyEngineAffineb): pass diff --git a/rrompy/parameter/parameter_list.py b/rrompy/parameter/parameter_list.py index b61c89b..ec0ccb6 100644 --- a/rrompy/parameter/parameter_list.py +++ b/rrompy/parameter/parameter_list.py @@ -1,232 +1,234 @@ # 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 . # import numpy as np from itertools import product as iterprod from copy import deepcopy as copy from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.utilities.base.types import Np2D __all__ = ['parameterList', 'emptyParameterList', 'checkParameterList'] -def checkParameterList(mu, npar = None, check_if_single : bool = False): +def checkParameterList(mu, npar = None, check_if_single : bool = False, + return_data : bool = False): if not isinstance(mu, (parameterList,)): mu = parameterList(mu, npar) else: if npar is not None: RROMPyAssert(mu.shape[1], npar, "Number of parameters") mu = copy(mu) if npar == 0: mu.reset((1, 0), mu.dtype) + if return_data: mu = mu.data if check_if_single: return mu, len(mu) <= 1 return mu -def checkParameter(mu, npar = None): - muL, wasPar = checkParameterList(mu, npar, True) +def checkParameter(mu, npar = None, return_data : bool = False): + muL, wasPar = checkParameterList(mu, npar, True, return_data) if not wasPar: - muL, wasPar = checkParameterList([mu], npar, True) + muL, wasPar = checkParameterList([mu], npar, True, return_data) if not wasPar: raise RROMPyException(("Only single parameter allowed. No " "parameter lists here.")) return muL def emptyParameterList(): return parameterList([[]]) def addMemberFromNumpyArray(self, fieldName): def objFunc(self, other): if not isinstance(other, (self.__class__,)): other = parameterList(other, self.shape[1]) return parameterList(getattr(np.ndarray, fieldName)(self.data, other.data)) setattr(self.__class__, fieldName, objFunc) def objIFunc(self, other): self.data = getattr(self.__class__, fieldName)(self, other).data setattr(self.__class__, "__i" + fieldName[2:], objIFunc) class parameterList: __all__ += [pre + post for pre, post in iterprod(["__", "__i"], ["add__", "sub__", "mul__", "div__", "truediv__", "floordiv__", "pow__"])] def __init__(self, data:Np2D, lengthCheck : int = None): if not hasattr(data, "__len__"): data = [data] elif isinstance(data, (self.__class__,)): data = data.data elif isinstance(data, (tuple,)): data = list(data) if (isinstance(data, (list,)) and len(data) > 0 and isinstance(data[0], (tuple,))): data = [list(x) for x in data] self.data = np.array(data, ndmin = 1, copy = 1) if self.data.ndim == 1: self.data = self.data[:, None] if np.size(self.data) > 0: self.data = self.data.reshape((len(self), -1)) if self.shape[0] * self.shape[1] == 0: lenEff = 0 if lengthCheck is None else lengthCheck self.reset((0, lenEff), self.dtype) if lengthCheck is not None: if lengthCheck != 1 and self.shape == (lengthCheck, 1): self.data = self.data.T RROMPyAssert(self.shape[1], lengthCheck, "Number of parameters") for fieldName in ["__add__", "__sub__", "__mul__", "__div__", "__truediv__", "__floordiv__", "__pow__"]: addMemberFromNumpyArray(self, fieldName) def __len__(self): return self.shape[0] def __str__(self): if len(self) == 0: selfstr = "[]" elif len(self) <= 3: selfstr = "[{}]".format(" ".join([str(x) for x in self.data])) else: selfstr = "[{} ..({}).. {}]".format(self[0], len(self) - 2, self[-1]) return selfstr def __repr__(self): return repr(self.data) @property def shape(self): return self.data.shape @property def size(self): return self.data.size @property def re(self): return parameterList(np.real(self.data)) @property def im(self): return parameterList(np.imag(self.data)) @property def abs(self): return parameterList(np.abs(self.data)) @property def angle(self): return parameterList(np.angle(self.data)) @property def conj(self): return parameterList(np.conj(self.data)) @property def dtype(self): return self.data.dtype def __getitem__(self, key): return self.data[key] def __call__(self, key, idx = None): if idx is None: return self.data[:, key] return self[key, idx] def __setitem__(self, key, value): - if isinstance(key, (tuple, list,)): + if isinstance(key, (tuple, list, np.ndarray)): RROMPyAssert(len(key), len(value), "Slice length") for k, val in zip(key, value): self[k] = val else: self.data[key] = value def __eq__(self, other): if not hasattr(other, "shape") or self.shape != other.shape: return False if isinstance(other, self.__class__): other = other.data return np.allclose(self.data, other) def __contains__(self, item): return next((x for x in self if np.allclose(x[0], item)), -1) != -1 def __iter__(self): return iter([parameterList([x]) for x in self.data]) def __copy__(self): return parameterList(self.data) def __deepcopy__(self, memo): return parameterList(copy(self.data, memo)) def __neg__(self): return parameterList(-self.data) def __pos__(self): return copy(self) def reset(self, size, dtype = complex): self.data = np.empty(size, dtype = dtype) self.data[:] = np.nan def insert(self, items, idx = None): if isinstance(items, self.__class__): items = items.data else: items = np.array(items, ndmin = 2) if len(self) == 0: self.data = parameterList(items).data elif idx is None: self.data = np.append(self.data, items, axis = 0) else: self.data = np.insert(self.data, idx, items, axis = 0) def append(self, items): self.insert(items) def pop(self, idx = -1): self.data = np.delete(self.data, idx, axis = 0) def find(self, item): if len(self) == 0: return None return next((j for j in range(len(self)) if np.allclose(self[j], item)), None) def findall(self, item): if len(self) == 0: return [] return [j for j in range(len(self)) if np.allclose(self[j], item)] def sort(self, overwrite = False, *args, **kwargs): dataT = np.array([tuple(x[0]) for x in self], dtype = [(str(j), self.dtype) for j in range(self.shape[1])]) sortedP = parameterList([list(x) for x in np.sort(dataT, *args, **kwargs)]) if overwrite: self.data = sortedP.data return sortedP def unique(self, overwrite = False, *args, **kwargs): dataT = np.array([tuple(x[0]) for x in self], dtype = [(str(j), self.dtype) for j in range(self.shape[1])]) uniqueT = np.unique(dataT, *args, **kwargs) if isinstance(uniqueT, (tuple,)): extraT = uniqueT[1:] uniqueT = uniqueT[0] else: extraT = () uniqueP = parameterList([list(x) for x in uniqueT]) if overwrite: self.data = uniqueP.data uniqueP = (uniqueP,) + extraT if len(uniqueP) == 1: return uniqueP[0] return uniqueP diff --git a/rrompy/parameter/parameter_sampling/segment/quadrature_sampler_total.py b/rrompy/parameter/parameter_sampling/segment/quadrature_sampler_total.py index f759cce..bbf3fdd 100644 --- a/rrompy/parameter/parameter_sampling/segment/quadrature_sampler_total.py +++ b/rrompy/parameter/parameter_sampling/segment/quadrature_sampler_total.py @@ -1,62 +1,62 @@ # 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 . # import numpy as np from scipy.special import binom, factorial as fact from .quadrature_sampler import QuadratureSampler from rrompy.utilities.base.types import paramList from rrompy.utilities.numerical import lowDiscrepancy from rrompy.parameter import checkParameterList __all__ = ['QuadratureSamplerTotal'] class QuadratureSamplerTotal(QuadratureSampler): """ Generator of quadrature sample points for total degree polynomial computations. """ def generatePoints(self, n:int, reorder : bool = True) -> paramList: """Array of sample points.""" d = self.npar n1d = int((fact(d) * n) ** (1. / d)) while binom(n1d + d - 1, d) > n: n1d -= 1 x = super().generatePoints(n1d ** d, reorder = False) nTot = n1d ** d indicesBase = np.zeros(nTot, dtype = int) if reorder: idxBase = ([y + 1 for y in lowDiscrepancy(n1d - 1, inverse = True)] + [0]) else: idxBase = list(range(n1d)) linearIdxs = np.array(idxBase) nleft, nright = 1, nTot for j in range(d): nright //= n1d kronIn = np.repeat(linearIdxs, nright) indicesBase += np.tile(kronIn, nleft) nleft *= n1d keepIdxs = np.zeros(nTot, dtype = bool) keepIdxs[indicesBase < n1d] = True - xmat = x.data[keepIdxs, :] + xmat = x[keepIdxs] if reorder: nx = len(xmat) fejerOrdering = [nx - 1] + lowDiscrepancy(nx - 1, inverse = True) xmat = xmat[fejerOrdering, :] return checkParameterList(xmat, d) diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py index 7c5a803..ed6b372 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py @@ -1,115 +1,115 @@ # 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 deepcopy as copy from itertools import product import numpy as np from rrompy.parameter.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds, sparseMap) from rrompy.utilities.base.types import List, DictAny, paramList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['SparseGridSampler'] class SparseGridSampler(GenericSampler): """Generator of sparse grid sample points.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", parameterMap : DictAny = 1.): super().__init__(lims = lims, parameterMap = parameterMap) self.kind = kind self.reset() def __str__(self) -> str: return "{}[{}_{}]_{}".format(self.name(), self.lims[0], self.lims[1], self.kind) @property def npoints(self): """Number of points.""" return len(self.points) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in [sk.split("_")[2] + extra for sk, extra in product(sparsekinds, ["", "-NOBOUNDARY"])]: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() self._noBoundary = "NOBOUNDARY" in self._kind def reset(self): limsE = self.mapParameterList(self.lims) centerEff = .5 * (limsE[0] + limsE[1]) self.points = self.mapParameterList(centerEff, "B") self.depth = np.zeros((1, self.npar), dtype = int) self.deltadepth = np.zeros(1, dtype = int) def refine(self, active : List[int] = None) -> List[int]: if active is None: active = np.arange(self.npoints) - limsX = self.mapParameterList(self.lims).data.T + limsX = self.mapParameterList(self.lims) newIdxs = [] for act in active: point, dpt = self.points[act], self.depth[act] exhausted = False while not exhausted: ddp = self.deltadepth[act] for jdelta in range(self.npar): for signdelta in [-1., 1.]: Pointj = sparseMap( self.mapParameterList(point[jdelta], idx = [jdelta])(0, 0), - limsX[jdelta], self.kind, False) + limsX(jdelta), self.kind, False) if not self._noBoundary and dpt[jdelta] == 1: Centerj = sparseMap( self.mapParameterList(self.points[0][jdelta], idx = [jdelta])(0, 0), - limsX[jdelta], self.kind, False) + limsX(jdelta), self.kind, False) gradj = Pointj - Centerj if signdelta * gradj > 0: continue pointj = copy(point) Pointj = Pointj + .5 ** (dpt[jdelta] + ddp + self._noBoundary) * signdelta pointj[jdelta] = self.mapParameterList( - sparseMap(Pointj, limsX[jdelta], self.kind), + sparseMap(Pointj, limsX(jdelta), self.kind), "B", [jdelta])(0, 0) dist = np.sum(np.abs(self.points.data - pointj.reshape(1, -1)), axis = 1) samePt = np.where(np.isclose(dist, 0.))[0] if len(samePt) > 0: if samePt[0] in newIdxs: exhausted = True continue newIdxs += [self.npoints] self.points.append(pointj) self.depth = np.append(self.depth, [dpt], 0) self.depth[-1, jdelta] += 1 + ddp self.deltadepth = np.append(self.deltadepth, [0]) exhausted = True self.deltadepth[act] += 1 return newIdxs def generatePoints(self, n:int, reorder = None) -> paramList: if self.npoints > n: self.reset() idx = np.arange(self.npoints) while self.npoints < n: idx = self.refine(idx) return self.points diff --git a/rrompy/sampling/engines/pod_engine.py b/rrompy/sampling/engines/pod_engine.py index 1f146bd..c617af0 100644 --- a/rrompy/sampling/engines/pod_engine.py +++ b/rrompy/sampling/engines/pod_engine.py @@ -1,137 +1,137 @@ # 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 . # import numpy as np from copy import deepcopy as copy from warnings import catch_warnings from rrompy.utilities.base.types import Np1D, Np2D, Tuple, HFEng, sampList -from rrompy.utilities.numerical import dot from rrompy.sampling import sampleList __all__ = ['PODEngine'] class PODEngine: """ POD engine for general matrix orthogonalization. """ def __init__(self, HFEngine:HFEng): self.HFEngine = HFEngine 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)) def GS(self, a:Np1D, Q:sampList, n : int = -1, is_state : bool = True) -> Tuple[Np1D, Np1D, bool]: """ Compute 1 Gram-Schmidt step with given projector. Args: a: vector to be projected; Q: orthogonal projection matrix; n: number of columns of Q to be considered; is_state: whether state-norm should be used. Returns: Resulting normalized vector, coefficients of a wrt the updated basis, whether computation is ill-conditioned. """ if n == -1: n = Q.shape[1] r = np.zeros((n + 1,), dtype = Q.dtype) if n > 0: + from rrompy.utilities.numerical import dot Q = Q[: n] for j in range(2): # twice is enough! nu = self.HFEngine.innerProduct(a, Q, is_state = is_state) a = a - dot(Q, nu) r[:-1] = r[:-1] + nu.flatten() r[-1] = self.HFEngine.norm(a, is_state = is_state) ill_cond = False with catch_warnings(record = True) as w: snr = np.abs(r[-1]) / np.linalg.norm(r) if len(w) > 0 or np.isclose(snr, 0.): ill_cond = True r[-1] = 1. a = a / r[-1] return a, r, ill_cond def generalizedQR(self, A:sampList, Q0 : sampList = None, only_R : bool = False, genTrials : int = 10, is_state : bool = True) -> Tuple[sampList, Np2D]: """ Compute generalized QR decomposition of a matrix through Householder method; see Trefethen, IMA J.N.A., 2010. Args: A: matrix to be decomposed; Q0(optional): initial orthogonal guess for Q; defaults to random; only_R(optional): whether to skip reconstruction of Q; defaults to False. genTrials(optional): number of trials of generation of linearly independent vector; defaults to 10. is_state(optional): whether state-norm should be used; defaults to True. Returns: Resulting (orthogonal and )upper-triangular factor(s). """ Nh, N = A.shape B = copy(A) V = sampleList(np.zeros(A.shape, dtype = A.dtype)) R = np.zeros((N, N), dtype = A.dtype) - Q = copy(V) if Q0 is None else copy(Q0) + Q = copy(V) if Q0 is None else sampleList(Q0) for k in range(N): a = B[k] R[k, k] = self.HFEngine.norm(a, is_state = is_state) if Q0 is None and k < Nh: for _ in range(genTrials): Q[k], _, illC = self.GS(np.random.randn(Nh), Q, k, is_state) if not illC: break else: illC = k >= Nh if illC: - if Q0 is not None or k < Nh: Q.data[:, k] = 0. + if Q0 is not None or k < Nh: Q[k] = 0. else: alpha = self.HFEngine.innerProduct(a, Q[k], is_state = is_state) if np.isclose(np.abs(alpha), 0.): s = 1. else: s = - alpha / np.abs(alpha) - Q.data[:, k] *= s + Q[k] = s * Q[k] V[k], _, _ = self.GS(R[k, k] * Q[k] - a, Q, k, is_state) J = np.arange(k + 1, N) vtB = self.HFEngine.innerProduct(B[J], V[k], is_state = is_state) - B.data[:, J] -= 2 * np.outer(V[k], vtB) + B[J] = (B[J] - 2 * np.outer(V[k], vtB)).T if not illC: R[k, J] = self.HFEngine.innerProduct(B[J], Q[k], is_state = is_state) - B.data[:, J] -= np.outer(Q[k], R[k, J]) + B[J] = (B[J] - np.outer(Q[k], R[k, J])).T if only_R: return R for k in range(min(Nh, N) - 1, -1, -1): J = np.arange(k, min(Nh, N)) vtQ = self.HFEngine.innerProduct(Q[J], V[k], is_state = is_state) - Q.data[:, J] -= 2 * np.outer(V[k], vtQ) + Q[J] = (Q[J] - 2 * np.outer(V[k], vtQ)).T return Q, R diff --git a/rrompy/sampling/engines/sampling_engine_standard.py b/rrompy/sampling/engines/sampling_engine_standard.py index 4ee3be4..b663ba4 100644 --- a/rrompy/sampling/engines/sampling_engine_standard.py +++ b/rrompy/sampling/engines/sampling_engine_standard.py @@ -1,358 +1,358 @@ # 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 numbers import Number from copy import deepcopy as copy import numpy as np from warnings import catch_warnings from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, List, paramVal, Any, paramList, sampList, Tuple, TupleAny, DictAny, FigHandle) from rrompy.utilities.base.data_structures import getNewFilename from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert, RROMPy_READY, RROMPy_FRAGILE) -from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.base.pickle_utilities import pickleDump, pickleLoad from rrompy.parameter import (emptyParameterList, checkParameter, checkParameterList) from rrompy.sampling import sampleList, emptySampleList from rrompy.utilities.parallel import bcast, masterCore __all__ = ['SamplingEngineStandard'] class SamplingEngineStandard: def __init__(self, HFEngine:HFEng, sample_state : bool = False, verbosity : int = 10, timestamp : bool = True, scaleFactor : Np1D = None): self.sample_state = sample_state self.verbosity = verbosity self.timestamp = timestamp vbMng(self, "INIT", "Initializing sampling engine of type {}.".format(self.name()), 10) self.HFEngine = HFEngine vbMng(self, "DEL", "Done initializing sampling engine.", 10) self.scaleFactor = scaleFactor 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 HFEngine(self): """Value of HFEngine. Its assignment resets history.""" return self._HFEngine @HFEngine.setter def HFEngine(self, HFEngine): self._HFEngine = HFEngine self.resetHistory() @property def scaleFactor(self): """Value of scaleFactor.""" return self._scaleFactor @scaleFactor.setter def scaleFactor(self, scaleFactor): if scaleFactor is None: scaleFactor = 1. if not hasattr(scaleFactor, "__len__"): scaleFactor = [scaleFactor] self._scaleFactor = scaleFactor def scaleDer(self, derIdx:Np1D): if not isinstance(self.scaleFactor, Number): RROMPyAssert(len(derIdx), len(self.scaleFactor), "Number of derivative indices") with catch_warnings(record = True) as w: res = np.prod(np.power(self.scaleFactor, derIdx)) if len(w) == 0: return res raise RROMPyException(("Error in computing derivative scaling " "factor: {}".format(str(w[-1].message)))) @property def feature_keys(self) -> TupleAny: return ["mus", "samples", "nsamples", "_derIdxs"] @property def feature_vals(self) -> DictAny: return {"mus":self.mus, "samples":self.samples, "nsamples":self.nsamples, "_derIdxs":self._derIdxs, "_scaleFactor":self.scaleFactor} def _mergeFeature(self, name:str, value:Any): if name in ["mus", "samples"]: getattr(self, name).append(value) elif name == "nsamples": self.nsamples += value elif name == "_derIdxs": self._derIdxs += value else: raise RROMPyException(("Invalid key {} in sampling engine " "merge.".format(name))) def store(self, filenameBase : str = "sampling_engine", forceNewFile : bool = True, local : bool = False) -> str: """Store sampling engine to file.""" filename = None if masterCore(): vbMng(self, "INIT", "Storing sampling engine to file.", 20) if forceNewFile: filename = getNewFilename(filenameBase, "pkl") else: filename = "{}.pkl".format(filenameBase) pickleDump(self.feature_vals, filename) vbMng(self, "DEL", "Done storing engine.", 20) if local: return filename = bcast(filename) return filename def load(self, filename:str, merge : bool = False): """Load sampling engine from file.""" if isinstance(filename, (list, tuple,)): self.load(filename[0], merge) for filen in filename[1 :]: self.load(filen, True) return vbMng(self, "INIT", "Loading stored sampling engine from file.", 20) datadict = pickleLoad(filename) for key in datadict: if key in self.feature_keys: if merge and key != "_scaleFactor": self._mergeFeature(key, datadict[key]) else: setattr(self, key, datadict[key]) self._mode = RROMPy_FRAGILE vbMng(self, "DEL", "Done loading stored engine.", 20) @property def projectionMatrix(self) -> Np2D: return self.samples.data def resetHistory(self): self._mode = RROMPy_READY self.samples = emptySampleList() self.nsamples = 0 self.mus = emptyParameterList() self._derIdxs = [] def setsample(self, u:sampList, overwrite : bool = False): if overwrite: self.samples[self.nsamples] = u else: if self.nsamples == 0: self.samples = sampleList(u) else: self.samples.append(u) def popSample(self): if hasattr(self, "nsamples") and self.nsamples > 1: if self.samples.shape[1] > self.nsamples: RROMPyWarning(("More than 'nsamples' memory allocated for " "samples. Popping empty sample column.")) self.nsamples += 1 self.nsamples -= 1 self.samples.pop() self.mus.pop() else: self.resetHistory() def preallocateSamples(self, u:sampList, mu:paramVal, n:int): self._mode = RROMPy_READY self.samples.reset((u.shape[0], n), u.dtype) self.samples[0] = u mu = checkParameter(mu, self.HFEngine.npar) self.mus.reset((n, self.HFEngine.npar)) self.mus[0] = mu[0] def postprocessu(self, u:sampList, overwrite : bool = False): self.setsample(u, overwrite) def postprocessuBulk(self): pass def solveLS(self, mu : paramList = [], RHS : sampList = None) -> sampList: """ Solve linear system. Args: mu: Parameter value. Returns: Solution of system. """ mu = checkParameterList(mu, self.HFEngine.npar) vbMng(self, "INIT", "Solving HF model for mu = {}.".format(mu), 15) u = self.HFEngine.solve(mu, RHS, return_state = self.sample_state) vbMng(self, "DEL", "Done solving HF model.", 15) return u def _getSampleConcurrence(self, mu:paramVal, previous:Np1D) -> sampList: RROMPyAssert(self._mode, message = "Cannot add samples.") if not (self.sample_state or self.HFEngine.isCEye): raise RROMPyException(("Derivatives of solution with non-scalar " "C not computable.")) + from rrompy.utilities.numerical import dot if len(previous) >= len(self._derIdxs): self._derIdxs += nextDerivativeIndices(self._derIdxs, len(self.scaleFactor), len(previous) + 1 - len(self._derIdxs)) derIdx = self._derIdxs[len(previous)] mu = checkParameter(mu, self.HFEngine.npar) samplesOld = self.samples(previous) RHS = self.scaleDer(derIdx) * self.HFEngine.b(mu, derIdx) for j, derP in enumerate(self._derIdxs[: len(previous)]): diffP = [x - y for (x, y) in zip(derIdx, derP)] if np.all([x >= 0 for x in diffP]): RHS -= self.scaleDer(diffP) * dot(self.HFEngine.A(mu, diffP), samplesOld[j]) return self.solveLS(mu, RHS = RHS) def nextSample(self, mu:paramVal, overwrite : bool = False, postprocess : bool = True) -> Np1D: RROMPyAssert(self._mode, message = "Cannot add samples.") mu = checkParameter(mu, self.HFEngine.npar) muidxs = self.mus.findall(mu[0]) if len(muidxs) > 0: u = self._getSampleConcurrence(mu, np.sort(muidxs)) else: u = self.solveLS(mu) if postprocess: self.postprocessu(u, overwrite = overwrite) else: self.setsample(u, overwrite) if overwrite: self.mus[self.nsamples] = mu[0] else: self.mus.append(mu) self.nsamples += 1 return self.samples[self.nsamples - 1] def iterSample(self, mus:paramList) -> sampList: mus = checkParameterList(mus, self.HFEngine.npar) vbMng(self, "INIT", "Starting sampling iterations.", 5) n = len(mus) if n <= 0: raise RROMPyException(("Number of samples must be positive.")) self.resetHistory() if len(mus.unique()) != n: for j in range(n): vbMng(self, "MAIN", "Computing sample {} / {}.".format(j + 1, n), 7) self.nextSample(mus[j], overwrite = (j > 0), postprocess = False) if n > 1 and j == 0: self.preallocateSamples(self.samples[0], mus[0], n) else: self.setsample(self.solveLS(mus), overwrite = False) self.mus = copy(mus) self.nsamples = n self.postprocessuBulk() vbMng(self, "DEL", "Finished sampling iterations.", 5) return self.samples def plotSamples(self, warpings : List[List[callable]] = None, name : str = "u", **kwargs) -> Tuple[List[FigHandle], List[str]]: """ Do some nice plots of the samples. Args: warpings(optional): Domain warping functions. name(optional): Name to be shown as title of the plots. Defaults to 'u'. Returns: Output filenames and figure handles. """ if warpings is None: warpings = [None] * self.nsamples figs = [None] * self.nsamples filesOut = [None] * self.nsamples for j in range(self.nsamples): pltOut = self.HFEngine.plot(self.samples[j], warpings[j], self.sample_state, "{}_{}".format(name, j), **kwargs) if isinstance(pltOut, (tuple,)): figs[j], filesOut[j] = pltOut else: figs[j] = pltOut if filesOut[0] is None: return figs return figs, filesOut def outParaviewSamples(self, warpings : List[List[callable]] = None, name : str = "u", filename : str = "out", times : Np1D = None, **kwargs) -> List[str]: """ Output samples to ParaView file. Args: warpings(optional): Domain warping functions. name(optional): Base name to be used for data output. filename(optional): Name of output file. times(optional): Timestamps. Returns: Output filenames. """ if warpings is None: warpings = [None] * self.nsamples if times is None: times = [0.] * self.nsamples filesOut = [None] * self.nsamples for j in range(self.nsamples): filesOut[j] = self.HFEngine.outParaview( self.samples[j], warpings[j], self.sample_state, "{}_{}".format(name, j), "{}_{}".format(filename, j), times[j], **kwargs) if filesOut[0] is None: return None return filesOut def outParaviewTimeDomainSamples(self, omegas : Np1D = None, warpings : List[List[callable]] = None, timeFinal : Np1D = None, periodResolution : List[int] = 20, name : str = "u", filename : str = "out", **kwargs) -> List[str]: """ Output samples to ParaView file, converted to time domain. Args: omegas(optional): frequencies. timeFinal(optional): final time of simulation. periodResolution(optional): number of time steps per period. name(optional): Base name to be used for data output. filename(optional): Name of output file. Returns: Output filename. """ if omegas is None: omegas = np.real(self.mus) if warpings is None: warpings = [None] * self.nsamples if not isinstance(timeFinal, (list, tuple,)): timeFinal = [timeFinal] * self.nsamples if not isinstance(periodResolution, (list, tuple,)): periodResolution = [periodResolution] * self.nsamples filesOut = [None] * self.nsamples for j in range(self.nsamples): filesOut[j] = self.HFEngine.outParaviewTimeDomain( self.samples[j], omegas[j], warpings[j], self.sample_state, timeFinal[j], periodResolution[j], "{}_{}".format(name, j), "{}_{}".format(filename, j), **kwargs) if filesOut[0] is None: return None return filesOut diff --git a/rrompy/sampling/sample_list.py b/rrompy/sampling/sample_list.py index 3366fea..ffd35fd 100644 --- a/rrompy/sampling/sample_list.py +++ b/rrompy/sampling/sample_list.py @@ -1,224 +1,224 @@ # 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 deepcopy as copy import numpy as np from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.utilities.base.types import Np1D, List -from rrompy.utilities.numerical import dot __all__ = ['emptySampleList', 'sampleList'] def emptySampleList(): return sampleList(np.empty((0, 0))) class sampleList: def __init__(self, data:List[Np1D], lengthCheck : int = None, deep : bool = True): if isinstance(data, (self.__class__,)): data = data.data if isinstance(data, (np.ndarray,)): self.data = copy(data) if deep else data if self.data.ndim <= 1: self.data.shape = (self.data.shape[0], 1) else: if not isinstance(data, (list,)): data = [data] self.data = np.empty((len(data[0]), len(data)), dtype = data[0].dtype) for j, par in enumerate(data): self[j] = copy(data[j]) if deep else data[j] if j == 0 and lengthCheck is None: lengthCheck = self.shape[0] RROMPyAssert(len(data[j]), lengthCheck, "Number of parameters") def __len__(self): return self.shape[1] def __str__(self): return str(self.data) def __repr__(self): return repr(self.data) @property def shape(self): return self.data.shape @property def re(self): return sampleList(np.real(self.data)) @property def im(self): return sampleList(np.imag(self.data)) @property def abs(self): return sampleList(np.abs(self.data)) @property def angle(self): return sampleList(np.angle(self.data)) def conj(self): return sampleList(np.conj(self.data)) @property def T(self): return sampleList(self.data.T) @property def H(self): return sampleList(self.data.T.conj()) @property def dtype(self): return self.data.dtype @dtype.setter def dtype(self, dtype): self.data.dtype = dtype def __getitem__(self, key): return self.data[:, key] def __call__(self, key): return sampleList(self.data[:, key]) def __setitem__(self, key, value): if isinstance(value, self.__class__): value = value.data - if isinstance(key, (tuple, list,)): + if isinstance(key, (tuple, list, np.ndarray)): RROMPyAssert(len(key), len(value), "Slice length") for k, val in zip(key, value): self[k] = val else: self.data[:, key] = value.flatten() def __iter__(self): return self.data.T.__iter__() def __eq__(self, other): if not hasattr(other, "shape") or self.shape != other.shape: return False if isinstance(other, self.__class__): fac = other.data else: fac = other return np.allclose(self.data, fac) def __ne__(self, other): return not self == other def __copy__(self): return sampleList(self.data) def __deepcopy__(self, memo): return sampleList(copy(self.data, memo)) def __add__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data + fac) def __iadd__(self, other): self.data = (self + other).data return self def __sub__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data - fac) def __isub__(self, other): self.data = (self - other).data return self def __mul__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data * fac) def __imul__(self, other): self.data = (self * other).data return self def __truediv__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(self.data / fac) def __idiv__(self, other): self.data = (self / other).data return self def __pow__(self, other): if isinstance(other, self.__class__): RROMPyAssert(self.shape, other.shape, "Sample shape") fac = other.data else: fac = other return sampleList(np.power(self.data, fac)) def __ipow__(self, other): self.data = (self ** other).data return self def __neg__(self): return sampleList(- self.data) def __pos__(self): return sampleList(self.data) def reset(self, size, dtype = np.complex): self.data = np.empty(size, dtype = dtype) self.data[:] = np.nan def append(self, items): if isinstance(items, self.__class__): fac = items.data else: fac = np.array(items, ndmin = 2) self.data = np.append(self.data, fac, axis = 1) def pop(self, idx = -1): self.data = np.delete(self.data, idx, axis = 1) def dot(self, other, sampleListOut : bool = True): + from rrompy.utilities.numerical import dot if isinstance(other, self.__class__): other = other.data try: prod = dot(self.data, other) except: prod = dot(other.T, self.data.T).T if sampleListOut: prod = sampleList(prod) return prod diff --git a/rrompy/solver/norm_utilities.py b/rrompy/solver/norm_utilities.py index 358d53c..eccef9d 100644 --- a/rrompy/solver/norm_utilities.py +++ b/rrompy/solver/norm_utilities.py @@ -1,91 +1,93 @@ # 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 abc import abstractmethod import numpy as np from numbers import Number from copy import deepcopy as copy from rrompy.utilities.base.types import Np1D, Np2D, DictAny from rrompy.utilities.numerical import dot as tdot, solve as tsolve +from rrompy.sampling.sample_list import sampleList +from rrompy.parameter.parameter_list import parameterList from .linear_solver import setupSolver from rrompy.utilities.exception_manager import RROMPyException __all__ = ['Np2DLike', 'Np2DLikeInv', 'Np2DLikeInvLowRank', 'normEngine'] class Np2DLike: @abstractmethod def dot(self, u:Np2D) -> Np2D: pass class Np2DLikeInv(Np2DLike): def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny): self.K, self.M = K, M self.MH = np.conj(M) if isinstance(self.M, Number) else M.T.conj() try: self.solver, self.solverArgs = setupSolver(solverType, solverArgs) except: self.solver, self.solverArgs = solverType, solverArgs def dot(self, u:Np2D) -> Np2D: return tdot(self.MH, tsolve(self.K, tdot(self.M, u), self.solver, self.solverArgs)).reshape(u.shape) @property def shape(self): if isinstance(self.M, Number): return self.K.shape return (self.MH.shape[0], self.M.shape[1]) class Np2DLikeInvLowRank(Np2DLike): def __init__(self, K:Np2D, M:Np2D, solverType:str, solverArgs:DictAny, rank:int, oversampling : int = 10, seed : int = 420): sizeO = K.shape[1] if hasattr(K, "shape") else M.shape[1] if rank > sizeO: raise RROMPyException(("Cannot select compressed rank larger than " "original size.")) if oversampling < 0: raise RROMPyException("Oversampling parameter must be positive.") HF = Np2DLikeInv(K, M, solverType, solverArgs) np.random.seed(seed) xs = np.random.randn(sizeO, rank + oversampling) samples = HF.dot(xs) try: Q, _ = np.linalg.qr(samples, mode = "reduced") R = HF.dot(Q).T.conj() # assuming HF (i.e. K) hermitian... U, s, Vh = np.linalg.svd(R, full_matrices = False) self.L = Q.dot(U[:, : rank]) * s[: rank] self.R = Vh[: rank, :] except np.linalg.LinAlgError as e: raise RROMPyException(e) def dot(self, u:Np2D) -> Np2D: return tdot(self.L, tdot(self.R, u)).reshape(u.shape) @property def shape(self): return (self.L.shape[0], self.R.shape[1]) class normEngine: def __init__(self, normMatrix:Np2D): self.normMatrix = copy(normMatrix) def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D: - if not isinstance(u, (np.ndarray,)): u = u.data - if not isinstance(v, (np.ndarray,)): v = v.data + if isinstance(u, (parameterList, sampleList)): u = u.data + if isinstance(v, (parameterList, sampleList)): v = v.data if onlyDiag: return np.sum(tdot(self.normMatrix, u) * v.conj(), axis = 0) return tdot(tdot(self.normMatrix, u).T, v.conj()).T def norm(self, u:Np2D) -> Np1D: return np.power(np.abs(self.innerProduct(u, u, onlyDiag = True)), .5) diff --git a/rrompy/utilities/expression/expression_evaluator.py b/rrompy/utilities/expression/expression_evaluator.py index e07cecc..69128fc 100644 --- a/rrompy/utilities/expression/expression_evaluator.py +++ b/rrompy/utilities/expression/expression_evaluator.py @@ -1,124 +1,125 @@ # 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 numbers import Number import numpy as np from copy import deepcopy as copy from .keys import (expressionKeysUnary, expressionKeysUnaryParam, expressionKeysBinary, expressionKeysBinaryParam) from rrompy.utilities.base.types import Tuple, TupleAny, paramList from rrompy.utilities.exception_manager import RROMPyException -from rrompy.parameter import parameterList, checkParameterList +from rrompy.sampling.sample_list import sampleList +from rrompy.parameter.parameter_list import parameterList, checkParameterList __all__ = ["expressionEvaluator"] def manageNotExpression(expr): if isinstance(expr, (str,)) and expr == "x": return None elif isinstance(expr, (Number,)): return expr - elif isinstance(expr, (parameterList,)): return expr.data + elif isinstance(expr, (parameterList, sampleList)): return expr.data else: try: return np.array(expr) except: raise RROMPyException(("Expression '{}' not " "recognized.").format(expr)) def expressionEvaluator(expr:TupleAny, x:paramList, force_shape : Tuple[int] = None): if not isinstance(x, (parameterList,)): x = checkParameterList(x) exprSimp = [None] * len(expr) for j, ex in enumerate(expr): if isinstance(ex, (tuple,)): exprSimp[j] = expressionEvaluator(ex, x) else: exprSimp[j] = ex z, zc = None, None pile, pilePars = [], [] j = -1 while j + 1 < len(exprSimp): j += 1 ex = exprSimp[j] if not isinstance(ex, (np.ndarray, parameterList, list, tuple,)): if ex in expressionKeysUnary.keys(): pile = pile + [ex] pilePars = pilePars + [None] continue if ex in expressionKeysUnaryParam.keys(): pile = pile + [ex] j += 1 if j >= len(exprSimp) or not isinstance(exprSimp[j], (dict,)): raise RROMPyException(("Parameters missing for unary " "operand '{}'.").format(ex)) pilePars = pilePars + [exprSimp[j]] continue if ex in expressionKeysBinary.keys(): if len(pile) > 0 or z is None or zc is not None: raise RROMPyException(("Binary operand '{}' must follow " "numerical expression.").format(ex)) zc = copy(z) pile = pile + [ex] pilePars = pilePars + [None] continue if ex in expressionKeysBinaryParam.keys(): if len(pile) > 0 or z is None or zc is not None: raise RROMPyException(("Binary operand '{}' must follow " "numerical expression.").format(ex)) zc = copy(z) pile = pile + [ex] j += 1 if j >= len(exprSimp) or not isinstance(exprSimp[j], (dict,)): raise RROMPyException(("Parameters missing for binary " "operand '{}'.").format(ex)) pilePars = pilePars + [exprSimp[j]] continue z = manageNotExpression(ex) - if z is None: z = checkParameterList(x).data + if z is None: z = checkParameterList(x, return_data = True) if len(pile) > 0: for pl, plp in zip(pile[::-1], pilePars[::-1]): if zc is None: if plp is None: z = expressionKeysUnary[pl](z) else: z = expressionKeysUnaryParam[pl](z, plp) else: if plp is None: z = expressionKeysBinary[pl](zc, z) else: z = expressionKeysBinaryParam[pl](zc, z, plp) zc, pile, pilePars = None, [], [] if len(pile) > 0: raise RROMPyException(("Missing numerical expression feeding into " "'{}'.").format(pile[-1])) if force_shape is not None: if hasattr(z, "__len__") and len(z) > 1: - if isinstance(z, (parameterList,)): z = z.data + if isinstance(z, (parameterList, sampleList)): z = z.data if isinstance(z, (list, tuple,)): z = np.array(z) if z.size == np.prod(force_shape): z = np.reshape(z, force_shape) else: zdim = len(z.shape) if z.shape != force_shape[: zdim]: raise RROMPyException(("Error in reshaping result: shapes " "{} and {} not compatible.").format( z.shape, force_shape)) else: z = np.tile(z, [1] * zdim + force_shape[zdim :]) else: if hasattr(z, "__len__"): z = z[0] z = z * np.ones(force_shape) return z diff --git a/rrompy/utilities/numerical/marginalize_poly_list.py b/rrompy/utilities/numerical/marginalize_poly_list.py index d36262c..48fc089 100644 --- a/rrompy/utilities/numerical/marginalize_poly_list.py +++ b/rrompy/utilities/numerical/marginalize_poly_list.py @@ -1,79 +1,79 @@ # 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 . # import numpy as np from scipy.sparse import csr from rrompy.utilities.base.types import Np1D, Np2D, ListAny from rrompy.utilities.base import freepar as fp from .hash_derivative import (hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.parameter import checkParameter __all__ = ['marginalizePolyList'] def marginalizePolyList(objs:ListAny, marginalVals : Np1D = [fp], zeroObj : Np2D = 0., recompress : bool = True) -> ListAny: res = [] freeLocations = [] fixedLocations = [] muFixed = [] if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] for i, m in enumerate(marginalVals): if m == fp: freeLocations += [i] else: fixedLocations += [i] muFixed += [m] - muFixed = checkParameter(muFixed, len(fixedLocations)) + muFixed = checkParameter(muFixed, len(fixedLocations), return_data = True) if zeroObj == "auto": if isinstance(objs[0], np.ndarray): zeroObj = np.zeros_like(objs[0]) elif isinstance(objs[0], csr.csr_matrix): d = objs[0].shape[0] zeroObj = csr.csr_matrix(([], [], np.zeros(d + 1)), shape = objs[0].shape, dtype = objs[0].dtype) else: zeroObj = 0. for j, obj in enumerate(objs): derjBase = hashI(j, len(marginalVals)) jNew = hashD([derjBase[i] for i in freeLocations]) derjFixed = [derjBase[i] for i in fixedLocations] - obj = np.prod((muFixed ** derjFixed).data) * obj + obj = np.prod(muFixed ** derjFixed) * obj if jNew >= len(res): for _ in range(len(res), jNew): res += [zeroObj] res += [obj] else: res[jNew] = res[jNew] + obj if recompress: for re in res[::-1]: try: if isinstance(re, np.ndarray): iszero = np.allclose(re, zeroObj, atol = 2 * np.finfo(re.dtype).eps) elif isinstance(re, csr.csr_matrix): iszero = re.nnz == 0 else: break if not iszero: break except: break res.pop() return res diff --git a/rrompy/utilities/numerical/tensor_la.py b/rrompy/utilities/numerical/tensor_la.py index 5901bf6..188768b 100644 --- a/rrompy/utilities/numerical/tensor_la.py +++ b/rrompy/utilities/numerical/tensor_la.py @@ -1,47 +1,53 @@ # 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 . # import numpy as np from numbers import Number +from rrompy.sampling.sample_list import sampleList +from rrompy.parameter.parameter_list import parameterList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['dot', 'solve'] def dot(u, v): if isinstance(u, Number) or isinstance(v, Number): return u * v + if isinstance(u, (parameterList, sampleList)): u = u.data + if isinstance(v, (parameterList, sampleList)): v = v.data if u.shape[-1] == v.shape[0]: if isinstance(u, np.ndarray): return np.tensordot(u, v, 1) else: return u.dot(v) M = u.shape[-1] N = v.shape[0] // M rshape = u.shape[: -2] + (N * u.shape[-2],) + v.shape[1 :] return u.dot(v.reshape(M, -1)).reshape(rshape) def solve(A, b, solver, kwargs): try: if isinstance(A, Number): return b / A + if isinstance(A, (parameterList, sampleList)): A = A.data + if isinstance(b, (parameterList, sampleList)): b = b.data if A.shape[-1] == b.shape[0]: return solver(A, b, kwargs) M = A.shape[-1] N = b.shape[0] // M rshape = A.shape[: -2] + (N * A.shape[-2],) + b.shape[1 :] return solver(A, b.reshape(M, -1), kwargs).reshape(rshape) except np.linalg.LinAlgError as e: raise RROMPyException(e) diff --git a/rrompy/utilities/poly_fitting/nearest_neighbor/val.py b/rrompy/utilities/poly_fitting/nearest_neighbor/val.py index f071267..d23f3ff 100644 --- a/rrompy/utilities/poly_fitting/nearest_neighbor/val.py +++ b/rrompy/utilities/poly_fitting/nearest_neighbor/val.py @@ -1,41 +1,41 @@ # 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 . # import numpy as np from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.parameter import checkParameterList __all__ = ['polyval'] def polyval(x:paramList, cL:Np2D, supportPoints:paramList, nNeighbors : int = 1, directionalWeights : Np1D = None) -> Np2D: - supportPoints = checkParameterList(supportPoints) + supportPoints = checkParameterList(supportPoints, return_data = True) if directionalWeights is None: directionalWeights = np.ones(supportPoints.shape[1]) npar = supportPoints.shape[1] - x = checkParameterList(x, npar) - muDiff = (np.repeat(supportPoints.data.reshape((len(supportPoints), 1, - npar)), len(x), axis = 1) - - x.data) * directionalWeights + x = checkParameterList(x, npar, return_data = True) + muDiff = (np.repeat(supportPoints.reshape((len(supportPoints), 1, + npar)), len(x), axis = 1) - x + ) * directionalWeights dist = (np.sum(np.abs(muDiff) ** 2., axis = 2) + np.finfo(float).eps ** 2.) ** -.5 if len(dist) > nNeighbors: iOut = np.argpartition(dist, - nNeighbors, axis = 0)[: - nNeighbors] np.put_along_axis(dist, iOut, 0., 0) dist /= np.linalg.norm(dist, axis = 0, ord = 1) return np.moveaxis(np.tensordot(dist.T, cL, 1), 0, -1) diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py index 8f3256c..5a913a9 100644 --- a/rrompy/utilities/poly_fitting/radial_basis/val.py +++ b/rrompy/utilities/poly_fitting/radial_basis/val.py @@ -1,54 +1,53 @@ # 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 deepcopy as copy import numpy as np from .kernel import kernels from rrompy.utilities.poly_fitting.polynomial.val import polyval as pvP from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.parameter import checkParameterList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['polyval'] def polyval(x:paramList, cG:Np2D, cL:Np2D, supportPoints:paramList, directionalWeights:Np1D, basis:str) -> Np2D: - x = checkParameterList(x) + x = checkParameterList(x, return_data = True) basisp, basisr = basis.split("_") c = pvP(x, cG, basisp) - x = x.data try: radialker = kernels[basisr.upper()] except: raise RROMPyException("Radial basis not recognized.") supportPoints = checkParameterList(supportPoints) csh = copy(c.shape) if len(csh) == 1: c = c.reshape(1, -1) radialVal = np.zeros((len(supportPoints), len(x))) xDiff2V, xDiff2I, xDiff2J = np.zeros(0), [], [] for i in range(len(supportPoints)): xiD2Loc = np.sum(np.abs((x - supportPoints[i]) * directionalWeights) ** 2., axis = 1) xiD2Good = np.where(xiD2Loc <= radialker.threshold)[0] xDiff2V = np.append(xDiff2V, xiD2Loc[xiD2Good]) xDiff2I += [i] * len(xiD2Good) xDiff2J += list(xiD2Good) radialVal[xDiff2I, xDiff2J] = radialker(xDiff2V, apply_threshold = False) c += np.tensordot(cL, radialVal, (0, 0)) if len(csh) == 1: c = c.flatten() return c