UQ.Quantify fails with scipy==1.9.1

[ReddTea]

Hello, I noticed that the same script failed for two different machines, I compared libraries, found out the scipy version was different. Forced the installation of scipy 1.9.1 on the machine where it worked, and it reproduces the same error.

The part of the code that produces the error is

UQ = un.UncertaintyQuantification(model=un_model,
                                   parameters=parameters,
                                   logger_level=u'critical',
                                   logger_filename=unsaveplace+'uncertainpy.log')

undata = UQ.quantify(seed=10,
                     method='pc',
                     plot=None,
                     pc_method='collocation',
                     logger_level=u'critical',
                     figure_folder=unsaveplace+'figures',
                     data_folder=unsaveplace+'data',
                     single=False)

Displaying

File ~/anaconda3/lib/python3.9/site-packages/uncertainpy/uncertainty.py:415, in UncertaintyQuantification.quantify(self, method, pc_method, rosenblatt, uncertain_parameters, polynomial_order, nr_collocation_nodes, quadrature_order, nr_pc_mc_samples, nr_mc_samples, allow_incomplete, seed, single, plot, figure_folder, figureformat, save, data_folder, filename, **custom_kwargs)
    397         data = self.polynomial_chaos_single(uncertain_parameters=uncertain_parameters,
    398                                             method=pc_method,
    399                                             rosenblatt=rosenblatt,
   (...)
    411                                             filename=filename,
    412                                             **custom_kwargs)
    414     else:
--> 415         data = self.polynomial_chaos(uncertain_parameters=uncertain_parameters,
    416                                      method=pc_method,
    417                                      rosenblatt=rosenblatt,
    418                                      polynomial_order=polynomial_order,
    419                                      nr_collocation_nodes=nr_collocation_nodes,
    420                                      quadrature_order=quadrature_order,
    421                                      nr_pc_mc_samples=nr_pc_mc_samples,
    422                                      allow_incomplete=allow_incomplete,
    423                                      seed=seed,
    424                                      plot=plot,
    425                                      figure_folder=figure_folder,
    426                                      figureformat=figureformat,
    427                                      save=save,
    428                                      data_folder=data_folder,
    429                                      filename=filename,
    430                                      **custom_kwargs)
    432 elif method.lower() == "mc":
    433     if single:

File ~/anaconda3/lib/python3.9/site-packages/uncertainpy/uncertainty.py:702, in UncertaintyQuantification.polynomial_chaos(self, method, rosenblatt, uncertain_parameters, polynomial_order, nr_collocation_nodes, quadrature_order, nr_pc_mc_samples, allow_incomplete, seed, plot, figure_folder, figureformat, save, data_folder, filename, **custom_kwargs)
    697 if len(uncertain_parameters) > 20:
    698     raise RuntimeWarning("The number of uncertain parameters is high."
    699                          + "The Monte-Carlo method might be faster.")
--> 702 self.data = self.uncertainty_calculations.polynomial_chaos(
    703     method=method,
    704     rosenblatt=rosenblatt,
    705     uncertain_parameters=uncertain_parameters,
    706     polynomial_order=polynomial_order,
    707     nr_collocation_nodes=nr_collocation_nodes,
    708     quadrature_order=quadrature_order,
    709     nr_pc_mc_samples=nr_pc_mc_samples,
    710     allow_incomplete=allow_incomplete,
    711     seed=seed,
    712     **custom_kwargs
    713     )
    715 self.data.backend = self.backend
    717 if filename is None:

File ~/anaconda3/lib/python3.9/site-packages/uncertainpy/core/uncertainty_calculations.py:1383, in UncertaintyCalculations.polynomial_chaos(self, method, rosenblatt, uncertain_parameters, polynomial_order, nr_collocation_nodes, quadrature_order, nr_pc_mc_samples, allow_incomplete, seed, **custom_kwargs)
   1374 # TODO add support for more methods here by using
   1375 # try:
   1376 #     getattr(self, method)
   (...)
   1379 
   1380 else:
   1381     raise ValueError("No polynomial chaos method with name {}".format(method))
-> 1383 data = self.analyse_PCE(U_hat, distribution, data, nr_samples=nr_pc_mc_samples)
   1385 data.seed = seed
   1387 return data

