Page MenuHomec4science

marginalize.py
No OneTemporary

File Metadata

Created
Wed, Jul 10, 17:13

marginalize.py

# 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 <http://www.gnu.org/licenses/>.
#
from numpy import array, polynomial as po
from copy import deepcopy as copy
from rrompy.utilities.base.types import Np1D, Np2D
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
__all__ = ['polymarginalize']
def polymarginalize(c:Np2D, basis:str, marginalVals : Np1D = [fp],
nMarginal : int = None) -> Np1D:
if not hasattr(c, "ndim"): c = array(c)
ndim = c.ndim
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
marginalVals = list(marginalVals)
try:
polyvalbase = {"CHEBYSHEV" : po.chebyshev.chebval,
"LEGENDRE" : po.legendre.legval,
"MONOMIAL" : po.polynomial.polyval}[basis.upper()]
except:
raise RROMPyException("Polynomial basis not recognized.")
RROMPyAssert(ndim, len(marginalVals), "Marginalized variables")
marginalDims = []
for j in range(len(marginalVals)):
if marginalVals[j] == fp:
marginalDims += [c.shape[j]]
if nMarginal is not None and len(marginalDims) != nMarginal:
raise RROMPyException(("Exactly {} 'freepar' entries in marginalVals "
"must be provided.").format(nMarginal))
cEff = [copy(c)]
for d in range(ndim):
if marginalVals[d] != fp:
for dj in range(len(cEff)):
cEff[dj] = polyvalbase(marginalVals[d], cEff[dj],
tensor = False)
else:
cEff2 = []
for dj in range(len(cEff)):
cEff2 += list(cEff[dj])
cEff = copy(cEff2)
return array(cEff).reshape(tuple(marginalDims))

Event Timeline