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