File ~/anaconda3/lib/python3.9/site-packages/uncertainpy/core/uncertainty_calculations.py:1032, in UncertaintyCalculations.analyse_PCE(self, U_hat, distribution, data, nr_samples)
   1030 if feature in U_hat:
   1031     data[feature].mean = cp.E(U_hat[feature], distribution)
-> 1032     data[feature].variance = cp.Var(U_hat[feature], distribution)
   1034     samples = distribution.sample(nr_samples, "M")
   1036     if len(data.uncertain_parameters) > 1:

File ~/anaconda3/lib/python3.9/site-packages/chaospy/descriptives/variance.py:43, in Var(poly, dist, **kws)
     41 poly = poly - E(poly, dist, **kws)
     42 poly = numpoly.square(poly)
---> 43 return E(poly, dist, **kws)

File ~/anaconda3/lib/python3.9/site-packages/chaospy/descriptives/expected.py:42, in E(poly, dist, **kws)
     39 if poly.isconstant():
     40     return poly.tonumpy()
---> 42 moments = dist.mom(poly.exponents.T, **kws)
     43 if len(dist) == 1:
     44     moments = moments[0]

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/distribution.py:668, in Distribution.mom(self, K, allow_approx, **kwargs)
    666 K = K.reshape(dim, size)
    667 try:
--> 668     out = [self._get_mom(kdata) for kdata in K.T]
    669     logger.debug("%s: moment calculated successfully", str(self))
    670 except chaospy.UnsupportedFeature:

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/distribution.py:668, in <listcomp>(.0)
    666 K = K.reshape(dim, size)
    667 try:
--> 668     out = [self._get_mom(kdata) for kdata in K.T]
    669     logger.debug("%s: moment calculated successfully", str(self))
    670 except chaospy.UnsupportedFeature:

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/distribution.py:693, in Distribution._get_mom(self, kdata)
    691     parameters = self.get_parameters(idx=None, cache={}, assert_numerical=False)
    692 assert "idx" not in parameters, (self, parameters)
--> 693 ret_val = float(self._mom(kdata, **parameters))
    694 assert not isinstance(ret_val, Distribution), (self, ret_val)
    695 self._mom_cache[tuple(kdata)] = ret_val

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/operators/joint.py:169, in J._mom(self, kloc, index)
    167 kloc = kloc[self._rotation]
    168 for unique_idx in numpy.unique(index[self._rotation]):
--> 169     output *= self._owners[unique_idx][1]._get_mom(kloc[index == unique_idx])
    170 return output

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/distribution.py:693, in Distribution._get_mom(self, kdata)
    691     parameters = self.get_parameters(idx=None, cache={}, assert_numerical=False)
    692 assert "idx" not in parameters, (self, parameters)
--> 693 ret_val = float(self._mom(kdata, **parameters))
    694 assert not isinstance(ret_val, Distribution), (self, ret_val)
    695 self._mom_cache[tuple(kdata)] = ret_val

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/shift_scale.py:109, in ShiftScaleDistribution._mom(self, kloc, dist, shift, scale)
    106 poly = numpoly.sum(scale * poly, axis=-1) + shift
    107 poly = numpoly.set_dimensions(numpoly.prod(poly**kloc), len(self))
    108 out = sum(
--> 109     [
    110         dist._get_mom(key) * coeff
    111         for key, coeff in zip(poly.exponents, poly.coefficients)
    112     ]
    113 )
    114 return out

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/shift_scale.py:110, in <listcomp>(.0)
    106 poly = numpoly.sum(scale * poly, axis=-1) + shift
    107 poly = numpoly.set_dimensions(numpoly.prod(poly**kloc), len(self))
    108 out = sum(
    109     [
--> 110         dist._get_mom(key) * coeff
    111         for key, coeff in zip(poly.exponents, poly.coefficients)
    112     ]
    113 )
    114 return out

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/baseclass/distribution.py:693, in Distribution._get_mom(self, kdata)
    691     parameters = self.get_parameters(idx=None, cache={}, assert_numerical=False)
    692 assert "idx" not in parameters, (self, parameters)
--> 693 ret_val = float(self._mom(kdata, **parameters))
    694 assert not isinstance(ret_val, Distribution), (self, ret_val)
    695 self._mom_cache[tuple(kdata)] = ret_val

File ~/anaconda3/lib/python3.9/site-packages/chaospy/distributions/collection/trunc_normal.py:50, in trunc_normal._mom(self, n, a, b, mu, sigma)
     49 def _mom(self, n, a, b, mu, sigma):
---> 50     return truncnorm.moment(int(n), a, b, loc=mu, scale=sigma)

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:1371, in rv_generic.moment(self, order, *args, **kwds)
   1369     mu, mu2, g1, g2 = self._stats(*shapes, **mdict)
   1370 val = np.empty(loc.shape)  # val needs to be indexed by loc
-> 1371 val[...] = _moment_from_stats(n, mu, mu2, g1, g2, self._munp, shapes)
   1373 # Convert to transformed  X = L + S*Y
   1374 # E[X^n] = E[(L+S*Y)^n] = L^n sum(comb(n, k)*(S/L)^k E[Y^k], k=0...n)
   1375 result = zeros(i0.shape)

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_distn_infrastructure.py:393, in _moment_from_stats(n, mu, mu2, g1, g2, moment_func, args)
    391         val = mu4+4*mu*mu3+6*mu*mu*mu2+mu*mu*mu*mu
    392 else:
--> 393     val = moment_func(n, *args)
    395 return val

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8141, in truncnorm_gen._munp(self, n, a, b)
   8138         moments.append(mk)
   8139     return moments[-1]
-> 8141 return _lazywhere((n >= 0) & (a == a) & (b == b), (n, a, b),
   8142                   np.vectorize(n_th_moment, otypes=[np.float64]),
   8143                   np.nan)

File ~/anaconda3/lib/python3.9/site-packages/scipy/_lib/_util.py:69, in _lazywhere(cond, arrays, f, fillvalue, f2)
     67 tcode = np.mintypecode([a.dtype.char for a in arrays])
     68 out = np.full(np.shape(arrays[0]), fill_value=fillvalue, dtype=tcode)
---> 69 np.place(out, cond, f(*temp))
     70 if f2 is not None:
     71     temp = tuple(np.extract(~cond, arr) for arr in arrays)

File ~/anaconda3/lib/python3.9/site-packages/numpy/lib/function_base.py:2163, in vectorize.__call__(self, *args, **kwargs)
   2160     vargs = [args[_i] for _i in inds]
   2161     vargs.extend([kwargs[_n] for _n in names])
-> 2163 return self._vectorize_call(func=func, args=vargs)

File ~/anaconda3/lib/python3.9/site-packages/numpy/lib/function_base.py:2246, in vectorize._vectorize_call(self, func, args)
   2243 # Convert args to object arrays first
   2244 inputs = [asanyarray(a, dtype=object) for a in args]
-> 2246 outputs = ufunc(*inputs)
   2248 if ufunc.nout == 1:
   2249     res = asanyarray(outputs, dtype=otypes[0])

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8127, in truncnorm_gen._munp.<locals>.n_th_moment(n, a, b)
   8122 def n_th_moment(n, a, b):
   8123     """
   8124     Returns n-th moment. Defined only if n >= 0.
   8125     Function cannot broadcast due to the loop over n
   8126     """
-> 8127     pA, pB = self._pdf([a, b], a, b)
   8128     probs = [pA, -pB]
   8129     moments = [0, 1]

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8047, in truncnorm_gen._pdf(self, x, a, b)
   8046 def _pdf(self, x, a, b):
-> 8047     return np.exp(self._logpdf(x, a, b))

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:8050, in truncnorm_gen._logpdf(self, x, a, b)
   8049 def _logpdf(self, x, a, b):
-> 8050     return _norm_logpdf(x) - _log_gauss_mass(a, b)

File ~/anaconda3/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:271, in _norm_logpdf(x)
    270 def _norm_logpdf(x):
--> 271     return -x**2 / 2.0 - _norm_pdf_logC

TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'int'

Could it be done so it is compatible with newer versions of scipy?

Cheers,
Pablo